Merge pull request #9 from noaa-oar-arl/develop

Merge Develop into Master

README.md

@@ -97,7 +97,11 @@ The variables listed here are set by the user in the NACC run script, and they a
 -   `IFVIIRS_GVF [default: False]`  
     Logical to choose if the vegetation fraction is taken from FV3 or VIIRS GVF.
     -   `False`: Use FV3 vegetation fraction
-    -   `True`: Use VIIRS GVF  
+    -   `True`: Use VIIRS GVF 
+-   `IFVIIRS_LAI [default: False]`
+    Logical to choose if the LAI is taken from FV3 or VIIRS.
+    -   `False`: Use FV3 LAI
+    -   `True`: Use VIIRS LAI
 -   `IFFENGSHA_DUST [default: False]`
     Logical to choose if the NOAA-ARL Fengsha Windblown dust is used in CMAQ:  based on input soil parameters in "GeoFile" and output in METCRO2D file.
     -   `False`: Do not input/output Fengsha WB Dust soil parameters
@@ -179,7 +183,7 @@ cd $NACC_home/serial/scripts
 |InMetFiles|netCDF (WRF or FV3-GFS)|List of WRF or FV3-GFS output files for input to NACC|required|
 |InSfcFiles|netCDF (FV3-GFS)|List of FV3-GFS output files for input to NACC|required (only FV3-GFS)|
 |InGeoFile|netCDF (WRFor FV3-GFS)|Output from WRF Geogrid processor|optional; only required if fractional land use, LAI, etc are not part of the WRF or FV3-GFS output.  Offline Pre-processed NOAA-ARL "geofiles" with LAI and LANDUSEF (based on 12-month climatological IGBP-MODIS) for the global GFSv16 Gaussian NetCDF Grid are available via FTP by request (Contact:  Patrick C. Campbell; Patrick.C.Campbell@noaa.gov)|
-|InVIIRSFile|netCDF (FV3-GFS)|Input from VIIRS GVF data |optional; only if global NetCDF VIIRS Input file is provided. Global ~ 4km VIIRS NetCDF Grid are available via FTP by request. (Contact:  Patrick C. Campbell; Patrick.C.Campbell@noaa.gov)|
+|InVIIRSFile|netCDF (FV3-GFS)|Input from VIIRS data |optional; only if global NetCDF VIIRS Input COARDS file is provided. Global ~ 4km VIIRS NetCDF Grid are available via FTP by request. (Contact:  Patrick C. Campbell; Patrick.C.Campbell@noaa.gov)|
 
 **Table 2. NACC output files**
 

parallel/scripts/run-nacc-fv3.ksh

@@ -24,9 +24,10 @@ export NODES=12
 APPL=aqm.t12z
 InMetDir=/gpfs/hps2/ptmp/Patrick.C.Campbell/NACC_FV3GFS16_runs/sens5_nacc_cmaq531_nofire_fengsha/com/aqm/prod/aqm.20190826
 InGeoDir=/gpfs/hps3/emc/naqfc/noscrub/Youhua.Tang/nwdev/NAQFC-WCOSS/fix
-InVIIRSDir=/gpfs/hps3/emc/naqfc/noscrub/Patrick.C.Campbell/viirs_gvf_test/grib2
-OutDir=/gpfs/hps2/ptmp/$USER/NACC-VIIRS-Test/output_fengsha_test
-ProgDir=/gpfs/hps3/emc/naqfc/noscrub/Patrick.C.Campbell/NACC/parallel/src
+InVIIRSDir_GVF=/gpfs/hps3/emc/naqfc/noscrub/Patrick.C.Campbell/viirs_gvf_test/grib2
+InVIIRSDir_LAI=/gpfs/hps3/emc/naqfc/noscrub/Patrick.C.Campbell/viirs_lai_test/
+OutDir=/gpfs/hps2/ptmp/$USER/NACC_test_lai/
+ProgDir=/gpfs/hps3/emc/naqfc/noscrub/Patrick.C.Campbell/NACC_viirs_lai/parallel/src
 
 if [ ! -s $InMetDir ]; then
   echo "No such input directory $InMetDir"
