Add wave-ice coupling to nuopc/cmeps driver

cicecore/cicedynB/general/ice_state.F90

@@ -192,6 +192,12 @@ subroutine alloc_state
       n_trcr_strata = 0
       nt_strata = 0
       trcr_base = c0
+      aicen = c0
+      aicen_init = c0
+      vicen = c0
+      vicen_init = c0
+      vsnon = c0
+      vsnon_init = c0
 
       end subroutine alloc_state
 

cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90

@@ -78,7 +78,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
@@ -94,10 +94,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)'
     !----------------------------------------------------
 
@@ -177,6 +179,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.

cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90

@@ -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
@@ -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
@@ -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)
@@ -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
@@ -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
@@ -581,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
@@ -594,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)

cicecore/drivers/nuopc/cmeps/ice_import_export.F90

@@ -3,12 +3,13 @@ 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_iavg, strocnyT_iavg
   use ice_flux           , only : alvdr, alidr, alvdf, alidf, Tref, Qref, Uref
@@ -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
   use ice_mesh_mod       , only : ocn_gridcell_frac
   use ice_boundary       , only : ice_HaloUpdate
   use ice_fileunits      , only : nu_diag, flush_fileunit
@@ -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
@@ -112,6 +116,7 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam
     character(char_len) :: stdname
     character(char_len) :: cvalue
     logical             :: flds_wiso         ! use case
+    logical             :: flds_wave         ! use case
     logical             :: isPresent, isSet
     character(len=*), parameter :: subname='(ice_import_export:ice_advertise_fields)'
     !-------------------------------------------------------------------------------
@@ -148,6 +153,17 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam
        write(nu_diag,*)'flds_wiso = ',flds_wiso
     end if
 
+    flds_wave = .false.
+    call NUOPC_CompAttributeGet(gcomp, name='wav_coupling_to_cice', value=cvalue, &
+         isPresent=isPresent, isSet=isSet, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    if (isPresent .and. isSet) then
+       read(cvalue,*) flds_wave
+    end if
+    if (my_task == master_task) then
+       write(nu_diag,*)'flds_wave = ',flds_wave
+    end if
+
     !-----------------
     ! advertise import fields
     !-----------------
@@ -192,6 +208,14 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam
     ! 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
+    if (flds_wave) then
+       call fldlist_add(fldsToIce_num, fldsToIce, 'Sw_elevation_spectrum', ungridded_lbound=1, &
+            ungridded_ubound=25)
+    end if
+
     do n = 1,fldsToIce_num
        call NUOPC_Advertise(importState, standardName=fldsToIce(n)%stdname, &
             TransferOfferGeomObject='will provide', rc=rc)
@@ -225,6 +249,10 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam
        call fldlist_add(fldsFrIce_num, fldsFrIce, 'ice_fraction_n', &
             ungridded_lbound=1, ungridded_ubound=ncat)
     end if
+    if (flds_wave) then
+       call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_thick'                 )
+       call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_floediam'              )
+    end if
 
     ! ice/atm fluxes computed by ice
     call fldlist_add(fldsFrIce_num, fldsFrIce, 'stress_on_air_ice_zonal'          )
@@ -292,7 +320,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(:)
@@ -403,11 +431,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
@@ -559,6 +586,29 @@ 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
@@ -845,7 +895,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                     ! indices
     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
@@ -859,7 +909,11 @@ 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
+    logical (kind=log_kind) :: tr_fsd
+    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(:,:)
@@ -877,6 +931,9 @@ 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_query_tracer_flags(tr_fsd_out=tr_fsd)
+
     call icepack_warnings_flush(nu_diag)
     if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
         file=u_FILE_u, line=__LINE__)
@@ -890,8 +947,10 @@ subroutine ice_export( exportState, rc )
     tauya(:,:,:) = c0
     tauxo(:,:,:) = c0
     tauyo(:,:,:) = c0
+    floediam(:,:,:) = c0
+    floethick(:,:,:) = c0
 
-    !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky, this_block, ilo, ihi, jlo, jhi)
+    !$OMP PARALLEL DO PRIVATE(iblk,i,j,k,workx,worky, this_block, ilo, ihi, jlo, jhi)
     do iblk = 1, nblocks
        this_block = get_block(blocks_ice(iblk),iblk)
        ilo = this_block%ilo
@@ -904,6 +963,27 @@ subroutine ice_export( exportState, rc )
              ! ice fraction
              ailohi(i,j,iblk) = min(aice(i,j,iblk), c1)
 
+             if (tr_fsd) then
+                ! floe thickness (m)
+                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)
+                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)
+             endif
+
              ! surface temperature
              Tsrf(i,j,iblk)  = Tffresh + trcr(i,j,1,iblk)     !Kelvin (original ???)
 
@@ -1054,6 +1134,22 @@ subroutine ice_export( exportState, rc )
     call state_setexport(exportState, 'Si_snowh' , input=tempfld , lmask=tmask, ifrac=ailohi, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
+    ! ------
+    ! optional floe diameter and ice thickness to wave
+    ! ------
+
+    ! Sea ice thickness (m)
+    if (State_FldChk(exportState, 'Si_thick')) then
+       call state_setexport(exportState, 'Si_thick' , input=floethick , lmask=tmask, ifrac=ailohi, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    end if
+
+    ! Sea ice floe diameter (m)
+    if (State_FldChk(exportState, 'Si_floediam')) then
+       call state_setexport(exportState, 'Si_floediam' , input=floediam , lmask=tmask, ifrac=ailohi, rc=rc)
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    end if
+
     ! ------
     ! ice/atm fluxes computed by ice
     ! ------