Skip to content

Add WW3 updates to CESM NUOPC driver #717

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
dabail10 wants to merge 8 commits into from
Closed

Add WW3 updates to CESM NUOPC driver #717

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
68b7634
Use internal CICE calendar parameters
dabail10 Apr 27, 2022
def5cdf
Revert "Use internal CICE calendar parameters"
dabail10 May 5, 2022
88b50af
Merge branch 'main' of https://github.com/ESCOMP/CICE into main
dabail10 May 11, 2022
eaa4bbd
NUOPC related updates for WW3 and leap years
dabail10 May 11, 2022
43e42de
Fix OMP PRIVATE of FSD variables and add tr_fsd if block
dabail10 May 11, 2022
d3a4d77
Replace t2u calls
dabail10 May 12, 2022
f99728e
Remove t2u/grid_average calls
dabail10 May 12, 2022
dbbb8e9
Additional fixes to the CAP to get it compiling again with Consortium…
dabail10 May 16, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 9 additions & 1 deletion cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ subroutine cice_init1()
use ice_grid , only: init_grid1, alloc_grid
use ice_domain , only: init_domain_blocks
use ice_arrays_column , only: alloc_arrays_column
use ice_arrays_column , only: wavefreq, dwavefreq, wave_spectrum
use ice_state , only: alloc_state
use ice_dyn_shared , only: alloc_dyn_shared
use ice_flux_bgc , only: alloc_flux_bgc
Expand Down Expand Up @@ -78,7 +79,7 @@ subroutine cice_init2()
use ice_calendar , only: dt, dt_dyn, istep, istep1, write_ic, init_calendar, calendar
use ice_communicate , only: my_task, master_task
use ice_diagnostics , only: init_diags
use ice_domain_size , only: ncat, nfsd
use ice_domain_size , only: ncat, nfsd, nfreq
use ice_dyn_eap , only: init_eap, alloc_dyn_eap
use ice_dyn_shared , only: kdyn, init_dyn
use ice_dyn_vp , only: init_vp
Expand All @@ -94,10 +95,12 @@ subroutine cice_init2()
use ice_restoring , only: ice_HaloRestore_init
use ice_timers , only: timer_total, init_ice_timers, ice_timer_start
use ice_transport_driver , only: init_transport
use ice_arrays_column , only: wavefreq, dwavefreq

logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers
logical(kind=log_kind) :: tr_iso, tr_fsd, wave_spec, tr_snow
character(len=char_len) :: snw_aging_table
real(kind=dbl_kind), dimension(25) :: wave_spectrum_profile ! hardwire for now
character(len=*), parameter :: subname = '(cice_init2)'
!----------------------------------------------------

Expand Down Expand Up @@ -177,6 +180,11 @@ subroutine cice_init2()
endif
endif

if (wave_spec) then
call icepack_init_wave(nfreq=nfreq, &
wave_spectrum_profile=wave_spectrum_profile, wavefreq=wavefreq, dwavefreq=dwavefreq)
end if

! Initialize shortwave components using swdn from previous timestep
! if restarting. These components will be scaled to current forcing
! in prep_radiation.
Expand Down
23 changes: 11 additions & 12 deletions cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ module ice_comp_nuopc
use ice_calendar , only : force_restart_now, write_ic
use ice_calendar , only : idate, mday, mmonth, myear, year_init
use ice_calendar , only : msec, dt, calendar, calendar_type, nextsw_cday, istep
use ice_calendar , only : ice_calendar_noleap, ice_calendar_gregorian
use ice_kinds_mod , only : dbl_kind, int_kind, char_len, char_len_long
use ice_fileunits , only : nu_diag, nu_diag_set, inst_index, inst_name
use ice_fileunits , only : inst_suffix, release_all_fileunits, flush_fileunit
Expand Down Expand Up @@ -80,9 +81,6 @@ module ice_comp_nuopc
character(len=*) , parameter :: orb_variable_year = 'variable_year'
character(len=*) , parameter :: orb_fixed_parameters = 'fixed_parameters'

character(len=*),parameter :: shr_cal_noleap = 'NO_LEAP'
character(len=*),parameter :: shr_cal_gregorian = 'GREGORIAN'

type(ESMF_Mesh) :: ice_mesh

