Merge branch 'CoLM-SYSU:master' into master

Makefile

@@ -77,6 +77,7 @@ OBJS_MKSRFDATA = \
 				  Aggregation_SoilParameters.o      \
 				  Aggregation_DBedrock.o            \
 				  Aggregation_Topography.o          \
+				  Aggregation_TopographyFactors.o   \
 				  Aggregation_Urban.o               \
 				  MOD_MeshFilter.o                  \
 				  MOD_RegionClip.o                  \
@@ -127,6 +128,7 @@ OBJS_BASIC =    \
 				 MOD_NdepData.o                 \
 				 MOD_FireData.o                 \
 				 MOD_OrbCoszen.o                \
+				 MOD_OrbCosazi.o                \
 				 MOD_3DCanopyRadiation.o        \
 				 MOD_Aerosol.o                  \
 				 MOD_SnowSnicar.o               \

main/MOD_OrbCosazi.F90

@@ -0,0 +1,69 @@
+#include <define.h>
+
+MODULE MOD_OrbCosazi
+
+!-----------------------------------------------------------------------
+   USE MOD_Precision
+   IMPLICIT NONE
+   SAVE
+
+! PUBLIC MEMBER FUNCTIONS:
+   PUBLIC :: orb_cosazi
+!-----------------------------------------------------------------------
+
+CONTAINS
+
+!-----------------------------------------------------------------------
+
+   FUNCTION orb_cosazi(calday, dlon, dlat, coszen)
+
+!-----------------------------------------------------------------------
+   USE MOD_Precision
+   IMPLICIT NONE
+
+   REAL(r8), intent(in) :: calday        !Julian cal day (1.xx to 365.xx)
+   REAL(r8), intent(in) :: dlat          !Centered latitude (radians)
+   REAL(r8), intent(in) :: dlon          !Centered longitude (radians)
+   REAL(r8), intent(in) :: coszen        !cosine of sun zenith angle
+   REAL(r8) :: orb_cosazi               !cosine of sun azimuth angle
+
+   ! --- Local variables ---
+   REAL(r8) declin                       !Solar declination (radians)
+   REAL(r8) eccf                         !Earth-sun distance factor (ie. (1/r)**2)
+   REAL(r8) lambm                        !Lambda m, mean long of perihelion (rad)
+   REAL(r8) lmm                          !Intermediate argument involving lambm
+   REAL(r8) lamb                         !Lambda, the earths long of perihelion
+   REAL(r8) invrho                       !Inverse normalized sun/earth distance
+   REAL(r8) sinl                         !Sine of lmm
+   REAL(r8) pi                           !3.14159265358979323846...
+   REAL(r8), parameter :: &
+          dayspy=365.0,              &!days per year
+          ve=80.5,                   &!Calday of vernal equinox assumes Jan 1 = calday 1
+          eccen=1.672393084E-2,      &!Eccentricity
+          obliqr=0.409214646,        &!Earths obliquity in radians
+          lambm0=-3.2625366E-2,      &!Mean long of perihelion at the vernal equinox (radians)
+          mvelpp=4.92251015           !moving vernal equinox longitude of
+                                         !perihelion plus pi (radians)
+   !-------------------------------------------------------------------------------
+
+      pi = 4.*atan(1.)
+      lambm = lambm0 + (calday - ve)*2.*pi/dayspy
+      lmm = lambm  - mvelpp
+
+      sinl = sin(lmm)
+      lamb = lambm + eccen*(2.*sinl + eccen*(1.25*sin(2.*lmm) &
+         + eccen*((13.0/12.0)*sin(3.*lmm) - 0.25*sinl)))
+      invrho = (1. + eccen*cos(lamb - mvelpp)) / (1. - eccen*eccen)
+
+      declin = asin(sin(obliqr)*sin(lamb))
+      eccf = invrho*invrho
+
+      orb_cosazi = (-1*cos(declin)*cos(calday*2.0*pi+dlon)- &
+         coszen*cos(dlat))/(sin(dlat)*sqrt(1-coszen*coszen))
+
+      IF (orb_cosazi<-1) orb_cosazi = -1
+      IF (orb_cosazi>1) orb_cosazi = 1
+
+   END FUNCTION orb_cosazi
+
+END MODULE MOD_OrbCosazi

main/MOD_OrbCoszen.F90

@@ -72,6 +72,9 @@ FUNCTION orb_coszen(calday,dlon,dlat)
       orb_coszen = sin(dlat)*sin(declin) &
                  - cos(dlat)*cos(declin)*cos(calday*2.0*pi+dlon)
 
