Add option to receive enthalpy fluxes via coupler

* Remove lrunoff_hflx and frunoff_hflx from the IOB type. These
are never used;

* mean_runoff_heat_flx and mean_calving_heat_flx are are never
advertized, so these have been deleted;

* If cesm_coupled=true, six new fields are imported:
  - heat_content_lprec
  - heat_content_fprec
  - heat_content_evap
  - heat_content_cond
  - heat_content_rofl
  - heat_content_rofi

* Add a new parameter (ENTHALPY_FROM_COUPLER) to control if the
enthalpy associated with mass entering/leaving the ocean is provided
via the coupler or calculated in MOM6.

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -654,8 +654,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   ocean_public%is_ocean_pe = .true.
   call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(restartfiles))
 
-  ! GMM, this call is not needed for NCAR. Check with EMC.
-  ! If this can be deleted, perhaps we should also delete ocean_model_flux_init
+  ! GMM, this call is not needed in CESM. Check with EMC if it can be deleted.
   call ocean_model_flux_init(ocean_state)
 
   call ocean_model_init_sfc(ocean_state, ocean_public)
@@ -680,9 +679,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
              Ice_ocean_boundary% ice_fraction (isc:iec,jsc:jec),    &
              Ice_ocean_boundary% u10_sqr (isc:iec,jsc:jec),         &
              Ice_ocean_boundary% p (isc:iec,jsc:jec),               &
-             Ice_ocean_boundary% lrunoff_hflx (isc:iec,jsc:jec),    &
-             Ice_ocean_boundary% frunoff_hflx (isc:iec,jsc:jec),    &
-             Ice_ocean_boundary% lrunoff (isc:iec,jsc:jec),       &
+             Ice_ocean_boundary% lrunoff (isc:iec,jsc:jec),         &
              Ice_ocean_boundary% frunoff (isc:iec,jsc:jec))
 
   Ice_ocean_boundary%u_flux          = 0.0
@@ -703,11 +700,25 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   Ice_ocean_boundary%ice_fraction    = 0.0
   Ice_ocean_boundary%u10_sqr         = 0.0
   Ice_ocean_boundary%p               = 0.0
-  Ice_ocean_boundary%lrunoff_hflx    = 0.0
-  Ice_ocean_boundary%frunoff_hflx    = 0.0
   Ice_ocean_boundary%lrunoff         = 0.0
   Ice_ocean_boundary%frunoff         = 0.0
 
+  if (cesm_coupled) then
+    allocate (Ice_ocean_boundary% hrain (isc:iec,jsc:jec),           &
+              Ice_ocean_boundary% hsnow (isc:iec,jsc:jec),           &
+              Ice_ocean_boundary% hrofl (isc:iec,jsc:jec),           &
+              Ice_ocean_boundary% hrofi (isc:iec,jsc:jec),           &
+              Ice_ocean_boundary% hevap (isc:iec,jsc:jec),           &
+              Ice_ocean_boundary% hcond (isc:iec,jsc:jec))
+
+    Ice_ocean_boundary%hrain           = 0.0
+    Ice_ocean_boundary%hsnow           = 0.0
+    Ice_ocean_boundary%hrofl           = 0.0
+    Ice_ocean_boundary%hrofi           = 0.0
+    Ice_ocean_boundary%hevap           = 0.0
+    Ice_ocean_boundary%hcond           = 0.0
+  endif
+
   call query_ocean_state(ocean_state, use_waves=use_waves, wave_method=wave_method)
   if (use_waves) then
     if (wave_method == "EFACTOR") then
@@ -756,9 +767,16 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fld_list_add(fldsToOcn_num, fldsToOcn, "So_duu10n"                  , "will provide") !-> wind^2 at 10m
   call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_fresh_water_to_ocean_rate", "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "net_heat_flx_to_ocn"        , "will provide")
