wrfv2_fire/phys/module_sf_noahmpdrv.F

   1 MODULE module_sf_noahmpdrv
   2
   3 !-------------------------------
   4   USE module_sf_noahmplsm
   5   USE module_sf_urban
   6   USE module_model_constants, ONLY : R_D, CP, XLF, XLV, RHOWATER, KARMAN
   7   USE module_sf_noahdrv, ONLY : SOIL_VEG_GEN_PARM
   8   USE module_sf_noah_seaice
   9   USE module_sf_noahlsm_glacial_only
  10   USE MODULE_RA_GFDLETA, ONLY: CAL_MON_DAY
  11 #ifdef WRF_CHEM
  12   USE module_data_gocart_dust
  13 #endif
  14 !-------------------------------
  15
  16 !
  17 CONTAINS
  18 !
  19   SUBROUTINE noahmplsm(DZ8W,QV3D,P8W3D,T3D,TSK,                      &
  20        HFX,QFX,LH,GRDFLX, QGH,GSW,SWDOWN,GLW,SMSTAV,SMSTOT, &
  21        SFCRUNOFF, UDRUNOFF,IVGTYP,ISLTYP,VEGFRA,     &
  22        ALBEDO,ALBBCK,ZNT,Z0,TMN,XLAND,XICE,XICE_THRESHOLD,ISICE,EMISS,EMBCK,   &
  23        SNOWC,QSFC,RAINBL,                            &
  24        num_soil_layers,DT,DZS,ITIMESTEP,             &
  25        SMOIS,TSLB,SNOW,CANWAT,                       &
  26        CHS,CHS2,CQS2,CPM,ROVCP,SR,chklowq,qz0,       & !H
  27        myj,RIB,frpcpn,                               &
  28        SH2O,SNOWH,                                   & !H
  29        U_PHY,V_PHY,                                  & !I
  30        COSZ_URB2D, XLAT_URB2D,                       & !I
  31        SNOALB,                                       & !I
  32        SNOTIME,ACSNOM,ACSNOW,                        & !O
  33        idveg    ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ,iopt_frz ,iopt_inf , &
  34        iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot,iopt_stc ,                     &
  35        isnowxy  ,tvxy     ,tgxy     ,canicexy ,                               &
  36        canliqxy ,eahxy    ,tahxy    ,cmxy     ,chxy     ,                     &
  37        fwetxy   ,sneqvoxy ,alboldxy ,qsnowxy  ,wslakexy ,zwtxy    ,waxy     , &
  38        wtxy     ,tsnoxy   ,zsnsoxy  ,snicexy  ,snliqxy  ,lfmassxy ,rtmassxy , &
  39        stmassxy ,woodxy   ,stblcpxy ,fastcpxy ,xlaixy   ,xsaixy   ,           &
  40        tradxy   ,tsxy     ,neexy    ,gppxy    ,nppxy    ,fvegxy   ,qinxy    , &
  41        runsfxy  ,runsbxy  ,ecanxy   ,edirxy   ,etranxy  ,fsaxy    ,firaxy   , &
  42        aparxy   ,psnxy    ,savxy    ,sagxy    ,                               &
  43        fsnoxy   ,YR       ,JULIAN   ,                                         &
  44        potevp,                                       & !O
  45
  46 !jref:start
  47        qcxy     ,pblhxy   ,isurban  ,iz0tlnd  ,dx     , & !I
  48        chstarxy ,t2mvxy   ,t2mbxy   ,rssunxy  ,rsshaxy, bgapxy ,wgapxy ,gapxy,   & !O
  49        tgvxy    ,tgbxy    ,q2mvxy   ,q2mbxy  ,shdmaxxy, chvxy  ,chbxy  ,         &
  50 !jref:end
  51
  52        ids,ide, jds,jde, kds,kde,                    &
  53        ims,ime, jms,jme, kms,kme,                    &
  54        its,ite, jts,jte, kts,kte                     )
  55 !----------------------------------------------------------------
  56     IMPLICIT NONE
  57 !----------------------------------------------------------------
  58 !----------------------------------------------------------------
  59 ! --- atmospheric (WRF generic) variables
  60 !-- DT          time step (seconds)
  61 !-- DZ8W        thickness of layers (m)
  62 !-- T3D         temperature (K)
  63 !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
  64 !-- P3D         3D pressure (Pa)
  65 !-- FLHC        exchange coefficient for heat (m/s)
  66 !-- FLQC        exchange coefficient for moisture (m/s)
  67 !-- PSFC        surface pressure (Pa)
  68 !-- XLAND       land mask (1 for land, 2 for water)
  69 !-- QGH         saturated mixing ratio at 2 meter
  70 !-- GSW         downward short wave flux at ground surface (W/m^2)
  71 !-- GLW         downward long wave flux at ground surface (W/m^2)
  72 !-- History variables
  73 !-- CANWAT      canopy moisture content (mm)
  74 !-- TSK         surface temperature (K)
  75 !-- TSLB        soil temp (k)
  76 !-- SMOIS       total soil moisture content (volumetric fraction)
  77 !-- SH2O        unfrozen soil moisture content (volumetric fraction)
  78 !                note: frozen soil moisture (i.e., soil ice) = SMOIS - SH2O
  79 !-- SNOWH       actual snow depth (m)
  80 !-- SNOW        liquid water-equivalent snow depth (m)
  81 !-- ALBEDO      time-varying surface albedo including snow effect (unitless fraction)
  82 !-- ALBBCK      background surface albedo (unitless fraction)
  83 !-- CHS          surface exchange coefficient for heat and moisture (m s-1);
  84 !-- CHS2        2m surface exchange coefficient for heat  (m s-1);
  85 !-- CQS2        2m surface exchange coefficient for moisture (m s-1);
  86 ! --- soil variables
  87 !-- num_soil_layers   the number of soil layers
  88 !-- ZS          depths of centers of soil layers   (m)
  89 !-- DZS         thicknesses of soil layers (m)
  90 !-- SLDPTH      thickness of each soil layer (m, same as DZS)
  91 !-- TMN         soil temperature at lower boundary (K)
  92 !-- SMCMAX      porosity, i.e. saturated value of soil moisture (volumetric)
  93 !-- NROOT       number of root layers, a function of veg type, determined
  94 !               in subroutine redprm.
  95 !-- SMSTAV      Soil moisture availability for evapotranspiration (
  96 !                   fraction between SMCWLT and SMCMXA)
  97 !-- SMSTOT      Total soil moisture content frozen+unfrozen) in the soil column (mm)
  98 ! --- snow variables
  99 !-- SNOWC       fraction snow coverage (0-1.0)
 100 ! --- vegetation variables
 101 !-- SNOALB      upper bound on maximum albedo over deep snow
 102 !-- Z0BRD       Background fixed roughness length (M)
 103 !-- Z0          Background vroughness length (M) as function
 104 !-- ZNT         Time varying roughness length (M) as function
 105 !-- ALBD(IVGTPK,ISN) background albedo reading from a table
 106 ! --- LSM output
 107 !-- HFX         upward heat flux at the surface (W/m^2)
 108 !-- QFX         upward moisture flux at the surface (kg/m^2/s)
 109 !-- LH          upward moisture flux at the surface (W m-2)
 110 !-- GRDFLX(I,J) ground heat flux (W m-2)
 111 !-- FDOWN       radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
 112 !----------------------------------------------------------------------------
 113 !-- EC          canopy water evaporation ((W m-2)
 114 !-- EDIR        direct soil evaporation (W m-2)
 115 !-- ESNOW       sublimation from (or deposition to if <0) snowpack (W m-2)
 116 !-- DEW         dewfall (or frostfall for t<273.15) (M)
 117 ! ----------------------------------------------------------------------
 118 !-- ETP         potential evaporation (W m-2)
 119 ! ----------------------------------------------------------------------
 120 !-- FLX1        precip-snow sfc (W m-2)
 121 !-- FLX2        freezing rain latent heat flux (W m-2)
 122 !-- FLX3        phase-change heat flux from snowmelt (W m-2)
 123 ! ----------------------------------------------------------------------
 124 !-- ACSNOM      snow melt (mm) (water equivalent)
 125 !-- ACSNOW      accumulated snow fall (mm) (water equivalent)
 126 !-- POTEVP      accumulated potential evaporation (W/m^2)
 127 !-- RIB         Bulk Richardson number from SFCLAY routine
 128 ! ----------------------------------------------------------------------
 129 !-- RUNOFF1     surface runoff (m s-1), not infiltrating the surface
 130 !-- RUNOFF2     subsurface runoff (m s-1), drainage out bottom of last
 131 !                  soil layer (baseflow)
 132 ! ----------------------------------------------------------------------
 133 !-- RC          canopy resistance (s m-1)
 134 !-- PC          plant coefficient (unitless fraction, 0-1) where PC*ETP = actual transp
 135 !-- EMISS       surface emissivity (between 0 and 1)
 136 !-- EMBCK       Background surface emissivity (between 0 and 1)
 137 !-- SHDMAX      Maximum vegetation fraction
 138
 139 !-- ROVCP       R/CP
 140 !               (R_d/R_v) (dimensionless)
 141 !-- ids         start index for i in domain
 142 !-- ide         end index for i in domain
 143 !-- jds         start index for j in domain
 144 !-- jde         end index for j in domain
 145 !-- kds         start index for k in domain
 146 !-- kde         end index for k in domain
 147 !-- ims         start index for i in memory
 148 !-- ime         end index for i in memory
 149 !-- jms         start index for j in memory
 150 !-- jme         end index for j in memory
 151 !-- kms         start index for k in memory
 152 !-- kme         end index for k in memory
 153 !-- its         start index for i in tile
 154 !-- ite         end index for i in tile
 155 !-- jts         start index for j in tile
 156 !-- jte         end index for j in tile
 157 !-- kts         start index for k in tile
 158 !-- kte         end index for k in tile
 159 !
 160 !-- SR          fraction of frozen precip (0.0 to 1.0)
 161 !----------------------------------------------------------------
 162
 163 ! IN only
 164
 165     INTEGER,  INTENT(IN   )   ::     ids,ide, jds,jde, kds,kde,  &
 166          &                           ims,ime, jms,jme, kms,kme,  &
 167          &                           its,ite, jts,jte, kts,kte
 168
 169     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
 170     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
 171     REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
 172     REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
 173
 174
 175     REAL,    DIMENSION( ims:ime, jms:jme )                     , &
 176          &   INTENT(IN   )    ::                            TMN, &
 177          &                                                XLAND, &
 178          &                                                 XICE, &
 179          &                                               VEGFRA, &
 180          &                                               SNOALB, &
 181          &                                                  GSW, &
 182          &                                               SWDOWN, &
 183          &                                                  GLW, &
 184          &                                                   Z0, &
 185          &                                               ALBBCK, &
 186          &                                               RAINBL, &
 187          &                                                EMBCK, &
 188          &                                                   SR
 189
 190
 191     REAL,    DIMENSION( ims:ime, kms:kme, jms:jme )            , &
 192          INTENT(IN   )    ::                           QV3D, &
 193          p8w3D, &
 194          DZ8W, &
 195          T3D
 196 !jref:start - changed to inout
 197     REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 198          INTENT(INOUT   )               ::               QGH,  &
 199          CHS,   &
 200          CPM
 201 !jref:end
 202
 203     INTEGER, DIMENSION( ims:ime, jms:jme )                     , &
 204          INTENT(IN   )    ::                         IVGTYP, &
 205          ISLTYP
 206
 207     INTEGER, INTENT(IN)       ::     num_soil_layers,ITIMESTEP
 208 !jref:start - xice_threshold
 209     REAL,     INTENT(IN   )   ::     DT,ROVCP,XICE_THRESHOLD
 210     INTEGER,  INTENT(IN   )   ::     ISICE
 211 !jref:end
 212
 213     REAL,     DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
 214
 215 ! IN and OUT
 216
 217     REAL,     DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
 218          &    INTENT(INOUT)   ::                          SMOIS, &
 219          &                                                 SH2O, &
 220          &                                                 TSLB
 221
 222     REAL,    DIMENSION( ims:ime, jms:jme )                     , &
 223          &   INTENT(INOUT)    ::                            TSK, &
 224          &                                                  HFX, &
 225          &                                                  QFX, &
 226          &                                                   LH, &
 227          &                                               GRDFLX, &
 228          &                                                 QSFC, &
 229          &                                                 CQS2, &
 230          &                                                 CHS2, &
 231          &                                                 SNOW, &
 232          &                                                SNOWC, &
 233          &                                                SNOWH, &
 234          &                                               CANWAT, &
 235          &                                               SMSTAV, &
 236          &                                               SMSTOT, &
 237          &                                            SFCRUNOFF, &
 238          &                                             UDRUNOFF, &
 239          &                                               ACSNOM, &
 240          &                                               ACSNOW, &
 241          &                                                EMISS, &
 242          &                                               POTEVP, &
 243          &                                                  RIB, &
 244          &                                               ALBEDO, &
 245          &                                                  ZNT
 246
 247     REAL,    DIMENSION( ims:ime, jms:jme )                     , &
 248          INTENT(OUT)    ::                        CHKLOWQ
 249     REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) ::        QZ0
 250
 251 !niuin:
 252 ! in
 253     INTEGER,  INTENT(IN) :: idveg     !dynamic vegetation (1 -> off ; 2 -> on) with opt_crs = 1
 254     INTEGER,  INTENT(IN) :: iopt_crs  !canopy stomatal resistance (1-> Ball-Berry; 2->Jarvis)
 255     INTEGER,  INTENT(IN) :: iopt_btr  !soil moisture factor for stomatal resistance (1-> Noah; 2-> CLM; 3-> SSiB)
 256     INTEGER,  INTENT(IN) :: iopt_run  !runoff and groundwater (1->SIMGM; 2->SIMTOP; 3->Schaake96; 4->BATS)
 257     INTEGER,  INTENT(IN) :: iopt_sfc  !surface layer drag coeff (CH & CM) (1->M-O; 2->Chen97)
 258     INTEGER,  INTENT(IN) :: iopt_frz  !supercooled liquid water (1-> NY06; 2->Koren99)
 259     INTEGER,  INTENT(IN) :: iopt_inf  !frozen soil permeability (1-> NY06; 2->Koren99)
 260     INTEGER,  INTENT(IN) :: iopt_rad  !radiation transfer (1->gap=F(3D,cosz); 2->gap=0; 3->gap=1-Fveg)
 261     INTEGER,  INTENT(IN) :: iopt_alb  !snow surface albedo (1->BATS; 2->CLASS)
 262     INTEGER,  INTENT(IN) :: iopt_snf  !rainfall & snowfall (1-Jordan91; 2->BATS; 3->Noah)
 263     INTEGER,  INTENT(IN) :: iopt_tbot !lower boundary of soil temperature (1->zero-flux; 2->Noah)
 264     INTEGER,  INTENT(IN) :: iopt_stc  !snow/soil temperature time scheme
 265
 266 ! in & out
 267     INTEGER, INTENT(IN) :: YR
 268     REAL,    INTENT(IN) :: JULIAN
 269     INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: isnowxy     !actual no. of snow layers
 270     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tvxy        !vegetation canopy temperature
 271     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tgxy        !ground surface temperature
 272     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canicexy    !canopy-intercepted ice (mm)
 273     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canliqxy    !canopy-intercepted liquid water (mm)
 274     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: eahxy       !canopy air vapor pressure (pa)
 275     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tahxy       !canopy air temperature (k)
 276     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: cmxy        !momentum drag coefficient
 277     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chxy        !sensible heat exchange coefficient
 278     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fwetxy      !wetted or snowed fraction of the canopy (-)
 279     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: sneqvoxy    !snow mass at last time step(mm h2o)
 280     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: alboldxy    !snow albedo at last time step (-)
 281     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: qsnowxy     !snowfall on the ground [mm/s]
 282     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wslakexy    !lake water storage [mm]
 283
 284     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: zwtxy       !water table depth [m]
 285     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: waxy        !water in the "aquifer" [mm]
 286     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wtxy        !groundwater storage [mm]
 287     REAL, DIMENSION(ims:ime,-2:num_soil_layers,jms:jme), INTENT(INOUT) :: zsnsoxy  !snow layer depth [m]
 288     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: tsnoxy   !snow temperature [K]
 289     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: snicexy  !snow layer ice [mm]
 290     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: snliqxy  !snow layer liquid water [mm]
 291     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: lfmassxy    !leaf mass [g/m2]
 292     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: rtmassxy    !mass of fine roots [g/m2]
 293     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stmassxy    !stem mass [g/m2]
 294     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: woodxy      !mass of wood (incl. woody roots) [g/m2]
 295     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stblcpxy    !stable carbon in deep soil [g/m2]
 296     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fastcpxy    !short-lived carbon, shallow soil [g/m2]
 297     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xlaixy      !leaf area index
 298     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xsaixy      !stem area index
 299 !jref:start
 300     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: SNOTIME      !snow age time
 301 !jref:end
 302
 303 !out
 304     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: tradxy      !surface radiative temperature (k)
 305     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: tsxy        !surface temperature (k)
 306     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: neexy       !net ecosys exchange (g/m2/s CO2)
 307     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: gppxy       !gross primary assimilation [g/m2/s C]
 308     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: nppxy       !net primary productivity [g/m2/s C]
 309     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: fvegxy      !greenness vegetation fraction [-]
 310     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: qinxy       !groundwater recharge [mm/s]
 311     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: runsfxy     !surface runoff [mm/s]
 312     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: runsbxy     !subsurface runoff [mm/s]
 313     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: ecanxy      !evaporation of intercepted water (mm/s)
 314     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: edirxy      !soil surface evaporation rate (mm/s]
 315     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: etranxy     !transpiration rate (mm/s)
 316     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: fsaxy       !total absorbed solar radiation (w/m2)
 317     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: firaxy      !total net longwave rad (w/m2) [+ to atm]
 318     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: aparxy      !photosyn active energy by canopy (w/m2)
 319     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: psnxy       !total photosynthesis (umol co2/m2/s) [+]
 320     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: savxy       !solar rad absorbed by veg. (w/m2)
 321     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: sagxy       !solar rad absorbed by ground (w/m2)
 322     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT)  :: fsnoxy      !snow cover fraction (-)
 323
 324 !jref:start
 325     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)  :: chstarxy    !effective ch
 326     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)  :: t2mvxy      !2m temperature of vegetation part
 327     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)  :: t2mbxy      !2m temperature of bare ground part
 328     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)  :: q2mvxy      !2m mixing ratio of vegetation part
 329     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)  :: q2mbxy      !2m mixing ratio of bare ground part
 330
 331     REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN)   :: qcxy        !cloud water mixing ratio
 332     REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN)   :: pblhxy      !Planetary boundary layer from sfclay
 333     INTEGER                          , INTENT(IN)   :: isurban
 334     INTEGER                          , INTENT(IN)   :: iz0tlnd
 335     REAL                             , INTENT(IN)   :: dx
 336     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: rssunxy        !sunlit leaf stomatal resistance (s/m)
 337     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: rsshaxy        !shaded leaf stomatal resistance (s/m)
 338     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: bgapxy        !between gap fraction
 339     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: wgapxy        !within gap fraction
 340     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: gapxy        !within gap fraction
 341     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: tgvxy
 342     REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: tgbxy
 343     REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN )::   shdmaxxy
 344     REAL, DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: chvxy        !sensible heat exchange coefficient vegetated
 345     REAL, DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: chbxy        !sensible heat exchange coefficient bare-ground
 346 !jref:end
 347
 348
 349 !niuout
 350
 351 ! Local variables (moved here from driver to make routine thread safe, 20031007 jm)
 352
 353     INTEGER :: YEARLEN
 354
 355     REAL  ::  ETP, SSOIL,EC, ESNOW,             &
 356          FLX1,FLX2,FLX3,DEW,FDOWN,RC,PC,FFROZP
 357
 358 !niuin
 359 !locals (prognostic):
 360
 361     INTEGER                             :: isnow        !actual no. of snow layers
 362     REAL, DIMENSION(-2:num_soil_layers) :: stc          !snow/soil tmperatures
 363     REAL, DIMENSION( 1:num_soil_layers) :: smc          !vol. soil moisture (m3/m3)
 364     REAL, DIMENSION( 1:num_soil_layers) :: smh2o        !vol. soil liquid water (m3/m3)
 365     REAL                                :: tv           !vegetation canopy temperature
 366     REAL                                :: tg           !ground surface temperature
 367     REAL                                :: canice       !canopy-intercepted ice (mm)
 368     REAL                                :: canliq       !canopy-intercepted liquid water (mm)
 369     REAL                                :: snowd        !snow depth (m)
 370     REAL                                :: swe          !snow water equivalent (mm)
 371     REAL                                :: eah          !canopy air vapor pressure (pa)
 372     REAL                                :: tah          !canopy air temperature (k)
 373     REAL                                :: cm           !momentum drag coefficient
 374     REAL                                :: ch           !sensible heat exchange coefficient
 375     REAL                                :: fwet         !wetted or snowed fraction of the canopy (-)
 376     REAL                                :: sneqvo       !snow mass at last time step(mm h2o)
 377     REAL                                :: albold       !snow albedo at last time step (-)
 378     REAL                                :: qsnow        !snowfall on the ground [mm/s]
 379     REAL                                :: wslake       !lake water storage [mm]
 380
 381     REAL                                :: zwt          !water table depth [m]
 382     REAL                                :: wa           !water in the "aquifer" [mm]
 383     REAL                                :: wt           !groundwater storage [mm]
 384     REAL, DIMENSION(-2:num_soil_layers) :: zsnso        !snow layer depth [m]
 385     REAL, DIMENSION(-2:              0) :: tsno         !snow temperature [K]
 386     REAL, DIMENSION(-2:              0) :: snice        !snow layer ice [mm]
 387     REAL, DIMENSION(-2:              0) :: snliq        !snow layer liquid water [mm]
 388     REAL                                :: lfmass       !leaf mass [g/m2]
 389     REAL                                :: rtmass       !mass of fine roots [g/m2]
 390     REAL                                :: stmass       !stem mass [g/m2]
 391     REAL                                :: wood         !mass of wood (incl. woody roots) [g/m2]
 392     REAL                                :: stblcp       !stable carbon in deep soil [g/m2]
 393     REAL                                :: fastcp       !short-lived carbon, shallow soil [g/m2]
 394     REAL                                :: plai         !leaf area index
 395     REAL                                :: psai         !stem area index
 396
 397 !jref:start
 398     REAL                                :: chstar2
 399     REAL                                :: cqstar2
 400     REAL                                :: chstar       !effective ch
 401     REAL                                :: tstar
 402     REAL                                :: t2mv         !2m temperature of vegetation part
 403     REAL                                :: t2mb         !2m temperature of bare ground part
 404     REAL                                :: q2mv         !2m mixing ratio of vegetation part
 405     REAL                                :: q2mb         !2m mixing ratio of bare ground part
 406     REAL                                :: qc           !
 407     REAL                                :: t2m
 408     REAL                                :: pblh
 409     REAL                                :: qsfc1d
 410     REAL, DIMENSION(ims:ime,jms:jme)    :: tstarxy      !effective skin temperature
 411     REAL, DIMENSION(ims:ime,jms:jme)    :: chstar2xy    !effective 2m exchange coefficients
 412     REAL                                :: rssun
 413     REAL                                :: rssha
 414     REAL                                :: bgap
 415     REAL                                :: wgap
 416     REAL                                :: gap
 417     REAL                                :: tgv
 418     REAL                                :: tgb
 419     REAL                                :: snowhk
 420     REAL                                :: snotime1
 421     REAL                                :: qv1d         !mixing ratio
 422     REAL                                :: dz8w1d
 423     REAL                                :: shdmax
 424     REAL                                :: chv          !sensible heat exchange coefficient vegetated
 425     REAL                                :: chb          !sensible heat exchange coefficient bare-ground
 426 !jref:end
 427
 428 !out (outputs)
 429
 430     REAL                                :: trad         !surface radiative temperature (k)
 431     REAL                                :: ts           !surface temperature (k)
 432     REAL                                :: nee          !net ecosys exchange (g/m2/s CO2)
 433     REAL                                :: gpp          !gross primary assimilation [g/m2/s C]
 434     REAL                                :: npp          !net primary productivity [g/m2/s C]
 435     REAL                                :: fveg         !greenness vegetation fraction [-]
 436     REAL                                :: qin          !groundwater recharge [mm/s]
 437     REAL                                :: runsf        !surface runoff [mm/s]
 438     REAL                                :: runsb        !subsurface runoff [mm/s]
 439     REAL                                :: ecan         !evaporation of intercepted water (mm/s)
 440     REAL                                :: esoil        !soil surface evaporation rate (mm/s]
 441     REAL                                :: etran        !transpiration rate (mm/s)
 442     REAL                                :: fsa          !total absorbed solar radiation (w/m2)
 443     REAL                                :: fira         !total net longwave rad (w/m2) [+ to atm]
 444     REAL                                :: fsh          !total sensible heat (w/m2) [+ to atm]
 445     REAL                                :: flh          !total latent heat (w/m2) [+ to atm]
 446     REAL                                :: apar         !photosyn active energy by canopy (w/m2)
 447     REAL                                :: psn          !total photosynthesis (umol co2/m2/s) [+]
 448     REAL                                :: sav          !solar rad absorbed by veg. (w/m2)
 449     REAL                                :: sag          !solar rad absorbed by ground (w/m2)
 450     REAL                                :: fsno         !snow cover fraction (-)
 451     REAL                                :: salb         !surface albedo (-)
 452     REAL                                :: errwat
 453     REAL                                :: qmelt
 454     REAL                                :: ponding
 455     REAL                                :: ponding1
 456     REAL                                :: ponding2
 457
 458 !local
 459
 460     real                                :: fsr          !total reflected solar radiation (w/m2)
 461     real                                :: fcev         !canopy evaporation heat (w/m2) [+ to atm]
 462     real                                :: fgev         !ground evaporation heat (w/m2) [+ to atm]
 463     real                                :: fctr         !transpiration heat flux (w/m2) [+ to atm]
 464     real, dimension(-2:              0) :: ficeold      !snow layer liquid water [mm]
 465
 466     INTEGER :: ILOC      !grid index
 467     INTEGER :: JLOC      !grid index
 468     INTEGER :: ISC       !soil color index
 469     INTEGER :: IST       !surface type 1-soil; 2-lake
 470
 471 !niuout
 472
 473     LOGICAL,    INTENT(IN   )    ::     myj,frpcpn
 474
 475 ! DECLARATIONS - LOGICAL
 476 ! ----------------------------------------------------------------------
 477     LOGICAL, PARAMETER :: LOCAL=.false.
 478     LOGICAL :: FRZGRA, SNOWNG
 479
 480     LOGICAL :: IPRINT
 481
 482 ! ----------------------------------------------------------------------
 483 ! DECLARATIONS - INTEGER
 484 ! ----------------------------------------------------------------------
 485     INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
 486     INTEGER :: NROOT
 487     INTEGER :: KZ ,K
 488     INTEGER :: NS
 489 ! ----------------------------------------------------------------------
 490 ! DECLARATIONS - REAL
 491 ! ----------------------------------------------------------------------
 492
 493     REAL  :: DQSDT2, LWDN, PRCP, PSFC, UU, VV, CO2AIR, O2AIR,         &
 494          &   Q2SAT,Q2SATI,SFCPRS,SFCTMP,SHDFAC,SNOALB1,               &
 495          &   SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ETA, ETA_KINEMATIC,         &
 496          &   EMBRD, FOLN, LAT,                                        &
 497          &   Z0K,RUNOFF1,RUNOFF2,SOLNET,E2SAT,SFCTSNO
 498
 499     REAL :: RIBB
 500     REAL ::  FDTW
 501
 502     REAL  :: EMISSI
 503
 504     REAL  :: SNCOVR,SNEQV,CHK,TH2
 505
 506     REAL  :: SMCMAX,SNOMLT,SOILM,SOILW,Q1,T1
 507
 508     REAL  :: Z0BRD
 509 !
 510     REAL  :: COSZ
 511 !
 512 !niu      REAL, DIMENSION(1:num_soil_layers)::  SLDPTH, STC,SMC,SWC
 513     REAL, DIMENSION(1:num_soil_layers)::  SLDPTH,SWC
 514
 515 !jref:start
 516     REAL, DIMENSION(1:num_soil_layers):: STCNEW
 517 !jref:end
 518 !
 519     REAL, DIMENSION(1:num_soil_layers) ::     ZSOIL, RTDIS
 520     REAL, PARAMETER  :: TRESH=.95E0, A2=17.67,A3=273.15,A4=29.65,   &
 521          T0=273.16E0, ELWV=2.50E6,  A23M4=A2*(A3-A4)
 522
 523     ! Used for calculating the 2-m Potential Temperature:
 524     REAL, PARAMETER                                       :: CAPA=R_D/CP
 525     REAL                                                  :: APELM
 526     REAL                                                  :: APES
 527     REAL                                                  :: SFCTH2
 528
 529 ! ----------------------------------------------------------------------
 530 ! ----------------------------------------------------------------------
 531
 532 ! MEK JUL2007
 533     FDTW=DT/(XLV*RHOWATER)
 534 ! debug printout
 535     IPRINT=.false.
 536
 537 !      SLOPETYP=2
 538     SLOPETYP=1
 539 !      SHDMIN=0.00
 540
 541     YEARLEN = 365
 542     if (mod(YR,4) == 0) then
 543        YEARLEN = 366
 544        if (mod(YR,100) == 0) then
 545           YEARLEN = 365
 546           if (mod(YR,400) == 0) then
 547              YEARLEN = 366
 548           endif
 549        endif
 550     endif
 551
 552     NSOIL=num_soil_layers
 553
 554     DO NS=1,NSOIL
 555        SLDPTH(NS)=DZS(NS)
 556     ENDDO
 557
 558     call noahmp_options(idveg     ,iopt_crs  ,iopt_btr  ,iopt_run  ,iopt_sfc  ,iopt_frz , &
 559          iopt_inf  ,iopt_rad  ,iopt_alb  ,iopt_snf  ,iopt_tbot, iopt_stc )
 560
 561     ISC = 4                  ! soil color: assuming a middle color category ?????????
 562
 563     ZSOIL(1) = -SLDPTH(1)    ! move out of  x-y do loops
 564     DO KZ = 2, NSOIL
 565        ZSOIL(KZ)         = -SLDPTH(KZ) + ZSOIL(KZ-1)
 566     END DO
 567     FOLN                 = 1.0
 568
 569 !niuout
 570
 571     DO J=jts,jte
 572
 573        IF(ITIMESTEP.EQ.1)THEN
 574           DO I=its,ite
 575
 576 !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
 577              IF((XLAND(I,J)-1.5).GE.0.)THEN
 578 ! check sea-ice point
 579                 IF(XICE(I,J).EQ.1..and.IPRINT)PRINT*,' sea-ice at water point, I=',I,  &
 580                      'J=',J
 581 !***   Open Water Case
 582                 SMSTAV(I,J)=1.0
 583                 SMSTOT(I,J)=1.0
 584                 DO NS=1,NSOIL
 585                    SMOIS(I,NS,J)=1.0
 586                    TSLB(I,NS,J)=273.16                                          !STEMP
 587                 ENDDO
 588              ELSE
 589                 IF(XICE(I,J).EQ.1.)THEN
 590 !***        SEA-ICE CASE
 591                    SMSTAV(I,J)=1.0
 592                    SMSTOT(I,J)=1.0
 593                    DO NS=1,NSOIL
 594                       SMOIS(I,NS,J)=1.0
 595                    ENDDO
 596                 ENDIF
 597              ENDIF
 598 !
 599           ENDDO
 600        ENDIF                                                               ! end of initialization over ocean
 601
 602
 603 !-----------------------------------------------------------------------
 604        DO I=its,ite
 605 ! surface pressure
 606           PSFC=P8w3D(i,1,j)
 607 ! pressure in middle of lowest layer
 608           SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
 609 ! convert from mixing ratio to specific humidity
 610           Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
 611 !
 612 !         Q2SAT=QGH(I,j)
 613           Q2SAT=QGH(I,J)/(1.0+QGH(I,J))        ! Q2SAT is sp humidity
 614 ! add check on myj=.true.
 615 !        IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
 616           IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
 617              CHKLOWQ(I,J)=0.
 618           ELSE
 619              CHKLOWQ(I,J)=1.
 620           ENDIF
 621
 622           SFCTMP=T3D(i,1,j)
 623           ZLVL=0.5*DZ8W(i,1,j)
 624
 625 !        TH2=SFCTMP+(0.0097545*ZLVL)
 626 ! calculate SFCTH2 via Exner function vs lapse-rate (above)
 627           APES=(1.E5/PSFC)**CAPA
 628           APELM=(1.E5/SFCPRS)**CAPA
 629           SFCTH2=SFCTMP*APELM
 630           TH2=SFCTH2/APES
 631 !
 632           EMISSI = EMISS(I,J)
 633 !          LWDN=GLW(I,J)*EMISSI
 634           LWDN=GLW(I,J)
 635 ! SOLDN is total incoming solar
 636           SOLDN=SWDOWN(I,J)
 637 ! GSW is net downward solar
 638 !        SOLNET=GSW(I,J)
 639 ! use mid-day albedo to determine net downward solar (no solar zenith angle correction)
 640           SOLNET=SOLDN*(1.-ALBEDO(I,J))
 641           PRCP=RAINBL(i,j)/DT
 642           VEGTYP=IVGTYP(I,J)
 643           SOILTYP=ISLTYP(I,J)
 644           SHDFAC=VEGFRA(I,J)/100.
 645           T1=TSK(I,J)
 646           CHK=CHS(I,J)
 647           SNOALB1=SNOALB(I,J)     !NEW
 648
 649 ! if "SR" present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
 650 ! SR from e.g. Ferrier microphysics
 651 ! otherwise define from 1st atmos level temperature
 652           IF(FRPCPN) THEN
 653              FFROZP=SR(I,J)
 654           ELSE
 655              IF (SFCTMP <=  273.15) THEN
 656                 FFROZP = 1.0
 657              ELSE
 658                 FFROZP = 0.0
 659              ENDIF
 660           ENDIF
 661 !***
 662           IF((XLAND(I,J)-1.5).GE.0.)THEN                                  ! begining of land/sea if block
 663 ! Open water points
 664           ELSE
 665 ! Land or sea-ice case
 666
 667              IF (XICE(I,J) .GT. 0.5) THEN
 668                 ICE=1
 669              ELSE
 670                 ICE=0
 671              ENDIF
 672              DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
 673
 674              IF(SNOW(I,J).GT.0.0)THEN
 675 ! snow on surface (use ice saturation properties)
 676                 SFCTSNO=SFCTMP
 677                 E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
 678                 Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
 679                 Q2SATI=Q2SATI/(1.0+Q2SATI)    ! spec. hum.
 680                 IF(T1 .GT. 273.15)THEN
 681 ! warm ground temps, weight the saturation between ice and water according to SNOWC
 682                    Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
 683                    DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
 684                 ELSE
 685 ! cold ground temps, use ice saturation only
 686                    Q2SAT=Q2SATI
 687                    DQSDT2=Q2SATI*6174./(SFCTSNO**2)
 688                 ENDIF
 689 ! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
 690                 IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
 691              ENDIF
 692
 693              IF(ICE.EQ.0)THEN
 694                 TBOT=TMN(I,J)
 695              ELSE
 696                 TBOT=271.16
 697              ENDIF
 698              IF(VEGTYP.EQ.25) SHDFAC=0.0000
 699              IF(VEGTYP.EQ.26) SHDFAC=0.0000
 700              IF(VEGTYP.EQ.27) SHDFAC=0.0000
 701              IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
 702                 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
 703                 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
 704                 SOILTYP=7
 705              ENDIF
 706
 707 !-------------------------------------------
 708              ALBBRD=ALBBCK(I,J)
 709              Z0BRD=Z0(I,J)
 710              EMBRD=EMBCK(I,J)
 711 !jref:start - check if this is correct!! Maybe snowd
 712              RIBB=RIB(I,J)
 713              SNOTIME1 = SNOTIME(I,J)
 714 !jref:end
 715
 716 !FEI: temporaray arrays above need to be changed later by using SI
 717
 718 !niu          DO 70 NS=1,NSOIL
 719 !niu            SMC(NS)=SMOIS(I,NS,J)
 720 !niu            STC(NS)=TSLB(I,NS,J)                                          !STEMP
 721 !niu            SWC(NS)=SH2O(I,NS,J)
 722 !niu   70     CONTINUE
 723 !
 724
 725              IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
 726                   IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
 727                 VEGTYP = ISURBAN
 728              ENDIF
 729
 730              IST    = 1
 731              IF(VEGTYP == 16) IST = 2 ! lake points
 732
 733              CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,SLDPTH,ZSOIL,NSOIL,ISURBAN)
 734
 735              UU                   = U_PHY(I,1,J)
 736              VV                   = V_PHY(I,1,J)
 737              CO2AIR               = 395.E-06 * SFCPRS       !partial pressure co2 (pa)
 738              O2AIR                = 0.209    * SFCPRS       !partial pressure  o2 (pa)
 739              COSZ                 = COSZ_URB2D(I,J)
 740              LAT                  = XLAT_URB2D(I,J)
 741
 742              isnow                = isnowxy (i,j)
 743              stc  (isnow+1:    0) = tsnoxy  (i,isnow+1:    0,j)
 744              stc  (      1:nsoil) = tslb    (i,      1:nsoil,j)
 745              smc (       1:nsoil) = smois   (i,      1:nsoil,j)
 746              smh2o(      1:nsoil) = sh2o    (i,      1:nsoil,j)
 747              tv                   = tvxy    (i,j)
 748              tg                   = tgxy    (i,j)
 749              canliq               = canliqxy(i,j)
 750              canice               = canicexy(i,j)
 751              snowd                = snowh   (i,j)
 752              swe                  = snow    (i,j)
 753              eah                  = eahxy   (i,j)
 754              tah                  = tahxy   (i,j)
 755              cm                   = cmxy    (i,j)
 756              ch                   = chxy    (i,j)
 757 !jref:start
 758              chstar               = chs     (i,j)
 759              chstar2              = chs2    (i,j)
 760              cqstar2              = cqs2    (i,j)
 761              tstar                = T1
 762              qc                   = qcxy    (i,j)
 763              pblh                 = pblhxy  (i,j)
 764              qsfc1d               = qsfc    (i,j)
 765              t2mv                 = t2mvxy  (i,j)
 766              t2mb                 = t2mbxy  (i,j)
 767              q2mv                 = q2mvxy  (i,j)
 768              q2mb                 = q2mbxy  (i,j)
 769              qv1d                 = qv3d    (i,1,j) ! seaice/glacial needs mixing ratio (q2k = specific hum).
 770              dz8w1d               = dz8w    (i,1,j)
 771              shdmax               = shdmaxxy (i,j)/100. !fraction
 772 !jref:end
 773
 774              fwet                 = fwetxy  (i,j)
 775              sneqvo               = sneqvoxy(i,j)
 776              albold               = alboldxy(i,j)
 777              qsnow                = qsnowxy (i,j)
 778              wslake               = wslakexy(i,j)
 779
 780              zwt                  = zwtxy   (i,j)
 781              wa                   = waxy    (i,j)
 782              wt                   = wtxy    (i,j)
 783              zsnso(isnow+1:nsoil) = zsnsoxy (i,isnow+1:nsoil,j)
 784              snice(isnow+1:    0) = snicexy (i,isnow+1:    0,j)
 785              snliq(isnow+1:    0) = snliqxy (i,isnow+1:    0,j)
 786
 787              lfmass               = lfmassxy(i,j)
 788              rtmass               = rtmassxy(i,j)
 789              stmass               = stmassxy(i,j)
 790              wood                 = woodxy  (i,j)
 791              stblcp               = stblcpxy(i,j)
 792              fastcp               = fastcpxy(i,j)
 793              plai                 = xlaixy  (i,j)
 794              psai                 = xsaixy  (i,j)
 795
 796              ficeold(isnow+1:0)   = snicexy(i,isnow+1:0,j) &
 797                   /(snicexy(i,isnow+1:0,j)+snliqxy(i,isnow+1:0,j))
 798
 799 ! glacial, seaice split - jref
 800
 801     IF ( XICE(I,J) >= XICE_THRESHOLD ) THEN
 802
 803        SH2O  (i,1:nsoil,j) = 1.0
 804        XLAIXY(i,j)         = 0.01
 805
 806        cycle ! Skip any processing at sea-ice points
 807
 808     ELSE IF ( VEGTYP == ISICE ) THEN
 809
 810        SNCOVR = SNOWC(I,J)
 811        swe = swe*0.001 !jref mm -> m
 812        if ( (swe.ne.0..AND.snowd.eq.0.).or.(snowd.le.swe) )THEN
 813              snowd= 5.*swe
 814        endif
 815
 816        CALL SFLX_GLACIAL(I,J,ISICE,FFROZP,DT,ZLVL,NSOIL,SLDPTH,   &    !C
 817             &    LWDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,              &    !F
 818             &    TH2,Q2SAT,DQSDT2,                                &    !I
 819             &    ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, &    !S
 820             &    tstar,STC(1:NSOIL),snowd,swe,salb,chstar,        &    !H
 821             &    ETA,fsh, ETA_KINEMATIC,FDOWN,                    &    !O
 822             &    ESNOW,DEW,                                       &    !O
 823             &    ETP,SSOIL,                                       &    !O
 824             &    FLX1,FLX2,FLX3,                                  &    !O
 825             &    SNOMLT,SNCOVR,                                   &    !O
 826             &    runsf,                                         &    !O
 827             &    Q1,                                              &    !D
 828             &    SNOTIME1,                                        &
 829             &    RIBB)
 830
 831        tgb    = sfctmp ! Bare ground temperature will be the surface temperature over glacial points.
 832        tgv    = 0.0    ! Temperature under vegetation undefined over glacial points.
 833        swe    = swe*1000.
 834        plai   = 0.01  ! Should make this zero?
 835        smc    = 1.00
 836        smh2o  = 1.00  ! Something else?
 837        runsb  = 0.00
 838        fgev   = ETA
 839        fcev   = 0.
 840        fctr   = 0.
 841        soilm  = 1.0   ! Something else?
 842 !       SMAV    = 1.00  ! Something else?
 843        SNOWC(I,J) = 1.0
 844
 845        QFX(I,J) = eta_kinematic
 846        POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
 847
 848        CHS2(I,J) = CQS2(I,J)
 849        IF ( Q1 .GT. QSFC(I,J) ) THEN
 850           CQS2(I,J) = CHS(I,J)
 851        ENDIF
 852
 853     ELSE
 854 !jref:end
 855
 856              nee = -1.E36
 857              npp = -1.E36
 858
 859 #if 0
 860              if ( I == 15 .and. J == 5 ) then
 861                 ! Intent (IN) or Intent (INOUT), but not Intent (OUT)
 862                 write(*,'("Before call to NOAHMP_SFLX, at point ", I8, I8)') i, j
 863                 write(*,'(10x, "ICE        = ",  I10   )') ICE
 864                 write(*,'(10x, "IST        = ",  I10   )') IST
 865                 write(*,'(10x, "VEGTYP     = ",  I10   )') VEGTYP
 866                 write(*,'(10x, "ISC        = ",  I10   )') ISC
 867                 write(*,'(10x, "NSOIL      = ",  I10   )') NSOIL
 868                 write(*,'(10x, "ZSOIL      = ", 7F20.10)') ZSOIL
 869                 write(*,'(10x, "DT         = ",  F20.10)') DT
 870                 write(*,'(10x, "QV1D       = ",  F20.10)') QV1D
 871                 write(*,'(10x, "SFCTMP     = ",  F20.10)') SFCTMP
 872                 write(*,'(10x, "UU         = ",  F20.10)') UU
 873                 write(*,'(10x, "VV         = ",  F20.10)') VV
 874                 write(*,'(10x, "SOLDN      = ",  F20.10)') SOLDN
 875                 write(*,'(10x, "LWDN       = ",  F20.10)') LWDN
 876                 write(*,'(10x, "PRCP       = ",  F20.10)') PRCP
 877                 write(*,'(10x, "ZLVL       = ",  F20.10)') ZLVL
 878                 write(*,'(10x, "CO2AIR     = ",  F20.10)') CO2AIR
 879                 write(*,'(10x, "O2AIR      = ",  F20.10)') O2AIR
 880                 write(*,'(10x, "COSZ       = ",  F20.10)') COSZ
 881                 write(*,'(10x, "TBOT       = ",  F20.10)') TBOT
 882                 write(*,'(10x, "FOLN       = ",  F20.10)') FOLN
 883                 write(*,'(10x, "SFCPRS     = ",  F20.10)') SFCPRS
 884                 write(*,'(10x, "SHDFAC     = ",  F20.10)') SHDFAC
 885                 write(*,'(10x, "LAT        = ",  F20.10)') LAT
 886                 write(*,'(10x, "DZ8W1D     = ",  F20.10)') DZ8W1D
 887                 write(*,'(10x, "EAH        = ",  F20.10)') EAH
 888                 write(*,'(10x, "TAH        = ",  F20.10)') TAH
 889                 write(*,'(10x, "FWET       = ",  F20.10)') FWET
 890                 write(*,'(10x, "FICEOLD    = ", 7F20.10)') FICEOLD
 891                 write(*,'(10x, "QSNOW      = ",  F20.10)') QSNOW
 892                 write(*,'(10x, "SNEQVO     = ",  F20.10)') SNEQVO
 893                 write(*,'(10x, "ISNOW      = ",  F20.10)') ISNOW
 894                 write(*,'(10x, "ZSNSO      = ", 7F20.10)') ZSNSO
 895                 write(*,'(10x, "CANLIQ     = ",  F20.10)') CANLIQ
 896                 write(*,'(10x, "CANICE     = ",  F20.10)') CANICE
 897                 write(*,'(10x, "SNOWD      = ",  F20.10)') SNOWD
 898                 write(*,'(10x, "SWE        = ",  F20.10)') SWE
 899                 write(*,'(10x, "SNICE      = ", 7F20.10)') SNICE
 900                 write(*,'(10x, "SNLIQ      = ", 7F20.10)') SNLIQ
 901                 write(*,'(10x, "TV         = ",  F20.10)') TV
 902                 write(*,'(10x, "TG         = ",  F20.10)') TG
 903                 write(*,'(10x, "STC        = ", 7F20.10)') STC
 904                 write(*,'(10x, "SMH2O      = ", 7F20.10)') SMH2O
 905                 write(*,'(10x, "SMC        = ", 7F20.10)') SMC
 906                 write(*,'(10x, "ZWT        = ",  F20.10)') ZWT
 907                 write(*,'(10x, "WA         = ",  F20.10)') WA
 908                 write(*,'(10x, "WT         = ",  F20.10)') WT
 909                 write(*,'(10x, "WSLAKE     = ",  F20.10)') WSLAKE
 910                 write(*,'(10x, "LFMASS     = ",  F20.10)') LFMASS
 911                 write(*,'(10x, "RTMASS     = ",  F20.10)') RTMASS
 912                 write(*,'(10x, "STMASS     = ",  F20.10)') STMASS
 913                 write(*,'(10x, "WOOD       = ",  F20.10)') WOOD
 914                 write(*,'(10x, "STBLCP     = ",  F20.10)') STBLCP
 915                 write(*,'(10x, "FASTCP     = ",  F20.10)') FASTCP
 916                 write(*,'(10x, "PLAI       = ",  F20.10)') PLAI
 917                 write(*,'(10x, "PSAI       = ",  F20.10)') PSAI
 918                 write(*,'(10x, "ALBOLD     = ",  F20.10)') ALBOLD
 919                 write(*,'(10x, "CM         = ",  F20.10)') CM
 920                 write(*,'(10x, "CH         = ",  F20.10)') CH
 921                 write(*,'(10x, "DX         = ",  F20.10)') DX
 922                 write(*,'(10x, "ISURBAN    = ",  I10   )') ISURBAN
 923                 write(*,'(10x, "IZ0TLND    = ",  I10   )') IZ0TLND
 924                 write(*,'(10x, "QC         = ",  F20.10)') QC
 925                 write(*,'(10x, "PBLH       = ",  F20.10)') PBLH
 926                 write(*,'(10x, "QSFC1D     = ",  F20.10)') QSFC1D
 927                 write(*,'(10x, "PSFC       = ",  F20.10)') PSFC
 928              endif
 929 #endif
 930
 931              CALL NOAHMP_SFLX (&
 932                   I       , J       , LAT     , YEARLEN , JULIAN  , COSZ    , & ! IN : Time/Space-related
 933                   DT      , DX      , DZ8W1D  , NSOIL   , ZSOIL   , 3       , & ! IN : Model configuration
 934                   SHDFAC  , SHDMAX  , VEGTYP  , ISURBAN , ICE     , IST     , & ! IN : Vegetation/Soil characteristics
 935                   ISC     ,                                                   & ! IN : Vegetation/Soil characteristics
 936                   IZ0TLND ,                                                   & ! IN : User options
 937                   SFCTMP  , SFCPRS  , PSFC    , UU      , VV      , QV1D    , & ! IN : Forcing
 938                   QC      , SOLDN   , LWDN    , PRCP    , TBOT    , CO2AIR  , & ! IN : Forcing
 939                   O2AIR   , FOLN    , FICEOLD , PBLH   ,                      & ! IN : Forcing
 940                   ZLVL    , ALBOLD  , SNEQVO  ,                               & ! IN/OUT :
 941                   STC     , SMH2O   , SMC     , TAH     , EAH     , FWET    , & ! IN/OUT :
 942                   CANLIQ  , CANICE  , TV      , TG      , QSFC1D  , QSNOW   , & ! IN/OUT :
 943                   ISNOW   , ZSNSO   , SNOWD   , SWE     , SNICE   , SNLIQ   , & ! IN/OUT :
 944                   ZWT     , WA      , WT      , WSLAKE  , LFMASS  , RTMASS  , & ! IN/OUT :
 945                   STMASS  , WOOD    , STBLCP  , FASTCP  , PLAI    , PSAI    , & ! IN/OUT :
 946                   CM      , CH      , CHSTAR  ,                               & ! IN/OUT :
 947                   FSA     , FSR     , FIRA    , FSH     , SSOIL   , FCEV    , & ! OUT :
 948                   FGEV    , FCTR    , ECAN    , ETRAN   , ESOIL   , TRAD    , & ! OUT :
 949                   TS      , TGB     , TGV     , T2MV    , T2MB    , TSTAR   , & ! OUT :
 950                   Q1      , Q2MV    , Q2MB    , RUNSF   , RUNSB   , APAR    , & ! OUT :
 951                   PSN     , SAV     , SAG     , FSNO    , NEE     , GPP     , & ! OUT :
 952                   NPP     , FVEG    , SALB    , QMELT   , PONDING , PONDING1, & ! OUT :
 953                   PONDING2, RSSUN   , RSSHA   , BGAP    , WGAP    , GAP     , & ! OUT :
 954                   ERRWAT  , CHV     , CHB     , EMISSI)                         ! OUT :
 955
 956 #if 0
 957             if ( I == 15 .and. J == 5 ) then
 958                 ! Intent (OUT) or Intent (INOUT), but not Intent (IN)
 959                 write(*,'("After  call to NOAHMP_SFLX, at point ", I8, I8)') i, j
 960                 write(*,'(10x, "ZLVL       = ", 7F20.10)') ZLVL
 961                 write(*,'(10x, "EAH        = ",  F20.10)') EAH
 962                 write(*,'(10x, "TAH        = ",  F20.10)') TAH
 963                 write(*,'(10x, "FWET       = ",  F20.10)') FWET
 964                 write(*,'(10x, "QSNOW      = ",  F20.10)') QSNOW
 965                 write(*,'(10x, "SNEQVO     = ",  F20.10)') SNEQVO
 966                 write(*,'(10x, "ISNOW      = ",  F20.10)') ISNOW
 967                 write(*,'(10x, "ZSNSO      = ", 7F20.10)') ZSNSO
 968                 write(*,'(10x, "CANLIQ     = ",  F20.10)') CANLIQ
 969                 write(*,'(10x, "CANICE     = ",  F20.10)') CANICE
 970                 write(*,'(10x, "SNOWD      = ",  F20.10)') SNOWD
 971                 write(*,'(10x, "SWE        = ",  F20.10)') SWE
 972                 write(*,'(10x, "SNICE      = ", 3F20.10)') SNICE
 973                 write(*,'(10x, "SNLIQ      = ", 3F20.10)') SNLIQ
 974                 write(*,'(10x, "TV         = ",  F20.10)') TV
 975                 write(*,'(10x, "TG         = ",  F20.10)') TG
 976                 write(*,'(10x, "STC        = ", 7F20.10)') STC
 977                 write(*,'(10x, "SMH2O      = ", 7F20.10)') SMH2O
 978                 write(*,'(10x, "SMC        = ", 7F20.10)') SMC
 979                 write(*,'(10x, "ZWT        = ",  F20.10)') ZWT
 980                 write(*,'(10x, "WA         = ",  F20.10)') WA
 981                 write(*,'(10x, "WT         = ",  F20.10)') WT
 982                 write(*,'(10x, "WSLAKE     = ",  F20.10)') WSLAKE
 983                 write(*,'(10x, "LFMASS     = ",  F20.10)') LFMASS
 984                 write(*,'(10x, "RTMASS     = ",  F20.10)') RTMASS
 985                 write(*,'(10x, "STMASS     = ",  F20.10)') STMASS
 986                 write(*,'(10x, "WOOD       = ",  F20.10)') WOOD
 987                 write(*,'(10x, "STBLCP     = ",  F20.10)') STBLCP
 988                 write(*,'(10x, "FASTCP     = ",  F20.10)') FASTCP
 989                 write(*,'(10x, "PLAI       = ",  F20.10)') PLAI
 990                 write(*,'(10x, "PSAI       = ",  F20.10)') PSAI
 991                 write(*,'(10x, "ALBOLD     = ",  F20.10)') ALBOLD
 992                 write(*,'(10x, "CM         = ",  F20.10)') CM
 993                 write(*,'(10x, "CH         = ",  F20.10)') CH
 994                 write(*,'(10x, "FSA        = ",  F20.10)') FSA
 995                 write(*,'(10x, "FSR        = ",  F20.10)') FSR
 996                 write(*,'(10x, "FIRA       = ",  F20.10)') FIRA
 997                 write(*,'(10x, "FSH        = ",  F20.10)') FSH
 998                 write(*,'(10x, "SSOIL      = ",  F20.10)') SSOIL
 999                 write(*,'(10x, "FCEV       = ",  F20.10)') FCEV
1000                 write(*,'(10x, "FGEV       = ",  F20.10)') FGEV
1001                 write(*,'(10x, "FCTR       = ",  F20.10)') FCTR
1002                 write(*,'(10x, "TRAD       = ",  F20.10)') TRAD
1003                 write(*,'(10x, "ECAN       = ",  F20.10)') ECAN
1004                 write(*,'(10x, "ETRAN      = ",  F20.10)') ETRAN
1005                 write(*,'(10x, "ESOIL      = ",  F20.10)') ESOIL
1006                 write(*,'(10x, "RUNSF      = ",  F20.10)') RUNSF
1007                 write(*,'(10x, "RUNSB      = ",  F20.10)') RUNSB
1008                 write(*,'(10x, "APAR       = ",  F20.10)') APAR
1009                 write(*,'(10x, "PSN        = ",  F20.10)') PSN
1010                 write(*,'(10x, "SAV        = ",  F20.10)') SAV
1011                 write(*,'(10x, "SAG        = ",  F20.10)') SAG
1012                 write(*,'(10x, "FSNO       = ",  F20.10)') FSNO
1013                 write(*,'(10x, "NEE        = ",  F20.10)') NEE
1014                 write(*,'(10x, "GPP        = ",  F20.10)') GPP
1015                 write(*,'(10x, "NPP        = ",  F20.10)') NPP
1016                 write(*,'(10x, "TS         = ",  F20.10)') TS
1017                 write(*,'(10x, "FVEG       = ",  F20.10)') FVEG
1018                 write(*,'(10x, "SALB       = ",  F20.10)') SALB
1019                 write(*,'(10x, "ERRWAT     = ",  F20.10)') ERRWAT
1020                 write(*,'(10x, "QMELT      = ",  F20.10)') QMELT
1021                 write(*,'(10x, "PONDING    = ",  F20.10)') PONDING
1022                 write(*,'(10x, "PONDING1   = ",  F20.10)') PONDING1
1023                 write(*,'(10x, "PONDING2   = ",  F20.10)') PONDING2
1024                 write(*,'(10x, "QSFC1D     = ",  F20.10)') QSFC1D
1025                 write(*,'(10x, "CHSTAR     = ",  F20.10)') CHSTAR
1026                 write(*,'(10x, "TSTAR      = ",  F20.10)') TSTAR
1027                 write(*,'(10x, "T2MV       = ",  F20.10)') T2MV
1028                 write(*,'(10x, "T2MB       = ",  F20.10)') T2MB
1029                 write(*,'(10x, "RSSUN      = ",  F20.10)') RSSUN
1030                 write(*,'(10x, "RSSHA      = ",  F20.10)') RSSHA
1031                 write(*,'(10x, "BGAP       = ",  F20.10)') BGAP
1032                 write(*,'(10x, "WGAP       = ",  F20.10)') WGAP
1033                 write(*,'(10x, "GAP        = ",  F20.10)') GAP
1034                 write(*,'(10x, "TGV        = ",  F20.10)') TGV
1035                 write(*,'(10x, "TGB        = ",  F20.10)') TGB
1036                 write(*,'(10x, "Q1         = ",  F20.10)') Q1
1037              endif
1038 #endif
1039
1040              !Q1           = eah * 0.622 / (SFCPRS - 0.378*eah)
1041
1042              chs2     (i,j)               = chstar2
1043              cqs2     (i,j)               = cqstar2
1044              QFX      (I,J)               = ecan + esoil + etran
1045              SNOWC    (I,J)               = fsno
1046
1047    ENDIF ! glacial, seaice split ends
1048 !jref:end
1049
1050              isnowxy  (i,j)               = isnow
1051              canliqxy (i,j)               = canliq
1052              canicexy (i,j)               = canice
1053              snowh    (i,j)               = snowd
1054              snow     (i,j)               = swe
1055              zsnsoxy  (i,isnow+1:nsoil,j) = zsnso (isnow+1:nsoil)
1056              tslb     (i,      1:nsoil,j) = stc   (      1:nsoil)
1057              tsnoxy   (i,isnow+1:    0,j) = stc   (isnow+1:    0)
1058              smois    (i,      1:nsoil,j) = smc   (      1:nsoil)
1059              sh2o     (i,      1:nsoil,j) = smh2o (      1:nsoil)
1060              snicexy  (i,isnow+1:    0,j) = snice (isnow+1:    0)
1061              snliqxy  (i,isnow+1:    0,j) = snliq (isnow+1:    0)
1062              tvxy     (i,j)               = tv
1063              tgxy     (i,j)               = tg
1064              zwtxy    (i,j)               = zwt
1065              waxy     (i,j)               = wa
1066              wtxy     (i,j)               = wt
1067              lfmassxy (i,j)               = lfmass
1068              rtmassxy (i,j)               = rtmass
1069              stmassxy (i,j)               = stmass
1070              woodxy   (i,j)               = wood
1071              stblcpxy (i,j)               = stblcp
1072              fastcpxy (i,j)               = fastcp
1073              xlaixy   (i,j)               = plai
1074              xsaixy   (i,j)               = psai
1075              emiss    (i,j)               = emissi
1076              eahxy    (i,j)               = eah
1077              tahxy    (i,j)               = tah
1078              fwetxy   (i,j)               = fwet
1079              sneqvoxy (i,j)               = sneqvo
1080              alboldxy (i,j)               = albold
1081              qsnowxy  (i,j)               = qsnow
1082              wslakexy (i,j)               = wslake
1083              cmxy     (i,j)               = cm
1084 !jref:start
1085              chxy     (i,j)               = chstar
1086              rssunxy  (i,j)               = rssun
1087              rsshaxy  (i,j)               = rssha
1088              bgapxy   (i,j)               = bgap
1089              wgapxy   (i,j)               = wgap
1090              gapxy    (i,j)               = gap
1091              tgvxy    (i,j)               = tgv
1092              tgbxy    (i,j)               = tgb
1093              chvxy    (i,j)               = chv
1094              chbxy    (i,j)               = chb
1095 !jref:end
1096
1097              !for output
1098
1099              runsfxy  (i,j)               = runsf
1100              runsbxy  (i,j)               = runsb
1101              ecanxy   (i,j)               = ecan
1102              edirxy   (i,j)               = esoil
1103              etranxy  (i,j)               = etran
1104              aparxy   (i,j)               = apar
1105              psnxy    (i,j)               = psn
1106              savxy    (i,j)               = sav
1107              sagxy    (i,j)               = sag
1108              fsnoxy   (i,j)               = fsno
1109              fsaxy    (i,j)               = fsa
1110              firaxy   (i,j)               = fira
1111              hfx      (i,j)               = fsh
1112              lh       (i,j)               = fcev + fgev + fctr
1113              grdflx   (i,j)               = ssoil
1114              tradxy   (i,j)               = trad
1115              tsxy     (i,j)               = ts
1116              neexy    (i,j)               = nee
1117              gppxy    (i,j)               = gpp
1118              nppxy    (i,j)               = npp
1119              fvegxy   (i,j)               = fveg
1120 !jref:4/21/2011
1121              t2mvxy   (i,j)               = t2mv
1122              t2mbxy   (i,j)               = t2mb
1123              q2mvxy   (i,j)               = q2mv
1124              q2mbxy   (i,j)               = q2mb
1125              chstarxy (i,j)               = chstar
1126              chs      (i,j)               = chstar
1127              tstarxy  (i,j)               = tstar
1128 !jref:4/21/2011
1129
1130              CANWAT(I,J)  = canliqxy (i,j) + canicexy (i,j)
1131              IF ( SALB > -999 ) THEN
1132                 ALBEDO(I,J)  = salb
1133              ENDIF
1134              TSK(I,J)     = tradxy   (i,j)
1135 !KWM          TSK(I,J)     = tstarxy  (i,j)
1136 !niu          POTEVP(I,J)  = ???
1137 !jref         CHS2(I,J)    = chxy     (i,j)
1138
1139              !IF (Q1.GT.QSFC(I,J)) THEN
1140              !   CQS2(I,J) = CHS(I,J)
1141              !END IF
1142              QSFC(I,J)    = Q1/(1.0-Q1)
1143 !jref: specific humidity to mixing ratio
1144              q2mvxy(i,j)  = q2mvxy(i,j)/(1.0-q2mvxy(i,j))
1145 !             IF (VEGTYP == ISURBAN) write(*,*) "IN SFCDRV: q2mb=",q2mb,"q2mbxy(i,j)=",q2mbxy(i,j)
1146              q2mbxy(i,j)  = q2mbxy(i,j)/(1.0-q2mbxy(i,j))
1147
1148 !***  DIAGNOSTICS
1149 !jref:start - THESE SHOULD BE LOOKED AT!!!
1150              SNOTIME(I,J) = SNOTIME1
1151              SMSTAV(I,J)=SOILW
1152              SMSTOT(I,J)=SOILM*1000.
1153 !         Convert the water unit into mm
1154              SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+runsfxy(i,j)*DT*1000.0
1155              UDRUNOFF(I,J)=UDRUNOFF(I,J)+runsbxy(i,j)*DT*1000.0
1156 !jref        SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
1157 !jref        UDRUNOFF(I,J)=UDRUNOFF(I,J)+(RUNOFF2+RUNOFF3)*DT*1000.0
1158 !jref:end
1159 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
1160              IF(FFROZP.GT.0.5)THEN
1161                 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
1162              ENDIF
1163              IF(SNOW(I,J).GT.0.)THEN
1164 !KWM                ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
1165              ENDIF
1166
1167           ENDIF                                                           ! endif of land-sea test
1168 !jref:start  make sure exchange coeff and TSK include water points
1169 !          IF((XLAND(I,J)-1.5).GE.0.)THEN                                  ! begining of land/sea if block
1170 !             chstar2xy(i,j) = CHS2(i,j)
1171 !             chstarxy(i,j) = CHS(i,j)
1172 !             tstarxy(i,j)  = T1 !TSK(i,j) test with T1
1173 !          ENDIF
1174 !jref:end
1175
1176
1177        ENDDO
1178
1179     ENDDO                                                                ! of J loop
1180 !------------------------------------------------------
1181   END SUBROUTINE noahmplsm
1182 !------------------------------------------------------
1183
1184   SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP ,           &
1185        TSLB , SMOIS , SH2O , DZS , FNDSOILW , FNDSNOWH ,                      &
1186        TSK, isnowxy , tvxy     ,tgxy     ,canicexy ,                          &
1187        canliqxy ,eahxy    ,tahxy    ,cmxy     ,chxy     ,                     &
1188        fwetxy   ,sneqvoxy ,alboldxy ,qsnowxy  ,wslakexy ,zwtxy    ,waxy     , &
1189        wtxy     ,tsnoxy   ,zsnsoxy  ,snicexy  ,snliqxy  ,lfmassxy ,rtmassxy , &
1190        stmassxy ,woodxy   ,stblcpxy ,fastcpxy ,xsaixy   , &
1191 !jref:start
1192        t2mvxy   ,t2mbxy   ,chstarxy ,            &
1193 !jref:end
1194        num_soil_layers, restart,                 &
1195        allowed_to_read ,                         &
1196        ids,ide, jds,jde, kds,kde,                &
1197        ims,ime, jms,jme, kms,kme,                &
1198        its,ite, jts,jte, kts,kte                 )
1199
1200 ! Initializing Canopy air temperature to 287 K seems dangerous to me [KWM].
1201
1202     INTEGER, INTENT(IN   )    ::     ids,ide, jds,jde, kds,kde,  &
1203          &                           ims,ime, jms,jme, kms,kme,  &
1204          &                           its,ite, jts,jte, kts,kte
1205
1206     INTEGER, INTENT(IN)       ::     num_soil_layers
1207
1208     LOGICAL, INTENT(IN)       ::     restart,                    &
1209          &                           allowed_to_read
1210
1211     REAL,    DIMENSION( num_soil_layers), INTENT(IN)    ::     DZS  ! Thickness of the soil layers [m]
1212
1213     REAL,    DIMENSION( ims:ime, num_soil_layers, jms:jme ) ,    &
1214          &   INTENT(INOUT)    ::     SMOIS,                      &
1215          &                           SH2O,                       &
1216          &                           TSLB
1217
1218     REAL,    DIMENSION( ims:ime, jms:jme ) ,                     &
1219          &   INTENT(INOUT)    ::     SNOW,                       &
1220          &                           SNOWH,                      &
1221          &                           CANWAT
1222
1223     INTEGER, DIMENSION( ims:ime, jms:jme ),                      &
1224          &   INTENT(IN)       ::     ISLTYP
1225
1226     LOGICAL, INTENT(IN)       ::     FNDSOILW,                   &
1227          &                           FNDSNOWH
1228
1229     REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: TSK         !skin temperature (k)
1230     INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: isnowxy     !actual no. of snow layers
1231     REAL, DIMENSION(ims:ime,-2:num_soil_layers,jms:jme), INTENT(INOUT) :: zsnsoxy  !snow layer depth [m]
1232     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: tsnoxy   !snow temperature [K]
1233     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: snicexy  !snow layer ice [mm]
1234     REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: snliqxy  !snow layer liquid water [mm]
1235     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tvxy        !vegetation canopy temperature
1236     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tgxy        !ground surface temperature
1237     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canicexy    !canopy-intercepted ice (mm)
1238     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canliqxy    !canopy-intercepted liquid water (mm)
1239     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: eahxy       !canopy air vapor pressure (pa)
1240     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tahxy       !canopy air temperature (k)
1241     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: cmxy        !momentum drag coefficient
1242     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chxy        !sensible heat exchange coefficient
1243     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fwetxy      !wetted or snowed fraction of the canopy (-)
1244     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: sneqvoxy    !snow mass at last time step(mm h2o)
1245     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: alboldxy    !snow albedo at last time step (-)
1246     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: qsnowxy     !snowfall on the ground [mm/s]
1247     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wslakexy    !lake water storage [mm]
1248     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: zwtxy       !water table depth [m]
1249     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: waxy        !water in the "aquifer" [mm]
1250     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wtxy        !groundwater storage [mm]
1251     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: lfmassxy    !leaf mass [g/m2]
1252     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: rtmassxy    !mass of fine roots [g/m2]
1253     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stmassxy    !stem mass [g/m2]
1254     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: woodxy      !mass of wood (incl. woody roots) [g/m2]
1255     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stblcpxy    !stable carbon in deep soil [g/m2]
1256     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fastcpxy    !short-lived carbon, shallow soil [g/m2]
1257     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xsaixy      !stem area index
1258
1259 !jref:start
1260     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: t2mvxy        !2m temperature vegetation part (k)
1261     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: t2mbxy        !2m temperature bare ground part (k)
1262     REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chstarxy    !effective exchange coefficient
1263 !jref:end
1264
1265
1266     REAL, DIMENSION(1:num_soil_layers)  :: ZSOIL      ! Depth of the soil layer bottom (m) from
1267     !                                                   the surface (negative)
1268
1269     REAL                      :: BX, SMCMAX, PSISAT, FREE
1270
1271     REAL, PARAMETER           :: BLIM  = 5.5
1272     REAL, PARAMETER           :: HLICE = 3.335E5
1273     REAL, PARAMETER           :: GRAV = 9.81
1274     REAL, PARAMETER           :: T0 = 273.15
1275
1276     INTEGER                   :: errflag
1277
1278     character(len=80) :: err_message
1279     character(len=4)  :: MMINSL
1280     character(len=*), intent(in) :: MMINLU
1281     MMINSL='STAS'
1282
1283     call read_mp_veg_parameters(trim(MMINLU))
1284
1285     !
1286     ! initialize three Noah LSM related tables
1287     !
1288     IF ( allowed_to_read ) THEN
1289        CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
1290        CALL  SOIL_VEG_GEN_PARM( MMINLU, MMINSL )
1291     ENDIF
1292
1293     IF( .NOT. restart ) THEN
1294
1295        itf=min0(ite,ide-1)
1296        jtf=min0(jte,jde-1)
1297
1298        errflag = 0
1299        DO j = jts,jtf
1300           DO i = its,itf
1301              IF ( ISLTYP( i,j ) .LT. 1 ) THEN
1302                 errflag = 1
1303                 WRITE(err_message,*)"module_sf_noahlsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
1304                 CALL wrf_message(err_message)
1305              ENDIF
1306           ENDDO
1307        ENDDO
1308        IF ( errflag .EQ. 1 ) THEN
1309           CALL wrf_error_fatal( "module_sf_noahlsm.F: lsminit: out of range value "// &
1310                "of ISLTYP. Is this field in the input?" )
1311        ENDIF
1312 #ifdef WRF_CHEM
1313        !
1314        ! need this parameter for dust parameterization in wrf/chem
1315        !
1316        do I=1,NSLTYPE
1317           porosity(i)=maxsmc(i)
1318        enddo
1319 #endif
1320
1321 ! initialize soil liquid water content SH2O
1322
1323 !  IF(.NOT.FNDSOILW) THEN
1324
1325 ! If no SWC, do the following
1326 !         PRINT *,'SOIL WATER NOT FOUND - VALUE SET IN LSMINIT'
1327        DO J = jts , jtf
1328           DO I = its , itf
1329              BX = BB(ISLTYP(I,J))
1330              SMCMAX = MAXSMC(ISLTYP(I,J))
1331              PSISAT = SATPSI(ISLTYP(I,J))
1332              IF ( ( bx > 0.0 ) .AND. ( smcmax > 0.0 ) .AND. ( psisat > 0.0 ) ) THEN
1333                 DO NS=1, num_soil_layers
1334
1335                    IF ( TSLB(I,NS,J) < 273.149 ) THEN
1336                       ! SH2O <= SMOIS for T < 273.149K (-0.001C)
1337
1338                       ! First guess of SH2O following explicit solution for
1339                       ! Flerchinger Eqn from Koren et al, JGR, 1999, Eqn 17
1340                       ! (KCOUNT=0 in function FRH2O).
1341                       BX = BB(ISLTYP(I,J))
1342                       SMCMAX = MAXSMC(ISLTYP(I,J))
1343                       PSISAT = SATPSI(ISLTYP(I,J))
1344                       IF ( BX >  BLIM ) BX = BLIM
1345                       FK=(( (HLICE/(GRAV*(-PSISAT))) *                              &
1346                            ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
1347                       FK = MAX(FK, 0.02)
1348                       SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
1349
1350                       ! Use iterative solution for liquid soil water content
1351                       ! using function FRH2O, with the initial guess for SH2O
1352                       ! from the above explicit first guess.
1353                       CALL FRH2O ( FREE , TSLB(I,NS,J) , SMOIS(I,NS,J) , SH2O(I,NS,J) )
1354                       SH2O(I,NS,J) = FREE
1355                    ELSE
1356                       ! SH2O = SMOIS ( for T => 273.149K (-0.001C)
1357                       SH2O(I,NS,J)=SMOIS(I,NS,J)
1358                    ENDIF
1359                 END DO
1360              ELSE
1361                 DO NS=1, num_soil_layers
1362                    SH2O(I,NS,J)=SMOIS(I,NS,J)
1363                 END DO
1364              ENDIF
1365           ENDDO
1366        ENDDO
1367 !  ENDIF
1368
1369
1370        !
1371        ! initialize physical snow height SNOWH
1372        !
1373        IF(.NOT.FNDSNOWH)THEN
1374           ! If no SNOWH do the following
1375           CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
1376           DO J = jts,jtf
1377              DO I = its,itf
1378                 SNOWH(I,J)=SNOW(I,J)*0.005               ! SNOW in mm and SNOWH in m
1379              ENDDO
1380           ENDDO
1381        ENDIF
1382
1383        DO J = jts,jtf
1384           DO I = its,itf
1385              tvxy       (I,J) = TSK(I,J)
1386              tgxy       (I,J) = TSK(I,J)
1387              CANWAT     (I,J) = 0.0
1388              canliqxy   (I,J) = CANWAT(I,J)
1389              canicexy   (I,J) = 0.
1390              eahxy      (I,J) = 2000.
1391              tahxy      (I,J) = 287.
1392 !jref:start
1393              t2mvxy     (I,J) = TSK(I,J)
1394              t2mbxy     (I,J) = TSK(I,J)
1395              chstarxy   (I,J) = 0.0
1396 !jref:end
1397
1398              cmxy       (I,J) = 0.0
1399              chxy       (I,J) = 0.0
1400              fwetxy     (I,J) = 0.0
1401              sneqvoxy   (I,J) = 0.0
1402              alboldxy   (I,J) = 0.65
1403              qsnowxy    (I,J) = 0.0
1404              wslakexy   (I,J) = 0.0
1405
1406              waxy       (I,J) = 4900.                                       !???
1407              wtxy       (I,J) = waxy(i,j)                                   !???
1408              zwtxy      (I,J) = (25. + 2.0) - waxy(i,j)/1000/0.2            !???
1409
1410              lfmassxy   (I,J) = 50.         !
1411              stmassxy   (I,J) = 50.0        !
1412              rtmassxy   (I,J) = 500.0       !
1413              woodxy     (I,J) = 500.0       !
1414              stblcpxy   (I,J) = 1000.0      !
1415              fastcpxy   (I,J) = 1000.0      !
1416              xsaixy     (I,J) = 0.1         !
1417
1418           enddo
1419        enddo
1420
1421        ! Given the soil layer thicknesses (in DZS), initialize the soil layer
1422        ! depths from the surface.
1423        ZSOIL(1)         = -DZS(1)          ! negative
1424        DO NS=2, num_soil_layers
1425           ZSOIL(NS)       = ZSOIL(NS-1) - DZS(NS)
1426        END DO
1427
1428        ! Initialize snow/soil layer arrays ZSNSOXY, TSNOXY, SNICEXY, SNLIQXY,
1429        ! and ISNOWXY
1430        CALL snow_init ( ims , ime , jms , jme , its , itf , jts , jtf , 3 , &
1431             &           num_soil_layers , zsoil , snow , tgxy , snowh ,     &
1432             &           zsnsoxy , tsnoxy , snicexy , snliqxy , isnowxy )
1433
1434     ENDIF
1435   END SUBROUTINE NOAHMP_INIT
1436
1437 !------------------------------------------------------------------------------------------
1438 !------------------------------------------------------------------------------------------
1439
1440   SUBROUTINE SNOW_INIT ( ims , ime , jms , jme , its , itf , jts , jtf ,                  &
1441        &                 NSNOW , NSOIL , ZSOIL , SWE , TGXY , SNODEP ,                    &
1442        &                 ZSNSOXY , TSNOXY , SNICEXY ,SNLIQXY , ISNOWXY )
1443 !------------------------------------------------------------------------------------------
1444 !   Initialize snow arrays for Noah-MP LSM, based in input SNOWDEP, NSNOW
1445 !   ISNOWXY is an index array, indicating the index of the top snow layer.  Valid indices
1446 !           for snow layers range from 0 (no snow) and -1 (shallow snow) to (-NSNOW)+1 (deep snow).
1447 !   TSNOXY holds the temperature of the snow layer.  Snow layers are initialized with
1448 !          temperature = ground temperature [?].  Snow-free levels in the array have value 0.0
1449 !   SNICEXY is the frozen content of a snow layer.  Initial estimate based on SNODEP and SWE
1450 !   SNLIQXY is the liquid content of a snow layer.  Initialized to 0.0
1451 !   ZNSNOXY is the layer depth from the surface.
1452 !------------------------------------------------------------------------------------------
1453     IMPLICIT NONE
1454 !------------------------------------------------------------------------------------------
1455     INTEGER, INTENT(IN)                              :: ims, ime, jms, jme
1456     INTEGER, INTENT(IN)                              :: its, itf, jts, jtf
1457     INTEGER, INTENT(IN)                              :: NSNOW
1458     INTEGER, INTENT(IN)                              :: NSOIL
1459     REAL,    INTENT(IN), DIMENSION(ims:ime, jms:jme) :: SWE
1460     REAL,    INTENT(IN), DIMENSION(ims:ime, jms:jme) :: SNODEP
1461     REAL,    INTENT(IN), DIMENSION(ims:ime, jms:jme) :: TGXY
1462     REAL,    INTENT(IN), DIMENSION(1:NSOIL)          :: ZSOIL
1463
1464     INTEGER, INTENT(OUT), DIMENSION(ims:ime, jms:jme)                :: ISNOWXY ! Top snow layer index
1465     REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:NSOIL,jms:jme) :: ZSNSOXY ! Snow/soil layer depth from surface [m]
1466     REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:    0,jms:jme) :: TSNOXY  ! Snow layer temperature [K]
1467     REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:    0,jms:jme) :: SNICEXY ! Snow layer ice content [mm]
1468     REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:    0,jms:jme) :: SNLIQXY ! snow layer liquid content [mm]
1469
1470 ! Local variables:
1471 !   DZSNO   holds the thicknesses of the various snow layers.
1472 !   DZSNOSO holds the thicknesses of the various soil/snow layers.
1473     INTEGER                           :: I,J,IZ
1474     REAL,   DIMENSION(-NSNOW+1:    0) :: DZSNO
1475     REAL,   DIMENSION(-NSNOW+1:NSOIL) :: DZSNSO
1476
1477 !------------------------------------------------------------------------------------------
1478
1479     DO J = jts , jtf
1480        DO I = its , itf
1481           IF ( SNODEP(I,J) < 0.025 ) THEN
1482              ISNOWXY(I,J) = 0
1483              DZSNO(-NSNOW+1:0) = 0.
1484           ELSE
1485              IF ( ( SNODEP(I,J) >= 0.025 ) .AND. ( SNODEP(I,J) <= 0.05 ) ) THEN
1486                 ISNOWXY(I,J)    = -1
1487                 DZSNO(0)  = SNODEP(I,J)
1488              ELSE IF ( ( SNODEP(I,J) > 0.05 ) .AND. ( SNODEP(I,J) <= 0.10 ) ) THEN
1489                 ISNOWXY(I,J)    = -2
1490                 DZSNO(-1) = SNODEP(I,J)/2.
1491                 DZSNO( 0) = SNODEP(I,J)/2.
1492              ELSE IF ( (SNODEP(I,J) > 0.10 ) .AND. ( SNODEP(I,J) <= 0.25 ) ) THEN
1493                 ISNOWXY(I,J)    = -2
1494                 DZSNO(-1) = 0.05
1495                 DZSNO( 0) = SNODEP(I,J) - DZSNO(-1)
1496              ELSE IF ( ( SNODEP(I,J) > 0.25 ) .AND. ( SNODEP(I,J) <= 0.35 ) ) THEN
1497                 ISNOWXY(I,J)    = -3
1498                 DZSNO(-2) = 0.05
1499                 DZSNO(-1) = 0.5*(SNODEP(I,J)-DZSNO(-2))
1500                 DZSNO( 0) = 0.5*(SNODEP(I,J)-DZSNO(-2))
1501              ELSE IF ( SNODEP(I,J) > 0.35 ) THEN
1502                 ISNOWXY(I,J)     = -3
1503                 DZSNO(-2) = 0.05
1504                 DZSNO(-1) = 0.10
1505                 DZSNO( 0) = SNODEP(I,J) - DZSNO(-1) - DZSNO(-2)
1506              ELSE
1507                 CALL wrf_error_fatal("Problem with the logic assigning snow layers.")
1508              END IF
1509           END IF
1510
1511           TSNOXY (I,-NSNOW+1:0,J) = 0.
1512           SNICEXY(I,-NSNOW+1:0,J) = 0.
1513           SNLIQXY(I,-NSNOW+1:0,J) = 0.
1514           DO IZ = ISNOWXY(I,J)+1 , 0
1515              TSNOXY(I,IZ,J)  = TGXY(I,J)  ! [k]
1516              SNLIQXY(I,IZ,J) = 0.00
1517              SNICEXY(I,IZ,J) = 1.00 * DZSNO(IZ) * (SWE(I,J)/SNODEP(I,J))  ! [kg/m3]
1518           END DO
1519
1520           ! Assign local variable DZSNSO, the soil/snow layer thicknesses, for snow layers
1521           DO IZ = ISNOWXY(I,J)+1 , 0
1522              DZSNSO(IZ) = -DZSNO(IZ)
1523           END DO
1524
1525           ! Assign local variable DZSNSO, the soil/snow layer thicknesses, for soil layers
1526           DZSNSO(1) = ZSOIL(1)
1527           DO IZ = 2 , NSOIL
1528              DZSNSO(IZ) = (ZSOIL(IZ) - ZSOIL(IZ-1))
1529           END DO
1530
1531           ! Assign ZSNSOXY, the layer depths, for soil and snow layers
1532           ZSNSOXY(I,ISNOWXY(I,J)+1,J) = DZSNSO(ISNOWXY(I,J)+1)
1533           DO IZ = ISNOWXY(I,J)+2 , NSOIL
1534              ZSNSOXY(I,IZ,J) = ZSNSOXY(I,IZ-1,J) + DZSNSO(IZ)
1535           ENDDO
1536
1537        END DO
1538     END DO
1539
1540   END SUBROUTINE SNOW_INIT
1541 !
1542 !------------------------------------------------------------------------------------------
1543 !------------------------------------------------------------------------------------------
1544 !
1545 END MODULE module_sf_noahmpdrv