+      IF (orb_coszen<0) orb_coszen = 0
+      IF (orb_coszen>1) orb_coszen = 1
+
    END FUNCTION orb_coszen
 
 END MODULE MOD_OrbCoszen

main/MOD_Vars_1DForcing.F90

@@ -30,6 +30,7 @@ MODULE MOD_Vars_1DForcing
    real(r8), allocatable :: forc_solsd (:) ! atm vis diffuse solar rad onto srf [W/m2]
    real(r8), allocatable :: forc_solld (:) ! atm nir diffuse solar rad onto srf [W/m2]
    real(r8), allocatable :: forc_frl   (:) ! atmospheric infrared (longwave) radiation [W/m2]
+   real(r8), allocatable :: forc_swrad (:) ! atmospheric shortwave radiation [W/m2]
    real(r8), allocatable :: forc_hgt_u (:) ! observational height of wind [m]
    real(r8), allocatable :: forc_hgt_t (:) ! observational height of temperature [m]
    real(r8), allocatable :: forc_hgt_q (:) ! observational height of humidity [m]
@@ -82,6 +83,7 @@ SUBROUTINE allocate_1D_Forcing
             allocate (forc_solsd  (numpatch) ) ! atm vis diffuse solar rad onto srf [W/m2]
             allocate (forc_solld  (numpatch) ) ! atm nir diffuse solar rad onto srf [W/m2]
             allocate (forc_frl    (numpatch) ) ! atmospheric infrared (longwave) radiation [W/m2]
+            allocate (forc_swrad  (numpatch) ) ! atmospheric shortwave radiation [W/m2]
             allocate (forc_hgt_u  (numpatch) ) ! observational height of wind [m]
             allocate (forc_hgt_t  (numpatch) ) ! observational height of temperature [m]
             allocate (forc_hgt_q  (numpatch) ) ! observational height of humidity [m]
@@ -133,6 +135,7 @@ SUBROUTINE deallocate_1D_Forcing ()
             deallocate ( forc_solsd  ) ! atm vis diffuse solar rad onto srf [W/m2]
             deallocate ( forc_solld  ) ! atm nir diffuse solar rad onto srf [W/m2]
             deallocate ( forc_frl    ) ! atmospheric infrared (longwave) radiation [W/m2]
+            deallocate ( forc_swrad  ) ! atmospheric shortwave radiation [W/m2]
             deallocate ( forc_hgt_u  ) ! observational height of wind [m]
             deallocate ( forc_hgt_t  ) ! observational height of temperature [m]
             deallocate ( forc_hgt_q  ) ! observational height of humidity [m]

main/MOD_Vars_Global.F90

@@ -54,6 +54,11 @@ MODULE MOD_Vars_Global
    integer, parameter :: nl_roof   = 10
    integer, parameter :: nl_wall   = 10
    integer, parameter :: nvegwcs   = 4  ! number of vegetation water potential nodes
+
+   ! used for downscaling
+   integer, parameter :: num_type = 4
+   integer, parameter :: num_zenith = 51
+   integer, parameter :: num_azimuth = 36
 
    ! bgc variables
    integer, parameter :: ndecomp_pools        = 7

main/MOD_Vars_TimeInvariants.F90

@@ -255,6 +255,14 @@ MODULE MOD_Vars_TimeInvariants
    real(r8) :: tcrit                            !critical temp. to determine rain or snow
    real(r8) :: wetwatmax                        !maximum wetland water (mm)
 
+   ! Used for downscaling
+   real(r8), allocatable    :: svf_patches (:)                                           ! sky view factor
+   real(r8), allocatable    :: cur_patches (:)                                           ! curvature
+   real(r8), allocatable    :: sf_lut_patches (:,:,:)                                    ! look up table of shadow factor of a patch
+   real(r8), allocatable    :: asp_type_patches        (:,:)                             ! topographic aspect of each character of one patch
+   real(r8), allocatable    :: slp_type_patches        (:,:)                             ! topographic slope of each character of one patch
+   real(r8), allocatable    :: area_type_patches       (:,:)                             ! area percentage of each character of one patch
+
 ! PUBLIC MEMBER FUNCTIONS:
    PUBLIC :: allocate_TimeInvariants
    PUBLIC :: deallocate_TimeInvariants
@@ -342,6 +350,14 @@ SUBROUTINE allocate_TimeInvariants ()
             allocate (ibedrock             (numpatch))
             allocate (topoelv              (numpatch))
             allocate (topostd              (numpatch))
+
+            ! Used for downscaling
+            allocate (svf_patches          (numpatch))
+            allocate (asp_type_patches     (num_type,numpatch))
+            allocate (slp_type_patches     (num_type,numpatch))
+            allocate (area_type_patches    (num_type,numpatch))
+            allocate (sf_lut_patches       (num_azimuth,num_zenith,numpatch))
+            allocate (cur_patches          (numpatch))
       ENDIF
 
 #if (defined LULC_IGBP_PFT || defined LULC_IGBP_PC)
@@ -384,7 +400,7 @@ SUBROUTINE READ_TimeInvariants (lc_year, casename, dir_restart)
    character(len=*), intent(in) :: dir_restart
 
    ! Local variables
-   character(len=256) :: file_restart, cyear
+   character(len=256) :: file_restart, cyear, lndname
 
       write(cyear,'(i4.4)') lc_year
       file_restart = trim(dir_restart) // '/const/' // trim(casename) //'_restart_const' // '_lc' // trim(cyear) // '.nc'
@@ -468,6 +484,15 @@ SUBROUTINE READ_TimeInvariants (lc_year, casename, dir_restart)
       CALL ncio_read_bcast_serial (file_restart, 'tcrit ', tcrit ) ! critical temp. to determine rain or snow
       CALL ncio_read_bcast_serial (file_restart, 'wetwatmax', wetwatmax) ! maximum wetland water (mm)
 
+      IF (DEF_USE_Forcing_Downscaling) THEN
+         CALL ncio_read_vector (file_restart, 'slp_type_patches', num_type, landpatch, slp_type_patches)
+         CALL ncio_read_vector (file_restart, 'svf_patches', landpatch, svf_patches)
+         CALL ncio_read_vector (file_restart, 'asp_type_patches', num_type, landpatch, asp_type_patches)
+         CALL ncio_read_vector (file_restart, 'area_type_patches', num_type, landpatch, area_type_patches)
+         CALL ncio_read_vector (file_restart, 'sf_lut_patches', num_azimuth, num_zenith, landpatch, sf_lut_patches)
+         CALL ncio_read_vector (file_restart, 'cur_patches', landpatch, cur_patches)
+       ENDIF
+
 #if (defined LULC_IGBP_PFT || defined LULC_IGBP_PC)
       file_restart = trim(dir_restart) // '/const/' // trim(casename) //'_restart_pft_const' // '_lc' // trim(cyear) // '.nc'
       CALL READ_PFTimeInvariants (file_restart)
@@ -546,6 +571,9 @@ SUBROUTINE WRITE_TimeInvariants (lc_year, casename, dir_restart)
       CALL ncio_define_dimension_vector (file_restart, landpatch, 'soilsnow', nl_soil-maxsnl)
       CALL ncio_define_dimension_vector (file_restart, landpatch, 'soil',     nl_soil)
       CALL ncio_define_dimension_vector (file_restart, landpatch, 'lake',     nl_lake)
+      CALL ncio_define_dimension_vector (file_restart, landpatch, 'type',     num_type)
+      CALL ncio_define_dimension_vector (file_restart, landpatch, 'azi',      num_azimuth)
+      CALL ncio_define_dimension_vector (file_restart, landpatch, 'zen',      num_zenith)
 
       CALL ncio_write_vector (file_restart, 'patchclass', 'patch', landpatch, patchclass)                            !
       CALL ncio_write_vector (file_restart, 'patchtype' , 'patch', landpatch, patchtype )                            !
@@ -610,6 +638,15 @@ SUBROUTINE WRITE_TimeInvariants (lc_year, casename, dir_restart)
 
       CALL ncio_write_vector (file_restart, 'topoelv', 'patch', landpatch, topoelv)
       CALL ncio_write_vector (file_restart, 'topostd', 'patch', landpatch, topostd)
+
+      IF (DEF_USE_Forcing_Downscaling) THEN
+         CALL ncio_write_vector (file_restart, 'svf_patches', 'patch', landpatch, svf_patches)
+         CALL ncio_write_vector (file_restart, 'cur_patches', 'patch', landpatch, cur_patches)
+         CALL ncio_write_vector (file_restart, 'slp_type_patches', 'type', num_type, 'patch', landpatch, slp_type_patches)
+         CALL ncio_write_vector (file_restart, 'asp_type_patches', 'type', num_type, 'patch', landpatch, asp_type_patches)
+         CALL ncio_write_vector (file_restart, 'area_type_patches', 'type', num_type, 'patch', landpatch, area_type_patches)
+         CALL ncio_write_vector (file_restart, 'sf_lut_patches', 'azi', num_azimuth, 'zen', num_zenith, 'patch', landpatch, sf_lut_patches)
+      ENDIF
 
 #ifdef USEMPI
       CALL mpi_barrier (p_comm_glb, p_err)
@@ -662,6 +699,7 @@ END SUBROUTINE WRITE_TimeInvariants
 
    SUBROUTINE deallocate_TimeInvariants ()
 
