Add spatially heterogeneous soil colour maps for NoahMP (based on feature/monan-757-NF)

src/core_atmosphere/Registry.xml

@@ -508,8 +508,10 @@
 #endif
 #ifdef DO_PHYSICS
 			<var name="isltyp"/>
+			<var name="isctyp"/>
 			<var name="ivgtyp"/>
 			<var name="mminlu"/>
+			<var name="mminsc"/>
 			<var name="isice_lu"/>
 			<var name="iswater_lu"/>
 			<var name="landmask"/>
@@ -984,6 +986,7 @@
 			<var name="pin"/>
 			<var name="ozmixm"/>
 			<var name="mminlu"/>
+			<var name="mminsc"/>
 			<var name="isice_lu"/>
 			<var name="iswater_lu"/>
 			<var name="xland"/>
@@ -4462,14 +4465,20 @@
                 <!-- ================================================================================================== -->
 
                 <var name="isltyp" type="integer" dimensions="nCells" units="unitless"
-                     description="dominant soil category"/>
+                     description="dominant soil texture category"/>
+
+                <var name="isctyp" type="integer" dimensions="nCells" units="unitless"
+                     description="dominant soil colour category"/>
 
                 <var name="ivgtyp" type="integer" dimensions="nCells" units="unitless"
                      description="dominant vegetation category"/>
 
                 <var name="mminlu" type="text" dimensions="" units="unitless"
                      description="land use classification"/>
 
+                <var name="mminsc" type="text" dimensions="" units="unitless"
+                     description="soil colour classification"/>
+
                 <var name="isice_lu" name_in_code="isice" type="integer" dimensions="" units="unitless" default_value="24"
                      description="Index category for snow/ice"/>
 

src/core_atmosphere/mpas_atm_core.F

@@ -409,6 +409,7 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
       integer :: nCells, nEdges, nVertices, nVertLevels
       integer :: thread
       character(len=StrKIND), pointer :: mminlu
+      character(len=StrKIND), pointer :: mminsc
 
       integer, pointer :: nThreads
       integer, dimension(:), pointer :: cellThreadStart, cellThreadEnd
@@ -562,6 +563,16 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
             write(mminlu,'(a)') 'USGS'
       end if
 
+      ! Before calling physics_init, ensure that mminsc contains the name of the soil colour dataset
+      call mpas_pool_get_array(sfc_input, 'mminsc', mminsc)
+      if (len_trim(mminsc) == 0) then
+            call mpas_log_write('****************************************************************')
+            call mpas_log_write('No information on soil colour dataset is available.')
+            call mpas_log_write('Assume that we are using ''DEFAULT_RAD_NOAH''.')
+            call mpas_log_write('****************************************************************')
+            write(mminsc,'(a)') 'DEFAULT_RAD_NOAH'
+      end if
+
 
       if (moist_physics) then
          !initialization of some input variables in registry:

src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F

@@ -139,6 +139,7 @@ subroutine allocate_lsm
  if(.not.allocated(cpm_p)       ) allocate(cpm_p(ims:ime,jms:jme)       )
  if(.not.allocated(cqs2_p)      ) allocate(cqs2_p(ims:ime,jms:jme)      )
  if(.not.allocated(isltyp_p)    ) allocate(isltyp_p(ims:ime,jms:jme)    )
+ if(.not.allocated(isctyp_p)    ) allocate(isctyp_p(ims:ime,jms:jme)    )
  if(.not.allocated(ivgtyp_p)    ) allocate(ivgtyp_p(ims:ime,jms:jme)    )
  if(.not.allocated(glw_p)       ) allocate(glw_p(ims:ime,jms:jme)       )
  if(.not.allocated(grdflx_p)    ) allocate(grdflx_p(ims:ime,jms:jme)    )
@@ -216,6 +217,7 @@ subroutine deallocate_lsm
  if(allocated(gsw_p)       ) deallocate(gsw_p       )
  if(allocated(hfx_p)       ) deallocate(hfx_p       )
  if(allocated(isltyp_p)    ) deallocate(isltyp_p    )
+ if(allocated(isctyp_p)    ) deallocate(isctyp_p    )
  if(allocated(ivgtyp_p)    ) deallocate(ivgtyp_p    )
  if(allocated(lai_p)       ) deallocate(lai_p       )
  if(allocated(lh_p)        ) deallocate(lh_p        )
@@ -279,7 +281,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
  character(len=StrKIND),pointer:: config_microp_scheme,    &
                                   config_convection_scheme
 
- integer,dimension(:),pointer:: isltyp,ivgtyp
+ integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp
 
  real(kind=RKIND),dimension(:),pointer  :: acsnom,acsnow,canwat,chs,chs2,chklowq,cpm,cqs2,glw,   &
                                            grdflx,gsw,hfx,lai,lh,noahres,potevp,qfx,qgh,qsfc,    &