integer :: nthrds ! Number of threads to use in this component
Expand Down Expand Up @@ -216,7 +214,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
type(ESMF_Time) :: stopTime ! Stop time
type(ESMF_Time) :: refTime ! Ref time
type(ESMF_TimeInterval) :: timeStep ! Model timestep
type(ESMF_Calendar) :: esmf_calendar ! esmf calendar
type(ESMF_CalKind_Flag) :: esmf_caltype ! esmf calendar type
integer :: start_ymd ! Start date (YYYYMMDD)
integer :: start_tod ! start time of day (s)
Expand Down Expand Up @@ -339,7 +336,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
call get_component_instance(gcomp, inst_suffix, inst_index, rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return

inst_name = "ICE"//trim(inst_suffix)
! inst_name = "ICE"//trim(inst_suffix)
inst_name = "ICE"

!----------------------------------------------------------------------------
! start cice timers
Expand Down Expand Up @@ -470,9 +468,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
call ESMF_TimeGet( currTime, calkindflag=esmf_caltype, rc=rc )
if (ChkErr(rc,__LINE__,u_FILE_u)) return
if (esmf_caltype == ESMF_CALKIND_NOLEAP) then
calendar_type = shr_cal_noleap
calendar_type = ice_calendar_noleap
else if (esmf_caltype == ESMF_CALKIND_GREGORIAN) then
calendar_type = shr_cal_gregorian
calendar_type = ice_calendar_gregorian
else
call abort_ice( subname//'ERROR:: bad calendar for ESMF' )
end if
Expand Down Expand Up @@ -549,7 +547,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
Expand Down Expand Up @@ -582,9 +579,11 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
end if
call icepack_query_parameters( tfrz_option_out=tfrz_option)
if (tfrz_option_driver /= tfrz_option) then
write(errmsg,'(a)') trim(subname)//'error: tfrz_option from driver '//trim(tfrz_option_driver)//&
' must be the same as tfrz_option from cice namelist '//trim(tfrz_option)
call abort_ice(trim(errmsg))
write(errmsg,'(a)') trim(subname)//'WARNING: tfrz_option from driver '//trim(tfrz_option_driver)//&
' is overwriting tfrz_option from cice namelist '//trim(tfrz_option)
write(nu_diag,*) trim(errmsg)
call icepack_warnings_flush(nu_diag)
call icepack_init_parameters(tfrz_option_in=tfrz_option_driver)
endif

! Flux convergence tolerance - always use the driver attribute value
Expand All @@ -595,7 +594,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
read(cvalue,*) atmiter_conv_driver
call icepack_query_parameters( atmiter_conv_out=atmiter_conv)
if (atmiter_conv_driver /= atmiter_conv) then
write(errmsg,'(a,d13.5,a,d13.5)') trim(subname)//'warning: atmiter_ from driver ',&
write(errmsg,'(a,d13.5,a,d13.5)') trim(subname)//'WARNING: atmiter_ from driver ',&
atmiter_conv_driver,' is overwritting atmiter_conv from cice namelist ',atmiter_conv
write(nu_diag,*) trim(errmsg)
call icepack_warnings_flush(nu_diag)
Expand Down
104 changes: 81 additions & 23 deletions cicecore/drivers/nuopc/cmeps/ice_import_export.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,15 @@ module ice_import_export
use ESMF
use NUOPC
use NUOPC_Model
use ice_kinds_mod , only : int_kind, dbl_kind, char_len, log_kind
use ice_kinds_mod , only : int_kind, dbl_kind, char_len, char_len_long, log_kind
use ice_constants , only : c0, c1, spval_dbl, radius
use ice_constants , only : field_loc_center, field_type_scalar, field_type_vector
use ice_blocks , only : block, get_block, nx_block, ny_block
use ice_domain , only : nblocks, blocks_ice, halo_info, distrb_info
use ice_domain_size , only : nx_global, ny_global, block_size_x, block_size_y, max_blocks, ncat
use ice_domain_size , only : nfreq, nfsd
use ice_exit , only : abort_ice
use ice_flux , only : strairxT, strairyT, strocnxT, strocnyT
use ice_flux , only : strairxt, strairyt, strocnxt, strocnyt
use ice_flux , only : alvdr, alidr, alvdf, alidf, Tref, Qref, Uref
use ice_flux , only : flat, fsens, flwout, evap, fswabs, fhocn, fswthru
use ice_flux , only : fswthru_vdr, fswthru_vdf, fswthru_idr, fswthru_idf
Expand All @@ -23,9 +24,10 @@ module ice_import_export
use ice_flux , only : fsnow, uocn, vocn, sst, ss_tltx, ss_tlty, frzmlt
use ice_flux , only : send_i2x_per_cat
use ice_flux , only : sss, Tf, wind, fsw
use ice_state , only : vice, vsno, aice, aicen_init, trcr
use ice_arrays_column , only : floe_rad_c, wave_spectrum
use ice_state , only : vice, vsno, aice, aicen_init, trcr, trcrn
use ice_grid , only : tlon, tlat, tarea, tmask, anglet, hm
use ice_grid , only : grid_type, grid_average_X2Y
use ice_grid , only : grid_type, t2ugrid_vector
use ice_mesh_mod , only : ocn_gridcell_frac
use ice_boundary , only : ice_HaloUpdate
use ice_fileunits , only : nu_diag, flush_fileunit
Expand All @@ -34,8 +36,10 @@ module ice_import_export
use ice_shr_methods , only : chkerr, state_reset
use icepack_intfc , only : icepack_warnings_flush, icepack_warnings_aborted
use icepack_intfc , only : icepack_query_parameters, icepack_query_tracer_flags
use icepack_intfc , only : icepack_query_tracer_indices
use icepack_intfc , only : icepack_liquidus_temperature
use icepack_intfc , only : icepack_sea_freezing_temperature
use icepack_parameters , only : puny, c2
use cice_wrapper_mod , only : t_startf, t_stopf, t_barrierf
#ifdef CESMCOUPLED
use shr_frz_mod , only : shr_frz_freezetemp
Expand Down Expand Up @@ -187,11 +191,16 @@ 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)

! the following are advertised but might not be connected if they are not advertised in the
! in the cmeps esmFldsExchange_xxx_mod.F90 that is model specific
! from wave
call fldlist_add(fldsToIce_num, fldsToIce, 'Sw_elevation_spectrum', ungridded_lbound=1, ungridded_ubound=25)

do n = 1,fldsToIce_num
call NUOPC_Advertise(importState, standardName=fldsToIce(n)%stdname, &
TransferOfferGeomObject='will provide', rc=rc)
Expand All @@ -214,6 +223,8 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam
call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_qref' )
call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_snowh' )
call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_u10' )
call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_thick' )
call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_floediam' )
call fldlist_add(fldsFrIce_num, fldsFrIce, 'inst_ice_vis_dir_albedo' )
call fldlist_add(fldsFrIce_num, fldsFrIce, 'inst_ice_ir_dir_albedo' )
call fldlist_add(fldsFrIce_num, fldsFrIce, 'inst_ice_vis_dif_albedo' )
Expand Down Expand Up @@ -292,7 +303,7 @@ subroutine ice_realize_fields(gcomp, mesh, flds_scalar_name, flds_scalar_num, rc
type(ESMF_State) :: exportState
type(ESMF_Field) :: lfield
integer :: numOwnedElements
integer :: i, j, iblk, n
integer :: i, j, iblk, n, k
integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain
type(block) :: this_block ! block information for current block
real(dbl_kind), allocatable :: mesh_areas(:)
Expand Down Expand Up @@ -403,11 +414,10 @@ subroutine ice_import( importState, rc )
! local variables
integer,parameter :: nflds=16
integer,parameter :: nfldv=6
integer :: i, j, iblk, n
integer :: i, j, iblk, n, k
integer :: ilo, ihi, jlo, jhi !beginning and end of physical domain
type(block) :: this_block ! block information for current block
real (kind=dbl_kind),allocatable :: aflds(:,:,:,:)
real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: work
real (kind=dbl_kind) :: workx, worky
real (kind=dbl_kind) :: MIN_RAIN_TEMP, MAX_SNOW_TEMP
real (kind=dbl_kind) :: Tffresh
Expand Down Expand Up @@ -559,6 +569,30 @@ subroutine ice_import( importState, rc )
end do
!$OMP END PARALLEL DO

! import wave elevation spectrum from wave (frequencies 1-25, assume that nfreq is 25)
if (State_FldChk(importState, 'Sw_elevation_spectrum')) then
if (nfreq /= 25) then
call abort_ice(trim(subname)//": ERROR nfreq not equal to 25 ")
end if
call state_getfldptr(importState, 'Sw_elevation_spectrum', fldptr=dataPtr2d, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return
do k = 1,nfreq
n = 0
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
do j = jlo, jhi
do i = ilo, ihi
n = n+1
wave_spectrum(i,j,k,iblk) = dataPtr2d(k,n)
end do
end do
end do
end do
end if


if ( State_fldChk(importState, 'Sa_ptem') .and. State_fldchk(importState,'air_density_height_lowest')) then
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
Expand Down Expand Up @@ -800,19 +834,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 t2ugrid_vector(uocn)

@DeniseWorthen DeniseWorthen May 12, 2022 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm getting a compile error for t2ugrid_vector, it has been removed from cicecore/cicedynB/infrastructure/ice_grid.F90 and replaced by the grid_averageX2Y interface it seems.

I'm not sure I understand the @apcraig comment; maybe this entire block should be removed.

@apcraig apcraig May 12, 2022

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

uocn, vocn, ss_tltx, ss_tlty are now assumed to be on the U grid and the averaging to the T grid is being done in CICE. It no longer needs to be done here. We are moving to a more general system in CICE where the location of the coupling/forcing variables are more clearly defined within CICE. This should be carefully reviewed in the context of CESM coupling to make sure it's working properly.

call t2ugrid_vector(vocn)
call t2ugrid_vector(ss_tltx)
call t2ugrid_vector(ss_tlty)
call t_stopf ('cice_imp_t2u')
end if

Expand Down Expand Up @@ -854,7 +879,7 @@ subroutine ice_export( exportState, rc )

! local variables
type(block) :: this_block ! block information for current block
integer :: i, j, iblk, n ! incides
integer :: i, j, iblk, n, k ! incides
integer :: n2 ! thickness category index
integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain
real (kind=dbl_kind) :: workx, worky ! tmps for converting grid
Expand All @@ -868,7 +893,10 @@ subroutine ice_export( exportState, rc )
real (kind=dbl_kind) :: tauxo (nx_block,ny_block,max_blocks) ! ice/ocean stress
real (kind=dbl_kind) :: tauyo (nx_block,ny_block,max_blocks) ! ice/ocean stress
real (kind=dbl_kind) :: ailohi(nx_block,ny_block,max_blocks) ! fractional ice area
real (kind=dbl_kind) :: floediam(nx_block,ny_block,max_blocks)
real (kind=dbl_kind) :: floethick(nx_block,ny_block,max_blocks) ! ice thickness
real (kind=dbl_kind) :: Tffresh
integer (kind=int_kind) :: nt_fsd
real (kind=dbl_kind), allocatable :: tempfld(:,:,:)
real (kind=dbl_kind), pointer :: dataptr_ifrac_n(:,:)
real (kind=dbl_kind), pointer :: dataptr_swpen_n(:,:)
Expand All @@ -886,6 +914,8 @@ subroutine ice_export( exportState, rc )
! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, &
! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit)

call icepack_query_tracer_indices(nt_fsd_out=nt_fsd)

call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=u_FILE_u, line=__LINE__)
Expand Down Expand Up @@ -913,6 +943,26 @@ subroutine ice_export( exportState, rc )
! ice fraction
ailohi(i,j,iblk) = min(aice(i,j,iblk), c1)

! ice thickness (m)

@DeniseWorthen DeniseWorthen May 11, 2022

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A comment about a preceding line (933) containing the OMP directive. Do additional variables need to be made private for the added floediam calculation?

@dabail10 dabail10 May 11, 2022

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch. Looks like k and floe_rad_c need to be added to the OMP PRIVATE statement. Anything else you noticed?

@DeniseWorthen DeniseWorthen May 11, 2022

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had the wave-ice coupling working in a branch, but I need to make sure it works if we don't have wave-coupling activated.

I also noticed the change at LN 802 in import export---is this related to the c-grid development?

@dabail10 dabail10 May 11, 2022

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes! The C-grid has been merged onto main.

@dabail10 dabail10 May 11, 2022

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, just found out that this code for the FSD needs to be in an if block. Probably the set_stateexport and so on needs to check for the fields being asked for as well.

if (aice(i,j,iblk) > puny) then
floethick(i,j,iblk) = vice(i,j,iblk) / aice(i,j,iblk)
else
floethick(i,j,iblk) = c0
end if

! floe diameter (m)
floediam(i,j,iblk) = c0
workx = c0
worky = c0
do n = 1, ncat
do k = 1, nfsd
workx = workx + floe_rad_c(k) * aicen_init(i,j,n,iblk) * trcrn(i,j,nt_fsd+k-1,n,iblk)
worky = worky + aicen_init(i,j,n,iblk) * trcrn(i,j,nt_fsd+k-1,n,iblk)
end do
end do
if (worky > c0) workx = c2*workx / worky
floediam(i,j,iblk) = MAX(c2*floe_rad_c(1),workx)

! surface temperature
Tsrf(i,j,iblk) = Tffresh + trcr(i,j,1,iblk) !Kelvin (original ???)

Expand Down Expand Up @@ -1037,6 +1087,14 @@ subroutine ice_export( exportState, rc )
call state_setexport(exportState, 'Si_qref' , input=Qref , lmask=tmask, ifrac=ailohi, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return

! Sea ice thickness (m)
call state_setexport(exportState, 'Si_thick' , input=floethick , lmask=tmask, ifrac=ailohi, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return

! Sea ice floe diameter (m)
call state_setexport(exportState, 'Si_floediam' , input=floediam , lmask=tmask, ifrac=ailohi, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return

! Snow volume
call state_setexport(exportState, 'mean_snow_volume' , input=vsno , lmask=tmask, ifrac=ailohi, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return
Expand Down