Modify NUOPC cap to accept separate glc runoff fluxes (#288)

* Modify NUOPC cap to accept separate glc runoff fluxes
* (1/2) Add separate fluxes for glc runoff. (2/2) Add heat content fields for lrunoff_glc and frunoff_glc.
* fix merge bugs and add more glc runoff diags
* enable glc runoff flux only if use_glc_runoff is present
* add ALLOW_GLC_RUNOFF_DIAGNOSTICS param to control whether to allow separate glacier runoff fluxes.

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -737,6 +737,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
              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% lrunoff_glc (isc:iec,jsc:jec),     &
+             Ice_ocean_boundary% frunoff_glc (isc:iec,jsc:jec),     &
+             Ice_ocean_boundary% hrofl_glc (isc:iec,jsc:jec),       &
+             Ice_ocean_boundary% hrofi_glc (isc:iec,jsc:jec),       &
              source=0.0)
 
     ! Needed for MARBL
@@ -797,6 +801,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Sa_pslv"        , "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofl"      , "will provide") !-> liquid runoff
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi"      , "will provide") !-> ice runoff
+  if (cesm_coupled) then
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "Forr_rofl_glc"  , "will provide") !-> liquid glc runoff
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "Forr_rofi_glc"  , "will provide") !-> frozen glc runoff
+  endif
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Si_ifrac"       , "will provide") !-> ice fraction
   call fld_list_add(fldsToOcn_num, fldsToOcn, "So_duu10n"      , "will provide") !-> wind^2 at 10m
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_meltw"     , "will provide")
@@ -808,6 +816,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hcond"     , "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hrofl"     , "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hrofi"     , "will provide")
+  if (cesm_coupled) then
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hrofl_glc" , "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hrofi_glc" , "will provide")
+  endif
 
   if (Ice_ocean_boundary%ice_ncat > 0) then
     call fld_list_add(fldsToOcn_num, fldsToOcn, "Sf_afracr", "will provide")
@@ -2855,6 +2867,34 @@ end subroutine shr_log_setLogUnit
 !!     <td></td>
 !! </tr>
 !! <tr>
+!!     <td>Forr_rofl_glc</td>
+!!     <td>kg m-2 s-1</td>
+!!     <td>runoff</td>
+!!     <td>mass flux of liquid glc runoff</td>
+!!     <td></td>
+!! </tr>
+!! <tr>
+!!     <td>Forr_rofi_glc</td>
+!!     <td>kg m-2 s-1</td>
+!!     <td>runoff</td>
+!!     <td>mass flux of frozen glc runoff</td>
+!!     <td></td>
+!! </tr>
+!! <tr>
+!!     <td>Foxx_hrofi_glc</td>
+!!     <td>W m-2</td>
+!!     <td>hrofi_glc</td>
+!!     <td>heat content (enthalpy) of frozen glc runoff</td>
+!!     <td></td>
+!! </tr>
+!! <tr>
+!!     <td>Foxx_hrofl_glc</td>
+!!     <td>W m-2</td>
+!!     <td>hrofl_glc</td>
+!!     <td>heat content (enthalpy) of liquid glc runoff</td>
+!!     <td></td>
+!! </tr>
+!! <tr>
 !!     <td>Fioi_salt</td>
 !!     <td>kg m-2 s-1</td>
 !!     <td>salt_flux</td>

config_src/drivers/nuopc_cap/mom_cap_methods.F90

@@ -216,6 +216,22 @@ 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
 