@@ -336,6 +338,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
  call mpas_pool_get_array(diag_physics,'znt'       ,znt       )
 
  call mpas_pool_get_array(sfc_input,'isltyp'    ,isltyp    )
+ call mpas_pool_get_array(sfc_input,'isctyp'    ,isctyp    )
  call mpas_pool_get_array(sfc_input,'ivgtyp'    ,ivgtyp    )
  call mpas_pool_get_array(sfc_input,'shdmin'    ,shdmin    )
  call mpas_pool_get_array(sfc_input,'shdmax'    ,shdmax    )
@@ -409,6 +412,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
     znt_p(i,j)        = znt(i)
 
     isltyp_p(i,j)     = isltyp(i)
+    isctyp_p(i,j)     = isctyp(i)
     ivgtyp_p(i,j)     = ivgtyp(i)
     shdmin_p(i,j)     = shdmin(i)
     shdmax_p(i,j)     = shdmax(i)
@@ -490,7 +494,7 @@ subroutine lsm_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
 !local pointers:
  character(len=StrKIND),pointer:: config_microp_scheme
 
- integer,dimension(:),pointer:: isltyp,ivgtyp
+ integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp
 
  real(kind=RKIND),dimension(:),pointer  :: acsnom,acsnow,canwat,chs,chs2,chklowq,cpm,cqs2,glw, &
                                            grdflx,gsw,hfx,lai,lh,noahres,potevp,qfx,qgh,qsfc,  &
@@ -544,6 +548,7 @@ subroutine lsm_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
  call mpas_pool_get_array(diag_physics,'znt'       ,znt       )
 
  call mpas_pool_get_array(sfc_input,'isltyp'    ,isltyp    )
+ call mpas_pool_get_array(sfc_input,'isctyp'    ,isctyp    )
  call mpas_pool_get_array(sfc_input,'ivgtyp'    ,ivgtyp    )
  call mpas_pool_get_array(sfc_input,'shdmin'    ,shdmin    )
  call mpas_pool_get_array(sfc_input,'shdmax'    ,shdmax    )

src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F

@@ -58,7 +58,7 @@ subroutine lsm_noahmp_fromMPAS(configs,mesh,diag,diag_physics,diag_physics_noahm
 
  integer:: i,its,ite
  integer:: n,ns,nsoil,nsnow,nzsnow
- integer,dimension(:),pointer:: isltyp,ivgtyp
+ integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp
 
  real(kind=RKIND),dimension(:),pointer:: latCell,lonCell
  real(kind=RKIND),dimension(:),pointer:: shdmax,shdmin,vegfra,tmn,xice,xland

src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F

@@ -91,7 +91,7 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input)
  character(len=StrKIND),pointer:: mminlu,mminsl
 
  integer,pointer:: nCells,nSoilLevels
- integer,dimension(:),pointer:: ivgtyp,isltyp
+ integer,dimension(:),pointer:: ivgtyp,isltyp,isctyp
 
  real(kind=RKIND),dimension(:),pointer:: snoalb,snow,snowh
  real(kind=RKIND),dimension(:,:),pointer:: tslb,smois,sh2o
@@ -118,6 +118,7 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input)
  call mpas_pool_get_dimension(mesh,'nSoilLevels',nSoilLevels)
 
  call mpas_pool_get_array(sfc_input,'isltyp', isltyp)
+ call mpas_pool_get_array(sfc_input,'isctyp', isctyp)
  call mpas_pool_get_array(sfc_input,'ivgtyp', ivgtyp)
  call mpas_pool_get_array(sfc_input,'sh2o'  , sh2o  )
  call mpas_pool_get_array(sfc_input,'smois' , smois )
@@ -149,6 +150,20 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input)
        call physics_error_fatal("module_sf_noahlsm.F: lsminit: out of range value "// &
                             "of ISLTYP. Is this field in the input?" )
 
+    ! Make sure all cells have reasonable soil colour indices.
+    errflag = 0
+    do iCell = 1, nCells
+       if(isctyp(iCell) < 1) then
+          errflag = 1
+          write(err_message,*) "module_sf_noahlsm.F: lsminit: out of range ISCTYP ", &
+                               iCell,isctyp(iCell),isltyp(iCell),ivgtyp(iCell)
+          call physics_message(err_message)
+       endif
+    end do
+    if (errflag == 1) &
+       call physics_error_fatal("module_sf_noahlsm.F: lsminit: out of range value "// &
+                                "of ISCTYP. Is this field in the input?" )
+
 !initializes soil liquid water content SH2O:
     do iCell = 1, nCells
        bx = bb(isltyp(iCell))

src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F

@@ -46,6 +46,7 @@ subroutine init_lsm_noahmp(configs,mesh,clock,diag_physics,diag_physics_noahmp,o
 
 !--- local variables and arrays:
  character(len=StrKIND),pointer:: mminlu
+ character(len=StrKIND),pointer:: mminsc
 
  integer:: ns
 
@@ -73,6 +74,10 @@ subroutine init_lsm_noahmp(configs,mesh,clock,diag_physics,diag_physics_noahmp,o
  call mpas_pool_get_array(sfc_input,'mminlu',mminlu)
  mpas_noahmp%llanduse = mminlu
 
+!--- read NoahmpTable.TBL:
+ call mpas_pool_get_array(sfc_input,'mminsc',mminsc)
+ mpas_noahmp%lsoilcol = mminsc
+
 !call mpas_log_write(' ')
 !call mpas_log_write('--- enter subroutine NoahmpReadTable:')
  call NoahmpReadTable(mpas_noahmp)
@@ -234,7 +239,7 @@ subroutine noahmp_init(configs,mesh,clock,diag_physics,diag_physics_noahmp,outpu
  logical,pointer:: urban_physics
 
  integer,pointer:: nsoilcomps
- integer,dimension(:),pointer:: isltyp,ivgtyp
+ integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp
  integer,dimension(:),pointer:: isnowxy
  integer,dimension(:),pointer:: irnumsi,irnummi,irnumfi
 
@@ -322,11 +327,14 @@ subroutine noahmp_init(configs,mesh,clock,diag_physics,diag_physics_noahmp,outpu
 !--- initialization of time-invariant surface variables needed in subroutine NoahmpInitMain:
  call mpas_pool_get_array(sfc_input,'dzs'   ,dzs   )
  call mpas_pool_get_array(sfc_input,'isltyp',isltyp)
+ call mpas_pool_get_array(sfc_input,'isctyp',isctyp)
  call mpas_pool_get_array(sfc_input,'ivgtyp',ivgtyp)
 
  do i = its, ite
     mpas_noahmp%isltyp(i) = isltyp(i)
+    mpas_noahmp%isctyp(i) = isctyp(i)
     mpas_noahmp%ivgtyp(i) = ivgtyp(i)
+
  enddo
  do ns = 1, nsoil
     mpas_noahmp%dzs(ns) = dzs(ns,its)

src/core_atmosphere/physics/mpas_atmphys_update_surface.F

@@ -151,7 +151,7 @@ subroutine physics_update_sst(dminfo,config_frac_seaice,mesh,sfc_input,diag_phys
 
  call mpas_pool_get_array(sfc_input,'isice'     ,isice     )
  call mpas_pool_get_array(sfc_input,'iswater'   ,iswater   )
- call mpas_pool_get_array(sfc_input,'isltyp'    ,isltyp    )    
+ call mpas_pool_get_array(sfc_input,'isltyp'    ,isltyp    )
  call mpas_pool_get_array(sfc_input,'ivgtyp'    ,ivgtyp    )
  call mpas_pool_get_array(sfc_input,'landmask'  ,landmask  )
  call mpas_pool_get_array(sfc_input,'vegfra'    ,vegfra    )

src/core_atmosphere/physics/mpas_atmphys_vars.F

@@ -882,6 +882,7 @@ module mpas_atmphys_vars
 
  integer,dimension(:,:),allocatable:: &
     isltyp_p,         &!dominant soil type category                                                            [-]
+    isctyp_p,         &!dominant soil colour category                                                            [-]
     ivgtyp_p           !dominant vegetation category                                                           [-]
 
  real(kind=RKIND),dimension(:),allocatable:: &

src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md

@@ -400,7 +400,7 @@
 
   - but in the file listed in the namelist as: HRLDAS_SETUP_FILE = "
   - The initialization fields are: SNOW,CANWAT,TSK,TSLB,SMOIS
-  - This file also contains the static grid/domain information: XLAT,XLONG,TMN,HGT,SEAICE,MAPFAC_MX,MAPFAC_MY,SHDMAX,SHDMIN,XLAND,IVGTYP,ISLTYP,DZS,ZS
+  - This file also contains the static grid/domain information: XLAT,XLONG,TMN,HGT,SEAICE,MAPFAC_MX,MAPFAC_MY,SHDMAX,SHDMIN,XLAND,IVGTYP,ISLTYP,ISCTYP,DZS,ZS
   - This file can also contains some optional fields: LAI  
   - NOTE: a WRF input file can be used as a HRLDAS_SETUP_FILE
 

src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90

@@ -62,7 +62,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO)
     ! config domain variable
     noahmp%config%domain%SurfaceType                 = 1
     noahmp%config%domain%NumSwRadBand                = 2
-    noahmp%config%domain%SoilColor                   = 4
+    noahmp%config%domain%SoilColor                   = NoahmpIO%ISCTYP(I)
     noahmp%config%domain%NumCropGrowStage            = 8
     noahmp%config%domain%FlagSoilProcess             = NoahmpIO%calculate_soil
     noahmp%config%domain%NumSoilTimeStep             = NoahmpIO%soil_update_steps
@@ -75,6 +75,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO)
     noahmp%config%domain%SoilTimeStep                = NoahmpIO%DTBL * NoahmpIO%soil_update_steps
     noahmp%config%domain%GridSize                    = NoahmpIO%DX
     noahmp%config%domain%LandUseDataName             = NoahmpIO%LLANDUSE
+    noahmp%config%domain%SoilColDataName             = NoahmpIO%LSOILCOL
     noahmp%config%domain%VegType                     = NoahmpIO%IVGTYP(I)
     noahmp%config%domain%CropType                    = NoahmpIO%CROPCAT(I)
     noahmp%config%domain%IndicatorIceSfc             = NoahmpIO%ICE

src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90

@@ -37,6 +37,7 @@ subroutine NoahmpIOVarFinalizeDefault(NoahmpIO)
     if ( allocated (NoahmpIO%zsoil)     ) deallocate ( NoahmpIO%zsoil               ) ! depth to soil interfaces [m]
     if ( allocated (NoahmpIO%ivgtyp)    ) deallocate ( NoahmpIO%ivgtyp              ) ! vegetation type
     if ( allocated (NoahmpIO%isltyp)    ) deallocate ( NoahmpIO%isltyp              ) ! soil type
+    if ( allocated (NoahmpIO%isctyp)    ) deallocate ( NoahmpIO%isctyp              ) ! soil colour class
     if ( allocated (NoahmpIO%vegfra)    ) deallocate ( NoahmpIO%vegfra              ) ! vegetation fraction []
     if ( allocated (NoahmpIO%tmn)       ) deallocate ( NoahmpIO%tmn                 ) ! deep soil temperature [K]
     if ( allocated (NoahmpIO%xland)     ) deallocate ( NoahmpIO%xland               ) ! =2 ocean; =1 land/seaice

src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90

@@ -41,6 +41,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
     if ( .not. allocated (NoahmpIO%zsoil)     ) allocate ( NoahmpIO%zsoil      (1:nsoil            )      ) ! depth to soil interfaces [m]
     if ( .not. allocated (NoahmpIO%ivgtyp)    ) allocate ( NoahmpIO%ivgtyp     (its:ite            )      ) ! vegetation type
     if ( .not. allocated (NoahmpIO%isltyp)    ) allocate ( NoahmpIO%isltyp     (its:ite            )      ) ! soil type
+    if ( .not. allocated (NoahmpIO%isctyp)    ) allocate ( NoahmpIO%isctyp     (its:ite            )      ) ! soil colour class
     if ( .not. allocated (NoahmpIO%vegfra)    ) allocate ( NoahmpIO%vegfra     (its:ite            )      ) ! vegetation fraction []
     if ( .not. allocated (NoahmpIO%tmn)       ) allocate ( NoahmpIO%tmn        (its:ite            )      ) ! deep soil temperature [K]
     if ( .not. allocated (NoahmpIO%xland)     ) allocate ( NoahmpIO%xland      (its:ite            )      ) ! =2 ocean; =1 land/seaice
@@ -469,6 +470,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
     NoahmpIO%ice             = undefined_int
     NoahmpIO%ivgtyp          = undefined_int
     NoahmpIO%isltyp          = undefined_int
+    NoahmpIO%isctyp          = undefined_int
     NoahmpIO%isnowxy         = undefined_int
     NoahmpIO%coszen          = undefined_real
     NoahmpIO%xlat            = undefined_real

src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90

@@ -68,6 +68,7 @@ module NoahmpIOVarType
     integer                                                ::  soil_update_steps   ! number of model time steps to update soil process
     integer,                allocatable, dimension(:)      ::  ivgtyp              ! vegetation type
     integer,                allocatable, dimension(:)      ::  isltyp              ! soil type
+    integer,                allocatable, dimension(:)      ::  isctyp              ! soil colour class
     real(kind=kind_noahmp), allocatable, dimension(:)      ::  coszen              ! cosine zenith angle
     real(kind=kind_noahmp), allocatable, dimension(:)      ::  xlat                ! latitude [rad]
     real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  dz8w                ! thickness of atmo layers [m]
@@ -516,6 +517,7 @@ module NoahmpIOVarType
 
     CHARACTER(LEN=256)                                     ::  mminsl  = 'STAS'    ! soil classification
     CHARACTER(LEN=256)                                     ::  llanduse            ! (=USGS, using USGS landuse classification)
+    CHARACTER(LEN=256)                                     ::  lsoilcol            ! (=Noah, using Noah's native soil colour classification)
 
 !------------------------------------------------------------------------
 ! Timing: