* In ./src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F, add…

…ed the initialization

  of variable sh2o (soil liquid water) prior to calling subroutine NoahmpInitMain in
  subroutine noahmp_init.

src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F

@@ -10,6 +10,7 @@ module mpas_atmphys_lsm_noahmpinit
  use mpas_log
  use mpas_pool_routines
 
+ use mpas_atmphys_constants,only: grav => gravity, t0 => svpt0
  use mpas_atmphys_utilities,only: physics_error_fatal
  use mpas_atmphys_vars,only     : mpas_noahmp
 
@@ -253,8 +254,12 @@ subroutine noahmp_init(configs,mesh,diag_physics,diag_physics_noahmp,output_noah
 !local variables and pointers:
  logical,pointer:: do_restart
  logical,parameter:: fndsnowh = .true.
+
  integer:: i,its,ite,ns,nsoil,nsnow,nzsnow
 
+ real(kind=RKIND),parameter:: hlice = 3.335E5
+ real(kind=RKIND):: bexp,fk,smcmax,psisat
+
 !-----------------------------------------------------------------------------------------------------------------
 !call mpas_log_write(' ')
 !call mpas_log_write('--- enter subroutine noahmp_init:')
@@ -390,6 +395,40 @@ subroutine noahmp_init(configs,mesh,diag_physics,diag_physics_noahmp,output_noah
  call mpas_pool_get_array(output_noahmp,'t2mvxy',t2mvxy  )
  call mpas_pool_get_array(output_noahmp,'qtdrain',qtdrain)
 
+!--- initialization of the soil liquid water content:
+ do i = its,ite
+    if(ivgtyp(i) == mpas_noahmp%isice_table .and. xice(i) .le. 0._RKIND) then
+       do ns = 1,nsoil
+          smois(ns,i) = 1._RKIND
+          sh2o(ns,i)  = 0._RKIND
+          tslb(ns,i)  = min(tslb(ns,i),263.15)
+       enddo
+    else
+       bexp   = mpas_noahmp%bexp_table(isltyp(i))
+       smcmax = mpas_noahmp%smcmax_table(isltyp(i))
+       psisat = mpas_noahmp%psisat_table(isltyp(i))
+
+       do ns = 1,nsoil
+          if(smois(ns,i) > smcmax) smois(ns,i) = smcmax
+       enddo
+       if(bexp.gt.0. .and. smcmax.gt.0. .and. psisat.gt.0.) then
+          do ns = 1,nsoil
+             if(tslb(ns,i) .lt. 273.149) then
+                fk = ( ((hlice/(grav*(-psisat)))*((tslb(ns,i)-t0)/tslb(ns,i)))**(-1/bexp) )*smcmax
+                fk = max(fk,0.02)
+                sh2o(ns,i) = min(fk,smois(n,i))
+             else
+                sh2o(ns,i) = smois(ns,i)
+             endif
+          enddo
+       else
+          do ns = 1,nsoil
+             sh2o(ns,i) = smois(ns,i)
+          enddo
+       endif
+    endif
+ enddo
+
 
  do i = its,ite
     mpas_noahmp%tmn(i)   = tmn(i)
@@ -419,9 +458,9 @@ subroutine noahmp_init(configs,mesh,diag_physics,diag_physics_noahmp,output_noah
     if(snow(i) .gt. 0._RKIND) snowc(i) = 1.
 
     do ns = 1,nsoil
-       mpas_noahmp%sh2o(i,ns)  = sh2o(ns,i)
-       mpas_noahmp%smois(i,ns) = smois(ns,i)
-       mpas_noahmp%tslb(i,ns)  = tslb(ns,i)
+       sh2o(ns,i)  = mpas_noahmp%sh2o(i,ns)
+       smois(ns,i) = mpas_noahmp%smois(i,ns)
+       tslb(ns,i)  = mpas_noahmp%tslb(i,ns)
    enddo
  enddo