Merge branch 'zhan391/precfactor' (PR #7361)

Dynamical adjustment on freshwater fluxes balance for E3SM runs

This pull request introduces a new module and function in the MPAS-Ocean
(MPAS-O) model to compute and apply a scaling factor on precipitation
fluxes. The goal of this feature is to help maintain freshwater balance
in ocean-only simulations by mitigating the spurious linear trend
observed in sea surface height (SSH) over time. The updates specifically
target uncoupled ocean configurations, including:
* G-case: ocean + sea ice with data atmosphere
* GP- or J-case: land + river + ocean + sea ice with data atmosphere
These configurations do not include a dynamically coupled atmosphere,
and thus lack the feedback mechanisms necessary for naturally
maintaining freshwater balance. This PR provides a workaround by
dynamically adjusting the precipitation flux.

[NML]
[BFB] - stealth feature, off by default

Author: jonbob

cime_config/tests.py

@@ -315,6 +315,7 @@
             "SMS_P480_Ld5.TL319_IcoswISC30E3r5.GMPAS-JRA1p5-DIB-PISMF-TMIX.mpaso-jra_1958",
             "PET_P480_Ld2.TL319_IcoswISC30E3r5.GMPAS-JRA1p5-DIB-PISMF-DSGR.mpaso-jra_1958",
             "PET_P480_Ld2.TL319_IcoswISC30E3r5.GMPAS-JRA1p5-DIB-PISMF-DSGR.mpaso-freshwater_tracers_jra_1958",
+            "SMS_P480_Ld5.TL319_IcoswISC30E3r5.GMPAS-JRA1p5.mpaso-precip_scaling_jra_1958",
             )
         },
 

components/mpas-ocean/bld/build-namelist

@@ -728,6 +728,17 @@ add_default($nl, 'config_use_sgr_opt_salt_prescribed');
 add_default($nl, 'config_sgr_temperature_prescribed');
 add_default($nl, 'config_sgr_salinity_prescribed');
 
+##################################
+# Namelist group: precip_scaling #
+##################################
+
+add_default($nl, 'config_use_precip_scaling');
+add_default($nl, 'config_precip_scaling_mode');
+add_default($nl, 'config_precip_scaling_initial_factor');
+add_default($nl, 'config_precip_scaling_constant_factor');
+add_default($nl, 'config_precip_scaling_history_years');
+add_default($nl, 'config_precip_scaling_write_to_logfile');
+
 ############################
 # Namelist group: coupling #
 ############################
@@ -1892,6 +1903,7 @@ my @groups = qw(run_modes
                 wave_coupling
                 gotm
                 forcing
+		precip_scaling
                 coupling
                 shortwaveradiation
                 self_attraction_loading

components/mpas-ocean/bld/build-namelist-group-list

@@ -15,6 +15,7 @@ my @groups = qw(run_modes
                 wave_coupling
                 gotm
                 forcing
+                precip_scaling
                 coupling
                 shortwaveradiation
                 self_attraction_loading

components/mpas-ocean/bld/build-namelist-section

@@ -234,6 +234,17 @@ add_default($nl, 'config_use_sgr_opt_salt_prescribed');
 add_default($nl, 'config_sgr_temperature_prescribed');
 add_default($nl, 'config_sgr_salinity_prescribed');
 
+##################################
+# Namelist group: precip_scaling #
+##################################
+
+add_default($nl, 'config_use_precip_scaling');
+add_default($nl, 'config_precip_scaling_mode');
+add_default($nl, 'config_precip_scaling_initial_factor');
+add_default($nl, 'config_precip_scaling_constant_factor');
+add_default($nl, 'config_precip_scaling_history_years');
+add_default($nl, 'config_precip_scaling_write_to_logfile');
+
 ############################
 # Namelist group: coupling #
 ############################

components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml

@@ -372,6 +372,14 @@
 <config_sgr_temperature_prescribed>0.0</config_sgr_temperature_prescribed>
 <config_sgr_salinity_prescribed>0.0</config_sgr_salinity_prescribed>
 
+<!-- precip_scaling -->
+<config_use_precip_scaling>.false.</config_use_precip_scaling>
+<config_precip_scaling_mode>'off'</config_precip_scaling_mode>
+<config_precip_scaling_initial_factor>1.0</config_precip_scaling_initial_factor>
+<config_precip_scaling_constant_factor>1.0</config_precip_scaling_constant_factor>
+<config_precip_scaling_history_years>1</config_precip_scaling_history_years>
+<config_precip_scaling_write_to_logfile>.false.</config_precip_scaling_write_to_logfile>
+
 <!-- coupling -->
 <config_remove_ais_river_runoff>.false.</config_remove_ais_river_runoff>
 <config_remove_ais_ice_runoff>.false.</config_remove_ais_ice_runoff>

components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml

@@ -1261,6 +1261,57 @@ Default: Defined in namelist_defaults.xml
 </entry>
 
 
+<!-- precip_scaling -->
+
+<entry id="config_use_precip_scaling" type="logical"
+	category="precip_scaling" group="precip_scaling">
+If true, precipitation fluxes will be scalled to ensure freshwater balance in FOSI runs referring to forced ocean-sea ice simulation where ocean model is driven by a data atmosphere with prescribed precipitation (rain + snow) fluxes.
+
+Valid values: .true. or .false.
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_precip_scaling_mode" type="char*1024"
+	category="precip_scaling" group="precip_scaling">
+Option for scaling precipitation fluxes to ensure freshwater balance in FOSI runs. config_precip_scaling_mode='constant' uses a prescribed constant factor applied to precipitation fluxes. config_precip_scaling_mode='time-dependent' enables the dynamical caclulation of scaling factors using annual and global mean precipitation (rain + snow) and sea surface height changes over one-year or multi-year period.
+
+Valid values: 'off', 'constant' or 'time-dependent'
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_precip_scaling_initial_factor" type="real"
+	category="precip_scaling" group="precip_scaling">
+Initial guess scaling factor on precipitation fluxes in FOSI runs when config_precip_scaling_mode = 'time-dependent'.
+
+Valid values: Non-negative number between 0 and 1.0. Default value of 1.0 means no scaling.
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_precip_scaling_constant_factor" type="real"
+	category="precip_scaling" group="precip_scaling">
+Prescribed constant scaling factor on precipitation fluxes in FOSI runs when config_precip_scaling_mode = 'constant'.
+
+Valid values: Non-negative number between 0 and 1.0. Default value of 1.0 means no scaling.
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_precip_scaling_history_years" type="integer"
+	category="precip_scaling" group="precip_scaling">
+The number of years over which scaling factor on precipitation fluxes is computed in FOSI runs. Default value is 1 year.
+
+Valid values: Positive integer greater than or equal to 1.
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_precip_scaling_write_to_logfile" type="logical"
+	category="precip_scaling" group="precip_scaling">
+Logical flag indicating whether information for precipitation flux scaling in FOSI runs will be written to the log file.
+
+Valid values: .true. or .false.
+Default: Defined in namelist_defaults.xml
+</entry>
+
+
 <!-- coupling -->
 
 <entry id="config_remove_ais_river_runoff" type="logical"

components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/precip_scaling_jra_1958/README

@@ -0,0 +1,10 @@
+This testdef is used to test a stealth feature in mpaso introduced by
+PR #7361. It simplies changes one mpaso namelist variable,
+    config_use_precip_scaling
+from its default value of .false. to .true.. It also sets associated
+variables,
+    config_precip_scaling_mode
+    config_precip_scaling_constant_factor
+for testing. This particular testdef also limmits the required JRA
+forcing data to a single year to reduce the number of files needed for
+testing.

components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/precip_scaling_jra_1958/shell_commands

@@ -0,0 +1,4 @@
+./xmlchange DATM_CLMNCEP_YR_START=1958
+./xmlchange DATM_CLMNCEP_YR_END=1958
+./xmlchange DROF_STRM_YR_START=1958
+./xmlchange DROF_STRM_YR_END=1958

components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/precip_scaling_jra_1958/user_nl_mpaso

@@ -0,0 +1,3 @@
+config_use_precip_scaling = .true.
+config_precip_scaling_mode = 'constant'
+config_precip_scaling_constant_factor = 0.95

components/mpas-ocean/driver/ocn_comp_mct.F

@@ -70,6 +70,7 @@ module ocn_comp_mct
    use ocn_submesoscale_eddies
    use ocn_eddy_parameterization_helpers
    use ocn_scaled_dismf
+   use ocn_scaled_sfwf
 !
 ! !PUBLIC MEMBER FUNCTIONS:
   implicit none
@@ -121,6 +122,8 @@ module ocn_comp_mct
    integer :: itimestep, &  ! time step number for MPAS
               ocn_cpl_dt    ! length of coupling interval in seconds - set by coupler/ESMF
 
+   real (kind=RKIND) :: precFactor 
+
 !=======================================================================
 
 contains
@@ -235,6 +238,12 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename )!{{{
     real (kind=RKIND), pointer :: &
        runningMeanRemovedIceRunoff ! the area integrated, running mean of removed ice runoff from the ocean
 
+    ! freshwater balance adjustment 
+    logical, pointer :: config_use_precip_scaling
+    character (len=StrKIND), pointer :: config_precip_scaling_mode
+    real (kind=RKIND), pointer :: SFWFScalingFactor ! scaling factor on freshwater fluxes 
+    real(kind=RKIND), pointer :: scalingFactor
+
 #ifdef HAVE_MOAB
     character*100 outfile, wopts
     integer :: ierrmb, numco, tagtype, tagindex, ent_type
@@ -833,6 +842,25 @@ end subroutine xml_stream_get_attributes
 
        ! initialize scaled data ice-shelf melt fluxes based on remove ice runoff
        call ocn_init_scaled_dismf(domain)
+
+       ! Initialize the precipitation scaling scheme 
+       call ocn_scaled_sfwf_init(domain)
+       ! Ensure a valid first guess at the initial time 
+       call mpas_pool_get_config(domain % configs, 'config_use_precip_scaling', config_use_precip_scaling)
+       if (config_use_precip_scaling) then
+          call mpas_pool_get_config(domain % configs, 'config_precip_scaling_mode', config_precip_scaling_mode)
+          ! independent of space so should be no need to loop over blocks
+          block_ptr => domain % blocklist
+          call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool)
+          call mpas_pool_get_array(forcingPool, "SFWFScalingFactor", SFWFScalingFactor)
+          if (trim(config_precip_scaling_mode) == 'time-dependent') then 
+            call mpas_pool_get_config(domain % configs, 'config_precip_scaling_initial_factor', scalingFactor)
+          else if (trim(config_precip_scaling_mode) == 'constant') then
+            call mpas_pool_get_config(domain % configs, 'config_precip_scaling_constant_factor', scalingFactor) 
+          end if
+          ! Broadcast the scalar to the field
+          SFWFScalingFactor = scalingFactor
+       end if
     end if
 
 !-----------------------------------------------------------------------
@@ -881,6 +909,10 @@ end subroutine xml_stream_get_attributes
 
       ! initialize scaled data ice-shelf melt fluxes based on remove ice runoff
        call ocn_init_scaled_dismf(domain)
+
+      ! initialize scaled data freshwater fluxes based on water balance relationship 
+       call ocn_scaled_sfwf_init(domain)
+
     end if
 
     call t_stopf ('mpaso_mct_init')
@@ -915,6 +947,16 @@ end subroutine xml_stream_get_attributes
     call init_glc_z_levels(domain, ierr_local)
     ierr = ior(ierr, ierr_local)
 
+    !configuration for freshwater conservation
+    call mpas_pool_get_config(domain % configs, 'config_use_precip_scaling', config_use_precip_scaling)
+    if (config_use_precip_scaling) then 
+      ! independent of space so should be no need to loop over blocks
+      block_ptr => domain % blocklist
+      call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool)
+      call mpas_pool_get_array(forcingPool, "SFWFScalingFactor", SFWFScalingFactor)
+      call seq_infodata_PutData(infodata, precip_fact=SFWFScalingFactor)
+    end if 
+
 !-----------------------------------------------------------------------
 !
 !   get initial state from driver
@@ -1038,6 +1080,10 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{
       real (kind=RKIND), pointer :: &
          runningMeanRemovedIceRunoff ! the area integrated, running mean of removed ice runoff from the ocean
 
+      ! freshwater balance adjustment 
+      logical, pointer :: config_use_precip_scaling
+      real (kind=RKIND), pointer :: SFWFScalingFactor ! scaling factor on freshwater fluxes 
+
 #ifdef HAVE_MOAB
 #ifdef MOABCOMP
     real(r8)                 :: difference
@@ -1286,6 +1332,9 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{
             call mpas_log_write('   Validating ocean state')
          endif
 
+        ! update scalling factors for freshwater fluxes 
+        call ocn_update_scaling_factor(domain, timeLevel=1)
+
          call ocn_validate_state(domain, timeLevel=1)
          if (debugOn) call mpas_log_write('   Completed validating ocean state')
 
@@ -1356,6 +1405,16 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{
          call seq_infodata_PutData(infodata, rmean_rmv_ice_runoff=runningMeanRemovedIceRunoff)
       end if
 
+      !configuration for freshwater conservation
+      call mpas_pool_get_config(domain % configs, 'config_use_precip_scaling', config_use_precip_scaling)
+      if (config_use_precip_scaling) then
+         ! independent of space so should be no need to loop over blocks
+         block_ptr => domain % blocklist
+         call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool)
+         call mpas_pool_get_array(forcingPool, "SFWFScalingFactor", SFWFScalingFactor)
+         call seq_infodata_PutData(infodata, precip_fact=SFWFScalingFactor)
+      end if
+
       ! Reset I/O logs
       call shr_file_setLogUnit (shrlogunit)
       call shr_file_setLogLevel(shrloglev)
@@ -1980,7 +2039,11 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
                                      ecosysSeaIceCoupling, &
                                      DMSSeaIceCoupling,    &
                                      MacroMoleculesSeaIceCoupling,      &
-                                     CFCAuxiliary
+                                     CFCAuxiliary, &
+                                     conservationCheckEnergyAMPool, &
+                                     conservationCheckMassAMPool, &
+                                     conservationCheckSaltAMPool, &
+                                     conservationCheckCarbonAMPool
 
    integer, pointer :: nCellsSolve
 
@@ -2092,6 +2155,9 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
 
    real (kind=RKIND) :: riverFactor
 
+   ! freshwater balance adjustment 
+   logical, pointer :: config_use_precip_scaling
+
 !-----------------------------------------------------------------------
 !
 !  zero out padded cells
@@ -2123,6 +2189,14 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
    call mpas_pool_get_config(domain % configs, 'config_remove_ais_ice_runoff', config_remove_ais_ice_runoff)
    call mpas_pool_get_config(domain % configs, 'config_cvmix_kpp_use_theory_wave', config_cvmix_kpp_use_theory_wave)
 
+   precFactor = 1.0_RKIND  ! Default to no scaling
+   !configuration for freshwater balance constrains
+   call mpas_pool_get_config(domain % configs, 'config_use_precip_scaling', config_use_precip_scaling)
+   if (config_use_precip_scaling) then
+      !extract scaling factors 
+      call seq_infodata_GetData(infodata, precip_fact=precFactor)
+   end if
+
    n = 0
    removedRiverRunoffFluxThisProc = 0.0_RKIND
    removedIceRunoffFluxThisProc = 0.0_RKIND
@@ -2379,7 +2453,7 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
         end if
 
         if ( rainFluxField % isActive ) then
-           rainFlux(i) = x2o_o % rAttr(index_x2o_Faxa_rain, n)
+           rainFlux(i) = x2o_o % rAttr(index_x2o_Faxa_rain, n) * precFactor
         end if
         if ( atmosphericPressureField % isActive ) then
            atmosphericPressure(i) = x2o_o % rAttr(index_x2o_Sa_pslv,   n)
@@ -3955,7 +4029,7 @@ subroutine ocn_import_moab(Eclock, errorCode)!{{{
            icebergHeatFlux(i) = x2o_om(n, index_x2o_Fioi_bergh)
         end if
         if ( snowFluxField % isActive ) then
-           snowFlux(i) = x2o_om(n, index_x2o_Faxa_snow)
+           snowFlux(i) = x2o_om(n, index_x2o_Faxa_snow) * precFactor
         end if
         if ( seaIceFreshWaterFluxField % isActive ) then
            seaIceFreshWaterFlux(i) = x2o_om(n, index_x2o_Fioi_meltw)