Adjust the initialization of xice for "cold" sea-surface temperatures

Setting xice to 1.0 for sea-surface temperatures colder than 271 K has proved to
produce too much sea ice at high latitudes. Furthermore, many datasets provide a
sea ice fraction field. Determination of using actual sea ice fraction input or
setting xice to 1.0 for sea-surface temperatures colder than a pre-defined value
in the case where fractional sea ice is not provided in the input is now moved
to the subroutine physics_init_seaice in module mpas_atmphys_initialize_real.F.

  -> In mpas_init_atm_cases.F, remove the lines that set the variable
     xice equals to 1.0 when sea-surface temperatures are colder than 271K.
     This proves to lead too much seaice at high latitudes.

  -> In Registry.xml, added the parameter config_tsk_seaice_threshold that has
     a default value equals 100. This parameter is in a new &physics namelist
     group that is hidden by default.

  -> In mpas_atmphys_initialize_real.F, replaces the threshold 271K with the
     parameter config_tsk_seaice_threshold. This leads to a decrease in values
     of xice equals to 1 at high latitude. It follows the correction made in
     the WRF release v3.5.1.

Author: ldfowler58

src/core_atmosphere/physics/mpas_atmphys_initialize_real.F

@@ -55,6 +55,10 @@ module mpas_atmphys_initialize_real
 ! * In subroutine physics_init_seaice, added the initialization of the annual maximum snow albedo over seaice
 !   points to 0.75.
 ! * Laura D. Fowler (laura@ucar.edu) / 2022-03-15).
+! * In subroutine physics_init_seaice, corrected the initialization of seaice points when the surface
+!   temperature was originally colder than 271K. We now use the seaice threshold config_tsk_seaice_threshold
+!   that has a default value set to 100K. This leads to decreased seaice at high latitudes.
+!   Laura D. Fowler (laura@ucar.edu) / 2022-03-25.
 
 
  contains
@@ -599,6 +603,7 @@ subroutine physics_init_seaice(mesh, input, dims, configs)
 
  real(kind=RKIND):: xice_threshold
  real(kind=RKIND):: mid_point_depth
+ real(kind=RKIND),pointer:: tsk_seaice_threshold
  real(kind=RKIND),dimension(:),pointer  :: vegfra
  real(kind=RKIND),dimension(:),pointer  :: seaice,snoalb,xice
  real(kind=RKIND),dimension(:),pointer  :: skintemp,tmn,xland
@@ -611,15 +616,15 @@ subroutine physics_init_seaice(mesh, input, dims, configs)
 !note that this threshold is also defined in module_physics_vars.F.It is defined here to avoid
 !adding "use module_physics_vars" since this subroutine is only used for the initialization of
 !a "real" forecast with $CORE = init_nhyd_atmos.
- real(kind=RKIND),parameter:: xice_tsk_threshold = 271.
  real(kind=RKIND),parameter:: total_depth        = 3.   ! 3-meter soil depth.
 
 !-----------------------------------------------------------------------------------------------------------------
 !call mpas_log_write('')
 !call mpas_log_write('--- enter physics_init_seaice:')
 
- call mpas_pool_get_config(configs, 'config_frac_seaice', config_frac_seaice)
- call mpas_pool_get_config(configs, 'config_landuse_data', config_landuse_data)
+ call mpas_pool_get_config(configs, 'config_frac_seaice'         , config_frac_seaice)
+ call mpas_pool_get_config(configs, 'config_tsk_seaice_threshold', tsk_seaice_threshold)
+ call mpas_pool_get_config(configs, 'config_landuse_data'        , config_landuse_data)
 
  call mpas_pool_get_dimension(dims, 'nCellsSolve', nCellsSolve)
  call mpas_pool_get_dimension(dims, 'nSoilLevels', nSoilLevels)
@@ -660,13 +665,14 @@ subroutine physics_init_seaice(mesh, input, dims, configs)
  endif
  call mpas_log_write('--- config_frac_seaice      : $l', logicArgs=(/config_frac_seaice/))
  call mpas_log_write('--- xice_threshold          : $r', realArgs=(/xice_threshold/))
+ call mpas_log_write('--- tsk_seaice_threshold    : $r', realArgs=(/tsk_seaice_threshold/))
 
 !convert seaice points to land points when the sea-ice fraction is greater than the
 !prescribed threshold:
  num_seaice_changes = 0
  do iCell = 1, nCellsSolve
     if(xice(iCell) .ge. xice_threshold .or. &
-      (landmask(iCell).eq.0 .and. skintemp(iCell).lt.xice_tsk_threshold)) then
+      (landmask(iCell).eq.0 .and. skintemp(iCell).lt.tsk_seaice_threshold)) then
 
        num_seaice_changes = num_seaice_changes + 1
        !... sea-ice points are converted to land points:

src/core_init_atmosphere/Registry.xml

@@ -278,7 +278,12 @@
                      units="-"
                      description="Whether to switch sea-ice threshold from 0.5 to 0.02"
                      possible_values="true or false"/>
+        </nml_record>
 
+        <nml_record name="physics" in_defaults="false">
+                <nml_option name="config_tsk_seaice_threshold" type="real" default_value="100." units="K"
+                     description="surface temperature threshold below which water points are set to seaice points"
+                     possible_values="Positive real values"/>
         </nml_record>
 
         <nml_record name="io" in_defaults="true">

src/core_init_atmosphere/mpas_init_atm_cases.F

@@ -4559,9 +4559,6 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
 
       if (allocated(maskslab)) deallocate(maskslab)
 
-      ! Freeze really cold ocean
-      where (sst < 271.0 .and. landmask == 0) xice = 1.0
-
       ! Limit XICE to values between 0 and 1. Although the input meteorological field is between 0.
       ! and 1., interpolation to the MPAS grid can yield values of XiCE less than 0. and greater
       ! than 1.: