Merge pull request CoLM-SYSU#296 from leelew/master

Add machine learning model for downscaling module

include/define.h

@@ -73,3 +73,5 @@
 #ifdef VectorInOneFileP
 #undef VectorInOneFileS
 #endif
+
+#define USESplitAI

main/CoLM.F90

@@ -132,7 +132,22 @@ PROGRAM CoLM
    integer*8 :: start_time, end_time, c_per_sec, time_used
 
 #ifdef USEMPI
+#ifdef USESplitAI
+      integer :: num_procs, my_rank, ierr, color, new_comm
+      CALL MPI_Init(ierr) ! Initialize MPI
+      CALL MPI_Comm_size(MPI_COMM_WORLD, num_procs, ierr) ! Get the total number of processes
+      CALL MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr) ! Get the rank of the current process
+      color = 1 ! The pyroot process will be in its own communicator
+      print*, 'before split I am process', my_rank, 'of', num_procs
+      CALL MPI_Comm_split(MPI_COMM_WORLD, color, my_rank, new_comm, ierr) ! Split the communicator
+      print*, 'after split I am process', my_rank, 'of', num_procs
+      CALL MPI_Comm_size(new_comm, num_procs, ierr) ! Get the total number of processes
+      CALL MPI_Comm_rank(new_comm, my_rank, ierr) ! Get the rank of the current process
+      print*,num_procs,"for CoLM"
+      CALL spmd_init (new_comm)
+#else
       CALL spmd_init ()
+#endif
 #endif
 
       CALL getarg (1, nlfile)

main/MOD_Forcing.F90

@@ -414,7 +414,9 @@ SUBROUTINE read_forcing (idate, dir_forcing)
    INTEGER  :: year, month, mday
    logical  :: has_u,has_v
    real solar, frl, prcp, tm, us, vs, pres, qm
-   real(r8) :: pco2m
+   real(r8) :: pco2m                        
+   real(r8), dimension(12, numpatch) :: spaceship !NOTE: 12 is the dimension size of spaceship
+   integer target_server, ierr
 
       IF (p_is_io) THEN
          !------------------------------------------------------------
@@ -709,8 +711,7 @@ SUBROUTINE read_forcing (idate, dir_forcing)
          CALL mg2p_forc%grid2part (forc_xy_us,      forc_us_grid   )
          CALL mg2p_forc%grid2part (forc_xy_vs,      forc_vs_grid   )
 
-         calday = calendarday(idate)
-         write(*,*) 'calday', calday
+         calday = calendarday(idate) 
 
          IF (p_is_worker) THEN
             DO np = 1, numpatch ! patches
@@ -780,10 +781,6 @@ SUBROUTINE read_forcing (idate, dir_forcing)
             ENDDO
          ENDIF
 
-         ! Conservation of short- and long- waves radiation in the grid of forcing
-         CALL mg2p_forc%normalize (forc_xy_solarin, forc_swrad_part)
-         CALL mg2p_forc%normalize (forc_xy_frl,     forc_frl_part  )
-
          ! mapping parts to patches
          CALL mg2p_forc%part2pset (forc_t_part,      forc_t     )
          CALL mg2p_forc%part2pset (forc_q_part,      forc_q     )
@@ -796,7 +793,61 @@ SUBROUTINE read_forcing (idate, dir_forcing)
          CALL mg2p_forc%part2pset (forc_us_part,     forc_us    )
          CALL mg2p_forc%part2pset (forc_vs_part,     forc_vs    )
 
