From 5c036cc442904dd8c5900e75eca45b53e5bf3827 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Tue, 7 Dec 2021 23:49:29 +0000 Subject: [PATCH 1/8] refactor seabed_stress. Bit for bit --- cicecore/cicedynB/dynamics/ice_dyn_eap.F90 | 44 ++++++++++--------- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 23 +++++----- cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 | 6 +-- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 9 ++-- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 24 +++++----- 5 files changed, 53 insertions(+), 53 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index 83374d4dd..759831eb7 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -395,30 +395,32 @@ subroutine eap (dt) !----------------------------------------------------------------- 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 - !$OMP PARALLEL DO PRIVATE(iblk) - 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 - - 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 + elseif ( seabed_stress_method == 'probabilistic' ) then + !$OMP PARALLEL DO PRIVATE(iblk) + do iblk = 1, nblocks - enddo - !$OMP END PARALLEL DO + 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)) + enddo + !$OMP END PARALLEL DO + + endif endif + do ksub = 1,ndte ! subcycling diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 8f3fc4910..6bc6f2ef8 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -332,29 +332,30 @@ subroutine evp (dt) !----------------------------------------------------------------- if (seabed_stress) then - - !$OMP PARALLEL DO PRIVATE(iblk) - do iblk = 1, nblocks - - if ( seabed_stress_method == 'LKD' ) 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 call ice_timer_start(timer_evp_2d) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 index c691453cb..4e2a0e2e6 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 23251b2d1..8fec2e259 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -742,13 +742,10 @@ subroutine stepu (nx_block, ny_block, & 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 + if (ksub == ndte) 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 - enddo ! ij end subroutine stepu diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 2f1285084..8e1159ff3 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 !----------------------------------------------------------------- From 4ab3476f4b6985c9ad8a661d319ec0359c479c04 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Sun, 12 Dec 2021 21:37:32 +0000 Subject: [PATCH 2/8] 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 --- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 8fec2e259..3a3d0dc6e 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -741,11 +741,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 - 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 + ! calculate seabed stress component for outputs + ! only needed on last iteration. + 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) enddo ! ij end subroutine stepu From 0122dd33b8e8f527bf5f565a1b8b46a337d02c5f Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Mon, 20 Dec 2021 12:23:21 +0000 Subject: [PATCH 3/8] changed capping from logical to numeric in order to remove if statement. Moved call to deformation out of loop --- cicecore/cicedynB/dynamics/ice_dyn_eap.F90 | 1 - cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 192 +++++++++++------- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 59 +++--- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 6 +- 4 files changed, 140 insertions(+), 118 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index 759831eb7..d4b602b70 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -473,7 +473,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), & diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 6bc6f2ef8..3da528e95 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 @@ -393,26 +393,25 @@ subroutine evp (dt) divu,rdg_conv,rdg_shear,shear,taubx,tauby ) else ! evp_algorithm == standard_2d (Standard CICE) + if (ndte > 1) then + do ksub = 1,ndte-1 ! subcycling - do ksub = 1,ndte ! subcycling - - !----------------------------------------------------------------- - ! stress tensor equation, total surface stress - !----------------------------------------------------------------- + !----------------------------------------------------------------- + ! stress tensor equation, total surface stress + !----------------------------------------------------------------- - !$TCXOMP PARALLEL DO PRIVATE(iblk,strtmp) - do iblk = 1, nblocks + !$TCXOMP PARALLEL DO PRIVATE(iblk,strtmp) + do iblk = 1, nblocks -! if (trim(yield_curve) == 'ellipse') then call stress (nx_block, ny_block, & - ksub, icellt(iblk), & + 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), & + tinyarea (:,:,iblk), & strength (:,:,iblk), & stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & @@ -420,54 +419,121 @@ subroutine evp (dt) 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 (:,:,:) ) + + !----------------------------------------------------------------- + ! momentum equation + !----------------------------------------------------------------- + + call stepu (nx_block, ny_block, & + icellu (iblk), Cdn_ocn (:,:,iblk), & + indxui (:,iblk), indxuj (:,iblk), & + aiu (:,:,iblk), strtmp (:,:,:), & + uocn (:,:,iblk), vocn (:,:,iblk), & + waterx (:,:,iblk), watery (:,:,iblk), & + forcex (:,:,iblk), forcey (:,:,iblk), & + umassdti (:,:,iblk), fm (:,:,iblk), & + uarear (:,:,iblk), & + strintx (:,:,iblk), strinty (:,:,iblk), & + taubx (:,:,iblk), tauby (:,:,iblk), & + uvel_init(:,:,iblk), vvel_init(:,:,iblk),& + uvel (:,:,iblk), vvel (:,:,iblk), & + Tbu (:,:,iblk)) + + enddo + !$TCXOMP END PARALLEL DO + + call stack_velocity_field(uvel, vvel, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, uvel, vvel) + + enddo ! subcycling + endif ! only subcycle when ndte larger than 1 + !$TCXOMP PARALLEL DO PRIVATE(iblk,strtmp) + do iblk = 1, nblocks + +! if (trim(yield_curve) == 'ellipse') then + 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 ! yield_curve + !----------------------------------------------------------------- + ! on last subcycle, save quantities for mechanical redistribution + !----------------------------------------------------------------- + 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) ) + !----------------------------------------------------------------- ! momentum equation !----------------------------------------------------------------- - 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), & - forcex (:,:,iblk), forcey (:,:,iblk), & - umassdti (:,:,iblk), fm (:,:,iblk), & - uarear (:,:,iblk), & - strintx (:,:,iblk), strinty (:,:,iblk), & - taubx (:,:,iblk), tauby (:,:,iblk), & - uvel_init(:,:,iblk), vvel_init(:,:,iblk),& - uvel (:,:,iblk), vvel (:,:,iblk), & - Tbu (:,:,iblk)) + call stepu (nx_block, ny_block, & + icellu (iblk), Cdn_ocn (:,:,iblk), & + indxui (:,iblk), indxuj (:,iblk), & + aiu (:,:,iblk), strtmp (:,:,:), & + uocn (:,:,iblk), vocn (:,:,iblk), & + waterx (:,:,iblk), watery (:,:,iblk), & + forcex (:,:,iblk), forcey (:,:,iblk), & + umassdti (:,:,iblk), fm (:,:,iblk), & + uarear (:,:,iblk), & + strintx (:,:,iblk), strinty (:,:,iblk), & + taubx (:,:,iblk), tauby (:,:,iblk), & + uvel_init(:,:,iblk), vvel_init(:,:,iblk),& + uvel (:,:,iblk), vvel (:,:,iblk), & + Tbu (:,:,iblk)) + + enddo + !$TCXOMP END PARALLEL DO + + call stack_velocity_field(uvel, vvel, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, uvel, vvel) - enddo - !$TCXOMP END PARALLEL DO - - call stack_velocity_field(uvel, vvel, fld2) - call ice_timer_start(timer_bound) - if (maskhalo_dyn) then - call ice_HaloUpdate (fld2, halo_info_mask, & - field_loc_NEcorner, field_type_vector) - else - call ice_HaloUpdate (fld2, halo_info, & - field_loc_NEcorner, field_type_vector) - endif - call ice_timer_stop(timer_bound) - call unstack_velocity_field(fld2, uvel, vvel) - - enddo ! subcycling endif ! evp_algorithm call ice_timer_stop(timer_evp_2d) deallocate(fld2) if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask) - ! Force symmetry across the tripole seam if (trim(grid_type) == 'tripole') then if (maskhalo_dyn) then @@ -574,14 +640,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, & @@ -589,15 +655,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) :: & @@ -616,7 +679,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) :: & @@ -624,12 +686,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 @@ -657,8 +713,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)' !----------------------------------------------------------------- @@ -666,7 +722,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) @@ -872,23 +927,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_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 3a3d0dc6e..a488346f8 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -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 @@ -1388,7 +1386,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 @@ -1401,42 +1399,33 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & ! 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 + 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 - -! endif + zetax2se = (c1+Ktens)*tmpcalcse ! southeast + rep_prsse = (c1-Ktens)*tmpcalcse*Deltase + etax2se = epp2i*zetax2se 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 8e1159ff3..efa47210c 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -1190,14 +1190,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 From 34b074c27929bb808f65835c9c802d16776b663b Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Wed, 22 Dec 2021 11:43:44 +0000 Subject: [PATCH 4/8] clean dyn_finish, correct intent(inout) to intent(in) for Cw, resse Cb in stepu, remove if from seabed_stress_LKD --- cicecore/cicedynB/dynamics/ice_dyn_eap.F90 | 2 - cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 2 - cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 43 +++++++++---------- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 2 - 4 files changed, 20 insertions(+), 29 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index d4b602b70..55022858b 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -541,8 +541,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 3da528e95..419495298 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -616,8 +616,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)) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index a488346f8..4d173c2dc 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -675,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,8 +741,8 @@ subroutine stepu (nx_block, ny_block, & ! calculate seabed stress component for outputs ! only needed on last iteration. - 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) + taubx(i,j) = -uvel(i,j)*Cb + tauby(i,j) = -vvel(i,j)*Cb enddo ! ij end subroutine stepu @@ -760,8 +760,6 @@ subroutine dyn_finish (nx_block, ny_block, & uvel, vvel, & uocn, vocn, & aiu, fm, & - strintx, strinty, & - strairx, strairy, & strocnx, strocny, & strocnxT, strocnyT) @@ -779,11 +777,7 @@ 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) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & strocnx , & ! ice-ocean stress, x-direction @@ -793,14 +787,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)' @@ -891,9 +886,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) :: & @@ -909,18 +905,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 +! if (hwu < threshold_hw) then + docalc_tbu = max(sign(c1,threshold_hw-hwu),c0) - 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 diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index efa47210c..71a4b497b 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -642,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)) From 22893ca645aec0b36f8a8b854eac583b12f72cbf Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Sat, 19 Feb 2022 17:44:26 +0000 Subject: [PATCH 5/8] Reolve conflicts after updating main --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 43 ++++++++++++---------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 2f30c1da7..bb3e5d512 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -454,7 +454,9 @@ subroutine evp (dt) enddo ! subcycling endif ! only subcycle when ndte larger than 1 - !$TCXOMP PARALLEL DO PRIVATE(iblk,strtmp) + + ! ksub = ndte last sub cycling step + !$OMP PARALLEL DO PRIVATE(iblk,strtmp) SCHEDULE(runtime) do iblk = 1, nblocks ! if (trim(yield_curve) == 'ellipse') then @@ -476,9 +478,11 @@ subroutine evp (dt) stress12_3(:,:,iblk), stress12_4(:,:,iblk), & strtmp (:,:,:) ) ! endif ! yield_curve - !----------------------------------------------------------------- - ! on last subcycle, save quantities for mechanical redistribution - !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! on last subcycle, save quantities for mechanical redistribution + !----------------------------------------------------------------- + call deformations (nx_block , ny_block , & icellt(iblk) , & indxti(:,iblk) , indxtj(:,iblk) , & @@ -491,9 +495,9 @@ subroutine evp (dt) rdg_conv(:,:,iblk) , rdg_shear(:,:,iblk) ) - !----------------------------------------------------------------- - ! momentum equation - !----------------------------------------------------------------- + !----------------------------------------------------------------- + ! momentum equation + !----------------------------------------------------------------- call stepu (nx_block, ny_block, & icellu (iblk), Cdn_ocn (:,:,iblk), & @@ -511,21 +515,20 @@ subroutine evp (dt) Tbu (:,:,iblk)) enddo - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO - call stack_velocity_field(uvel, vvel, fld2) - call ice_timer_start(timer_bound) - if (maskhalo_dyn) then - call ice_HaloUpdate (fld2, halo_info_mask, & - field_loc_NEcorner, field_type_vector) - else - call ice_HaloUpdate (fld2, halo_info, & - field_loc_NEcorner, field_type_vector) - endif - call ice_timer_stop(timer_bound) - call unstack_velocity_field(fld2, uvel, vvel) + call stack_velocity_field(uvel, vvel, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, uvel, vvel) - enddo ! subcycling endif ! evp_algorithm call ice_timer_stop(timer_evp_2d) From c833f90cd60ce9dea5a06ed1c6a8cddffe7c74ed Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Mon, 21 Feb 2022 17:39:34 +0000 Subject: [PATCH 6/8] modified environment for Freya to accomodate for additional OMP commands --- configuration/scripts/machines/env.freya_gnu | 3 ++- configuration/scripts/machines/env.freya_intel | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) 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 From 7376807288e776f64e25aa8a87058b978203be2e Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Mon, 21 Feb 2022 20:14:43 +0000 Subject: [PATCH 7/8] 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 --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 20 +++++++++---------- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 18 +++++++++++++---- 2 files changed, 24 insertions(+), 14 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index bb3e5d512..481075381 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -483,16 +483,16 @@ subroutine evp (dt) ! on last subcycle, save quantities for mechanical redistribution !----------------------------------------------------------------- - 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) ) + 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) ) !----------------------------------------------------------------- diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index d9d7a64cc..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 @@ -760,6 +760,8 @@ subroutine dyn_finish (nx_block, ny_block, & uvel, vvel, & uocn, vocn, & aiu, fm, & +! strintx, strinty, & +! strairx, strairy, & strocnx, strocny, & strocnxT, strocnyT) @@ -778,6 +780,10 @@ subroutine dyn_finish (nx_block, ny_block, & 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 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & strocnx , & ! ice-ocean stress, x-direction @@ -905,8 +911,8 @@ 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 = max(sign(c1,threshold_hw-hwu),c0) + 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)) @@ -1395,7 +1401,8 @@ 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 tmpcalcne = capping *(strength/max(Deltane, tinyarea))+ & (c1-capping)* strength/ (Deltane+ tinyarea) tmpcalcnw = capping *(strength/max(Deltanw, tinyarea))+ & @@ -1420,6 +1427,9 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & zetax2se = (c1+Ktens)*tmpcalcse ! southeast rep_prsse = (c1-Ktens)*tmpcalcse*Deltase etax2se = epp2i*zetax2se + ! else + + ! endif end subroutine viscous_coeffs_and_rep_pressure From e312c5ba06611e224a6104ef765b5b8324ef7beb Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Tue, 22 Feb 2022 19:29:36 +0000 Subject: [PATCH 8/8] move deformation to subcycling --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 175 +++++++-------------- 1 file changed, 59 insertions(+), 116 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 481075381..d54a73dd4 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -390,77 +390,18 @@ subroutine evp (dt) divu,rdg_conv,rdg_shear,shear,taubx,tauby ) else ! evp_algorithm == standard_2d (Standard CICE) - if (ndte > 1) then - do ksub = 1,ndte-1 ! subcycling - !----------------------------------------------------------------- - ! stress tensor equation, total surface stress - !----------------------------------------------------------------- + do ksub = 1,ndte ! subcycling - !$OMP PARALLEL DO PRIVATE(iblk,strtmp) SCHEDULE(runtime) - do iblk = 1, nblocks - - 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 (:,:,:) ) - - !----------------------------------------------------------------- - ! momentum equation - !----------------------------------------------------------------- - - call stepu (nx_block, ny_block, & - icellu (iblk), Cdn_ocn (:,:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & - aiu (:,:,iblk), strtmp (:,:,:), & - uocn (:,:,iblk), vocn (:,:,iblk), & - waterx (:,:,iblk), watery (:,:,iblk), & - forcex (:,:,iblk), forcey (:,:,iblk), & - umassdti (:,:,iblk), fm (:,:,iblk), & - uarear (:,:,iblk), & - strintx (:,:,iblk), strinty (:,:,iblk), & - taubx (:,:,iblk), tauby (:,:,iblk), & - uvel_init(:,:,iblk), vvel_init(:,:,iblk),& - uvel (:,:,iblk), vvel (:,:,iblk), & - Tbu (:,:,iblk)) - - enddo - !$OMP END PARALLEL DO - - call stack_velocity_field(uvel, vvel, fld2) - call ice_timer_start(timer_bound) - if (maskhalo_dyn) then - call ice_HaloUpdate (fld2, halo_info_mask, & - field_loc_NEcorner, field_type_vector) - else - call ice_HaloUpdate (fld2, halo_info, & - field_loc_NEcorner, field_type_vector) - endif - call ice_timer_stop(timer_bound) - call unstack_velocity_field(fld2, uvel, vvel) - - enddo ! subcycling - endif ! only subcycle when ndte larger than 1 + !----------------------------------------------------------------- + ! stress tensor equation, total surface stress + !----------------------------------------------------------------- - ! ksub = ndte last sub cycling step - !$OMP PARALLEL DO PRIVATE(iblk,strtmp) SCHEDULE(runtime) - do iblk = 1, nblocks + !$OMP PARALLEL DO PRIVATE(iblk,strtmp) SCHEDULE(runtime) + do iblk = 1, nblocks -! if (trim(yield_curve) == 'ellipse') then - call stress (nx_block, ny_block, & +! if (trim(yield_curve) == 'ellipse') then + call stress (nx_block, ny_block, & icellt(iblk), & indxti (:,iblk), indxtj (:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & @@ -477,64 +418,66 @@ subroutine evp (dt) stress12_1(:,:,iblk), stress12_2(:,:,iblk), & stress12_3(:,:,iblk), stress12_4(:,:,iblk), & strtmp (:,:,:) ) -! endif ! yield_curve - - !----------------------------------------------------------------- - ! on last subcycle, save quantities for mechanical redistribution - !----------------------------------------------------------------- - - 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) ) - - - !----------------------------------------------------------------- - ! momentum equation - !----------------------------------------------------------------- - - call stepu (nx_block, ny_block, & - icellu (iblk), Cdn_ocn (:,:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & - aiu (:,:,iblk), strtmp (:,:,:), & - uocn (:,:,iblk), vocn (:,:,iblk), & - waterx (:,:,iblk), watery (:,:,iblk), & - forcex (:,:,iblk), forcey (:,:,iblk), & - umassdti (:,:,iblk), fm (:,:,iblk), & - uarear (:,:,iblk), & - strintx (:,:,iblk), strinty (:,:,iblk), & - taubx (:,:,iblk), tauby (:,:,iblk), & - uvel_init(:,:,iblk), vvel_init(:,:,iblk),& - uvel (:,:,iblk), vvel (:,:,iblk), & - Tbu (:,:,iblk)) +! 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 + !----------------------------------------------------------------- + + call stepu (nx_block, ny_block, & + icellu (iblk), Cdn_ocn (:,:,iblk), & + indxui (:,iblk), indxuj (:,iblk), & + aiu (:,:,iblk), strtmp (:,:,:), & + uocn (:,:,iblk), vocn (:,:,iblk), & + waterx (:,:,iblk), watery (:,:,iblk), & + forcex (:,:,iblk), forcey (:,:,iblk), & + umassdti (:,:,iblk), fm (:,:,iblk), & + uarear (:,:,iblk), & + strintx (:,:,iblk), strinty (:,:,iblk), & + taubx (:,:,iblk), tauby (:,:,iblk), & + uvel_init(:,:,iblk), vvel_init(:,:,iblk),& + uvel (:,:,iblk), vvel (:,:,iblk), & + Tbu (:,:,iblk)) - enddo - !$OMP END PARALLEL DO + enddo + !$OMP END PARALLEL DO - call stack_velocity_field(uvel, vvel, fld2) - call ice_timer_start(timer_bound) - if (maskhalo_dyn) then - call ice_HaloUpdate (fld2, halo_info_mask, & - field_loc_NEcorner, field_type_vector) - else - call ice_HaloUpdate (fld2, halo_info, & - field_loc_NEcorner, field_type_vector) - endif - call ice_timer_stop(timer_bound) - call unstack_velocity_field(fld2, uvel, vvel) + call stack_velocity_field(uvel, vvel, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, uvel, vvel) + enddo ! subcycling endif ! evp_algorithm call ice_timer_stop(timer_evp_2d) deallocate(fld2) if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask) + ! Force symmetry across the tripole seam if (trim(grid_type) == 'tripole') then if (maskhalo_dyn) then