Merge branch 'MPAS-v8.0/sf_noah' into develop (PR #1049)

This merge updates the Noah land surface model to match what is available in the
WRF v4.5 release.

Specifically, this merge makes the following changes:

 * The variables swddir, swddni, and swddif from module_ra_rrtmg_sw.F are now
   output; these variables are only used for the urban option BEP+BEM and not
   the Noah LSM itself.

 * In the namelist.atmosphere file, the namelist option config_lsm_scheme now
   uses "sf_noah" rather than "noah" to select the Noah LSM.

 * All code modules needed to run the Noah LSM have been updated to match what
   is available in the WRF v4.5 release. As part of this update, various code
   used only by urban options has also been introduced.

Author: mgduda

src/core_atmosphere/Registry.xml

@@ -2669,6 +2669,15 @@
                 <var name="acswuptc" type="real" dimensions="nCells Time" units="W m^{-2}"
                      description="accumulated clear-sky upward top-of-atmosphere shortwave radiation flux"/>
 
+                <var name="swddir" type="real" dimensions="nCells Time" units="W m^{-2}"
+                     description="shortwave surface downward direct irradiance"/>
+
+                <var name="swddni" type="real" dimensions="nCells Time" units="W m^{-2}"
+                     description="shortwave surface downward direct normal irradiance"/>
+
+                <var name="swddif" type="real" dimensions="nCells Time" units="W m^{-2}"
+                     description="shortwave surface downward diffuse irradiance"/>
+
                 <var name="gsw" type="real" dimensions="nCells Time" units="W m^{-2}"
                      description="net surface shortwave radiation flux"/>
 

src/core_atmosphere/physics/checkout_data_files.sh

@@ -23,7 +23,7 @@
 ################################################################################
 
 
-mpas_vers="7.0"
+mpas_vers="8.0"
 
 github_org="MPAS-Dev"   # GitHub organization where the MPAS-Data repository is found.
                         # For physics development, it can be helpful for a developer

src/core_atmosphere/physics/mpas_atmphys_control.F

@@ -74,6 +74,8 @@ module mpas_atmphys_control
 !   Laura D. Fowler (laura@ucar.edu) / 2016-12-30.
 ! * added the option sf_monin_obukhov_rev to run the revised surface layer scheme with the YSU PBL scheme.
 !   Laura D. Fowler (laura@ucar.edu) / 2023-05-15.
+! * replaced the option "noah" with "sf_noah" to run the NOAH land surface scheme.
+!   Laura D. Fowler (laura@ucar.edu) / 2022-02-18.
 
 
  contains
@@ -129,7 +131,7 @@ subroutine physics_namelist_check(configs)
     if (trim(config_radt_sw_scheme)    == 'suite') config_radt_sw_scheme    = 'rrtmg_sw'
     if (trim(config_radt_cld_scheme)   == 'suite') config_radt_cld_scheme   = 'cld_fraction'
     if (trim(config_sfclayer_scheme)   == 'suite') config_sfclayer_scheme   = 'sf_monin_obukhov'
-    if (trim(config_lsm_scheme)        == 'suite') config_lsm_scheme        = 'noah'
+    if (trim(config_lsm_scheme)        == 'suite') config_lsm_scheme        = 'sf_noah'
 
  else if (trim(config_physics_suite) == 'convection_permitting') then
 
@@ -141,7 +143,7 @@ subroutine physics_namelist_check(configs)
     if (trim(config_radt_sw_scheme)    == 'suite') config_radt_sw_scheme    = 'rrtmg_sw'
     if (trim(config_radt_cld_scheme)   == 'suite') config_radt_cld_scheme   = 'cld_fraction'
     if (trim(config_sfclayer_scheme)   == 'suite') config_sfclayer_scheme   = 'sf_mynn'
-    if (trim(config_lsm_scheme)        == 'suite') config_lsm_scheme        = 'noah'
+    if (trim(config_lsm_scheme)        == 'suite') config_lsm_scheme        = 'sf_noah'
 
  else if (trim(config_physics_suite) == 'none') then
 
@@ -286,7 +288,7 @@ subroutine physics_namelist_check(configs)
                              'set config_sfclayer_scheme different than off')
 
  elseif(.not. (config_lsm_scheme .eq. 'off ' .or. &
-               config_lsm_scheme .eq. 'noah')) then
+               config_lsm_scheme .eq. 'sf_noah')) then
 
     write(mpas_err_message,'(A,A10)') 'illegal value for land surface scheme: ', &
           trim(config_lsm_scheme)