-         ! divide fractions of downscaled shortwave radiation
+#ifndef SinglePoint
+         IF (trim(DEF_DS_precipitation_adjust_scheme) == 'III') THEN 
+            ! Sisi Chen, Lu Li, Yongjiu Dai et al., 2024, JGR
+            ! Using MPI to pass the forcing variable field to Python to accomplish precipitation downscaling
+            IF (p_is_worker) THEN
+               spaceship(1,1:numpatch) = forc_topo
+               spaceship(2,1:numpatch) = forc_t
+               spaceship(3,1:numpatch) = forc_pbot
+               spaceship(4,1:numpatch) = forc_q
+               spaceship(5,1:numpatch) = forc_frl
+               spaceship(6,1:numpatch) = forc_swrad
+               spaceship(7,1:numpatch) = forc_us
+               spaceship(8,1:numpatch) = forc_vs
+               spaceship(9,1:numpatch) = INT(calday)
+               spaceship(10,1:numpatch) = patchlatr
+               spaceship(11,1:numpatch) = patchlonr
+
+               target_server = p_iam_glb/5+p_np_glb
+               CALL MPI_SEND(spaceship,12*numpatch,MPI_REAL8,target_server,0,MPI_COMM_WORLD,ierr)
+               CALL MPI_RECV(forc_prc,numpatch,MPI_REAL8,target_server,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
+
+               forc_prl = forc_prc/3600*2/3._r8
+               forc_prc = forc_prc/3600*1/3._r8
+            ENDIF
+
+            ! mapping forc_prl to forc_prl_part, forc_prc to forc_prc_part
+            IF (p_is_worker) THEN
+               DO np = 1, numpatch ! patches
+                  DO ipart = 1, mg2p_forc%npart(np) ! part loop of each patch
+                     IF (mg2p_forc%areapart(np)%val(ipart) == 0.) CYCLE
+
+                     forc_prl_part(np)%val(ipart) = forc_prl(np)
+                     forc_prc_part(np)%val(ipart) = forc_prc(np)
+
+                  ENDDO
+               ENDDO
+            ENDIF
+
+            ! Conservation of convective and large scale precipitation in the grid of forcing
+            CALL mg2p_forc%normalize (forc_xy_prc, forc_prc_part)
+            CALL mg2p_forc%normalize (forc_xy_prl, forc_prl_part)
+
+            ! mapping parts to patches
+            CALL mg2p_forc%part2pset (forc_prc_part, forc_prc)
+            CALL mg2p_forc%part2pset (forc_prl_part, forc_prl)
+         ENDIF
+
+         ! Conservation of short- and long- waves radiation in the grid of forcing
+         CALL mg2p_forc%normalize (forc_xy_solarin, forc_swrad_part)
+         CALL mg2p_forc%normalize (forc_xy_frl,     forc_frl_part  )
+         CALL mg2p_forc%part2pset (forc_frl_part,    forc_frl   )
+         CALL mg2p_forc%part2pset (forc_swrad_part,  forc_swrad )
+#endif
+
+         ! divide fractions of downscaled shortwave radiation 
          IF (p_is_worker) THEN
             DO j = 1, numpatch
                   a = forc_swrad(j)

main/MOD_ForcingDownscaling.F90

@@ -119,7 +119,7 @@ SUBROUTINE downscale_forcings (&
    IMPLICIT NONE
 
    integer,  parameter :: S = 1370            ! solar constant (W/m**2)
-   real(r8), parameter :: thr = 85*PI/180      ! threshold of ??
+   real(r8), parameter :: thr = 85*PI/180      ! threshold of zenith angle
 
    ! ARGUMENTS:
    logical,  intent(in) :: glaciers            ! true: glacier column (itypwat = 3)
@@ -132,7 +132,7 @@ SUBROUTINE downscale_forcings (&
    real(r8), intent(in) :: svf_c                                    ! sky view factor
    real(r8), intent(in) :: cur_c                                    ! curvature
    real(r8), intent(in) :: sf_lut_c   (1:num_azimuth,1:num_zenith)  ! look up table of shadow mask of a patch
-   real(r8), intent(in) :: asp_type_c (1:num_type)                  ! topographic aspect of each type of one patch
+   real(r8), intent(in) :: asp_type_c (1:num_type)                  ! topographic aspect of each type of one patch(rad)
    real(r8), intent(in) :: slp_type_c (1:num_type)                  ! topographic slope of each character of one patch
    real(r8), intent(in) :: area_type_c(1:num_type)                  ! area percentage of each character of one patch
 
@@ -223,7 +223,7 @@ SUBROUTINE downscale_forcings (&
       ! save
       forc_t_c    = tbot_c
       forc_th_c   = thbot_c
-      forc_q_c    = qbot_c
+      forc_q_c    = qbot_c  
       forc_pbot_c = pbot_c
       forc_rho_c  = rhos_c
 
@@ -265,22 +265,15 @@ SUBROUTINE downscale_forcings (&
          ! Liston, G. E. and Elder, K.: A meteorological distribution system
          ! for high-resolution terrestrial modeling (MicroMet), J. Hydrometeorol., 7, 217-234, 2006.
          ! Equation (33) and Table 1: chi range from January to December:
-         ! [0.35,0.35,0.35,0.30,0.25,0.20,0.20,0.20,0.20,0.25,0.30,0.35] (1/m)
+         ! [0.35,0.35,0.35,0.30,0.25,0.20,0.20,0.20,0.20,0.25,0.30,0.35] (1/km)
 
-         delta_prc_c = forc_prc_g * 2.0*0.27e-3*(forc_topo_c - forc_topo_g) &
-            /(1.0 - 0.27e-3*(forc_topo_c - forc_topo_g))
+         delta_prc_c = forc_prc_g *1.0*0.27*(forc_topo_c - forc_topo_g) &
+            /(1.0 - 0.27*(forc_topo_c - forc_topo_g))
          forc_prc_c = forc_prc_g + delta_prc_c   ! large scale precipitation [mm/s]
 
-         delta_prl_c = forc_prl_g * 2.0*0.27e-3*(forc_topo_c - forc_topo_g) &
-            /(1.0 - 0.27e-3*(forc_topo_c - forc_topo_g))
+         delta_prl_c = forc_prl_g *1.0*0.27*(forc_topo_c - forc_topo_g) &
+            /(1.0 - 0.27*(forc_topo_c - forc_topo_g))
          forc_prl_c = forc_prl_g + delta_prl_c   ! large scale precipitation [mm/s]
-
-      ELSEIF (trim(DEF_DS_precipitation_adjust_scheme) == 'III') THEN 
-         ! Mei, Y., Maggioni, V., Houser, P., Xue, Y., & Rouf, T. (2020). A nonparametric statistical 
-         ! technique for spatial downscaling of precipitation over High Mountain Asia. Water Resources Research, 
-         ! 56, e2020WR027472. https://doi.org/ 10.1029/2020WR027472
-
-         ! We implement this scheme in MOD_Forcing.F90
       END IF
 
       IF (forc_prl_c < 0) THEN
@@ -334,7 +327,7 @@ SUBROUTINE downscale_wind(forc_us_g, forc_vs_g, &
 
       ! calculate wind direction
       IF (forc_us_g == 0.) THEN
-         wind_dir  = 0
+         wind_dir  = PI/2
       ELSE
          wind_dir  = atan(forc_vs_g /forc_us_g)
       ENDIF
@@ -344,7 +337,7 @@ SUBROUTINE downscale_wind(forc_us_g, forc_vs_g, &
 
       ! compute the slope in the direction of the wind
       DO i = 1, num_type
-         wind_dir_slp(i) = slp_type_c(i)*cos(wind_dir-asp_type_c(i)*PI/180)
+         wind_dir_slp(i) = slp_type_c(i)*cos(wind_dir-asp_type_c(i))
       ENDDO
 
       ! compute wind speed ajustment
@@ -482,7 +475,7 @@ SUBROUTINE downscale_shortwave( &
    IMPLICIT NONE
 
    integer,  parameter :: S = 1370                               ! solar constant (W/m**2)
-   real(r8), parameter :: thr = 85*PI/180                        ! threshold of ??
+   real(r8), parameter :: thr = 85*PI/180                        ! threshold of zenith
    real(r8), parameter :: shortwave_downscaling_limit = 0.5_r8   ! relative limit for how much shortwave downscaling can be done (unitless)
 
    ! ARGUMENTS:
@@ -519,7 +512,7 @@ SUBROUTINE downscale_shortwave( &
    real(r8) :: svf, balb
 
    real(r8) :: diff_swrad_g, beam_swrad_g              ! diffuse and beam radiation
-   real(r8) :: diff_swrad_c, beam_swrad_c, refl_swrad_c! downscaled diffuse and beam radiation 
+   real(r8) :: diff_swrad_c, beam_swrad_c, refl_swrad_c! downscaled diffuse, beam radiation and reflect radiation 
    real(r8) :: beam_swrad_type (1:num_type)            ! beam radiation of one characterized patch
    real(r8) :: refl_swrad_type (1:num_type)            ! reflect radiation of one characterized patch
    real(r8) :: tcf_type        (1:num_type)            ! terrain configure factor
@@ -539,6 +532,7 @@ SUBROUTINE downscale_shortwave( &
       idx_zen = INT(zen_deg*num_zenith/90)
       IF (idx_azi==0) idx_azi = 1
       IF (idx_zen==0) idx_zen = 1
+      IF (idx_zen>num_zenith) idx_zen = num_zenith !constrain the upper boundary of zenith angle to 90 deg
 
       sf_c = sf_lut_c(idx_azi, idx_zen)
       IF (sf_c<0) sf_c = 0
@@ -549,8 +543,8 @@ SUBROUTINE downscale_shortwave( &
       toa_swrad = S*(rt_R**2)*coszen
 
       ! calculate clearness index
-      IF (toa_swrad.le.0) THEN
-         clr_idx = 1
+      IF (toa_swrad.eq.0) THEN
+         clr_idx = 0
       ELSE
          clr_idx = forc_swrad_g/toa_swrad
       ENDIF
@@ -589,7 +583,7 @@ SUBROUTINE downscale_shortwave( &
          ! calculate the cosine of solar illumination angle, cos(θ), 
          ! ranging between −1 and 1, indicates if the sun is below or 
          ! above the local horizon (note that values lower than 0 are set to 0 indicate self shadow)
-         cosill_type(i) = cos(slp_type_c(i))+tan(zen_rad)*sin(slp_type_c(i))*cos(asp_type_c(i)*PI/180)
+         cosill_type(i) = cos(slp_type_c(i))+tan(zen_rad)*sin(slp_type_c(i))*cos(asp_type_c(i))
          IF (cosill_type(i)>1) cosill_type(i) = 1
          IF (cosill_type(i)<0) cosill_type(i) = 0
 
@@ -633,7 +627,7 @@ SUBROUTINE downscale_shortwave( &
                forc_swrad_g * (1._r8 + shortwave_downscaling_limit))
       forc_swrad_c = max(forc_swrad_c, &
                forc_swrad_g * (1._r8 - shortwave_downscaling_limit))
-      ! for normalize
+      ! Ensure that the denominator is not 0 during shortwave normalization
       IF (forc_swrad_c==0.) forc_swrad_c = 0.0001
 
    END SUBROUTINE downscale_shortwave

main/MOD_OrbCoszen.F90

@@ -72,9 +72,6 @@ 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

mksrfdata/Aggregation_TopographyFactors.F90

@@ -36,8 +36,8 @@ SUBROUTINE Aggregation_TopographyFactors ( &
 
    ! local variables:
    ! ---------------------------------------------------------------
-   CHARACTER(len=256) :: landdir, lndname, cyear
-   CHARACTER(len=3)   :: sdir
+   CHARACTER(len=256) :: landdir, lndname, cyear                                                                   
+   CHARACTER(len=3)   :: sdir, sdir1
 
    TYPE (block_data_real8_2d) :: slp_grid    ! slope
    TYPE (block_data_real8_2d) :: asp_grid    ! aspect
@@ -83,7 +83,7 @@ SUBROUTINE Aggregation_TopographyFactors ( &
    REAL(r8) :: zenith_angle(num_zenith)    ! sine of sun zenith angle (divided by num_zenith part)
 
    ! local variables
-   INTEGER :: ipatch, i, ps, pe, type, a, z, count_pixels, num_pixels
+   INTEGER :: ipatch, i, ps, pe, type, a, z, count_pixels, num_pixels, j  
 
 #ifdef SrfdataDiag
    INTEGER :: typpatch(N_land_classification+1), ityp  ! number of land classification
@@ -209,11 +209,11 @@ SUBROUTINE Aggregation_TopographyFactors ( &
          sum_area_one = sum(area_one, mask = area_one>0)
 
          DO a = 1, num_azimuth
-            ! terrain elevation angle at each azimuth
-            tea_f_one(:) = tea_f_azi_one(a,:)
-            tea_b_one(:) = tea_b_azi_one(a,:)
-
             DO z = 1, num_zenith
+               ! terrain elevation angle at each azimuth
+               tea_f_one(:) = tea_f_azi_one(a,:)
+               tea_b_one(:) = tea_b_azi_one(a,:) 
+
                ! count the pixels which are not missing value
                count_pixels = 0
 
@@ -280,19 +280,19 @@ SUBROUTINE Aggregation_TopographyFactors ( &
 
          DO i = 1, num_pixels
             ! Define the south slope, north slope, abrupt slope and gentle lope of target pixel
-            IF ((asp_one(i).ge.0 .and. asp_one(i).le.90) .or. (asp_one(i).ge.270 .and. asp_one(i).le.360).and.(slp_one(i).ge.15*pi/180)) THEN  ! north abrupt slope
-                 type = 1
-            ELSE IF ((asp_one(i).ge.0 .and. asp_one(i).le.90) .or. (asp_one(i).ge.270 .and. asp_one(i).le.360).and.(slp_one(i)<15*pi/180)) THEN  ! north gentle slope
+            IF ((asp_one(i).ge.0 .and. asp_one(i).le.90*pi/180) .or. (asp_one(i).ge.270*pi/180 .and. asp_one(i).le.360*pi/180).and.(slp_one(i).ge.15*pi/180)) THEN  ! north abrupt slope
+                 type = 1 
+            ELSE IF ((asp_one(i).ge.0 .and. asp_one(i).le.90*pi/180) .or. (asp_one(i).ge.270*pi/180 .and. asp_one(i).le.360*pi/180).and.(slp_one(i)<15*pi/180)) THEN  ! north gentle slope
                  type = 2
-            ELSE IF ((asp_one(i).gt.90) .and. (asp_one(i).lt.270) .and. (slp_one(i).ge.15*pi/180)) THEN  ! south abrupt slope
-                 type = 3
-            ELSE IF ((asp_one(i).gt.90) .and. (asp_one(i).lt.270) .and. (slp_one(i).lt.15*pi/180)) THEN  ! south gentle slope
-                 type = 4
+            ELSE IF ((asp_one(i).gt.90*pi/180) .and. (asp_one(i).lt.270*pi/180) .and. (slp_one(i).ge.15*pi/180)) THEN  ! south abrupt slope
+                 type = 3 
+            ELSE IF ((asp_one(i).gt.90*pi/180) .and. (asp_one(i).lt.270*pi/180) .and. (slp_one(i).lt.15*pi/180)) THEN  ! south gentle slope
+                 type = 4  
             ELSE ! missing value=-9999
                  cycle
             END IF
 
-            IF ((area_one(i)>0).and.(area_one(i)<sum_area_one)) THEN      ! quality control
+            IF ((area_one(i)>0).and.(area_one(i)<=sum_area_one)) THEN      ! quality control
                   area_type_one(type,i) = area_one(i)
                   asp_type_one (type,i) = asp_one(i)*area_one(i)
                   slp_type_one (type,i) = slp_one(i)*area_one(i)
@@ -381,24 +381,39 @@ SUBROUTINE Aggregation_TopographyFactors ( &
 
 #ifdef SrfdataDiag
    typpatch = (/(ityp, ityp = 0, N_land_classification)/)
-   lndname  = trim(dir_model_landdata) // '/diag/topo_factor_' // trim(cyear) // '.nc'
 
    ! only write the first type of slope and aspect at patches
-   CALL srfdata_map_and_write (slp_type_patches(1,:), landpatch%settyp, typpatch, m_patch2diag, &
-      -1.0e36_r8, lndname, 'slp', compress = 1, write_mode = 'one')
-   CALL srfdata_map_and_write (asp_type_patches(1,:), landpatch%settyp, typpatch, m_patch2diag, &
-      -1.0e36_r8, lndname, 'asp', compress = 1, write_mode = 'one')
+   lndname  = trim(dir_model_landdata) // '/diag/topo_factor_slp_' // trim(cyear) // '.nc'
+   DO i = 1, num_type
+      write(sdir,'(I0)') i
+      CALL srfdata_map_and_write (slp_type_patches(i,:), landpatch%settyp, typpatch, m_patch2diag, &
+         -1.0e36_r8, lndname, 'slp_'//trim(sdir), compress = 1, write_mode = 'one')
+   ENDDO
+
+   lndname  = trim(dir_model_landdata) // '/diag/topo_factor_asp_' // trim(cyear) // '.nc'
+   DO i = 1, num_type
+      write(sdir,'(I0)') i
+      CALL srfdata_map_and_write (asp_type_patches(i,:), landpatch%settyp, typpatch, m_patch2diag, &
+         -1.0e36_r8, lndname, 'asp_'//trim(sdir), compress = 1, write_mode = 'one')
+   ENDDO
 
+   lndname  = trim(dir_model_landdata) // '/diag/topo_factor_svf_' // trim(cyear) // '.nc'
    CALL srfdata_map_and_write (svf_patches, landpatch%settyp, typpatch, m_patch2diag, &
       -1.0e36_r8, lndname, 'svf', compress = 1, write_mode = 'one')
+
+   lndname  = trim(dir_model_landdata) // '/diag/topo_factor_cur_' // trim(cyear) // '.nc'
    CALL srfdata_map_and_write (cur_patches, landpatch%settyp, typpatch, m_patch2diag, &
       -1.0e36_r8, lndname, 'cur', compress = 1, write_mode = 'one')
 
-   ! only write part of lut to save cache
-   DO i = 1, 11
-      write(sdir,'(I0)') i
-      CALL srfdata_map_and_write (sf_lut_patches(1,i,:), landpatch%settyp, typpatch, m_patch2diag, &
-         -1.0e36_r8, lndname, 'sf_'//trim(sdir), compress = 1, write_mode = 'one')
+   lndname  = trim(dir_model_landdata) // '/diag/topo_factor_sf_lut_' // trim(cyear) // '.nc'
+
+   DO j = 1, num_azimuth
+      DO i = 1, num_zenith 
+         write(sdir,'(I0)') j
+         write(sdir1,'(I0)') i
+         CALL srfdata_map_and_write (sf_lut_patches(j,i,:), landpatch%settyp, typpatch, m_patch2diag, &
+            -1.0e36_r8, lndname, 'sf_'//trim(sdir)//'_'//trim(sdir1), compress = 1, write_mode = 'one')
+      ENDDO
    ENDDO
 #endif
 #else