From cf65c3987106a878664499f8d97debe55e91854a Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Wed, 17 Nov 2021 19:29:06 +0000 Subject: [PATCH 01/13] Changes in ice_dyn* are interpolation of uvelE/vvelN to B grid. ice_transport files are changed in a way so that velocities are interpolated to b grid for depature point function and kept at E or N grid possible. --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 5 ++--- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 1 + 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index fc6d8eed0..b7d3e5b51 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -914,10 +914,9 @@ subroutine evp (dt) field_loc_NEcorner, field_type_vector) call grid_average_X2Y('U2TF',strocnxT) ! shift call grid_average_X2Y('U2TF',strocnyT) -! shift velocity components from CD grid locations (N, E) to B grid location (U) for transport if (grid_system == 'CD') then - call grid_average_X2Y('E2US',uvelE,uvel) - call grid_average_X2Y('N2US',vvelN,vvel) + call grid_average_X2Y('E2UF',uvelE,uvel) + call grid_average_X2Y('N2UF',vvelN,vvel) endif call ice_timer_stop(timer_dynamics) ! dynamics diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 2c1b23032..03010d590 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -657,6 +657,7 @@ subroutine implicit_solver (dt) field_loc_NEcorner, field_type_vector) call grid_average_X2Y('U2TF',strocnxT) ! shift call grid_average_X2Y('U2TF',strocnyT) +<<<<<<< HEAD ! shift velocity components from CD grid locations (N, E) to B grid location (U) for transport ! commented out in order to focus on EVP for now within the cdgrid ! should be used when routine is ready From 89e401a95a79c49ab459d63c40b26967b48fa1ab Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Wed, 17 Nov 2021 20:47:36 +0000 Subject: [PATCH 02/13] changed according to comments. changed average from F to S and. commented out in vp and eap --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 5 +++-- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 1 - 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index b7d3e5b51..fc6d8eed0 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -914,9 +914,10 @@ subroutine evp (dt) field_loc_NEcorner, field_type_vector) call grid_average_X2Y('U2TF',strocnxT) ! shift call grid_average_X2Y('U2TF',strocnyT) +! shift velocity components from CD grid locations (N, E) to B grid location (U) for transport if (grid_system == 'CD') then - call grid_average_X2Y('E2UF',uvelE,uvel) - call grid_average_X2Y('N2UF',vvelN,vvel) + call grid_average_X2Y('E2US',uvelE,uvel) + call grid_average_X2Y('N2US',vvelN,vvel) endif call ice_timer_stop(timer_dynamics) ! dynamics diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 03010d590..2c1b23032 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -657,7 +657,6 @@ subroutine implicit_solver (dt) field_loc_NEcorner, field_type_vector) call grid_average_X2Y('U2TF',strocnxT) ! shift call grid_average_X2Y('U2TF',strocnyT) -<<<<<<< HEAD ! shift velocity components from CD grid locations (N, E) to B grid location (U) for transport ! commented out in order to focus on EVP for now within the cdgrid ! should be used when routine is ready From 24d98dbe8ac207ae078f63b6a8601879a06ff8ee Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Wed, 17 Nov 2021 19:29:06 +0000 Subject: [PATCH 03/13] Changes in ice_dyn* are interpolation of uvelE/vvelN to B grid. ice_transport files are changed in a way so that velocities are interpolated to b grid for depature point function and kept at E or N grid possible. --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index fc6d8eed0..b7d3e5b51 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -914,10 +914,9 @@ subroutine evp (dt) field_loc_NEcorner, field_type_vector) call grid_average_X2Y('U2TF',strocnxT) ! shift call grid_average_X2Y('U2TF',strocnyT) -! shift velocity components from CD grid locations (N, E) to B grid location (U) for transport if (grid_system == 'CD') then - call grid_average_X2Y('E2US',uvelE,uvel) - call grid_average_X2Y('N2US',vvelN,vvel) + call grid_average_X2Y('E2UF',uvelE,uvel) + call grid_average_X2Y('N2UF',vvelN,vvel) endif call ice_timer_stop(timer_dynamics) ! dynamics From 8849523757d2bca631a7df73bccfd29753e99247 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Wed, 17 Nov 2021 20:47:36 +0000 Subject: [PATCH 04/13] changed according to comments. changed average from F to S and. commented out in vp and eap --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index b7d3e5b51..fc6d8eed0 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -914,9 +914,10 @@ subroutine evp (dt) field_loc_NEcorner, field_type_vector) call grid_average_X2Y('U2TF',strocnxT) ! shift call grid_average_X2Y('U2TF',strocnyT) +! shift velocity components from CD grid locations (N, E) to B grid location (U) for transport if (grid_system == 'CD') then - call grid_average_X2Y('E2UF',uvelE,uvel) - call grid_average_X2Y('N2UF',vvelN,vvel) + call grid_average_X2Y('E2US',uvelE,uvel) + call grid_average_X2Y('N2US',vvelN,vvel) endif call ice_timer_stop(timer_dynamics) ! dynamics From ae409640020e5ea605ef7f27f4229ef9036dd708 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Fri, 19 Nov 2021 19:07:09 +0000 Subject: [PATCH 05/13] Create local uvmCD that is projected into umaskCD and icellu. Added function to fint zeta, eta and replacement pressure at U --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 44 +++++-- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 122 ++++++++++++++---- cicecore/cicedynB/infrastructure/ice_grid.F90 | 25 ++-- 3 files changed, 144 insertions(+), 47 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index fc6d8eed0..65d240adb 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -1015,7 +1015,7 @@ 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)' @@ -1024,7 +1024,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) @@ -1054,16 +1053,37 @@ 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 (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 diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 3a3f4b154..9e43b14c6 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -605,6 +605,7 @@ subroutine dyn_prep2 (nx_block, ny_block, & !----------------------------------------------------------------- icellu = 0 + if (grid_system == 'B') then do j = jlo, jhi do i = ilo, ihi @@ -637,7 +638,40 @@ subroutine dyn_prep2 (nx_block, ny_block, & vvel_init(i,j) = vvel(i,j) enddo enddo + else + do j = jlo, jhi + do i = ilo, ihi + + ! ice extent mask (U-cells) +! iceumask_old(i,j) = iceumask(i,j) ! save +! iceumask(i,j) = (umaskCD(i,j)) .and. (aiu (i,j) > a_min) & +! .and. (umass(i,j) > m_min) + iceumask(i,j) = umaskCD(i,j) ! SHOULD BE EXTENDED TO CHECK FOR ICE + if (iceumask(i,j)) then + icellu = icellu + 1 + indxui(icellu) = i + indxuj(icellu) = j + + ! initialize velocity for new ice points to ocean sfc current + if (.not. iceumask_old(i,j)) then + uvel(i,j) = uocn(i,j) + vvel(i,j) = vocn(i,j) + endif + else + ! set velocity and stresses to zero for masked-out points + uvel(i,j) = c0 + vvel(i,j) = c0 + strintx(i,j) = c0 + strinty(i,j) = c0 + strocnx(i,j) = c0 + strocny(i,j) = c0 + endif + uvel_init(i,j) = uvel(i,j) + vvel_init(i,j) = vvel(i,j) + enddo + enddo + endif !----------------------------------------------------------------- ! Define variables for momentum equation !----------------------------------------------------------------- @@ -1899,7 +1933,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 @@ -1913,18 +1947,14 @@ 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 + tmpcalcne = capping*(strength/max(Deltane, tinyarea))+ & + (1-capping)*strength/(Deltane + tinyarea) & + tmpcalcnw = capping*(strength/max(Deltanw, tinyarea))+ & + (1-capping)*strength/(Deltanw + tinyarea) & + tmpcalcsw = capping*(strength/max(Deltasw, tinyarea))+ & + (1-capping)*strength/(Deltasw + tinyarea) + tmpcalcse = capping*(strength/max(Deltase,tinyarea)) + & + (1-capping)*strength/(Deltase + tinyarea) zetax2ne = (c1+Ktens)*tmpcalcne ! northeast rep_prsne = (c1-Ktens)*tmpcalcne*Deltane @@ -1961,6 +1991,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 +2002,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 +2015,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, maskeT_01, & + maskT_11, maskeT_10, & + Tarea_00, Tarea_01, & + Tarea_11, Tarea_10, & + deltaU, & + zetax2U, etax2U, & + rep_prsU + + + real (kind=dbl_kind), intent(in):: & + zetax2_00,zetax2_10,zetax2_11,zetax2_01, & + etax2_00, etax2_10, etax2_11, etax2_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(i,j )*Tarea(i,j) + & + maskT(i+1,j)*Tarea(i+1,j) + & + maskT(i+1,j+1)*Tarea(i+1,j+1) + & + maskT(i,j+1)*Tarea(i,j+1) + + zetax2U = (maskT(i,j)*Tarea(i,j) *zetax2T(i,j) + & + maskT(i+1,j)*Tarea(i+1,j) *zetax2T(i+1,j) + & + maskT(i,j+1)*Tarea(i,j+1) *zetax2T(i,j+1) + & + maskT(i+1,j+1)*Tarea(i+1,j+1)*zetax2T(i+1,j+1))/Totarea + + etax2U = (maskT(i,j) * etax2T(i,j) + & + maskT(i+1,j) * etax2T(i+1,j) + & + maskT(i,j+1) * etax2T(i,j+1) + & + maskT(i+1,j+1)* etax2T(i+1,j+1))/Totarea + + rep_prsU = (c1-ktens)/(1+Ktens)*zetax2U*deltaU + + end subroutine viscous_coeffs_and_rep_pressure_T2U !======================================================================= diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index 688b7155e..6bd828fef 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -153,7 +153,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) - land if one point is one npm , & ! land/boundary mask (N-cell) epm , & ! land/boundary mask (E-cell) kmt ! ocean topography mask for bathymetry (T-cell) @@ -165,7 +165,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) (calc point if all point is ocean) + umaskCD, & ! land/boundary mask, velocity (U-cell) (calc point if one point is ocean) nmask , & ! land/boundary mask, (N-cell) emask , & ! land/boundary mask, (E-cell) lmask_n, & ! northern hemisphere mask @@ -235,12 +236,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 @@ -1911,7 +1913,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 @@ -1922,6 +1924,9 @@ subroutine makemask real (kind=dbl_kind) :: & puny + real (kind=dbl_kind) :: & + uvmCD + type (block) :: & this_block ! block information for current block @@ -1985,14 +1990,18 @@ 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. + uvmCD = (hm(i,j, iblk)+hm(i+1,j, iblk)+ + hm(i,j+1,iblk)+hm(i+1,j+1,iblk)) + if (uvmCD > 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 From 3855b8c017079404c04946e0cc03fd2b748d219d Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Fri, 19 Nov 2021 19:28:19 +0000 Subject: [PATCH 06/13] Added call to routine --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 19 +++++++++---------- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 1 + 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 65d240adb..3b8e21688 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -1511,17 +1511,16 @@ subroutine stress_U (nx_block, ny_block, & ! viscous coefficients and replacement pressure at T point !----------------------------------------------------------------- -! COMING SOON!!! -! 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 9e43b14c6..bc238a589 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 From 004aa9cbd801dcdd7033e7db0faae808e31adb15 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Fri, 19 Nov 2021 20:58:42 +0000 Subject: [PATCH 07/13] bux fixes --- cicecore/cicedynB/dynamics/ice_dyn_eap.F90 | 4 +- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 26 +++++----- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 47 ++++++++++--------- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 8 ++-- cicecore/cicedynB/infrastructure/ice_grid.F90 | 4 +- 5 files changed, 44 insertions(+), 45 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index e5a89b118..b0b47c5cb 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -133,7 +133,7 @@ subroutine eap (dt) stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 - use ice_grid, only: tmask, umask, dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, & + use ice_grid, only: tmask, umask, umaskCD, dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, & tarear, uarear, grid_average_X2Y!, grid_system commented out until implementation of cd-grid use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, & aice_init, aice0, aicen, vicen, strength !, uvelE, vvelN grid_system commented out until implementation of cd-grid @@ -301,7 +301,7 @@ subroutine eap (dt) indxui (:,iblk), indxuj (:,iblk), & aiu (:,:,iblk), umass (:,:,iblk), & umassdti (:,:,iblk), fcor_blk (:,:,iblk), & - umask (:,:,iblk), & + umask (:,:,iblk), umaskCD (:,:,iblk), & uocn (:,:,iblk), vocn (:,:,iblk), & strairx (:,:,iblk), strairy (:,:,iblk), & ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), & diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 3b8e21688..2de428784 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, tinyarea, grid_average_X2Y, & + tarear, uarear, tinyarea, grid_average_X2Y, tarea, & grid_type, grid_system use ice_state, only: aice, vice, vsno, uvel, vvel, uvelN, vvelN, & uvelE, vvelE, divu, shear, & @@ -349,7 +349,7 @@ subroutine evp (dt) indxui (:,iblk), indxuj (:,iblk), & aiu (:,:,iblk), umass (:,:,iblk), & umassdti (:,:,iblk), fcor_blk (:,:,iblk), & - umask (:,:,iblk), & + umask (:,:,iblk), umaskCD (:,:,iblk), & uocn (:,:,iblk), vocn (:,:,iblk), & strairx (:,:,iblk), strairy (:,:,iblk), & ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), & @@ -412,8 +412,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), nmask (:,:,iblk), & uocnN (:,:,iblk), vocnN (:,:,iblk), & strairxN (:,:,iblk), strairyN (:,:,iblk), & ss_tltxN (:,:,iblk), ss_tltyN (:,:,iblk), & @@ -445,8 +445,8 @@ subroutine evp (dt) indxti (:,iblk), indxtj (:,iblk), & indxei (:,iblk), indxej (:,iblk), & aiE (:,:,iblk), emass (:,:,iblk), & - emassdti (:,:,iblk), fcorE_blk (:,:,iblk), & - emask (:,:,iblk), & + emassdti (:,:,iblk), fcorE_blk (:,:,iblk),& + emask (:,:,iblk), emask (:,:,iblk),& uocnE (:,:,iblk), vocnE (:,:,iblk), & strairxE (:,:,iblk), strairyE (:,:,iblk), & ss_tltxE (:,:,iblk), ss_tltyE (:,:,iblk), & @@ -691,7 +691,7 @@ subroutine evp (dt) 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)) @@ -1015,7 +1015,7 @@ subroutine stress (nx_block, ny_block, & str12ew, str12we, str12ns, str12sn , & strp_tmp, strm_tmp, tmp - real(kind=dbl_kind,parameter :: capping = c1 ! of the viscous coef + real(kind=dbl_kind),parameter :: capping = c1 ! of the viscous coef character(len=*), parameter :: subname = '(stress)' @@ -1339,15 +1339,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) @@ -1431,7 +1430,7 @@ subroutine stress_U (nx_block, ny_block, & dxU, dyU, & ratiodxN, ratiodxNr, & ratiodyE, ratiodyEr, & - epm, npm, uvm, & + epm, npm, hm, uvm, & zetax2T, etax2T, & stresspU, stressmU, & stress12U ) @@ -1465,6 +1464,7 @@ subroutine stress_U (nx_block, ny_block, & ratiodyEr, & ! -dyE(i,j)/dyE(i,j+1) for BCs 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 diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index bc238a589..f145037f0 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -453,7 +453,7 @@ subroutine dyn_prep2 (nx_block, ny_block, & indxui, indxuj, & aiu, umass, & umassdti, fcor, & - umask, & + umask, umaskCD, & uocn, vocn, & strairx, strairy, & ss_tltx, ss_tlty, & @@ -490,7 +490,7 @@ subroutine dyn_prep2 (nx_block, ny_block, & indxuj ! compressed index in j-direction logical (kind=log_kind), dimension (nx_block,ny_block), intent(in) :: & - umask ! land/boundary mask, thickness (U-cell) + umask, umaskCD ! land/boundary mask, thickness (U-cell) integer (kind=int_kind), dimension (nx_block,ny_block), intent(in) :: & icetmask ! ice extent mask (T-cell) @@ -1949,9 +1949,9 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & ! if (trim(yield_curve) == 'ellipse') then tmpcalcne = capping*(strength/max(Deltane, tinyarea))+ & - (1-capping)*strength/(Deltane + tinyarea) & + (1-capping)*strength/(Deltane + tinyarea) tmpcalcnw = capping*(strength/max(Deltanw, tinyarea))+ & - (1-capping)*strength/(Deltanw + tinyarea) & + (1-capping)*strength/(Deltanw + tinyarea) tmpcalcsw = capping*(strength/max(Deltasw, tinyarea))+ & (1-capping)*strength/(Deltasw + tinyarea) tmpcalcse = capping*(strength/max(Deltase,tinyarea)) + & @@ -2029,18 +2029,18 @@ 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, maskeT_01, & - maskT_11, maskeT_10, & + maskT_00, maskT_01, & + maskT_11, maskT_10, & Tarea_00, Tarea_01, & Tarea_11, Tarea_10, & deltaU, & zetax2U, etax2U, & - rep_prsU + rep_prsU) real (kind=dbl_kind), intent(in):: & - zetax2_00,zetax2_10,zetax2_11,zetax2_01, & - etax2_00, etax2_10, etax2_11, etax2_01, & ! 2 x visous coeffs, replacement pressure + 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 @@ -2049,25 +2049,26 @@ subroutine viscous_coeffs_and_rep_pressure_T2U (zetax2T_00, zetax2T_01, & ! local variables - real (kind=dbl_kind), :: Totarea + 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(i,j )*Tarea(i,j) + & - maskT(i+1,j)*Tarea(i+1,j) + & - maskT(i+1,j+1)*Tarea(i+1,j+1) + & - maskT(i,j+1)*Tarea(i,j+1) + Totarea = maskT_00*Tarea_00 + & + maskT_10*Tarea_10 + & + maskT_11*Tarea_11 + & + maskT_01*Tarea_01 - zetax2U = (maskT(i,j)*Tarea(i,j) *zetax2T(i,j) + & - maskT(i+1,j)*Tarea(i+1,j) *zetax2T(i+1,j) + & - maskT(i,j+1)*Tarea(i,j+1) *zetax2T(i,j+1) + & - maskT(i+1,j+1)*Tarea(i+1,j+1)*zetax2T(i+1,j+1))/Totarea - - etax2U = (maskT(i,j) * etax2T(i,j) + & - maskT(i+1,j) * etax2T(i+1,j) + & - maskT(i,j+1) * etax2T(i,j+1) + & - maskT(i+1,j+1)* etax2T(i+1,j+1))/Totarea + 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 diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 2c1b23032..bbbcb67c5 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -199,7 +199,7 @@ subroutine implicit_solver (dt) stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 - use ice_grid, only: tmask, umask, dxt, dyt, cxp, cyp, cxm, cym, & + use ice_grid, only: tmask, umask, umaskCD, dxt, dyt, cxp, cyp, cxm, cym, & tarear, grid_type, grid_average_X2Y !, grid_system commented out until implementation of c grid use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, & aice_init, aice0, aicen, vicen, strength!, uvelE, vvelN ommented out until implementation of c grid @@ -366,7 +366,7 @@ subroutine implicit_solver (dt) indxui (:,iblk), indxuj (:,iblk), & aiu (:,:,iblk), umass (:,:,iblk), & umassdti (:,:,iblk), fcor_blk (:,:,iblk), & - umask (:,:,iblk), & + umask (:,:,iblk), umaskCD (:,:,iblk), & uocn (:,:,iblk), vocn (:,:,iblk), & strairx (:,:,iblk), strairy (:,:,iblk), & ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), & @@ -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 6bd828fef..a69e16f46 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -1913,7 +1913,7 @@ end subroutine primary_grid_lengths_HTE subroutine makemask - use ice_constants, only: c0, p5, c1p5& + use ice_constants, only: c0, p5, c1p5, & field_loc_center, field_loc_NEcorner, field_type_scalar, & field_loc_Nface, field_loc_Eface @@ -1999,7 +1999,7 @@ subroutine makemask 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. - uvmCD = (hm(i,j, iblk)+hm(i+1,j, iblk)+ + uvmCD = (hm(i,j, iblk)+hm(i+1,j, iblk)+ & hm(i,j+1,iblk)+hm(i+1,j+1,iblk)) if (uvmCD > c1p5) umaskCD(i,j,iblk) = .true. if (npm(i,j,iblk) > p5) nmask(i,j,iblk) = .true. From fcd291b95b3f20158b83aac13efa6b1b81dbfae9 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Sat, 20 Nov 2021 15:36:04 +0000 Subject: [PATCH 08/13] last bugfixes. Code compiles and is bit for bit with code before. Only quick suite has been run. --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 30 ++++++++++--------- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 8 ++--- cicecore/cicedynB/infrastructure/ice_grid.F90 | 13 ++++++-- 3 files changed, 30 insertions(+), 21 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 2de428784..7c305547b 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -694,7 +694,7 @@ subroutine evp (dt) hm (:,:,iblk), uvm (:,:,iblk), & zetax2T (:,:,iblk), etax2T (:,:,iblk), & stresspU (:,:,iblk), stressmU (:,:,iblk), & - stress12U (:,:,iblk)) + stress12U (:,:,iblk), tarea (:,:,iblk)) call step_vel (nx_block, ny_block, & ! E point icelle (iblk), Cdn_ocn (:,:,iblk), & @@ -951,7 +951,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, viscous_coeffs_and_rep_pressure_T integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -1433,10 +1433,10 @@ subroutine stress_U (nx_block, ny_block, & epm, npm, hm, uvm, & zetax2T, etax2T, & stresspU, stressmU, & - stress12U ) + stress12U, Tarea ) - 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 @@ -1465,6 +1465,7 @@ subroutine stress_U (nx_block, ny_block, & epm , & ! E-cell mask npm , & ! E-cell mask hm , & ! T-cell mask + Tarea , & ! area of T-cell uvm , & ! U-cell mask zetax2T , & ! 2*zeta at the T point etax2T ! 2*eta at the T point @@ -1512,15 +1513,16 @@ subroutine stress_U (nx_block, ny_block, & !----------------------------------------------------------------- - 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) + 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 f145037f0..1d95b2ba8 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -1949,13 +1949,13 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & ! if (trim(yield_curve) == 'ellipse') then tmpcalcne = capping*(strength/max(Deltane, tinyarea))+ & - (1-capping)*strength/(Deltane + tinyarea) + (c1-capping)*strength/(Deltane + tinyarea) tmpcalcnw = capping*(strength/max(Deltanw, tinyarea))+ & - (1-capping)*strength/(Deltanw + tinyarea) + (c1-capping)*strength/(Deltanw + tinyarea) tmpcalcsw = capping*(strength/max(Deltasw, tinyarea))+ & - (1-capping)*strength/(Deltasw + tinyarea) + (c1-capping)*strength/(Deltasw + tinyarea) tmpcalcse = capping*(strength/max(Deltase,tinyarea)) + & - (1-capping)*strength/(Deltase + tinyarea) + (c1-capping)*strength/(Deltase + tinyarea) zetax2ne = (c1+Ktens)*tmpcalcne ! northeast rep_prsne = (c1-Ktens)*tmpcalcne*Deltane diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index a69e16f46..fd7717b88 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -1995,13 +1995,20 @@ subroutine makemask umaskCD(:,:,iblk) = .false. nmask(:,:,iblk) = .false. emask(:,:,iblk) = .false. + do j = jlo-nghost, jhi + do i = ilo-nghost, ihi +!!!!!WARNING THIS SUM is not calculating uvmCD mask on the northern and eastern halo +! Loop had to be separate as plus 1 indexes exceed array size + uvmCD = (hm(i,j, iblk)+hm(i+1,j, iblk)+ & + hm(i,j+1,iblk)+hm(i+1,j+1,iblk)) + if (uvmCD > c1p5) umaskCD(i,j,iblk) = .true. + enddo + enddo + 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. - uvmCD = (hm(i,j, iblk)+hm(i+1,j, iblk)+ & - hm(i,j+1,iblk)+hm(i+1,j+1,iblk)) - if (uvmCD > 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 From 710967f9d4f36619de51de25db2d3c311e15031d Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Sun, 21 Nov 2021 19:46:29 +0000 Subject: [PATCH 09/13] First move towards merge with merge of commit 39 --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 57 ++++++--------- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 57 +++++++-------- cicecore/cicedynB/infrastructure/ice_grid.F90 | 71 ++++++++++++------- .../infrastructure/ice_restart_driver.F90 | 4 ++ .../scripts/options/set_nml.boxnodyn | 31 +++++++- 5 files changed, 128 insertions(+), 92 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 7c305547b..2c4df2ac2 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -102,7 +102,7 @@ subroutine evp (dt) dxe, dxn, dxt, dxu, dye, dyn, dyt, dyu, & ratiodxN, ratiodxNr, ratiodyE, ratiodyEr, & dxhy, dyhx, cxp, cyp, cxm, cym, & - tarear, uarear, tinyarea, grid_average_X2Y, tarea, & + 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, & @@ -643,9 +643,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), & @@ -1016,7 +1016,7 @@ subroutine stress (nx_block, ny_block, & strp_tmp, strm_tmp, tmp real(kind=dbl_kind),parameter :: capping = c1 ! of the viscous coef - + character(len=*), parameter :: subname = '(stress)' !----------------------------------------------------------------- @@ -1052,17 +1052,6 @@ 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, & @@ -1339,8 +1328,8 @@ subroutine stress_T (nx_block, ny_block, & divT, tensionT, shearT, DeltaT, & ! strain rates at T point rep_prsT ! replacement pressure at T point - real(kind=dbl_kind),parameter :: capping = c1 ! of the viscous coef - + real(kind=dbl_kind), parameter :: capping = c1 ! of the viscous coef + character(len=*), parameter :: subname = '(stress_T)' !----------------------------------------------------------------- @@ -1435,8 +1424,8 @@ subroutine stress_U (nx_block, ny_block, & stresspU, stressmU, & stress12U, Tarea ) - use ice_dyn_shared, only: strain_rates_U, & - viscous_coeffs_and_rep_pressure_T2U + 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 @@ -1449,14 +1438,14 @@ 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 @@ -1512,17 +1501,15 @@ subroutine stress_U (nx_block, ny_block, & ! viscous coefficients and replacement pressure at T point !----------------------------------------------------------------- - - 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) + 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 1d95b2ba8..a999b611b 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -1034,7 +1034,7 @@ subroutine dyn_finish (nx_block, ny_block, & if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - if (present(strocnxT)) then + if (present(strocnxT) .and. present(strocnyT)) then do j = 1, ny_block do i = 1, nx_block @@ -1070,7 +1070,7 @@ subroutine dyn_finish (nx_block, ny_block, & ! strocnx(i,j) = -(strairx(i,j) + strintx(i,j)) ! strocny(i,j) = -(strairy(i,j) + strinty(i,j)) - if (present(strocnxT)) then + if (present(strocnxT) .and. present(strocnyT)) then ! Prepare to convert to T grid ! divide by aice for coupling @@ -1947,35 +1947,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 - 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 + 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 @@ -2031,8 +2026,8 @@ subroutine viscous_coeffs_and_rep_pressure_T2U (zetax2T_00, zetax2T_01, & etax2T_11, etax2T_10, & maskT_00, maskT_01, & maskT_11, maskT_10, & - Tarea_00, Tarea_01, & - Tarea_11, Tarea_10, & + tarea_00, tarea_01, & + tarea_11, tarea_10, & deltaU, & zetax2U, etax2U, & rep_prsU) @@ -2042,7 +2037,7 @@ subroutine viscous_coeffs_and_rep_pressure_T2U (zetax2T_00, zetax2T_01, & 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, & + tarea_00, tarea_10, tarea_11, tarea_01, & deltaU real (kind=dbl_kind), intent(out):: zetax2U, etax2U, rep_prsU diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index fd7717b88..4bb2710b8 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -78,6 +78,8 @@ module ice_grid earea , & ! area of E-cell (m^2) tarear , & ! 1/tarea uarear , & ! 1/uarea + narear , & ! 1/narea + earear , & ! 1/earea tinyarea,& ! puny*tarea tarean , & ! area of NH T-cells tareas , & ! area of SH T-cells @@ -165,8 +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) (calc point if all point is ocean) - umaskCD, & ! land/boundary mask, velocity (U-cell) (calc point if one point is ocean) + 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 @@ -209,6 +211,8 @@ subroutine alloc_grid earea (nx_block,ny_block,max_blocks), & ! area of E-cell (m^2) tarear (nx_block,ny_block,max_blocks), & ! 1/tarea uarear (nx_block,ny_block,max_blocks), & ! 1/uarea + narear (nx_block,ny_block,max_blocks), & ! 1/narea + earear (nx_block,ny_block,max_blocks), & ! 1/earea tinyarea (nx_block,ny_block,max_blocks), & ! puny*tarea tarean (nx_block,ny_block,max_blocks), & ! area of NH T-cells tareas (nx_block,ny_block,max_blocks), & ! area of SH T-cells @@ -242,7 +246,7 @@ subroutine alloc_grid 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) + 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 @@ -499,6 +503,16 @@ subroutine init_grid2 else uarear(i,j,iblk) = c0 ! possible on boundaries endif + if (narea(i,j,iblk) > c0) then + narear(i,j,iblk) = c1/narea(i,j,iblk) + else + narear(i,j,iblk) = c0 ! possible on boundaries + endif + if (earea(i,j,iblk) > c0) then + earear(i,j,iblk) = c1/earea(i,j,iblk) + else + earear(i,j,iblk) = c0 ! possible on boundaries + endif tinyarea(i,j,iblk) = puny*tarea(i,j,iblk) enddo enddo @@ -1924,7 +1938,7 @@ subroutine makemask real (kind=dbl_kind) :: & puny - real (kind=dbl_kind) :: & + real (kind=dbl_kind), dimension(:,:,:), allocatable :: & uvmCD type (block) :: & @@ -1949,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 @@ -1965,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 @@ -1973,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, & @@ -1990,30 +2009,32 @@ subroutine makemask jhi = this_block%jhi ! needs to cover halo (no halo update for logicals) - tmask(:,:,iblk) = .false. - umask(:,:,iblk) = .false. + tmask(:,:,iblk) = .false. + umask(:,:,iblk) = .false. umaskCD(:,:,iblk) = .false. - nmask(:,:,iblk) = .false. - emask(:,:,iblk) = .false. - do j = jlo-nghost, jhi - do i = ilo-nghost, ihi -!!!!!WARNING THIS SUM is not calculating uvmCD mask on the northern and eastern halo -! Loop had to be separate as plus 1 indexes exceed array size - uvmCD = (hm(i,j, iblk)+hm(i+1,j, iblk)+ & - hm(i,j+1,iblk)+hm(i+1,j+1,iblk)) - if (uvmCD > c1p5) umaskCD(i,j,iblk) = .true. - enddo - enddo - + 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 !----------------------------------------------------------------- @@ -2047,6 +2068,8 @@ subroutine makemask enddo ! iblk !$OMP END PARALLEL DO + deallocate(uvmCD) + end subroutine makemask !======================================================================= @@ -2875,7 +2898,7 @@ real(kind=dbl_kind) function grid_neighbor_min(field, i, j, grid_location) resul case('N') mini = min(field(i,j), field(i,j+1)) case default - call abort_ice(subname // ' unkwown grid_location: ' // grid_location) + call abort_ice(subname // ' unknown grid_location: ' // grid_location) end select end function grid_neighbor_min @@ -2907,7 +2930,7 @@ real(kind=dbl_kind) function grid_neighbor_max(field, i, j, grid_location) resul case('N') maxi = max(field(i,j), field(i,j+1)) case default - call abort_ice(subname // ' unkwown grid_location: ' // grid_location) + call abort_ice(subname // ' unknown grid_location: ' // grid_location) end select end function grid_neighbor_max diff --git a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 b/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 index 9a5a75bea..11836c073 100644 --- a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 +++ b/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 @@ -381,6 +381,9 @@ subroutine restartfile (ice_ic) call read_restart_field(nu_restart,0,stress12_4,'ruf8', & 'stress12_4',1,diag,field_loc_center,field_type_scalar) ! stress12_4 +! tcraig, comment these out now to allow restarts from B grid file +! this will affect exact restart when we get to that point +#if (1 == 0) if (grid_system == 'CD') then call read_restart_field(nu_restart,0,stresspT,'ruf8', & 'stresspT' ,1,diag,field_loc_center,field_type_scalar) ! stresspT @@ -395,6 +398,7 @@ subroutine restartfile (ice_ic) call read_restart_field(nu_restart,0,stress12U,'ruf8', & 'stress12U',1,diag,field_loc_center,field_type_scalar) ! stress12U endif +#endif if (trim(grid_type) == 'tripole') then call ice_HaloUpdate_stress(stressp_1, stressp_3, halo_info, & diff --git a/configuration/scripts/options/set_nml.boxnodyn b/configuration/scripts/options/set_nml.boxnodyn index 35a224c92..8e5d4a692 100644 --- a/configuration/scripts/options/set_nml.boxnodyn +++ b/configuration/scripts/options/set_nml.boxnodyn @@ -8,6 +8,30 @@ dumpfreq_n = 2 histfreq = 'd','x','x','x','x' histfreq_n = 2,1,1,1,1 f_aice = 'd' +f_uvel = 'd' +f_vvel = 'd' +f_uatm = 'd' +f_vatm = 'd' +f_uocn = 'd' +f_vocn = 'd' +f_strairx = 'd' +f_strairy = 'd' +f_strtltx = 'd' +f_strtlty = 'd' +f_strcorx = 'd' +f_strcory = 'd' +f_strocnx = 'd' +f_strocny = 'd' +f_strintx = 'd' +f_strinty = 'd' +f_taubx = 'd' +f_tauby = 'd' +f_strength = 'd' +f_divu = 'd' +f_shear = 'd' +f_sig1 = 'd' +f_sig2 = 'd' +f_sigP = 'd' kcatbound = 0 ew_boundary_type = 'open' ns_boundary_type = 'open' @@ -20,10 +44,13 @@ tr_pond_lvl = .false. tr_aero = .false. kitd = 0 ktherm = 0 -kdyn = 0 +kdyn = 1 revised_evp = .false. -kstrength = 0 +kstrength = 1 krdg_partic = 1 krdg_redist = 1 +seabed_stress = .true. +atm_data_type = 'calm' +ocn_data_type = 'calm' shortwave = 'ccsm3' albedo_type = 'constant' From 9b1ec14dfce17a34ebbd8edbc53b247f3d1c2760 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Sun, 21 Nov 2021 20:17:57 +0000 Subject: [PATCH 10/13] update ice_dyn_evp --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 2c4df2ac2..6ca39f84f 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -690,11 +690,12 @@ subroutine evp (dt) dxU (:,:,iblk), dyU (:,:,iblk), & ratiodxN (:,:,iblk), ratiodxNr (:,:,iblk), & ratiodyE (:,:,iblk), ratiodyEr (:,:,iblk), & + tarea (:,:,iblk), & epm (:,:,iblk), npm (:,:,iblk), & hm (:,:,iblk), uvm (:,:,iblk), & zetax2T (:,:,iblk), etax2T (:,:,iblk), & stresspU (:,:,iblk), stressmU (:,:,iblk), & - stress12U (:,:,iblk), tarea (:,:,iblk)) + stress12U (:,:,iblk)) call step_vel (nx_block, ny_block, & ! E point icelle (iblk), Cdn_ocn (:,:,iblk), & @@ -951,7 +952,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, viscous_coeffs_and_rep_pressure_T + 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 @@ -1417,12 +1418,13 @@ subroutine stress_U (nx_block, ny_block, & uvelU, vvelU, & dxE, dyN, & dxU, dyU, & + tarea, & ratiodxN, ratiodxNr, & ratiodyE, ratiodyEr, & epm, npm, hm, uvm, & zetax2T, etax2T, & stresspU, stressmU, & - stress12U, Tarea ) + stress12U ) use ice_dyn_shared, only: strain_rates_U, & viscous_coeffs_and_rep_pressure_T2U @@ -1447,14 +1449,14 @@ subroutine stress_U (nx_block, ny_block, & dyN , & ! height of N-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 - Tarea , & ! area of T-cell uvm , & ! U-cell mask zetax2T , & ! 2*zeta at the T point etax2T ! 2*eta at the T point @@ -1629,7 +1631,7 @@ subroutine div_stress (nx_block, ny_block, & + (c1/dyE_N(i,j)) * ( (dyT_U(i,j)**2) * stress12(i,j) & -(dyT_U(i-1,j)**2)*stress12(i-1,j) ) ) case default - call abort_ice(subname // ' unkwown grid_location: ' // grid_location) + call abort_ice(subname // ' unknown grid_location: ' // grid_location) end select From f5d7d89d1a3ef318ea2b01c61db873e1d47ff22b Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Sun, 21 Nov 2021 22:51:25 +0000 Subject: [PATCH 11/13] It is now bit for bit. Number of mask is reduced from 2 to one. Instead an additional if else gridsystem is introduced. subroutine evp should be refactored --- cicecore/cicedynB/dynamics/ice_dyn_eap.F90 | 4 +- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 93 ++++++++++++------- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 49 ++-------- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 4 +- 4 files changed, 75 insertions(+), 75 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index b0b47c5cb..e5a89b118 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -133,7 +133,7 @@ subroutine eap (dt) stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 - use ice_grid, only: tmask, umask, umaskCD, dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, & + use ice_grid, only: tmask, umask, dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, & tarear, uarear, grid_average_X2Y!, grid_system commented out until implementation of cd-grid use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, & aice_init, aice0, aicen, vicen, strength !, uvelE, vvelN grid_system commented out until implementation of cd-grid @@ -301,7 +301,7 @@ subroutine eap (dt) indxui (:,iblk), indxuj (:,iblk), & aiu (:,:,iblk), umass (:,:,iblk), & umassdti (:,:,iblk), fcor_blk (:,:,iblk), & - umask (:,:,iblk), umaskCD (:,:,iblk), & + umask (:,:,iblk), & uocn (:,:,iblk), vocn (:,:,iblk), & strairx (:,:,iblk), strairy (:,:,iblk), & ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), & diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 6ca39f84f..9ba8f5b41 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -341,35 +341,66 @@ 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), 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)) + + if (grid_system == trim('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)) + else + 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 @@ -413,7 +444,7 @@ subroutine evp (dt) indxni (:,iblk), indxnj (:,iblk), & aiN (:,:,iblk), nmass (:,:,iblk), & nmassdti (:,:,iblk), fcorN_blk (:,:,iblk),& - nmask (:,:,iblk), nmask (:,:,iblk), & + nmask (:,:,iblk), & uocnN (:,:,iblk), vocnN (:,:,iblk), & strairxN (:,:,iblk), strairyN (:,:,iblk), & ss_tltxN (:,:,iblk), ss_tltyN (:,:,iblk), & @@ -446,7 +477,7 @@ subroutine evp (dt) indxei (:,iblk), indxej (:,iblk), & aiE (:,:,iblk), emass (:,:,iblk), & emassdti (:,:,iblk), fcorE_blk (:,:,iblk),& - emask (:,:,iblk), emask (:,:,iblk),& + emask (:,:,iblk), & uocnE (:,:,iblk), vocnE (:,:,iblk), & strairxE (:,:,iblk), strairyE (:,:,iblk), & ss_tltxE (:,:,iblk), ss_tltyE (:,:,iblk), & diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index a999b611b..bcd7ed247 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -453,7 +453,7 @@ subroutine dyn_prep2 (nx_block, ny_block, & indxui, indxuj, & aiu, umass, & umassdti, fcor, & - umask, umaskCD, & + umask, & uocn, vocn, & strairx, strairy, & ss_tltx, ss_tlty, & @@ -490,7 +490,7 @@ subroutine dyn_prep2 (nx_block, ny_block, & indxuj ! compressed index in j-direction logical (kind=log_kind), dimension (nx_block,ny_block), intent(in) :: & - umask, umaskCD ! land/boundary mask, thickness (U-cell) + umask ! land/boundary mask, thickness (U-cell) integer (kind=int_kind), dimension (nx_block,ny_block), intent(in) :: & icetmask ! ice extent mask (T-cell) @@ -606,48 +606,18 @@ subroutine dyn_prep2 (nx_block, ny_block, & !----------------------------------------------------------------- icellu = 0 - if (grid_system == 'B') then + do j = jlo, jhi do i = ilo, ihi - - ! ice extent mask (U-cells) iceumask_old(i,j) = iceumask(i,j) ! save - iceumask(i,j) = (umask(i,j)) .and. (aiu (i,j) > a_min) & - .and. (umass(i,j) > m_min) - - if (iceumask(i,j)) then - icellu = icellu + 1 - indxui(icellu) = i - indxuj(icellu) = j - - ! initialize velocity for new ice points to ocean sfc current - if (.not. iceumask_old(i,j)) then - uvel(i,j) = uocn(i,j) - vvel(i,j) = vocn(i,j) - endif - else - ! set velocity and stresses to zero for masked-out points - uvel(i,j) = c0 - vvel(i,j) = c0 - strintx(i,j) = c0 - strinty(i,j) = c0 - strocnx(i,j) = c0 - strocny(i,j) = c0 + 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) + else ! ice mask shpuld be applied to cd grid. For now it is not implemented. + iceumask(i,j) = umask(i,j) endif - uvel_init(i,j) = uvel(i,j) - vvel_init(i,j) = vvel(i,j) - enddo - enddo - else - do j = jlo, jhi - do i = ilo, ihi - - ! ice extent mask (U-cells) -! iceumask_old(i,j) = iceumask(i,j) ! save -! iceumask(i,j) = (umaskCD(i,j)) .and. (aiu (i,j) > a_min) & -! .and. (umass(i,j) > m_min) - iceumask(i,j) = umaskCD(i,j) ! SHOULD BE EXTENDED TO CHECK FOR ICE if (iceumask(i,j)) then icellu = icellu + 1 indxui(icellu) = i @@ -672,7 +642,6 @@ subroutine dyn_prep2 (nx_block, ny_block, & vvel_init(i,j) = vvel(i,j) enddo enddo - endif !----------------------------------------------------------------- ! Define variables for momentum equation !----------------------------------------------------------------- diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index bbbcb67c5..5bbef02a4 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -199,7 +199,7 @@ subroutine implicit_solver (dt) stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 - use ice_grid, only: tmask, umask, umaskCD, dxt, dyt, cxp, cyp, cxm, cym, & + use ice_grid, only: tmask, umask, dxt, dyt, cxp, cyp, cxm, cym, & tarear, grid_type, grid_average_X2Y !, grid_system commented out until implementation of c grid use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, & aice_init, aice0, aicen, vicen, strength!, uvelE, vvelN ommented out until implementation of c grid @@ -366,7 +366,7 @@ subroutine implicit_solver (dt) indxui (:,iblk), indxuj (:,iblk), & aiu (:,:,iblk), umass (:,:,iblk), & umassdti (:,:,iblk), fcor_blk (:,:,iblk), & - umask (:,:,iblk), umaskCD (:,:,iblk), & + umask (:,:,iblk), & uocn (:,:,iblk), vocn (:,:,iblk), & strairx (:,:,iblk), strairy (:,:,iblk), & ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), & From 65d46c08111e24e90cedb4e055c91d970277a943 Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Tue, 23 Nov 2021 19:24:37 +0000 Subject: [PATCH 12/13] moved trim --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 9ba8f5b41..da1840f4d 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -342,7 +342,7 @@ subroutine evp (dt) jlo = this_block%jlo jhi = this_block%jhi - if (grid_system == trim('B')) then + if (trim(grid_system) == 'B') then call dyn_prep2 (nx_block, ny_block, & ilo, ihi, jlo, jhi, & icellt(iblk), icellu(iblk), & From 82c90b2143a339cf942ebab8ab3190d2bcc3eecb Mon Sep 17 00:00:00 2001 From: TillRasmussen Date: Tue, 23 Nov 2021 23:57:33 +0000 Subject: [PATCH 13/13] Based on Quicksuite: B grid bit for bit. CD grid passes 13 of 16. 2 fails on 2 test. 1 is pending. Comments from Phillipe also included --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 14 +++++++------- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 15 ++++++--------- cicecore/cicedynB/infrastructure/ice_grid.F90 | 4 ++-- 3 files changed, 15 insertions(+), 18 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index da1840f4d..e1cf14e3a 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -342,7 +342,7 @@ subroutine evp (dt) jlo = this_block%jlo jhi = this_block%jhi - if (trim(grid_system) == 'B') then + if (trim(grid_system) == 'B') then call dyn_prep2 (nx_block, ny_block, & ilo, ihi, jlo, jhi, & icellt(iblk), icellu(iblk), & @@ -371,7 +371,8 @@ subroutine evp (dt) uvel_init (:,:,iblk), vvel_init (:,:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & Tbu (:,:,iblk)) - else + + elseif (trim(grid_system) == 'CD') then call dyn_prep2 (nx_block, ny_block, & ilo, ihi, jlo, jhi, & icellt(iblk), icellu(iblk), & @@ -400,7 +401,7 @@ subroutine evp (dt) uvel_init (:,:,iblk), vvel_init (:,:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & Tbu (:,:,iblk)) - endif + endif !----------------------------------------------------------------- ! ice strength @@ -444,7 +445,7 @@ subroutine evp (dt) indxni (:,iblk), indxnj (:,iblk), & aiN (:,:,iblk), nmass (:,:,iblk), & nmassdti (:,:,iblk), fcorN_blk (:,:,iblk),& - nmask (:,:,iblk), & + nmask (:,:,iblk), & uocnN (:,:,iblk), vocnN (:,:,iblk), & strairxN (:,:,iblk), strairyN (:,:,iblk), & ss_tltxN (:,:,iblk), ss_tltyN (:,:,iblk), & @@ -677,7 +678,6 @@ subroutine evp (dt) !----------------------------------------------------------------- ! momentum equation !----------------------------------------------------------------- - call stepu (nx_block, ny_block, & icellu (iblk), Cdn_ocn (:,:,iblk), & indxui (:,iblk), indxuj (:,iblk), & @@ -719,15 +719,15 @@ 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), & - tarea (:,:,iblk), & epm (:,:,iblk), npm (:,:,iblk), & hm (:,:,iblk), uvm (:,:,iblk), & zetax2T (:,:,iblk), etax2T (:,:,iblk), & stresspU (:,:,iblk), stressmU (:,:,iblk), & stress12U (:,:,iblk)) - + call step_vel (nx_block, ny_block, & ! E point icelle (iblk), Cdn_ocn (:,:,iblk), & indxei (:,iblk), indxej (:,iblk), & diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index bcd7ed247..d4622cb6e 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -610,13 +610,13 @@ subroutine dyn_prep2 (nx_block, ny_block, & do j = jlo, jhi do i = ilo, ihi iceumask_old(i,j) = iceumask(i,j) ! save - if (grid_system == 'B') then ! include ice mask. +! 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) & + iceumask(i,j) = (umask(i,j)) .and. (aiu (i,j) > a_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 +! 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 @@ -920,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 @@ -2023,7 +2021,6 @@ subroutine viscous_coeffs_and_rep_pressure_T2U (zetax2T_00, zetax2T_01, & 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 + & diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index 4bb2710b8..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) - land if one point is one + 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) @@ -2018,7 +2018,7 @@ subroutine makemask 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 (uvmCD(i,j,iblk) > c1p5) umaskCD(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