Cleaned code and updated comments

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/GEOS_IrrigationGridComp.F90

@@ -560,11 +560,10 @@ subroutine INITIALIZE (GC,IMPORT, EXPORT, CLOCK, RC )
     call MAPL_GetPointer(EXPORT, FLOODRATE,     'FLOODRATE',    ALLOC=.true., RC=STATUS) ; VERIFY_(STATUS)
 
     ! Update IRRIGFRAC and PADDYFRAC for applications that are run on regular tiles in which IRRIGFRAC and PADDYFRAC in BCs are fractions.
-    ! The irrigation model would run on tiles whose IRRIGFRAC + PADDYFRAC > IRRIG_THRES (defult is 0.5) assuming the larger 
-    ! of the two fractions is the dominant surface type.
+    ! The irrigation model would run on tiles whose IRRIGFRAC + PADDYFRAC > IRRIG_THRES (default is 0.5).
 
     where (IRRIGFRAC + PADDYFRAC > IM%IRRIG_THRES)
-!To assign the entire cell to the largest fraction remove comments 568-574
+!To assign the entire cell to the largest fraction remove comments 568-574,
         ! where (PADDYFRAC >= IRRIGFRAC)
         !    PADDYFRAC = 1.
         !    IRRIGFRAC = 0.
@@ -597,8 +596,7 @@ subroutine INITIALIZE (GC,IMPORT, EXPORT, CLOCK, RC )
 
     ! Scale computed SPRINKLERRATE, DRIPRATE, and FLOODRATE to the total
     ! irrigated tile fraction before exporting to land models.
-    ! Since revised IRRIGFRAC, and PADDYFRAC in subtiling mode are 0. or 1., below scaling
-    ! has no effect in that mode.
+
 
     SPRINKLERRATE = SPRINKLERRATE*(IRRIGFRAC + PADDYFRAC)
     DRIPRATE      = DRIPRATE     *(IRRIGFRAC + PADDYFRAC)
@@ -836,9 +834,7 @@ subroutine RUN (GC,IMPORT, EXPORT, CLOCK, RC )
 
     ! Scale computed SPRINKLERRATE, DRIPRATE, and FLOODRATE to the total
     ! irrigated tile fraction before exporting to land models.
-    ! Since revised IRRIGFRAC, and PADDYFRAC in subtiling mode are 0. or 1., below scaling
-    ! has no effect in that mode. 
-
+
     SPRINKLERRATE = SPRINKLERRATE*(IRRIGFRAC + PADDYFRAC)
     DRIPRATE      = DRIPRATE     *(IRRIGFRAC + PADDYFRAC)
     FLOODRATE     = FLOODRATE    *(IRRIGFRAC + PADDYFRAC)

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSirrigation_GridComp/irrigation_model.F90

@@ -21,22 +21,17 @@ MODULE IRRIGATION_MODULE
   ! Drip irrigation method calculation is similar to that of sprinkler, albeit the drip irrigation method assumes a 0% water loss.
   !
   ! March 21, 2021 (Sarith Mahanama) - First Version
+  ! June 25, 2024 (Stefano Casirati) - Second Version
   !
   ! (1) EXPORTS - MODEL OUTPUTS TO THE LAND MODEL (IRRIGATION RATES):
   !    1) SPRINKLERRATE [kg m-2 s-1]
   !    2) DRIPRATE [kg m-2 s-1]
   !    3) FLOODRATE [kg m-2 s-1]    
   !
-  ! (2) COMPUTATIONAL TILES:
-  !    During land BCs generation, land tiles where irrigated crops or paddy is present were further subdivided to a
-  !    mosaic of upto 3 fractions: i) non-irrigated land, ii) irrigated crop, or iii) paddy.
-  !    The model treats each fraction as a separate computational tile and runs on each individual fraction with own parameters and prognostics.
-  !    All fractions inherited model and soil parameters from the main land tile that they belong to. A special treatment of setting BF3 to a high
-  !    value (25.) was applied to paddy/crop tiles to account for the uniquely flat nature of farmlands. Vegetation characteristics and vegetation dynamic
-  !    parameters for irrigated crop and paddy tiles were taken from the nearest grass or crops land tile. 
-  !    During tiling and BCs data preparation, computed fractional coverages for land tiles were also adjusted
-  !    to reflect each computational tile under the land grid component represents entirely one of the 3 irrigation surface types: a non-irrigated land,
-  !    OR a irrigated-crops OR a paddy tile.
+  ! (2) IRRIGATED AND PADDY TILES:
+  ! During land BC's generation, the fraction of irrigated crops and paddy is set to zero if their sum is below an irrigation threshold.
+  ! Irrigated fractions can be irrigated with sprinkler, drip, and flood, while paddy fractions can only be irrigated using the flood irrigation method.
+  ! Vegetation characteristics and vegetation dynamic parameters for irrigated crops and paddy tiles were taken from the nearest grass or cropland tile.
   !
   ! (3) MODEL INPUTS:
   !     SMWP    : rootzone soil moisture content at wilting point [mm] 
