* In ./src/core_atmosphere/physics and ./src/core_atmosphere/physics/physics_wrf,

  updated the Noah land surface scheme to that from WRF-v4.5. note that to run
  this updated version of the Noah LSM, we also need to use the updated versions
  of GENPARM.TBL, LANDUSE.TBL, SOILPARM.TBL, and VEGPARM.TBL. These updated
  tables are obtained from the v8.0 tag in the MPAS-Data repository.

Author: ldfowler58

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_driver_lsm.F

@@ -165,6 +165,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)       )
@@ -184,9 +186,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)   )
@@ -253,6 +252,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       )
@@ -272,9 +273,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   )
@@ -317,8 +315,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_albedo_seaice,sfc_emibck,sfc_emiss, &
-                                           sfcrunoff,smstav,smstot,snotime,snopcx,sr,udrunoff,   &
-                                           z0,znt
+                                           sfcrunoff,smstav,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
@@ -365,6 +363,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              )
@@ -439,6 +439,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)
@@ -472,13 +474,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
 
@@ -838,18 +836,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               &

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
+!=================================================================================================================
+
+
+
+

src/core_atmosphere/physics/mpas_atmphys_landuse.F

@@ -210,7 +210,7 @@ subroutine landuse_init_forMPAS(dminfo,julday,mesh,configs,diag_physics,sfc_inpu
                                      therin(ic,is),scfx(ic,is),sfhc(ic,is)
        enddo
 !      do ic = 1, lucats
-!         call mpas_log_write('$i $r $r $r $r $r $r $r $r', intArgs=(/ic/), &
+!         call mpas_log_write('$i $r $r $r $r $r $r $r', intArgs=(/ic/), &
 !                realArgs=(/albd(ic,is),slmo(ic,is),sfem(ic,is),sfz0(ic,is), &
 !                therin(ic,is),scfx(ic,is),sfhc(ic,is)/))
 !      enddo

src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F

@@ -215,11 +215,13 @@ subroutine soil_veg_gen_parm(dminfo,mminlu,mminsl)
  character(len=*),intent(inout):: mminlu, mminsl
 
 !local variables:
- character*128:: mess , message
+ character*128:: mess,message
+ character*128:: astring
 
- integer,parameter:: open_ok = 0
+ integer,parameter:: open_ok  = 0
+ integer,parameter:: loop_max = 10
  integer:: lumatch,iindex,lc,num_slope
- integer:: read_unit, istat
+ integer:: istat,loop_count,read_unit
 
 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
 !ALBBCK: SFC albedo (in percentage)
@@ -261,8 +263,10 @@ subroutine soil_veg_gen_parm(dminfo,mminlu,mminsl)
                                 'failure opening VEGPARM.TBL')
 
     lumatch=0
+
+    loop_count = 0
+    read(read_unit,fmt='(A)',end=2002) astring
     find_lutype : do while (lumatch == 0)
-       read(read_unit,*,end=2002)
        read(read_unit,*,end=2002) lutype
        read(read_unit,*) lucats,iindex
 
@@ -272,10 +276,17 @@ subroutine soil_veg_gen_parm(dminfo,mminlu,mminsl)
           call physics_message(mess)
           lumatch=1
        else
+          loop_count = loop_count + 1
           call physics_message('    skipping over lutype = ' // trim ( lutype ))
-          do lc = 1, lucats+20
-             read(read_unit,*)
-          enddo
+
+          find_vegetation_parameter_flag: do
+             read(read_unit,fmt='(A)',end=2002) astring
+             if(astring(1:21) .eq. 'Vegetation Parameters') then
+                exit find_vegetation_parameter_flag
+             elseif(loop_count .ge. loop_max) then
+                call physics_error_fatal('too many loops in VEGPARM.TBL')
+             endif
+          enddo find_vegetation_parameter_flag
        endif
     enddo find_lutype
 
@@ -307,6 +318,7 @@ subroutine soil_veg_gen_parm(dminfo,mminlu,mminsl)
                      albedomintbl(lc),albedomaxtbl(lc),z0mintbl(lc),z0maxtbl(lc),ztopvtbl(lc), &
                      zbotvtbl(lc)
        enddo
+
        read (read_unit,*)
        read (read_unit,*)topt_data
        read (read_unit,*)
@@ -322,11 +334,27 @@ subroutine soil_veg_gen_parm(dminfo,mminlu,mminsl)
        read (read_unit,*)
        read (read_unit,*)
        read (read_unit,*)
-       read (read_unit,*)low_density_residential
+       read (read_unit,*)lcz_1
+       read (read_unit,*)
+       read (read_unit,*)lcz_2
+       read (read_unit,*)
+       read (read_unit,*)lcz_3
+       read (read_unit,*)
+       read (read_unit,*)lcz_4
+       read (read_unit,*)
+       read (read_unit,*)lcz_5
+       read (read_unit,*)
+       read (read_unit,*)lcz_6
+       read (read_unit,*)
+       read (read_unit,*)lcz_7
+       read (read_unit,*)
+       read (read_unit,*)lcz_8
+       read (read_unit,*)
+       read (read_unit,*)lcz_9
        read (read_unit,*)
-       read (read_unit,*)high_density_residential
+       read (read_unit,*)lcz_10
        read (read_unit,*)
-       read (read_unit,*)high_intensity_industrial
+       read (read_unit,*)lcz_11
     endif
 
     2002 continue
@@ -364,9 +392,17 @@ subroutine soil_veg_gen_parm(dminfo,mminlu,mminsl)
  DM_BCAST_REAL(rsmax_data)
  DM_BCAST_INTEGER(bare)
  DM_BCAST_INTEGER(natural)
- DM_BCAST_INTEGER(low_density_residential)
- DM_BCAST_INTEGER(high_density_residential)
- DM_BCAST_INTEGER(high_intensity_industrial)
+ DM_BCAST_INTEGER(lcz_1)
+ DM_BCAST_INTEGER(lcz_2)
+ DM_BCAST_INTEGER(lcz_3)
+ DM_BCAST_INTEGER(lcz_4)
+ DM_BCAST_INTEGER(lcz_5)
+ DM_BCAST_INTEGER(lcz_6)
+ DM_BCAST_INTEGER(lcz_7)
+ DM_BCAST_INTEGER(lcz_8)
+ DM_BCAST_INTEGER(lcz_9)
+ DM_BCAST_INTEGER(lcz_10)
+ DM_BCAST_INTEGER(lcz_11)
 
 !call mpas_log_write(' LUTYPE  = '//trim(lutype))
 !call mpas_log_write(' LUCATS  = $i',intArgs=(/lucats/))
@@ -379,16 +415,24 @@ subroutine soil_veg_gen_parm(dminfo,mminlu,mminsl)
 !call mpas_log_write(' RSMAX_DATA  = $r',realArgs=(/rsmax_data/))
 !call mpas_log_write(' BARE        = $i',intArgs=(/bare/))
 !call mpas_log_write(' NATURAL     = $i',intArgs=(/natural/))
-!call mpas_log_write(' LOW_DENSITY_RESIDENTIAL   = $i , intArgs=(/low_density_residential/))
-!call mpas_log_write(' HIGH_DENSITY_RESIDENTIAL  = $i , intArgs=(/high_density_residential/))
-!call mpas_log_write(' HIGH_DENSITY_INDUSTRIAL   = $i , intArgs=(/high_density_industrial/))
+!call mpas_log_write(' LCZ_1       = $i', intArgs=(/lcz_1/))
+!call mpas_log_write(' LCZ_2       = $i', intArgs=(/lcz_2/))
+!call mpas_log_write(' LCZ_3       = $i', intArgs=(/lcz_3/))
+!call mpas_log_write(' LCZ_4       = $i', intArgs=(/lcz_4/))
+!call mpas_log_write(' LCZ_5       = $i', intArgs=(/lcz_5/))
+!call mpas_log_write(' LCZ_6       = $i', intArgs=(/lcz_6/))
+!call mpas_log_write(' LCZ_7       = $i', intArgs=(/lcz_7/))
+!call mpas_log_write(' LCZ_8       = $i', intArgs=(/lcz_8/))
+!call mpas_log_write(' LCZ_9       = $i', intArgs=(/lcz_9/))
+!call mpas_log_write(' LCZ_10      = $i', intArgs=(/lcz_10/))
+!call mpas_log_write(' LCZ_11      = $i', intArgs=(/lcz_11/))
 !call mpas_log_write('')
 !do lc = 1, lucats
 !   call mpas_log_write('$i $r $r $r $r $r $r $r $r $r $r $r $r $r $r $r $r $r', intArgs=(/lc/),     &
 !                realArgs=(/shdtbl(lc),float(nrotbl(lc)),rstbl(lc),rgltbl(lc),hstbl(lc),snuptbl(lc), &
 !                maxalb(lc),laimintbl(lc),laimaxtbl(lc),emissmintbl(lc),emissmaxtbl(lc),             &
 !                albedomintbl(lc),albedomaxtbl(lc),z0mintbl(lc),z0maxtbl(lc),ztopvtbl(lc),           &
-!                zbottvtbl(lc)/))
+!                zbotvtbl(lc)/))
 !enddo
  call mpas_log_write('    end read VEGPARM.TBL')
 

src/core_atmosphere/physics/mpas_atmphys_vars.F

@@ -649,7 +649,7 @@ module mpas_atmphys_vars
 !=================================================================================================================
 
  logical,parameter:: &
-    ua_phys=.false.    !option to activate UA Noah changes: a different snow-cover physics in the land-surface
+    ua_phys = .false.  !option to activate UA Noah changes: a different snow-cover physics in the land-surface
                        !scheme. That option is not currently implemented in MPAS.
 
  integer,parameter:: &
@@ -660,6 +660,9 @@ module mpas_atmphys_vars
  integer,parameter:: &
     fasdas = 0         !for WRF surface data assimilation system (not used in MPAS).
 
+ integer,parameter:: &
+    nurb = 1           !generic dimension for all dimensions needed to run the urban physics.
+
  integer,public:: &
     sf_surface_physics !used to define the land surface scheme by a number instead of name. It
                        !is only needed in module_ra_rrtmg_sw.F to define the spectral surface
@@ -728,13 +731,6 @@ module mpas_atmphys_vars
     frc_urb_p,        &!urban fraction                                                                         [-]
     ust_urb_p          !urban u* in similarity theory                                                        [m/s]
 
-!.. arrays needed in the argument list in the call to the Noah LSM hydrology model: note that these arrays are
-!.. initialized to zero since we do not run a hydrology model:
- real(kind=RKIND),dimension(:,:),allocatable:: &
-    infxsrt_p,        &!timestep infiltration excess                                                          [mm]
-    sfcheadrt_p,      &!surface water detph                                                                   [mm]
-    soldrain_p         !soil column drainage                                                                  [mm]
-
 !=================================================================================================================
 !.. variables and arrays related to surface characteristics:
 !=================================================================================================================

src/core_atmosphere/physics/physics_wrf/Makefile

@@ -7,6 +7,7 @@ dummy:
 
 OBJS = \
 	libmassv.o                     \
+	module_bep_bem_helper.o        \
 	module_bl_gwdo.o               \
 	module_bl_mynn.o               \
 	module_bl_ysu.o                \
@@ -69,9 +70,11 @@ module_ra_rrtmg_sw.o: \
 	module_ra_rrtmg_vinterp.o
 
 module_sf_bep.o: \
+	module_bep_bem_helper.o \
 	module_sf_urban.o
 
 module_sf_bep_bem.o: \
+	module_bep_bem_helper.o \
 	module_sf_bem.o \
 	module_sf_urban.o
 

src/core_atmosphere/physics/physics_wrf/module_bep_bem_helper.F

@@ -0,0 +1,4 @@
+MODULE module_bep_bem_helper
+   integer, save :: nurbm  ! Maximum number of urban classes    
+CONTAINS
+END MODULE module_bep_bem_helper