@@ -357,7 +359,7 @@ subroutine physics_registry_init(mesh,configs,sfc_input)
 
     lsm_select: select case(trim(config_lsm_scheme))
 
-       case("noah")
+       case("sf_noah")
        !initialize the thickness of the soil layers for the Noah scheme:
           do iCell = 1, nCells
              dzs(1,iCell) = 0.10_RKIND

src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F

@@ -96,6 +96,8 @@ module mpas_atmphys_driver_lsm
 ! * the initialization of the variable albsi_p is switched from sfc_albedo_seaice (which is originally
 !   initialized to albbck to seaice_albedo_default. Note that albsi_p is not used if seaice_albedo_opt = 0.
 !   Laura D. Fowler (laura@ucar.edu) / 2020-05-10.
+! * replaced the option "noah" with "sf_noah" to run the NOAH land surface scheme.
+!   Laura D. Fowler (laura@ucar.edu) / 2022-02-18.
 
 
 !
@@ -166,6 +168,8 @@ subroutine allocate_lsm(config_frac_seaice)
  if(.not.allocated(snowc_p)     ) allocate(snowc_p(ims:ime,jms:jme)     )
  if(.not.allocated(snowh_p)     ) allocate(snowh_p(ims:ime,jms:jme)     )
  if(.not.allocated(sr_p)        ) allocate(sr_p(ims:ime,jms:jme)        )
+ if(.not.allocated(swddir_p)    ) allocate(swddir_p(ims:ime,jms:jme)    )
+ if(.not.allocated(swddif_p)    ) allocate(swddif_p(ims:ime,jms:jme)    )
  if(.not.allocated(swdown_p)    ) allocate(swdown_p(ims:ime,jms:jme)    )
  if(.not.allocated(tmn_p)       ) allocate(tmn_p(ims:ime,jms:jme)       )
  if(.not.allocated(tsk_p)       ) allocate(tsk_p(ims:ime,jms:jme)       )
@@ -185,9 +189,6 @@ subroutine allocate_lsm(config_frac_seaice)
  if(.not.allocated(frc_urb_p)   ) allocate(frc_urb_p(ims:ime,jms:jme)   )
  if(.not.allocated(ust_urb_p)   ) allocate(ust_urb_p(ims:ime,jms:jme)   )
  if(.not.allocated(utype_urb_p) ) allocate(utype_urb_p(ims:ime,jms:jme) )
- if(.not.allocated(infxsrt_p)   ) allocate(infxsrt_p(ims:ime,jms:jme)   )
- if(.not.allocated(sfcheadrt_p) ) allocate(sfcheadrt_p(ims:ime,jms:jme) )
- if(.not.allocated(soldrain_p)  ) allocate(soldrain_p(ims:ime,jms:jme)  )
 
  if(config_frac_seaice) then
     if(.not.allocated(tsk_sea)   ) allocate(tsk_sea(ims:ime,jms:jme)   )
@@ -254,6 +255,8 @@ subroutine deallocate_lsm(config_frac_seaice)
  if(allocated(snowc_p)     ) deallocate(snowc_p     )
  if(allocated(snowh_p)     ) deallocate(snowh_p     )
  if(allocated(sr_p)        ) deallocate(sr_p        )
+ if(allocated(swddir_p)    ) deallocate(swddir_p    )
+ if(allocated(swddif_p)    ) deallocate(swddif_p    )
  if(allocated(swdown_p)    ) deallocate(swdown_p    )
  if(allocated(tmn_p)       ) deallocate(tmn_p       )
  if(allocated(tsk_p)       ) deallocate(tsk_p       )
@@ -273,9 +276,6 @@ subroutine deallocate_lsm(config_frac_seaice)
  if(allocated(frc_urb_p)   ) deallocate(frc_urb_p   )
  if(allocated(ust_urb_p)   ) deallocate(ust_urb_p   )
  if(allocated(utype_urb_p) ) deallocate(utype_urb_p )
- if(allocated(infxsrt_p)   ) deallocate(infxsrt_p   )
- if(allocated(sfcheadrt_p) ) deallocate(sfcheadrt_p )
- if(allocated(soldrain_p)  ) deallocate(soldrain_p  )
 
  if(config_frac_seaice) then
     if(allocated(chs_sea)   ) deallocate(chs_sea   )
@@ -318,7 +318,8 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
  real(kind=RKIND),dimension(:),pointer  :: acsnom,acsnow,canwat,chs,chs2,chklowq,cpm,cqs2,glw,   &
                                            grdflx,gsw,hfx,lai,lh,noahres,potevp,qfx,qgh,qsfc,    &
                                            br,sfc_albedo,sfc_emibck,sfc_emiss,sfcrunoff,smstav,  &
-                                           smstot,snotime,snopcx,sr,udrunoff,z0,znt
+                                           smstot,snotime,snopcx,sr,swddif,swddir,udrunoff,      &
+                                           z0,znt
  real(kind=RKIND),dimension(:),pointer  :: shdmin,shdmax,snoalb,sfc_albbck,snow,snowc,snowh,tmn, &
                                            skintemp,vegfra,xice,xland
  real(kind=RKIND),dimension(:),pointer  :: t2m,th2m,q2
@@ -364,6 +365,8 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
  call mpas_pool_get_array(diag_physics,'smstot'           ,smstot           )
  call mpas_pool_get_array(diag_physics,'snotime'          ,snotime          )
  call mpas_pool_get_array(diag_physics,'snopcx'           ,snopcx           )
+ call mpas_pool_get_array(diag_physics,'swddif'           ,swddif           )
+ call mpas_pool_get_array(diag_physics,'swddir'           ,swddir           )
  call mpas_pool_get_array(diag_physics,'udrunoff'         ,udrunoff         )
  call mpas_pool_get_array(diag_physics,'z0'               ,z0               )
  call mpas_pool_get_array(diag_physics,'znt'              ,znt              )
@@ -438,6 +441,8 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
     smstot_p(i,j)     = smstot(i)
     snotime_p(i,j)    = snotime(i)
     snopcx_p(i,j)     = snopcx(i)
+    swddif_p(i,j)     = swddif(i)
+    swddir_p(i,j)     = swddir(i)
     udrunoff_p(i,j)   = udrunoff(i)
     z0_p(i,j)         = z0(i)
     znt_p(i,j)        = znt(i)
@@ -471,13 +476,9 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
     !initialization of arrays to run the Noah LSM urban parameterization (not currently
     frc_urb_p(i,j)    = 0._RKIND
     ust_urb_p(i,j)    = 0._RKIND
-    utype_urb_p(i,j)  = low_density_residential
+    utype_urb_p(i,j)  = 0
 
-    !initialization of arrays to run the Noah LSM hydrological parameterization (not currently
-    !implemented in MPAS):
-    infxsrt_p(i,j)    = 0._RKIND
-    sfcheadrt_p(i,j)  = 0._RKIND
-    soldrain_p(i,j)   = 0._RKIND
+    if(swddir_p(i,j).gt.0.) call mpas_log_write('--- sw: $i $r $r',intArgs=(/i/),realArgs=(/swddir_p(i,j),swddif_p(i,j)/))
  enddo
  enddo
 
@@ -762,7 +763,7 @@ subroutine init_lsm(dminfo,mesh,configs,diag_physics,sfc_input)
 
  lsm_select: select case (trim(lsm_scheme))
 
-    case ("noah")
+    case ("sf_noah")
        call noah_init_forMPAS(dminfo,mesh,configs,diag_physics,sfc_input)
 
     case default
@@ -809,8 +810,8 @@ subroutine driver_lsm(itimestep,configs,mesh,diag_physics,sfc_input,its,ite)
 !call to land-surface scheme:
  lsm_select: select case (trim(lsm_scheme))
 
-    case("noah")
-       call mpas_timer_start('Noah')
+    case("sf_noah")
+       call mpas_timer_start('sf_noah')
        call lsm( &
                 dz8w        = dz_p        , p8w3d       = pres2_hyd_p  , t3d       = t_p          , &
                 qv3d        = qv_p        , xland       = xland_p      , xice      = xice_p       , &
@@ -837,18 +838,21 @@ subroutine driver_lsm(itimestep,configs,mesh,diag_physics,sfc_input,its,ite)
                 opt_thcnd   = opt_thcnd   , ua_phys     = ua_phys      , flx4_2d   = flxsnow_p    , &
                 fvb_2d      = fvbsnow_p   , fbur_2d     = fbursnow_p   , fgsn_2d   = fgsnsnow_p   , &
                 utype_urb2d = utype_urb_p , frc_urb2d   = frc_urb_p    , ust_urb2d = ust_urb_p    , &
-                infxsrt     = infxsrt_p   , sfcheadrt   = sfcheadrt_p  , soldrain  = soldrain_p   , &
-                fasdas      = fasdas      , julian      = 0            , julyr     = 0            , &
+                swddir      = swddir_p    , swddif      = swddif_p     , fasdas    = fasdas       , &
+                julian      = 0           , julyr       = 0            ,                            &
+                num_soil_layers  = num_soils         ,                                              &
                 xice_threshold   = xice_threshold    ,                                              &
                 usemonalb        = config_sfc_albedo ,                                              &
-                mminlu           = mminlu ,                                                         &
-                num_soil_layers  = num_soils ,                                                      &
-                num_roof_layers  = num_soils ,                                                      &
-                num_wall_layers  = num_soils ,                                                      &
-                num_road_layers  = num_soils ,                                                      &
-                num_urban_layers = num_soils ,                                                      &
-                num_urban_hi     = num_soils ,                                                      &
-                sf_urban_physics = sf_urban_physics ,                                               &
+                mminlu           = mminlu            ,                                              &
+                sf_urban_physics = sf_urban_physics  ,                                              &
+                num_roof_layers  = nurb , num_wall_layers = nurb  ,                                 &
+                num_road_layers  = nurb , num_urban_hi    = nurb  ,                                 &
+                num_urban_ndm    = nurb , urban_map_zrd   = nurb  ,                                 &
+                urban_map_zwd    = nurb , urban_map_gd    = nurb  ,                                 &
+                urban_map_zd     = nurb , urban_map_zdf   = nurb  ,                                 &
+                urban_map_bd     = nurb , urban_map_wd    = nurb  ,                                 &
+                urban_map_gbd    = nurb , urban_map_fbd   = nurb  ,                                 &
+                urban_map_zgrd   = nurb ,                                                           &
                 ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde ,             &
                 ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme ,             &
                 its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte               &
@@ -895,7 +899,7 @@ subroutine driver_lsm(itimestep,configs,mesh,diag_physics,sfc_input,its,ite)
                 ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme ,           &
                 its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte             &
                     )
-       call mpas_timer_stop('Noah')
+       call mpas_timer_stop('sf_noah')
 
 
     case default

src/core_atmosphere/physics/mpas_atmphys_driver_lsm_shared.F

@@ -0,0 +1,65 @@
+! Copyright (c) 2013,  Los Alamos National Security, LLC (LANS)
+! and the University Corporation for Atmospheric Research (UCAR).
+!
+! Unless noted otherwise source code is licensed under the BSD license.
+! Additional copyright and license information can be found in the LICENSE file
+! distributed with this code, or at http://mpas-dev.github.com/license.html
+!
+!=================================================================================================================
+ module mpas_atmphys_driver_lsm_shared
+ use mpas_kind_types
+
+
+ implicit none
+ private
+ public:: correct_tsk_over_seaice
+
+
+ contains
+
+
+!=================================================================================================================
+ subroutine correct_tsk_over_seaice(ims,ime,jms,jme,its,ite,jts,jte,xice_thresh,xice,tsk,tsk_sea,tsk_ice)
+!=================================================================================================================
+
+!input arguments:
+ integer,intent(in):: ims,ime,its,ite,jms,jme,jts,jte
+ real(kind=RKIND),intent(in):: xice_thresh
+ real(kind=RKIND),intent(in),dimension(ims:ime,jms:jme):: tsk,xice
+
+!inout arguments:
+ real(kind=RKIND),intent(inout),dimension(ims:ime,jms:jme):: tsk_sea,tsk_ice
+
+!local variables:
+ integer:: i,j
+
+!-----------------------------------------------------------------------------------------------------------------
+
+!initialize the local sea-surface temperature and local sea-ice temperature to the local surface
+!temperature:
+ do j = jts,jte
+ do i = its,ite
+    tsk_sea(i,j) = tsk(i,j)
+    tsk_ice(i,j) = tsk(i,j)
+
+    if(xice(i,j).ge.xice_thresh .and. xice(i,j).le.1._RKIND) then
+       !over sea-ice grid cells, limit sea-surface temperatures to temperatures warmer than 271.4:
+       tsk_sea(i,j) = max(tsk_sea(i,j),271.4_RKIND)
+
+       !over sea-ice grid cells, avoids unphysically too cold sea-ice temperatures for grid cells
+       !with small sea-ice fractions:
+       if(xice(i,j).lt.0.2_RKIND .and. tsk_ice(i,j).lt.253.15_RKIND) tsk_ice(i,j) = 253.15_RKIND
+       if(xice(i,j).lt.0.1_RKIND .and. tsk_ice(i,j).lt.263.15_RKIND) tsk_ice(i,j) = 263.15_RKIND
+    endif
+ enddo
+ enddo
+
+ end subroutine correct_tsk_over_seaice
+
+!=================================================================================================================
+ end module mpas_atmphys_driver_lsm_shared
+!=================================================================================================================
+
+
+
+