From a59fa474f049d44608ca41f6d83e989cf9807cff Mon Sep 17 00:00:00 2001 From: TRasmussen <33480590+TillRasmussen@users.noreply.github.com> Date: Wed, 23 Feb 2022 06:37:02 +0100 Subject: [PATCH] seabed stress - remove if statements (#673) * refactor seabed_stress. Bit for bit * Removed if statement from stepu. Results are binary identical, however taubx and tauby is updated on all iterations instead of just the last one. Not used within iteration * changed capping from logical to numeric in order to remove if statement. Moved call to deformation out of loop * clean dyn_finish, correct intent(inout) to intent(in) for Cw, resse Cb in stepu, remove if from seabed_stress_LKD * Reolve conflicts after updating main * modified environment for Freya to accomodate for additional OMP commands * Requested changes after review. Only changed in seabed stress and not bit for bit if cor=0.0 added legacy comment in ice_dyn_finish * move deformation to subcycling --- cicecore/cicedynB/dynamics/ice_dyn_eap.F90 | 35 +++-- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 119 ++++++++--------- cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 | 6 +- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 120 ++++++++---------- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 32 ++--- configuration/scripts/machines/env.freya_gnu | 3 +- .../scripts/machines/env.freya_intel | 1 + 7 files changed, 142 insertions(+), 174 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index 5324a9a8f..eaaa54cd1 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -399,29 +399,29 @@ subroutine eap (dt) !----------------------------------------------------------------- if (seabed_stress) then + if ( seabed_stress_method == 'LKD' ) then + !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) + do iblk = 1, nblocks + call seabed_stress_factor_LKD (nx_block, ny_block, & + icellu (iblk), & + indxui(:,iblk), indxuj(:,iblk), & + vice(:,:,iblk), aice(:,:,iblk), & + hwater(:,:,iblk), Tbu(:,:,iblk)) + enddo + !$OMP END PARALLEL DO - !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) - do iblk = 1, nblocks - - if ( seabed_stress_method == 'LKD' ) then - - call seabed_stress_factor_LKD (nx_block, ny_block, & - icellu (iblk), & - indxui(:,iblk), indxuj(:,iblk), & - vice(:,:,iblk), aice(:,:,iblk), & - hwater(:,:,iblk), Tbu(:,:,iblk)) - - elseif ( seabed_stress_method == 'probabilistic' ) then + elseif ( seabed_stress_method == 'probabilistic' ) then + !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) + do iblk = 1, nblocks call seabed_stress_factor_prob (nx_block, ny_block, & icellt(iblk), indxti(:,iblk), indxtj(:,iblk), & icellu(iblk), indxui(:,iblk), indxuj(:,iblk), & aicen(:,:,:,iblk), vicen(:,:,:,iblk), & hwater(:,:,iblk), Tbu(:,:,iblk)) - endif - - enddo - !$OMP END PARALLEL DO + enddo + !$OMP END PARALLEL DO + endif endif do ksub = 1,ndte ! subcycling @@ -476,7 +476,6 @@ subroutine eap (dt) call stepu (nx_block, ny_block, & icellu (iblk), Cdn_ocn (:,:,iblk), & indxui (:,iblk), indxuj (:,iblk), & - ksub, & aiu (:,:,iblk), strtmp (:,:,:), & uocn (:,:,iblk), vocn (:,:,iblk), & waterx (:,:,iblk), watery (:,:,iblk), & @@ -546,8 +545,6 @@ subroutine eap (dt) uvel (:,:,iblk), vvel (:,:,iblk), & uocn (:,:,iblk), vocn (:,:,iblk), & aiu (:,:,iblk), fm (:,:,iblk), & - strintx (:,:,iblk), strinty (:,:,iblk), & - strairx (:,:,iblk), strairy (:,:,iblk), & strocnx (:,:,iblk), strocny (:,:,iblk), & strocnxT(:,:,iblk), strocnyT(:,:,iblk)) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 1e5cdecff..d54a73dd4 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -98,7 +98,7 @@ subroutine evp (dt) use ice_dyn_evp_1d, only: ice_dyn_evp_1d_copyin, ice_dyn_evp_1d_kernel, & ice_dyn_evp_1d_copyout use ice_dyn_shared, only: evp_algorithm, stack_velocity_field, unstack_velocity_field - + use ice_dyn_shared, only: deformations real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -330,29 +330,29 @@ subroutine evp (dt) !----------------------------------------------------------------- if (seabed_stress) then - - !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) - do iblk = 1, nblocks - - if ( seabed_stress_method == 'LKD' ) then - + if ( seabed_stress_method == 'LKD' ) then + !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) + do iblk = 1, nblocks call seabed_stress_factor_LKD (nx_block, ny_block, & icellu (iblk), & indxui(:,iblk), indxuj(:,iblk), & vice(:,:,iblk), aice(:,:,iblk), & hwater(:,:,iblk), Tbu(:,:,iblk)) + enddo + !$OMP END PARALLEL DO - elseif ( seabed_stress_method == 'probabilistic' ) then - + elseif ( seabed_stress_method == 'probabilistic' ) then + !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) + do iblk = 1, nblocks call seabed_stress_factor_prob (nx_block, ny_block, & icellt(iblk), indxti(:,iblk), indxtj(:,iblk), & icellu(iblk), indxui(:,iblk), indxuj(:,iblk), & aicen(:,:,:,iblk), vicen(:,:,:,iblk), & hwater(:,:,iblk), Tbu(:,:,iblk)) - endif - - enddo - !$OMP END PARALLEL DO + enddo + !$OMP END PARALLEL DO + + endif endif call ice_timer_start(timer_evp_2d) @@ -401,27 +401,41 @@ subroutine evp (dt) do iblk = 1, nblocks ! if (trim(yield_curve) == 'ellipse') then - call stress (nx_block, ny_block, & - ksub, icellt(iblk), & - indxti (:,iblk), indxtj (:,iblk), & - uvel (:,:,iblk), vvel (:,:,iblk), & - dxt (:,:,iblk), dyt (:,:,iblk), & - dxhy (:,:,iblk), dyhx (:,:,iblk), & - cxp (:,:,iblk), cyp (:,:,iblk), & - cxm (:,:,iblk), cym (:,:,iblk), & - tarear (:,:,iblk), tinyarea (:,:,iblk), & - strength (:,:,iblk), & - stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & - stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & - stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & - stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & - stress12_1(:,:,iblk), stress12_2(:,:,iblk), & - stress12_3(:,:,iblk), stress12_4(:,:,iblk), & - shear (:,:,iblk), divu (:,:,iblk), & - rdg_conv (:,:,iblk), rdg_shear (:,:,iblk), & - strtmp (:,:,:) ) -! endif ! yield_curve + call stress (nx_block, ny_block, & + icellt(iblk), & + indxti (:,iblk), indxtj (:,iblk), & + uvel (:,:,iblk), vvel (:,:,iblk), & + dxt (:,:,iblk), dyt (:,:,iblk), & + dxhy (:,:,iblk), dyhx (:,:,iblk), & + cxp (:,:,iblk), cyp (:,:,iblk), & + cxm (:,:,iblk), cym (:,:,iblk), & + tinyarea (:,:,iblk), & + strength (:,:,iblk), & + stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & + stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & + stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & + stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & + stress12_1(:,:,iblk), stress12_2(:,:,iblk), & + stress12_3(:,:,iblk), stress12_4(:,:,iblk), & + strtmp (:,:,:) ) +! endif + !----------------------------------------------------------------- + ! on last subcycle, save quantities for mechanical redistribution + !----------------------------------------------------------------- + if (ksub == ndte) then + call deformations (nx_block, ny_block , & + icellt(iblk) , & + indxti(:,iblk) , indxtj(:,iblk) , & + uvel(:,:,iblk) , vvel(:,:,iblk) , & + dxt(:,:,iblk) , dyt(:,:,iblk) , & + cxp(:,:,iblk) , cyp(:,:,iblk) , & + cxm(:,:,iblk) , cym(:,:,iblk) , & + tarear(:,:,iblk) , & + shear(:,:,iblk) , divu(:,:,iblk) , & + rdg_conv(:,:,iblk), rdg_shear(:,:,iblk) ) + endif + !----------------------------------------------------------------- ! momentum equation !----------------------------------------------------------------- @@ -429,7 +443,6 @@ subroutine evp (dt) call stepu (nx_block, ny_block, & icellu (iblk), Cdn_ocn (:,:,iblk), & indxui (:,iblk), indxuj (:,iblk), & - ksub, & aiu (:,:,iblk), strtmp (:,:,:), & uocn (:,:,iblk), vocn (:,:,iblk), & waterx (:,:,iblk), watery (:,:,iblk), & @@ -547,8 +560,6 @@ subroutine evp (dt) uvel (:,:,iblk), vvel (:,:,iblk), & uocn (:,:,iblk), vocn (:,:,iblk), & aiu (:,:,iblk), fm (:,:,iblk), & - strintx (:,:,iblk), strinty (:,:,iblk), & - strairx (:,:,iblk), strairy (:,:,iblk), & strocnx (:,:,iblk), strocny (:,:,iblk), & strocnxT(:,:,iblk), strocnyT(:,:,iblk)) @@ -571,14 +582,14 @@ end subroutine evp ! author: Elizabeth C. Hunke, LANL subroutine stress (nx_block, ny_block, & - ksub, icellt, & + icellt, & indxti, indxtj, & uvel, vvel, & dxt, dyt, & dxhy, dyhx, & cxp, cyp, & cxm, cym, & - tarear, tinyarea, & + tinyarea, & strength, & stressp_1, stressp_2, & stressp_3, stressp_4, & @@ -586,15 +597,12 @@ subroutine stress (nx_block, ny_block, & stressm_3, stressm_4, & stress12_1, stress12_2, & stress12_3, stress12_4, & - shear, divu, & - rdg_conv, rdg_shear, & str ) use ice_dyn_shared, only: strain_rates, deformations, viscous_coeffs_and_rep_pressure integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - ksub , & ! subcycling step icellt ! no. of cells where icetmask = 1 integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & @@ -613,7 +621,6 @@ subroutine stress (nx_block, ny_block, & cxp , & ! 1.5*HTN - 0.5*HTS cym , & ! 0.5*HTE - 1.5*HTW cxm , & ! 0.5*HTN - 1.5*HTS - tarear , & ! 1/tarea tinyarea ! puny*tarea real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & @@ -621,12 +628,6 @@ subroutine stress (nx_block, ny_block, & stressm_1, stressm_2, stressm_3, stressm_4 , & ! sigma11-sigma22 stress12_1,stress12_2,stress12_3,stress12_4 ! sigma12 - real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & - shear , & ! strain rate II component (1/s) - divu , & ! strain rate I component, velocity divergence (1/s) - rdg_conv , & ! convergence term for ridging (1/s) - rdg_shear ! shear term for ridging (1/s) - real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(out) :: & str ! stress combinations @@ -654,8 +655,8 @@ subroutine stress (nx_block, ny_block, & str12ew, str12we, str12ns, str12sn , & strp_tmp, strm_tmp, tmp - logical :: capping ! of the viscous coef - + real(kind=dbl_kind),parameter :: capping = c1 ! of the viscous coef + character(len=*), parameter :: subname = '(stress)' !----------------------------------------------------------------- @@ -663,7 +664,6 @@ subroutine stress (nx_block, ny_block, & !----------------------------------------------------------------- str(:,:,:) = c0 - capping = .true. ! could be later included in ice_in do ij = 1, icellt i = indxti(ij) @@ -869,23 +869,6 @@ subroutine stress (nx_block, ny_block, & enddo ! ij - !----------------------------------------------------------------- - ! on last subcycle, save quantities for mechanical redistribution - !----------------------------------------------------------------- - if (ksub == ndte) then - call deformations (nx_block , ny_block , & - icellt , & - indxti , indxtj , & - uvel , vvel , & - dxt , dyt , & - cxp , cyp , & - cxm , cym , & - tarear , & - shear , divu , & - rdg_conv , rdg_shear ) - - endif - end subroutine stress !======================================================================= diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 index 99e0aee85..43bbe41a0 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 @@ -869,10 +869,8 @@ subroutine stepu_last(NA_len, rhow, lb, ub, Cw, aiu, uocn, vocn, & vvel(iw) = (cca * cc2 - ccb * cc1) / ab2 ! calculate seabed stress component for outputs - if (seabed_stress) then - taubx(iw) = -uvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0) - tauby(iw) = -vvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0) - end if + taubx(iw) = -uvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0) + tauby(iw) = -vvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0) end do #ifdef _OPENACC diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 8946e644a..76d0caf41 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -176,7 +176,7 @@ subroutine init_dyn (dt) if (trim(coriolis) == 'constant') then fcor_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s else if (trim(coriolis) == 'zero') then - fcor_blk(i,j,iblk) = 0.0 + fcor_blk(i,j,iblk) = c0 else fcor_blk(i,j,iblk) = c2*omega*sin(ULAT(i,j,iblk)) ! 1/s endif @@ -627,7 +627,6 @@ end subroutine dyn_prep2 subroutine stepu (nx_block, ny_block, & icellu, Cw, & indxui, indxuj, & - ksub, & aiu, str, & uocn, vocn, & waterx, watery, & @@ -642,8 +641,7 @@ subroutine stepu (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu, & ! total count when iceumask is true - ksub ! subcycling iteration + icellu ! total count when iceumask is true integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & indxui , & ! compressed index in i-direction @@ -677,7 +675,7 @@ subroutine stepu (nx_block, ny_block, & taubx , & ! seabed stress, x-direction (N/m^2) tauby ! seabed stress, y-direction (N/m^2) - real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & Cw ! ocean-ice neutral drag coefficient ! local variables @@ -741,14 +739,10 @@ subroutine stepu (nx_block, ny_block, & uvel(i,j) = (cca*cc1 + ccb*cc2) / ab2 ! m/s vvel(i,j) = (cca*cc2 - ccb*cc1) / ab2 - ! calculate seabed stress component for outputs - if (ksub == ndte) then ! on last subcycling iteration - if ( seabed_stress ) then - taubx(i,j) = -uvel(i,j)*Tbu(i,j) / (sqrt(uold**2 + vold**2) + u0) - tauby(i,j) = -vvel(i,j)*Tbu(i,j) / (sqrt(uold**2 + vold**2) + u0) - endif - endif - + ! calculate seabed stress component for outputs + ! only needed on last iteration. + taubx(i,j) = -uvel(i,j)*Cb + tauby(i,j) = -vvel(i,j)*Cb enddo ! ij end subroutine stepu @@ -766,8 +760,8 @@ subroutine dyn_finish (nx_block, ny_block, & uvel, vvel, & uocn, vocn, & aiu, fm, & - strintx, strinty, & - strairx, strairy, & +! strintx, strinty, & +! strairx, strairy, & strocnx, strocny, & strocnxT, strocnyT) @@ -785,11 +779,11 @@ subroutine dyn_finish (nx_block, ny_block, & uocn , & ! ocean current, x-direction (m/s) vocn , & ! ocean current, y-direction (m/s) aiu , & ! ice fraction on u-grid - fm , & ! Coriolis param. * mass in U-cell (kg/s) - strintx , & ! divergence of internal ice stress, x (N/m^2) - strinty , & ! divergence of internal ice stress, y (N/m^2) - strairx , & ! stress on ice by air, x-direction - strairy ! stress on ice by air, y-direction + fm ! Coriolis param. * mass in U-cell (kg/s) +! strintx , & ! divergence of internal ice stress, x (N/m^2) +! strinty , & ! divergence of internal ice stress, y (N/m^2) +! strairx , & ! stress on ice by air, x-direction +! strairy ! stress on ice by air, y-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & strocnx , & ! ice-ocean stress, x-direction @@ -799,14 +793,15 @@ subroutine dyn_finish (nx_block, ny_block, & strocnxT, & ! ice-ocean stress, x-direction strocnyT ! ice-ocean stress, y-direction + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & + Cw ! ocean-ice neutral drag coefficient + ! local variables integer (kind=int_kind) :: & i, j, ij real (kind=dbl_kind) :: vrel, rhow - real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & - Cw ! ocean-ice neutral drag coefficient character(len=*), parameter :: subname = '(dyn_finish)' @@ -897,9 +892,10 @@ subroutine seabed_stress_factor_LKD (nx_block, ny_block, & Tbu ! seabed stress factor (N/m^2) real (kind=dbl_kind) :: & - au, & ! concentration of ice at u location - hu, & ! volume per unit area of ice at u location (mean thickness, m) - hwu, & ! water depth at u location (m) + au , & ! concentration of ice at u location + hu , & ! volume per unit area of ice at u location (mean thickness, m) + hwu , & ! water depth at u location (m) + docalc_tbu, & ! logical as real (C0,C1) decides whether c0 is 0 or hcu ! critical thickness at u location (m) integer (kind=int_kind) :: & @@ -915,18 +911,19 @@ subroutine seabed_stress_factor_LKD (nx_block, ny_block, & hwu = min(hwater(i,j),hwater(i+1,j),hwater(i,j+1),hwater(i+1,j+1)) - if (hwu < threshold_hw) then + docalc_tbu = merge(c1,c0,hwu < threshold_hw) + - au = max(aice(i,j),aice(i+1,j),aice(i,j+1),aice(i+1,j+1)) - hu = max(vice(i,j),vice(i+1,j),vice(i,j+1),vice(i+1,j+1)) + au = max(aice(i,j),aice(i+1,j),aice(i,j+1),aice(i+1,j+1)) + hu = max(vice(i,j),vice(i+1,j),vice(i,j+1),vice(i+1,j+1)) - ! 1- calculate critical thickness - hcu = au * hwu / k1 + ! 1- calculate critical thickness + hcu = au * hwu / k1 - ! 2- calculate seabed stress factor - Tbu(i,j) = k2 * max(c0,(hu - hcu)) * exp(-alphab * (c1 - au)) + ! 2- calculate seabed stress factor + Tbu(i,j) = docalc_tbu*k2 * max(c0,(hu - hcu)) * exp(-alphab * (c1 - au)) - endif +! endif enddo ! ij @@ -1392,7 +1389,7 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & real (kind=dbl_kind), intent(in):: & Deltane, Deltanw, Deltasw, Deltase ! Delta at each corner - logical , intent(in):: capping + real(kind=dbl_kind) , intent(in):: capping real (kind=dbl_kind), intent(out):: & zetax2ne, zetax2nw, zetax2sw, zetax2se, & ! zetax2 at each corner @@ -1404,43 +1401,38 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & tmpcalcne, tmpcalcnw, tmpcalcsw, tmpcalcse ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code - -! if (trim(yield_curve) == 'ellipse') then - - if (capping) then - tmpcalcne = strength/max(Deltane,tinyarea) - tmpcalcnw = strength/max(Deltanw,tinyarea) - tmpcalcsw = strength/max(Deltasw,tinyarea) - tmpcalcse = strength/max(Deltase,tinyarea) - else - tmpcalcne = strength/(Deltane + tinyarea) - tmpcalcnw = strength/(Deltanw + tinyarea) - tmpcalcsw = strength/(Deltasw + tinyarea) - tmpcalcse = strength/(Deltase + tinyarea) - endif - zetax2ne = (c1+Ktens)*tmpcalcne ! northeast - rep_prsne = (c1-Ktens)*tmpcalcne*Deltane - etax2ne = epp2i*zetax2ne + ! if (trim(yield_curve) == 'ellipse') then + tmpcalcne = capping *(strength/max(Deltane, tinyarea))+ & + (c1-capping)* strength/ (Deltane+ tinyarea) + tmpcalcnw = capping *(strength/max(Deltanw, tinyarea))+ & + (c1-capping)* strength/ (Deltanw+ tinyarea) + tmpcalcsw = capping *(strength/max(Deltasw, tinyarea))+ & + (c1-capping)* strength/ (Deltasw+ tinyarea) + tmpcalcse = capping *(strength/max(Deltase, tinyarea))+ & + (c1-capping)* strength/ (Deltase+ tinyarea) + + zetax2ne = (c1+Ktens)*tmpcalcne ! northeast + rep_prsne = (c1-Ktens)*tmpcalcne*Deltane + etax2ne = epp2i*zetax2ne - zetax2nw = (c1+Ktens)*tmpcalcnw ! northwest - rep_prsnw = (c1-Ktens)*tmpcalcnw*Deltanw - etax2nw = epp2i*zetax2nw + zetax2nw = (c1+Ktens)*tmpcalcnw ! northwest + rep_prsnw = (c1-Ktens)*tmpcalcnw*Deltanw + etax2nw = epp2i*zetax2nw - zetax2sw = (c1+Ktens)*tmpcalcsw ! southwest - rep_prssw = (c1-Ktens)*tmpcalcsw*Deltasw - etax2sw = epp2i*zetax2sw + zetax2sw = (c1+Ktens)*tmpcalcsw ! southwest + rep_prssw = (c1-Ktens)*tmpcalcsw*Deltasw + etax2sw = epp2i*zetax2sw - zetax2se = (c1+Ktens)*tmpcalcse ! southeast - rep_prsse = (c1-Ktens)*tmpcalcse*Deltase - etax2se = epp2i*zetax2se - -! else + zetax2se = (c1+Ktens)*tmpcalcse ! southeast + rep_prsse = (c1-Ktens)*tmpcalcse*Deltase + etax2se = epp2i*zetax2se + ! else -! endif + ! endif end subroutine viscous_coeffs_and_rep_pressure - + !======================================================================= ! Load velocity components into array for boundary updates diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 6fa6ab108..a9da7e300 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -440,33 +440,35 @@ subroutine implicit_solver (dt) !----------------------------------------------------------------- ! seabed stress factor Tbu (Tbu is part of Cb coefficient) !----------------------------------------------------------------- - - if (seabed_stress) then - - !$OMP PARALLEL DO PRIVATE(iblk) - do iblk = 1, nblocks - - if ( seabed_stress_method == 'LKD' ) then + if (seabed_stress) then + if ( seabed_stress_method == 'LKD' ) then + !$OMP PARALLEL DO PRIVATE(iblk) + do iblk = 1, nblocks call seabed_stress_factor_LKD (nx_block, ny_block, & icellu (iblk), & indxui(:,iblk), indxuj(:,iblk), & vice(:,:,iblk), aice(:,:,iblk), & hwater(:,:,iblk), Tbu(:,:,iblk)) + enddo + !$OMP END PARALLEL DO - elseif ( seabed_stress_method == 'probabilistic' ) then + elseif ( seabed_stress_method == 'probabilistic' ) then + !$OMP PARALLEL DO PRIVATE(iblk) + do iblk = 1, nblocks call seabed_stress_factor_prob (nx_block, ny_block, & icellt(iblk), indxti(:,iblk), indxtj(:,iblk), & icellu(iblk), indxui(:,iblk), indxuj(:,iblk), & aicen(:,:,:,iblk), vicen(:,:,:,iblk), & hwater(:,:,iblk), Tbu(:,:,iblk)) - endif + enddo + !$OMP END PARALLEL DO - enddo - !$OMP END PARALLEL DO + endif endif + !----------------------------------------------------------------- ! calc size of problem (ntot) and allocate solution vector !----------------------------------------------------------------- @@ -640,8 +642,6 @@ subroutine implicit_solver (dt) uvel (:,:,iblk), vvel (:,:,iblk), & uocn (:,:,iblk), vocn (:,:,iblk), & aiu (:,:,iblk), fm (:,:,iblk), & - strintx (:,:,iblk), strinty (:,:,iblk), & - strairx (:,:,iblk), strairy (:,:,iblk), & strocnx (:,:,iblk), strocny (:,:,iblk), & strocnxT(:,:,iblk), strocnyT(:,:,iblk)) @@ -1188,14 +1188,10 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & stressp_1, stressp_2, stressp_3, stressp_4 , & strp_tmp - logical :: capping ! of the viscous coeff + real(kind=dbl_kind),parameter :: capping = c0 ! of the viscous coef character(len=*), parameter :: subname = '(calc_zeta_dPr)' - ! Initialize - - capping = .false. - ! Initialize stPr, zetax2 and etax2 to zero ! (for cells where icetmask is false) stPr = c0 diff --git a/configuration/scripts/machines/env.freya_gnu b/configuration/scripts/machines/env.freya_gnu index b655d6dd0..2681e1318 100755 --- a/configuration/scripts/machines/env.freya_gnu +++ b/configuration/scripts/machines/env.freya_gnu @@ -8,7 +8,7 @@ endif if ("$inp" != "-nomodules") then source /opt/modules/default/init/csh # Initialize modules for csh - Clear environment +# Clear environment module rm PrgEnv-intel module rm PrgEnv-cray module rm PrgEnv-gnu @@ -37,3 +37,4 @@ setenv ICE_MACHINE_ACCT P0000000 setenv ICE_MACHINE_QUEUE "development" setenv ICE_MACHINE_BLDTHRDS 18 setenv ICE_MACHINE_QSTAT "qstat " +setenv OMP_STACKSIZE 64M diff --git a/configuration/scripts/machines/env.freya_intel b/configuration/scripts/machines/env.freya_intel index dcbc1f8ba..4b45cd9e7 100755 --- a/configuration/scripts/machines/env.freya_intel +++ b/configuration/scripts/machines/env.freya_intel @@ -36,3 +36,4 @@ setenv ICE_MACHINE_ACCT P0000000 setenv ICE_MACHINE_QUEUE "development" setenv ICE_MACHINE_BLDTHRDS 18 setenv ICE_MACHINE_QSTAT "qstat " +setenv OMP_STACKSIZE 64M