Merge branch 'CoLM-SYSU:master' into master

main/MOD_Forcing.F90

@@ -487,7 +487,7 @@ SUBROUTINE read_forcing (idate, dir_forcing)
                         cosz = orb_coszen(calday, gforc%rlon(ilon), gforc%rlat(ilat))
                         cosz = max(0.001, cosz)
                         ! 10/24/2024, yuan: deal with time log with backward or foreward
-                        IF (timelog(ivar) == 'foreward') THEN
+                        IF (trim(timelog(ivar)) == 'foreward') THEN
                            forcn(ivar)%blk(ib,jb)%val(i,j) = &
                               cosz / avgcos%blk(ib,jb)%val(i,j) * forcn_LB(ivar)%blk(ib,jb)%val(i,j)
                         ELSE

main/MOD_GroundFluxes.F90

@@ -32,18 +32,20 @@ SUBROUTINE GroundFluxes (zlnd, zsno, hu, ht, hq, hpbl, &
 ! Original author : Yongjiu Dai, 09/15/1999; 08/30/2002
 !
 ! REVISIONS:
-! Hua Yuan, 09/2019: removed sigf to be consistant with PFT runs, removed fsena,
-!                    fevpa, renamed z0ma to z0m
-! Shaofeng Liu, 05/2023: add option to call moninobuk_leddy, the LargeEddy
-!                        surface turbulence scheme (LZD2022);
-!                        make a proper update of um.
+! 09/2019, Hua Yuan: removed sigf to be consistant with PFT runs, removed
+!          fsena, fevpa, renamed z0ma to z0m.
+!
+! 05/2023, Shaofeng Liu: add option to call moninobuk_leddy, the LargeEddy
+!          surface turbulence scheme (LZD2022); make a proper update of um.
+!
 !=======================================================================
 
    USE MOD_Precision
    USE MOD_Const_Physical, only: cpair,vonkar,grav
    USE MOD_FrictionVelocity
    USE mod_namelist, only: DEF_USE_CBL_HEIGHT,DEF_RSS_SCHEME
    USE MOD_TurbulenceLEddy
+
    IMPLICIT NONE
 
 !----------------------- Dummy argument --------------------------------

main/MOD_LAIReadin.F90

@@ -69,7 +69,7 @@ SUBROUTINE LAI_readin (year, time, dir_landdata)
 
 #ifdef SinglePoint
 #ifndef URBAN_MODEL
-      iyear = findloc_ud(SITE_LAI_year == min(DEF_LAI_END_YEAR, max(DEF_LAI_START_YEAR,year) )
+      iyear = findloc_ud(SITE_LAI_year == min(DEF_LAI_END_YEAR, max(DEF_LAI_START_YEAR,year)))
       IF (.not. DEF_LAI_MONTHLY) THEN
          itime = (time-1)/8 + 1
       ENDIF

main/MOD_Lake.F90

@@ -1839,16 +1839,16 @@ SUBROUTINE snowwater_lake ( USE_Dynamic_Lake, &
          a = wliq_soisno(j)/(dz_soisno(j)*denh2o) + wice_soisno(j)/(dz_soisno(j)*denice)
 
          IF (a < porsl(j)) THEN
-            wliq_soisno(j) = max( 0., (porsl(j)*dz_soisno(j) - wice_soisno(j)/denice)*denh2o )
-            wice_soisno(j) = max( 0., (porsl(j)*dz_soisno(j) - wliq_soisno(j)/denh2o)*denice )
+            wliq_soisno(j) = max(0., (porsl(j)*dz_soisno(j) - wice_soisno(j)/denice)*denh2o )
+            wice_soisno(j) = max(0., (porsl(j)*dz_soisno(j) - wliq_soisno(j)/denh2o)*denice )
          ELSE
             wliq_soisno(j) = max(0., wliq_soisno(j) - (a - porsl(j))*denh2o*dz_soisno(j) )
-            wice_soisno(j) = max( 0., (porsl(j)*dz_soisno(j) - wliq_soisno(j)/denh2o)*denice )
+            wice_soisno(j) = max(0., (porsl(j)*dz_soisno(j) - wliq_soisno(j)/denh2o)*denice )
          ENDIF
 
          IF (wliq_soisno(j) > porsl(j)*denh2o*dz_soisno(j)) THEN
-             wliq_soisno(j) = porsl(j)*denh2o*dz_soisno(j)
-             wice_soisno(j) = 0.0
+            wliq_soisno(j) = porsl(j)*denh2o*dz_soisno(j)
+            wice_soisno(j) = 0.0
          ENDIF
 
          dw_soil = dw_soil - wliq_soisno(j) - wice_soisno(j)
@@ -1878,7 +1878,16 @@ SUBROUTINE snowwater_lake ( USE_Dynamic_Lake, &
          lake_icefrac(1) = min(max(lake_icefrac(1), 0.), 1.)
 
          dz_lake(nl_lake) = dz_lake(nl_lake) + dw_soil/1.e3
-         dz_lake(nl_lake) = max(dz_lake(nl_lake), 0.)
+
+         IF (dz_lake(nl_lake) < 0.) THEN
+            j = nl_lake
+            DO WHILE (dz_lake(j) < 0.)
+               IF (j > 1) dz_lake(j-1) = dz_lake(j-1) + dz_lake(j)
+               dz_lake(j) = 0.
+               j = j - 1
+               IF (j == 0) EXIT
+            ENDDO 
+         ENDIF
 
          CALL adjust_lake_layer (nl_lake, dz_lake, t_lake, lake_icefrac)
 

main/MOD_LeafTemperature.F90

@@ -105,6 +105,7 @@ SUBROUTINE LeafTemperature ( &
    USE MOD_CanopyLayerProfile
    USE MOD_TurbulenceLEddy
    USE MOD_AssimStomataConductance
+   USE MOD_UserSpecifiedForcing, only: HEIGHT_mode
    USE MOD_Vars_TimeInvariants, only: patchclass
    USE MOD_Const_LC, only: z0mr, displar
    USE MOD_PlantHydraulic, only :PlantHydraulicStress_twoleaf, getvegwp_twoleaf
@@ -304,6 +305,9 @@ SUBROUTINE LeafTemperature ( &
 
    real(r8) :: &
         displa,     &! displacement height [m]
+        hu_,        &! adjusted observational height of wind [m]
+        ht_,        &! adjusted observational height of temperature [m]
+        hq_,        &! adjusted observational height of humidity [m]
         zldis,      &! reference height "minus" zero displacement heght [m]
         zii,        &! convective boundary layer height [m]
         z0mv,       &! roughness length, momentum [m]
@@ -518,7 +522,23 @@ SUBROUTINE LeafTemperature ( &
       dth = thm - taf
       dqh = qm  - qaf
       dthv  = dth*(1.+0.61*qm) + 0.61*th*dqh
-      zldis = hu - displa
+
+      hu_ = hu; ht_ = ht; hq_ = hq;
+
+      IF (trim(HEIGHT_mode) == 'absolute') THEN
+#ifndef SingPoint
+         ! to ensure the obs height >= htop+10.
+         hu_ = max(htop+10., hu)
+         ht_ = max(htop+10., ht)
+         hq_ = max(htop+10., hq)
+#endif
+      ELSE ! relative height
+         hu_ = htop + hu
+         ht_ = htop + ht
+         hq_ = htop + hq
+      ENDIF
+
+      zldis = hu_ - displa
 
       IF(zldis <= 0.0) THEN
          write(6,*) 'the obs height of u less than the zero displacement heght'
@@ -550,11 +570,11 @@ SUBROUTINE LeafTemperature ( &
 ! Evaluate stability-dependent variables using moz from prior iteration
          IF (rd_opt == 3) THEN
             IF (DEF_USE_CBL_HEIGHT) THEN
-              CALL moninobukm_leddy(hu,ht,hq,displa,z0mv,z0hv,z0qv,obu,um, &
+              CALL moninobukm_leddy(hu_,ht_,hq_,displa,z0mv,z0hv,z0qv,obu,um, &
                                     displasink,z0mv,hpbl,ustar,fh2m,fq2m, &
                                     htop,fmtop,fm,fh,fq,fht,fqt,phih)
             ELSE
-              CALL moninobukm(hu,ht,hq,displa,z0mv,z0hv,z0qv,obu,um, &
+              CALL moninobukm(hu_,ht_,hq_,displa,z0mv,z0hv,z0qv,obu,um, &
                               displasink,z0mv,ustar,fh2m,fq2m, &
                               htop,fmtop,fm,fh,fq,fht,fqt,phih)
             ENDIF
@@ -564,10 +584,10 @@ SUBROUTINE LeafTemperature ( &
             raw = 1./(vonkar/(fq-fqt)*ustar)
          ELSE
             IF (DEF_USE_CBL_HEIGHT) THEN
-               CALL moninobuk_leddy(hu,ht,hq,displa,z0mv,z0hv,z0qv,obu,um,hpbl, &
+               CALL moninobuk_leddy(hu_,ht_,hq_,displa,z0mv,z0hv,z0qv,obu,um,hpbl, &
                                     ustar,fh2m,fq2m,fm10m,fm,fh,fq)
             ELSE
-               CALL moninobuk(hu,ht,hq,displa,z0mv,z0hv,z0qv,obu,um,&
+               CALL moninobuk(hu_,ht_,hq_,displa,z0mv,z0hv,z0qv,obu,um,&
                               ustar,fh2m,fq2m,fm10m,fm,fh,fq)
             ENDIF
             ! Aerodynamic resistance
@@ -918,7 +938,7 @@ SUBROUTINE LeafTemperature ( &
            um = max(ur,.1)
          ELSE
            IF (DEF_USE_CBL_HEIGHT) THEN !//TODO: Shaofeng, 2023.05.18
-             zii = max(5.*hu,hpbl)
+             zii = max(5.*hu_,hpbl)
            ENDIF !//TODO: Shaofeng, 2023.05.18
            wc = (-grav*ustar*thvstar*zii/thv)**(1./3.)
           wc2 = beta*beta*(wc*wc)

main/MOD_LeafTemperaturePC.F90

@@ -91,14 +91,14 @@ SUBROUTINE  LeafTemperaturePC ( &
 !
 ! !REVISIONS:
 !
-!  01/2021, Xingjie Lu and Nan Wei: added plant hydraulic process interface
+!  01/2021, Xingjie Lu and Nan Wei: added plant hydraulic process interface.
 !
-!  01/2021, Nan Wei: added interaction btw prec and canopy
+!  01/2021, Nan Wei: added interaction btw prec and canopy.
 !
 !  05/2023, Shaofeng Liu: add option to call moninobuk_leddy, the LargeEddy
 !           surface turbulence scheme (LZD2022); make a proper update of um.
 !
-!  04/2024, Hua Yuan: add option to account for vegetation snow process
+!  04/2024, Hua Yuan: add option to account for vegetation snow process.
 !
 !=======================================================================
 
@@ -114,6 +114,8 @@ SUBROUTINE  LeafTemperaturePC ( &
    USE MOD_AssimStomataConductance
    USE MOD_PlantHydraulic, only: PlantHydraulicStress_twoleaf
    USE MOD_Ozone, only: CalcOzoneStress
+   USE MOD_UserSpecifiedForcing, only: HEIGHT_mode
+
    IMPLICIT NONE
 
 !-----------------------Arguments---------------------------------------
@@ -325,6 +327,9 @@ SUBROUTINE  LeafTemperaturePC ( &
         rootfr(nl_soil,ps:pe) ! root fraction
 
    real(r8) :: &
+        hu_,           &! adjusted observational height of wind [m]
+        ht_,           &! adjusted observational height of temperature [m]
+        hq_,           &! adjusted observational height of humidity [m]
         zldis,         &! reference height "minus" zero displacement heght [m]
         zii,           &! convective boundary layer height [m]
         z0mv,          &! roughness length, momentum [m]
@@ -901,7 +906,23 @@ SUBROUTINE  LeafTemperaturePC ( &
       dth   = thm - taf(toplay)
       dqh   = qm  - qaf(toplay)
       dthv  = dth*(1.+0.61*qm) + 0.61*th*dqh
-      zldis = hu - displa_lays(3)
+
+      hu_ = hu; ht_ = ht; hq_ = hq;
+
+      IF (trim(HEIGHT_mode) == 'absolute') THEN
+#ifndef SingPoint
+         ! to ensure the obs height >= htop+10.
+         hu_ = max(htop_lay(toplay)+10., hu)
+         ht_ = max(htop_lay(toplay)+10., ht)
+         hq_ = max(htop_lay(toplay)+10., hq)
+#endif
+      ELSE ! relative height
+         hu_ = htop_lay(toplay) + hu
+         ht_ = htop_lay(toplay) + ht
+         hq_ = htop_lay(toplay) + hq
+      ENDIF
+
+      zldis = hu_ - displa_lays(3)
 
       IF(zldis <= 0.0) THEN
          write(6,*) 'the obs height of u less than the zero displacement heght'
@@ -935,11 +956,11 @@ SUBROUTINE  LeafTemperaturePC ( &
 ! Evaluate stability-dependent variables using moz from prior iteration
 
          IF (DEF_USE_CBL_HEIGHT) THEN
-            CALL moninobukm_leddy(hu,ht,hq,displa_lays(toplay),z0mv,z0hv,z0qv,obu,um, &
+            CALL moninobukm_leddy(hu_,ht_,hq_,displa_lays(toplay),z0mv,z0hv,z0qv,obu,um, &
                                   displa_lay(toplay),z0m_lay(toplay),hpbl,ustar,fh2m,fq2m, &
                                   htop_lay(toplay),fmtop,fm,fh,fq,fht,fqt,phih)
          ELSE
-            CALL moninobukm(hu,ht,hq,displa_lays(toplay),z0mv,z0hv,z0qv,obu,um, &
+            CALL moninobukm(hu_,ht_,hq_,displa_lays(toplay),z0mv,z0hv,z0qv,obu,um, &
                             displa_lay(toplay),z0m_lay(toplay),ustar,fh2m,fq2m, &
                             htop_lay(toplay),fmtop,fm,fh,fq,fht,fqt,phih)
          ENDIF
@@ -1659,7 +1680,7 @@ SUBROUTINE  LeafTemperaturePC ( &
             um = max(ur,.1)
          ELSE
             IF (DEF_USE_CBL_HEIGHT) THEN !//TODO: Shaofeng, 2023.05.18
-              zii = max(5.*hu,hpbl)
+              zii = max(5.*hu_,hpbl)
             ENDIF !//TODO: Shaofeng, 2023.05.18
             wc = (-grav*ustar*thvstar*zii/thv)**(1./3.)
             wc2 = beta*beta*(wc*wc)

main/MOD_UserSpecifiedForcing.F90

@@ -45,6 +45,8 @@ MODULE MOD_UserSpecifiedForcing
    character(len=256) :: dataset
 
    logical  :: solarin_all_band   ! whether solar radiation in all bands is available
+
+   character(len=256) :: HEIGHT_mode ! observation height mode
    real(r8) :: HEIGHT_V           ! observation height of wind speed
    real(r8) :: HEIGHT_T           ! observation height of air temperature
    real(r8) :: HEIGHT_Q           ! observation height of specific humidity
@@ -110,6 +112,7 @@ SUBROUTINE init_user_specified_forcing
       allocate (tintalgo (NVAR))
 
       solarin_all_band = DEF_forcing%solarin_all_band ! whether solar radiation in all bands is available
+      HEIGHT_mode      = DEF_forcing%HEIGHT_mode      ! observation height mode
       HEIGHT_V         = DEF_forcing%HEIGHT_V         ! observation height of wind speed
       HEIGHT_T         = DEF_forcing%HEIGHT_T         ! observation height of air temperature
       HEIGHT_Q         = DEF_forcing%HEIGHT_Q         ! observation height of specific humidity

main/MOD_Vars_TimeVariables.F90

@@ -884,7 +884,7 @@ FUNCTION save_to_restart (idate, deltim, itstamp, ptstamp) result(rwrite)
       ENDSELECT
 
       IF (rwrite) THEN
-         rwrite = (ptstamp <= itstamp)
+         rwrite = ((ptstamp <= itstamp) .or. isendofyear(idate,deltim))
       ENDIF
 
    END FUNCTION save_to_restart

main/URBAN/MOD_Urban_Flux.F90

@@ -111,6 +111,7 @@ SUBROUTINE UrbanOnlyFlux ( &
    USE MOD_Const_Physical, only: cpair,vonkar,grav,hvap
    USE MOD_FrictionVelocity
    USE MOD_CanopyLayerProfile
+   USE MOD_UserSpecifiedForcing, only: HEIGHT_mode
    IMPLICIT NONE
 
 !----------------------- Dummy argument --------------------------------
@@ -262,9 +263,9 @@ SUBROUTINE UrbanOnlyFlux ( &
         numlay         ! available layer number
 
    real(r8) :: &
-        huu,          &! observational height of wind [m]
-        htu,          &! observational height of temperature [m]
-        hqu,          &! observational height of humidity [m]
+        hu_,          &! adjusted observational height of wind [m]
+        ht_,          &! adjusted observational height of temperature [m]
+        hq_,          &! adjusted observational height of humidity [m]
         ktop,         &! K value at a specific height
         utop,         &! u value at a specific height
         fht,          &! integral of profile function for heat at the top layer
@@ -509,12 +510,22 @@ SUBROUTINE UrbanOnlyFlux ( &
       dqh  =  qm - qaf(2)
       dthv = dth*(1.+0.61*qm) + 0.61*th*dqh
 
-      ! to ensure the obs height >= hroof+10.
-      huu = max(hroof+10., hu)
-      htu = max(hroof+10., ht)
-      hqu = max(hroof+10., hq)
+      hu_ = hu; ht_ = ht; hq_ = hq;
 
-      zldis = huu - displa
+      IF (trim(HEIGHT_mode) == 'absolute') THEN
+#ifndef SingPoint
+         ! to ensure the obs height >= hroof+10.
+         hu_ = max(hroof+10., hu)
+         ht_ = max(hroof+10., ht)
+         hq_ = max(hroof+10., hq)
+#endif
+      ELSE ! relative height
+         hu_ = hroof + hu
+         ht_ = hroof + ht
+         hq_ = hroof + hq
+      ENDIF
+
+      zldis = hu_ - displa
 
       IF (zldis <= 0.0) THEN
          write(6,*) 'the obs height of u less than the zero displacement heght'
@@ -538,7 +549,7 @@ SUBROUTINE UrbanOnlyFlux ( &
 
          !NOTE: displat=hroof, z0mt=0, are set for roof
          ! fmtop is calculated at the same height of fht, fqt
-         CALL moninobukm(huu,htu,hqu,displa,z0m,z0h,z0q,obu,um, &
+         CALL moninobukm(hu_,ht_,hq_,displa,z0m,z0h,z0q,obu,um, &
               hroof,0.,ustar,fh2m,fq2m,hroof,fmtop,fm,fh,fq,fht,fqt,phih)
 
 ! Aerodynamic resistance
@@ -885,6 +896,7 @@ SUBROUTINE  UrbanVegFlux ( &
    USE MOD_FrictionVelocity
    USE MOD_CanopyLayerProfile
    USE MOD_AssimStomataConductance
+   USE MOD_UserSpecifiedForcing, only: HEIGHT_mode
    IMPLICIT NONE
 
 !-----------------------Arguments---------------------------------------
@@ -1150,9 +1162,9 @@ SUBROUTINE  UrbanVegFlux ( &
         numlay         ! available layer number
 
    real(r8) :: &
-        huu,          &! observational height of wind [m]
-        htu,          &! observational height of temperature [m]
-        hqu,          &! observational height of humidity [m]
+        hu_,          &! adjusted observational height of wind [m]
+        ht_,          &! adjusted observational height of temperature [m]
+        hq_,          &! adjusted observational height of humidity [m]
         ktop,         &! K value at a specific height
         utop,         &! u value at a specific height
         fht,          &! integral of profile function for heat at the top layer
@@ -1487,12 +1499,22 @@ SUBROUTINE  UrbanVegFlux ( &
       dqh =  qm - qaf(2)
       dthv = dth*(1.+0.61*qm) + 0.61*th*dqh
 
-      ! To ensure the obs height >= hroof+10.
-      huu = max(hroof+10., hu)
-      htu = max(hroof+10., ht)
-      hqu = max(hroof+10., hq)
+      hu_ = hu; ht_ = ht; hq_ = hq;
+
+      IF (trim(HEIGHT_mode) == 'absolute') THEN
+#ifndef SingPoint
+         ! to ensure the obs height >= hroof+10.
+         hu_ = max(hroof+10., hu)
+         ht_ = max(hroof+10., ht)
+         hq_ = max(hroof+10., hq)
+#endif
+      ELSE ! relative height
+         hu_ = hroof + hu
+         ht_ = hroof + ht
+         hq_ = hroof + hq
+      ENDIF
 
-      zldis = huu - displa
+      zldis = hu_ - displa
 
       IF (zldis <= 0.0) THEN
          write(6,*) 'the obs height of u less than the zero displacement heght'
@@ -1517,7 +1539,7 @@ SUBROUTINE  UrbanVegFlux ( &
 !-----------------------------------------------------------------------
 ! Evaluate stability-dependent variables using moz from prior iteration
 
-         CALL moninobukm(huu,htu,hqu,displa,z0m,z0h,z0q,obu,um, &
+         CALL moninobukm(hu_,ht_,hq_,displa,z0m,z0h,z0q,obu,um, &
               hroof,0.,ustar,fh2m,fq2m,hroof,fmtop,fm,fh,fq,fht,fqt,phih)
 
 ! Aerodynamic resistance