diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 1e308b65a..aefb7679b 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -42,7 +42,8 @@ module ice_dyn_evp use ice_constants, only: c0, p027, p055, p111, p166, & p222, p25, p333, p5, c1 use ice_dyn_shared, only: stepu, dyn_prep1, dyn_prep2, dyn_finish, & - ndte, yield_curve, ecci, denom1, arlx1i, fcor_blk, uvel_init, vvel_init, & + ndte, yield_curve, ecci, denom1, arlx1i, fcor_blk, fcorE_blk, fcorN_blk, & + uvel_init, vvel_init, uvelE_init, vvelE_init, uvelN_init, vvelN_init, & seabed_stress_factor_LKD, seabed_stress_factor_prob, seabed_stress_method, & seabed_stress, Ktens, revp use ice_fileunits, only: nu_diag @@ -84,14 +85,22 @@ subroutine evp (dt) strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, & strtltx, strtlty, strocnx, strocny, strintx, strinty, taubx, tauby, & strocnxT, strocnyT, strax, stray, & + strairxN, strairyN, uocnN, vocnN, ss_tltxN, ss_tltyN, icenmask, fmN, & + strtltxN, strtltyN, strocnxN, strocnyN, strintxN, strintyN, taubxN, taubyN, & + straxN, strayN, TbN, & + strairxE, strairyE, uocnE, vocnE, ss_tltxE, ss_tltyE, iceemask, fmE, & + strtltxE, strtltyE, strocnxE, strocnyE, strintxE, strintyE, taubxE, taubyE, & + straxE, strayE, TbE, & Tbu, hwater, & 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, nmask, emask, dxt, dyt, & + dxhy, dyhx, cxp, cyp, cxm, cym, & tarear, uarear, tinyarea, grid_average_X2Y, & grid_type, grid_system - use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, & + use ice_state, only: aice, vice, vsno, uvel, vvel, uvelN, vvelN, & + uvelE, vvelE, divu, shear, & aice_init, aice0, aicen, vicen, strength use ice_timers, only: timer_dynamics, timer_bound, & ice_timer_start, ice_timer_stop, timer_evp_1d, timer_evp_2d @@ -112,11 +121,17 @@ subroutine evp (dt) integer (kind=int_kind), dimension(max_blocks) :: & icellt , & ! no. of cells where icetmask = 1 + icelln , & ! no. of cells where icenmask = 1 + icelle , & ! no. of cells where iceemask = 1 icellu ! no. of cells where iceumask = 1 integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks) :: & indxti , & ! compressed index in i-direction indxtj , & ! compressed index in j-direction + indxei , & ! compressed index in i-direction + indxej , & ! compressed index in j-direction + indxni , & ! compressed index in i-direction + indxnj , & ! compressed index in j-direction indxui , & ! compressed index in i-direction indxuj ! compressed index in j-direction @@ -130,6 +145,24 @@ subroutine evp (dt) umass , & ! total mass of ice and snow (u grid) umassdti ! mass of U-cell/dte (kg/m^2 s) + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & + waterxN , & ! for ocean stress calculation, x (m/s) + wateryN , & ! for ocean stress calculation, y (m/s) + forcexN , & ! work array: combined atm stress and ocn tilt, x + forceyN , & ! work array: combined atm stress and ocn tilt, y + aiN , & ! ice fraction on N-grid + nmass , & ! total mass of ice and snow (N grid) + nmassdti ! mass of N-cell/dte (kg/m^2 s) + + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & + waterxE , & ! for ocean stress calculation, x (m/s) + wateryE , & ! for ocean stress calculation, y (m/s) + forcexE , & ! work array: combined atm stress and ocn tilt, x + forceyE , & ! work array: combined atm stress and ocn tilt, y + aiE , & ! ice fraction on E-grid + emass , & ! total mass of ice and snow (E grid) + emassdti ! mass of E-cell/dte (kg/m^2 s) + real (kind=dbl_kind), allocatable :: fld2(:,:,:,:) real (kind=dbl_kind), dimension(nx_block,ny_block,8):: & @@ -207,6 +240,12 @@ subroutine evp (dt) strairx (:,:,iblk), strairy (:,:,iblk), & tmass (:,:,iblk), icetmask(:,:,iblk)) + if (grid_system == 'CD') then + strairxN(:,:,iblk) = strairxT(:,:,iblk) + strairyN(:,:,iblk) = strairyT(:,:,iblk) + strairxE(:,:,iblk) = strairxT(:,:,iblk) + strairyE(:,:,iblk) = strairyT(:,:,iblk) + endif enddo ! iblk !$OMP END PARALLEL DO @@ -222,6 +261,12 @@ subroutine evp (dt) call grid_average_X2Y('T2UF',tmass,umass) call grid_average_X2Y('T2UF',aice_init, aiu) + if (grid_system == 'CD') then + call grid_average_X2Y('T2EF',tmass,emass) + call grid_average_X2Y('T2EF',aice_init, aie) + call grid_average_X2Y('T2NF',tmass,nmass) + call grid_average_X2Y('T2NF',aice_init, aie) + endif !---------------------------------------------------------------- ! Set wind stress to values supplied via NEMO or other forcing ! This wind stress is rotated on u grid and multiplied by aice @@ -243,6 +288,29 @@ subroutine evp (dt) call grid_average_X2Y('T2UF',strairy) endif + if (grid_system == 'CD') then + + if (.not. calc_strair) then + strairxN(:,:,:) = strax(:,:,:) + strairyN(:,:,:) = stray(:,:,:) + strairxE(:,:,:) = strax(:,:,:) + strairyE(:,:,:) = stray(:,:,:) + else + call ice_HaloUpdate (strairxN, halo_info, & + field_loc_center, field_type_vector) + call ice_HaloUpdate (strairyN, halo_info, & + field_loc_center, field_type_vector) + call ice_HaloUpdate (strairxE, halo_info, & + field_loc_center, field_type_vector) + call ice_HaloUpdate (strairyE, halo_info, & + field_loc_center, field_type_vector) + call grid_average_X2Y('T2NF',strairxN) + call grid_average_X2Y('T2NF',strairyN) + call grid_average_X2Y('T2EF',strairxE) + call grid_average_X2Y('T2EF',strairyE) + endif + + endif ! tcraig, tcx, threading here leads to some non-reproducbile results and failures in icepack_ice_strength ! need to do more debugging !$TCXOMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) @@ -307,6 +375,100 @@ subroutine evp (dt) enddo ! iblk !$TCXOMP END PARALLEL DO + if (grid_system == 'CD') then + + !$TCXOMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) + do iblk = 1, nblocks + + !----------------------------------------------------------------- + ! more preparation for dynamics on N grid + !----------------------------------------------------------------- + + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + 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), icelln(iblk), & + indxti (:,iblk), indxtj (:,iblk), & + indxni (:,iblk), indxnj (:,iblk), & + aiN (:,:,iblk), nmass (:,:,iblk), & + nmassdti (:,:,iblk), fcorN_blk (:,:,iblk), & + nmask (:,:,iblk), & + uocnN (:,:,iblk), vocnN (:,:,iblk), & + strairxN (:,:,iblk), strairyN (:,:,iblk), & + ss_tltxN (:,:,iblk), ss_tltyN (:,:,iblk), & + icetmask (:,:,iblk), icenmask (:,:,iblk), & + fmN (:,:,iblk), dt, & + strtltxN (:,:,iblk), strtltyN (:,:,iblk), & + strocnxN (:,:,iblk), strocnyN (:,:,iblk), & + strintxN (:,:,iblk), strintyN (:,:,iblk), & + taubxN (:,:,iblk), taubyN (:,:,iblk), & + waterxN (:,:,iblk), wateryN (:,:,iblk), & + forcexN (:,:,iblk), forceyN (:,:,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), & + uvelN_init (:,:,iblk), vvelN_init (:,:,iblk), & + uvelN (:,:,iblk), vvelN (:,:,iblk), & + TbN (:,:,iblk)) + + enddo ! iblk + !$TCXOMP END PARALLEL DO + + !$TCXOMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) + do iblk = 1, nblocks + + !----------------------------------------------------------------- + ! more preparation for dynamics on N grid + !----------------------------------------------------------------- + + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + 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), icelle(iblk), & + indxti (:,iblk), indxtj (:,iblk), & + indxei (:,iblk), indxej (:,iblk), & + aiE (:,:,iblk), emass (:,:,iblk), & + emassdti (:,:,iblk), fcorE_blk (:,:,iblk), & + emask (:,:,iblk), & + uocnE (:,:,iblk), vocnE (:,:,iblk), & + strairxE (:,:,iblk), strairyE (:,:,iblk), & + ss_tltxE (:,:,iblk), ss_tltyE (:,:,iblk), & + icetmask (:,:,iblk), icenmask (:,:,iblk), & + fmE (:,:,iblk), dt, & + strtltxE (:,:,iblk), strtltyE (:,:,iblk), & + strocnxE (:,:,iblk), strocnyE (:,:,iblk), & + strintxE (:,:,iblk), strintyE (:,:,iblk), & + taubxE (:,:,iblk), taubyE (:,:,iblk), & + waterxE (:,:,iblk), wateryE (:,:,iblk), & + forcexE (:,:,iblk), forceyE (:,:,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), & + uvelE_init (:,:,iblk), vvelE_init (:,:,iblk), & + uvelE (:,:,iblk), vvelE (:,:,iblk), & + TbE (:,:,iblk)) + + enddo ! iblk + !$TCXOMP END PARALLEL DO + + endif ! grid_system + call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 1d27ff792..2f569d62f 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -84,10 +84,22 @@ module ice_dyn_shared real (kind=dbl_kind), allocatable, public :: & fcor_blk(:,:,:) ! Coriolis parameter (1/s) + real (kind=dbl_kind), allocatable, public :: & + fcorE_blk(:,:,:), & ! Coriolis parameter at E points (1/s) + fcorN_blk(:,:,:) ! Coriolis parameter at N points (1/s) + real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & uvel_init, & ! x-component of velocity (m/s), beginning of timestep vvel_init ! y-component of velocity (m/s), beginning of timestep + real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & + uvelN_init, & ! x-component of velocity (m/s), beginning of timestep + vvelN_init ! y-component of velocity (m/s), beginning of timestep + + real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & + uvelE_init, & ! x-component of velocity (m/s), beginning of timestep + vvelE_init ! y-component of velocity (m/s), beginning of timestep + ! ice isotropic tensile strength parameter real (kind=dbl_kind), public :: & Ktens ! T=Ktens*P (tensile strength: see Konig and Holland, 2010) @@ -121,6 +133,16 @@ subroutine alloc_dyn_shared stat=ierr) if (ierr/=0) call abort_ice('(alloc_dyn_shared): Out of memory') + if (grid_system == 'CD') then + allocate( & + uvelE_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep + vvelE_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep + uvelN_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep + vvelN_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep + stat=ierr) + if (ierr/=0) call abort_ice('(alloc_dyn_shared): Out of memory') + endif + end subroutine alloc_dyn_shared !======================================================================= @@ -138,7 +160,7 @@ subroutine init_dyn (dt) stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 use ice_state, only: uvel, vvel, uvelE, vvelE, uvelN, vvelN, divu, shear - use ice_grid, only: ULAT + use ice_grid, only: ULAT, NLAT, ELAT real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -161,6 +183,11 @@ subroutine init_dyn (dt) allocate(fcor_blk(nx_block,ny_block,max_blocks)) + if (grid_system == 'CD') then + allocate(fcorE_blk(nx_block,ny_block,max_blocks)) + allocate(fcorN_blk(nx_block,ny_block,max_blocks)) + endif + !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks do j = 1, ny_block @@ -170,10 +197,10 @@ subroutine init_dyn (dt) uvel(i,j,iblk) = c0 ! m/s vvel(i,j,iblk) = c0 ! m/s if (grid_system == 'CD') then ! extra velocity variables - uvelE = c0 - vvelE = c0 - uvelN = c0 - vvelN = c0 + uvelE(i,j,iblk) = c0 + vvelE(i,j,iblk) = c0 + uvelN(i,j,iblk) = c0 + vvelN(i,j,iblk) = c0 endif ! strain rates @@ -191,6 +218,21 @@ subroutine init_dyn (dt) fcor_blk(i,j,iblk) = c2*omega*sin(ULAT(i,j,iblk)) ! 1/s endif + if (grid_system == 'CD') then + + if (trim(coriolis) == 'constant') then + fcorE_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s + fcorN_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s + else if (trim(coriolis) == 'zero') then + fcorE_blk(i,j,iblk) = 0.0 + fcorN_blk(i,j,iblk) = 0.0 + else + fcorE_blk(i,j,iblk) = c2*omega*sin(ELAT(i,j,iblk)) ! 1/s + fcorN_blk(i,j,iblk) = c2*omega*sin(NLAT(i,j,iblk)) ! 1/s + endif + + endif + ! stress tensor, kg/s^2 stressp_1 (i,j,iblk) = c0 stressp_2 (i,j,iblk) = c0 diff --git a/cicecore/cicedynB/general/ice_flux.F90 b/cicecore/cicedynB/general/ice_flux.F90 index 1b17c130d..65c21f6d9 100644 --- a/cicecore/cicedynB/general/ice_flux.F90 +++ b/cicecore/cicedynB/general/ice_flux.F90 @@ -40,12 +40,24 @@ module ice_flux ! in from atmos (if .not.calc_strair) strax , & ! wind stress components (N/m^2) stray , & ! + straxE , & ! wind stress components (N/m^2) + strayE , & ! + straxN , & ! wind stress components (N/m^2) + strayN , & ! ! in from ocean uocn , & ! ocean current, x-direction (m/s) vocn , & ! ocean current, y-direction (m/s) + uocnE , & ! ocean current, x-direction (m/s) + vocnE , & ! ocean current, y-direction (m/s) + uocnN , & ! ocean current, x-direction (m/s) + vocnN , & ! ocean current, y-direction (m/s) ss_tltx , & ! sea surface slope, x-direction (m/m) ss_tlty , & ! sea surface slope, y-direction + ss_tltxE, & ! sea surface slope, x-direction (m/m) + ss_tltyE, & ! sea surface slope, y-direction + ss_tltxN, & ! sea surface slope, x-direction (m/m) + ss_tltyN, & ! sea surface slope, y-direction hwater , & ! water depth for seabed stress calc (landfast ice) ! out to atmosphere @@ -128,6 +140,14 @@ module ice_flux dimension (:,:,:), allocatable, public :: & iceumask ! ice extent mask (U-cell) + logical (kind=log_kind), & + dimension (:,:,:), allocatable, public :: & + icenmask ! ice extent mask (N-cell) + + logical (kind=log_kind), & + dimension (:,:,:), allocatable, public :: & + iceemask ! ice extent mask (E-cell) + ! internal real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & @@ -565,6 +585,12 @@ subroutine alloc_flux if (grid_system == "CD") & allocate( & + straxN (nx_block,ny_block,max_blocks), & ! wind stress components (N/m^2) + strayN (nx_block,ny_block,max_blocks), & ! + uocnN (nx_block,ny_block,max_blocks), & ! ocean current, x-direction (m/s) + vocnN (nx_block,ny_block,max_blocks), & ! ocean current, y-direction (m/s) + ss_tltxN (nx_block,ny_block,max_blocks), & ! sea surface slope, x-direction (m/m) + ss_tltyN (nx_block,ny_block,max_blocks), & ! sea surface slope, y-direction taubxN (nx_block,ny_block,max_blocks), & ! seabed stress (x) at N points (N/m^2) taubyN (nx_block,ny_block,max_blocks), & ! seabed stress (y) at N points (N/m^2) strairxN (nx_block,ny_block,max_blocks), & ! stress on ice by air, x-direction at N points @@ -575,6 +601,13 @@ subroutine alloc_flux strtltyN (nx_block,ny_block,max_blocks), & ! stress due to sea surface slope, y-direction at N points strintxN (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, x at N points (N/m^2) strintyN (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, y at N points (N/m^2) + icenmask (nx_block,ny_block,max_blocks), & ! ice extent mask (N-cell) + straxE (nx_block,ny_block,max_blocks), & ! wind stress components (N/m^2) + strayE (nx_block,ny_block,max_blocks), & ! + uocnE (nx_block,ny_block,max_blocks), & ! ocean current, x-direction (m/s) + vocnE (nx_block,ny_block,max_blocks), & ! ocean current, y-direction (m/s) + ss_tltxE (nx_block,ny_block,max_blocks), & ! sea surface slope, x-direction (m/m) + ss_tltyE (nx_block,ny_block,max_blocks), & ! sea surface slope, y-direction taubxE (nx_block,ny_block,max_blocks), & ! seabed stress (x) at E points (N/m^2) taubyE (nx_block,ny_block,max_blocks), & ! seabed stress (y) at E points (N/m^2) strairxE (nx_block,ny_block,max_blocks), & ! stress on ice by air, x-direction at E points @@ -585,6 +618,7 @@ subroutine alloc_flux strtltyE (nx_block,ny_block,max_blocks), & ! stress due to sea surface slope, y-direction at E points strintxE (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, x at E points (N/m^2) strintyE (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, y at E points (N/m^2) + iceemask (nx_block,ny_block,max_blocks), & ! ice extent mask (E-cell) stat=ierr) if (ierr/=0) call abort_ice('(alloc_flux): Out of memory')