-  !These are not currently used and changing requires a nuopc dictionary change
-  !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_heat_flx"        , "will provide")
-  !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx"       , "will provide")
+
+  if (cesm_coupled) then
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "heat_content_lprec", "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "heat_content_fprec", "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "heat_content_evap" , "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "heat_content_cond" , "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "heat_content_rofl" , "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "heat_content_rofi" , "will provide")
+  endif
+
   if (use_waves) then
     if (wave_method == "EFACTOR") then
       call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult"                 , "will provide")
@@ -1663,7 +1681,8 @@ subroutine ModelAdvance(gcomp, rc)
      !---------------
 
      if(profile_memory) call ESMF_VMLogMemInfo("Entering MOM update_ocean_model: ")
-     call update_ocean_model(Ice_ocean_boundary, ocean_state, ocean_public, Time, Time_step_coupled)
+     call update_ocean_model(Ice_ocean_boundary, ocean_state, ocean_public, Time, Time_step_coupled, &
+                             cesm_coupled)
      if(profile_memory) call ESMF_VMLogMemInfo("Leaving MOM update_ocean_model: ")
 
      !---------------
@@ -2388,7 +2407,7 @@ end subroutine shr_file_getLogUnit
 !! infrastructure when it's time for MOM to advance in time. During this subroutine, there is a
 !! call into the MOM update routine:
 !!
-!!      call update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_public, Time, Time_step_coupled)
+!!      call update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_public, Time, Time_step_coupled, cesm_coupled)
 !!
 !! Priori to the call to `update_ocean_model()`, the cap performs these steps
 !! - the `Time` and `Time_step_coupled` parameters, based on FMS types, are derived from the incoming ESMF clock
@@ -2499,13 +2518,6 @@ end subroutine shr_file_getLogUnit
 !!     <td></td>
 !! </tr>
 !! <tr>
-!!     <td>mean_calving_heat_flx</td>
-!!     <td>W m-2</td>
-!!     <td>calving_hflx</td>
-!!     <td>heat flux, relative to 0C, of frozen land water into ocean</td>
-!!     <td></td>
-!! </tr>
-!! <tr>
 !!     <td>mean_calving_rate</td>
 !!     <td>kg m-2 s-1</td>
 !!     <td>calving</td>
@@ -2576,10 +2588,45 @@ end subroutine shr_file_getLogUnit
 !!     <td></td>
 !! </tr>
 !! <tr>
-!!     <td>mean_runoff_heat_flx</td>
+!!     <td>heat_content_lprec</td>
+!!     <td>W m-2</td>
+!!     <td>hrain</td>
+!!     <td>heat content (enthalpy) of liquid water entering the ocean</td>
+!!     <td></td>
+!! </tr>
+!! <tr>
+!!     <td>heat_content_fprec</td>
+!!     <td>W m-2</td>
+!!     <td>hsnow</td>
+!!     <td>heat content (enthalpy) of frozen water entering the ocean</td>
+!!     <td></td>
+!! </tr>
+!! <tr>
+!!     <td>heat_content_evap</td>
+!!     <td>W m-2</td>
+!!     <td>hevap</td>
+!!     <td>heat content (enthalpy) of water leaving the ocean</td>
+!!     <td></td>
+!! </tr>
+!! <tr>
+!!     <td>heat_content_cond</td>
+!!     <td>W m-2</td>
+!!     <td>hcond</td>
+!!     <td>heat content (enthalpy) of liquid water entering the ocean due to condensation</td>
+!!     <td></td>
+!! </tr>
+!! <tr>
+!!     <td>heat_content_rofl</td>
+!!     <td>W m-2</td>
+!!     <td>hrofl</td>
+!!     <td>heat content (enthalpy) of liquid runoff</td>
+!!     <td></td>
+!! </tr>
+!! <tr>
+!!     <td>heat_content_rofi</td>
 !!     <td>W m-2</td>
-!!     <td>runoff_hflx</td>
-!!     <td>heat flux, relative to 0C, of liquid land water into ocean</td>
+!!     <td>hrofi</td>
+!!     <td>heat content (enthalpy) of frozen runoff</td>
 !!     <td></td>
 !! </tr>
 !! <tr>