@@ -59,7 +60,8 @@ cat>namelist.mcip<<!
   file_mm    = '$InMetDir/gfs.t12z.atmf','.nc'
   file_sfc   = '$InMetDir/gfs.t12z.sfcf','.nc'
   file_geo   = '$InGeoDir/gfs.t12z.geo.08.nc'
-  file_viirs_gvf = '$InVIIRSDir/GVF-WKL-GLB_v2r3_j01_s20200824_e20200830_c202008311235100.grib2.nc'
+  file_viirs_gvf = '$InVIIRSDir_GVF/GVF-WKL-GLB_v2r3_j01_s20200824_e20200830_c202008311235100.grib2.nc'
+  file_viirs_lai = '$InVIIRSDir_LAI/VIIRS_VNP15A2H.001_20190829.nc'
   ioform     =  1
  &END
 
@@ -73,7 +75,8 @@ cat>namelist.mcip<<!
   lwout      =  1
   luvbout    =  1
   ifdiag_pbl = .FALSE.
-  ifviirs_gvf = .FALSE. 
+  ifviirs_gvf = .TRUE. 
+  ifviirs_lai = .TRUE.
   iffengsha_dust = .FALSE.
   mcip_start = "2019-08-26-12:00:00.0000"
   mcip_end   = "2019-08-29-13:00:00.0000"

parallel/src/alloc_met.f90

@@ -202,7 +202,7 @@ SUBROUTINE alloc_met
     !test
   ENDIF
 
-  IF ( iflai ) THEN  ! leaf area index available
+  IF ( iflai ) THEN  ! leaf area index available (Either from model or VIIRS)
     ALLOCATE ( lai    (ix, jy) )
   ENDIF
 
@@ -241,9 +241,9 @@ SUBROUTINE alloc_met
     ALLOCATE ( veg    (ix, jy) )
   ENDIF
 
-   IF ( ifveg_viirs ) THEN  ! VIIRS GVF, i.e.,vegetation fraction available
+  IF ( ifveg_viirs ) THEN  ! VIIRS GVF, i.e.,vegetation fraction available
     ALLOCATE ( veg    (ix, jy) )
-   ENDIF
+  ENDIF
 
   IF ( ifwr ) THEN  ! canopy wetness available
     ALLOCATE ( wr     (ix, jy) )

parallel/src/coord_mod.f90

@@ -189,8 +189,9 @@ MODULE coord
 !FV3GFS
   REAL, SAVE,    ALLOCATABLE   :: fv3lat(:),fv3lon(:), xindex(:,:), yindex(:,:), xuindex(:,:), yuindex(:,:), &
     xvindex(:,:),yvindex(:,:),xdindex(:,:),ydindex(:,:)
-!VIIRS (Optional for GVF--> veg frac currently)
-  REAL, SAVE,    ALLOCATABLE   :: viirslat(:),viirslon(:),xindex_viirs(:,:),yindex_viirs(:,:)
+!VIIRS (Optional for GVF--> veg frac or LAI currently)
+  REAL, SAVE,    ALLOCATABLE   :: viirslat_gvf(:),viirslon_gvf(:),xindex_viirs_gvf(:,:),yindex_viirs_gvf(:,:)
+  REAL, SAVE,    ALLOCATABLE   :: viirslat_lai(:),viirslon_lai(:),xindex_viirs_lai(:,:),yindex_viirs_lai(:,:)
 !-------------------------------------------------------------------------------
 ! COORDNAM_GD:
 ! The coordinate system name used for I/O-API description and GRIDDESC.

parallel/src/files_mod.f90

@@ -69,6 +69,7 @@ MODULE files
   CHARACTER(LEN=256), PARAMETER     :: file_nml   = 'namelist.mcip'
   CHARACTER(LEN=256)                :: file_geo
   CHARACTER(LEN=256)                :: file_viirs_gvf
+  CHARACTER(LEN=256)                :: file_viirs_lai
 
   CHARACTER(LEN=16),  PARAMETER     :: gridbdy2d  = 'GRID_BDY_2D     '
   CHARACTER(LEN=16),  PARAMETER     :: gridcro2d  = 'GRID_CRO_2D     '

parallel/src/mcipparm_mod.f90

@@ -265,12 +265,14 @@ MODULE mcipparm
   LOGICAL            :: iftkef          ! TKE (if exists) on full levels?
   LOGICAL            :: ifveg           ! vegetation fraction in input file?
   LOGICAL            :: ifveg_viirs     ! using vegetation fraction in VIIRS?
+  LOGICAL            :: iflai_viirs     ! using LAI in VIIRS?
   LOGICAL            :: ifw10m          ! 10-m wind components in input file?
   LOGICAL            :: ifwr            ! canopy wetness in input file?
   LOGICAL            :: ifznt           ! roughness length in input file?
   LOGICAL            :: needseaice      ! sea ice in input file?
   LOGICAL            :: ifdiag_pbl      ! use diagnostic PBL?
   LOGICAL            :: ifviirs_gvf     ! use NAQFC VIIRS GVF?
+  LOGICAL            :: ifviirs_lai     ! use NAQFC VIIRS LAI?
   LOGICAL            :: iffengsha_dust  ! use NAQFC Fengsha Windblown Dust?
 
 !-------------------------------------------------------------------------------

parallel/src/rdfv3.f90

@@ -197,6 +197,7 @@ SUBROUTINE rdfv3 (mcip_now,nn)
   INTEGER, SAVE                     :: cdfid, cdfid2
   INTEGER                           :: cdfidg
   INTEGER                           :: cdfid_vgvf
+  INTEGER                           :: cdfid_vlai
   INTEGER                           :: dimid
   INTEGER                           :: dimids     ( nf90_max_var_dims )
   REAL,    SAVE,      ALLOCATABLE   :: dum1d      ( : )
@@ -661,7 +662,7 @@ END SUBROUTINE windrotation
   if(.not.allocated(atmp)) allocate(atmp(ncols_x,nrows_x))
   if(.not.allocated(utmp)) allocate(utmp(ncols_x+1,nrows_x+1))
 
-   IF ( ifveg_viirs ) THEN !If using VIIRS GVF for vegetation fraction-->Need to allocate array for VIIRS
+   IF ( ifveg_viirs .OR. iflai_viirs ) THEN !If using VIIRS GVF or LAI ->Need to allocate array for VIIRS
     IF ( .NOT. ALLOCATED ( dum2d_viirs   ) )  &
     ALLOCATE ( dum2d_viirs   (met_nx_viirs, met_ny_viirs)      )        ! 2D
 
@@ -685,9 +686,14 @@ END SUBROUTINE windrotation
     yuindex(ncols_x+1,nrows_x+1),xvindex(ncols_x+1,nrows_x+1),yvindex(ncols_x+1,nrows_x+1), &
     xdindex(ncols_x+1,nrows_x+1),ydindex(ncols_x+1,nrows_x+1))
 
-   IF ( ifveg_viirs ) THEN !If using VIIRS GVF for vegetation fraction
-    allocate(xindex_viirs(ncols_x,nrows_x), yindex_viirs(ncols_x,nrows_x))
+   IF ( ifveg_viirs ) THEN !If using VIIRS GVF 
+    allocate(xindex_viirs_gvf(ncols_x,nrows_x), yindex_viirs_gvf(ncols_x,nrows_x))
    ENDIF
+
+    IF ( iflai_viirs) THEN !If using VIIRS LAI
+    allocate(xindex_viirs_lai(ncols_x,nrows_x), yindex_viirs_lai(ncols_x,nrows_x))
+   ENDIF
+
 
     ! Compute distance from origin (at reflat, standlon) to domain center, and
     ! store in MET_XXCTR and MET_YYCTR.  Then calculate latitude, longitude,
@@ -755,12 +761,34 @@ END SUBROUTINE windrotation
 
              mapcrs(i,j) = mapfac_lam (latcrs(i,j), met_tru1, met_tru2)
 
-             call getxyindex(latcrs(i,j),loncrs(i,j),xindex_viirs(i,j),yindex_viirs(i,j),viirslat,viirslon,met_nx_viirs,met_ny_viirs)
+             call getxyindex(latcrs(i,j),loncrs(i,j),xindex_viirs_gvf(i,j),yindex_viirs_gvf(i,j),viirslat_gvf,viirslon_gvf,met_nx_viirs,met_ny_viirs)
 
            ENDDO
          ENDDO
         ENDIF
 
+
+        IF ( iflai_viirs ) THEN !If using VIIRS LAI for vegetation fraction
+         xoff=-1.5
+         yoff=-1.5
+         DO j = 1, nrows_x
+           DO i = 1, ncols_x
+
+             xxin = met_xxctr + (FLOAT(i) + xoff) * met_resoln
+             yyin = met_yyctr + (FLOAT(j) + yoff) * met_resoln
+
+             CALL xy2ll_lam (xxin, yyin, met_tru1, met_tru2, met_proj_clon,  &
+                             met_ref_lat, latcrs(i,j), loncrs(i,j))
+
+             mapcrs(i,j) = mapfac_lam (latcrs(i,j), met_tru1, met_tru2)
+
+             call getxyindex(latcrs(i,j),loncrs(i,j),xindex_viirs_lai(i,j),yindex_viirs_lai(i,j),viirslat_lai,viirslon_lai,met_nx_viirs,met_ny_viirs)
+
+           ENDDO
+         ENDDO
+        ENDIF
+
+
         IF ( .NOT. gotfaces ) THEN  ! get lat, lon, map-scale factor on faces
 
           xoff = -2.0  ! U-face: no offset in X from dot-point center value
@@ -1701,8 +1729,8 @@ END SUBROUTINE windrotation
     ENDIF
    ENDIF !Fengsha WB dust variables
 
-    IF ( iflai ) THEN
-      IF ( iflaiwrfout ) THEN  ! leaf area index in FV3 history file
+  IF ( iflai ) THEN
+   IF ( iflaiwrfout ) THEN  ! leaf area index in FV3 history file
         CALL get_var_2d_real_cdf (cdfid2, 'LAI', dum2d, it, rcode)
         IF ( rcode == nf90_noerr ) THEN
            call myinterp(dum2d,met_nx,met_ny,atmp,xindex,yindex,ncols_x,nrows_x,1)
@@ -1717,7 +1745,7 @@ END SUBROUTINE windrotation
           WRITE (*,f9400) TRIM(pname), 'LAI', TRIM(nf90_strerror(rcode))
           CALL graceful_stop (pname)
         ENDIF
-      ELSE  ! leaf area index in GEOGRID file from WPS
+   ELSE  ! leaf area index in GEOGRID file from WPS
        flg = file_geo
         rcode = nf90_open (flg, nf90_nowrite, cdfidg)
         IF ( rcode /= nf90_noerr ) THEN
@@ -1743,8 +1771,33 @@ END SUBROUTINE windrotation
           WRITE (*,f9950) TRIM(pname)
           CALL graceful_stop (pname)
         ENDIF
+   ENDIF
+   IF ( iflai_viirs ) THEN !Instead Using VIIRS LAI
+    !Open VIIRS LAI File
+    flg = file_viirs_lai
+    rcode = nf90_open (flg, nf90_nowrite, cdfid_vlai)
+    IF ( rcode /= nf90_noerr ) THEN
+       WRITE (*,f9910) TRIM(pname)
+       CALL graceful_stop (pname)
+    ENDIF
+    CALL get_var_2d_real_cdf (cdfid_vlai, 'lai', dum2d_viirs, it, rcode)
+      IF ( rcode == nf90_noerr ) THEN
+         call myinterp(dum2d_viirs,met_nx_viirs,met_ny_viirs,atmp,xindex_viirs_lai,yindex_viirs_lai,ncols_x,nrows_x,1)
+        lai(1:ncols_x,1:nrows_x) = atmp(1:ncols_x,1:nrows_x)
+        where(lai.lt.0) lai=0.
+        WRITE (*,f6000) 'lai   ', lai(lprt_metx, lprt_mety), ' from VIIRS LAI'
+      ELSE
+        WRITE (*,f9440) TRIM(pname), 'lai', TRIM(nf90_strerror(rcode))
+        CALL graceful_stop (pname)
       ENDIF
+    rcode = nf90_close (cdfid_vlai)
+    IF ( rcode /= nf90_noerr ) THEN
+       WRITE (*,f9960) TRIM(pname)
+       CALL graceful_stop (pname)
     ENDIF
