*Updates in FPMix and Stokes Most

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]
@@ -321,7 +325,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
                                press=.true., fix_accum_bug=.not.CS%ustar_gustless_bug, &
                                cfc=CS%use_CFC, marbl=CS%use_marbl_tracers, hevap=CS%enthalpy_cpl, &
                                tau_mag=.true., ice_ncat=IOB%ice_ncat)
-    !call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed)
+    call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_nir_dir,isd,ied,jsd,jed)
@@ -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
@@ -734,7 +762,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     call safe_alloc_ptr(forces%p_surf,isd,ied,jsd,jed)
     call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed)
-    !call safe_alloc_ptr(forces%omega_w2x,isd,ied,jsd,jed)
+    call safe_alloc_ptr(forces%omega_w2x,isd,ied,jsd,jed)
 
     if (CS%rigid_sea_ice) then
       call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed)
@@ -895,7 +923,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
       forces%tau_mag(i,j) = gustiness + G%mask2dT(i,j) * sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2)
       forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * &
                                sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2))
-      !forces%omega_w2x(i,j) = atan(tauy_at_h(i,j), taux_at_h(i,j))
+      forces%omega_w2x(i,j) = atan(tauy_at_h(i,j), taux_at_h(i,j))
     enddo ; enddo
     call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1)
   else ! C-grid wind stresses.
@@ -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)

src/ALE/MOM_regridding.F90

@@ -200,7 +200,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
   character(len=80)  :: string, string2, varName ! Temporary strings
   character(len=40)  :: coord_units, coord_res_param ! Temporary strings
   character(len=MAX_PARAM_LENGTH) :: param_name
-  character(len=200) :: inputdir, fileName
+  character(len=200) :: inputdir, fileName, longString
   character(len=320) :: message ! Temporary strings
   character(len=12) :: expected_units, alt_units ! Temporary strings
   logical :: tmpLogical, do_sum, main_parameters
@@ -680,7 +680,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
 
     ! Optionally specify maximum thicknesses for each layer, enforced by moving
     ! the interface below a layer downward.
-    call get_param(param_file, mdl, "MAX_LAYER_THICKNESS_CONFIG", string, &
+    call get_param(param_file, mdl, "MAX_LAYER_THICKNESS_CONFIG", longString, &
                    "Determines how to specify the maximum layer thicknesses.\n"//&
                    "Valid options are:\n"//&
                    " NONE        - there are no maximum layer thicknesses\n"//&
@@ -692,26 +692,26 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
                    default='NONE')
     message = "The list of maximum thickness for each layer."
     allocate(h_max(ke))
-    if ( trim(string) == "NONE") then
+    if ( trim(longString) == "NONE") then
       ! Do nothing.
-    elseif ( trim(string) ==  "PARAM") then
+    elseif ( trim(longString) ==  "PARAM") then
       call get_param(param_file, mdl, "MAX_LAYER_THICKNESS", h_max, &
                    trim(message), units="m", fail_if_missing=.true., scale=GV%m_to_H)
       call set_regrid_max_thickness(CS, h_max)
-    elseif (index(trim(string),'FILE:')==1) then
-      if (string(6:6)=='.' .or. string(6:6)=='/') then
+    elseif (index(trim(longString),'FILE:')==1) then
+      if (longString(6:6)=='.' .or. longString(6:6)=='/') then
         ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
-        fileName = trim( extractWord(trim(string(6:80)), 1) )
+        fileName = trim( extractWord(trim(longString(6:200)), 1) )
       else
         ! Otherwise assume we should look for the file in INPUTDIR
-        fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) )
+        fileName = trim(inputdir) // trim( extractWord(trim(longString(6:200)), 1) )
       endif
       if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
-        "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")")
+        "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(longString)//")")
 
-      varName = trim( extractWord(trim(string(6:)), 2) )
+      varName = trim( extractWord(trim(longString(6:)), 2) )
       if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
-        "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
+        "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(longString)//")")
       if (len_trim(varName)==0) then
         if (field_exists(fileName,'h_max')) then; varName = 'h_max'
         elseif (field_exists(fileName,'dz_max')) then; varName = 'dz_max'
@@ -723,14 +723,14 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
       call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
                  trim(message), units=coordinateUnits(coord_mode))
       call set_regrid_max_thickness(CS, h_max, GV%m_to_H)
-    elseif (index(trim(string),'FNC1:')==1) then
-      call dz_function1( trim(string(6:)), h_max )
+    elseif (index(trim(longString),'FNC1:')==1) then
+      call dz_function1( trim(longString(6:)), h_max )
       call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
                  trim(message), units=coordinateUnits(coord_mode))
       call set_regrid_max_thickness(CS, h_max, GV%m_to_H)
     else
       call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
-        "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string))
+        "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(longString))
     endif
     deallocate(h_max)
   endif