config_src/drivers/nuopc_cap/mom_cap_methods.F90

@@ -217,17 +217,53 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
        isc, iec, jsc, jec, ice_ocean_boundary%frunoff, areacor=med2mod_areacor, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
-  ! heat content of lrunoff
-  ice_ocean_boundary%lrunoff_hflx(:,:) = 0._ESMF_KIND_R8
-  call state_getimport(importState, 'mean_runoff_heat_flx',  &
-       isc, iec, jsc, jec, ice_ocean_boundary%lrunoff_hflx, areacor=med2mod_areacor, rc=rc)
-  if (ChkErr(rc,__LINE__,u_FILE_u)) return
+  !----
+  ! Enthalpy terms (only in CESM)
+  !----
+  if (cesm_coupled) then
+    !----
+    ! enthalpy from liquid precipitation (hrain)
+    !----
+    call state_getimport(importState, 'heat_content_lprec', &
+         isc, iec, jsc, jec, ice_ocean_boundary%hrain, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    !----
+    ! enthalpy from frozen precipitation (hsnow)
+    !----
+    call state_getimport(importState, 'heat_content_fprec', &
+         isc, iec, jsc, jec, ice_ocean_boundary%hsnow, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    !----
+    ! enthalpy from liquid runoff (hrofl)
+    !----
+    call state_getimport(importState, 'heat_content_rofl', &
+         isc, iec, jsc, jec, ice_ocean_boundary%hrofl, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    !----
+    ! enthalpy from frozen runoff (hrofi)
+    !----
+    call state_getimport(importState, 'heat_content_rofi', &
+         isc, iec, jsc, jec, ice_ocean_boundary%hrofi, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    !----
+    ! enthalpy from evaporation (hevap)
+    !----
+    call state_getimport(importState, 'heat_content_evap', &
+         isc, iec, jsc, jec, ice_ocean_boundary%hevap, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    !----
+    ! enthalpy from condensation (hcond)
+    !----
+    call state_getimport(importState, 'heat_content_cond', &
+         isc, iec, jsc, jec, ice_ocean_boundary%hcond, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
-  ! heat content of frunoff
-  ice_ocean_boundary%frunoff_hflx(:,:) = 0._ESMF_KIND_R8
-  call state_getimport(importState, 'mean_calving_heat_flx',  &
-       isc, iec, jsc, jec, ice_ocean_boundary%frunoff_hflx, areacor=med2mod_areacor, rc=rc)
-  if (ChkErr(rc,__LINE__,u_FILE_u)) return
+  endif
 
   !----
   ! salt flux from ice

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -459,7 +459,7 @@ end subroutine ocean_model_init
 !! storing the new ocean properties in Ocean_state.
 subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
                               time_start_update, Ocean_coupling_time_step, &
-                              update_dyn, update_thermo, Ocn_fluxes_used)
+                              cesm_coupled, update_dyn, update_thermo, Ocn_fluxes_used)
   type(ice_ocean_boundary_type), &
                      intent(in)    :: Ice_ocean_boundary !< A structure containing the
                                               !! various forcing fields coming from the ice.
@@ -474,6 +474,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
   type(time_type),   intent(in)    :: time_start_update  !< The time at the beginning of the update step.
   type(time_type),   intent(in)    :: Ocean_coupling_time_step !< The amount of time over
                                               !! which to advance the ocean.
+  logical,           intent(in)    :: cesm_coupled !< Flag to check if coupled with cesm
   logical, optional, intent(in)    :: update_dyn !< If present and false, do not do updates
                                               !! due to the ocean dynamics.
   logical, optional, intent(in)    :: update_thermo !< If present and false, do not do updates
@@ -523,7 +524,6 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
 
   do_dyn = .true. ; if (present(update_dyn)) do_dyn = update_dyn
   do_thermo = .true. ; if (present(update_thermo)) do_thermo = update_thermo
-
   ! This is benign but not necessary if ocean_model_init_sfc was called or if
   ! OS%sfc_state%tr_fields was spawned in ocean_model_init.  Consider removing it.
   is = OS%grid%isc ; ie = OS%grid%iec ; js = OS%grid%jsc ; je = OS%grid%jec
@@ -690,7 +690,12 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
   call mech_forcing_diags(OS%forces, dt_coupling, OS%grid, OS%Time, OS%diag, OS%forcing_CSp%handles)
 
   if (OS%fluxes%fluxes_used) then
-    call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%grid, OS%US, OS%Time, OS%diag, OS%forcing_CSp%handles)
+    if (cesm_coupled) then
+      call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%grid, OS%US, OS%Time, OS%diag, &
+                               OS%forcing_CSp%handles, enthalpy=.true.)
+    else
+      call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%grid, OS%US, OS%Time, OS%diag, OS%forcing_CSp%handles)
+    endif
   endif
 
 ! Translate state into Ocean.

config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90

@@ -82,6 +82,8 @@ module MOM_surface_forcing_nuopc
                                 !! pressure limited by max_p_surf instead of the
                                 !! full atmospheric pressure.  The default is true.
   logical :: use_CFC            !< enables the MOM_CFC_cap tracer package.
+  logical :: enthalpy_cpl       !< Controls if enthalpy terms are provided by the coupler or computed
+                                !! internally.
   real :: gust_const            !< constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa]
   logical :: read_gust_2d       !< If true, use a 2-dimensional gustiness supplied
                                 !! from an input file.