+  ! liquid glc runoff
+  if ( associated(ice_ocean_boundary%lrunoff_glc) ) then
+    ice_ocean_boundary%lrunoff_glc (:,:) = 0._ESMF_KIND_R8
+    call state_getimport(importState, 'Forr_rofl_glc',  &
+         isc, iec, jsc, jec, ice_ocean_boundary%lrunoff_glc, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+  endif
+
+  ! frozen glc runoff
+  if ( associated(ice_ocean_boundary%frunoff_glc) ) then
+    ice_ocean_boundary%frunoff_glc (:,:) = 0._ESMF_KIND_R8
+    call state_getimport(importState, 'Forr_rofi_glc',  &
+         isc, iec, jsc, jec, ice_ocean_boundary%frunoff_glc, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+  endif
+
   !----
   ! Enthalpy terms
   !----
@@ -256,6 +272,23 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
   end if
 
+  !----
+  ! enthalpy from liquid glc runoff (hrofl_glc)
+  !----
+  if ( associated(ice_ocean_boundary%hrofl_glc) ) then
+    call state_getimport(importState, 'Foxx_hrofl_glc', isc, iec, jsc, jec, &
+         ice_ocean_boundary%hrofl_glc, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+  end if
+
+  !----
+  ! enthalpy from frozen glc runoff (hrofi_glc)
+  !----
+  if ( associated(ice_ocean_boundary%hrofi_glc) ) then
+    call state_getimport(importState, 'Foxx_hrofi_glc', isc, iec, jsc, jec, &
+         ice_ocean_boundary%hrofi_glc, areacor=med2mod_areacor, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+  end if
   !----
   ! enthalpy from evaporation (hevap)
   !----

config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90

@@ -164,6 +164,8 @@ module MOM_surface_forcing_nuopc
 type, public :: ice_ocean_boundary_type
   real, pointer, dimension(:,:) :: lrunoff           =>NULL() !< liquid runoff [kg/m2/s]
   real, pointer, dimension(:,:) :: frunoff           =>NULL() !< ice runoff [kg/m2/s]
+  real, pointer, dimension(:,:) :: lrunoff_glc       =>NULL() !< liquid glc runoff via rof [kg/m2/s]
+  real, pointer, dimension(:,:) :: frunoff_glc       =>NULL() !< frozen glc runoff via rof [kg/m2/s]
   real, pointer, dimension(:,:) :: u_flux            =>NULL() !< i-direction wind stress [Pa]
   real, pointer, dimension(:,:) :: v_flux            =>NULL() !< j-direction wind stress [Pa]
   real, pointer, dimension(:,:) :: t_flux            =>NULL() !< sensible heat flux [W/m2]
@@ -183,6 +185,8 @@ module MOM_surface_forcing_nuopc
   real, pointer, dimension(:,:) :: mass_berg         =>NULL() !< mass of icebergs(kg/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(:,:) :: hrofl_glc         =>NULL() !< heat content from liquid glc runoff [W/m2]
+  real, pointer, dimension(:,:) :: hrofi_glc         =>NULL() !< heat content from frozen glc 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]
@@ -494,6 +498,16 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
       fluxes%frunoff(i,j) = kg_m2_s_conversion * IOB%frunoff(i-i0,j-j0) * G%mask2dT(i,j)
     endif
 
+    ! add liquid glc runoff flux via rof
+    if (associated(IOB%lrunoff_glc)) then
+      fluxes%lrunoff_glc(i,j) = kg_m2_s_conversion * IOB%lrunoff_glc(i-i0,j-j0) * G%mask2dT(i,j)
+    endif
+
+    ! ice glc runoff flux via rof
+    if (associated(IOB%frunoff_glc)) then
+      fluxes%frunoff_glc(i,j) = kg_m2_s_conversion * IOB%frunoff_glc(i-i0,j-j0) * G%mask2dT(i,j)
+    endif
+
     if (associated(IOB%ustar_berg)) &
       fluxes%ustar_berg(i,j) = US%m_to_Z*US%T_to_s * IOB%ustar_berg(i-i0,j-j0) * G%mask2dT(i,j)
 
@@ -531,6 +545,13 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
       fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * &
           IOB%frunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion
     endif
+    ! notice minus sign since frunoff_glc is positive into the ocean
+    if (associated(IOB%frunoff_glc)) then
+      fluxes%latent(i,j)              = fluxes%latent(i,j) - &
+          IOB%frunoff_glc(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion
+      fluxes%latent_frunoff_glc_diag(i,j) = fluxes%latent_frunoff_glc_diag(i,j) - G%mask2dT(i,j) * &
+          IOB%frunoff_glc(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion
+    endif
     if (associated(IOB%q_flux)) then
       fluxes%latent(i,j)           = fluxes%latent(i,j) + &
           IOB%q_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_vapor
@@ -572,6 +593,12 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
 
        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)
+
+       if (associated(IOB%hrofl_glc)) &
+        fluxes%heat_content_lrunoff_glc(i,j) = US%W_m2_to_QRZ_T * IOB%hrofl_glc(i-i0,j-j0) * G%mask2dT(i,j)
+
+       if (associated(IOB%hrofi_glc)) &
+        fluxes%heat_content_frunoff_glc(i,j) = US%W_m2_to_QRZ_T * IOB%hrofi_glc(i-i0,j-j0) * G%mask2dT(i,j)
     endif
 
     ! sea ice fraction [nondim]
@@ -633,7 +660,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     do j=js,je ; do i=is,ie
       net_FW(i,j) = US%RZ_T_to_kg_m2s * &
                     (((fluxes%lprec(i,j)   + fluxes%fprec(i,j) + fluxes%seaice_melt(i,j)) + &
-                      (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + &
+                      (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j) + &
+                       fluxes%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j))) + &
                       (fluxes%evap(i,j)    + fluxes%vprec(i,j)) ) * US%L_to_m**2*G%areaT(i,j)
       net_FW2(i,j) = net_FW(i,j) / (US%L_to_m**2*G%areaT(i,j))
     enddo ; enddo
@@ -1133,7 +1161,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
   ! Local variables
   real :: utide  ! The RMS tidal velocity [Z T-1 ~> m s-1].
   type(directories)  :: dirs
-  logical            :: new_sim, iceberg_flux_diags, fix_ustar_gustless_bug
+  logical            :: new_sim, iceberg_flux_diags, glc_runoff_diags, fix_ustar_gustless_bug
   logical :: test_value  ! This is used to determine whether a logical parameter is being set explicitly.
   logical :: explicit_bug, explicit_fix ! These indicate which parameters are set explicitly.
   type(time_type)    :: Time_frc
@@ -1431,8 +1459,13 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
                  "If true, makes available diagnostics of fluxes from icebergs "//&
                  "as seen by MOM6.", default=.false.)
 
+  call get_param(param_file, mdl, "ALLOW_GLC_RUNOFF_DIAGNOSTICS", glc_runoff_diags, &
+                 "If true, makes available diagnostics of separate glacier runoff fluxes"//&
+                 "as seen by MOM6.", default=.false.)
+
   call register_forcing_type_diags(Time, diag, US, CS%use_temperature, CS%handles, &
-                                   use_berg_fluxes=iceberg_flux_diags, use_waves=use_waves, use_cfcs=CS%use_CFC)
+                                   use_berg_fluxes=iceberg_flux_diags, use_waves=use_waves, &
+                                   use_cfcs=CS%use_CFC, use_glc_runoff=glc_runoff_diags)
 
   call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, &
                  "If true, allows flux adjustments to specified via the "//&
@@ -1541,6 +1574,8 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
   chks = field_chksum( iobt%fprec          ) ; if (root) write(outunit,100) 'iobt%fprec          ', chks
   chks = field_chksum( iobt%lrunoff        ) ; if (root) write(outunit,100) 'iobt%lrunoff        ', chks
   chks = field_chksum( iobt%frunoff        ) ; if (root) write(outunit,100) 'iobt%frunoff        ', chks
+  chks = field_chksum( iobt%lrunoff_glc    ) ; if (root) write(outunit,100) 'iobt%lrunoff_glc    ', chks
+  chks = field_chksum( iobt%frunoff_glc    ) ; if (root) write(outunit,100) 'iobt%frunoff_glc    ', chks
   chks = field_chksum( iobt%p              ) ; if (root) write(outunit,100) 'iobt%p              ', chks
   if (associated(iobt%ice_fraction)) then
     chks = field_chksum( iobt%ice_fraction ) ; if (root) write(outunit,100) 'iobt%ice_fraction   ', chks
@@ -1631,6 +1666,12 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
   if (associated(iobt%hcond)) then
     chks = field_chksum( iobt%hcond  ) ; if (root) write(outunit,100) 'iobt%hcond      ', chks
   endif
+  if (associated(iobt%hrofl_glc)) then
+    chks = field_chksum( iobt%hrofl_glc  ) ; if (root) write(outunit,100) 'iobt%hrofl_glc      ', chks
+  endif
+  if (associated(iobt%hrofl_glc)) then
+    chks = field_chksum( iobt%hrofl_glc  ) ; if (root) write(outunit,100) 'iobt%hrofl_glc      ', chks
+  endif
 
 100 FORMAT("   CHECKSUM::",A20," = ",Z20)
 110 FORMAT("   CHECKSUM::",A30," = ",Z20)