From 471c010e4b56ca0c4f57011ac65afc096c212edc Mon Sep 17 00:00:00 2001 From: Denise Worthen Date: Thu, 23 Jun 2022 14:47:08 -0400 Subject: [PATCH] add Cgrid-related fixes for nuopc/cmeps (#728) * replace save_init with step_prep in CICE_RunMod * fixes for cgrid repro * remove added haloupdates * baselines pass with these extra halo updates removed * change F->S for ocean velocities and tilts * fix debug failure when grid_ice=C * compiling in debug mode using -init=snan,arrays requires initialization of variables * respond to review comments * remove inserted whitespace for uvelE,N and vvelE,N Co-authored-by: apcraig Co-authored-by: David Bailey Co-authored-by: Elizabeth Hunke Co-authored-by: Mariana Vertenstein Co-authored-by: Minsuk Ji <57227195+MinsukJi-NOAA@users.noreply.github.com> Co-authored-by: Tony Craig Co-authored-by: Philippe Blain --- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 7 +- cicecore/cicedynB/general/ice_flux.F90 | 91 ++++++++++++------- cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 | 2 +- .../drivers/nuopc/cmeps/ice_comp_nuopc.F90 | 1 - .../drivers/nuopc/cmeps/ice_import_export.F90 | 21 ++--- 5 files changed, 71 insertions(+), 51 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 722022fad..a30cc1b1c 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -192,7 +192,7 @@ 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, & @@ -311,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 diff --git a/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 index 6f145ab0e..79066e82a 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 @@ -129,7 +129,7 @@ subroutine ice_step use ice_restoring, only: restore_ice, ice_HaloRestore use ice_step_mod, only: prep_radiation, step_therm1, step_therm2, & update_state, step_dyn_horiz, step_dyn_ridge, step_radiation, & - biogeochemistry, save_init, step_dyn_wave, step_snow + biogeochemistry, step_prep, step_dyn_wave, step_snow use ice_timers, only: ice_timer_start, ice_timer_stop, & timer_diags, timer_column, timer_thermo, timer_bound, & timer_hist, timer_readwrite diff --git a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 index a9d71e479..8920ea386 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 @@ -549,7 +549,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) scol_valid = (scol_mask == 1) if (.not. scol_valid) then - write(6,*)'DEBUG: i am here' ! Advertise fields call ice_advertise_fields(gcomp, importState, exportState, flds_scalar_name, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index dbdf5c07d..7bfc53f45 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -187,9 +187,9 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam ! in the cmeps esmFldsExchange_xxx_mod.F90 that is model specific ! from atm - black carbon deposition fluxes (3) call fldlist_add(fldsToIce_num, fldsToIce, 'Faxa_bcph', ungridded_lbound=1, ungridded_ubound=3) - ! from atm - wet dust deposition frluxes (4 sizes) + ! from atm - wet dust deposition fluxes (4 sizes) call fldlist_add(fldsToIce_num, fldsToIce, 'Faxa_dstwet', ungridded_lbound=1, ungridded_ubound=4) - ! from - atm dry dust deposition frluxes (4 sizes) + ! from atm - dry dust deposition fluxes (4 sizes) call fldlist_add(fldsToIce_num, fldsToIce, 'Faxa_dstdry', ungridded_lbound=1, ungridded_ubound=4) do n = 1,fldsToIce_num @@ -800,19 +800,10 @@ subroutine ice_import( importState, rc ) if (.not.prescribed_ice) then call t_startf ('cice_imp_t2u') - call ice_HaloUpdate(uocn, halo_info, field_loc_center, field_type_scalar) - call ice_HaloUpdate(vocn, halo_info, field_loc_center, field_type_scalar) - call ice_HaloUpdate(ss_tltx, halo_info, field_loc_center, field_type_scalar) - call ice_HaloUpdate(ss_tlty, halo_info, field_loc_center, field_type_scalar) - ! tcraig, moved to dynamics for consistency - !work = uocn - !call grid_average_X2Y('F',work,'T',uocn,'U') - !work = vocn - !call grid_average_X2Y('F',work,'T',vocn,'U') - !work = ss_tltx - !call grid_average_X2Y('F',work,'T',ss_tltx,'U') - !work = ss_tlty - !call grid_average_X2Y('F',work,'T',ss_tlty,'U') + call ice_HaloUpdate(uocn, halo_info, field_loc_center, field_type_vector) + call ice_HaloUpdate(vocn, halo_info, field_loc_center, field_type_vector) + call ice_HaloUpdate(ss_tltx, halo_info, field_loc_center, field_type_vector) + call ice_HaloUpdate(ss_tlty, halo_info, field_loc_center, field_type_vector) call t_stopf ('cice_imp_t2u') end if