diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 367adebcd..2ecc6e5f6 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -192,16 +192,12 @@ subroutine init_dyn (dt) use ice_blocks, only: nx_block, ny_block use ice_domain, only: nblocks, halo_dynbundle use ice_domain_size, only: max_blocks - use ice_flux, only: rdg_conv, rdg_shear, iceumask, & + use ice_flux, only: rdg_conv, rdg_shear, iceumask, iceemask, icenmask, & stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4, & stresspT, stressmT, stress12T, & - stresspU, stressmU, stress12U, & - taubxE, taubyE, strairxE, strairyE, strocnxE, strocnyE, & - strtltxE, strtltyE, strintxE, strintyE, iceEmask, fmE, TbE, & - taubxN, taubyN, strairxN, strairyN, strocnxN, strocnyN, & - strtltxN, strtltyN, strintxN, strintyN, iceNmask, fmN, TbN + stresspU, stressmU, stress12U use ice_state, only: uvel, vvel, uvelE, vvelE, uvelN, vvelN, divu, shear use ice_grid, only: ULAT, NLAT, ELAT, tarea @@ -250,35 +246,8 @@ subroutine init_dyn (dt) if (grid_ice == 'CD' .or. grid_ice == 'C') then ! extra velocity variables uvelE (i,j,iblk) = c0 vvelE (i,j,iblk) = c0 - taubxE (i,j,iblk) = c0 - taubyE (i,j,iblk) = c0 - strairxE(i,j,iblk) = c0 - strairyE(i,j,iblk) = c0 - strocnxE(i,j,iblk) = c0 - strocnyE(i,j,iblk) = c0 - strtltxE(i,j,iblk) = c0 - strtltyE(i,j,iblk) = c0 - strintxE(i,j,iblk) = c0 - strintyE(i,j,iblk) = c0 - iceEmask(i,j,iblk) = c0 - fmE (i,j,iblk) = c0 - TbE (i,j,iblk) = c0 - uvelN (i,j,iblk) = c0 vvelN (i,j,iblk) = c0 - taubxN (i,j,iblk) = c0 - taubyN (i,j,iblk) = c0 - strairxN(i,j,iblk) = c0 - strairyN(i,j,iblk) = c0 - strocnxN(i,j,iblk) = c0 - strocnyN(i,j,iblk) = c0 - strtltxN(i,j,iblk) = c0 - strtltyN(i,j,iblk) = c0 - strintxN(i,j,iblk) = c0 - strintyN(i,j,iblk) = c0 - iceNmask(i,j,iblk) = c0 - fmN (i,j,iblk) = c0 - TbN (i,j,iblk) = c0 endif ! strain rates @@ -342,7 +311,10 @@ subroutine init_dyn (dt) ! ice extent mask on velocity points iceumask(i,j,iblk) = .false. - + if (grid_ice == 'CD' .or. grid_ice == 'C') then + iceemask(i,j,iblk) = .false. + icenmask(i,j,iblk) = .false. + end if enddo ! i enddo ! j enddo ! iblk diff --git a/cicecore/cicedynB/general/ice_flux.F90 b/cicecore/cicedynB/general/ice_flux.F90 index 312891f95..18727b63e 100644 --- a/cicecore/cicedynB/general/ice_flux.F90 +++ b/cicecore/cicedynB/general/ice_flux.F90 @@ -39,7 +39,7 @@ module ice_flux real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & - ! in from atmos (if .not.calc_strair) + ! in from atmos (if .not.calc_strair) strax , & ! wind stress components (N/m^2), on grid_atm_dynu stray , & ! on grid_atm_dynv @@ -48,7 +48,7 @@ module ice_flux vocn , & ! ocean current, y-direction (m/s), on grid_ocn_dynv ss_tltx , & ! sea surface slope, x-direction (m/m), on grid_ocn_dynu ss_tlty , & ! sea surface slope, y-direction, on grid_ocn_dynv - hwater , & ! water depth for seabed stress calc (landfast ice) + hwater , & ! water depth for seabed stress calc (landfast ice) ! out to atmosphere strairxT, & ! stress on ice by air, x-direction at T points, computed in icepack @@ -103,7 +103,7 @@ module ice_flux dvirdgdt, & ! rate of ice volume ridged (m/s) opening ! rate of opening due to divergence/shear (1/s) - real (kind=dbl_kind), & + real (kind=dbl_kind), & dimension (:,:,:,:), allocatable, public :: & ! ridging diagnostics in categories dardg1ndt, & ! rate of area loss by ridging ice (1/s) @@ -114,7 +114,7 @@ module ice_flux ardgn, & ! fractional area of ridged ice vrdgn, & ! volume of ridged ice araftn, & ! rafting ice area - vraftn, & ! rafting ice volume + vraftn, & ! rafting ice volume aredistn, & ! redistribution function: fraction of new ridge area vredistn ! redistribution function: fraction of new ridge volume @@ -178,7 +178,7 @@ module ice_flux ! NOTE: when in CICE_IN_NEMO mode, these are gridbox mean fields, ! not per ice area. When in standalone mode, these are per ice area. - real (kind=dbl_kind), & + real (kind=dbl_kind), & dimension (:,:,:,:), allocatable, public :: & fsurfn_f , & ! net flux to top surface, excluding fcondtop fcondtopn_f, & ! downward cond flux at top surface (W m-2) @@ -201,7 +201,7 @@ module ice_flux Tf , & ! freezing temperature (C) qdp , & ! deep ocean heat flux (W/m^2), negative upward hmix , & ! mixed layer depth (m) - daice_da ! data assimilation concentration increment rate + daice_da ! data assimilation concentration increment rate ! (concentration s-1)(only used in hadgem drivers) ! out to atmosphere (if calc_Tsfc) @@ -247,8 +247,8 @@ module ice_flux dimension(:,:,:,:), allocatable, public :: & albcnt ! counter for zenith angle - ! out to ocean - ! (Note CICE_IN_NEMO does not use these for coupling. + ! out to ocean + ! (Note CICE_IN_NEMO does not use these for coupling. ! It uses fresh_ai,fsalt_ai,fhocn_ai and fswthru_ai) real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & fpond , & ! fresh water flux to ponds (kg/m^2/s) @@ -280,7 +280,7 @@ module ice_flux snoicen ! snow-ice formation in category n (m) real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: & - keffn_top ! effective thermal conductivity of the top ice layer + keffn_top ! effective thermal conductivity of the top ice layer ! on categories (W/m^2/K) ! quantities passed from ocean mixed layer to atmosphere @@ -324,7 +324,7 @@ module ice_flux frz_onset, &! day of year that freezing begins (congel or frazil) frazil_diag ! frazil ice growth diagnostic (m/step-->cm/day) - real (kind=dbl_kind), & + real (kind=dbl_kind), & dimension (:,:,:,:), allocatable, public :: & fsurfn, & ! category fsurf fcondtopn,& ! category fcondtop @@ -339,7 +339,7 @@ module ice_flux ! As above but these remain grid box mean values i.e. they are not ! divided by aice at end of ice_dynamics. These are used in ! CICE_IN_NEMO for coupling and also for generating - ! ice diagnostics and history files as these are more accurate. + ! ice diagnostics and history files as these are more accurate. ! (The others suffer from problem of incorrect values at grid boxes ! that change from an ice free state to an icy state.) @@ -369,12 +369,12 @@ module ice_flux rside , & ! fraction of ice that melts laterally fside , & ! lateral heat flux (W/m^2) fsw , & ! incoming shortwave radiation (W/m^2) - coszen , & ! cosine solar zenith angle, < 0 for sun below horizon + coszen , & ! cosine solar zenith angle, < 0 for sun below horizon rdg_conv, & ! convergence term for ridging (1/s) rdg_shear ! shear term for ridging (1/s) real (kind=dbl_kind), dimension(:,:,:,:), allocatable, public :: & - salinz ,& ! initial salinity profile (ppt) + salinz ,& ! initial salinity profile (ppt) Tmltz ! initial melting temperature (^oC) !======================================================================= @@ -383,7 +383,7 @@ module ice_flux !======================================================================= ! -! Allocate space for all variables +! Allocate space for all variables ! subroutine alloc_flux @@ -393,12 +393,12 @@ subroutine alloc_flux allocate( & strax (nx_block,ny_block,max_blocks), & ! wind stress components (N/m^2) - stray (nx_block,ny_block,max_blocks), & ! + stray (nx_block,ny_block,max_blocks), & ! uocn (nx_block,ny_block,max_blocks), & ! ocean current, x-direction (m/s) vocn (nx_block,ny_block,max_blocks), & ! ocean current, y-direction (m/s) ss_tltx (nx_block,ny_block,max_blocks), & ! sea surface slope, x-direction (m/m) ss_tlty (nx_block,ny_block,max_blocks), & ! sea surface slope, y-direction - hwater (nx_block,ny_block,max_blocks), & ! water depth for seabed stress calc (landfast ice) + hwater (nx_block,ny_block,max_blocks), & ! water depth for seabed stress calc (landfast ice) strairxT (nx_block,ny_block,max_blocks), & ! stress on ice by air, x-direction strairyT (nx_block,ny_block,max_blocks), & ! stress on ice by air, y-direction strocnxT (nx_block,ny_block,max_blocks), & ! ice-ocean stress, x-direction @@ -544,7 +544,7 @@ subroutine alloc_flux rside (nx_block,ny_block,max_blocks), & ! fraction of ice that melts laterally fside (nx_block,ny_block,max_blocks), & ! lateral melt rate (W/m^2) fsw (nx_block,ny_block,max_blocks), & ! incoming shortwave radiation (W/m^2) - coszen (nx_block,ny_block,max_blocks), & ! cosine solar zenith angle, < 0 for sun below horizon + coszen (nx_block,ny_block,max_blocks), & ! cosine solar zenith angle, < 0 for sun below horizon rdg_conv (nx_block,ny_block,max_blocks), & ! convergence term for ridging (1/s) rdg_shear (nx_block,ny_block,max_blocks), & ! shear term for ridging (1/s) dardg1ndt (nx_block,ny_block,ncat,max_blocks), & ! rate of area loss by ridging ice (1/s) @@ -555,7 +555,7 @@ subroutine alloc_flux ardgn (nx_block,ny_block,ncat,max_blocks), & ! fractional area of ridged ice vrdgn (nx_block,ny_block,ncat,max_blocks), & ! volume of ridged ice araftn (nx_block,ny_block,ncat,max_blocks), & ! rafting ice area - vraftn (nx_block,ny_block,ncat,max_blocks), & ! rafting ice volume + vraftn (nx_block,ny_block,ncat,max_blocks), & ! rafting ice volume aredistn (nx_block,ny_block,ncat,max_blocks), & ! redistribution function: fraction of new ridge area vredistn (nx_block,ny_block,ncat,max_blocks), & ! redistribution function: fraction of new ridge volume fsurfn_f (nx_block,ny_block,ncat,max_blocks), & ! net flux to top surface, excluding fcondtop @@ -575,7 +575,7 @@ subroutine alloc_flux flatn (nx_block,ny_block,ncat,max_blocks), & ! category latent heat flux albcnt (nx_block,ny_block,max_blocks,max_nstrm), & ! counter for zenith angle snwcnt (nx_block,ny_block,max_blocks,max_nstrm), & ! counter for snow - salinz (nx_block,ny_block,nilyr+1,max_blocks), & ! initial salinity profile (ppt) + salinz (nx_block,ny_block,nilyr+1,max_blocks), & ! initial salinity profile (ppt) Tmltz (nx_block,ny_block,nilyr+1,max_blocks), & ! initial melting temperature (^oC) stat=ierr) if (ierr/=0) call abort_ice('(alloc_flux): Out of memory') @@ -719,7 +719,7 @@ subroutine init_coupler_flux fcondtopn_f(:,:,:,:) = c0 ! conductive heat flux (W/m^2) flatn_f (:,:,:,:) = -1.0_dbl_kind ! latent heat flux (W/m^2) fsensn_f (:,:,:,:) = c0 ! sensible heat flux (W/m^2) - endif ! + endif ! fiso_atm (:,:,:,:) = c0 ! isotope deposition rate (kg/m2/s) faero_atm (:,:,:,:) = c0 ! aerosol deposition rate (kg/m2/s) @@ -762,7 +762,7 @@ subroutine init_coupler_flux flat (:,:,:) = c0 fswabs (:,:,:) = c0 fswint_ai(:,:,:) = c0 - flwout (:,:,:) = -stefan_boltzmann*Tffresh**4 + flwout (:,:,:) = -stefan_boltzmann*Tffresh**4 ! in case atm model diagnoses Tsfc from flwout evap (:,:,:) = c0 evaps (:,:,:) = c0 @@ -816,7 +816,7 @@ subroutine init_coupler_flux coszen (:,:,:) = c0 ! Cosine of the zenith angle fsw (:,:,:) = c0 ! shortwave radiation (W/m^2) - scale_factor(:,:,:) = c1 ! shortwave scaling factor + scale_factor(:,:,:) = c1 ! shortwave scaling factor wind (:,:,:) = sqrt(uatm(:,:,:)**2 & + vatm(:,:,:)**2) ! wind speed, (m/s) Cdn_atm(:,:,:) = (vonkar/log(zref/iceruf)) & @@ -986,8 +986,8 @@ subroutine init_history_therm snowfrac (:,:,:) = c0 frazil_diag (:,:,:) = c0 - ! drag coefficients are computed prior to the atmo_boundary call, - ! during the thermodynamics section + ! drag coefficients are computed prior to the atmo_boundary call, + ! during the thermodynamics section Cdn_ocn(:,:,:) = dragio Cdn_atm(:,:,:) = (vonkar/log(zref/iceruf)) & * (vonkar/log(zref/iceruf)) ! atmo drag for RASM @@ -1023,6 +1023,7 @@ end subroutine init_history_therm subroutine init_history_dyn use ice_state, only: aice, vice, trcr, strength + use ice_grid, only: grid_ice logical (kind=log_kind) :: & tr_iage @@ -1074,6 +1075,32 @@ subroutine init_history_dyn aredistn (:,:,:,:) = c0 vredistn (:,:,:,:) = c0 + if (grid_ice == "CD" .or. grid_ice == "C") then + taubxE (:,:,:) = c0 + taubyE (:,:,:) = c0 + strocnxE (:,:,:) = c0 + strocnyE (:,:,:) = c0 + strairxE (:,:,:) = c0 + strairyE (:,:,:) = c0 + strtltxE (:,:,:) = c0 + strtltyE (:,:,:) = c0 + strintxE (:,:,:) = c0 + strintyE (:,:,:) = c0 + fmE (:,:,:) = c0 + TbE (:,:,:) = c0 + taubxN (:,:,:) = c0 + taubyN (:,:,:) = c0 + strocnxN (:,:,:) = c0 + strocnyN (:,:,:) = c0 + strairxN (:,:,:) = c0 + strairyN (:,:,:) = c0 + strtltxN (:,:,:) = c0 + strtltyN (:,:,:) = c0 + strintxN (:,:,:) = c0 + strintyN (:,:,:) = c0 + fmN (:,:,:) = c0 + TbN (:,:,:) = c0 + end if end subroutine init_history_dyn !======================================================================= @@ -1166,8 +1193,8 @@ subroutine scale_fluxes (nx_block, ny_block, & ! zsalinity fluxes real (kind=dbl_kind), dimension(nx_block,ny_block), intent(inout) :: & - fzsal , & ! salt flux to ocean with prognositic salinity (kg/m2/s) - fzsal_g ! Gravity drainage salt flux to ocean (kg/m2/s) + fzsal , & ! salt flux to ocean with prognositic salinity (kg/m2/s) + fzsal_g ! Gravity drainage salt flux to ocean (kg/m2/s) ! isotopes real (kind=dbl_kind), dimension(nx_block,ny_block,icepack_max_iso), & @@ -1221,8 +1248,8 @@ subroutine scale_fluxes (nx_block, ny_block, & alidr (i,j) = alidr (i,j) * ar alvdf (i,j) = alvdf (i,j) * ar alidf (i,j) = alidf (i,j) * ar - fzsal (i,j) = fzsal (i,j) * ar - fzsal_g (i,j) = fzsal_g (i,j) * ar + fzsal (i,j) = fzsal (i,j) * ar + fzsal_g (i,j) = fzsal_g (i,j) * ar flux_bio (i,j,:) = flux_bio (i,j,:) * ar faero_ocn(i,j,:) = faero_ocn(i,j,:) * ar if (present(Qref_iso )) Qref_iso (i,j,:) = Qref_iso (i,j,:) * ar @@ -1251,10 +1278,10 @@ subroutine scale_fluxes (nx_block, ny_block, & fswthru_idf (i,j) = c0 alvdr (i,j) = c0 ! zero out albedo where ice is absent alidr (i,j) = c0 - alvdf (i,j) = c0 + alvdf (i,j) = c0 alidf (i,j) = c0 - fzsal (i,j) = c0 - fzsal_g (i,j) = c0 + fzsal (i,j) = c0 + fzsal_g (i,j) = c0 flux_bio (i,j,:) = c0 faero_ocn(i,j,:) = c0 if (present(Qref_iso )) Qref_iso (i,j,:) = c0 @@ -1265,7 +1292,7 @@ subroutine scale_fluxes (nx_block, ny_block, & enddo ! j ! Scale fluxes for history output - if (present(fsurf) .and. present(fcondtop) ) then + if (present(fsurf) .and. present(fcondtop) ) then do j = 1, ny_block do i = 1, nx_block