@@ -181,8 +183,12 @@ module MOM_surface_forcing_nuopc
   real, pointer, dimension(:,:) :: ustar_berg        =>NULL() !< frictional velocity beneath icebergs [m/s]
   real, pointer, dimension(:,:) :: area_berg         =>NULL() !< area covered by icebergs[m2/m2]
   real, pointer, dimension(:,:) :: mass_berg         =>NULL() !< mass of icebergs(kg/m2)
-  real, pointer, dimension(:,:) :: lrunoff_hflx      =>NULL() !< heat content of liquid runoff [W/m2]
-  real, pointer, dimension(:,:) :: frunoff_hflx      =>NULL() !< heat content of frozen runoff [W/m2]
+  real, pointer, dimension(:,:) :: hrofl             =>NULL() !< heat content from liquid runoff [W/m2]
+  real, pointer, dimension(:,:) :: hrofi             =>NULL() !< heat content from frozen runoff [W/m2]
+  real, pointer, dimension(:,:) :: hrain             =>NULL() !< heat content from liquid precipitation [W/m2]
+  real, pointer, dimension(:,:) :: hsnow             =>NULL() !< heat content from frozen precipitation [W/m2]
+  real, pointer, dimension(:,:) :: hevap             =>NULL() !< heat content from evaporation [W/m2]
+  real, pointer, dimension(:,:) :: hcond             =>NULL() !< heat content from condensation [W/m2]
   real, pointer, dimension(:,:) :: p                 =>NULL() !< pressure of overlying ice and atmosphere
                                                               !< on ocean surface [Pa]
   real, pointer, dimension(:,:) :: ice_fraction      =>NULL() !< fractional ice area [nondim]
@@ -304,7 +310,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
   if (fluxes%dt_buoy_accum < 0) then
     call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., &
                                press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug, &
-                               cfc=CS%use_CFC)
+                               cfc=CS%use_CFC, hevap=CS%enthalpy_cpl)
 
     call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
@@ -487,13 +493,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     if (associated(IOB%mass_berg)) &
       fluxes%mass_berg(i,j) = US%m_to_Z*US%kg_m3_to_R * IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j)
 
-    if (associated(IOB%lrunoff_hflx)) &
-      fluxes%heat_content_lrunoff(i,j) = US%W_m2_to_QRZ_T * IOB%lrunoff_hflx(i-i0,j-j0) * G%mask2dT(i,j)
-
-    if (associated(IOB%frunoff_hflx)) &
-      fluxes%heat_content_frunoff(i,j) = US%W_m2_to_QRZ_T * kg_m2_s_conversion * &
-                                         IOB%frunoff_hflx(i-i0,j-j0) * G%mask2dT(i,j)
-
     if (associated(IOB%lw_flux)) &
       fluxes%lw(i,j) = US%W_m2_to_QRZ_T * IOB%lw_flux(i-i0,j-j0) * G%mask2dT(i,j)
 