+   USE MOD_Namelist, only: DEF_USE_Forcing_Downscaling 
    USE MOD_SPMD_Task
    USE MOD_LandPatch, only: numpatch
 
@@ -736,6 +774,15 @@ SUBROUTINE deallocate_TimeInvariants ()
             deallocate (topoelv        )
             deallocate (topostd        )
 
+            IF (DEF_USE_Forcing_Downscaling) THEN
+               deallocate(slp_type_patches  )
+               deallocate(svf_patches       )
+               deallocate(asp_type_patches  )
+               deallocate(area_type_patches )
+               deallocate(sf_lut_patches    )
+               deallocate(cur_patches       )
+            ENDIF
+
          ENDIF
       ENDIF
 
@@ -758,7 +805,7 @@ SUBROUTINE check_TimeInvariants ()
 
    USE MOD_SPMD_Task
    USE MOD_RangeCheck
-   USE MOD_Namelist, only : DEF_USE_BEDROCK
+   USE MOD_Namelist, only : DEF_USE_BEDROCK, DEF_USE_Forcing_Downscaling
 
    IMPLICIT NONE
 
@@ -817,6 +864,15 @@ SUBROUTINE check_TimeInvariants ()
       CALL check_vector_data ('topostd      [m]     ', topostd     ) !
       CALL check_vector_data ('BVIC        [-]      ', BVIC        ) !
 
+      IF (DEF_USE_Forcing_Downscaling) THEN
+         CALL check_vector_data ('slp_type_patches     [rad] ' , slp_type_patches)      ! slope
+         CALL check_vector_data ('svf_patches          [-] '   , svf_patches)           ! sky view factor
+         CALL check_vector_data ('asp_type_patches     [rad] ' , asp_type_patches)      ! aspect
+         CALL check_vector_data ('area_type_patches    [-] '   , area_type_patches)     ! area percent
+         CALL check_vector_data ('cur_patches          [-]'    , cur_patches )
+         CALL check_vector_data ('sf_lut_patches       [-] '   , sf_lut_patches)        ! shadow mask
+      ENDIF
+
 #ifdef USEMPI
       CALL mpi_barrier (p_comm_glb, p_err)
 #endif

mkinidata/MOD_Initialize.F90

@@ -103,7 +103,7 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, &
    character(len=256) :: fsoildat
    character(len=256) :: fsnowdat
    character(len=256) :: fcndat
-   character(len=256) :: ftopo
+   character(len=256) :: ftopo, lndname
    type(grid_type) :: gsoil
    type(grid_type) :: gsnow
    type(grid_type) :: gcn
@@ -383,7 +383,32 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, &
       ftopo = trim(dir_landdata)//'/topography/'//trim(cyear)//'/topostd_patches.nc'
       CALL ncio_read_vector (ftopo, 'topostd_patches', landpatch, topostd)
 #endif
-
+! ......................................
+! 1.5 Initialize topography factor data
+! ......................................
+#ifdef SinglePoint
+      slp_type_patches(:,1) = SITE_slp_type
+      asp_type_patches(:,1) = SITE_asp_type
+      area_type_patches(:,1) = SITE_area_type
+      svf_patches(:) = SITE_svf
+      cur_patches(:) = SITE_cur
+      sf_lut_patches(:,:,1) = SITE_sf_lut
+#else
+      IF (DEF_USE_Forcing_Downscaling) THEN
+         lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/slp_type_patches.nc'             ! slope
+         CALL ncio_read_vector (lndname, 'slp_type_patches', num_type, landpatch, slp_type_patches)
+         lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/svf_patches.nc'               ! sky view factor
+         CALL ncio_read_vector (lndname, 'svf_patches', landpatch, svf_patches)
+         lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/asp_type_patches.nc'            ! aspect
+         CALL ncio_read_vector (lndname, 'asp_type_patches', num_type, landpatch, asp_type_patches)
+         lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/area_type_patches.nc'         ! area percent
+         CALL ncio_read_vector (lndname, 'area_type_patches', num_type, landpatch, area_type_patches)
+         lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/sf_lut_patches.nc'        ! shadow mask
+         CALL ncio_read_vector (lndname, 'sf_lut_patches', num_azimuth, num_zenith, landpatch, sf_lut_patches)
+         lndname = trim(DEF_dir_landdata) // '/topography/'//trim(cyear)//'/cur_patches.nc'               ! curvature
+         CALL ncio_read_vector (lndname, 'cur_patches', landpatch, cur_patches)
+       ENDIF
+#endif
 ! ................................
 ! 1.6 Initialize TUNABLE constants
 ! ................................