diff --git a/lis/metforcing/usaf/AGRMET_forcingMod.F90 b/lis/metforcing/usaf/AGRMET_forcingMod.F90 index 4a707ac24..220ca1563 100644 --- a/lis/metforcing/usaf/AGRMET_forcingMod.F90 +++ b/lis/metforcing/usaf/AGRMET_forcingMod.F90 @@ -302,7 +302,9 @@ module AGRMET_forcingMod ! contains the previous and next forcing values, respectively ! \end{description} ! -! !USES: +! !USES: + use LIS_constantsMod, only : LIS_CONST_PATH_LEN + implicit none PRIVATE @@ -719,6 +721,12 @@ module AGRMET_forcingMod ! EMK Add GFS-to-GALWEM bias correction integer :: back_bias_corr real, allocatable :: pcp_back_bias_ratio(:,:) + ! EMK Add NRT bias correction toward IMERG-Final Run + ! (back_bias_corr == 2) + character(LIS_CONST_PATH_LEN) :: gfs_nrt_bias_ratio_file + character(LIS_CONST_PATH_LEN) :: galwem_nrt_bias_ratio_file + real, allocatable :: gfs_nrt_bias_ratio(:,:) + real, allocatable :: galwem_nrt_bias_ratio(:,:) integer :: pcp_back_bias_ratio_month end type agrmet_type_dec @@ -747,7 +755,7 @@ module AGRMET_forcingMod subroutine init_AGRMET(findex) ! !USES: - use LIS_coreMod, only : LIS_rc, LIS_domain + use LIS_coreMod, only : LIS_rc use LIS_timeMgrMod, only : LIS_update_timestep use LIS_logMod, only : LIS_logunit, LIS_endrun use LIS_snowMod, only : LIS_snow_setup @@ -827,6 +835,21 @@ subroutine init_AGRMET(findex) integer :: istat1 integer :: ftn + external :: readcrd_agrmet + external :: AGRMET_readmask + external :: AGRMET_readterrain + external :: AGRMET_readpcpcntm + external :: AGRMET_read_sfcalccntm + external :: polarToLatLon + external :: bilinear_interp_input_withgrid + external :: neighbor_interp_input_withgrid + external :: AGRMET_read_pcpclimodata + external :: gfs_reset_interp_input + external :: galwem_reset_interp_input + external :: conserv_interp_input + external :: bilinear_interp_input + external :: upscaleByAveraging_input + allocate(agrmet_struc(LIS_rc%nnest)) call readcrd_agrmet() diff --git a/lis/metforcing/usaf/USAF_bratsethMod.F90 b/lis/metforcing/usaf/USAF_bratsethMod.F90 index b6dedf826..883e79acf 100644 --- a/lis/metforcing/usaf/USAF_bratsethMod.F90 +++ b/lis/metforcing/usaf/USAF_bratsethMod.F90 @@ -628,6 +628,8 @@ subroutine USAF_addSSMIObsData(this,imax,jmax,ra_tmp,nest) integer :: count_good_ssmi real, parameter :: FILL = -9999.0 + external :: polarToLatLon + net = "SSMI" platform = "SSMI" sigmaOSqr = agrmet_struc(nest)%bratseth_precip_ssmi_sigma_o_sqr @@ -774,6 +776,8 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex) character(len=10) :: yyyymmddhh integer :: c, r + external :: AGRMET_julhr_date10 + TRACE_ENTER("bratseth_getBackNWP") rc = 0 @@ -801,6 +805,8 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex) ! EMK...Fetch bias ratio if requested if (agrmet_struc(nest)%back_bias_corr .eq. 1) then call USAF_pcpBackBiasRatio_s2s(nest, yyyymmddhh) + else if (agrmet_struc(nest)%back_bias_corr .eq. 2) then + call USAF_pcpBackBiasRatio_nrt(nest, yyyymmddhh) end if write(LIS_logunit,*) & @@ -824,6 +830,20 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex) ! Use GALWEM Bratseth settings if we have the data if (rc .eq. 0) then pcp_src(k) = 'GALWEM' + + ! Apply bias correction to background field + if (agrmet_struc(nest)%back_bias_corr .eq. 2) then + write(LIS_logunit,*) & + '[INFO] Applying IMERG-based bias correction to GALWEM precip' + do r = 1, LIS_rc%gnr(nest) + do c = 1, LIS_rc%gnc(nest) + fg_data_glb(c,r) = & + fg_data_glb(c,r) * & + agrmet_struc(nest)%galwem_nrt_bias_ratio(c,r) + end do + end do + end if + end if endif @@ -851,9 +871,10 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex) end if ! Apply bias correction to background field - if (agrmet_struc(nest)%back_bias_corr .eq. 1) then + if (agrmet_struc(nest)%back_bias_corr .eq. 1 .and. & + rc .eq. 0) then write(LIS_logunit,*) & - '[INFO] Applying GALWEM-based bias correction to GFS precip' + '[INFO] Applying GALWEM-based bias correction to GFS precip' do r = 1, LIS_rc%gnr(nest) do c = 1, LIS_rc%gnc(nest) fg_data_glb(c,r) = & @@ -861,6 +882,19 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex) agrmet_struc(nest)%pcp_back_bias_ratio(c,r) end do end do + + else if (agrmet_struc(nest)%back_bias_corr .eq. 2 .and. & + rc .eq. 0) then + ! New IMERG option for NRT + write(LIS_logunit,*) & + '[INFO] Applying IMERG-based bias correction to GFS precip' + do r = 1, LIS_rc%gnr(nest) + do c = 1, LIS_rc%gnc(nest) + fg_data_glb(c,r) = & + fg_data_glb(c,r) * & + agrmet_struc(nest)%gfs_nrt_bias_ratio(c,r) + end do + end do end if ! Use GFS Bratseth settings if we have the data @@ -946,6 +980,8 @@ subroutine USAF_getSSMIObsData(nest,j6hr,use_twelve,precip3,precip6, & integer :: local_diff_grid integer :: ii,jj,kk,count_good_ssmi ! EMK TEST + external :: agrmet_ssmiprec_filename + TRACE_ENTER("bratseth_getSSMI") routine_name = 'USAF_getSSMIObsData' @@ -1135,6 +1171,9 @@ subroutine USAF_getGeoPrecipObsData(nest,j6hr,use_twelve,precip3,precip6,& logical, external :: is_geo_corrupted real, allocatable :: gest_temp(:,:,:) + external :: agrmet_geoprec_filename + external :: polarToLatLon + TRACE_ENTER("bratseth_getGeopPrcp") net = "GEOPRECIP" platform = "GEOPRECIP" @@ -1546,7 +1585,6 @@ subroutine USAF_analyzePrecip(precipAll,nest,back,hourindex,mrgp,precipOBA) use AGRMET_forcingMod, only: agrmet_struc use LIS_coreMod, only: LIS_rc, LIS_ews_halo_ind, LIS_ewe_halo_ind, & LIS_nss_halo_ind, LIS_nse_halo_ind, LIS_localPet - use LIS_logMod, only: LIS_logunit use USAF_OBAMod, only: OBA, createOBA,assignOBA ! Defaults @@ -1565,8 +1603,6 @@ subroutine USAF_analyzePrecip(precipAll,nest,back,hourindex,mrgp,precipOBA) real, allocatable :: invDataDensities(:) real, allocatable :: sumObsEstimates(:) integer :: npasses - integer :: r - character(len=10) :: new_name,type real :: convergeThresh real :: sigmaBSqr integer :: good_obs @@ -1688,6 +1724,9 @@ subroutine interpBack2ObsPts(this,nest,imax,jmax,back, & integer, external :: get_fieldpos real, parameter :: FILL = MISSING + external :: compute_grid_coord + external :: bilinear_interp + ! See if we have any observations available nobs = this%nobs if (nobs .eq. 0) return @@ -2765,6 +2804,11 @@ subroutine fldbld_precip_nwp(nest,findex,julhr,src,fc_hr, & integer :: rc2 logical :: first_time + external :: getAVNfilename + external :: AGRMET_getGALWEMfilename + external :: AGRMET_fg2lis_precip + external :: interp_galwem_first_guess + rc = 0 found = .false. @@ -2974,7 +3018,6 @@ subroutine check_grib_file(gribfile,yr1,mo1,da1,hr1,found, & #if (defined USE_GRIBAPI) use grib_api #endif - use LIS_coreMod, only: LIS_masterproc use LIS_logMod, only : LIS_logunit, LIS_abort, LIS_alert, & LIS_verify, LIS_endrun use LIS_mpiMod @@ -4093,7 +4136,9 @@ end function is_gauge logical function is_stn(net) implicit none character(len=10), intent(in) :: net + character(len=10) :: net_local logical :: answer + net_local = net answer = .true. is_stn = answer end function is_stn @@ -4367,8 +4412,7 @@ end subroutine destroy_hash2d subroutine build_hash2d(this,nest,imax,jmax,hash2d) ! Imports - use LIS_coreMod, only: LIS_rc, LIS_domain, LIS_localPet - use LIS_logMod, only: LIS_logunit + use LIS_coreMod, only: LIS_rc use LIS_mpiMod ! Defaults @@ -4382,26 +4426,21 @@ subroutine build_hash2d(this,nest,imax,jmax,hash2d) type(hash_list), allocatable, intent(out) :: hash2d(:,:) ! Local variables - real :: dlat, dlon, ctrlat, ctrlon integer :: nobs integer, allocatable :: cols(:), rows(:) - integer :: pet, pet_incr - logical :: found - integer :: j,r,c,gindex - integer :: ierr nobs = this%nobs if (nobs .eq. 0) return ! Get the columns and rows of the observations in the LIS grid - call find_LIS_cols_rows(this,nest,cols,rows) + call find_LIS_cols_rows(this, nest, cols, rows) ! Loop through the collected col, row information and add to the ! hash2d table. - call init_hash2d(LIS_rc%gnc(nest),LIS_rc%gnr(nest),hash2d) + call init_hash2d(LIS_rc%gnc(nest), LIS_rc%gnr(nest), hash2d) imax = LIS_rc%gnc(nest) ! Output jmax = LIS_rc%gnr(nest) ! Output - call insert_hash2d_array(imax,jmax,hash2d,nobs,cols,rows) + call insert_hash2d_array(imax, jmax, hash2d, nobs, cols, rows) ! do j = 1,nobs ! c = cols(j) ! if (c .eq. 0) cycle @@ -4910,12 +4949,10 @@ end subroutine find_LIS_cols_rows subroutine USAF_waterQC(this,nest,silent_rejects) ! Imports - use LIS_coreMod, only: LIS_rc use LIS_LMLCMod, only: LIS_LMLC use LIS_logMod, only: LIS_logunit use LIS_mpiMod - ! Defaults implicit none @@ -5440,9 +5477,9 @@ subroutine USAF_analyzeScreen(screenObs,nest,back,sigmaBSqr, & real, allocatable :: invDataDensities(:) real, allocatable :: sumObsEstimates(:) integer :: npasses - character(len=10) :: new_name,type + character(len=10) :: new_name, type real :: convergeThresh - integer :: i,j + integer :: j TRACE_ENTER("bratseth_analyzeScrn") ! Initialize analysis field with the background first guess. This will @@ -5614,6 +5651,8 @@ subroutine USAF_getCMORPHObsData(nest,j6hr,use_twelve, & integer :: ftn, ios integer :: i,j,k + external :: cmorfile_agrmet + TRACE_ENTER("bratseth_getCMORPH") net = "CMORPH" platform = "CMORPH" @@ -6397,14 +6436,13 @@ subroutine USAF_setBratsethScreenStats(src,n) end subroutine USAF_setBratsethScreenStats - ! Read bias ratio for background field + ! Read bias ratio for background field (S2S version, which uses + ! data in the LDT parameter file) subroutine USAF_pcpBackBiasRatio_s2s(n, yyyymmddhh) ! Imports use AGRMET_forcingMod, only: agrmet_struc - use LIS_coreMod, only: LIS_rc, LIS_localPet, & - LIS_ews_halo_ind, LIS_ewe_halo_ind, & - LIS_nss_halo_ind, LIS_nse_halo_ind + use LIS_coreMod, only: LIS_rc use LIS_logMod, only: LIS_logunit, LIS_endrun, LIS_verify #if ( defined USE_NETCDF3 || defined USE_NETCDF4 ) use netcdf @@ -6468,4 +6506,107 @@ subroutine USAF_pcpBackBiasRatio_s2s(n, yyyymmddhh) #endif end subroutine USAF_pcpBackBiasRatio_s2s -end module USAF_bratsethMod + + ! Read bias ratio for background field (NRT version, which uses + ! data in standalone netCDF file) + subroutine USAF_pcpBackBiasRatio_nrt(n, yyyymmddhh) + + ! Imports + use AGRMET_forcingMod, only: agrmet_struc + use LIS_coreMod, only: LIS_rc!, LIS_localPet, & + !LIS_ews_halo_ind, LIS_ewe_halo_ind, & + !LIS_nss_halo_ind, LIS_nse_halo_ind + + use LIS_logMod, only: LIS_logunit, LIS_endrun, LIS_verify +#if ( defined USE_NETCDF3 || defined USE_NETCDF4 ) + use netcdf +#endif + + ! Defaults + implicit none + + ! Arguments + integer, intent(in) :: n + character(len=10), intent(in) :: yyyymmddhh + + ! Locals + integer :: imonth + logical :: file_exists + integer :: ncid, ppt_ratio_Id + +#if ( defined USE_NETCDF3 || defined USE_NETCDF4 ) + + read(yyyymmddhh(5:6),'(i2.2)') imonth + + ! Accumulations ending at 00Z first of month are part of the + ! *previous* month. So decrement imonth by one in that situation, + ! and reset to 12 (December) at the year change. + if (agrmet_struc(n)%pcp_back_bias_ratio_month .ne. 0) then + if (yyyymmddhh(9:10) .eq. '00' .and. & + yyyymmddhh(7:8) .eq. '01') then + imonth = imonth - 1 + if (imonth .eq. 0) then + imonth = 12 + end if + end if + end if + if (imonth .eq. agrmet_struc(n)%pcp_back_bias_ratio_month .and. & + agrmet_struc(n)%pcp_back_bias_ratio_month .ne. 0) return + + agrmet_struc(n)%pcp_back_bias_ratio_month = imonth + + ! First, update the GFS ratios + inquire(file=agrmet_struc(n)%gfs_nrt_bias_ratio_file, & + exist=file_exists) + if (.not. file_exists) then + write(LIS_logunit,*) '[ERR] Cannot find ', & + trim(agrmet_struc(n)%gfs_nrt_bias_ratio_file) + call LIS_endrun() + end if + call LIS_verify( & + nf90_open(path=agrmet_struc(n)%gfs_nrt_bias_ratio_file, & + mode=NF90_NOWRITE, ncid=ncid), & + 'Error in nf90_open in USAF_pcpBackBiasRatio_nrt') + call LIS_verify(nf90_inq_varid(ncid, "biasRatio", ppt_ratio_Id), & + 'biasRatio field not found in bias file') + call LIS_verify(nf90_get_var(ncid, ppt_ratio_Id, & + agrmet_struc(n)%gfs_nrt_bias_ratio, & + start=(/1,1,imonth/), & + count=(/LIS_rc%gnc(n),LIS_rc%gnr(n),1/)), & + 'n90_get_var failed for biasRatio in LIS param file') + write(LIS_logunit,*) & + '[INFO] Read GFS bias ratios for month ', imonth, & + ' from ', trim(agrmet_struc(n)%gfs_nrt_bias_ratio_file) + call LIS_verify(nf90_close(ncid),& + 'nf90_close failed in USAF_pcpBackBiasRatio_nrt') + + ! Repeat for GALWEM ratios + inquire(file=agrmet_struc(n)%galwem_nrt_bias_ratio_file, & + exist=file_exists) + if (.not. file_exists) then + write(LIS_logunit,*) '[ERR] Cannot find ', & + trim(agrmet_struc(n)%galwem_nrt_bias_ratio_file) + call LIS_endrun() + end if + call LIS_verify( & + nf90_open(path=agrmet_struc(n)%galwem_nrt_bias_ratio_file, & + mode=NF90_NOWRITE, ncid=ncid), & + 'Error in nf90_open in USAF_pcpBackBiasRatio_nrt') + call LIS_verify(nf90_inq_varid(ncid, "biasRatio", ppt_ratio_Id), & + 'biasRatio field not found in bias file') + call LIS_verify(nf90_get_var(ncid, ppt_ratio_Id, & + agrmet_struc(n)%galwem_nrt_bias_ratio, & + start=(/1,1,imonth/), & + count=(/LIS_rc%gnc(n),LIS_rc%gnr(n),1/)),& + 'n90_get_var failed for biasRatio in LIS param file') + write(LIS_logunit,*) & + '[INFO] Read GALWEM bias ratios for month ', imonth, & + ' from ', trim(agrmet_struc(n)%galwem_nrt_bias_ratio_file) + call LIS_verify(nf90_close(ncid),& + 'nf90_close failed in USAF_pcpBackBiasRatio_nrt') + +#endif + + end subroutine USAF_pcpBackBiasRatio_nrt + + end module USAF_bratsethMod diff --git a/lis/metforcing/usaf/readcrd_agrmet.F90 b/lis/metforcing/usaf/readcrd_agrmet.F90 index 05452a3a4..47635b67f 100644 --- a/lis/metforcing/usaf/readcrd_agrmet.F90 +++ b/lis/metforcing/usaf/readcrd_agrmet.F90 @@ -40,7 +40,7 @@ subroutine readcrd_agrmet() use LIS_logMod, only : LIS_logunit, LIS_verify, LIS_abort, & LIS_endrun #if (defined SPMD) - use LIS_mpiMod, only: LIS_MPI_COMM + use LIS_mpiMod #endif use LIS_pluginIndices, only : LIS_agrmetrunId use AGRMET_forcingMod, only : agrmet_struc @@ -62,6 +62,7 @@ subroutine readcrd_agrmet() integer, external :: LIS_create_subdirs integer :: tmp_imerg_plp_thresh integer :: ierr + logical :: use_nrt_bias_files ! EMK call ESMF_ConfigFindLabel(LIS_config,"AGRMET forcing directory:",rc=rc) do n=1,LIS_rc%nnest @@ -1091,16 +1092,59 @@ subroutine readcrd_agrmet() call LIS_verify(rc, & "[ERR] AGRMET PPT Background bias correction option: not specified in config file") if (agrmet_struc(n)%back_bias_corr .lt. 0 .or. & - agrmet_struc(n)%back_bias_corr .gt. 1) then + agrmet_struc(n)%back_bias_corr .gt. 2) then call LIS_verify(rc, & - "[ERR] AGRMET PPT Background bias correction option: bad value in config file, set 0 or 1") + "[ERR] AGRMET PPT Background bias correction option: bad value in config file, set 0, 1, or 2") end if if (agrmet_struc(n)%back_bias_corr .eq. 1) then allocate(agrmet_struc(n)%pcp_back_bias_ratio(LIS_rc%gnc(n),LIS_rc%gnr(n))) agrmet_struc(n)%pcp_back_bias_ratio = 1. agrmet_struc(n)%pcp_back_bias_ratio_month = 0 + else if (agrmet_struc(n)%back_bias_corr .eq. 2) then + + allocate(agrmet_struc(n)%gfs_nrt_bias_ratio( & + LIS_rc%gnc(n),LIS_rc%gnr(n))) + allocate(agrmet_struc(n)%galwem_nrt_bias_ratio( & + LIS_rc%gnc(n),LIS_rc%gnr(n))) + agrmet_struc(n)%gfs_nrt_bias_ratio = 1. + agrmet_struc(n)%galwem_nrt_bias_ratio = 1. + agrmet_struc(n)%pcp_back_bias_ratio_month = 0 end if - enddo ! n + end do + + ! EMK Add support for NRT bias files for GFS and GALWEM + use_nrt_bias_files = .false. + do n = 1, LIS_rc%nnest + if (agrmet_struc(n)%back_bias_corr .eq. 2) then + use_nrt_bias_files = .true. + exit + end if + end do + + if (use_nrt_bias_files) then + call ESMF_ConfigFindLabel(LIS_config, & + "AGRMET PPT GFS NRT bias file:", rc=rc) + call LIS_verify(rc, & + "[ERR] AGRMET PPT GFS NRT bias file: not specified in config file") + + do n = 1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config, & + agrmet_struc(n)%gfs_nrt_bias_ratio_file, rc=rc) + call LIS_verify(rc, & + "[ERR] AGRMET PPT GFS NRT bias file: not specified in config file") + end do + + call ESMF_ConfigFindLabel(LIS_config, & + "AGRMET PPT GALWEM NRT bias file:", rc=rc) + call LIS_verify(rc, & + "[ERR] AGRMET PPT GALWEM NRT bias file: not specified in config file") + do n = 1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config, & + agrmet_struc(n)%galwem_nrt_bias_ratio_file, rc=rc) + call LIS_verify(rc, & + "[ERR] AGRMET PPT GALWEM NRT bias file: not specified in config file") + end do + end if do n=1,LIS_rc%nnest agrmet_struc(n)%radProcessInterval = 1 diff --git a/lvt/configs/lvt.config.adoc b/lvt/configs/lvt.config.adoc index 9ba94fea5..0c3f15a93 100644 --- a/lvt/configs/lvt.config.adoc +++ b/lvt/configs/lvt.config.adoc @@ -118,6 +118,7 @@ LVT output methodology: "2d gridspace" | "`ANSA SWE`" | ANSA SWE retrievals | "`CPC precipitation`" | CPC unified precipitation product | "`GPM IMERG`" | IMERG precipitation data product +| "`GPM IMERG Monthly`" | IMERG-Monthly precipitation product | "`USGS streamflow`" | USGS streamflow observations | "`Naturalized streamflow`" | Naturalized streamflow estimates | "`FLUXNET MTE`" | Gridded FLUXNET MTE data from MPI @@ -1276,6 +1277,21 @@ IMERG version: V06B IMERG product: final .... +==== GPM IMERG Monthly precipitation data + +`IMERG monthly data directory:` specifies the location of the GPM IMERG Monthly precipitation product + +`IMERG monthly version:` specifies the version of the GPM IMERG Monthly precipitation product. Most current version is V07A. + +`IMERG monthly product:` specifies the GPM IMERG Monthly precipitation product. Current option is: final. (early and late will be added later) + +.Example _lvt.config_ entry +.... +IMERG monthly data directory: input/IMERGF_V07A_Monthly +IMERG monthly version: V07A +IMERG monthly product: final +.... + ==== ESA CCI soil moisture data `ESA CCI soil moisture data directory:` specifies the location of the ESA CCI soil moisture data diff --git a/lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90 b/lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90 new file mode 100644 index 000000000..9a98a1c91 --- /dev/null +++ b/lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90 @@ -0,0 +1,170 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.4 +! +! Copyright (c) 2022 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +! +! NOTE: Currently only V07A IMERG Final Run data are supported +module IMERG_monthly_dataMod + + ! Defaults + implicit none + private + + ! Public routines + public :: IMERG_monthly_datainit + + ! Public types + public :: imergmonthlydata + + type, public :: imergmonthlydatadec + character*100 :: odir + character*10 :: imergver, imergprd + real :: datares + real, allocatable :: rlat(:) + real, allocatable :: rlon(:) + + ! This is only used with upscale averaging + integer, allocatable :: n11(:) + + ! These are only used with budget interpolation + integer, allocatable :: n112(:,:) + integer, allocatable :: n122(:,:) + integer, allocatable :: n212(:,:) + integer, allocatable :: n222(:,:) + real, allocatable :: w112(:,:) + real, allocatable :: w122(:,:) + real, allocatable :: w212(:,:) + real, allocatable :: w222(:,:) + + integer :: nc + integer :: nr + end type imergmonthlydatadec + + type(imergmonthlydatadec), allocatable :: imergmonthlydata(:) + +contains + + subroutine IMERG_monthly_datainit(i) + + ! Imports + use ESMF + use LVT_coreMod, only: LVT_rc, LVT_config, LVT_isAtAFinerResolution + use LVT_logMod, only: LVT_logunit, LVT_verify, LVT_endrun + + ! Defaults + implicit none + + ! Arguments + integer, intent(in) :: i + + ! Local variables + integer :: status + real :: gridDesci(50) + + ! External routines + external :: conserv_interp_input + external :: upscaleByAveraging_input + + if (.not. allocated(imergmonthlydata)) then + allocate(imergmonthlydata(LVT_rc%nDataStreams)) + endif + + ! Get top level IMERG data directory + call ESMF_ConfigGetAttribute(LVT_Config, imergmonthlydata(i)%odir, & + label='IMERG monthly data directory:', rc=status) + call LVT_verify(status, 'IMERG monthly data directory: not defined') + + ! Get IMERG product + call ESMF_ConfigFindLabel(LVT_config,"IMERG monthly product:", & + rc=status) + call LVT_verify(status, 'IMERG monthly product: not defined') + call ESMF_ConfigGetAttribute(LVT_config, & + imergmonthlydata(i)%imergprd, rc=status) + if (imergmonthlydata(i)%imergprd .ne. 'final') then + write(LVT_logunit,*)'[ERR] Invalid IMERG monthly product specified' + write(LVT_logunit,*)"[ERR] Currently only 'final' is supported" + call LVT_endrun() + end if + + ! Get IMERG version + call ESMF_ConfigFindLabel(LVT_config, "IMERG monthly version:", & + rc=status) + call LVT_verify(status, 'IMERG monthly version: not defined') + call ESMF_ConfigGetAttribute(LVT_config, & + imergmonthlydata(i)%imergver, rc=status) + if (imergmonthlydata(i)%imergver .ne. 'V07A') then + write(LVT_logunit,*)'[ERR] Invalid IMERG monthly version specified' + write(LVT_logunit,*)"[ERR] Currently only 'V07A' is supported" + call LVT_endrun() + end if + + ! Allocate arrays on LVT grid + allocate(imergmonthlydata(i)%rlat(LVT_rc%lnc*LVT_rc%lnr)) + allocate(imergmonthlydata(i)%rlon(LVT_rc%lnc*LVT_rc%lnr)) + + ! Set IMERG grid and map projection information. + gridDesci(:) = 0 + gridDesci(1) = 0 ! Lat/lon + gridDesci(2) = 3600 ! Points along latitude circle + gridDesci(3) = 1800 ! Points along longitude circle + gridDesci(4) = -89.95 ! Latitude of first grid point + gridDesci(5) = -179.95 ! Longitude of first grid point + gridDesci(6) = 128 + gridDesci(7) = 89.95 ! Latitude of last grid point + gridDesci(8) = 179.95 ! Longitude of last grid point + gridDesci(9) = 0.1 ! Longitudinal direction increment + gridDesci(10) = 0.1 ! Latitude direction increment + gridDesci(20) = 64 + + ! Set up interpolation data + imergmonthlydata(i)%datares = 0.1 + imergmonthlydata(i)%nc = 3600 + imergmonthlydata(i)%nr = 1800 + + ! Use budget-bilinear interpolation if IMERG resolution (0.1 deg) + ! is coarser than the analysis grid. Use upscale averaging if + ! IMERG is finer than the analysis grid. + if (LVT_isAtAFinerResolution(imergmonthlydata(i)%datares)) then + + ! Used only with budget interpolation + allocate(imergmonthlydata(i)%n112(LVT_rc%lnc*LVT_rc%lnr,25)) + allocate(imergmonthlydata(i)%n122(LVT_rc%lnc*LVT_rc%lnr,25)) + allocate(imergmonthlydata(i)%n212(LVT_rc%lnc*LVT_rc%lnr,25)) + allocate(imergmonthlydata(i)%n222(LVT_rc%lnc*LVT_rc%lnr,25)) + allocate(imergmonthlydata(i)%w112(LVT_rc%lnc*LVT_rc%lnr,25)) + allocate(imergmonthlydata(i)%w122(LVT_rc%lnc*LVT_rc%lnr,25)) + allocate(imergmonthlydata(i)%w212(LVT_rc%lnc*LVT_rc%lnr,25)) + allocate(imergmonthlydata(i)%w222(LVT_rc%lnc*LVT_rc%lnr,25)) + imergmonthlydata(i)%n112 = 0 + imergmonthlydata(i)%n122 = 0 + imergmonthlydata(i)%n212 = 0 + imergmonthlydata(i)%n222 = 0 + imergmonthlydata(i)%w112 = 0 + imergmonthlydata(i)%w122 = 0 + imergmonthlydata(i)%w212 = 0 + imergmonthlydata(i)%w222 = 0 + call conserv_interp_input(gridDesci, LVT_rc%gridDesc, & + LVT_rc%lnc*LVT_rc%lnr, & + imergmonthlydata(i)%rlat, imergmonthlydata(i)%rlon, & + imergmonthlydata(i)%n112, imergmonthlydata(i)%n122, & + imergmonthlydata(i)%n212, imergmonthlydata(i)%n222, & + imergmonthlydata(i)%w112, imergmonthlydata(i)%w122, & + imergmonthlydata(i)%w212, imergmonthlydata(i)%w222) + else + + ! Used only with upscale averaging + allocate(imergmonthlydata(i)%n11(imergmonthlydata(i)%nc * & + imergmonthlydata(i)%nr)) + imergmonthlydata(i)%n11 = 0 + call upscaleByAveraging_input(gridDesci, LVT_rc%gridDesc,& + imergmonthlydata(i)%nc*imergmonthlydata(i)%nr, & + LVT_rc%lnc*LVT_rc%lnr, & + imergmonthlydata(i)%n11) + end if + end subroutine IMERG_monthly_datainit +end module IMERG_monthly_dataMod diff --git a/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90 b/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90 new file mode 100644 index 000000000..b0da4236e --- /dev/null +++ b/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90 @@ -0,0 +1,332 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.4 +! +! Copyright (c) 2022 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- + +#include "LVT_misc.h" + +!NOTE: Currently only V07A IMERG Final Run monthly data are supported. +subroutine readIMERGmonthlydata(source) + + ! Imports + use ESMF + use IMERG_monthly_dataMod, only: imergmonthlydata + use LVT_coreMod, only: LVT_rc, LVT_isAtAFinerResolution + use LVT_histDataMod, only: LVT_logSingleDataStreamVar, & + LVT_MOC_totalprecip + use LVT_logMod, only: LVT_verify, LVT_logunit + use LVT_timeMgrMod, only: LVT_calendar + + ! Defaults + implicit none + + ! Arguments + integer, intent(in) :: source + + ! Local variables + character*255 :: filename + logical :: file_exists + real :: prcp_in(imergmonthlydata(source)%nc, & + imergmonthlydata(source)%nr) + real :: prcp_in1(imergmonthlydata(source)%nc * & + imergmonthlydata(source)%nr) + logical*1 :: lb(imergmonthlydata(source)%nc * & + imergmonthlydata(source)%nr) + logical*1 :: lo(LVT_rc%lnc*LVT_rc%lnr) + real :: prcp(LVT_rc%lnc*LVT_rc%lnr) + real :: prcp_final(LVT_rc%lnc, LVT_rc%lnr) + integer :: t, c, r + integer :: iret, ireaderr + real :: currTime + logical :: alarmCheck + integer :: yr1, mo1, da1, hr1, mn1, ss1 + integer :: yr2, mo2, da2, hr2, mn2, ss2 + type(ESMF_Time) :: time1 + type(ESMF_Time) :: time2 + type(ESMF_TimeInterval) :: lis_ts + integer :: status + integer :: days_in_month + + ! External routines + external :: create_IMERG_monthly_filename + external :: read_imergmonthly_hdf + external :: conserv_interp + external :: upscaleByAveraging + + ! Initialize variables + prcp = LVT_rc%udef + prcp_final = LVT_rc%udef + currTime = float(LVT_rc%dhr(source))*3600 + & + 60*LVT_rc%dmn(source) + LVT_rc%dss(source) + + ! Read IMERG at beginning of the month + alarmCheck = .false. + if (LVT_rc%dda(source) == 1 .and. & + LVT_rc%dhr(source) == 0 .and. & + LVT_rc%dmn(source) == 0 .and. & + LVT_rc%dss(source) == 0) alarmCheck = .true. + + yr1 = LVT_rc%dyr(source) + mo1 = LVT_rc%dmo(source) + da1 = LVT_rc%dda(source) + hr1 = LVT_rc%dhr(source) + mn1 = LVT_rc%dmn(source) + ss1 = 0 + call ESMF_TimeSet(time1, yy=yr1, mm=mo1, dd=da1, & + h=hr1, m=mn1, s=ss1, calendar=LVT_calendar, rc=status) + call LVT_verify(status) + + ! Use previous month file. + call ESMF_TimeIntervalSet(lis_ts, mm=1, calendar=LVT_calendar, & + rc=status) + call LVT_verify(status) + time2 = time1 - lis_ts + call ESMF_TimeGet(time2, yy=yr2, mm=mo2, dd=da2, & + h=hr2, m=mn2, s=ss2, calendar=LVT_calendar, & + rc=status) + call LVT_verify(status) + + if (alarmCheck) then + call create_IMERG_monthly_filename(source, yr2, mo2, filename) + inquire(file=trim(filename), exist=file_exists) + + if (file_exists) then + write(LVT_logunit,*) '[INFO] Reading IMERG data ',trim(filename) + call read_imergmonthly_hdf(filename, & + imergmonthlydata(source)%nc, & + imergmonthlydata(source)%nr, prcp_in, ireaderr) + if (ireaderr .eq. 0) then + + ! Use budget-bilinear interpolation if IMERG data are at + ! coarser resolution than the analysis grid; otherwise, use + ! upscale averaging. + prcp_in1 = LVT_rc%udef + lb = .false. + t = 1 + do r = 1, imergmonthlydata(source)%nr + do c = 1,imergmonthlydata(source)%nc + prcp_in1(t) = prcp_in(c,r) + if (prcp_in1(t) .ge. 0) then + lb(t) = .true. + end if + t = t + 1 + end do ! c + end do ! r + + if (LVT_isAtAFinerResolution( & + imergmonthlydata(source)%datares)) then + call conserv_interp(LVT_rc%gridDesc, lb, prcp_in1, & + lo, prcp, & + (imergmonthlydata(source)%nc * & + imergmonthlydata(source)%nr), & + (LVT_rc%lnc*LVT_rc%lnr), & + imergmonthlydata(source)%rlat, & + imergmonthlydata(source)%rlon, & + imergmonthlydata(source)%w112, & + imergmonthlydata(source)%w122, & + imergmonthlydata(source)%w212, & + imergmonthlydata(source)%w222, & + imergmonthlydata(source)%n112, & + imergmonthlydata(source)%n122, & + imergmonthlydata(source)%n212, & + imergmonthlydata(source)%n222, & + LVT_rc%udef, iret) + else + call upscaleByAveraging( & + (imergmonthlydata(source)%nc * & + imergmonthlydata(source)%nr), & + (LVT_rc%lnc*LVT_rc%lnr), LVT_rc%udef, & + imergmonthlydata(source)%n11, lb, & + prcp_in1, lo, prcp) + endif + write(LVT_logunit,*) '[INFO] Finished processing ', & + trim(filename) + else + write(LVT_logunit,*) '[ERR] Read error with IMERG file ', & + trim(filename) + prcp = LVT_rc%udef + endif + else + write(LVT_logunit,*) '[ERR] Missing IMERG file ', trim(filename) + prcp = LVT_rc%udef + end if ! file_exists + + do r = 1, LVT_rc%lnr + do c = 1, LVT_rc%lnc + prcp_final(c,r) = prcp(c+(r-1)*LVT_rc%lnc) + end do ! c + end do ! r + end if ! alarmCheck + + ! Convert mm/hr to kg/m2s. + do r = 1, LVT_rc%lnr + do c = 1, LVT_rc%lnc + if (prcp_final(c,r) .ge. 0) then + prcp_final(c,r) = prcp_final(c,r) / 3600. + else + prcp_final(c,r) = LVT_rc%udef + end if + end do ! c + end do ! r + call LVT_logSingleDataStreamVar(LVT_MOC_totalprecip, source, & + prcp_final, vlevel=1, units='kg/m2s') + + ! Now convert from kg/m2s to kg/m2. Need to calculate length of + ! month in seconds. + lis_ts = time1 - time2 + call ESMF_TimeIntervalGet(lis_ts, d=days_in_month, rc=status) + call LVT_verify(status) + do r = 1, LVT_rc%lnr + do c = 1, LVT_rc%lnc + if (prcp_final(c,r) >= 0) then + prcp_final(c,r) = & + prcp_final(c,r) * days_in_month * 86400 ! kg/m2 + else + prcp_final(c,r) = LVT_rc%udef + endif + enddo ! c + enddo ! r + call LVT_logSingleDataStreamVar(LVT_MOC_totalprecip, source, & + prcp_final, vlevel=1, units='kg/m2') + +end subroutine readIMERGmonthlydata + +subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr) + + ! Imports + use LVT_coreMod, only: LVT_rc + use LVT_logMod, only: LVT_logunit + +#if (defined USE_HDF5) + use HDF5 +#endif + + ! Defaults + implicit none + + ! Arguments + character(len=*), intent(in) :: filename + integer, intent(in) :: col, row + integer, intent(out) :: ireaderr + real, intent(out) :: precipout(col,row) + + ! Local variables + integer :: xsize, ysize + character(len=40) :: dsetname = '/Grid/precipitation' + real :: precipin(row,col) + integer :: istatus + integer :: i, j +#if (defined USE_HDF5) + integer(HSIZE_T), dimension(2) :: dims + integer(HID_T) :: fileid, dsetid +#endif + + precipin = LVT_rc%udef + precipout = LVT_rc%udef + +#if (defined USE_HDF5) + xsize = col + ysize = row + dims(1) = xsize + dims(2) = ysize + + ireaderr = 0 + + ! Open Fortran interface + call h5open_f(istatus) + if (istatus .ne. 0) then + write(LVT_logunit,*) 'Error opening HDF5 fortran interface' + ireaderr = istatus + return + end if + + ! Open HDF5 file + call h5fopen_f(filename, H5F_ACC_RDONLY_F, fileid, istatus) + if (istatus .ne. 0) then + write(LVT_logunit,*) 'Error opening IMERG file', trim(filename) + ireaderr = istatus + call h5close_f(istatus) ! Close HDF5 interface + return + end if + + ! Open precip dataset + call h5dopen_f(fileid, dsetname, dsetid, istatus) + if (istatus .ne. 0) then + write(LVT_logunit,*) 'Error opening IMERG dataset', trim(dsetname) + ireaderr = istatus + call h5fclose_f(fileid, istatus) ! Close HDF5 file + call h5close_f(istatus) ! Close HDF5 interface + return + end if + + ! Read dataset + call h5dread_f(dsetid, H5T_NATIVE_REAL, precipin, dims, istatus) + if (istatus .ne. 0) then + write(LVT_logunit,*) 'Error reading IMERG dataset', trim(dsetname) + ireaderr = istatus + call h5dclose_f(dsetid, istatus) ! Close dataset + call h5fclose_f(fileid, istatus) ! Close HDF5 file + call h5close_f(istatus) ! Close HDF5 interface + return + end if + + ! Put the real(1:,1:) on the precipout(0:,0:) + ! precipin is (ysize,xsize) starting at (lon=-179.9,lat=-89.9) + precipout(1:xsize,1:ysize) = transpose(precipin) + + ! Clean up. Since the data have already been copied, we will assume + ! any HDF5 errors from this point will have no impact on the data, and + ! therefore we will ignore the status codes. + call h5dclose_f(dsetid, istatus) ! Close HDF5 dataset + call h5fclose_f(fileid, istatus) ! Close HDF5 file + call h5close_f(istatus) ! Close HDF5 interface + +#endif + +end subroutine read_imergmonthly_hdf + +subroutine create_IMERG_monthly_filename(source, yr, mo, filename) + + ! Imports + use IMERG_monthly_dataMod, only: imergmonthlydata + use LVT_logMod, only: LVT_logunit, LVT_endrun + + ! Defaults + implicit none + + ! Arguments + integer, intent(in) :: source, yr, mo + character(len=*), intent(out) :: filename + + ! Local variables + character*4 :: cyr + character*2 :: cmo + character*100 :: fstem, fext + character*255 :: odir + character*4 :: imVer + + write(cyr, '(I4.4)') yr + write(cmo, '(I2.2)') mo + odir = trim(imergmonthlydata(source)%odir) + imVer = imergmonthlydata(source)%imergver + + ! FIXME: Add support for Early and Late Runs + if (imergmonthlydata(source)%imergprd == 'final') then + fstem = '/3B-MO.MS.MRG.3IMERG.' + fext = '.HDF5' + else + write(LVT_logunit,*) "[ERR] Invalid IMERG product option was chosen." + write(LVT_logunit,*) "[ERR] Please choose either 'final'." + call LVT_endrun() + endif + + filename = trim(odir) // "/" // cyr // trim(fstem) // & + cyr // cmo // "01-S000000-E235959." // cmo // "." // & + trim(imVer) // fext + +end subroutine create_IMERG_monthly_filename diff --git a/lvt/make/Filepath b/lvt/make/Filepath index d0034d33b..a1e42cc42 100644 --- a/lvt/make/Filepath +++ b/lvt/make/Filepath @@ -1 +1 @@ -dirs := . ../main ../core ../metrics ../plugins ../interp ../lib/bil ../domains/UTM ../domains/lambert ../domains/latlon ../datastreams/template ../datastreams/LISout ../datastreams/LISDAobs ../datastreams/ISCCP_Tskin ../datastreams/CEOP ../datastreams/SCAN ../datastreams/SCANGMAO ../datastreams/NASMD ../datastreams/GHCN ../datastreams/SNOTEL ../datastreams/SURFRAD ../datastreams/WGPBMR ../datastreams/LSWG_Tb ../datastreams/FMI_SWE ../datastreams/FMI_SNWD ../datastreams/CMC_SNWD ../datastreams/SNODAS ../datastreams/NASA_AMSRE_sm ../datastreams/LPRM_AMSRE_sm ../datastreams/AMMA ../datastreams/Ameriflux ../datastreams/ARM ../datastreams/SMOSREX ../datastreams/AGRMET ../datastreams/GlobSnow ../datastreams/SNODEPmetobs ../datastreams/SNODEP ../datastreams/MOD10A1 ../datastreams/MODSCAG ../datastreams/MOD10A1V6 ../datastreams/ANSA_SNWD ../datastreams/ANSA_SWE ../datastreams/CPC_PRCP ../datastreams/USGS_streamflow_gridded ../datastreams/USGS_streamflow ../datastreams/Natural_streamflow ../datastreams/ISMN ../datastreams/FLUXNETmte ../datastreams/FLUXNET2015 ../datastreams/FLUXNET2015_NC ../datastreams/MOD16A2 ../datastreams/UWET ../datastreams/ARSsm ../datastreams/NLDAS2 ../datastreams/ALEXI ../datastreams/ALEXIesi ../datastreams/GRACE ../datastreams/simGRACE ../datastreams/USGS_GWwell ../datastreams/PBOH2O ../datastreams/SMOS_L2sm ../datastreams/SMOS_NESDIS ../datastreams/SMOS_L1TB ../datastreams/SMOS_CATDS_L3sm ../datastreams/GCOMW_AMSR2L3sm ../datastreams/GCOMW_AMSR2L3snd ../datastreams/SMOPS ../datastreams/MODIS_LST ../datastreams/GLERL_HydroData ../datastreams/JULESdata ../datastreams/JULES2Ddata ../datastreams/ESACCI_sm ../datastreams/GIMMS_AVHRR_NDVI ../datastreams/GIMMS_MODIS_NDVI ../datastreams/GLDAS1 ../datastreams/GLDAS2 ../datastreams/MERRA2 ../datastreams/MERRA2asm ../datastreams/MERRA-Land ../datastreams/SSEBop ../datastreams/GRDC ../datastreams/GOES_LST ../datastreams/ERAinterimLand ../datastreams/SMAPsm ../datastreams/SMAPvod ../datastreams/LPRMvod ../datastreams/SMAPvwc ../datastreams/SMAP_L3TB ../datastreams/SMAPTB ../datastreams/GOME2_SIF ../datastreams/LVTbenchmarkOUT ../datastreams/LIS6out ../datastreams/LISDAdiag ../datastreams/Daymet ../datastreams/CMORPH ../datastreams/CHIRPSv2 ../datastreams/3B42V7 ../datastreams/USCRNsm ../datastreams/GLEAM ../datastreams/USDM ../datastreams/LVTpercentile ../datastreams/IMD_PRCP ../datastreams/APHRO_PRCP ../datastreams/GLASSlai ../datastreams/MODISsportLAI ../datastreams/FLUXCOM ../datastreams/HAR ../datastreams/OCO2_SIF ../datastreams/ECMWFforc ../datastreams/GDASforc ../datastreams/ASO_SWE ../runmodes/DataComp ../runmodes/557post ../runmodes/LISpost ../runmodes/DAstats ../runmodes/optUE ../runmodes/Benchmarking ../runmodes/USAFSIpost ../training/LinearRegression/ ../datastreams/GLASSalbedo ../datastreams/IMERG ../datastreams/UA_SNOW ../datastreams/OzFlux ../datastreams/jasmin ../datastreams/MCD15A2H ../datastreams/ERA5 ../datastreams/FluxSat_GPP ../datastreams/THySM/ ../datastreams/GRUNrunoff ../datastreams/UASMAP/ ../datastreams/COAMPSout/ ../datastreams/SMAP_E_OPL/ +dirs := . ../main ../core ../metrics ../plugins ../interp ../lib/bil ../domains/UTM ../domains/lambert ../domains/latlon ../datastreams/template ../datastreams/LISout ../datastreams/LISDAobs ../datastreams/ISCCP_Tskin ../datastreams/CEOP ../datastreams/SCAN ../datastreams/SCANGMAO ../datastreams/NASMD ../datastreams/GHCN ../datastreams/SNOTEL ../datastreams/SURFRAD ../datastreams/WGPBMR ../datastreams/LSWG_Tb ../datastreams/FMI_SWE ../datastreams/FMI_SNWD ../datastreams/CMC_SNWD ../datastreams/SNODAS ../datastreams/NASA_AMSRE_sm ../datastreams/LPRM_AMSRE_sm ../datastreams/AMMA ../datastreams/Ameriflux ../datastreams/ARM ../datastreams/SMOSREX ../datastreams/AGRMET ../datastreams/GlobSnow ../datastreams/SNODEPmetobs ../datastreams/SNODEP ../datastreams/MOD10A1 ../datastreams/MODSCAG ../datastreams/MOD10A1V6 ../datastreams/ANSA_SNWD ../datastreams/ANSA_SWE ../datastreams/CPC_PRCP ../datastreams/USGS_streamflow_gridded ../datastreams/USGS_streamflow ../datastreams/Natural_streamflow ../datastreams/ISMN ../datastreams/FLUXNETmte ../datastreams/FLUXNET2015 ../datastreams/FLUXNET2015_NC ../datastreams/MOD16A2 ../datastreams/UWET ../datastreams/ARSsm ../datastreams/NLDAS2 ../datastreams/ALEXI ../datastreams/ALEXIesi ../datastreams/GRACE ../datastreams/simGRACE ../datastreams/USGS_GWwell ../datastreams/PBOH2O ../datastreams/SMOS_L2sm ../datastreams/SMOS_NESDIS ../datastreams/SMOS_L1TB ../datastreams/SMOS_CATDS_L3sm ../datastreams/GCOMW_AMSR2L3sm ../datastreams/GCOMW_AMSR2L3snd ../datastreams/SMOPS ../datastreams/MODIS_LST ../datastreams/GLERL_HydroData ../datastreams/JULESdata ../datastreams/JULES2Ddata ../datastreams/ESACCI_sm ../datastreams/GIMMS_AVHRR_NDVI ../datastreams/GIMMS_MODIS_NDVI ../datastreams/GLDAS1 ../datastreams/GLDAS2 ../datastreams/MERRA2 ../datastreams/MERRA2asm ../datastreams/MERRA-Land ../datastreams/SSEBop ../datastreams/GRDC ../datastreams/GOES_LST ../datastreams/ERAinterimLand ../datastreams/SMAPsm ../datastreams/SMAPvod ../datastreams/LPRMvod ../datastreams/SMAPvwc ../datastreams/SMAP_L3TB ../datastreams/SMAPTB ../datastreams/GOME2_SIF ../datastreams/LVTbenchmarkOUT ../datastreams/LIS6out ../datastreams/LISDAdiag ../datastreams/Daymet ../datastreams/CMORPH ../datastreams/CHIRPSv2 ../datastreams/3B42V7 ../datastreams/USCRNsm ../datastreams/GLEAM ../datastreams/USDM ../datastreams/LVTpercentile ../datastreams/IMD_PRCP ../datastreams/APHRO_PRCP ../datastreams/GLASSlai ../datastreams/MODISsportLAI ../datastreams/FLUXCOM ../datastreams/HAR ../datastreams/OCO2_SIF ../datastreams/ECMWFforc ../datastreams/GDASforc ../datastreams/ASO_SWE ../runmodes/DataComp ../runmodes/557post ../runmodes/LISpost ../runmodes/DAstats ../runmodes/optUE ../runmodes/Benchmarking ../runmodes/USAFSIpost ../training/LinearRegression/ ../datastreams/GLASSalbedo ../datastreams/IMERG ../datastreams/IMERG_monthly ../datastreams/UA_SNOW ../datastreams/OzFlux ../datastreams/jasmin ../datastreams/MCD15A2H ../datastreams/ERA5 ../datastreams/FluxSat_GPP ../datastreams/THySM/ ../datastreams/GRUNrunoff ../datastreams/UASMAP/ ../datastreams/COAMPSout/ ../datastreams/SMAP_E_OPL/ diff --git a/lvt/plugins/LVT_datastream_pluginMod.F90 b/lvt/plugins/LVT_datastream_pluginMod.F90 index 1016c4e7b..8a5302abc 100644 --- a/lvt/plugins/LVT_datastream_pluginMod.F90 +++ b/lvt/plugins/LVT_datastream_pluginMod.F90 @@ -172,6 +172,7 @@ subroutine LVT_datastream_plugin use GDASforc_dataMod, only : GDASforc_datainit use ASOSWE_obsMod, only : ASOSWE_obsinit use IMERG_dataMod, only : IMERG_datainit + use IMERG_monthly_dataMod, only : IMERG_monthly_datainit use UASNOW_obsMod, only : UASNOW_obsinit use OzFlux_obsMod, only : OzFlux_obsinit use JASMINsm_obsMod, only : JASMINsm_obsInit @@ -284,6 +285,7 @@ subroutine LVT_datastream_plugin external readGDASforcdata external readASOSWEObs external readIMERGdata + external readIMERGmonthlydata external readUASNOWObs external readOzFluxObs external readJASMINsmobs @@ -700,7 +702,12 @@ subroutine LVT_datastream_plugin readASOSWEObs) call registerobssetup(trim(LVT_IMERGdataId)//char(0), IMERG_datainit) - call registerobsread(trim(LVT_IMERGdataId)//char(0) , readIMERGdata) + call registerobsread(trim(LVT_IMERGdataId)//char(0), readIMERGdata) + + call registerobssetup(trim(LVT_IMERGmonthlydataId)//char(0), & + IMERG_monthly_datainit) + call registerobsread(trim(LVT_IMERGmonthlydataId)//char(0), & + readIMERGmonthlydata) call registerobssetup(trim(LVT_UASNOWdataId)//char(0), UASNOW_obsinit) call registerobsread(trim(LVT_UASNOWdataId)//char(0) , readUASNOWObs) diff --git a/lvt/plugins/LVT_pluginIndices.F90 b/lvt/plugins/LVT_pluginIndices.F90 index e188f6197..e596e60eb 100644 --- a/lvt/plugins/LVT_pluginIndices.F90 +++ b/lvt/plugins/LVT_pluginIndices.F90 @@ -262,6 +262,8 @@ module LVT_pluginIndices character*50, public, parameter :: LVT_GDASdataId = "GDAS" character*50, public, parameter :: LVT_ASOSWEdataId = "ASO SWE" character*50, public, parameter :: LVT_IMERGdataId = "GPM IMERG" + character*50, public, parameter :: LVT_IMERGMonthlydataId = & + "GPM IMERG Monthly" character*50, public, parameter :: LVT_UASNOWdataId = "UA SNOW" character*50, public, parameter :: LVT_ozFluxdataId = "OzFlux" character*50, public, parameter :: LVT_JASMINsmobsId = "JASMIN soil moisture" diff --git a/lvt/utils/usaf/nafpa/biasratio.cfg b/lvt/utils/usaf/nafpa/biasratio.cfg new file mode 100644 index 000000000..c8ea41fa9 --- /dev/null +++ b/lvt/utils/usaf/nafpa/biasratio.cfg @@ -0,0 +1,17 @@ +# +# Entries for calc_gfsgalwem_biasratios_multiyear.py +# + +[Input] +gfsdir: /discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/GFS_NAFPA_Monthly_all + +galwemdir: /discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/GALWEM_NAFPA_Monthly_all + +imergdir: /discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/IMERGF_V07A_NAFPA_Monthly + +[Output] +outdir: biasratios + + + + diff --git a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py new file mode 100755 index 000000000..4d2e12b45 --- /dev/null +++ b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py @@ -0,0 +1,355 @@ +#!/usr/bin/env python3 + +""" +SCRIPT: calc_gfsgalwem_biasratios_multiyear.py + +Calculates bias ratio fields using monthly IMERG-FR V07A, GFS, and GALWEM +data, all already interpolated to the NAFPA grid and summed to monthly +totals via LVT. Outputs 12 fields total, each covering a month based +on multiple years of data. + +REVISION HISTORY: +04 Oct 2023: Eric Kemp: Initial specification. +06 Oct 2023: Eric Kemp: Added log10 transform to output for plotting. + Refactored to please pylint. +""" + +import configparser +import datetime +import os +import sys + +# Disable false alarm in pylint (Dataset is not visible in netCDF4 module +# because it is not implemented in python) +# pylint: disable=no-name-in-module +from netCDF4 import Dataset as nc4_dataset +# pylint: enable=no-name-in-module +import numpy as np + +def _usage(): + """Print usage message for this script.""" + print(f"Usage: {sys.argv[0]} CFGFILE BACKSOURCE STARTDATE ENDDATE") + print(" CFGFILE is path to config file") + print(" BACKSOURCE is GFS or GALWEM") + print(" STARTDATE is start date of calculation (YYYYMM)") + print(" ENDDATE is end date of calculateion (YYYYMM)") + +def _process_cmd_line(): + """Process command line arguments.""" + if len(sys.argv) != 5: + _usage() + sys.exit(1) + + cfgfile = sys.argv[1] + if not os.path.exists(cfgfile): + print(f"[ERR] {cfgfile} does not exist!") + sys.exit(1) + + backsource = sys.argv[2] + if backsource not in ["GFS", "GALWEM"]: + print(f"[ERR] {backsource} is not a valid background source") + print("[ERR] Only GFS and GALWEM recognized") + sys.exit(1) + + yyyymm = sys.argv[3] + year = int(yyyymm[0:4]) + month = int(yyyymm[4:6]) + startdate = datetime.datetime(year, month, 1) + + yyyymm = sys.argv[4] + year = int(yyyymm[0:4]) + month = int(yyyymm[4:6]) + enddate = datetime.datetime(year, month, 1) + + if startdate > enddate: + print("[ERR] STARTDATE is beyond ENDDATE!") + sys.exit(1) + + return cfgfile, backsource, startdate, enddate + +def _process_cfg_file(cfgfile, backsource): + """Processes config file for this script.""" + config = configparser.ConfigParser() + config.read(cfgfile) + + if backsource == "GFS": + backindir = config.get('Input', 'gfsdir') + elif backsource == "GALWEM": + backindir = config.get('Input', 'galwemdir') + imergdir = config.get('Input', 'imergdir') + outdir = config.get('Output', 'outdir') + + return backindir, imergdir, outdir + +def _create_arrays(): + """Allocate and zero out numpy arrays""" + nlat = 1920 + nlon = 2560 + nmon = 12 + sum_blended = np.zeros([nmon,nlat,nlon]) + sum_back = np.zeros([nmon,nlat,nlon]) + precip_ratio = np.zeros([nmon,nlat,nlon]) + return sum_blended, sum_back, precip_ratio + +def _read_imerg_file(imergdir, nextdate): + """Read the interpolated monthly IMERG file.""" + filename = f"{imergdir}/SUM_TS." + filename += f"{nextdate.year:04d}{nextdate.month:02d}" + filename += f"{nextdate.day:02d}0000.d01.nc" + ncid_imerg = nc4_dataset(filename, mode='r', \ + format="NETCDF4_CLASSIC") + precip_imerg = ncid_imerg.variables['TotalPrecip'][:,:] + del ncid_imerg + return precip_imerg + +def _read_imerg_file_latlondims(imergdir, startdate): + """Read the interpolated monthly IMERG file.""" + nextdate = _get_nextdate(startdate) + filename = f"{imergdir}/SUM_TS." + filename += f"{nextdate.year:04d}{nextdate.month:02d}" + filename += f"{nextdate.day:02d}0000.d01.nc" + ncid_imerg = nc4_dataset(filename, mode='r', \ + format="NETCDF4_CLASSIC") + lats_imerg = ncid_imerg.variables["latitude"][:,:] + lons_imerg = ncid_imerg.variables["longitude"][:,:] + north_south = ncid_imerg.dimensions['north_south'].size + east_west = ncid_imerg.dimensions['east_west'].size + del ncid_imerg + return lats_imerg, lons_imerg, north_south, east_west + +def _read_background_file(backindir, nextdate): + """Read interpolated monthly background file.""" + filename = f"{backindir}/SUM_TS." + filename += f"{nextdate.year:04d}{nextdate.month:02d}" + filename += f"{nextdate.day:02d}0000.d01.nc" + ncid_back = nc4_dataset(filename, mode='r', \ + format="NETCDF4_CLASSIC") + precip_back = ncid_back.variables['TotalPrecip'][:,:] + del ncid_back + return precip_back + +def _calc_precip_blended(lats_imerg, precip_imerg, precip_back): + """Calculate weighted blend of IMERG and background precip""" + + # Use IMERG from 40S to 71N; use linear tapers from 40S to 60S, + # and 51N to 71N; and don't use IMERG south of 60S and north of + # 71N. Rationale: No IR and little gauge data is available + # south of 60S, so we don't expect IMERG to be useful here for + # bias correction. In the northern hemisphere, there is no IR data + # north of 60N, but there is good GPCC gauge density in Scandinavia + # up to 71N. Finally, the 20 degree latitude linear tapers mimicks + # MERRA-2 usage of CPCU gauge analyses. + precip_imerg_weights = np.ones(np.shape(lats_imerg)) + precip_imerg_weights = np.where(lats_imerg < -60, + 0, + precip_imerg_weights) + precip_imerg_weights = np.where( (lats_imerg >= -60) & + (lats_imerg < -40), + ( (lats_imerg + 60.) / 20.), + precip_imerg_weights) + precip_imerg_weights = np.where( (lats_imerg > 51) & + (lats_imerg <= 71), + ( (71. - lats_imerg) / 20.), + precip_imerg_weights) + precip_imerg_weights = np.where(lats_imerg > 71, + 0, + precip_imerg_weights) + + # Conservative alternative: Linear taper from 40N to 60N, and + # screen out everything north of 60N. We'll do this as a backup + # if we find unphysical IMERG patterns north of 60N + #precip_imerg_weights = np.where( (lats_imerg > 40) & + # (lats_imerg <= 60), + # ( (60. - lats_imerg) / 20.), + # precip_imerg_weights) + #precip_imerg_weights = np.where(lats_imerg > 60, + # 0, + # precip_imerg_weights) + + precip_blended = precip_imerg_weights[:,:]*precip_imerg[:,:] + \ + (1. - precip_imerg_weights[:,:])*precip_back[:,:] + + return precip_blended + +def _get_nextdate(curdate): + """Determine next date (first of next month)""" + # LVT output files have data from the prior month, so we must + # advance one month to find the appropriate file, e.g., data + # for May 2020 will be in SUM_TS.202006010000.d01.nc. Python's + # timedelta object isn't smart enough to jump a whole month, so + # we loop through each day instead. + timedelta = datetime.timedelta(days=1) + nextdate = curdate + timedelta + while nextdate.day != 1: + nextdate += timedelta + return nextdate + +def _calc_biasratios(imergdir, backindir, startdate, enddate, lats_imerg): + """Calculate the bias ratios from the input data files.""" + + sum_blended, sum_back, precip_ratio = _create_arrays() + + # Loop through each month + curdate = startdate + while curdate <= enddate: + + # LVT output files have data from the prior month, so we must + # advance one month to find the appropriate file, e.g., data + # for May 2020 will be in SUM_TS.202006010000.d01.nc. + nextdate = _get_nextdate(curdate) + + # First, the IMERG file + precip_imerg = \ + _read_imerg_file(imergdir, nextdate) + + # Next, the background data (GFS or GALWEM) + precip_back = _read_background_file(backindir, nextdate) + + # Calculate blended precipitation (weighted average of IMERG + # and background, varying by latitude). + precip_blended = \ + _calc_precip_blended(lats_imerg, precip_imerg, precip_back) + + # Updated precip sums. Add trace precipitation every month to + # prevent undefined ratios in deserts. + sum_blended[(curdate.month-1),:,:] += precip_blended[:,:] + 0.05 + sum_back[(curdate.month-1),:,:] += precip_back[:,:] + 0.05 + + # Move on to next month + curdate = nextdate + + # Finish calculation + precip_ratio[:,:,:] = sum_blended[:,:,:] / sum_back[:,:,:] + precip_ratio[:,:,:] = np.where(precip_ratio == 0, 1, precip_ratio) + + return precip_ratio + +def _create_output_filename(outdir, backsource, startdate, enddate): + """Create output netCDF file name""" + filename = f"{outdir}" + filename += f"/{backsource}_pcp_biasratios_" + filename += f"{startdate.year:04d}{startdate.month:02d}_" + filename += f"{enddate.year:04d}{enddate.month:02d}.nc" + return filename + +def _create_latitude(rootgrp): + """Create latitude dataset in output file""" + latitude = rootgrp.createVariable("latitude", "f4", \ + ("north_south", "east_west",)) + latitude.units = "degree_north" + latitude.standard_name = "latitude" + latitude.long_name = "latitude" + latitude.scale_factor = np.float32("1.") + latitude.add_offset = np.float32("0.") + latitude.missing_value = np.float32("-9999.") + return latitude + +def _create_longitude(rootgrp): + """Create longitude dataset in output file.""" + longitude = rootgrp.createVariable("longitude", "f4", \ + ("north_south", "east_west",)) + longitude.units = "degree_east" + longitude.standard_name = "longitude" + longitude.long_name = "longitude" + longitude.scale_factor = np.float32("1.") + longitude.add_offset = np.float32("0.") + longitude.missing_value = np.float32("-9999.") + return longitude + +def _create_bias_ratio(rootgrp, backsource): + """Create the biasRatio dataset in the output file.""" + bias_ratio = rootgrp.createVariable("biasRatio", "f4", \ + ("months", "north_south", "east_west",)) + bias_ratio.units = "-" + bias_ratio.long_name = \ + f"bias_ratio_for_{backsource}_precipitation" + bias_ratio.scale_factor = np.float32("1.") + bias_ratio.add_offset = np.float32("0.") + bias_ratio.missing_value = np.float32("-9999.") + return bias_ratio + +def _create_log10_bias_ratio(rootgrp, backsource): + """Create the log10BiasRatio dataset in the output file.""" + log10_bias_ratio = rootgrp.createVariable("log10BiasRatio", "f4", \ + ("months", "north_south", "east_west",)) + log10_bias_ratio.units = "-" + log10_bias_ratio.long_name = \ + f"log10_bias_ratio_for_{backsource}_precipitation" + log10_bias_ratio.scale_factor = np.float32("1.") + log10_bias_ratio.add_offset = np.float32("0.") + log10_bias_ratio.missing_value = np.float32("-9999.") + return log10_bias_ratio + +def _write_biasratios(args): + """Write out bias ratios to netCDF file""" + + os.makedirs(args['outdir'], exist_ok=True) + + now = datetime.datetime.utcnow() + + history = "created on date: " + history += f"{now.year:04d}-{now.month:02d}-{now.day:02d}" + history += f"T{now.hour:02d}:{now.minute:02d}:{now.second:02d}" + + rootgrp = nc4_dataset(args['outfile'], "w", format="NETCDF4") + rootgrp.missing_value = np.float32("-9999.") + rootgrp.title = f"Monthly bias ratio for IMERG / {args['backsource']}" + rootgrp.institution = "NASA GSFC" + + rootgrp.history = f"created on date: {history}" + rootgrp.comment = "website: http://lis.gsfc.nasa.gov/" + rootgrp.MAP_PROJECTION = "EQUIDISTANT CYLINDRICAL" + rootgrp.SOUTH_WEST_CORNER_LAT = np.float32("-89.95312") + rootgrp.SOUTH_WEST_CORNER_LON = np.float32("-179.9297") + rootgrp.DX = np.float32("0.140625") + rootgrp.DY = np.float32("0.09375") + + # Define dimensions + rootgrp.createDimension("months", 12) + rootgrp.createDimension("north_south", \ + args['north_south']) + rootgrp.createDimension("east_west", \ + args['east_west']) + + # Define output variables + latitude = _create_latitude(rootgrp) + longitude = _create_longitude(rootgrp) + bias_ratio = _create_bias_ratio(rootgrp, args['backsource']) + log10_bias_ratio = _create_log10_bias_ratio(rootgrp, + args['backsource']) + + latitude[:,:] = args['lats_imerg'][:,:] + longitude[:,:] = args['lons_imerg'][:,:] + bias_ratio[:,:,:] = args['precip_ratio'][:,:,:] + log10_bias_ratio[:,:,:] = np.log10(args['precip_ratio'][:,:,:]) + + rootgrp.close() + +def _main(): + """Main driver""" + + cfgfile, backsource, startdate, enddate = _process_cmd_line() + backindir, imergdir, outdir = \ + _process_cfg_file(cfgfile, backsource) + lats_imerg, lons_imerg, north_south, east_west = \ + _read_imerg_file_latlondims(imergdir, startdate) + precip_ratio = \ + _calc_biasratios(imergdir, backindir, startdate, enddate, \ + lats_imerg) + outfile = _create_output_filename(outdir, backsource, + startdate, enddate) + # To satisfy pylint.... + args = { + "outfile" : outfile, + "outdir" : outdir, + "backsource" : backsource, + "north_south" : north_south, + "east_west" : east_west, + "lats_imerg" : lats_imerg, + "lons_imerg" : lons_imerg, + "precip_ratio" : precip_ratio, + } + _write_biasratios(args) + +if __name__ == "__main__": + _main()