@@ -544,6 +543,25 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + &
                      fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j)
 
+    ! enthalpy terms
+    if (associated(IOB%hrofl)) &
+      fluxes%heat_content_lrunoff(i,j) = US%W_m2_to_QRZ_T * IOB%hrofl(i-i0,j-j0) * G%mask2dT(i,j)
+
+    if (associated(IOB%hrofi)) &
+      fluxes%heat_content_frunoff(i,j) = US%W_m2_to_QRZ_T * IOB%hrofi(i-i0,j-j0) * G%mask2dT(i,j)
+
+    if (associated(IOB%hrain)) &
+      fluxes%heat_content_lprec(i,j)   = US%W_m2_to_QRZ_T * IOB%hrain(i-i0,j-j0) * G%mask2dT(i,j)
+
+    if (associated(IOB%hsnow)) &
+      fluxes%heat_content_fprec(i,j)   = US%W_m2_to_QRZ_T * IOB%hsnow(i-i0,j-j0) * G%mask2dT(i,j)
+
+     if (associated(IOB%hevap)) &
+      fluxes%heat_content_evap(i,j)    = US%W_m2_to_QRZ_T * IOB%hevap(i-i0,j-j0) * G%mask2dT(i,j)
+
+     if (associated(IOB%hcond)) &
+      fluxes%heat_content_cond(i,j)    = US%W_m2_to_QRZ_T * IOB%hcond(i-i0,j-j0) * G%mask2dT(i,j)
+
     ! sea ice fraction [nondim]
     if (associated(IOB%ice_fraction) .and. associated(fluxes%ice_fraction)) &
          fluxes%ice_fraction(i,j) = G%mask2dT(i,j) * IOB%ice_fraction(i-i0,j-j0)
@@ -1201,6 +1219,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
   call get_param(param_file, mdl, "USE_CFC_CAP", CS%use_CFC, &
                  default=.false., do_not_log=.true.)
 
+  call get_param(param_file, mdl, "ENTHALPY_FROM_COUPLER", CS%enthalpy_cpl, &
+                 default=.false., do_not_log=.true.)
+
   if (restore_salt) then
     call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, &
                  "The constant that relates the restoring surface fluxes to the relative "//&
@@ -1507,6 +1528,26 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
     chks = field_chksum( iobt%mass_berg  ) ; if (root) write(outunit,100) 'iobt%mass_berg      ', chks
   endif
 
+  ! enthalpy
+  if (associated(iobt%hrofl)) then
+    chks = field_chksum( iobt%hrofl  ) ; if (root) write(outunit,100) 'iobt%hrofl      ', chks
+  endif
+  if (associated(iobt%hrofi)) then
+    chks = field_chksum( iobt%hrofi  ) ; if (root) write(outunit,100) 'iobt%hrofi      ', chks
+  endif
+  if (associated(iobt%hrain)) then
+    chks = field_chksum( iobt%hrain  ) ; if (root) write(outunit,100) 'iobt%hrain      ', chks
+  endif
+  if (associated(iobt%hsnow)) then
+    chks = field_chksum( iobt%hsnow  ) ; if (root) write(outunit,100) 'iobt%hsnow      ', chks
+  endif
+  if (associated(iobt%hevap)) then
+    chks = field_chksum( iobt%hevap  ) ; if (root) write(outunit,100) 'iobt%hevap      ', chks
+  endif
+  if (associated(iobt%hcond)) then
+    chks = field_chksum( iobt%hcond  ) ; if (root) write(outunit,100) 'iobt%hcond      ', chks
+  endif
+
 100 FORMAT("   CHECKSUM::",A20," = ",Z20)
 
   call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%')