Add forcing downscaling module

The same with #272 PR but modify on the newest edition of CoLM

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,68 @@
+#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

@@ -15,6 +15,7 @@ MODULE MOD_Vars_PFTimeInvariants
 
    USE MOD_Precision
    USE MOD_Vars_Global
+
    IMPLICIT NONE
    SAVE
 
@@ -23,9 +24,6 @@ MODULE MOD_Vars_PFTimeInvariants
    real(r8), allocatable :: pftfrac     (:)    !PFT fractional cover
    real(r8), allocatable :: htop_p      (:)    !canopy top height [m]
    real(r8), allocatable :: hbot_p      (:)    !canopy bottom height [m]
-#ifdef CROP
-   real(r8), allocatable :: cropfrac    (:)    !Crop fractional cover
-#endif
 
 ! PUBLIC MEMBER FUNCTIONS:
    PUBLIC :: allocate_PFTimeInvariants
@@ -50,7 +48,6 @@ SUBROUTINE allocate_PFTimeInvariants
    ! --------------------------------------------------------------------
 
    USE MOD_SPMD_Task
-   USE MOD_LandPatch, only : numpatch
    USE MOD_LandPFT,   only : numpft
    USE MOD_Precision
    IMPLICIT NONE
@@ -61,9 +58,6 @@ SUBROUTINE allocate_PFTimeInvariants
             allocate (pftfrac       (numpft))
             allocate (htop_p        (numpft))
             allocate (hbot_p        (numpft))
-#ifdef CROP
-            allocate (cropfrac    (numpatch))
-#endif
          ENDIF
       ENDIF
 
@@ -72,7 +66,6 @@ END SUBROUTINE allocate_PFTimeInvariants
    SUBROUTINE READ_PFTimeInvariants (file_restart)
 
    USE MOD_NetCDFVector
-   USE MOD_LandPatch
    USE MOD_LandPFT
    IMPLICIT NONE
 
@@ -82,17 +75,13 @@ SUBROUTINE READ_PFTimeInvariants (file_restart)
       CALL ncio_read_vector (file_restart, 'pftfrac ', landpft, pftfrac ) !
       CALL ncio_read_vector (file_restart, 'htop_p  ', landpft, htop_p  ) !
       CALL ncio_read_vector (file_restart, 'hbot_p  ', landpft, hbot_p  ) !
-#ifdef CROP
-      CALL ncio_read_vector (file_restart, 'cropfrac ', landpatch, cropfrac) !
-#endif
 
    END SUBROUTINE READ_PFTimeInvariants
 
    SUBROUTINE WRITE_PFTimeInvariants (file_restart)
 
    USE MOD_NetCDFVector
    USE MOD_LandPFT
-   USE MOD_LandPatch
    USE MOD_Namelist
    USE MOD_Vars_Global
    IMPLICIT NONE
@@ -111,11 +100,6 @@ SUBROUTINE WRITE_PFTimeInvariants (file_restart)
       CALL ncio_write_vector (file_restart, 'htop_p  ', 'pft', landpft, htop_p  , compress) !
       CALL ncio_write_vector (file_restart, 'hbot_p  ', 'pft', landpft, hbot_p  , compress) !
 
-#ifdef CROP
-      CALL ncio_define_dimension_vector (file_restart, landpatch, 'patch')
-      CALL ncio_write_vector (file_restart, 'cropfrac', 'patch', landpatch, cropfrac, compress) !
-#endif
-
    END SUBROUTINE WRITE_PFTimeInvariants
 
    SUBROUTINE deallocate_PFTimeInvariants
@@ -131,9 +115,6 @@ SUBROUTINE deallocate_PFTimeInvariants
             deallocate (pftfrac )
             deallocate (htop_p  )
             deallocate (hbot_p  )
-#ifdef CROP
-            deallocate (cropfrac)
-#endif
          ENDIF
       ENDIF
 
@@ -143,13 +124,18 @@ END SUBROUTINE deallocate_PFTimeInvariants
    SUBROUTINE check_PFTimeInvariants ()
 
    USE MOD_RangeCheck
+#ifdef CROP
+   USE MOD_LandPatch
+#endif
    IMPLICIT NONE
 
       CALL check_vector_data ('pftfrac', pftfrac) !
       CALL check_vector_data ('htop_p ', htop_p ) !
       CALL check_vector_data ('hbot_p ', hbot_p ) !
 #ifdef CROP
-      CALL check_vector_data ('cropfrac', cropfrac) !
+      IF (landpatch%has_shared) THEN
+         CALL check_vector_data ('pct crop', landpatch%pctshared) !
+      ENDIF
 #endif
 
    END SUBROUTINE check_PFTimeInvariants
@@ -255,6 +241,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 +336,13 @@ 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)
@@ -382,13 +383,11 @@ SUBROUTINE READ_TimeInvariants (lc_year, casename, dir_restart)
    integer         , intent(in) :: lc_year
    character(len=*), intent(in) :: casename
    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'
-
       CALL ncio_read_vector (file_restart, 'patchclass',   landpatch, patchclass)          !
       CALL ncio_read_vector (file_restart, 'patchtype' ,   landpatch, patchtype )          !
       CALL ncio_read_vector (file_restart, 'patchmask' ,   landpatch, patchmask )          !
@@ -468,6 +467,26 @@ 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)
 
+      !-------------------------------------------------------------------------------------------------
+      ! Used for downscaling
+      !-------------------------------------------------------------------------------------------------
+      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
+
+
+
 #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)
@@ -664,6 +683,7 @@ SUBROUTINE deallocate_TimeInvariants ()
 
    USE MOD_SPMD_Task
    USE MOD_LandPatch, only: numpatch
+   USE MOD_Namelist
 
    IMPLICIT NONE
 
@@ -736,6 +756,14 @@ 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    )
+         ENDIF
+
          ENDIF
       ENDIF
 
@@ -758,7 +786,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
 
    IMPLICIT NONE
 
@@ -817,6 +845,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 ('sf_lut_patches      [-] '    , sf_lut_patches)       ! shadow mask
+ENDIF
+
 #ifdef USEMPI
       CALL mpi_barrier (p_comm_glb, p_err)
 #endif