+   ENDIF
+
+  ENDIF
 
   IF ( ifwr ) THEN
     CALL get_var_2d_real_cdf (cdfid2, 'cnwat', dum2d, it, rcode)
@@ -1768,7 +1821,7 @@ END SUBROUTINE windrotation
     ENDIF
     CALL get_var_2d_real_cdf (cdfid_vgvf, 'VEG_surface', dum2d_viirs, it, rcode)
       IF ( rcode == nf90_noerr ) THEN
-         call myinterp(dum2d_viirs,met_nx_viirs,met_ny_viirs,atmp,xindex_viirs,yindex_viirs,ncols_x,nrows_x,1)
+         call myinterp(dum2d_viirs,met_nx_viirs,met_ny_viirs,atmp,xindex_viirs_gvf,yindex_viirs_gvf,ncols_x,nrows_x,1)
         veg(1:ncols_x,1:nrows_x) = atmp(1:ncols_x,1:nrows_x)*0.01
         where(veg.gt.100 .OR. veg.lt.0) veg=0.
         WRITE (*,f6000) 'veg   ', veg(lprt_metx, lprt_mety), 'fraction (from VIIRS GVF)'

parallel/src/readnml.f90

@@ -107,11 +107,11 @@ SUBROUTINE readnml (ctmlays)
   INTEGER                           :: nrowsin
   CHARACTER(LEN=16),  PARAMETER     :: pname      = 'READNML'
 
-  NAMELIST /filenames/   file_gd, file_mm, file_sfc, file_geo, file_viirs_gvf, ioform
+  NAMELIST /filenames/   file_gd, file_mm, file_sfc, file_geo, file_viirs_gvf, file_viirs_lai, ioform
 
   NAMELIST /userdefs/    inmetmodel, dx_in, dy_in, met_cen_lat_in, met_cen_lon_in, &
                          lpv, lwout, luvbout, ifdiag_pbl, ifviirs_gvf, &
-                         iffengsha_dust,    &
+                         ifviirs_lai, iffengsha_dust,    &
                          eradm, mcip_start, mcip_end, ntimes, intvl,  &
                          coordnam, grdnam, ctmlays,           &
                          btrim, lprt_col, lprt_row,           &
@@ -219,6 +219,7 @@ SUBROUTINE readnml (ctmlays)
   file_geo       = " "
   file_sfc(:)    = " "
   file_viirs_gvf = " "
+  file_viirs_lai = " "
 !-------------------------------------------------------------------------------
 ! Set default value for user-selected model (2 = WRF, 3 = FV3).
   inmetmodel = 2
@@ -254,6 +255,7 @@ SUBROUTINE readnml (ctmlays)
 
   ifdiag_pbl  = .false.  ! To re-diagnose the PBLH
   ifviirs_gvf = .false.  ! To use NAQFC VIIRS GVF input
+  ifviirs_lai = .false.  ! To use NAQFC VIIRS LAI input
   iffengsha_dust = .false.  ! To use NAQFC Fengsha Windblown Dust
 !-------------------------------------------------------------------------------
 ! Set default value for earth radius in meters (ERADM).  The default value is
@@ -368,6 +370,8 @@ SUBROUTINE readnml (ctmlays)
   file_viirs_gvf = TRIM( ADJUSTL(file_viirs_gvf) )
   IF ( file_viirs_gvf(1:7) == "no_file" ) file_viirs_gvf = " "
 
+  file_viirs_lai = TRIM( ADJUSTL(file_viirs_lai) )
+  IF ( file_viirs_lai(1:7) == "no_file" ) file_viirs_lai = " "
 !-------------------------------------------------------------------------------
 ! Verify values of user-defined options.
 !-------------------------------------------------------------------------------

parallel/src/setgriddefs.f90

@@ -762,6 +762,13 @@ SUBROUTINE setgriddefs
   ENDIF
   WRITE (*,f6150) 'VIIRS GVF-VEGETATION FRACTION', TRIM(yesno)
 
+  IF ( iflai_viirs ) THEN
+    yesno = ''
+  ELSE
+    yesno = 'NOT'
+  ENDIF
+  WRITE (*,f6150) 'VIIRS LAI', TRIM(yesno)
+
   IF ( ifw10m ) THEN
     yesno = ''
   ELSE