@@ -61,8 +56,10 @@ MODULE IRRIGATION_MODULE
   !          iii)3: Flood irrigation on entire tile.
   !
   !     IRRIG_TRIGGER: 0 SPECIFIC INPUTS:
-  !          IRRIGFRAC    : fraction of tile covered by irrigated crops (per Section 2, values will be 0. or 1.)
-  !          PADDYFRAC    : fraction of tile covered by paddy (per Section 2, values will be 0. or 1.)
+  !          IRRIGFRAC    : fraction of tile covered by irrigated crops (values between 0 and 1 (if IRRIGFRAC + PADDYFRAC > Irrigation
+  !                         Threshold)
+  !          PADDYFRAC    : fraction of tile covered by paddy (values between 0 and 1 (if IRRIGFRAC + PADDYFRAC > Irrigation
+  !                         Threshold)
   !          SPRINKLERFR  : fraction of tile equipped for sprinkler irrigation
   !          DRIPFR       : fraction of tile equipped for drip irrigation
   !          FLOODFR      : fraction of tile equipped for flood irrigation
@@ -102,12 +99,6 @@ MODULE IRRIGATION_MODULE
   !     The second dimensions of 2D arrays is for different crop fractions i.e. the second dimension is 2 for above
   !     IRRIG_TRIGGER: 0 to separately store irrigation rates in irrigated crop and paddy fractions.
   !     It would be 26 for IRRIG_TRIGGER: 1.
-  !     Note also that runnning the irrigation model on subtiling mode (specific irrigated crop and paddy tiles with their own land prognostics)
-  !     is preferred (Section 2) to running the irrigation model on fractions of typical land tiles. For the subtiling mode with own soil moisture
-  !     prognostics, IRRIGFRAC, PADDYFRAC and CROPIRRIGFRAC fractions have been adjusted to represent irrigation type of the subtile in question. i.e.  
-  !     IRRIFRAC has been set to 1. on irrigated-crop subtiles; PADDYFRAC and CROPIRRIGFRAC(N,3) are set to 1. on paddy subtiles; and  SUM of
-  !     CROPIRRIGFRAC(N,:) excluding the 3rd element is set to 1. on irrigated crop tiles.
-  !
   !     The crop calendar implemetation (IRRIG_TRIGGER: 1) computes SPRINKLERRATE, DRIPRATE, and FLOODRATE as weighted averages of irrigation rates from
   !     all active crops in SRATE, DRATE and FRATE arrays.
 

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/module_irrig_params.F90

@@ -2,6 +2,8 @@
 #define ASSERT_(A)   if(.not.A)then;print *,'Error:',__FILE__,__LINE__;stop;endif
 
 module module_irrig_params
+! Version 1 - Sarith Mahanama
+! Version 2 - Stefano Casirati - LAI min-max obtained from LAI climatology boundary conditions   
 
   use rmTinyCatchParaMod,  ONLY : RegridRaster,regridrasterreal
   use process_hres_data,   ONLY : get_country_codes
@@ -143,7 +145,7 @@ subroutine create_irrig_params (nc, nr, gfile)
     allocate (GIAFRAC  (NTILES))
     call ReadProcess_GIA (NC, NR, NTILES, tile_id, GIAFRAC)
 
-    ! (4) Process GRIPC and MCD15A2H LAI
+    ! (4) Process GRIPC and LAI (Min-Max)
     ! ----------------------------------
 
     allocate (IGRIPC  (NTILES))
@@ -329,8 +331,8 @@ SUBROUTINE OpenFile
              date_time_values(1),'-',date_time_values(2),'-',date_time_values(3),'at', &
              date_time_values(5),':',date_time_values(6),':',date_time_values(7)
 
-        status = NF_PUT_ATT_TEXT(NCOutID, NF_GLOBAL, 'CreatedBy', LEN_TRIM('Sarith Mahanama'),   &
-             'Sarith Mahanama') ; VERIFY_(STATUS) 
+        status = NF_PUT_ATT_TEXT(NCOutID, NF_GLOBAL, 'CreatedBy', LEN_TRIM('Sarith Mahanama, Stefano Casirati'),   &
+             'Sarith Mahanama, Stefano Casirati') ; VERIFY_(STATUS) 
         status = NF_PUT_ATT_TEXT(NCOutID, NF_GLOBAL, 'Contact'   , LEN_TRIM('sarith.p.mahanama@nasa.gov'),                  &
              'sarith.p.mahanama@nasa.gov')
         status = NF_PUT_ATT_TEXT(NCOutID, NF_GLOBAL, 'Date'   , LEN_TRIM(time_stamp),trim(time_stamp))
@@ -387,7 +389,7 @@ SUBROUTINE MergeData (NTILES)
         ALLOCATE (DRIP     (1:NTILES))
         CALL ReadProcess_IMethod (NTILES, sprinkler, drip, flood)        
 
-        ! MERGING PROCEDURE
+        ! MERGING PROCEDURE  
         ! =================
         !                                                                                                                                  CROP CALENDAR
         ! GIA                 GIA-GRIPC                           GIA-GRIPC-MIRCA                                                          -------------
@@ -1271,50 +1273,8 @@ SUBROUTINE ReadProcess_GRIPC (NC, NR, NTILES, tile_id,IGRIPC, RGRIPC, PGRIPC, NG
       real,allocatable :: dum,yr,mn,dy,nt
       logical                                 :: write_lai = .false.
 
-    ! if (write_lai) then
-
-         ! Process LAI first
-         ! -----------------
-
-    !     status = NF_CREATE ('land/irrigation/lai_clim_min_max/v1/MCD15A2H.006_LAI_climMinMax.nc4' , NF_NETCDF4, NCIDW) ; VERIFY_(STATUS)
-    !     status = NF_DEF_DIM(NCIDW, 'lon'        , NX_LAI, xid)    ; VERIFY_(STATUS)
-    !     status = NF_DEF_DIM(NCIDW, 'lat'        , NY_LAI, yid)    ; VERIFY_(STATUS)
-    !     status = NF_DEF_VAR(NCIDW, 'LAIMIN'   , NF_FLOAT, 2 ,(/xid,yid/), vid) ; VERIFY_(STATUS)
-    !     status = NF_DEF_VAR(NCIDW, 'LAIMAX'   , NF_FLOAT, 2 ,(/xid,yid/), vid) ; VERIFY_(STATUS)
-    !     status = NF_ENDDEF(NCIDW)
-
-    !     allocate (Lai_clim  (1:NX_LAI, 1: NY_LAI))
-    !     allocate (clim_min  (1:NX_LAI, 1: NY_LAI))
-    !     allocate (clim_max  (1:NX_LAI, 1: NY_LAI))
-    !     allocate (clim_lai  (1:NX_LAI, 1: NY_LAI))
-    !     clim_max = -9999.
-    !     clim_min =  9999.
-
-    !     DO DOY = 1, 361, 8
-    !        write (DDD, '(i3.3)') DOY
-    !        print *,trim(LAI_file)//DDD//'.nc4'
-    !        status = NF_OPEN(trim(LAI_file)//DDD//'.nc4',NF_NOWRITE, ncid); VERIFY_(STATUS)
-    !        STATUS = NF_GET_VARA_INT2 (NCID,VarID(NCID,'Lai_500m'),(/1,1/),(/NX_LAI,NY_LAI/),Lai_clim) ; VERIFY_(STATUS)
-
-    !        clim_lai = -9999.
-    !        where ((Lai_clim >=0).and.(Lai_clim <= 100))
-    !           clim_lai = lai_clim * 0.1
-    !        end where
-    !        where ((clim_lai >= 0.) .and. (clim_lai > clim_max))
-    !           clim_max = clim_lai
-    !        end where
-    !        where ((clim_lai >= 0.) .and. (clim_lai < clim_min))
-    !           clim_min = clim_lai
-    !        end where
-    !        STATUS = NF_CLOSE (NCID)
-    !     END DO
-
-     !    STATUS = NF_PUT_VARA_REAL (NCIDW,VarID(NCIDW,'LAIMIN'),(/1,1/),(/NX_LAI,NY_LAI/),clim_min) ; VERIFY_(STATUS)
-     !    STATUS = NF_PUT_VARA_REAL (NCIDW,VarID(NCIDW,'LAIMAX'),(/1,1/),(/NX_LAI,NY_LAI/),clim_max) ; VERIFY_(STATUS)
-     !    STATUS = NF_CLOSE (NCIDW)
-     !   deallocate (lai_clim, clim_lai, clim_min, clim_max)
-     ! endif
-
+    !V2: Min Max LAI from LAI Climatology for consistence    
+
       allocate( var_in(NX_gripc,NY_gripc)) 
       var_in = UNDEF
 
@@ -1350,15 +1310,7 @@ SUBROUTINE ReadProcess_GRIPC (NC, NR, NTILES, tile_id,IGRIPC, RGRIPC, PGRIPC, NG
       allocate (lai (NTILES))
       LAI_MIN = -9999.
       LAI_MAX = -9999.
-!!!! Qui lo devo leggere .........  Ma se lo creo da lai.dat, Avra' dimensioni come le NTILES
-      ! Read Min/Max LAI
- !!!!     allocate (clim_min  (1:NX_LAI, 1: NY_LAI))
- !!!!     allocate (clim_max  (1:NX_LAI, 1: NY_LAI))
- !     status = NF_OPEN('/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/irrigation/lai_clim_min_max/v1/MCD15A2H.006_LAI_climMinMax.nc4', NF_NOWRITE, NCIDW); VERIFY_(STATUS)
-  !    STATUS = NF_GET_VARA_REAL (NCIDW,VarID(NCIDW,'LAIMIN'),(/1,1/),(/NX_LAI,NY_LAI/),clim_min)  ; VERIFY_(STATUS)
-  !    STATUS = NF_GET_VARA_REAL (NCIDW,VarID(NCIDW,'LAIMAX'),(/1,1/),(/NX_LAI,NY_LAI/),clim_max)  ; VERIFY_(STATUS)
-  !    STATUS = NF_CLOSE (NCIDW)
-
+
       do j = 1,NY_gripc
          do i =  1,NX_gripc
             if((iraster (i,j) >=1).and.(iraster (i,j) <=ntiles)) then
@@ -1368,19 +1320,6 @@ SUBROUTINE ReadProcess_GRIPC (NC, NR, NTILES, tile_id,IGRIPC, RGRIPC, PGRIPC, NG
                if (var_in(i,j) == 3) PGRIPC(iraster (i,j)) = PGRIPC(iraster (i,j)) + 1.
                if (var_in(i,j) == 4) NGRIPC(iraster (i,j)) = NGRIPC(iraster (i,j)) + 1.
 
-   !            if(var_in(i,j) == 2) then
-    !              if(clim_min(i,j) < 10.) then
-     !                if(LAI_MIN(iraster (i,j)) < 0.) LAI_MIN(iraster (i,j)) = 0.
-      !               LAI_MIN (iraster (i,j)) = LAI_MIN (iraster (i,j)) + clim_min(i,j)
-       !              min_cnt (iraster (i,j)) = min_cnt (iraster (i,j)) + 1.
-        !          endif
-
-   !               if(clim_max(i,j) >= 0.) then
-    !                 if(LAI_MAX(iraster (i,j)) < 0.) LAI_MAX(iraster (i,j)) = 0.
-     !                LAI_MAX (iraster (i,j)) = LAI_MAX (iraster (i,j)) + clim_max(i,j)
-      !               max_cnt (iraster (i,j)) = max_cnt (iraster (i,j)) + 1.
-       !           endif
-        !       endif
             endif
          end do
       end do
@@ -1390,19 +1329,6 @@ SUBROUTINE ReadProcess_GRIPC (NC, NR, NTILES, tile_id,IGRIPC, RGRIPC, PGRIPC, NG
       PGRIPC = PGRIPC / tot_cnt
       NGRIPC = NGRIPC / tot_cnt
 
-   !   where (max_cnt > 0.) LAI_MAX = LAI_MAX / max_cnt
-   !   where (min_cnt > 0.) LAI_MIN = LAI_MIN / min_cnt
-
-      ! Fill LAI gaps 
-      ! -------------
-      !      print *, 'START LAI gap filling',COUNT(LAI_MAX >=0.), MAXVAL(LAI_MAX),count (IGRIPC > 0.)
-   !   DO I = 1, NTILES
-   !      IF((IGRIPC (I) > 0.) .AND. (LAI_MAX(I) < 0.)) THEN
-   !         j = getNeighbor (I,lai_in=LAI_MAX)
-   !         LAI_MIN (I) = LAI_MIN (j)
-   !         LAI_MAX (I) = LAI_MAX (j)
-   !      ENDIF
-   !   END DO
    allocate(dum)
    allocate(yr)
    allocate(mn)
@@ -1429,7 +1355,7 @@ SUBROUTINE ReadProcess_GRIPC (NC, NR, NTILES, tile_id,IGRIPC, RGRIPC, PGRIPC, NG
 
       deallocate (var_in, iraster, min_cnt, max_cnt, tot_cnt, LAI_MIN, LAI_MAX, lai, yr, mn, dy, dum, nt)
 
-      print *,'DONE PROCESSING GRIPC and MCD15A2H LAI'
+      print *,'DONE PROCESSING GRIPC and LAI'
 
     END SUBROUTINE ReadProcess_GRIPC