Explicit urban vars

README.md

@@ -1,39 +1,5 @@
-# noahmp
-Noah-MP Community Repository
+# Noah-MP version 4.7 release
 
 
-This is the official Noah-MP unified Github repository for code downloading and contribution. Note that Noah-MP is a community model contributed by the whole Noah-MP scientific community. For maintenance and release of this GitHub, please contact: Cenlin He (cenlinhe@ucar.edu) and Fei Chen (feichen@ucar.edu).
+This is the official Noah-MP code version 4.7 consistent with that released in WRF v4.7. Note that WRF v4.7 GitHub code is directly connected to this Noah-MP GitHub through the submodule mechanism. 
 
-Some changes have been made to the structure of archiving the stand-alone version of Noah-MP/HRLDAS code in the Github repository. Now, it separately archives the core Noah-MP source code, Noah-MP driver, and parameter tables in this unified Noah-MP Github repository and the rest of the HRLDAS related files (e.g., module_sf_urban.F, etc.) in another unified HRLDAS Github repository (https://github.com/NCAR/hrldas). The HRLDAS Github repo is already linked to this unified core Noah-MP code repository. This new archiving structure will allow different host model platforms/systems (e.g., HRLDAS, WRF, UFS, NWM, LIS, etc.) to connect to the core Noah-MP source code and develop their own host model drivers. 
-
-Model developers can make code development based on the develop branch and create pull request to the develop branch. The pull request will be reviewed by Noah-MP model release team and if everything looks good, the new code development will be merged to the develop branch. Eventually, the updates in the develop branch will be merged to the master branch for official Noah-MP model release.
-
-Branch structures of this Noah-MP repository:
-
-1. "master" branch: most stable & recent version, updated whenever there are bug fixes or major model update/release;
-
-2. "develop" branch: used for continuous NoahMP development, keep being updated by including bug fixes and code updates (e.g., new physics options, processes, etc.); 
-
-3. "develop_no_gecros" branch: same as the "develop" branch, except for excluding the gecros crop section (to be consistent with recent Noah-MP changes for National Water Model (NWM));
-
-4. "release-v4.3-WRF" branch: version archive of the stable version consistent with the WRFv4.3 release;
-
-5. "release-v4.3-NWM" branch: version archive, same as the "release-v4.3-WRF" branch, except for excluding the gecros crop section (to be consistent with recent Noah-MP changes for National Water Model (NWM));
-
-Some suggestions for model developers to contribute to Noah-MP code through the Github repository (typical procedures):
-
-1. Step (1) Create a fork of this official NoahMP repository to your own Github account; 
-
-2. Step (2) Make code updates based on the "develop" branch of the forked repository under your own account; 
-
-3. Step (3) Finalize and test the code updates you make; 
-
-4. Step (4) Submit a pull request for your model updates from your own forked Github repository to the "develop" branch of the official repository;
-
-5. Step (5) The Noah-MP release team reviews and tests the model updates in the submitted pull request and discusses with the developer if there is any problem; 
-
-6. Step (6) The Noah-MP release team confirms the pull request and merges the updated code to the "develop" and "develop_no_gecros" branches in the official repository;
-
-7. Step (7) The Noah-MP release team will merge the updated "develop" branch to the master branch and version-release branch during the annual model release.
-
-Note: If updates are made to both the core NoahMP source codes and other HRLDAS files, two separate pull requests need to be submitted to this NoahMP repository and the HRLDAS repository (https://github.com/NCAR/hrldas), respectively, regarding the changes in the Noah-MP code files and other HRLDAS-related files. This could be done by using the same titles for the two pull requests to ensure that the submitted code changes are handled together by the release team in the two repositories.

drivers/wrf/module_sf_noahmpdrv.F

@@ -645,7 +645,6 @@ SUBROUTINE noahmplsm(ITIMESTEP,        YR,   JULIAN,   COSZIN,XLAT,XLONG, & ! IN
     IPRINT    =  .false.                     ! debug printout
 
 ! for using soil update timestep difference from noahmp main timestep
-    calculate_soil = .false.
     soil_update_steps = nint(soiltstep/DT)  ! 3600 = 1 hour
     soil_update_steps = max(soil_update_steps,1)
     if ( soil_update_steps == 1 ) then
@@ -673,7 +672,7 @@ SUBROUTINE noahmplsm(ITIMESTEP,        YR,   JULIAN,   COSZIN,XLAT,XLONG, & ! IN
      end if
     end if
 
-    if (mod(itimestep,soil_update_steps) == 0) calculate_soil = .true.
+    calculate_soil = mod(itimestep,soil_update_steps) == 0
 ! end soil timestep
 
     YEARLEN = 365                            ! find length of year for phenology (also S Hemisphere)
@@ -1056,7 +1055,7 @@ SUBROUTINE noahmplsm(ITIMESTEP,        YR,   JULIAN,   COSZIN,XLAT,XLONG, & ! IN
                                 FSA,     FSR,    FIRA,     FSH,    FGEV,   SSOIL, & ! OUT : 
                                TRAD,   ESOIL,   RUNSF,   RUNSB,     SAG,    SALB, & ! OUT :
                               QSNBOT,PONDING,PONDING1,PONDING2,    T2MB,    Q2MB, & ! OUT :
-			      EMISSI,  FPICE,    CHB2,   QMELT                    & ! OUT :
+			      EMISSI,  FPICE,    CHB2,   QMELT,   HCPCT            & ! OUT :
 #ifdef WRF_HYDRO
                               , sfcheadrt(i,j)                                      &
 #endif
@@ -2004,7 +2003,7 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP ,   IVGTYP, XLAT
     call read_mp_crop_parameters()
     call read_tiledrain_parameters()
     call read_mp_optional_parameters()
-    if(iopt_irr  >= 1) call read_mp_irrigation_parameters()
+    call read_mp_irrigation_parameters(iopt_irr)
 
     IF( .NOT. restart ) THEN
 
@@ -2181,7 +2180,11 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP ,   IVGTYP, XLAT
 
 	     lai        (I,J) = max(lai(i,j),0.05)             ! at least start with 0.05 for arbitrary initialization (v3.7)
              xsaixy     (I,J) = max(0.1*lai(I,J),0.05)         ! MB: arbitrarily initialize SAI using input LAI (v3.7)
-             masslai = 1000. / max(SLA_TABLE(IVGTYP(I,J)),1.0) ! conversion from lai to mass  (v3.7)
+             if (urbanpt_flag) then
+                 masslai = 1000. / max(SLA_TABLE(NATURAL_TABLE),1.0) ! conversion from lai to mass  (v3.7)
+             else
+                 masslai = 1000. / max(SLA_TABLE(IVGTYP(I,J)),1.0) ! conversion from lai to mass  (v3.7)
+             endif 
              lfmassxy   (I,J) = lai(i,j)*masslai               ! use LAI to initialize (v3.7)
              masssai = 1000. / 3.0                             ! conversion from lai to mass (v3.7)
              stmassxy   (I,J) = xsaixy(i,j)*masssai            ! use SAI to initialize (v3.7)
@@ -3035,6 +3038,9 @@ SUBROUTINE noahmp_urban(sf_urban_physics,   NSOIL,         IVGTYP,  ITIMESTEP,
                           dgr_urb3d,       dg_urb3d,      lfr_urb3d,  lfg_urb3d,            & !RMS
                            lp_urb2d,       hi_urb2d,       lb_urb2d,  hgt_urb2d,            & !H multi-layer urban
                            mh_urb2d,     stdh_urb2d,       lf_urb2d,                        & !SLUCM
+                         lf_urb2d_s,       z0_urb2d,         vegfra,                        & !SLUCM
+                         albr_urb2d,     albb_urb2d,     albg_urb2d,                        & !I urban explicit radiation PK2025
+                         epsr_urb2d,     epsb_urb2d,     epsg_urb2d,                        & !I urban explicit radiation PK2025
                              th_phy,            rho,          p_phy,        ust,            & !I multi-layer urban
                                 gmt,         julday,          xlong,       xlat,            & !I multi-layer urban
                             a_u_bep,        a_v_bep,        a_t_bep,    a_q_bep,            & !O multi-layer urban
@@ -3172,6 +3178,15 @@ SUBROUTINE noahmp_urban(sf_urban_physics,   NSOIL,         IVGTYP,  ITIMESTEP,
      REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TGRL_URB3D
      REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: SMR_URB3D
 
+!  Explicit radiation parameters PK2025
+     real :: albr_urb,albb_urb,albg_urb,epsr_urb,epsb_urb,epsg_urb
+     real, optional, dimension( ims:ime, jms:jme ), intent(in) :: albr_urb2d
+     real, optional, dimension( ims:ime, jms:jme ), intent(in) :: albb_urb2d
+     real, optional, dimension( ims:ime, jms:jme ), intent(in) :: albg_urb2d
+     real, optional, dimension( ims:ime, jms:jme ), intent(in) :: epsr_urb2d
+     real, optional, dimension( ims:ime, jms:jme ), intent(in) :: epsb_urb2d
+     real, optional, dimension( ims:ime, jms:jme ), intent(in) :: epsg_urb2d
+
 
 ! state variable surface_driver <--> lsm <--> urban
 
@@ -3239,6 +3254,11 @@ SUBROUTINE noahmp_urban(sf_urban_physics,   NSOIL,         IVGTYP,  ITIMESTEP,
      REAL :: hgt_urb
      REAL, DIMENSION(4) :: lf_urb
 
+! Distributed aerodynamics parameters
+     REAL :: lf_urb_s
+     REAL :: z0_urb
+     REAL :: vegfrac
+
 ! Local variables
 
      INTEGER :: I,J,K
@@ -3277,6 +3297,9 @@ SUBROUTINE noahmp_urban(sf_urban_physics,   NSOIL,         IVGTYP,  ITIMESTEP,
      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ),                     INTENT(IN   ) :: mh_urb2d
      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ),                     INTENT(IN   ) :: stdh_urb2d
      REAL, OPTIONAL, DIMENSION( ims:ime, 4, jms:jme ),                  INTENT(IN   ) :: lf_urb2d
+     REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ),                     INTENT(IN   ) :: lf_urb2d_s
+     REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ),                     INTENT(IN   ) :: z0_urb2d
+     REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ),                     INTENT(IN   ) :: vegfra
 
      REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ),    INTENT(INOUT) :: trb_urb4d
      REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ),    INTENT(INOUT) :: tw1_urb4d
@@ -3472,6 +3495,18 @@ SUBROUTINE noahmp_urban(sf_urban_physics,   NSOIL,         IVGTYP,  ITIMESTEP,
     IF (I.EQ.73.AND.J.EQ.125)THEN
       CHECK = 1
     END IF
+! Distributed aerodynamics
+    lf_urb_s = lf_urb2d_s(I, J)
+    z0_urb = z0_urb2d(I, J)
+    vegfrac = vegfra(I, J) / 100
+
+! Explicit radiation parameters PK2025
+    albr_urb = albr_urb2d(i,j)
+    albb_urb = albb_urb2d(i,j)
+    albg_urb = albg_urb2d(i,j)
+    epsr_urb = epsr_urb2d(i,j)
+    epsb_urb = epsb_urb2d(i,j)
+    epsg_urb = epsg_urb2d(i,j)
 
 ! Call urban
 
@@ -3493,10 +3528,14 @@ SUBROUTINE noahmp_urban(sf_urban_physics,   NSOIL,         IVGTYP,  ITIMESTEP,
             CMR_URB,    CHR_URB,    CMC_URB,    CHC_URB,                               &
             U10_URB,    V10_URB,    TH2_URB,     Q2_URB,                               & ! O
             UST_URB,     mh_urb,   stdh_urb,     lf_urb,   lp_urb,                     & ! 0
-            hgt_urb,    frc_urb,     lb_urb,      check, CMCR_URB,TGR_URB,             & ! H
+            hgt_urb,    frc_urb,     lb_urb,                                           & ! H
+           albr_urb,   albb_urb,   albg_urb,                                           & ! I
+           epsr_urb,   epsb_urb,   epsg_urb,                                           & ! I
+              check,   CMCR_URB,    TGR_URB,                                                   & ! H
            TGRL_URB,    SMR_URB,   CMGR_URB,   CHGR_URB,   jmonth,                     & ! H
           DRELR_URB,  DRELB_URB,                                                       & ! H
-          DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB )
+          DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB,                               &
+          lf_urb_s, z0_urb, vegfrac)
 
     TS_URB2D(I,J) = TS_URB
 
@@ -3505,7 +3544,7 @@ SUBROUTINE noahmp_urban(sf_urban_physics,   NSOIL,         IVGTYP,  ITIMESTEP,
     QFX(I,J)      = FRC_URB2D(I,J) * LH_KINEMATIC_URB &
                        + (1-FRC_URB2D(I,J))* QFX(I,J)                                  ![kg/m/m/s] 
     LH(I,J)       = FRC_URB2D(I,J) * LH_URB  + (1-FRC_URB2D(I,J)) * LH(I,J)            ![W/m/m]   
-    GRDFLX(I,J)   = FRC_URB2D(I,J) * (G_URB) + (1-FRC_URB2D(I,J)) * GRDFLX(I,J)        ![W/m/m]  
+    GRDFLX(I,J)   = FRC_URB2D(I,J) * (G_URB*(-1.0)) + (1-FRC_URB2D(I,J)) * GRDFLX(I,J) ![W/m/m] positive: downward direction
     TSK(I,J)      = FRC_URB2D(I,J) * TS_URB  + (1-FRC_URB2D(I,J)) * TSK(I,J)           ![K]    
 !    Q1            = QSFC(I,J)/(1.0+QSFC(I,J))                                         
 !    Q1            = FRC_URB2D(I,J) * QS_URB  + (1-FRC_URB2D(I,J)) * Q1                 ![-]
@@ -3734,7 +3773,7 @@ SUBROUTINE noahmp_urban(sf_urban_physics,   NSOIL,         IVGTYP,  ITIMESTEP,
 
 ! rename *_urb to sh_urb2d,lh_urb2d,g_urb2d,rn_urb2d
 
-    grdflx(i,j)   = (1.-frc_urb2d(i,j))*grdflx_rural(i,j)+ frc_urb2d(i,j)*grdflx_urb(i,j)
+    grdflx(i,j)   = (1.-frc_urb2d(i,j))*grdflx_rural(i,j)+ frc_urb2d(i,j)*grdflx_urb(i,j)*(-1.0) ! positive: downward direction
     qfx(i,j)      = (1.-frc_urb2d(i,j))*qfx_rural(i,j)   + qfx_urb(i,j)
     lh(i,j)       = qfx(i,j)*xlv
     hfx(i,j)      = hfx_urb(i,j)                         + (1-frc_urb2d(i,j))*hfx_rural(i,j)      ![W/m/m]

src/module_sf_noahmp_glacier.F

@@ -112,7 +112,7 @@ SUBROUTINE NOAHMP_GLACIER (&
                    FSA     ,FSR     ,FIRA    ,FSH     ,FGEV    ,SSOIL   , & ! OUT : 
                    TRAD    ,EDIR    ,RUNSRF  ,RUNSUB  ,SAG     ,ALBEDO  , & ! OUT :
                    QSNBOT  ,PONDING ,PONDING1,PONDING2,T2M     ,Q2E     , & ! OUT :
-                   EMISSI,  FPICE,    CH2B   ,QMELT                       & ! OUT :
+                   EMISSI,  FPICE,    CH2B   ,QMELT   ,HCPCT              & ! OUT :
 #ifdef WRF_HYDRO
                    , sfcheadrt                                            &
 #endif
@@ -191,6 +191,8 @@ SUBROUTINE NOAHMP_GLACIER (&
   REAL                           , INTENT(OUT)   :: EMISSI
   REAL                           , INTENT(OUT)   :: FPICE
   REAL                           , INTENT(OUT)   :: CH2B
+  REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT)   :: HCPCT  !heat capacity [j/m3/k]
+  REAL                           , INTENT(OUT)   :: QMELT  !internal pack melt due to phase change [mm/s]
 
 ! local
   INTEGER                                        :: IZ     !do-loop index
@@ -209,7 +211,6 @@ SUBROUTINE NOAHMP_GLACIER (&
   REAL                                           :: QDEW   !ground surface dew rate [mm/s]
   REAL                                           :: QVAP   !ground surface evap. rate [mm/s]
   REAL                                           :: LATHEA !latent heat [j/kg]
-  REAL, INTENT(OUT)                              :: QMELT  !internal pack melt due to phase change [mm/s]
   REAL                                           :: SWDOWN !downward solar [w/m2]
   REAL                                           :: BEG_WB !beginning water for error check
   REAL                                           :: ZBOT = -8.0 
@@ -245,7 +246,8 @@ SUBROUTINE NOAHMP_GLACIER (&
                          TAUSS  ,QSFC   ,                                 & !inout
                          IMELT  ,SNICEV ,SNLIQV ,EPORE  ,QMELT  ,PONDING, & !out
 		         SAG    ,FSA    ,FSR    ,FIRA   ,FSH    ,FGEV   , & !out
-		         TRAD   ,T2M    ,SSOIL  ,LATHEA ,Q2E    ,EMISSI, CH2B )   !out
+		         TRAD   ,T2M    ,SSOIL  ,LATHEA ,Q2E    ,EMISSI,  & !out
+                         CH2B   ,HCPCT  )   !out
 
     SICE = MAX(0.0, SMC - SH2O)   
     SNEQVO  = SNEQV
@@ -358,7 +360,7 @@ SUBROUTINE ENERGY_GLACIER (NSNOW  ,NSOIL  ,ISNOW  ,DT     ,QSNOW  ,RHOAIR , & !i
                              TAUSS  ,QSFC   ,                                 & !inout
                              IMELT  ,SNICEV ,SNLIQV ,EPORE  ,QMELT  ,PONDING, & !out
                              SAG    ,FSA    ,FSR    ,FIRA   ,FSH    ,FGEV   , & !out
-                             TRAD   ,T2M    ,SSOIL  ,LATHEA ,Q2E    ,EMISSI, CH2B )   !out
+                             TRAD   ,T2M    ,SSOIL  ,LATHEA ,Q2E,EMISSI, CH2B, HCPCT )   !out
 
 ! --------------------------------------------------------------------------------------------------
 ! --------------------------------------------------------------------------------------------------
@@ -426,6 +428,7 @@ SUBROUTINE ENERGY_GLACIER (NSNOW  ,NSOIL  ,ISNOW  ,DT     ,QSNOW  ,RHOAIR , & !i
   REAL                              , INTENT(OUT)   :: Q2E
   REAL                              , INTENT(OUT)   :: EMISSI
   REAL                              , INTENT(OUT)   :: CH2B   !sensible heat conductance, canopy air to ZLVL air (m/s)
+  REAL, DIMENSION(-NSNOW+1:NSOIL)   , INTENT(OUT)   :: HCPCT  !heat capacity [j/m3/k]
 
 
 ! local
@@ -438,7 +441,6 @@ SUBROUTINE ENERGY_GLACIER (NSNOW  ,NSOIL  ,ISNOW  ,DT     ,QSNOW  ,RHOAIR , & !i
   REAL                                              :: FIRE   !emitted IR (w/m2)
   REAL, DIMENSION(-NSNOW+1:NSOIL)                   :: FACT   !temporary used in phase change
   REAL, DIMENSION(-NSNOW+1:NSOIL)                   :: DF     !thermal conductivity [w/m/k]
-  REAL, DIMENSION(-NSNOW+1:NSOIL)                   :: HCPCT  !heat capacity [j/m3/k]
   REAL                                              :: GAMMA  !psychrometric constant (pa/k)
   REAL                                              :: RHSUR  !raltive humidity in surface soil/snow air space (-)
 
@@ -711,26 +713,25 @@ SUBROUTINE RADIATION_GLACIER (DT      ,TG      ,SNEQVO  ,SNEQV   ,COSZ    , & !i
 ! compute snow albedo only if cosz >0
 ! to be consistent with the main NoahMP code, currently do not include snow aging when sun is not present
 ! this needs more future work
-  if (COSZ > 0.0) then
 
-     ! snow age
+  ! snow age (allow nighttime aging)
 
-     CALL SNOW_AGE_GLACIER (DT,TG,SNEQVO,SNEQV,TAUSS,FAGE)
+  CALL SNOW_AGE_GLACIER (DT,TG,SNEQVO,SNEQV,TAUSS,FAGE)
 
+  if (cosz > 0.0) then ! only compute albedo when there is sunlight
      IF(OPT_ALB == 1) &
         CALL SNOWALB_BATS_GLACIER (NBAND,COSZ,FAGE,ALBSND,ALBSNI)
      IF(OPT_ALB == 2) THEN
         CALL SNOWALB_CLASS_GLACIER(NBAND,QSNOW,DT,ALB,ALBOLD,ALBSND,ALBSNI)
         ALBOLD = ALB
      END IF
-
   endif
 
 ! zero summed solar fluxes
 
-   SAG = 0.
-   FSA = 0.
-   FSR = 0.
+   SAG = 0.0
+   FSA = 0.0
+   FSR = 0.0
 
    FSNO = 0.0
    IF(SNEQV > 0.0) FSNO = 1.0