diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 10a756de6..0c9e1126e 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -98,11 +98,11 @@ subroutine evp (dt) stress12_1, stress12_2, stress12_3, stress12_4, & stresspT, stressmT, stress12T, & stresspU, stressmU, stress12U - use ice_grid, only: tmask, umask, nmask, emask, uvm, epm, npm, & + use ice_grid, only: hm, tmask, umask, umaskCD, nmask, emask, uvm, epm, npm, & dxe, dxn, dxt, dxu, dye, dyn, dyt, dyu, & ratiodxN, ratiodxNr, ratiodyE, ratiodyEr, & dxhy, dyhx, cxp, cyp, cxm, cym, & - tarear, uarear, earea, narea, tinyarea, grid_average_X2Y, & + tarear, uarear, earea, narea, tinyarea, grid_average_X2Y, tarea, & grid_type, grid_system use ice_state, only: aice, vice, vsno, uvel, vvel, uvelN, vvelN, & uvelE, vvelE, divu, shear, & @@ -341,35 +341,67 @@ subroutine evp (dt) ihi = this_block%ihi jlo = this_block%jlo jhi = this_block%jhi - - call dyn_prep2 (nx_block, ny_block, & - ilo, ihi, jlo, jhi, & - icellt(iblk), icellu(iblk), & - indxti (:,iblk), indxtj (:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & - aiu (:,:,iblk), umass (:,:,iblk), & - umassdti (:,:,iblk), fcor_blk (:,:,iblk), & - umask (:,:,iblk), & - uocn (:,:,iblk), vocn (:,:,iblk), & - strairx (:,:,iblk), strairy (:,:,iblk), & - ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), & - icetmask (:,:,iblk), iceumask (:,:,iblk), & - fm (:,:,iblk), dt, & - strtltx (:,:,iblk), strtlty (:,:,iblk), & - strocnx (:,:,iblk), strocny (:,:,iblk), & - strintx (:,:,iblk), strinty (:,:,iblk), & - taubx (:,:,iblk), tauby (:,:,iblk), & - waterx (:,:,iblk), watery (:,:,iblk), & - forcex (:,:,iblk), forcey (:,:,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), & - uvel_init (:,:,iblk), vvel_init (:,:,iblk), & - uvel (:,:,iblk), vvel (:,:,iblk), & - Tbu (:,:,iblk)) + + if (trim(grid_system) == 'B') then + call dyn_prep2 (nx_block, ny_block, & + ilo, ihi, jlo, jhi, & + icellt(iblk), icellu(iblk), & + indxti (:,iblk), indxtj (:,iblk), & + indxui (:,iblk), indxuj (:,iblk), & + aiu (:,:,iblk), umass (:,:,iblk), & + umassdti (:,:,iblk), fcor_blk (:,:,iblk), & + umask (:,:,iblk), & + uocn (:,:,iblk), vocn (:,:,iblk), & + strairx (:,:,iblk), strairy (:,:,iblk), & + ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), & + icetmask (:,:,iblk), iceumask (:,:,iblk), & + fm (:,:,iblk), dt, & + strtltx (:,:,iblk), strtlty (:,:,iblk), & + strocnx (:,:,iblk), strocny (:,:,iblk), & + strintx (:,:,iblk), strinty (:,:,iblk), & + taubx (:,:,iblk), tauby (:,:,iblk), & + waterx (:,:,iblk), watery (:,:,iblk), & + forcex (:,:,iblk), forcey (:,:,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), & + uvel_init (:,:,iblk), vvel_init (:,:,iblk), & + uvel (:,:,iblk), vvel (:,:,iblk), & + Tbu (:,:,iblk)) + + elseif (trim(grid_system) == 'CD') then + call dyn_prep2 (nx_block, ny_block, & + ilo, ihi, jlo, jhi, & + icellt(iblk), icellu(iblk), & + indxti (:,iblk), indxtj (:,iblk), & + indxui (:,iblk), indxuj (:,iblk), & + aiu (:,:,iblk), umass (:,:,iblk), & + umassdti (:,:,iblk), fcor_blk (:,:,iblk), & + umaskCD (:,:,iblk), & + uocn (:,:,iblk), vocn (:,:,iblk), & + strairx (:,:,iblk), strairy (:,:,iblk), & + ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), & + icetmask (:,:,iblk), iceumask (:,:,iblk), & + fm (:,:,iblk), dt, & + strtltx (:,:,iblk), strtlty (:,:,iblk), & + strocnx (:,:,iblk), strocny (:,:,iblk), & + strintx (:,:,iblk), strinty (:,:,iblk), & + taubx (:,:,iblk), tauby (:,:,iblk), & + waterx (:,:,iblk), watery (:,:,iblk), & + forcex (:,:,iblk), forcey (:,:,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), & + uvel_init (:,:,iblk), vvel_init (:,:,iblk), & + uvel (:,:,iblk), vvel (:,:,iblk), & + Tbu (:,:,iblk)) + endif !----------------------------------------------------------------- ! ice strength @@ -412,8 +444,8 @@ subroutine evp (dt) indxti (:,iblk), indxtj (:,iblk), & indxni (:,iblk), indxnj (:,iblk), & aiN (:,:,iblk), nmass (:,:,iblk), & - nmassdti (:,:,iblk), fcorN_blk (:,:,iblk), & - nmask (:,:,iblk), & + nmassdti (:,:,iblk), fcorN_blk (:,:,iblk),& + nmask (:,:,iblk), & uocnN (:,:,iblk), vocnN (:,:,iblk), & strairxN (:,:,iblk), strairyN (:,:,iblk), & ss_tltxN (:,:,iblk), ss_tltyN (:,:,iblk), & @@ -445,7 +477,7 @@ subroutine evp (dt) indxti (:,iblk), indxtj (:,iblk), & indxei (:,iblk), indxej (:,iblk), & aiE (:,:,iblk), emass (:,:,iblk), & - emassdti (:,:,iblk), fcorE_blk (:,:,iblk), & + emassdti (:,:,iblk), fcorE_blk (:,:,iblk),& emask (:,:,iblk), & uocnE (:,:,iblk), vocnE (:,:,iblk), & strairxE (:,:,iblk), strairyE (:,:,iblk), & @@ -643,10 +675,9 @@ subroutine evp (dt) rdg_conv (:,:,iblk), rdg_shear (:,:,iblk), & strtmp (:,:,:) ) - !----------------------------------------------------------------- - ! momentum equation - !----------------------------------------------------------------- - + !----------------------------------------------------------------- + ! momentum equation + !----------------------------------------------------------------- call stepu (nx_block, ny_block, & icellu (iblk), Cdn_ocn (:,:,iblk), & indxui (:,iblk), indxuj (:,iblk), & @@ -688,13 +719,14 @@ subroutine evp (dt) uvel (:,:,iblk), vvel (:,:,iblk), & dxE (:,:,iblk), dyN (:,:,iblk), & dxU (:,:,iblk), dyU (:,:,iblk), & + tarea (:,:,iblk), & ratiodxN (:,:,iblk), ratiodxNr (:,:,iblk), & ratiodyE (:,:,iblk), ratiodyEr (:,:,iblk), & epm (:,:,iblk), npm (:,:,iblk), & - uvm (:,:,iblk), & + hm (:,:,iblk), uvm (:,:,iblk), & zetax2T (:,:,iblk), etax2T (:,:,iblk), & stresspU (:,:,iblk), stressmU (:,:,iblk), & - stress12U (:,:,iblk)) + stress12U (:,:,iblk)) call div_stress (nx_block, ny_block, & ! E point ksub, icelle(iblk), & @@ -977,7 +1009,7 @@ subroutine stress (nx_block, ny_block, & rdg_conv, rdg_shear, & str ) - use ice_dyn_shared, only: strain_rates, deformations, viscous_coeffs_and_rep_pressure + use ice_dyn_shared, only: strain_rates, deformations, viscous_coeffs_and_rep_pressure_T integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -1041,8 +1073,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)' !----------------------------------------------------------------- @@ -1050,7 +1082,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) @@ -1079,17 +1110,27 @@ subroutine stress (nx_block, ny_block, & !----------------------------------------------------------------- ! viscous coefficients and replacement pressure !----------------------------------------------------------------- - - call viscous_coeffs_and_rep_pressure (strength(i,j), tinyarea(i,j),& - Deltane, Deltanw, & - Deltasw, Deltase, & - zetax2ne, zetax2nw, & - zetax2sw, zetax2se, & - etax2ne, etax2nw, & - etax2sw, etax2se, & - rep_prsne, rep_prsnw, & - rep_prssw, rep_prsse, & - capping) + + call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),& + Deltane, zetax2ne, & + etax2ne, rep_prsne, & + capping) + + call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),& + Deltanw, zetax2nw, & + etax2nw, rep_prsnw, & + capping) + + call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),& + Deltasw, zetax2sw, & + etax2sw, rep_prssw, & + capping) + + call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),& + Deltase, zetax2se, & + etax2se, rep_prsse, & + capping) + !----------------------------------------------------------------- ! the stresses ! kg/s^2 @@ -1345,15 +1386,14 @@ subroutine stress_T (nx_block, ny_block, & divT, tensionT, shearT, DeltaT, & ! strain rates at T point rep_prsT ! replacement pressure at T point - logical :: capping ! of the viscous coef - + real(kind=dbl_kind), parameter :: capping = c1 ! of the viscous coef + character(len=*), parameter :: subname = '(stress_T)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- - capping = .true. ! could be later included in ice_in do ij = 1, icellt i = indxti(ij) @@ -1435,15 +1475,16 @@ subroutine stress_U (nx_block, ny_block, & uvelU, vvelU, & dxE, dyN, & dxU, dyU, & + tarea, & ratiodxN, ratiodxNr, & ratiodyE, ratiodyEr, & - epm, npm, uvm, & + epm, npm, hm, uvm, & zetax2T, etax2T, & stresspU, stressmU, & - stress12U ) + stress12U ) - use ice_dyn_shared, only: strain_rates_U!, & - ! viscous_coeffs_and_rep_pressure_U + use ice_dyn_shared, only: strain_rates_U, & + viscous_coeffs_and_rep_pressure_T2U integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -1456,21 +1497,23 @@ subroutine stress_U (nx_block, ny_block, & real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & uvelE , & ! x-component of velocity (m/s) at the E point - vvelE , & ! y-component of velocity (m/s) at the N point - uvelN , & ! x-component of velocity (m/s) at the E point + vvelE , & ! y-component of velocity (m/s) at the E point + uvelN , & ! x-component of velocity (m/s) at the N point vvelN , & ! y-component of velocity (m/s) at the N point uvelU , & ! x-component of velocity (m/s) at the U point vvelU , & ! y-component of velocity (m/s) at the U point - dxE , & ! width of E-cell through the middle (m) + dxE , & ! width of E-cell through the middle (m) dyN , & ! height of N-cell through the middle (m) - dxU , & ! width of U-cell through the middle (m) + dxU , & ! width of U-cell through the middle (m) dyU , & ! height of U-cell through the middle (m) - ratiodxN , & ! -dxN(i+1,j)/dxN(i,j) for BCs - ratiodxNr, & ! -dxN(i,j)/dxN(i+1,j) for BCs - ratiodyE , & ! -dyE(i,j+1)/dyE(i,j) for BCs - ratiodyEr, & ! -dyE(i,j)/dyE(i,j+1) for BCs + tarea , & ! area of T-cell (m^2) + ratiodxN , & ! -dxN(i+1,j)/dxN(i,j) factor for BCs across coastline + ratiodxNr, & ! -dxN(i,j)/dxN(i+1,j) factor for BCs across coastline + ratiodyE , & ! -dyE(i,j+1)/dyE(i,j) factor for BCs across coastline + ratiodyEr, & ! -dyE(i,j)/dyE(i,j+1) factor for BCs across coastline epm , & ! E-cell mask npm , & ! E-cell mask + hm , & ! T-cell mask uvm , & ! U-cell mask zetax2T , & ! 2*zeta at the T point etax2T ! 2*eta at the T point @@ -1517,19 +1560,15 @@ subroutine stress_U (nx_block, ny_block, & ! viscous coefficients and replacement pressure at T point !----------------------------------------------------------------- -! COMING SOON!!! - zetax2U = c0 - etax2U = c0 - rep_prsU = c0 -! call viscous_coeffs_and_rep_pressure_U (zetax2T(i,j), zetax2T(i,j+1), & -! zetax2T(i+1,j+1),zetax2T(i+1,j), & -! etax2T(i,j), etax2T(i,j+1), & -! etax2T(i+1,j+1), etax2T(i+1,j), & -! hm(i,j), hm(i,j+1), & -! hm(i+1,j+1), hm(i+1,j), & -! tarea(i,j), tarea(i,j+1), & -! tarea(i+1,j+1), tarea(i+1,j), & -! DeltaU ) + call viscous_coeffs_and_rep_pressure_T2U (zetax2T(i ,j ), zetax2T(i ,j+1), & + zetax2T(i+1,j+1), zetax2T(i+1,j ), & + etax2T (i ,j ), etax2T (i ,j+1), & + etax2T (i+1,j+1), etax2T (i+1,j ), & + hm (i ,j ), hm (i ,j+1), & + hm (i+1,j+1), hm (i+1,j ), & + tarea (i ,j ), tarea (i ,j+1), & + tarea (i+1,j+1), tarea (i+1,j ), & + DeltaU,zetax2U, etax2U, rep_prsU) !----------------------------------------------------------------- ! the stresses ! kg/s^2 diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index dd2e45237..d4622cb6e 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -31,6 +31,7 @@ module ice_dyn_shared strain_rates, strain_rates_T, strain_rates_U, & viscous_coeffs_and_rep_pressure, & viscous_coeffs_and_rep_pressure_T, & + viscous_coeffs_and_rep_pressure_T2U, & stack_velocity_field, unstack_velocity_field ! namelist parameters @@ -605,13 +606,17 @@ subroutine dyn_prep2 (nx_block, ny_block, & !----------------------------------------------------------------- icellu = 0 + do j = jlo, jhi do i = ilo, ihi - - ! ice extent mask (U-cells) iceumask_old(i,j) = iceumask(i,j) ! save +! if (grid_system == 'B') then ! include ice mask. + ! ice extent mask (U-cells) iceumask(i,j) = (umask(i,j)) .and. (aiu (i,j) > a_min) & - .and. (umass(i,j) > m_min) + .and. (umass(i,j) > m_min) +! else ! ice mask shpuld be applied to cd grid. For now it is not implemented. +! iceumask(i,j) = umask(i,j) +! endif if (iceumask(i,j)) then icellu = icellu + 1 @@ -637,7 +642,6 @@ subroutine dyn_prep2 (nx_block, ny_block, & vvel_init(i,j) = vvel(i,j) enddo enddo - !----------------------------------------------------------------- ! Define variables for momentum equation !----------------------------------------------------------------- @@ -916,17 +920,15 @@ subroutine step_vel (nx_block, ny_block, & ccb = fm(i,j) + sign(c1,fm(i,j)) * vrel * sinw ! kg/m^2 s - ab2 = cca**2 + ccb**2 + ab2 = cca**2 + ccb**2 ! compute the velocity components cc1 = strintx(i,j) + forcex(i,j) + taux & + massdti(i,j)*(brlx*uold + revp*uvel_init(i,j)) cc2 = strinty(i,j) + forcey(i,j) + tauy & + massdti(i,j)*(brlx*vold + revp*vvel_init(i,j)) - 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 .and. seabed_stress) then ! on last subcycling iteration taubx(i,j) = -uvel(i,j)*Tb(i,j) / ccc @@ -1899,7 +1901,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 @@ -1912,39 +1914,30 @@ 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 @@ -1961,6 +1954,7 @@ end subroutine viscous_coeffs_and_rep_pressure ! Lemieux, J. F. et al. (2016). Improving the simulation of landfast ice ! by combining tensile strength and a parameterization for grounded ridges. ! J. Geophys. Res. Oceans, 121, 7354-7368. +! capping must be 1 (c1) for evp and 0 for vp solver subroutine viscous_coeffs_and_rep_pressure_T (strength, tinyarea, & Delta , zetax2 , & @@ -1971,9 +1965,7 @@ subroutine viscous_coeffs_and_rep_pressure_T (strength, tinyarea, & strength, tinyarea real (kind=dbl_kind), intent(in):: & - Delta - - logical, intent(in):: capping + Delta, capping real (kind=dbl_kind), intent(out):: & zetax2, etax2, rep_prs ! 2 x visous coeffs, replacement pressure @@ -1986,23 +1978,62 @@ subroutine viscous_coeffs_and_rep_pressure_T (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 - tmpcalc = strength/max(Delta,tinyarea) - else - tmpcalc = strength/(Delta + tinyarea) - endif - + tmpcalc = capping *(strength/max(Delta,tinyarea))+ & + (c1-capping)*(strength/(Delta + tinyarea)) zetax2 = (c1+Ktens)*tmpcalc rep_prs = (c1-Ktens)*tmpcalc*Delta etax2 = epp2i*zetax2 -! else + end subroutine viscous_coeffs_and_rep_pressure_T -! endif - end subroutine viscous_coeffs_and_rep_pressure_T + subroutine viscous_coeffs_and_rep_pressure_T2U (zetax2T_00, zetax2T_01, & + zetax2T_11, zetax2T_10, & + etax2T_00, etax2T_01, & + etax2T_11, etax2T_10, & + maskT_00, maskT_01, & + maskT_11, maskT_10, & + tarea_00, tarea_01, & + tarea_11, tarea_10, & + deltaU, & + zetax2U, etax2U, & + rep_prsU) + + + real (kind=dbl_kind), intent(in):: & + zetax2T_00,zetax2T_10,zetax2T_11,zetax2T_01, & + etax2T_00, etax2T_10, etax2T_11, etax2T_01, & ! 2 x visous coeffs, replacement pressure + maskT_00, maskT_10, maskT_11, maskT_01, & + tarea_00, tarea_10, tarea_11, tarea_01, & + deltaU + + real (kind=dbl_kind), intent(out):: zetax2U, etax2U, rep_prsU + + ! local variables + + real (kind=dbl_kind) :: & + Totarea + + character(len=*), parameter :: subname = '(viscous_coeffs_and_rep_pressure_T2U)' + + ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code + Totarea = maskT_00*Tarea_00 + & + maskT_10*Tarea_10 + & + maskT_11*Tarea_11 + & + maskT_01*Tarea_01 + zetax2U = (maskT_00*Tarea_00 *zetax2T_00 + & + maskT_10*Tarea_10 *zetax2T_10 + & + maskT_11*Tarea_11 *zetax2T_11 + & + maskT_01*Tarea_01 *zetax2T_01)/Totarea + + etax2U = (maskT_00*Tarea_00 *etax2T_00 + & + maskT_10*Tarea_10 *etax2T_10 + & + maskT_11*Tarea_11 *etax2T_11 + & + maskT_01*Tarea_01 *etax2T_01)/Totarea + + rep_prsU = (c1-ktens)/(1+Ktens)*zetax2U*deltaU + + end subroutine viscous_coeffs_and_rep_pressure_T2U !======================================================================= diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 2c1b23032..5bbef02a4 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -1202,13 +1202,11 @@ 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) diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index 8212def06..9dfc7b27b 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -155,7 +155,7 @@ module ice_grid real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & hm , & ! land/boundary mask, thickness (T-cell) bm , & ! task/block id - uvm , & ! land/boundary mask, velocity (U-cell) + uvm , & ! land/boundary mask, velocity (U-cell) npm , & ! land/boundary mask (N-cell) epm , & ! land/boundary mask (E-cell) kmt ! ocean topography mask for bathymetry (T-cell) @@ -167,7 +167,8 @@ module ice_grid logical (kind=log_kind), & dimension (:,:,:), allocatable, public :: & tmask , & ! land/boundary mask, thickness (T-cell) - umask , & ! land/boundary mask, velocity (U-cell) + umask , & ! land/boundary mask, velocity (U-cell) (1 if all surrounding T cells are ocean) + umaskCD, & ! land/boundary mask, velocity (U-cell) (1 if at least two surrounding T cells are ocean) nmask , & ! land/boundary mask, (N-cell) emask , & ! land/boundary mask, (E-cell) lmask_n, & ! northern hemisphere mask @@ -239,12 +240,13 @@ subroutine alloc_grid yyav (nx_block,ny_block,max_blocks), & ! mean T-cell value of yy hm (nx_block,ny_block,max_blocks), & ! land/boundary mask, thickness (T-cell) bm (nx_block,ny_block,max_blocks), & ! task/block id - uvm (nx_block,ny_block,max_blocks), & ! land/boundary mask, velocity (U-cell) + uvm (nx_block,ny_block,max_blocks), & ! land/boundary mask, velocity (U-cell) - water in case of all water point npm (nx_block,ny_block,max_blocks), & ! land/boundary mask (N-cell) epm (nx_block,ny_block,max_blocks), & ! land/boundary mask (E-cell) kmt (nx_block,ny_block,max_blocks), & ! ocean topography mask for bathymetry (T-cell) tmask (nx_block,ny_block,max_blocks), & ! land/boundary mask, thickness (T-cell) umask (nx_block,ny_block,max_blocks), & ! land/boundary mask, velocity (U-cell) + umaskCD (nx_block,ny_block,max_blocks), & ! land/boundary mask, velocity (U-cell) nmask (nx_block,ny_block,max_blocks), & ! land/boundary mask (N-cell) emask (nx_block,ny_block,max_blocks), & ! land/boundary mask (E-cell) lmask_n (nx_block,ny_block,max_blocks), & ! northern hemisphere mask @@ -1925,7 +1927,7 @@ end subroutine primary_grid_lengths_HTE subroutine makemask - use ice_constants, only: c0, p5, & + use ice_constants, only: c0, p5, c1p5, & field_loc_center, field_loc_NEcorner, field_type_scalar, & field_loc_Nface, field_loc_Eface @@ -1936,6 +1938,9 @@ subroutine makemask real (kind=dbl_kind) :: & puny + real (kind=dbl_kind), dimension(:,:,:), allocatable :: & + uvmCD + type (block) :: & this_block ! block information for current block @@ -1958,6 +1963,7 @@ subroutine makemask !----------------------------------------------------------------- bm = c0 + allocate(uvmCD(nx_block,ny_block,max_blocks)) !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) do iblk = 1, nblocks @@ -1974,6 +1980,8 @@ subroutine makemask npm(i,j,iblk) = min (hm(i,j, iblk), hm(i,j+1,iblk)) epm(i,j,iblk) = min (hm(i,j, iblk), hm(i+1,j,iblk)) bm(i,j,iblk) = my_task + iblk/100.0_dbl_kind + uvmCD(i,j,iblk) = (hm(i,j, iblk)+hm(i+1,j, iblk) & + + hm(i,j+1,iblk)+hm(i+1,j+1,iblk)) enddo enddo enddo @@ -1982,6 +1990,8 @@ subroutine makemask call ice_timer_start(timer_bound) call ice_HaloUpdate (uvm, halo_info, & field_loc_NEcorner, field_type_scalar) + call ice_HaloUpdate (uvmCD, halo_info, & + field_loc_NEcorner, field_type_scalar) call ice_HaloUpdate (npm, halo_info, & field_loc_Nface, field_type_scalar) call ice_HaloUpdate (epm, halo_info, & @@ -1999,19 +2009,32 @@ subroutine makemask jhi = this_block%jhi ! needs to cover halo (no halo update for logicals) - tmask(:,:,iblk) = .false. - umask(:,:,iblk) = .false. - nmask(:,:,iblk) = .false. - emask(:,:,iblk) = .false. + tmask(:,:,iblk) = .false. + umask(:,:,iblk) = .false. + umaskCD(:,:,iblk) = .false. + nmask(:,:,iblk) = .false. + emask(:,:,iblk) = .false. do j = jlo-nghost, jhi+nghost do i = ilo-nghost, ihi+nghost - if ( hm(i,j,iblk) > p5) tmask(i,j,iblk) = .true. - if (uvm(i,j,iblk) > p5) umask(i,j,iblk) = .true. - if (npm(i,j,iblk) > p5) nmask(i,j,iblk) = .true. - if (epm(i,j,iblk) > p5) emask(i,j,iblk) = .true. + if ( hm(i,j,iblk) > p5 ) tmask (i,j,iblk) = .true. + if (uvm(i,j,iblk) > p5 ) umask (i,j,iblk) = .true. + if (uvmCD(i,j,iblk) > c1p5) umaskCD(i,j,iblk) = .true. + if (npm(i,j,iblk) > p5 ) nmask (i,j,iblk) = .true. + if (epm(i,j,iblk) > p5 ) emask (i,j,iblk) = .true. enddo enddo + enddo + !$OMP END PARALLEL DO + + !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + !----------------------------------------------------------------- ! create hemisphere masks !----------------------------------------------------------------- @@ -2045,6 +2068,8 @@ subroutine makemask enddo ! iblk !$OMP END PARALLEL DO + deallocate(uvmCD) + end subroutine makemask !=======================================================================