diff --git a/ldt/configs/ldt.config.adoc b/ldt/configs/ldt.config.adoc index 4f72b49c1f..be9a3887af 100644 --- a/ldt/configs/ldt.config.adoc +++ b/ldt/configs/ldt.config.adoc @@ -182,6 +182,7 @@ endif::devonly[] |"`TRMM 3B42V6`" |TRMM-based 3B42 V6 rainfall |"`TRMM 3B42V7`" |TRMM-based 3B42 V7 rainfall |"`ERA5`" |ERA5 reanalysis +|"`ERA5CDS`" |ERA5 reanalysis from Climate Data Store |=== .Example _ldt.config_ entry @@ -3292,6 +3293,27 @@ ERA5 forcing tile to grid mapping file: ../ERA5/mapping.nc ERA5 forcing terrain height file: ../ERA5/era5_elev.nc .... +==== ERA5CDS Reanalysis Forcing based parameter inputs + +`ERA5CDS forcing directory:` specifies the location of the puglic ERA5 Reanalysis forcing data files. + +`ERA5CDS forcing terrain height file:` specifies the ERA5CDS reanalysis elevation file. + +`ERA5CDS use lowest model level forcing:` specifies the ERA5CDS reanalysis level of instantaneous fields. +[cols="<,<",] +|=== +|Value |Description + +|"`0`" |Use 2m temperature, 2m dew point temperature, and 10m u/v wind +|"`1`" |Use model lowest level temperature, specific humidity, and u/v wind +|=== + +.Example _ldt.config_ entry +.... +ERA5CDS forcing directory: ./ERA5/single_hourly +ERA5CDS forcing terrain height file: ../ERA5/era5_elev.nc +ERA5CDS use lowest model level forcing: 1 +.... ==== MERRA-2 Forcing based parameter inputs diff --git a/ldt/interp/zterp.F90 b/ldt/interp/zterp.F90 index 974e7c3ad7..20e32b3d9a 100644 --- a/ldt/interp/zterp.F90 +++ b/ldt/interp/zterp.F90 @@ -26,13 +26,16 @@ ! 4/30/15 Hiroko Beaudoing; Changed average solar zenith angle ! calculation for iflag=0 case. During sunrise/sunset, weight1 ! becomes huge due to extremely small avgangle, resulting in -! unrealdttically large SWdown flux at the instance. This fix +! unrealistically large SWdown flux at the instance. This fix ! accounts for time interval consisting of fewer than 37 ! above-zero czavgdata. ! 9/24/15 Hiroko Beaudoing; Modified above fix to be applied only when ! less than 50% of 37 intervals (arbitary) are day-time. This ! modification is needed to properly interpolate ECMWF data, ! where daily peak is too low due to the nature of dataset. +! 3/4/26 Hiroko Beaudoing; Additional fix for iflag=0 case, when weight1 +! becomes huge due to extremely small avgangle, resulting in +! unrealistically large SWdown flux at the instance. ! ! !INTERFACE: subroutine zterp (iflag,lat,lon,btime,etime, & @@ -229,6 +232,11 @@ subroutine zterp (iflag,lat,lon,btime,etime, & else weight1=(czmodel/avgangle) endif + if ( ldt%zterp_correction ) then + if (weight1.gt.5.0) then + weight1=1.0 + endif + endif endif !-------------------------------------------------------------------- ! Second case: IFLAG=1 (instantaneous input and output) diff --git a/ldt/make/Filepath b/ldt/make/Filepath index 3a469d2e06..bc9e8a7faf 100644 --- a/ldt/make/Filepath +++ b/ldt/make/Filepath @@ -1 +1 @@ -dirs := . ../main ../core ../interp ../plugins ../lib/bil ../lib/xmrg ../params/mask ../params/landcover ../params/soils ../params/topo ../params/albedo ../params/gfrac ../params/LAISAI ../params/slopetype ../params/tbot ../params/pet ../params/climateparms ../params/crop_parms ../params/irrigation ../params/catchment ../params/CLM2 ../params/GeoWRSI ../params/FLAKE ../params/HYMAP ../params/RAPID ../params/JULES ../params/Mosaic ../params/Noah ../params/AquaCrop ../params/RDHM ../params/RUC ../params/SACHTET ../params/Snow17 ../params/SnowModel ../params/SiB2 ../params/VIC ../params/glacier ../runmodes/ANNproc ../runmodes/DApreproc ../runmodes/LSMparamproc ../runmodes/EnsRstpreproc ../runmodes/NUWRFpreproc ../runmodes/Metforcproc ../runmodes/MetTimeDScale ../runmodes/climoRstproc ../runmodes/rstTransformProc ../runmodes/StatDscaleMetForc ../runmodes/OPTUEproc/ ../domains/latlon ../domains/gaussian ../domains/hrap ../domains/lambert ../domains/polar ../domains/merc ../domains/easev2 ../DAobs/MCD15A2H_LAI ../DAobs/LISlsmSMobs ../DAobs/LISlsmPrecipobs ../DAobs/LISlsmTEFFobs ../DAobs/synthetic_sm ../DAobs/NASA_AMSRE_sm ../DAobs/LPRM_AMSREsm ../DAobs/LPRMvod ../DAobs/ESACCI_sm ../DAobs/WindSat_sm ../DAobs/GRACE_tws ../DAobs/GRACEQL_tws ../DAobs/SMOPS ../DAobs/ASCAT_TUW ../DAobs/SMOS_L2sm ../DAobs/SMOS_NESDIS ../DAobs/GCOMW_AMSR2L3snd ../DAobs/GCOMW_AMSR2L3sm ../DAobs/Aquarius_L2sm ../DAobs/simGRACE_JPL ../DAobs/SMMR_SNWD ../DAobs/SSMI_SNWD ../DAobs/ANSA_SNWD ../DAobs/NASA_SMAPsm ../DAobs/SMOS_NRTNN_L2sm ../DAobs/NASA_SMAPvod ../DAobs/GLASSlai ../DAobs/THySM ../DAobs/GEOS_FP_TEFFobs ../DAobs/SMAP_E_OPLsm ../DAobs/VIIRS_GVF ../DAobs/CDFS_GVF ../RESTRICTED/usaf/usaf ../metforcing/agrradps ../metforcing/cmap ../metforcing/cmorph ../metforcing/ecmwf ../metforcing/gdas ../metforcing/geos5fcst ../metforcing/gfs ../metforcing/gldas ../metforcing/gswp1 ../metforcing/gswp2 ../metforcing/merra2 ../metforcing/era5 ../metforcing/nam242 ../metforcing/narr ../metforcing/princeton ../metforcing/nldas2 ../metforcing/WRFoutv2 ../metforcing/WRFAKdom ../metforcing/stg2 ../metforcing/stg4 ../metforcing/3B42RTV7 ../metforcing/3B42V6 ../metforcing/3B42V7 ../metforcing/RFE2Daily ../metforcing/RFE2gdas ../metforcing/chirps2 ../MetforcScale ../ANNdata/LISlsmSM ../ANNdata/synthetic_sm ../ANNdata/LPRM_AMSREsm ../ANNdata/MODIS_LST ../ANNdata/MOD10A1 ../ANNdata/GHCNsnwd ../ANNdata/GCOMW_AMSR2_TB ../statDscale/BayesianMerging ../statDscale/Climo ../params/CLM45 ../USAFSI ../SNIP ../SMAP_E_OPL ../runmodes/USAFSI ../runmodes/SNIP ../runmodes/LISHydropreproc ../runmodes/SMAP_E_OPL ../params/Crocus ../runmodes/obsSim ../obsSim/NatureRun/LISout/ ../obsSim/OSSEmask/LISout/ ../obsSim/OSSEmask/AMSR2/ ../obsSim/OSSEmask/TSMM/ ../obsSim/OSSEmask/MODIS/ ../obsSim/OSSEmask/Sentinel1A/ ../DAobs/WSF_sm ../WSF_OPL ../runmodes/WSF_OPL +dirs := . ../main ../core ../interp ../plugins ../lib/bil ../lib/xmrg ../params/mask ../params/landcover ../params/soils ../params/topo ../params/albedo ../params/gfrac ../params/LAISAI ../params/slopetype ../params/tbot ../params/pet ../params/climateparms ../params/crop_parms ../params/irrigation ../params/catchment ../params/CLM2 ../params/GeoWRSI ../params/FLAKE ../params/HYMAP ../params/RAPID ../params/JULES ../params/Mosaic ../params/Noah ../params/AquaCrop ../params/RDHM ../params/RUC ../params/SACHTET ../params/Snow17 ../params/SnowModel ../params/SiB2 ../params/VIC ../params/glacier ../runmodes/ANNproc ../runmodes/DApreproc ../runmodes/LSMparamproc ../runmodes/EnsRstpreproc ../runmodes/NUWRFpreproc ../runmodes/Metforcproc ../runmodes/MetTimeDScale ../runmodes/climoRstproc ../runmodes/rstTransformProc ../runmodes/StatDscaleMetForc ../runmodes/OPTUEproc/ ../domains/latlon ../domains/gaussian ../domains/hrap ../domains/lambert ../domains/polar ../domains/merc ../domains/easev2 ../DAobs/MCD15A2H_LAI ../DAobs/LISlsmSMobs ../DAobs/LISlsmPrecipobs ../DAobs/LISlsmTEFFobs ../DAobs/synthetic_sm ../DAobs/NASA_AMSRE_sm ../DAobs/LPRM_AMSREsm ../DAobs/LPRMvod ../DAobs/ESACCI_sm ../DAobs/WindSat_sm ../DAobs/GRACE_tws ../DAobs/GRACEQL_tws ../DAobs/SMOPS ../DAobs/ASCAT_TUW ../DAobs/SMOS_L2sm ../DAobs/SMOS_NESDIS ../DAobs/GCOMW_AMSR2L3snd ../DAobs/GCOMW_AMSR2L3sm ../DAobs/Aquarius_L2sm ../DAobs/simGRACE_JPL ../DAobs/SMMR_SNWD ../DAobs/SSMI_SNWD ../DAobs/ANSA_SNWD ../DAobs/NASA_SMAPsm ../DAobs/SMOS_NRTNN_L2sm ../DAobs/NASA_SMAPvod ../DAobs/GLASSlai ../DAobs/THySM ../DAobs/GEOS_FP_TEFFobs ../DAobs/SMAP_E_OPLsm ../DAobs/VIIRS_GVF ../DAobs/CDFS_GVF ../RESTRICTED/usaf/usaf ../metforcing/agrradps ../metforcing/cmap ../metforcing/cmorph ../metforcing/ecmwf ../metforcing/gdas ../metforcing/geos5fcst ../metforcing/gfs ../metforcing/gldas ../metforcing/gswp1 ../metforcing/gswp2 ../metforcing/merra2 ../metforcing/era5 ../metforcing/nam242 ../metforcing/narr ../metforcing/princeton ../metforcing/nldas2 ../metforcing/WRFoutv2 ../metforcing/WRFAKdom ../metforcing/stg2 ../metforcing/stg4 ../metforcing/3B42RTV7 ../metforcing/3B42V6 ../metforcing/3B42V7 ../metforcing/RFE2Daily ../metforcing/RFE2gdas ../metforcing/chirps2 ../MetforcScale ../ANNdata/LISlsmSM ../ANNdata/synthetic_sm ../ANNdata/LPRM_AMSREsm ../ANNdata/MODIS_LST ../ANNdata/MOD10A1 ../ANNdata/GHCNsnwd ../ANNdata/GCOMW_AMSR2_TB ../statDscale/BayesianMerging ../statDscale/Climo ../params/CLM45 ../USAFSI ../SNIP ../SMAP_E_OPL ../runmodes/USAFSI ../runmodes/SNIP ../runmodes/LISHydropreproc ../runmodes/SMAP_E_OPL ../params/Crocus ../runmodes/obsSim ../obsSim/NatureRun/LISout/ ../obsSim/OSSEmask/LISout/ ../obsSim/OSSEmask/AMSR2/ ../obsSim/OSSEmask/TSMM/ ../obsSim/OSSEmask/MODIS/ ../obsSim/OSSEmask/Sentinel1A/ ../DAobs/WSF_sm ../WSF_OPL ../runmodes/WSF_OPL ../metforcing/era5cds diff --git a/ldt/metforcing/era5cds/era5cds_forcingMod.F90 b/ldt/metforcing/era5cds/era5cds_forcingMod.F90 new file mode 100644 index 0000000000..9b1c3a946d --- /dev/null +++ b/ldt/metforcing/era5cds/era5cds_forcingMod.F90 @@ -0,0 +1,339 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LDT_misc.h" +module era5cds_forcingMod +!BOP +! !MODULE: era5cds_forcingMod +! +! !DESCRIPTION: +! This module contains variables and data structures that are used +! for the implementation of the ERA5 forcing data from Climate Data Store. +! The data is global 0.25 degree dataset in latlon +! projection, and at 1 hourly intervals. The derived +! data type {\tt era5cds\_struc} +! includes the variables that specify the runtime options, and the +! weights and neighbor information to be used for spatial interpolation. +! They are described below: +! \begin{description} +! \item[ncold] +! Number of columns (along the east west dimension) for the input data +! \item[nrold] +! Number of rows (along the north south dimension) for the input data +! \item[nmif] +! Number of forcing variables in the ECMWF data +! \item[era5cdstime1] +! The nearest, previous 1 hour instance of the incoming +! data (as a real time). +! \item[era5cdstime2] +! The nearest, next 1 hour instance of the incoming +! data (as a real time). +! \item[era5cdsdir] +! Directory containing the input data +! \item[mi] +! Number of points in the input grid +! \item[n111,n121,n211,n221] +! Arrays containing the neighbor information of the input grid +! for each grid point in LDT, for bilinear interpolation. +! \item[w111,w121,w211,w221] +! Arrays containing the weights of the input grid +! for each grid point in LDT, for bilinear interpolation. +! \item[n122,n122,n212,n222] +! Arrays containing the neighbor information of the input grid +! for each grid point in LDT, for conservative interpolation. +! \item[w112,w122,w212,w222] +! Arrays containing the weights of the input grid +! for each grid point in LDT, for conservative interpolation. +! \item[n113] +! Arrays containing the neighbor information of the input grid +! for each grid point in LDT, for n. neighbor interpolation. +! \item[findtime1, findtime2] +! boolean flags to indicate which time is to be read for +! temporal interpolation. +! \end{description} +! +! !USES: + use LDT_constantsMod, only : LDT_CONST_PATH_LEN + implicit none + + PRIVATE +!----------------------------------------------------------------------------- +! !PUBLIC MEMBER FUNCTIONS: +!----------------------------------------------------------------------------- + public :: init_era5cds !defines the native resolution of + !the input data +!----------------------------------------------------------------------------- +! !PUBLIC TYPES: +!----------------------------------------------------------------------------- + public :: era5cds_struc + +!EOP + type, public :: era5cds_type_dec + + real :: ts + integer :: ncold, nrold + character(len=LDT_CONST_PATH_LEN) :: era5cdsdir !ERA5 Forcing Directory + character(len=LDT_CONST_PATH_LEN) :: era5cdshgt_file + real*8 :: era5cdstime1,era5cdstime2 + logical :: reset_flag + integer :: mon + + integer :: mi + integer, allocatable :: n111(:) + integer, allocatable :: n121(:) + integer, allocatable :: n211(:) + integer, allocatable :: n221(:) + real, allocatable :: w111(:),w121(:) + real, allocatable :: w211(:),w221(:) + + integer, allocatable :: n112(:,:) + integer, allocatable :: n122(:,:) + integer, allocatable :: n212(:,:) + integer, allocatable :: n222(:,:) + real, allocatable :: w112(:,:),w122(:,:) + real, allocatable :: w212(:,:),w222(:,:) + integer, allocatable :: n113(:) + + integer :: findtime1, findtime2 + logical :: startFlag, dayFlag + real*8 :: validstart + + real, allocatable :: tair(:,:) + real, allocatable :: qair(:,:) + real, allocatable :: uwind(:,:) + real, allocatable :: vwind(:,:) + real, allocatable :: ps(:,:) + real, allocatable :: rainf(:,:) + real, allocatable :: crainf(:,:) + real, allocatable :: swd(:,:) + real, allocatable :: lwd(:,:) + + real, allocatable :: prev_rainf(:,:) + real, allocatable :: prev_crainf(:,:) + real, allocatable :: prev_swd(:,:) + real, allocatable :: prev_lwd(:,:) + + integer :: nmif + integer :: uselml + integer :: tdimsize + + real*8 :: ringtime + + end type era5cds_type_dec + + type(era5cds_type_dec), allocatable :: era5cds_struc(:) + +contains + +!BOP +! +! !ROUTINE: init_era5cds +! \label{init_era5cds} +! +! !REVISION HISTORY: +! 23 Dec 2019: Sujay Kumar, initial code +! 17 Apr 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: + subroutine init_era5cds(findex) + +! !USES: + use LDT_coreMod, only : LDT_rc, LDT_domain + use LDT_logMod, only : LDT_logunit, LDT_endrun + use LDT_timeMgrMod + use map_utils, only : proj_latlon + +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none +! !AGRUMENTS: + integer, intent(in) :: findex +! +! !DESCRIPTION: +! Defines the native resolution of the input forcing for ERA5 +! data. The grid description arrays are based on the decoding +! schemes used by NCEP and followed in the LDT interpolation +! schemes (see Section~\ref{interp}). +! +! The routines invoked are: +! \begin{description} +! \item[readcrd\_era5cds](\ref{readcrd_era5cds}) \newline +! reads the runtime options specified for ERA5 data +! \item[bilinear\_interp\_input](\ref{bilinear_interp_input}) \newline +! computes the neighbor, weights for bilinear interpolation +! \item[conserv\_interp\_input](\ref{conserv_interp_input}) \newline +! computes the neighbor, weights for conservative interpolation +! \end{description} +!EOP + + integer :: updoy, yr1,mo1,da1,hr1,mn1,ss1 + real :: upgmt + integer :: n + integer :: ftn + real :: gridDesc(20) + + allocate(era5cds_struc(LDT_rc%nnest)) + + write(LDT_logunit,*) "Reading the ERA5CDS forcing data" + +!- Read in config ERA5CDS inputs: + call readcrd_era5cds(findex) + + do n=1, LDT_rc%nnest + era5cds_struc(n)%ts = 3600 !hour + call LDT_update_timestep(LDT_rc, n, era5cds_struc(n)%ts) + enddo + + era5cds_struc%reset_flag = .false. + era5cds_struc(:)%nmif = 9 + + ! Metforcing and parameter grid info: + LDT_rc%met_proj(findex) = "latlon" + + era5cds_struc(:)%ncold = 1440 + era5cds_struc(:)%nrold = 720 + gridDesc = 0 + gridDesc(1) = 0 + gridDesc(2) = era5cds_struc(1)%ncold + gridDesc(3) = era5cds_struc(1)%nrold + gridDesc(4) = -89.875 + gridDesc(5) = -179.875 + gridDesc(6) = 128 + gridDesc(7) = 89.875 + gridDesc(8) = 179.875 + gridDesc(9) = 0.25 + gridDesc(10) = 0.25 + gridDesc(20) = 0 + + LDT_rc%met_nc(findex) = era5cds_struc(1)%ncold + LDT_rc%met_nr(findex) = era5cds_struc(1)%nrold + LDT_rc%met_gridDesc(findex,1:20) = gridDesc(1:20) + + !- If only processing parameters, then return to main routine calls ... + if( LDT_rc%runmode == "LSM parameter processing" ) return + + LDT_rc%met_nf(findex) = 9 + LDT_rc%met_ts(findex) = 3600 + LDT_rc%met_zterp(findex) = .true. + + do n=1,LDT_rc%nnest + + era5cds_struc(n)%findtime1 = 0 + era5cds_struc(n)%findtime2 = 0 + + era5cds_struc(n)%mi = era5cds_struc(n)%ncold*era5cds_struc(n)%nrold + + ! ERA5 accumulation data starts at 7z on 1940-01-01 + yr1 = 1940 + mo1 = 01 + da1 = 01 + hr1 = 07 + mn1 = 0; ss1 = 0 + call LDT_date2time( era5cds_struc(n)%validstart,updoy,upgmt,yr1,mo1,da1,hr1,mn1,ss1 ) + + if( gridDesc(9) == LDT_rc%gridDesc(n,9) .and. & + gridDesc(10) == LDT_rc%gridDesc(n,10).and. & + LDT_rc%gridDesc(n,1) == proj_latlon .and. & + LDT_rc%met_gridtransform(findex) .ne. "neighbor" ) then + write(LDT_logunit,*) "[ERR] The ERA5CDS 0.25 deg grid was selected for the" + write(LDT_logunit,*) " LDT run domain; however, 'bilinear', 'budget-bilinear'," + write(LDT_logunit,*) " or some other unknown option was selected to spatially" + write(LDT_logunit,*) " downscale the grid, which will cause errors during runtime." + write(LDT_logunit,*) "Program stopping ..." + call LDT_endrun() + endif + + ! Setting up weights for Interpolation + select case( LDT_rc%met_gridtransform(findex) ) + case("bilinear") + allocate(era5cds_struc(n)%n111(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%n121(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%n211(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%n221(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%w111(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%w121(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%w211(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%w221(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + call bilinear_interp_input(n, gridDesc(:),& + era5cds_struc(n)%n111,era5cds_struc(n)%n121,& + era5cds_struc(n)%n211,era5cds_struc(n)%n221,& + era5cds_struc(n)%w111,era5cds_struc(n)%w121,& + era5cds_struc(n)%w211,era5cds_struc(n)%w221) + + case("budget-bilinear") + allocate(era5cds_struc(n)%n111(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%n121(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%n211(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%n221(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%w111(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%w121(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%w211(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(era5cds_struc(n)%w221(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + call bilinear_interp_input(n, gridDesc(:),& + era5cds_struc(n)%n111,era5cds_struc(n)%n121,& + era5cds_struc(n)%n211,era5cds_struc(n)%n221,& + era5cds_struc(n)%w111,era5cds_struc(n)%w121,& + era5cds_struc(n)%w211,era5cds_struc(n)%w221) + + allocate(era5cds_struc(n)%n112(LDT_rc%lnc(n)*LDT_rc%lnr(n),25)) + allocate(era5cds_struc(n)%n122(LDT_rc%lnc(n)*LDT_rc%lnr(n),25)) + allocate(era5cds_struc(n)%n212(LDT_rc%lnc(n)*LDT_rc%lnr(n),25)) + allocate(era5cds_struc(n)%n222(LDT_rc%lnc(n)*LDT_rc%lnr(n),25)) + allocate(era5cds_struc(n)%w112(LDT_rc%lnc(n)*LDT_rc%lnr(n),25)) + allocate(era5cds_struc(n)%w122(LDT_rc%lnc(n)*LDT_rc%lnr(n),25)) + allocate(era5cds_struc(n)%w212(LDT_rc%lnc(n)*LDT_rc%lnr(n),25)) + allocate(era5cds_struc(n)%w222(LDT_rc%lnc(n)*LDT_rc%lnr(n),25)) + call conserv_interp_input(n, gridDesc(:),& + era5cds_struc(n)%n112,era5cds_struc(n)%n122,& + era5cds_struc(n)%n212,era5cds_struc(n)%n222,& + era5cds_struc(n)%w112,era5cds_struc(n)%w122,& + era5cds_struc(n)%w212,era5cds_struc(n)%w222) + + case("neighbor") + allocate(era5cds_struc(n)%n113(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + call neighbor_interp_input(n, gridDesc(:),& + era5cds_struc(n)%n113) + + case default + write(LDT_logunit,*) '[ERR] Interpolation option '// & + trim(LDT_rc%met_gridtransform(findex))//& + ' for ERA5CDS forcing is not supported' + call LDT_endrun() + end select + + call LDT_registerAlarm("ERA5CDS forcing alarm",& + 86400.0,86400.0) + era5cds_struc(n)%startFlag = .true. + era5cds_struc(n)%dayFlag = .true. + + era5cds_struc(n)%mon = -1 + allocate(era5cds_struc(n)%tair(LDT_rc%lnc(n)*LDT_rc%lnr(n),745)) + allocate(era5cds_struc(n)%qair(LDT_rc%lnc(n)*LDT_rc%lnr(n),745)) + allocate(era5cds_struc(n)%uwind(LDT_rc%lnc(n)*LDT_rc%lnr(n),745)) + allocate(era5cds_struc(n)%vwind(LDT_rc%lnc(n)*LDT_rc%lnr(n),745)) + allocate(era5cds_struc(n)%ps(LDT_rc%lnc(n)*LDT_rc%lnr(n),745)) + allocate(era5cds_struc(n)%rainf(LDT_rc%lnc(n)*LDT_rc%lnr(n),745)) + allocate(era5cds_struc(n)%crainf(LDT_rc%lnc(n)*LDT_rc%lnr(n),745)) + allocate(era5cds_struc(n)%swd(LDT_rc%lnc(n)*LDT_rc%lnr(n),745)) + allocate(era5cds_struc(n)%lwd(LDT_rc%lnc(n)*LDT_rc%lnr(n),745)) + + allocate(era5cds_struc(n)%prev_rainf(LDT_rc%lnc(n)*LDT_rc%lnr(n),7)) + allocate(era5cds_struc(n)%prev_crainf(LDT_rc%lnc(n)*LDT_rc%lnr(n),7)) + allocate(era5cds_struc(n)%prev_swd(LDT_rc%lnc(n)*LDT_rc%lnr(n),7)) + allocate(era5cds_struc(n)%prev_lwd(LDT_rc%lnc(n)*LDT_rc%lnr(n),7)) + + enddo ! End nest loop + + + end subroutine init_era5cds +end module era5cds_forcingMod + diff --git a/ldt/metforcing/era5cds/finalize_era5cds.F90 b/ldt/metforcing/era5cds/finalize_era5cds.F90 new file mode 100644 index 0000000000..8f295456f9 --- /dev/null +++ b/ldt/metforcing/era5cds/finalize_era5cds.F90 @@ -0,0 +1,75 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! !ROUTINE: finalize_era5cds +! \label{finalize_era5cds} +! +! !REVISION HISTORY: +! 23 Dec 2019: Sujay Kumar, initial code +! 17 Apr 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: + +subroutine finalize_era5cds(findex) + +! !USES: + use LDT_coreMod, only : LDT_rc + use era5cds_forcingMod, only : era5cds_struc +! +! !DESCRIPTION: +! Routine to cleanup ERA5 forcing related memory allocations. +! +!EOP + implicit none + + integer :: findex + integer :: n + + do n=1,LDT_rc%nnest + select case( LDT_rc%met_gridtransform(findex) ) + case( "bilinear" ) + deallocate(era5cds_struc(n)%n111) + deallocate(era5cds_struc(n)%n121) + deallocate(era5cds_struc(n)%n211) + deallocate(era5cds_struc(n)%n221) + deallocate(era5cds_struc(n)%w111) + deallocate(era5cds_struc(n)%w121) + deallocate(era5cds_struc(n)%w211) + deallocate(era5cds_struc(n)%w221) + + case( "budget-bilinear" ) + deallocate(era5cds_struc(n)%n111) + deallocate(era5cds_struc(n)%n121) + deallocate(era5cds_struc(n)%n211) + deallocate(era5cds_struc(n)%n221) + deallocate(era5cds_struc(n)%w111) + deallocate(era5cds_struc(n)%w121) + deallocate(era5cds_struc(n)%w211) + deallocate(era5cds_struc(n)%w221) + deallocate(era5cds_struc(n)%n112) + deallocate(era5cds_struc(n)%n122) + deallocate(era5cds_struc(n)%n212) + deallocate(era5cds_struc(n)%n222) + deallocate(era5cds_struc(n)%w112) + deallocate(era5cds_struc(n)%w122) + deallocate(era5cds_struc(n)%w212) + deallocate(era5cds_struc(n)%w222) + + case( "neighbor" ) + deallocate(era5cds_struc(n)%n113) + + case( "average" ) + deallocate(era5cds_struc(n)%n111) + end select + enddo + deallocate(era5cds_struc) + +end subroutine finalize_era5cds diff --git a/ldt/metforcing/era5cds/get_era5cds.F90 b/ldt/metforcing/era5cds/get_era5cds.F90 new file mode 100644 index 0000000000..bbf4d0c419 --- /dev/null +++ b/ldt/metforcing/era5cds/get_era5cds.F90 @@ -0,0 +1,320 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! +! !ROUTINE: get_era5cds +! \label{get_era5cds} +! +! +! !REVISION HISTORY: +! 23 dec 2019: Sujay Kumar, initial code +! 17 apr 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: +subroutine get_era5cds(n, findex) +! !USES: + use LDT_coreMod + use LDT_timeMgrMod + use LDT_logMod + use LDT_constantsMod, only : LDT_CONST_PATH_LEN + use LDT_metforcingMod + use era5cds_forcingMod + + implicit none + +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: findex +! +! !DESCRIPTION: +! Opens, reads, and interpolates 1-hourly ERA5 forcing from +! the Climate Data Store. +! +! The ERA5 forcing data are organized into monthly files, where each +! file contains 24 one-hourly records over days in a month per instantaneous +! or accumulation forcing fields. Entire month of data is read in at the +! beginning of month. + +! In general, metforcing readers read the forcing data before the current +! time, referred to as bookend1, and after the current time, referred to as +! bookend2. Then the readers temporally interpolate between bookend1 and +! bookend2. +! +! +! The arguments are: +! \begin{description} +! \item[n] +! index of the nest +! \item[findex] +! forcing dataset index +! \end{description} +! +! The routines invoked are: +! \begin{description} +! \item[LDT\_tick](\ref{LDT_tick}) \newline +! call to advance or retract time +! \item[era5cdsfiles](\ref{era5cdsfiles}) \newline +! Puts together appropriate file name for 1 hour intervals +! \item[read\_era5cds](\ref{read_era5cds}) \newline +! call to read the ERA5 data and perform spatial interpolation +! \end{description} +!EOP + integer :: order + integer :: ferror + character(len=LDT_CONST_PATH_LEN) :: instfilename, lmlfilename + character(len=LDT_CONST_PATH_LEN) :: accfilename, prevaccfilename + integer :: c, r,kk,f,try + integer :: yr1, mo1, da1, hr1, mn1, ss1, doy1 + integer :: yr2, mo2, da2, hr2, mn2, ss2, doy2 + real*8 :: time1, time2, timenow + real*8 :: dtime1, dtime2 + real :: gmt1, gmt2 + real :: ts1, ts2 + + integer :: movetime ! Flag to move bookend2 files to bookend1 + logical :: retrieve_file + +! _________________________________________________________ + + if( LDT_rc%nts(n).gt.3600 ) then ! > 1-hr timestep + write(LDT_logunit,*) '[ERR] When running LDT with ERA5CDS, the clock ' + write(LDT_logunit,*) '[ERR] should run with a timestep less than or ' + write(LDT_logunit,*) '[ERR] equal to one hour.' + call LDT_endrun() + endif + if(LDT_rc%ts.gt.3600) then + write(LDT_logunit,*) '[ERR] The model timestep is > forcing data timestep' + write(LDT_logunit,*) '[ERR] LDT does not support this mode currently' + write(LDT_logunit,*) '[ERR] Program stopping ...' + call LDT_endrun() + endif + + era5cds_struc(n)%findtime1 = 0 + era5cds_struc(n)%findtime2 = 0 + movetime = 0 + retrieve_file = .false. + +!=== Determine Required ERA5 Data Times (The previous hour and the future hour) + yr1=LDT_rc%yr + mo1=LDT_rc%mo + da1=LDT_rc%da + hr1=LDT_rc%hr + mn1=LDT_rc%mn + ss1=0 + ts1=0 + call LDT_tick(timenow,doy1,gmt1,yr1,mo1,da1,hr1,mn1,ss1,ts1) + + if( timenow < era5cds_struc(n)%validstart ) then + write(LDT_logunit,*) '[ERR] ERA5CDS forecast begins at 7z on 1940/01/01' + write(LDT_logunit,*) '[ERR] Starting hour should be set to 7z ' + write(LDT_logunit,*) '[ERR] in lis.config file. Stopping ... ' + call LDT_endrun() + endif + + yr1 = LDT_rc%yr ! previous hour + mo1=LDT_rc%mo + da1=LDT_rc%da + hr1=LDT_rc%hr + mn1=0 + ss1=0 + ts1=0 + call LDT_tick(time1,doy1,gmt1,yr1,mo1,da1,hr1,mn1,ss1,ts1) + + yr2=LDT_rc%yr !next hour + mo2=LDT_rc%mo + da2=LDT_rc%da + hr2=LDT_rc%hr + mn2=0 + ss2=0 + ts2=60*60 + call LDT_tick(time2,doy2,gmt2,yr2,mo2,da2,hr2,mn2,ss2,ts2) + + ! First time step, read in data and set metdata1 + if(LDT_rc%tscount(n).eq.1 .or.LDT_rc%rstflag(n).eq.1 ) then + era5cds_struc(n)%findtime1=1 + era5cds_struc(n)%findtime2=1 + movetime=0 + LDT_rc%rstflag(n) = 0 + era5cds_struc(n)%mon = mo1 ! current month + endif + + ! Flag date/time when it is the switch of a month: + if( era5cds_struc(n)%mon .ne. mo2 ) then + era5cds_struc(n)%findtime1 = 1 + era5cds_struc(n)%mon = mo2 + + yr1 = yr2 + mo1 = mo2 + da1 = da2 + hr1 = hr2 + endif + + ! Flag date/time for when to open and read next ERA5 file: + if( LDT_get_nstep(LDT_rc,n) == 1 & + .or. (era5cds_struc(n)%findtime1 == 1 ) ) then + retrieve_file = .true. + era5cds_struc(n)%findtime1=0 + endif + + if(timenow.ge.era5cds_struc(n)%era5cdstime2) then + movetime = 1 + era5cds_struc(n)%findtime2 = 1 + endif + + if(movetime.eq.1) then + era5cds_struc(n)%era5cdstime1=era5cds_struc(n)%era5cdstime2 + do f=1,LDT_rc%met_nf(findex) + do c=1,LDT_rc%ngrid(n) + LDT_forc(n,findex)%metdata1(f,c)= LDT_forc(n,findex)%metdata2(f,c) + enddo + enddo + endif !end of movetime=1 + + if( retrieve_file ) then + !- Obtaining ERA5 File: + kk= 1 + if (LDT_get_nstep(LDT_rc,n) == 1) then + order = 1 + else + order = 2 + era5cds_struc(n)%findtime2 = 0 + endif + call era5cdsfiles(n,kk,findex,era5cds_struc(n)%era5cdsdir, & + yr1, mo1, da1, hr1, instfilename, accfilename, lmlfilename, & + prevaccfilename) + write(LDT_logunit,*) '[INFO] opening Bookend1 ',LDT_get_nstep(LDT_rc,n),hr1 + call read_era5cds(n, kk, order, yr1, mo1, da1, hr1, retrieve_file,& + findex, instfilename, accfilename, lmlfilename, & + prevaccfilename, ferror) + + if(ferror.ge.1) then !successfully retrieved forcing data + era5cds_struc(n)%era5cdstime1=time1 + retrieve_file = .false. + if(era5cds_struc(n)%findtime2==0) era5cds_struc(n)%era5cdstime2=time2 + else !ferror still=0 + write(LDT_logunit,*)'[ERR] ERA5CDS data missing file 1' + call LDT_endrun() + endif + endif + + if(era5cds_struc(n)%findtime2.eq.1) then + ! just need to assign metdata2 + kk= 1 + order = 2 + call era5cdsfiles(n,kk,findex,era5cds_struc(n)%era5cdsdir, & + yr2, mo2, da2, hr2, instfilename, accfilename, lmlfilename, & + prevaccfilename) + write(LDT_logunit,*) '[INFO] using Bookend2 ',LDT_get_nstep(LDT_rc,n),hr2 + call read_era5cds(n, kk, order, yr2, mo2, da2, hr2, retrieve_file, & + findex, instfilename, accfilename, lmlfilename, & + prevaccfilename, ferror) + + if(ferror.ge.1) then !successfully retrieved forcing data + era5cds_struc(n)%era5cdstime2=time2 + else + write(LDT_logunit,*)'[ERR] ERA5CDS data missing file 2' + call LDT_endrun() + endif + endif + +! should findex = 1 always? + if( LDT_rc%met_ecor(findex).eq."lapse-rate") then + call read_era5cdselev_ldtproc(n, findex, 0) + endif + +end subroutine get_era5cds + + +!BOP +! !ROUTINE: era5cdsfiles +! \label{era5cdsfiles} +! +! !INTERFACE: +subroutine era5cdsfiles(n, kk, findex, era5cdsdir, yr, mo, da, hr, & + instfilename, accfilename, lmlfilename, prevaccfilename) + +! !USES: + use LDT_coreMod + use LDT_logMod + use LDT_timeMgrMod + + implicit none +! !ARGUMENTS: + integer :: n + integer :: kk + integer :: findex + character(len=*), intent(in) :: era5cdsdir + integer, intent(in) :: yr,mo,da,hr + character(len=*), intent(out) :: instfilename, accfilename, lmlfilename + character(len=*), intent(out) :: prevaccfilename + +! !DESCRIPTION: +! This subroutine puts together ERA5 file names for +! daily netcdf files +! +! The arguments are: +! \begin{description} +! \item[era5cdsdir] +! Name of the ERA5 directory +! \item[yr] +! year +! \item[mo] +! month +! \item[da] +! day of month +! \item[fname] +! name of the timestamped ERA5 file +! \end{description} +! +!EOP + + character*4 :: cyear + character*2 :: cmonth + character*8 :: cdate + character*2 :: chour + integer :: mn, ss + real*8 :: time + integer :: doy + real :: gmt + integer :: iyr,imo,ida,ihr,imn,iss,ts,idoy + real :: igmt + real*8 :: itime + character :: ciyr*4, cimo*2, cida*2 + + mn = 0 + ss = 0 + + write(unit=cyear, fmt='(i4.4)') yr + write(unit=cmonth,fmt='(i2.2)') mo + write(unit=cdate,fmt='(i2.2)') da + write(unit=chour,fmt='(i2.2)') hr + + instfilename = trim(era5cdsdir)//'/'//cyear//'/era5_forcing_'//cyear//cmonth//'_instant.nc' + lmlfilename = trim(era5cdsdir)//'/'//cyear//'/era5_forcing_'//cyear//cmonth//'_ml_instant.nc' + accfilename = trim(era5cdsdir)//'/'//cyear//'/era5_forcing_'//cyear//cmonth//'_accum.nc' + ! accum file starts from 7z on the 1st of current month and ends on 6z on + ! the 1st of next month, and need previous months + !== roll back for last day of last month + iyr=yr; imo=mo; ida=da + ihr=hr; imn=0; iss=0 + ts = -24*60*60*da + call LDT_tick(itime,idoy,igmt,iyr,imo,ida,ihr,imn,iss,real(ts)) + write(ciyr, '(i4.4)') iyr + write(cimo, '(i2.2)') imo + write(cida, '(i2.2)') ida + if ( imo /= mo ) then + prevaccfilename = trim(era5cdsdir)//'/'//ciyr//'/era5_forcing_'//ciyr//cimo//'_accum.nc' + else + prevaccfilename = "none" + endif + +end subroutine era5cdsfiles + diff --git a/ldt/metforcing/era5cds/read_era5cds.F90 b/ldt/metforcing/era5cds/read_era5cds.F90 new file mode 100644 index 0000000000..04e0e795e7 --- /dev/null +++ b/ldt/metforcing/era5cds/read_era5cds.F90 @@ -0,0 +1,983 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LDT_misc.h" +!BOP +! +! !ROUTINE: read_era5cds +! \label{read_era5cds} +! +! !REVISION HISTORY: +! 23 dec 2019: Sujay Kumar, initial code +! 15 apr 2025: Hiroko Beaudoing, adopted ERA5 routines for the public CDS +! data format +! 19 jun 2026: Hiroko Beaudoing, updated to read in monthly data. +! Upscaling and Downscaling are not supported. +! +! !INTERFACE: +subroutine read_era5cds(n, kk,order, year, month, day, hour, read_flag, findex,& + instfile, avgfile, lmlfile, prevavgfile, ferror) +! !USES: + use LDT_coreMod, only : LDT_rc, LDT_domain, LDT_masterproc + use LDT_logMod + use LDT_FORC_AttributesMod + use LDT_metforcingMod, only : LDT_forc + use era5cds_forcingMod, only : era5cds_struc +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: kk + integer, intent(in) :: order + integer, intent(in) :: year + integer, intent(in) :: month + integer, intent(in) :: day + integer, intent(in) :: hour + logical, intent(in) :: read_flag + integer, intent(in) :: findex + character(len=*), intent(in) :: instfile, avgfile, lmlfile, prevavgfile + integer, intent(out) :: ferror + +! +! !DESCRIPTION: +! Mothly ERA5 data is read in at the beginning of simulation and/or month based +! on read_flag flag. For the given time, assigns appropriate forcing fields, +! transforms into 9 LIS forcing parameters and interpolates to the LIS domain. +! +! ERA5 model output variables used to force LIS are provided in 3 files: +! 1) inst, Instantaneous values, available every hour. +! 2) accum, Time integrated values (accumulations) over the past hour. +! +! ERA5 FORCING VARIABLES: +! 1. T2M inst, Atmospheric temperature (K) (2m) +! 2. D2M inst, Dewpoint temperature (K) (2m) +! 3. SSRD accum, surface solar radiation downwards [J m**-2] +! 4. STRD accum, surface thermal radiation downwards [J m**-2] +! 5. U10M inst, zonal wind [m/s] (10m) +! 6. V10M inst, meridional wind [m/s] (10m) +! 7. SP inst, surface pressure [Pa] +! 8. LSP accum, large scale precipitation [m] +! 9. CP accum, convective precipitation [m] +! ERA5 LOWEST MODEL LEVEL VARIABLES: +! 1. T inst, air temperature (K) +! 2. Q inst, specific_humidity (kg kg**-1) +! 5. U inst, zonal wind [m/s] +! 6. V inst, meridional wind [m/s] +! +! also inst file includes: +! Z: Geopotential (m**2 s**-2) +! +! NOTE 1: be aware that ECMWF outputs large-scale and convective precipitation +! separately. For total precipitation, need to sum the two fields, +! LSP+CP=TP. +! NOTE 2: SW & LW flux accumulations from time1 are interpolated via new zterp. +! NOTE 3: accum files starts at 7z on the 1st day of month and ends at 6z +! on the 1st day of next month. +! NOTE 4: Current month includes data up to 5 days ago. If the simulation end +! time is past avaiable data, forcing fields are assigned to the last +! available data with [WARN] message. LSM fields will be invalid. +! +! The arguments are: +! \begin{description} +! \item[order] +! flag indicating which data to be read (order=1, initial read, +! order=2, assign the next 1 hourly instance ) +! \item[n] +! index of the nest +! \item[name] +! name of the 1 hour ERA5 analysis file +! \item[tscount] +! time step count +! \item[ferror] +! return error code (1 indicates success) +! \end{description} +! +! The routines invoked are: +! \begin{description} +! \item[era5grid\_2\_lisgrid](\ref{era5grid_2_lisgrid}) \newline +! reverse south-north direction and shift 180 deg east-west +! \item[interp\_era5cds\_var](\ref{interp_era5cds_var}) \newline +! spatially interpolate the forcing data using interpolation of choice +! \item[assign\_processed\_era5cdsf](\ref{assign_processed_era5cdsf}) \newline +! assigns interpolated forcing data to the module data structure +! \end{description} +!EOP + +!--------------------------------------------------------------------------- + integer :: ftn + integer :: tmpId, qId, uwindId, vwindId, lwdId, psId, rainfId, crainfId + integer :: swdId, timeId + integer :: c,r,t,k,l,i,j,ll + integer :: tindex,atindex + integer :: mo,rec_size, prev_rec_size + integer :: start_time_index + logical :: file_exists + real :: missingValue + integer :: dimid, dim_len + integer :: cyear, cmonth + character(len=8) :: date + character(len=10) :: time + character(len=5) :: zone + integer, dimension(8) :: current + + real, allocatable :: tair(:,:,:) + real, allocatable :: qair(:,:,:) + real, allocatable :: swd(:,:,:) + real, allocatable :: lwd(:,:,:) + real, allocatable :: uwind(:,:,:) + real, allocatable :: vwind(:,:,:) + real, allocatable :: ps(:,:,:) + real, allocatable :: rainf(:,:,:) + real, allocatable :: crainf(:,:,:) + real, allocatable :: data4d(:,:,:,:) + real, allocatable :: var1(:,:,:) + + real, allocatable, dimension(:,:) :: datalis + + integer, parameter :: n_last_steps = 7 + integer :: days(12) + data days /31,28,31,30,31,30,31,31,30,31,30,31/ + character(len=*), parameter :: dim_name = "valid_time" + real, parameter :: epsln = 0.622 !Constants' values taken from + real, parameter :: A = 2.53E8, B=5.42E3 !Roger&Yau, A Short Course + !in Cloud Physics, pp.12-17 + ! A in kPa, B in K. + real :: p_kPa + + external :: era5grid_2_lisgrid + external :: interp_era5cds_var + external :: assign_processed_era5cdsff +! __________________________________________________________________________ + + ferror = 0 ! 1 success + tindex = (day - 1)*24 + hour + 1 + atindex = tindex - n_last_steps + + if(read_flag) then +#if (defined USE_NETCDF4) + + if((mod(year,4) .eq. 0 .and. mod(year, 100).ne.0) &!leap year + .or.(mod(year,400) .eq.0)) then + days(2) = 29 + else + days(2) = 28 + endif + + mo = LDT_rc%lnc(n)*LDT_rc%lnr(n) + rec_size = days(month)*24 + era5cds_struc(n)%tdimsize = rec_size + call date_and_time(date, time, zone, current) + cyear = current(1) + cmonth = current(2) + ! present month has incomplete data record + if (cmonth == month .and. cyear == year) then + call LDT_verify(nf90_open(path=trim(instfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + ! Get the ID for the dimension named "var_time" + call LDT_verify(nf90_inq_dimid(ftn, dim_name, dimid), & + 'nf90_inq_dimid failed in read_era5cds') + ! Get the length (size) of that dimension + call LDT_verify(nf90_inquire_dimension(ftn, dimid, len = dim_len), & + 'nf90_inquire_dimension failed in read_era5cds') + call LDT_verify(nf90_close(ftn), & + 'failed to close file in read_era5cds inst in lml') + rec_size = dim_len + write(LDT_logunit,*)'[INFO] in current month, reset rec_size =',rec_size + era5cds_struc(n)%tdimsize = rec_size + endif + + allocate(datalis(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold)) + !=== If using the lowest model level forcing (inst and mlinst): + !=== "t","q","u","v", and "sp " ==== + if(era5cds_struc(n)%uselml.eq.1) then + inquire(file=lmlfile,exist=file_exists) + if(file_exists) then + era5cds_struc(n)%ps = LDT_rc%udef + era5cds_struc(n)%tair = LDT_rc%udef + era5cds_struc(n)%qair = LDT_rc%udef + era5cds_struc(n)%uwind = LDT_rc%udef + era5cds_struc(n)%vwind = LDT_rc%udef + write(LDT_logunit,*)'[INFO] Reading ERA5 file (bookend,', order,' -',trim(instfile), ')' + + call LDT_verify(nf90_open(path=trim(instfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + + allocate(ps(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + + call LDT_verify(nf90_inq_varid(ftn,'sp',psId), & + 'nf90_inq_varid failed for psurf in read_era5cds') + call LDT_verify(nf90_get_var(ftn,psId, ps),& + 'nf90_get_var failed for ps in read_era5cds') + call LDT_verify(nf90_get_att(ftn,psId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for sp FillValue') + call LDT_verify(nf90_close(ftn), & + 'failed to close file in read_era5cds inst in lml') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,ps(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%ps(:,l)) + enddo + deallocate(ps) + + write(LDT_logunit,*)'[INFO] Reading ERA5 file (bookend,', order,' -',trim(lmlfile), ')' + call LDT_verify(nf90_open(path=trim(lmlfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + + ! model level data are in 4-dimentional + allocate(data4d(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,1,rec_size)) + call LDT_verify(nf90_inq_varid(ftn,'t',tmpId), & + 'nf90_inq_varid failed for t in read_era5cds') + call LDT_verify(nf90_get_var(ftn,tmpId, data4d),& + 'nf90_get_var failed for t in read_era5cds') + call LDT_verify(nf90_get_att(ftn,tmpId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for t FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,data4d(:,:,1,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%tair(:,l)) + enddo + + call LDT_verify(nf90_inq_varid(ftn,'q',qId), & + 'nf90_inq_varid failed for q in read_era5cds') + call LDT_verify(nf90_get_var(ftn,qId, data4d),& + 'nf90_get_var failed for q in read_era5cds') + call LDT_verify(nf90_get_att(ftn,qId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for q FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,data4d(:,:,1,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%qair(:,l)) + enddo + + call LDT_verify(nf90_inq_varid(ftn,'u',uwindId), & + 'nf90_inq_varid failed for u in read_era5cds') + call LDT_verify(nf90_get_var(ftn,uwindId, data4d),& + 'nf90_get_var failed for u in read_era5cds') + call LDT_verify(nf90_get_att(ftn,uwindId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for u FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,data4d(:,:,1,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%uwind(:,l)) + enddo + + call LDT_verify(nf90_inq_varid(ftn,'v',vwindId), & + 'nf90_inq_varid failed for v in read_era5cds') + call LDT_verify(nf90_get_var(ftn,vwindId, data4d),& + 'nf90_get_var failed for v in read_era5cds') + call LDT_verify(nf90_get_att(ftn,vwindId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for v FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,data4d(:,:,1,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%vwind(:,l)) + enddo + deallocate(data4d) + + call LDT_verify(nf90_close(ftn), & + 'failed to close file in read_era5cds lml') + ferror = 1 + else + write(LDT_logunit,*) '[ERR] ',trim(lmlfile)//' does not exist' + call LDT_endrun() + endif + + else !=== If using 2m and 10m forcing (inst): + !=== "t2m","d2m","u10","v10","sp " ==== + + inquire(file=instfile,exist=file_exists) + if(file_exists) then + era5cds_struc(n)%ps = LDT_rc%udef + era5cds_struc(n)%tair = LDT_rc%udef + era5cds_struc(n)%qair = LDT_rc%udef + era5cds_struc(n)%uwind = LDT_rc%udef + era5cds_struc(n)%vwind = LDT_rc%udef + + write(LDT_logunit,*)'[INFO] Reading ERA5 2m10m (bookend,', order,' -',trim(instfile), ')' + + call LDT_verify(nf90_open(path=trim(instfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + + allocate(ps(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + + call LDT_verify(nf90_inq_varid(ftn,'sp',psId), & + 'nf90_inq_varid failed for psurf in read_era5cds') + call LDT_verify(nf90_get_var(ftn,psId, ps),& + 'nf90_get_var failed for ps in read_era5cds') + call LDT_verify(nf90_get_att(ftn,psId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for sp FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,ps(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%ps(:,l)) + enddo + deallocate(ps) + allocate(qair(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LDT_verify(nf90_inq_varid(ftn,'d2m',qId), & + 'nf90_inq_varid failed for d2m in read_era5cds') + call LDT_verify(nf90_get_var(ftn,qId, qair),& + 'nf90_get_var failed for d2m in read_era5cds') + call LDT_verify(nf90_get_att(ftn,qId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for d2m FillValue') + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,qair(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%qair(:,l)) + enddo + deallocate(qair) + !----------------------------------------------------------------- + ! Calculate specific humidity from Dew point Temp. & Sfc Pressure. + ! Approximate: q~epsilon*e/p, p~p_sfc, e=e_s(Td), e_s=A*exp(-B/T) + !----------------------------------------------------------------- + do l=1,rec_size + do r=1,LDT_rc%lnr(n) + do c=1,LDT_rc%lnc(n) + k=c+(r-1)*LDT_rc%lnc(n) + if(era5cds_struc(n)%ps(k,l).ne.LDT_rc%udef.and.& + era5cds_struc(n)%qair(k,l).ne.LDT_rc%udef) then + p_kPa = era5cds_struc(n)%ps(k,l) / 1000.0 + era5cds_struc(n)%qair(k,l) = epsln*(A * EXP(-B/era5cds_struc(n)%qair(k,l))) / p_kPa + endif + enddo + enddo + enddo + + allocate(tair(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LDT_verify(nf90_inq_varid(ftn,'t2m',tmpId), & + 'nf90_inq_varid failed for t2m in read_era5cds') + call LDT_verify(nf90_get_var(ftn,tmpId, tair),& + 'nf90_get_var failed for t2m in read_era5cds') + call LDT_verify(nf90_get_att(ftn,tmpId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for t2m FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,tair(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%tair(:,l)) + enddo + deallocate(tair) + + allocate(uwind(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LDT_verify(nf90_inq_varid(ftn,'u10',uwindId), & + 'nf90_inq_varid failed for u10 in read_era5cds') + call LDT_verify(nf90_get_var(ftn,uwindId, uwind),& + 'nf90_get_var failed for u10 in read_era5cds') + call LDT_verify(nf90_get_att(ftn,uwindId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for u10 FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,uwind(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%uwind(:,l)) + enddo + deallocate(uwind) + + allocate(vwind(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LDT_verify(nf90_inq_varid(ftn,'v10',vwindId), & + 'nf90_inq_varid failed for v10 in read_era5cds') + call LDT_verify(nf90_get_var(ftn,vwindId, vwind),& + 'nf90_get_var failed for v10 in read_era5cds') + call LDT_verify(nf90_get_att(ftn,vwindId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for v10 FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,vwind(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%vwind(:,l)) + enddo + deallocate(vwind) + + call LDT_verify(nf90_close(ftn), & + 'failed to close file in read_era5cds inst') + ferror = 1 + else + write(LDT_logunit,*) '[ERR] ',trim(instfile)//' does not exist' + call LDT_endrun() + endif + + endif !=== era5cds_struc(n)%uselml + + !=== Read precipitation and radiation fields in accum file + !=== "ssrd","strd","lsp","cp" ==== + inquire(file=avgfile,exist=file_exists) + if(file_exists) then + era5cds_struc(n)%rainf = LDT_rc%udef + era5cds_struc(n)%crainf = LDT_rc%udef + era5cds_struc(n)%swd = LDT_rc%udef + era5cds_struc(n)%lwd = LDT_rc%udef + + write(LDT_logunit,*)'[INFO] Reading ERA5 file (bookend,', order,' -',trim(avgfile), ')' + call LDT_verify(nf90_open(path=trim(avgfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + + allocate(rainf(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LDT_verify(nf90_inq_varid(ftn,'lsp',rainfId), & + 'nf90_inq_varid failed for lsp in read_era5cds') + call LDT_verify(nf90_get_var(ftn,rainfId, rainf),& + 'nf90_get_var failed for lsp in read_era5cds') + call LDT_verify(nf90_get_att(ftn,rainfId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for lsp FillValue') + + ! large scale rainfall: lsp [m] -> mm/s + do l=1,rec_size + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(rainf(i,j,l).ne.missingValue) then + rainf(i,j,l) = rainf(i,j,l)*1000.0/(1.0*60*60) + endif + enddo + enddo + enddo + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,rainf(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.true.,& + era5cds_struc(n)%rainf(:,l)) + enddo + deallocate(rainf) + + allocate(crainf(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LDT_verify(nf90_inq_varid(ftn,'cp',crainfId), & + 'nf90_inq_varid failed for cp in read_era5cds') + call LDT_verify(nf90_get_var(ftn,crainfId, crainf),& + 'nf90_get_var failed for cp in read_era5cds') + call LDT_verify(nf90_get_att(ftn,crainfId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for cp FillValue') + + ! convective rainfall: cp [m] -> mm/s + do l=1,rec_size + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(crainf(i,j,l).ne.missingValue) then + crainf(i,j,l) = crainf(i,j,l)*1000.0/(1.0*60*60) + endif + enddo + enddo + enddo + + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,crainf(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.true.,& + era5cds_struc(n)%crainf(:,l)) + enddo + deallocate(crainf) + + !----------------------------------------------------------------- + ! rain = lsp + cp [mm/s] + !----------------------------------------------------------------- + do l=1,rec_size + do r=1,LDT_rc%lnr(n) + do c=1,LDT_rc%lnc(n) + k=c+(r-1)*LDT_rc%lnc(n) + if(era5cds_struc(n)%rainf(k,l).ne.LDT_rc%udef.and.& + era5cds_struc(n)%crainf(k,l).ne.LDT_rc%udef) then + era5cds_struc(n)%rainf(k,l) = era5cds_struc(n)%rainf(k,l) + & + era5cds_struc(n)%crainf(k,l) + endif + enddo + enddo + enddo + + allocate(swd(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LDT_verify(nf90_inq_varid(ftn,'ssrd',swdId), & + 'nf90_inq_varid failed for ssrd in read_era5cds') + call LDT_verify(nf90_get_var(ftn,swdId, swd),& + 'nf90_get_var failed for swd in read_era5cds') + call LDT_verify(nf90_get_att(ftn,swdId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for swd FillValue') + + ! swd - convert accumulated field to rate + do l=1,rec_size + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(swd(i,j,l).ne.missingValue) then + swd(i,j,l) = swd(i,j,l)/(1.0*60*60) + endif + enddo + enddo + enddo + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,swd(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%swd(:,l)) + enddo + deallocate(swd) + + allocate(lwd(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LDT_verify(nf90_inq_varid(ftn,'strd',lwdId), & + 'nf90_inq_varid failed for strd in read_era5cds') + call LDT_verify(nf90_get_var(ftn,lwdId, lwd),& + 'nf90_get_var failed for lwd in read_era5cds') + call LDT_verify(nf90_get_att(ftn,lwdId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for lwd FillValue') + + ! lwd - convert accumulated field to rate + do l=1,rec_size + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(lwd(i,j,l).ne.missingValue) then + lwd(i,j,l) = lwd(i,j,l)/(1.0*60*60) + endif + enddo + enddo + enddo + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,lwd(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%lwd(:,l)) + enddo + deallocate(lwd) + + call LDT_verify(nf90_close(ftn), & + 'failed to close avg file in read_era5cds') + ferror = 1 + else + write(LDT_logunit,*) '[ERR] ',trim(avgfile)//' does not exist' + call LDT_endrun() + + endif + !=== Read previous precipitation and radiation fields in accum file + !=== populate prev* arrays + inquire(file=prevavgfile,exist=file_exists) + if(file_exists) then + ! backward fill: read in last 7 time steps of prevavgfile + write(LDT_logunit,*)'[INFO] Reading prev ERA5 file (bookend,', order,' -',trim(prevavgfile), ')' + call LDT_verify(nf90_open(path=trim(prevavgfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + + call LDT_verify(nf90_inq_dimid(ftn,'valid_time',timeId), & + 'nf90_inq_dimid failed for timeId in read_era5cds') + call LDT_verify(nf90_inquire_dimension(ftn,timeId,len=prev_rec_size), & + 'nf90_inquire_dimension failed for timeId in read_era5cds') + start_time_index = prev_rec_size - n_last_steps + 1 + + allocate(var1(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,n_last_steps)) + call LDT_verify(nf90_inq_varid(ftn,'lsp',rainfId), & + 'nf90_inq_varid failed for lsp in read_era5cds') + call LDT_verify(nf90_get_var(ftn,rainfId, var1, & + start=(/1,1,start_time_index/), & + count=(/era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,n_last_steps/)),& + 'nf90_get_var failed for lsp in read_era5cds') + call LDT_verify(nf90_get_att(ftn,rainfId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for lsp FillValue') + + ! lsp [m] -> mm/s + do l=1,n_last_steps + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(var1(i,j,l).ne.missingValue) then + var1(i,j,l) = var1(i,j,l)*1000.0/(1.0*60*60) + endif + enddo + enddo + enddo + + do l=1,n_last_steps + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,var1(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.true.,& + era5cds_struc(n)%prev_rainf(:,l)) + enddo + + call LDT_verify(nf90_inq_varid(ftn,'cp',crainfId), & + 'nf90_inq_varid failed for cp in read_era5cds') + call LDT_verify(nf90_get_var(ftn,crainfId, var1, & + start=(/1,1,start_time_index/), & + count=(/era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,n_last_steps/)),& + 'nf90_get_var failed for cp in read_era5cds') + + ! cp [m] -> mm/s + do l=1,n_last_steps + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(var1(i,j,l).ne.missingValue) then + var1(i,j,l) = var1(i,j,l)*1000.0/(1.0*60*60) + endif + enddo + enddo + enddo + + do l=1,n_last_steps + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,var1(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.true.,& + era5cds_struc(n)%prev_crainf(:,l)) + enddo + + !----------------------------------------------------------------- + ! rain = lsp + cp [mm/s] + !----------------------------------------------------------------- + do l=1,n_last_steps + do r=1,LDT_rc%lnr(n) + do c=1,LDT_rc%lnc(n) + k=c+(r-1)*LDT_rc%lnc(n) + if(era5cds_struc(n)%prev_rainf(k,l).ne.LDT_rc%udef.and.& + era5cds_struc(n)%prev_crainf(k,l).ne.LDT_rc%udef) then + era5cds_struc(n)%prev_rainf(k,l) = & + era5cds_struc(n)%prev_rainf(k,l) + & + era5cds_struc(n)%prev_crainf(k,l) + endif + enddo + enddo + enddo + + call LDT_verify(nf90_inq_varid(ftn,'ssrd',swdId), & + 'nf90_inq_varid failed for ssrd in read_era5cds') + call LDT_verify(nf90_get_var(ftn,swdId, var1, & + start=(/1,1,start_time_index/), & + count=(/era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,n_last_steps/)),& + 'nf90_get_var failed for swd in read_era5cds') + call LDT_verify(nf90_get_att(ftn,swdId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for swd FillValue') + + ! swd - convert accumulated field to rate + do l=1,n_last_steps + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(var1(i,j,l).ne.missingValue) then + var1(i,j,l) = var1(i,j,l)/(1.0*60*60) + endif + enddo + enddo + enddo + do l=1,n_last_steps + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,var1(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%prev_swd(:,l)) + enddo + + call LDT_verify(nf90_inq_varid(ftn,'strd',lwdId), & + 'nf90_inq_varid failed for strd in read_era5cds') + call LDT_verify(nf90_get_var(ftn,lwdId, var1, & + start=(/1,1,start_time_index/), & + count=(/era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,n_last_steps/)),& + 'nf90_get_var failed for lwd in read_era5cds') + call LDT_verify(nf90_get_att(ftn,lwdId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for lwd FillValue') + + ! lwd - convert accumulated field to rate + do l=1,n_last_steps + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(var1(i,j,l).ne.missingValue) then + var1(i,j,l) = var1(i,j,l)/(1.0*60*60) + endif + enddo + enddo + enddo + do l=1,n_last_steps + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,var1(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.false.,& + era5cds_struc(n)%prev_lwd(:,l)) + enddo + deallocate(var1) + + call LDT_verify(nf90_close(ftn), & + 'failed to close avg file in read_era5cds') + else ! prevavgfile exists + ! at start of ERA5 in Jan 1940, no data for 0-6z + write(LDT_logunit,*) '[INFO] ',trim(prevavgfile)//' not available!' + do l = 1, n_last_steps + ll = rec_size - n_last_steps + l + era5cds_struc(n)%prev_rainf(:,l) = LDT_rc%udef + era5cds_struc(n)%prev_crainf(:,l) = LDT_rc%udef + era5cds_struc(n)%prev_swd(:,l) = LDT_rc%udef + era5cds_struc(n)%prev_lwd(:,l) = LDT_rc%udef + enddo + endif + + deallocate(datalis) +#endif + endif !if(read_flag) + + ! make sure current time is within valid data + if (era5cds_struc(n)%tdimsize < tindex) then + write(LDT_logunit,*) '[WARN] Current month passed available data, filling with last data' + tindex = era5cds_struc(n)%tdimsize + atindex = era5cds_struc(n)%tdimsize - n_last_steps + endif !(rec_size < tindex) + + call assign_processed_era5cdsf(n,kk,findex,order,1,era5cds_struc(n)%tair(:,tindex)) + call assign_processed_era5cdsf(n,kk,findex,order,2,era5cds_struc(n)%qair(:,tindex)) + call assign_processed_era5cdsf(n,kk,findex,order,5,era5cds_struc(n)%uwind(:,tindex)) + call assign_processed_era5cdsf(n,kk,findex,order,6,era5cds_struc(n)%vwind(:,tindex)) + call assign_processed_era5cdsf(n,kk,findex,order,7,era5cds_struc(n)%ps(:,tindex)) + if ( atindex .le. 0 ) then + atindex = atindex + n_last_steps + call assign_processed_era5cdsf(n,kk,findex,order,3,era5cds_struc(n)%prev_swd(:,atindex)) + call assign_processed_era5cdsf(n,kk,findex,order,4,era5cds_struc(n)%prev_lwd(:,atindex)) + call assign_processed_era5cdsf(n,kk,findex,order,8,era5cds_struc(n)%prev_rainf(:,atindex)) + call assign_processed_era5cdsf(n,kk,findex,order,9,era5cds_struc(n)%prev_crainf(:,atindex)) + else + call assign_processed_era5cdsf(n,kk,findex,order,3,era5cds_struc(n)%swd(:,atindex)) + call assign_processed_era5cdsf(n,kk,findex,order,4,era5cds_struc(n)%lwd(:,atindex)) + call assign_processed_era5cdsf(n,kk,findex,order,8,era5cds_struc(n)%rainf(:,atindex)) + call assign_processed_era5cdsf(n,kk,findex,order,9,era5cds_struc(n)%crainf(:,atindex)) + endif + if ( .not. read_flag .and. ferror == 0 ) ferror = 1 + +end subroutine read_era5cds + +!BOP +! +! !ROUTINE: interp_era5cds_var +! \label{interp_era5cds_var} +! +! !INTERFACE: +subroutine interp_era5cds_var(n,findex, month, input_var, missingValue, & + pcp_flag, output_var) + +! !USES: + use LDT_coreMod + use LDT_logMod + use era5cds_forcingMod, only : era5cds_struc +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + implicit none + +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: findex + integer, intent(in) :: month + real, intent(in) :: input_var(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold) + real, intent(in) :: missingValue + logical, intent(in) :: pcp_flag + real, intent(out) :: output_var(LDT_rc%lnc(n)*LDT_rc%lnr(n)) + +! +! !DESCRIPTION: +! This subroutine spatially interpolates a ERA5CDS field +! to the LDT running domain +! +!EOP + + integer :: t,c,r,k,iret + integer :: doy + integer :: ftn + real :: f (era5cds_struc(n)%ncold*era5cds_struc(n)%nrold) + logical*1 :: lb(era5cds_struc(n)%mi) + logical*1 :: lo(LDT_rc%lnc(n)*LDT_rc%lnr(n)) + integer :: input_size + integer :: input_nc, input_nr + integer :: count1,nrec +! _____________________________________________________________ + + input_size = era5cds_struc(n)%mi + output_var = LDT_rc%udef + lo = .true. + input_nc = era5cds_struc(n)%ncold + input_nr = era5cds_struc(n)%nrold + + lb = .false. + do r=1,era5cds_struc(n)%nrold + do c=1,era5cds_struc(n)%ncold + k= c+(r-1)*era5cds_struc(n)%ncold + f(k) = input_var(c,r) + if(f(k) .ne. missingValue) then + lb(k) = .true. + endif + enddo + enddo + + if(pcp_flag.and.& + trim(LDT_rc%met_gridtransform(findex)).eq."budget-bilinear") then + + call conserv_interp(LDT_rc%gridDesc(n,:),lb,f,lo,& + output_var, & + era5cds_struc(n)%mi,LDT_rc%lnc(n)*LDT_rc%lnr(n),& + LDT_domain(n)%lat, LDT_domain(n)%lon,& + era5cds_struc(n)%w112,era5cds_struc(n)%w122,& + era5cds_struc(n)%w212,era5cds_struc(n)%w222,& + era5cds_struc(n)%n112,era5cds_struc(n)%n122,& + era5cds_struc(n)%n212,era5cds_struc(n)%n222,& + LDT_rc%udef, iret) + + elseif(trim(LDT_rc%met_gridtransform(findex)).eq."bilinear".or.& + trim(LDT_rc%met_gridtransform(findex)).eq."budget-bilinear") then + + call bilinear_interp(LDT_rc%gridDesc(n,:),lb,f,lo,& + output_var, & + era5cds_struc(n)%mi,LDT_rc%lnc(n)*LDT_rc%lnr(n), & + LDT_domain(n)%lat, LDT_domain(n)%lon,& + era5cds_struc(n)%w111,era5cds_struc(n)%w121,& + era5cds_struc(n)%w211,era5cds_struc(n)%w221,& + era5cds_struc(n)%n111,era5cds_struc(n)%n121,& + era5cds_struc(n)%n211,era5cds_struc(n)%n221,& + LDT_rc%udef, iret) + + elseif(trim(LDT_rc%met_gridtransform(findex)).eq."neighbor") then + call neighbor_interp(LDT_rc%gridDesc(n,:),lb,f,lo,& + output_var,era5cds_struc(n)%mi,& + LDT_rc%lnc(n)*LDT_rc%lnr(n),& + LDT_domain(n)%lat, LDT_domain(n)%lon,& + era5cds_struc(n)%n113,LDT_rc%udef,iret) + else + write(LDT_logunit,*) '[ERR] Spatial interpolation option '//& + trim(LDT_rc%met_gridtransform(findex))//& + ' not supported for ERA5CDS' + call LDT_endrun() + endif + + +end subroutine interp_era5cds_var +!BOP +! +! !ROUTINE: era5grid_2_lisgrid +! +! !DESCRIPTION: +! Changes grid_data from ECMWF data convention to GLDAS convention +! +! ERA5: North-to-South around Greenwich Meridian +! Global grid. Data are written as flat binary from "upper left to +! lower right" starting at 0.25-degree grid point center coordinates: +! 0.125E,89.875N and going to 0.125W,89.875S. +! +! LIS: South-to-North around Date Line +! Full global grid. Starts at the southernmost latitude and date line, +! going east and then north. +! +! !REVISION HISTORY: +! 10 Apr 2002: Urszula Jambor; Code adapted from +! ecmwfgrid_2_grid2catgrid, by R. Reichle +! 21 Apr 2025: Hiroko Beaudoing; modified berggrid_2_gldasgrid.F90 +! to work for odd ny. +! +! !INTERFACE: +subroutine era5grid_2_lisgrid( nx, ny, grid_data, out_data ) +!EOP + implicit none + + integer, intent(in) :: nx, ny + real, intent(in), dimension(nx,ny) :: grid_data + real, intent(out), dimension(nx*ny):: out_data + real, dimension(nx,ny):: lis_data + + integer :: i, j, m, n, c + real :: tmp, tmp_data1(nx) + + ! ------------------------------------------------------------------ + ! some checks + + if ((nx /= 1440) .or. (ny /= 720)) then + write (*,*) 'era5grid_2_gldasgrid(): This routine has only been' + write (*,*) 'checked for nx=1440 and ny=720. STOPPING' + stop + end if + if ((mod(nx,2) /= 0)) then + write (*,*) 'era5grid_2_gldasgrid(): This routine can only work' + write (*,*) 'for even nx. STOPPING.' + stop + end if + + !------------------------------------------------------------------- + + do j=1,ny + + ! swap latitude bands (North-to-South becomes South-to-North) + n = ny-j+1 + tmp_data1 = grid_data(:,j) + + do i=1,nx/2 + + ! shift longitudes (wrapping around Greenwhich Meridian becomes + ! wrapping around Date Line) + m = i + nx/2 + tmp = tmp_data1(i) + tmp_data1(i) = tmp_data1(m) + tmp_data1(m) = tmp + + end do + + lis_data(:,n) = tmp_data1 + + end do + + ! return output in 1D + c=0 + do j=1,ny + do i=1,nx + c=c+1 + out_data(c) = lis_data(i,j) + end do + end do + +end subroutine era5grid_2_lisgrid + +!BOP +! +! !ROUTINE: assign_processed_era5cdsf +! \label{assign_processed_era5cdsf} +! +! !INTERFACE: +subroutine assign_processed_era5cdsf(n,kk,findex,order,var_index,era5data) +! !USES: + use LDT_coreMod, only : LDT_rc, LDT_domain + use LDT_metforcingMod, only : LDT_forc +! +! !DESCRIPTION: +! This routine assigns the interpolated ERA5 forcing data +! to the module data structures to be used later for +! time interpolation +! +!EOP + implicit none + + integer, intent(in) :: n + integer, intent(in) :: kk + integer, intent(in) :: findex + integer, intent(in) :: order + integer, intent(in) :: var_index + real, intent(in) :: era5data(LDT_rc%lnc(n)*LDT_rc%lnr(n)) + + + integer :: c,r + + do r=1,LDT_rc%lnr(n) + do c=1,LDT_rc%lnc(n) + if(LDT_domain(n)%gindex(c,r).ne.-1) then + if(order.eq.1) then + LDT_forc(n,findex)%metdata1(var_index,& + LDT_domain(n)%gindex(c,r)) = & + era5data(c+(r-1)*LDT_rc%lnc(n)) + elseif(order.eq.2) then + LDT_forc(n,findex)%metdata2(var_index,& + LDT_domain(n)%gindex(c,r)) = & + era5data(c+(r-1)*LDT_rc%lnc(n)) + endif + endif + enddo + enddo +end subroutine assign_processed_era5cdsf diff --git a/ldt/metforcing/era5cds/read_era5cds_elev.F90 b/ldt/metforcing/era5cds/read_era5cds_elev.F90 new file mode 100644 index 0000000000..5f1764a057 --- /dev/null +++ b/ldt/metforcing/era5cds/read_era5cds_elev.F90 @@ -0,0 +1,150 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LDT_misc.h" +!BOP +! !ROUTINE: read_era5cds_elev +! \label{read_era5cds_elev} +! +! !REVISION HISTORY: +! +! 23 dec 2019: Sujay Kumar; Added ERA5 terrain height reader +! 17 apr 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: +subroutine read_era5cds_elev( n, findex, era5cdselev, elevdiff ) + +! !USES: + use LDT_coreMod, only : LDT_rc, LDT_domain + use LDT_metforcingMod, only : LDT_forc + use era5cds_forcingMod, only : era5cds_struc + use LDT_logMod, only : LDT_logunit, LDT_verify, & + LDT_endrun + use LDT_fileIOMod, only : LDT_transform_paramgrid +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none + +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: findex +!- Terrain height will be set to run domain: + real, intent(inout) :: era5cdselev(LDT_rc%lnc(n),LDT_rc%lnr(n),1) + real, intent(inout) :: elevdiff(LDT_rc%met_nc(findex), LDT_rc%met_nr(findex)) + +! !DESCRIPTION: +! +! Opens, reads, and interpolates ERA5 model elevation to the LDT +! grid. The data will be used to perform any topographical +! adjustments to the forcing. +! +! The arguments are: +! \begin{description} +! \item[n] +! index of the nest +! \item[findex] +! index of the forcing dataset selected +! \end{description} +! +! The routines invoked are: +! \begin{description} +! \end{description} +!EOP + logical :: file_exists + integer :: ftn_const + integer :: i,c,r,k,iret + + integer :: phisId + real :: read_phis(era5cds_struc(n)%ncold, era5cds_struc(n)%nrold) + + ! Grid transform fields: + integer :: inpts, outpts + real :: elev1d(era5cds_struc(n)%ncold*era5cds_struc(n)%nrold) + logical*1 :: lb(era5cds_struc(n)%ncold*era5cds_struc(n)%nrold) + real :: elev_regrid(LDT_rc%lnc(n)*LDT_rc%lnr(n)) + logical*1 :: lb_regrid(LDT_rc%lnc(n)*LDT_rc%lnr(n)) + +! _____________________________________________________________________________ + +#if (defined USE_NETCDF3) + write(LDT_logunit,*) "ERR: ERA5CDS terrain height reader requires NetCDF4" + call LDT_endrun() +#endif + + era5cdselev = LDT_rc%udef + elevdiff = LDT_rc%udef + + inquire(file = trim(era5cds_struc(n)%era5cdshgt_file), exist=file_exists) + if(.not. file_exists) then + write(LDT_logunit,*) "[ERR] The ERA5CDS terrain height file ",& + trim(era5cds_struc(n)%era5cdshgt_file)," is not found." + write(LDT_logunit,*) "Program stopping ..." + call LDT_endrun + endif + +! ------------------------------------------------------------------- +! Open and Read-in Forcing Terrain Hght File - Bring to LIS run domain +! ------------------------------------------------------------------- + write(LDT_logunit,*) "[INFO] Reading the ERA5CDS terrain height file: ", & + trim(era5cds_struc(n)%era5cdshgt_file) + +#if (defined USE_NETCDF4) + + ! Open the ERA5 surface elevation field (PHIS): + call LDT_verify(nf90_open(path=trim(era5cds_struc(n)%era5cdshgt_file), & + mode=NF90_NOWRITE, ncid=ftn_const), & + 'nf90_open failed in read_era5cds_elev') + + call LDT_verify(nf90_inq_varid(ftn_const,'elev',phisId), & + 'nf90_inq_varid failed for elev in read_era5cds_elev') + + ! Reading in PHIS field: + call LDT_verify(nf90_get_var(ftn_const, phisId, read_phis), & + 'nf90_get_var failed for elev in read_era5cds_elev') + call LDT_verify(nf90_close(ftn_const), & + 'nf90_close failed for elev in read_era5cds_elev') + + ! Initialize arrays for grid transformation: + inpts = LDT_rc%met_nc(findex)*LDT_rc%met_nr(findex) + outpts = LDT_rc%lnr(n)*LDT_rc%lnc(n) + lb = .true. + lb_regrid = .true. + elev1d = -9999.0 + + ! Convert 2D to 1D array for the interplation call: + do r = 1, LDT_rc%met_nr(findex) + do c = 1, LDT_rc%met_nc(findex) + k= c+(r-1)*LDT_rc%met_nc(findex) + elev1d(k) = read_phis(c,r) + if ( elev1d(k) == -9999.0 ) then + elev1d(k) = LDT_rc%udef + lb(k) = .false. + endif + enddo + enddo + + ! Interp elevation field to output field: + call LDT_transform_paramgrid(n, LDT_rc%met_gridtransform_parms(findex), & + LDT_rc%met_gridDesc(findex,:), inpts, 1, elev1d, lb, & + outpts, elev_regrid, lb_regrid ) + + i = 0 + do r = 1, LDT_rc%lnr(n) + do c = 1, LDT_rc%lnc(n) + i = i + 1 + era5cdselev(c,r,1) = elev_regrid(i) + end do + end do + + +#endif +end subroutine read_era5cds_elev diff --git a/ldt/metforcing/era5cds/read_era5cdselev_ldtproc.F90 b/ldt/metforcing/era5cds/read_era5cdselev_ldtproc.F90 new file mode 100644 index 0000000000..a357d81d9c --- /dev/null +++ b/ldt/metforcing/era5cds/read_era5cdselev_ldtproc.F90 @@ -0,0 +1,70 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.8 +! +! Copyright (c) 2026 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! !ROUTINE: read_era5cdselev_ldtproc +! \label{read_era5cdselev_ldtproc} +! +! !REVISION HISTORY: +! +! 17Dec2004; Sujay Kumar; Initial Specificaton +! +! !INTERFACE: +subroutine read_era5cdselev_ldtproc(n, findex, change) +! !USES: + use LDT_constantsMod, only : LDT_CONST_PATH_LEN + use LDT_coreMod, only : LDT_rc, LDT_domain + use LDT_metforcingMod, only : LDT_forc + use LDT_fileIOMod, only : LDT_read_param + use LDT_logMod, only : LDT_logunit, LDT_endrun + use era5cds_forcingMod, only : era5cds_struc + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: findex + integer, intent(in) :: change + +! !DESCRIPTION: +! +! Opens, reads, and interpolates ERA5CDS model elevation to the LDT +! grid. The data will be used to perform any topographical +! adjustments to the forcing. +! +! The arguments are: +! \begin{description} +! \item[n] +! index of the nest +! \end{description} +!EOP + integer :: c,r,line1,line2,nc_dom,line + integer :: glnc, glnr + real :: go(LDT_rc%lnc(n),LDT_rc%lnr(n)) + + if ( trim(LDT_rc%met_ecor(findex)).eq."lapse-rate")then + write(LDT_logunit,*) " Reading the ERA5CDS elevation: " + if ( change == 0 ) then ! constant for all period + call LDT_read_param(n,"ELEV_ERA5CDS",go) + else + write(LDT_logunit,*) 'ERR: invalid update request. ', change, & + 'is not for ERA5CDS elev' + call LDT_endrun + endif + + do r=1,LDT_rc%lnr(n) + do c=1,LDT_rc%lnc(n) + if(LDT_domain(n)%gindex(c,r).ne.-1) then + LDT_forc(n,findex)%modelelev(LDT_domain(n)%gindex(c,r)) = go(c,r) + endif + enddo + enddo + + endif + +end subroutine read_era5cdselev_ldtproc diff --git a/ldt/metforcing/era5cds/readcrd_era5cds.F90 b/ldt/metforcing/era5cds/readcrd_era5cds.F90 new file mode 100644 index 0000000000..9cfab595e4 --- /dev/null +++ b/ldt/metforcing/era5cds/readcrd_era5cds.F90 @@ -0,0 +1,72 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! +! !ROUTINE: readcrd_era5cds +! \label{readcrd_era5cds} +! +! !REVISION HISTORY: +! 23 Dec 2019: Sujay Kumar, initial code +! 17 Apr 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: +subroutine readcrd_era5cds() +! !USES: + use ESMF + use LDT_coreMod, only : LDT_rc, LDT_config + use LDT_logMod + use era5cds_forcingMod, only : era5cds_struc +! +! !DESCRIPTION: +! +! This routine reads the options specific to ERA5 forcing +! (downloaded from the Climate Data) from the LDT configuration file. +! +!EOP + implicit none + + integer :: n,t,rc + + call ESMF_ConfigFindLabel(LDT_config,"ERA5CDS forcing directory:",rc=rc) + do n=1,LDT_rc%nnest + call ESMF_ConfigGetAttribute(LDT_config,era5cds_struc(n)%era5cdsdir,& + rc=rc) + call LDT_verify(rc,& + 'ERA5CDS forcing directory: not defined') + enddo + + call ESMF_ConfigFindLabel(LDT_config,"ERA5CDS forcing terrain height file:",rc=rc) + do n=1,LDT_rc%nnest + call ESMF_ConfigGetAttribute(LDT_config,era5cds_struc(n)%era5cdshgt_file,& + rc=rc) + call LDT_verify(rc,& + 'ERA5CDS forcing terrain height file: not defined') + enddo + + call ESMF_ConfigFindLabel(LDT_config,"ERA5CDS use lowest model level forcing:",rc=rc) + do n=1,LDT_rc%nnest + call ESMF_ConfigGetAttribute(LDT_config,era5cds_struc(n)%uselml,rc=rc) + call LDT_verify(rc,& + 'ERA5CDS use lowest model level forcing: not defined') + enddo + + do n=1,LDT_rc%nnest + write(LDT_logunit,*) '[INFO] Using ERA5CDS forcing' + write(LDT_logunit,*) '[INFO] ERA5CDS forcing directory: ',& + trim(era5cds_struc(n)%era5cdsDIR) + write(LDT_logunit,*) '[INFO] ERA5CDS use lowest model level forcing: ',& + era5cds_struc(n)%uselml + + era5cds_struc(n)%era5cdstime1 = 3000.0 + era5cds_struc(n)%era5cdstime2 = 0.0 + + enddo +end subroutine readcrd_era5cds diff --git a/ldt/metforcing/era5cds/reset_era5cds.F90 b/ldt/metforcing/era5cds/reset_era5cds.F90 new file mode 100644 index 0000000000..0501dd8a5d --- /dev/null +++ b/ldt/metforcing/era5cds/reset_era5cds.F90 @@ -0,0 +1,41 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! !MODULE: reset_era5cds +! \label{reset_era5cds} +! +! !REVISION HISTORY: +! 23 dec 2019: Sujay Kumar, initial code +! 17 apr 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: +subroutine reset_era5cds +! !USES: + use LDT_coreMod, only : LDT_rc + use LDT_timeMgrMod, only : LDT_date2time + use era5cds_forcingMod +! +! !DESCRIPTION: +! Routine to cleanup allocated structures for era5cds forcing. +! +!EOP + implicit none + integer :: n + + do n=1,LDT_rc%nnest + era5cds_struc(n)%startFlag = .true. + era5cds_struc(n)%dayFlag = .true. + era5cds_struc(n)%era5cdstime1 = 3000.0 + era5cds_struc(n)%era5cdstime2 = 0.0 + era5cds_struc(n)%ringtime = 0.0 + era5cds_struc(n)%reset_flag = .true. + enddo +end subroutine reset_era5cds diff --git a/ldt/metforcing/era5cds/timeinterp_era5cds.F90 b/ldt/metforcing/era5cds/timeinterp_era5cds.F90 new file mode 100644 index 0000000000..824a4e3c6b --- /dev/null +++ b/ldt/metforcing/era5cds/timeinterp_era5cds.F90 @@ -0,0 +1,227 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! !ROUTINE: timeinterp_era5cds +! \label{timeinterp_era5cds} +! +! !REVISION HISTORY: +! 23 dec 2019: Sujay Kumar, initial code +! 17 Apr 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: + +subroutine timeinterp_era5cds(n,findex) + +! !USES: + use ESMF + use LDT_coreMod, only : LDT_rc, LDT_domain, LDT_localPet + use LDT_metforcingMod, only : LDT_forc, LDT_FORC_Base_State + use LDT_FORC_AttributesMod + use LDT_constantsMod, only : LDT_CONST_SOLAR + use LDT_timeMgrMod + use LDT_logMod, only : LDT_logunit, LDT_verify, LDT_endrun + use era5cds_forcingMod, only : era5cds_struc + + implicit none + +! !ARGUMENTS: + integer, intent(in):: n + integer, intent(in):: findex +! +! !DESCRIPTION: +! Temporally interpolates the forcing data to the current model +! timestep. Downward shortwave radiation is interpolated using a +! zenith-angled based approach. Precipitation and downward longwave +! are not temporally interpolated, and the previous 1 hourly value +! is used. All other variables are linearly interpolated between +! the 1 hourly blocks. +! +! The routines invoked are: +! \begin{description} +! \item[LDT\_time2date](\ref{LDT_time2date}) \newline +! converts the time to a date format +! \item[zterp](\ref{zterp}) \newline +! zenith-angle based interpolation +! \end{description} +!EOP + integer :: t,zdoy,k,kk + integer :: index1 + real :: wt1,wt2,czb,cze,czm,gmt1,gmt2 + real :: zw1,zw2 + integer :: bdoy,byr,bmo,bda,bhr,bmn + real*8 :: btime + integer :: status + type(ESMF_Field) :: tmpField,q2Field,uField,vField,swdField,lwdField + type(ESMF_Field) :: psurfField,pcpField,cpcpField + real,pointer :: tmp(:),q2(:),uwind(:),vwind(:) + real,pointer :: swd(:),lwd(:),psurf(:),pcp(:),cpcp(:) + +!----------------------------------------------------------------------- +! Interpolate Data in Time +!----------------------------------------------------------------------- + + btime=era5cds_struc(n)%era5cdstime1 + call LDT_time2date(btime,bdoy,gmt1,byr,bmo,bda,bhr,bmn) + btime=era5cds_struc(n)%era5cdstime2 + call LDT_time2date(btime,bdoy,gmt2,byr,bmo,bda,bhr,bmn) + + wt1=(era5cds_struc(n)%era5cdstime2-LDT_rc%time)/ & + (era5cds_struc(n)%era5cdstime2-era5cds_struc(n)%era5cdstime1) + wt2=1.0-wt1 + + call ESMF_StateGet(LDT_FORC_Base_State(n,findex),LDT_FORC_Tair%varname(1),tmpField,& + rc=status) + call LDT_verify(status, 'Error: Enable Tair in the forcing variables list') + + call ESMF_StateGet(LDT_FORC_Base_State(n,findex),LDT_FORC_Qair%varname(1),q2Field,& + rc=status) + call LDT_verify(status, 'Error: Enable Qair in the forcing variables list') + + call ESMF_StateGet(LDT_FORC_Base_State(n,findex),LDT_FORC_SWdown%varname(1),swdField,& + rc=status) + call LDT_verify(status, 'Error: Enable SWdown in the forcing variables list') + + call ESMF_StateGet(LDT_FORC_Base_State(n,findex),LDT_FORC_LWdown%varname(1),lwdField,& + rc=status) + call LDT_verify(status, 'Error: Enable LWdown in the forcing variables list') + + call ESMF_StateGet(LDT_FORC_Base_State(n,findex),LDT_FORC_Wind_E%varname(1),uField,& + rc=status) + call LDT_verify(status, 'Error: Enable Wind_E in the forcing variables list') + + call ESMF_StateGet(LDT_FORC_Base_State(n,findex),LDT_FORC_Wind_N%varname(1),vField,& + rc=status) + call LDT_verify(status, 'Error: Enable Wind_N in the forcing variables list') + + call ESMF_StateGet(LDT_FORC_Base_State(n,findex),LDT_FORC_Psurf%varname(1),psurfField,& + rc=status) + call LDT_verify(status, 'Error: Enable Psurf in the forcing variables list') + + call ESMF_StateGet(LDT_FORC_Base_State(n,findex),LDT_FORC_Rainf%varname(1),pcpField,& + rc=status) + call LDT_verify(status, 'Error: Enable Rainf in the forcing variables list') + + call ESMF_StateGet(LDT_FORC_Base_State(n,findex),LDT_FORC_CRainf%varname(1),cpcpField,& + rc=status) + call LDT_verify(status, 'Error: Enable CRainf in the forcing variables list') + + call ESMF_FieldGet(tmpField,localDE=0,farrayPtr=tmp,rc=status) + call LDT_verify(status) + + call ESMF_FieldGet(q2Field,localDE=0,farrayPtr=q2,rc=status) + call LDT_verify(status) + + call ESMF_FieldGet(swdField,localDE=0,farrayPtr=swd,rc=status) + call LDT_verify(status) + + call ESMF_FieldGet(lwdField,localDE=0,farrayPtr=lwd,rc=status) + call LDT_verify(status) + + call ESMF_FieldGet(uField,localDE=0,farrayPtr=uwind,rc=status) + call LDT_verify(status) + + call ESMF_FieldGet(vField,localDE=0,farrayPtr=vwind,rc=status) + call LDT_verify(status) + + call ESMF_FieldGet(psurfField,localDE=0,farrayPtr=psurf,rc=status) + call LDT_verify(status) + + call ESMF_FieldGet(pcpField,localDE=0,farrayPtr=pcp,rc=status) + call LDT_verify(status) + + call ESMF_FieldGet(cpcpField,localDE=0,farrayPtr=cpcp,rc=status) + call LDT_verify(status) + +!----------------------------------------------------------------------- +! SW down +!----------------------------------------------------------------------- + do t=1,LDT_rc%ntiles(n) + index1 = LDT_domain(n)%tile(t)%index + zdoy=LDT_rc%doy + call zterp(0,LDT_domain(n)%grid(index1)%lat,& + LDT_domain(n)%grid(index1)%lon,& + gmt1,gmt2,LDT_rc%gmt,zdoy,zw1,zw2,czb,cze,czm,LDT_rc) + +! past hour average + if (LDT_forc(n,findex)%metdata1(3,index1).ne.LDT_rc%udef) then + + swd(t) = zw1 * LDT_forc(n,findex)%metdata1(3,index1) + + if (swd(t).gt.LDT_CONST_SOLAR) then + write(unit=LDT_logunit,fmt=*)'[WARN] sw radiation too high in ERA5!!' + write(unit=LDT_logunit,fmt=*)'[WARN] it is',swd(t), zw1 + write(unit=LDT_logunit,fmt=*)'[WARN] at ',t,LDT_domain(n)%grid(index1)%lat,LDT_domain(n)%grid(index1)%lon + write(unit=LDT_logunit,fmt=*)'[WARN] era5cdsdata1=',LDT_forc(n,findex)%metdata1(3,index1) + write(unit=LDT_logunit,fmt=*)'[WARN] era5cdsdata2=',LDT_forc(n,findex)%metdata2(3,index1) + swd(t) = LDT_CONST_SOLAR + write(unit=LDT_logunit,fmt=*)'[WARN] forcing set to ',swd(t) + endif + endif + + if ((swd(t).ne.LDT_rc%udef).and.(swd(t).lt.0)) then + ! tiny negative swd appear in the high latitudes sometimes + ! WARN in stead of ERR and endrun -- 2/2/2026 HKB + swd(t) = 0.0 + write(LDT_logunit,*) & + '[WARN] timeinterp_era5cds -- forcing not udef but lt0' + write(LDT_logunit,*)'[WARN] timeinterp_era5cds -- ', & + t,swd(t),LDT_forc(n,findex)%metdata2(3,index1), & + ' (',LDT_localPet,')' + !call LDT_endrun + endif + + enddo + +!----------------------------------------------------------------------- +! precip variable - constant rate over the ERA5 hour +!----------------------------------------------------------------------- + do t=1,LDT_rc%ntiles(n) + index1 = LDT_domain(n)%tile(t)%index + pcp(t) = LDT_forc(n,findex)%metdata1(8,index1) + cpcp(t) = LDT_forc(n,findex)%metdata1(9,index1) + + if ( pcp(t) < 0 ) then + pcp(t) = 0 + endif + if ( cpcp(t) < 0 ) then + cpcp(t) = 0 + endif + enddo + +!----------------------------------------------------------------------- +! LW down +!----------------------------------------------------------------------- + do t=1,LDT_rc%ntiles(n) + index1 = LDT_domain(n)%tile(t)%index + lwd(t) = LDT_forc(n,findex)%metdata1(4,index1) + enddo + +!----------------------------------------------------------------------- +! Linearly interpolate everything else +!----------------------------------------------------------------------- + + do t=1,LDT_rc%ntiles(n) + index1 = LDT_domain(n)%tile(t)%index + + tmp(t) = LDT_forc(n,findex)%metdata1(1,index1)*wt1 + & + LDT_forc(n,findex)%metdata2(1,index1)*wt2 + q2(t) = LDT_forc(n,findex)%metdata1(2,index1)*wt1 + & + LDT_forc(n,findex)%metdata2(2,index1)*wt2 + uwind(t) = LDT_forc(n,findex)%metdata1(5,index1)*wt1+& + LDT_forc(n,findex)%metdata2(5,index1)*wt2 + vwind(t) = LDT_forc(n,findex)%metdata1(6,index1)*wt1+& + LDT_forc(n,findex)%metdata2(6,index1)*wt2 + psurf(t) = LDT_forc(n,findex)%metdata1(7,index1)*wt1 + & + LDT_forc(n,findex)%metdata2(7,index1)*wt2 + enddo + +end subroutine timeinterp_era5cds + diff --git a/ldt/plugins/LDT_metforcing_pluginMod.F90 b/ldt/plugins/LDT_metforcing_pluginMod.F90 index c320eb6658..5762406166 100644 --- a/ldt/plugins/LDT_metforcing_pluginMod.F90 +++ b/ldt/plugins/LDT_metforcing_pluginMod.F90 @@ -101,6 +101,7 @@ subroutine LDT_metforcing_plugin use gfs_forcingMod use merra2_forcingMod use era5_forcingMod + use era5cds_forcingMod use gswp1_forcingMod use gswp2_forcingMod #if ( defined MF_AGRMET ) @@ -197,6 +198,11 @@ subroutine LDT_metforcing_plugin external finalize_era5 external reset_era5 + external get_era5cds + external timeinterp_era5cds + external finalize_era5cds + external reset_era5cds + external get_agrradps external timeinterp_agrradps external finalize_agrradps @@ -398,6 +404,13 @@ subroutine LDT_metforcing_plugin call registerresetmetforc(trim(LDT_ERA5Id)//char(0),reset_ERA5) call registerfinalmetforc(trim(LDT_ERA5Id)//char(0),finalize_ERA5) +! - ERA5 Reanalysis Forcing from Climate Data Store: + call registerinitmetforc(trim(LDT_era5cdsId)//char(0),init_era5cds) + call registerretrievemetforc(trim(LDT_era5cdsId)//char(0),get_era5cds) + call registertimeinterpmetforc(trim(LDT_era5cdsId)//char(0),timeinterp_era5cds) + call registerresetmetforc(trim(LDT_era5cdsId)//char(0),reset_era5cds) + call registerfinalmetforc(trim(LDT_era5cdsId)//char(0),finalize_era5cds) + ! - WRFv2 Analysis Forcing: call registerinitmetforc(trim(LDT_wrfoutv2Id)//char(0),init_WRFoutv2) diff --git a/ldt/plugins/LDT_param_pluginMod.F90 b/ldt/plugins/LDT_param_pluginMod.F90 index 41bd922bd9..3f4a8012f3 100644 --- a/ldt/plugins/LDT_param_pluginMod.F90 +++ b/ldt/plugins/LDT_param_pluginMod.F90 @@ -1293,6 +1293,7 @@ subroutine LDT_forcingparams_plugin external read_ecmwf_elev external read_merra2_elev external read_era5_elev + external read_era5cds_elev external read_wrfoutv2_elev external read_wrfak_elev ! external read_geos5_elev @@ -1328,6 +1329,10 @@ subroutine LDT_forcingparams_plugin call registerreadforcelev(trim(LDT_era5Id)//char(0),& read_era5_elev) +!- ERA5 Climate Data Store forcing: + call registerreadforcelev(trim(LDT_era5cdsId)//char(0),& + read_era5cds_elev) + !- WRFoutv2 forcing: call registerreadforcelev(trim(LDT_wrfoutv2Id)//char(0),& read_WRFoutv2_elev) diff --git a/ldt/plugins/LDT_pluginIndices.F90 b/ldt/plugins/LDT_pluginIndices.F90 index 9119c37025..0545b1246a 100644 --- a/ldt/plugins/LDT_pluginIndices.F90 +++ b/ldt/plugins/LDT_pluginIndices.F90 @@ -149,6 +149,7 @@ module LDT_pluginIndices character*50, public, parameter :: LDT_princetonId = "PRINCETON" character*50, public, parameter :: LDT_merra2Id = "MERRA2" character*50, public, parameter :: LDT_era5Id = "ERA5" + character*50, public, parameter :: LDT_era5cdsId = "ERA5CDS" character*50, public, parameter :: LDT_gswp1Id = "GSWP1" character*50, public, parameter :: LDT_gswp2Id = "GSWP2" character*50, public, parameter :: LDT_nldas2Id = "NLDAS2" diff --git a/lis/configs/lis.config.adoc b/lis/configs/lis.config.adoc index d419f28cd6..bd07c34b2d 100644 --- a/lis/configs/lis.config.adoc +++ b/lis/configs/lis.config.adoc @@ -391,6 +391,7 @@ endif::devonly[] |"`AWAP`" | AWAP precipitation data |"`GDAS T1534`" | NCEP-specific GDAS T1534 forcing data |"`ERA5`" | ERA5 reanalysis meteorology +|"`ERA5CDS`" | ERA5 reanalysis meteorology from Climate Data Store |"`AWRAL processed forcing`" | AWRA-L processed forcing |"`PLUMBER2`" | PLUMBER2 forcing |"`GDDP`" | NEX GDDP forcing @@ -6259,6 +6260,31 @@ ERA5 forcing directory: ./ERA5/ ERA5 forcing tile to grid mapping file: ../ERA5/mapping.nc .... +[[sssec_forcings_ERA5CDS,ERA5CDS]] +==== ERA5CDS + +`ERA5CDS forcing directory:` specifies the location of +the ERA5 from Climate Data Store forcing files. + + +`ERA5CDS surface altitude file:` specifies the ERA5CDS reanalysis elevation file. + +`ERA5CDS use lowest model level forcing:` specifies whether to use the +model lowest level or 2m/10m diagnosed wind fields. The 2m/10m file is used +if this model level forcing option is turned off. +Acceptable values are: + +|==== +|Value | Description + +|0 | Use 2m temperature, 2m dew point temperature, and 10m u/v wind +|1 | Use model lowest level temperature, specific humidity, and u/v wind +|==== +.Example _lis.config_ entry +.... +ERA5CDS forcing directory: ./ERA5/single_hourly +ERA5CDS surface altitude file: ../ERA5/era5cds_elev.nc +ERA5CDS use lowest model level forcing: 1 [[sssec_forcings_gswp1,GSWP1]] ==== GSWP1 diff --git a/lis/interp/zterp.F90 b/lis/interp/zterp.F90 index 996f24fd46..f2d1f7556f 100644 --- a/lis/interp/zterp.F90 +++ b/lis/interp/zterp.F90 @@ -33,6 +33,9 @@ ! less than 50% of 37 intervals (arbitary) are day-time. This ! modification is needed to properly interpolate ECMWF data, ! where daily peak is too low due to the nature of dataset. +! 3/4/26 Hiroko Beaudoing; Additional fix for iflag=0 case, when weight1 +! becomes huge due to extremely small avgangle, resulting in +! unrealistically large SWdown flux at the instance. ! ! !INTERFACE: subroutine zterp (iflag,lat,lon,btime,etime, & @@ -131,6 +134,10 @@ subroutine zterp (iflag,lat,lon,btime,etime, & integer zoneetime,zonebtime,zonembtime integer zone,juliane,julianmb,juliantemp,i integer acnt ! counter for cavgdata > 0 (daytime) + + external :: LIS_localtime + external :: coszenith + !-------------------------------------------------------------------- ! This section contains hardwired data that will be supplied by main program. ! These values were chosen arbitrarily and exist simply to check the @@ -229,6 +236,11 @@ subroutine zterp (iflag,lat,lon,btime,etime, & else weight1=(czmodel/avgangle) endif + if ( lis%zterp_correction ) then + if (weight1.gt.5.0) then + weight1=1.0 + endif + endif endif !-------------------------------------------------------------------- ! Second case: IFLAG=1 (instantaneous input and output) @@ -420,7 +432,7 @@ subroutine coszenith (lon,latd,lhour,zone,julian,czenith,dec,omega) czenith=sin(dec)*sin(llat)+cos(dec)*cos(llat)*cos(omega) czenith=amax1(0.,czenith) if ( czenith > 1.0 ) then - write(LIS_logunit,fmt='(a,e16.10,a)') 'czenith > 1 (',czenith,& + write(LIS_logunit,fmt='(a,e17.10,a)') 'czenith > 1 (',czenith,& '). Resetting to 1.' czenith = 1.0 endif diff --git a/lis/make/default.cfg b/lis/make/default.cfg index d368ef7de7..4376a0ae6e 100755 --- a/lis/make/default.cfg +++ b/lis/make/default.cfg @@ -297,6 +297,11 @@ enabled: True macro: MF_ERA5 path: metforcing/era5 +[ERA5CDS] +enabled: True +macro: MF_ERA5CDS +path: metforcing/era5cds + [CMAP] enabled: True macro: MF_CMAP diff --git a/lis/metforcing/era5cds/era5cds_forcingMod.F90 b/lis/metforcing/era5cds/era5cds_forcingMod.F90 new file mode 100644 index 0000000000..04c35c0cf3 --- /dev/null +++ b/lis/metforcing/era5cds/era5cds_forcingMod.F90 @@ -0,0 +1,433 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LIS_misc.h" +module era5cds_forcingMod +!BOP +! !MODULE: era5cds_forcingMod +! +! !DESCRIPTION: +! This module contains variables and data structures that are used +! for the implementation of the ERA5 forcing data from Climate Data Store. +! The data is global 0.25 degree dataset in latlon +! projection, and at 1 hourly intervals. The derived +! data type {\tt era5cds\_struc} +! includes the variables that specify the runtime options, and the +! weights and neighbor information to be used for spatial interpolation. +! They are described below: +! \begin{description} +! \item[ncold] +! Number of columns (along the east west dimension) for the input data +! \item[nrold] +! Number of rows (along the north south dimension) for the input data +! \item[nmif] +! Number of forcing variables in the ECMWF data +! \item[era5cdstime1] +! The nearest, previous 1 hour instance of the incoming +! data (as a real time). +! \item[era5cdstime2] +! The nearest, next 1 hour instance of the incoming +! data (as a real time). +! \item[era5cdsdir] +! Directory containing the input data +! \item[mi] +! Number of points in the input grid +! \item[n111,n121,n211,n221] +! Arrays containing the neighbor information of the input grid +! for each grid point in LIS, for bilinear interpolation. +! \item[w111,w121,w211,w221] +! Arrays containing the weights of the input grid +! for each grid point in LIS, for bilinear interpolation. +! \item[n122,n122,n212,n222] +! Arrays containing the neighbor information of the input grid +! for each grid point in LIS, for conservative interpolation. +! \item[w112,w122,w212,w222] +! Arrays containing the weights of the input grid +! for each grid point in LIS, for conservative interpolation. +! \item[n113] +! Arrays containing the neighbor information of the input grid +! for each grid point in LIS, for n. neighbor interpolation. +! \item[findtime1, findtime2] +! boolean flags to indicate which time is to be read for +! temporal interpolation. +! \end{description} +! +! !USES: + + use LIS_constantsMod, only : LIS_CONST_PATH_LEN + + implicit none + + PRIVATE +!----------------------------------------------------------------------------- +! !PUBLIC MEMBER FUNCTIONS: +!----------------------------------------------------------------------------- + public :: init_era5cds !defines the native resolution of + !the input data +!----------------------------------------------------------------------------- +! !PUBLIC TYPES: +!----------------------------------------------------------------------------- + public :: era5cds_struc + +!EOP + type, public :: era5cds_type_dec + + real :: ts + integer :: ncold, nrold + character(len=LIS_CONST_PATH_LEN) :: era5cdsdir !ERA5CDS Forcing Directory + character(len=LIS_CONST_PATH_LEN) :: era5cdsalt_file ! ERA5CDS altitude file + real*8 :: era5cdstime1,era5cdstime2 + character*50 :: met_interp + logical :: reset_flag + integer :: mon + + integer :: mi + integer, allocatable :: n111(:) + integer, allocatable :: n121(:) + integer, allocatable :: n211(:) + integer, allocatable :: n221(:) + real, allocatable :: w111(:),w121(:) + real, allocatable :: w211(:),w221(:) + + integer, allocatable :: n112(:,:) + integer, allocatable :: n122(:,:) + integer, allocatable :: n212(:,:) + integer, allocatable :: n222(:,:) + real, allocatable :: w112(:,:),w122(:,:) + real, allocatable :: w212(:,:),w222(:,:) + integer, allocatable :: n113(:) + ! Subset parameters for "none" reprojection type: + integer, allocatable :: lat_line(:,:) + integer, allocatable :: lon_line(:,:) + real :: subset_gridDesc(20) + integer :: subset_nc, subset_nr + + integer :: findtime1, findtime2 + logical :: startFlag, dayFlag + real*8 :: validstart + + real, allocatable :: tair(:,:) + real, allocatable :: qair(:,:) + real, allocatable :: uwind(:,:) + real, allocatable :: vwind(:,:) + real, allocatable :: ps(:,:) + real, allocatable :: rainf(:,:) + real, allocatable :: crainf(:,:) + real, allocatable :: swd(:,:) + real, allocatable :: lwd(:,:) + + real, allocatable :: prev_rainf(:,:) + real, allocatable :: prev_crainf(:,:) + real, allocatable :: prev_swd(:,:) + real, allocatable :: prev_lwd(:,:) + + integer :: nvars + integer :: uselml + integer :: tdimsize + + real*8 :: ringtime + + integer :: nIter, st_iterid,en_iterid + + real, allocatable :: metdata1(:,:,:) + real, allocatable :: metdata2(:,:,:) + + end type era5cds_type_dec + + type(era5cds_type_dec), allocatable :: era5cds_struc(:) + +contains + +!BOP +! +! !ROUTINE: init_era5cds +! \label{init_era5cds} +! +! !REVISION HISTORY: +! 23 Dec 2019: Sujay Kumar, initial code +! 04 Mar 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: + subroutine init_era5cds(findex) + +! !USES: + use LIS_coreMod + use LIS_timeMgrMod + use LIS_logMod + use LIS_spatialDownscalingMod, only : LIS_init_pcpclimo_native + use LIS_forecastMod + use LIS_gridmappingMod, only : LIS_RunDomainPts +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none +! !AGRUMENTS: + integer, intent(in) :: findex +! +! !DESCRIPTION: +! Defines the native resolution of the input forcing for ERA5CDS +! data. The grid description arrays are based on the decoding +! schemes used by NCEP and followed in the LIS interpolation +! schemes (see Section~\ref{interp}). +! +! The routines invoked are: +! \begin{description} +! \item[readcrd\_era5cds](\ref{readcrd_era5cds}) \newline +! reads the runtime options specified for ERA5CDS data +! \item[bilinear\_interp\_input](\ref{bilinear_interp_input}) \newline +! computes the neighbor, weights for bilinear interpolation +! \item[conserv\_interp\_input](\ref{conserv_interp_input}) \newline +! computes the neighbor, weights for conservative interpolation +! \end{description} +!EOP + real :: gridDesci(LIS_rc%nnest,50) + integer :: updoy, yr1,mo1,da1,hr1,mn1,ss1 + real :: upgmt + integer :: n + integer :: glpnc, glpnr + + external :: readcrd_era5cds + external :: bilinear_interp_input + external :: conserv_interp_input + external :: neighbor_interp_input + external :: upscaleByAveraging_input + external :: read_era5cds_elev + + allocate(era5cds_struc(LIS_rc%nnest)) + + do n=1, LIS_rc%nnest + + era5cds_struc(n)%ncold = 1440 + era5cds_struc(n)%nrold = 720 + era5cds_struc(n)%mon = -1 + + allocate(era5cds_struc(n)%tair(LIS_rc%lnc(n)*LIS_rc%lnr(n),745)) + allocate(era5cds_struc(n)%qair(LIS_rc%lnc(n)*LIS_rc%lnr(n),745)) + allocate(era5cds_struc(n)%uwind(LIS_rc%lnc(n)*LIS_rc%lnr(n),745)) + allocate(era5cds_struc(n)%vwind(LIS_rc%lnc(n)*LIS_rc%lnr(n),745)) + allocate(era5cds_struc(n)%ps(LIS_rc%lnc(n)*LIS_rc%lnr(n),745)) + allocate(era5cds_struc(n)%rainf(LIS_rc%lnc(n)*LIS_rc%lnr(n),745)) + allocate(era5cds_struc(n)%crainf(LIS_rc%lnc(n)*LIS_rc%lnr(n),745)) + allocate(era5cds_struc(n)%swd(LIS_rc%lnc(n)*LIS_rc%lnr(n),745)) + allocate(era5cds_struc(n)%lwd(LIS_rc%lnc(n)*LIS_rc%lnr(n),745)) + + ! save accumulation fields + ! previous month's 17-23Z on last day = 0-6Z on 1st next month + allocate(era5cds_struc(n)%prev_rainf(LIS_rc%lnc(n)*LIS_rc%lnr(n),7)) + allocate(era5cds_struc(n)%prev_crainf(LIS_rc%lnc(n)*LIS_rc%lnr(n),7)) + allocate(era5cds_struc(n)%prev_swd(LIS_rc%lnc(n)*LIS_rc%lnr(n),7)) + allocate(era5cds_struc(n)%prev_lwd(LIS_rc%lnc(n)*LIS_rc%lnr(n),7)) + + enddo + + call readcrd_era5cds() + LIS_rc%met_nf(findex) = 9 + + era5cds_struc%reset_flag = .false. + + do n=1, LIS_rc%nnest + era5cds_struc(n)%ts = 3600 !hourly + call LIS_update_timestep(LIS_rc, n, era5cds_struc(n)%ts) + enddo + + gridDesci = 0 + LIS_rc%met_proj(findex) = "latlon" + + do n=1,LIS_rc%nnest + gridDesci(n,1) = 0 + gridDesci(n,2) = era5cds_struc(n)%ncold + gridDesci(n,3) = era5cds_struc(n)%nrold + gridDesci(n,4) = -89.875 + gridDesci(n,5) = -179.875 + gridDesci(n,6) = 128 + gridDesci(n,7) = 89.875 + gridDesci(n,8) = 179.875 + gridDesci(n,9) = 0.25 + gridDesci(n,10) = 0.25 + gridDesci(n,20) = 0 + + ! ERA5 accumulation data starts at 7z on 1940-01-01 + yr1 = 1940 + mo1 = 01 + da1 = 01 + hr1 = 07 + mn1 = 0; ss1 = 0 + call LIS_date2time( era5cds_struc(n)%validstart,updoy,upgmt, & + yr1,mo1,da1,hr1,mn1,ss1 ) + + era5cds_struc(n)%mi = era5cds_struc(n)%ncold*era5cds_struc(n)%nrold + + ! Check resolution and set up weights for Interpolation + if ( LIS_isatAfinerResolution(n,gridDesci(n,9)) ) then + + era5cds_struc(n)%met_interp = LIS_rc%met_interp(findex) + + write(LIS_logunit,*) '[INFO] The ERA5CDS forcing resolution is ' // & + ' coarser than the running domain.' + write(LIS_logunit,*) ' Interpolating with the ' // & + trim(era5cds_struc(n)%met_interp) // ' method.' + ! Setting up weights for Interpolation + if(trim(LIS_rc%met_interp(findex)).eq."bilinear") then + allocate(era5cds_struc(n)%n111(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%n121(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%n211(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%n221(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%w111(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%w121(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%w211(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%w221(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + call bilinear_interp_input(n, gridDesci(n,:),& + era5cds_struc(n)%n111,era5cds_struc(n)%n121,& + era5cds_struc(n)%n211,era5cds_struc(n)%n221,& + era5cds_struc(n)%w111,era5cds_struc(n)%w121,& + era5cds_struc(n)%w211,era5cds_struc(n)%w221) + + elseif(trim(LIS_rc%met_interp(findex)).eq."budget-bilinear") then + allocate(era5cds_struc(n)%n111(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%n121(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%n211(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%n221(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%w111(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%w121(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%w211(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + allocate(era5cds_struc(n)%w221(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + call bilinear_interp_input(n, gridDesci(n,:),& + era5cds_struc(n)%n111,era5cds_struc(n)%n121,& + era5cds_struc(n)%n211,era5cds_struc(n)%n221,& + era5cds_struc(n)%w111,era5cds_struc(n)%w121,& + era5cds_struc(n)%w211,era5cds_struc(n)%w221) + + allocate(era5cds_struc(n)%n112(LIS_rc%lnc(n)*LIS_rc%lnr(n),25)) + allocate(era5cds_struc(n)%n122(LIS_rc%lnc(n)*LIS_rc%lnr(n),25)) + allocate(era5cds_struc(n)%n212(LIS_rc%lnc(n)*LIS_rc%lnr(n),25)) + allocate(era5cds_struc(n)%n222(LIS_rc%lnc(n)*LIS_rc%lnr(n),25)) + allocate(era5cds_struc(n)%w112(LIS_rc%lnc(n)*LIS_rc%lnr(n),25)) + allocate(era5cds_struc(n)%w122(LIS_rc%lnc(n)*LIS_rc%lnr(n),25)) + allocate(era5cds_struc(n)%w212(LIS_rc%lnc(n)*LIS_rc%lnr(n),25)) + allocate(era5cds_struc(n)%w222(LIS_rc%lnc(n)*LIS_rc%lnr(n),25)) + call conserv_interp_input(n, gridDesci(n,:),& + era5cds_struc(n)%n112,era5cds_struc(n)%n122,& + era5cds_struc(n)%n212,era5cds_struc(n)%n222,& + era5cds_struc(n)%w112,era5cds_struc(n)%w122,& + era5cds_struc(n)%w212,era5cds_struc(n)%w222) + + elseif(trim(LIS_rc%met_interp(findex)).eq."neighbor") then + allocate(era5cds_struc(n)%n113(LIS_rc%lnc(n)*LIS_rc%lnr(n))) + call neighbor_interp_input(n, gridDesci(n,:),& + era5cds_struc(n)%n113) + + endif + ! Running domain at 0.25 degree, no need for interpolation + else if (gridDesci(n,9) == LIS_rc%gridDesc(n,9) .and. & + gridDesci(n,10) == LIS_rc%gridDesc(n,10).and. & + LIS_rc%gridDesc(n,1) == 0 ) then + write(LIS_logunit,*) '[INFO] ERA5CDS and LIS resolutions match ' // & + 'no spatial trainsorm of input forcing (only subset may occur). ' + era5cds_struc(n)%met_interp = "none" + call LIS_RunDomainPts( n, LIS_rc%met_proj(findex), gridDesci(n,:), & + glpnc, glpnr, era5cds_struc(n)%subset_nc, & + era5cds_struc(n)%subset_nr, era5cds_struc(n)%subset_gridDesc, & + era5cds_struc(n)%lat_line, era5cds_struc(n)%lon_line ) + + else + era5cds_struc(n)%met_interp = LIS_rc%met_upscale(findex) + write(LIS_logunit,*) '[INFO] The ERA5CDS forcing resolution is finer ' // & + 'than the running domain.' + write(LIS_logunit,*) ' Upscaling with the ' // & + trim(era5cds_struc(n)%met_interp) // ' method.' + + select case( era5cds_struc(n)%met_interp ) + case( "average" ) + allocate(era5cds_struc(n)%n111(era5cds_struc(n)%mi)) + + call upscaleByAveraging_input(gridDesci, & + LIS_rc%gridDesc(n,:), & + era5cds_struc(n)%mi, & + LIS_rc%lnc(n)*LIS_rc%lnr(n), & + era5cds_struc(n)%n111) + case default + write(LIS_logunit,*) '[ERR] Interpolation option '// & + trim(LIS_rc%met_interp(findex))//& + ' for ERA5CDS forcing is not supported' + call LIS_endrun() + end select + endif + + call LIS_registerAlarm("ERA5CDS forcing alarm",& + 86400.0,86400.0) + era5cds_struc(n)%startFlag = .true. + era5cds_struc(n)%dayFlag = .true. + + era5cds_struc(n)%nvars = 9 + + ! Forecast mode: + if(LIS_rc%forecastMode.eq.1) then + + if(mod(LIS_rc%nensem(n),& + LIS_forecast_struc(1)%niterations).ne.0) then + write(LIS_logunit,*) & + '[ERR] The number of ensembles must be a multiple' + write(LIS_logunit,*) '[ERR] of the number of iterations ' + write(LIS_logunit,*) '[ERR] nensem = ',LIS_rc%nensem(n) + write(LIS_logunit,*) '[ERR] niter = ', & + LIS_forecast_struc(1)%niterations + call LIS_endrun() + endif + + era5cds_struc(n)%st_iterid = LIS_forecast_struc(1)%st_iterId + era5cds_struc(n)%en_iterId = LIS_forecast_struc(1)%niterations + era5cds_struc(n)%nIter = LIS_forecast_struc(1)%niterations + + allocate(era5cds_struc(n)%metdata1(LIS_forecast_struc(1)%niterations,& + LIS_rc%met_nf(findex),& + LIS_rc%ngrid(n))) + allocate(era5cds_struc(n)%metdata2(LIS_forecast_struc(1)%niterations,& + LIS_rc%met_nf(findex),& + LIS_rc%ngrid(n))) + + ! Regular retrospective or non-forecast mode: + else + + era5cds_struc(n)%st_iterid = 1 + era5cds_struc(n)%en_iterId = 1 + era5cds_struc(n)%nIter = 1 + + allocate(era5cds_struc(n)%metdata1(1,LIS_rc%met_nf(findex),& + LIS_rc%ngrid(n))) + allocate(era5cds_struc(n)%metdata2(1,LIS_rc%met_nf(findex),& + LIS_rc%ngrid(n))) + + endif + + era5cds_struc(n)%metdata1 = 0 + era5cds_struc(n)%metdata2 = 0 + + era5cds_struc(n)%findtime1 = 0 + era5cds_struc(n)%findtime2 = 0 + + ! Set up precipitation climate downscaling: + if(LIS_rc%pcp_downscale(findex).ne.0) then + call LIS_init_pcpclimo_native(n,findex,& + era5cds_struc(n)%ncold,& + era5cds_struc(n)%nrold) + endif + + if ( LIS_rc%met_ecor(findex) == "lapse-rate" .or. & + LIS_rc%met_ecor(findex) == "lapse-rate and slope-aspect" ) then + + call read_era5cds_elev(n,findex) + endif + + enddo ! End nest loop + + end subroutine init_era5cds +end module era5cds_forcingMod + diff --git a/lis/metforcing/era5cds/finalize_era5cds.F90 b/lis/metforcing/era5cds/finalize_era5cds.F90 new file mode 100644 index 0000000000..98322a6137 --- /dev/null +++ b/lis/metforcing/era5cds/finalize_era5cds.F90 @@ -0,0 +1,77 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! !ROUTINE: finalize_era5cds +! \label{finalize_era5cds} +! +! !REVISION HISTORY: +! 23 Dec 2019: Sujay Kumar, initial code +! 16 Apr 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: + +subroutine finalize_era5cds(findex) + + ! !USES: + use LIS_coreMod, only : LIS_rc + use era5cds_forcingMod, only : era5cds_struc +! +! !DESCRIPTION: +! Routine to cleanup ERA5 forcing related memory allocations. +! +!EOP + implicit none + + integer :: findex + integer :: n + + do n=1,LIS_rc%nnest + select case( LIS_rc%met_interp(findex) ) + case( "bilinear" ) + if (allocated(era5cds_struc(n)%n111)) deallocate(era5cds_struc(n)%n111) + if (allocated(era5cds_struc(n)%n121)) deallocate(era5cds_struc(n)%n121) + if (allocated(era5cds_struc(n)%n211)) deallocate(era5cds_struc(n)%n211) + if (allocated(era5cds_struc(n)%n221)) deallocate(era5cds_struc(n)%n221) + if (allocated(era5cds_struc(n)%w111)) deallocate(era5cds_struc(n)%w111) + if (allocated(era5cds_struc(n)%w121)) deallocate(era5cds_struc(n)%w121) + if (allocated(era5cds_struc(n)%w211)) deallocate(era5cds_struc(n)%w211) + if (allocated(era5cds_struc(n)%w221)) deallocate(era5cds_struc(n)%w221) + + case( "budget-bilinear" ) + if (allocated(era5cds_struc(n)%n111)) deallocate(era5cds_struc(n)%n111) + if (allocated(era5cds_struc(n)%n121)) deallocate(era5cds_struc(n)%n121) + if (allocated(era5cds_struc(n)%n211)) deallocate(era5cds_struc(n)%n211) + if (allocated(era5cds_struc(n)%n221)) deallocate(era5cds_struc(n)%n221) + if (allocated(era5cds_struc(n)%w111)) deallocate(era5cds_struc(n)%w111) + if (allocated(era5cds_struc(n)%w121)) deallocate(era5cds_struc(n)%w121) + if (allocated(era5cds_struc(n)%w211)) deallocate(era5cds_struc(n)%w211) + if (allocated(era5cds_struc(n)%w221)) deallocate(era5cds_struc(n)%w221) + if (allocated(era5cds_struc(n)%n112)) deallocate(era5cds_struc(n)%n112) + if (allocated(era5cds_struc(n)%n122)) deallocate(era5cds_struc(n)%n122) + if (allocated(era5cds_struc(n)%n212)) deallocate(era5cds_struc(n)%n212) + if (allocated(era5cds_struc(n)%n222)) deallocate(era5cds_struc(n)%n222) + if (allocated(era5cds_struc(n)%w112)) deallocate(era5cds_struc(n)%w112) + if (allocated(era5cds_struc(n)%w122)) deallocate(era5cds_struc(n)%w122) + if (allocated(era5cds_struc(n)%w212)) deallocate(era5cds_struc(n)%w212) + if (allocated(era5cds_struc(n)%w222)) deallocate(era5cds_struc(n)%w222) + + case( "neighbor" ) + if (allocated(era5cds_struc(n)%n113)) deallocate(era5cds_struc(n)%n113) + + case( "average" ) + if (allocated(era5cds_struc(n)%n111)) deallocate(era5cds_struc(n)%n111) + end select + if (allocated(era5cds_struc(n)%lat_line)) deallocate(era5cds_struc(n)%lat_line) + if (allocated(era5cds_struc(n)%lon_line)) deallocate(era5cds_struc(n)%lon_line) + enddo + deallocate(era5cds_struc) + +end subroutine finalize_era5cds diff --git a/lis/metforcing/era5cds/get_era5cds.F90 b/lis/metforcing/era5cds/get_era5cds.F90 new file mode 100644 index 0000000000..6199614fe6 --- /dev/null +++ b/lis/metforcing/era5cds/get_era5cds.F90 @@ -0,0 +1,354 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! +! !ROUTINE: get_era5cds +! \label{get_era5cds} +! +! +! !REVISION HISTORY: +! 23 dec 2019: Sujay Kumar, initial code +! 04 Mar 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: +subroutine get_era5cds(n, findex) +! !USES: + use LIS_coreMod + use LIS_timeMgrMod + use LIS_logMod + use LIS_metforcingMod + use era5cds_forcingMod + use LIS_constantsMod, only: LIS_CONST_PATH_LEN + + implicit none + +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: findex +! +! !DESCRIPTION: +! Opens, reads, and interpolates 1-hourly ERA5 forcing from +! the Climate Data Store. +! +! The ERA5 forcing data are organized into monthly files, where each +! file contains 24 one-hourly records over days in a month per instantaneous +! or accumulation forcing fields. Entire month of data is read in at the +! beginning of month. +! +! In general, metforcing readers read the forcing data before the current +! time, referred to as bookend1, and after the current time, referred to as +! bookend2. Then the readers temporally interpolate between bookend1 and +! bookend2. +! +! +! The arguments are: +! \begin{description} +! \item[n] +! index of the nest +! \item[findex] +! forcing dataset index +! \end{description} +! +! The routines invoked are: +! \begin{description} +! \item[LIS\_tick](\ref{LIS_tick}) \newline +! call to advance or retract time +! \item[era5cdsfiles](\ref{era5cdsfiles}) \newline +! Puts together appropriate file name for 1 hour intervals +! \item[read\_era5cds](\ref{read_era5cds}) \newline +! call to read the ERA5 data and perform spatial interpolation +! \end{description} +!EOP + integer :: order + integer :: ferror + character(len=LIS_CONST_PATH_LEN) :: instfilename, lmlfilename + character(len=LIS_CONST_PATH_LEN) :: accfilename, prevaccfilename + integer :: c,kk,f + integer :: yr1, mo1, da1, hr1, mn1, ss1, doy1 + integer :: yr2, mo2, da2, hr2, mn2, ss2, doy2 + real*8 :: time1, time2, timenow + real :: gmt1, gmt2 + real :: ts1, ts2 + + integer :: movetime ! Flag to move bookend2 files to bookend1 + logical :: retrieve_file + + external :: era5cdsfiles + external :: read_era5cds + external :: LISrun + +! _________________________________________________________ + + if( LIS_rc%nts(n).gt.3600 ) then ! > 1-hr timestep + write(LIS_logunit,*) '[ERR] When running LIS with ERA5CDS, the clock ' + write(LIS_logunit,*) '[ERR] should run with a timestep less than or ' + write(LIS_logunit,*) '[ERR] equal to one hour.',LIS_rc%nts(n) + call LIS_endrun() + endif + if(LIS_rc%ts.gt.3600) then + write(LIS_logunit,*) '[ERR] The model timestep is > forcing data timestep' + write(LIS_logunit,*) '[ERR] LIS does not support this mode currently' + write(LIS_logunit,*) '[ERR] Program stopping ...',LIS_rc%ts + call LIS_endrun() + endif + + era5cds_struc(n)%findtime1 = 0 + era5cds_struc(n)%findtime2 = 0 + movetime = 0 + retrieve_file = .false. + +!=== Determine Required ERA5 Data Times (The previous hour and the future hour) + yr1=LIS_rc%yr + mo1=LIS_rc%mo + da1=LIS_rc%da + hr1=LIS_rc%hr + mn1=LIS_rc%mn + ss1=0 + ts1=0 + call LIS_tick(timenow,doy1,gmt1,yr1,mo1,da1,hr1,mn1,ss1,ts1) + + if( timenow < era5cds_struc(n)%validstart ) then + write(LIS_logunit,*) '[ERR] ERA5CDS forecast begins at 7z on 1940/01/01' + write(LIS_logunit,*) '[ERR] Starting hour should be set to 7z ' + write(LIS_logunit,*) '[ERR] in lis.config file. Stopping ... ' + call LIS_endrun() + endif + + yr1 = LIS_rc%yr !previous hour + mo1=LIS_rc%mo + da1=LIS_rc%da + hr1=LIS_rc%hr + mn1=0 + ss1=0 + ts1=0 + call LIS_tick(time1,doy1,gmt1,yr1,mo1,da1,hr1,mn1,ss1,ts1) + + yr2=LIS_rc%yr !next hour + mo2=LIS_rc%mo + da2=LIS_rc%da + hr2=LIS_rc%hr + mn2=0 + ss2=0 + ts2=60*60 + call LIS_tick(time2,doy2,gmt2,yr2,mo2,da2,hr2,mn2,ss2,ts2) + + ! First time step, read in data and set metdata1 + if(LIS_rc%tscount(n).eq.1 .or.LIS_rc%rstflag(n).eq.1 ) then + era5cds_struc(n)%findtime1=1 + era5cds_struc(n)%findtime2=1 + movetime=0 + LIS_rc%rstflag(n) = 0 + era5cds_struc(n)%mon = mo1 ! current month + endif + + ! Flag date/time when it is the switch of a month: + if( era5cds_struc(n)%mon .ne. mo2 ) then + era5cds_struc(n)%findtime1 = 1 + era5cds_struc(n)%mon = mo2 + ! advance filename to next hour + yr1 = yr2 + mo1 = mo2 + da1 = da2 + hr1 = hr2 + endif + + ! Flag date/time for when to open and read next ERA5 file: + if( LIS_get_nstep(LIS_rc,n) == 1 & + .or. (era5cds_struc(n)%findtime1 == 1 ) ) then + retrieve_file = .true. + era5cds_struc(n)%findtime1=0 + endif + + if(timenow.ge.era5cds_struc(n)%era5cdstime2) then + movetime = 1 + era5cds_struc(n)%findtime2 = 1 + endif + + if(movetime.eq.1) then + era5cds_struc(n)%era5cdstime1=era5cds_struc(n)%era5cdstime2 + do f=1,LIS_rc%met_nf(findex) + do c=1,LIS_rc%ngrid(n) + era5cds_struc(n)%metdata1(:,f,c)=era5cds_struc(n)%metdata2(:,f,c) + enddo + enddo + endif !end of movetime=1 + + if( retrieve_file ) then + !- Obtaining ERA5 File: + do kk= era5cds_struc(n)%st_iterid, era5cds_struc(n)%en_iterid + if (LIS_get_nstep(LIS_rc,n) == 1) then + order = 1 + else + order = 2 + era5cds_struc(n)%findtime2 = 0 + endif + call era5cdsfiles(n,kk,findex,era5cds_struc(n)%era5cdsdir, & + yr1, mo1, da1, hr1, instfilename, accfilename, lmlfilename, & + prevaccfilename) + write(LIS_logunit,*) '[INFO] opening Bookend1 ',LIS_get_nstep(LIS_rc,n),hr1 + call read_era5cds(n, kk, order, yr1, mo1, da1, hr1, retrieve_file,& + findex, instfilename, accfilename, lmlfilename, & + prevaccfilename, ferror) + enddo + + if(ferror.ge.1) then !successfully retrieved forcing data + era5cds_struc(n)%era5cdstime1=time1 + retrieve_file = .false. + if(era5cds_struc(n)%findtime2==0) era5cds_struc(n)%era5cdstime2=time2 + else !ferror still=0 + write(LIS_logunit,*)'[ERR] ERA5CDS data missing file 1' + call LISrun() + endif + endif + + if(era5cds_struc(n)%findtime2.eq.1) then + ! just need to assign metdata2 + do kk= era5cds_struc(n)%st_iterid, era5cds_struc(n)%en_iterid + order = 2 + call era5cdsfiles(n,kk,findex,era5cds_struc(n)%era5cdsdir, & + yr2, mo2, da2, hr2, instfilename, accfilename, lmlfilename, & + prevaccfilename) + write(LIS_logunit,*) '[INFO] using Bookend2 ',LIS_get_nstep(LIS_rc,n),hr2 + call read_era5cds(n, kk, order, yr2, mo2, da2, hr2, retrieve_file, & + findex, instfilename, accfilename, lmlfilename, & + prevaccfilename, ferror) + end do + + if(ferror.ge.1) then !successfully retrieved forcing data + era5cds_struc(n)%era5cdstime2=time2 + else + write(LIS_logunit,*)'[ERR] ERA5CDS data missing file 2' + call LIS_endrun() + endif + endif + +end subroutine get_era5cds + + +!BOP +! !ROUTINE: era5cdsfiles +! \label{era5cdsfiles} +! +! !INTERFACE: +subroutine era5cdsfiles(n, kk, findex, era5cdsdir, yr, mo, da, hr, & + instfilename, accfilename, lmlfilename, prevaccfilename) + +! !USES: + use LIS_coreMod + use LIS_logMod + use LIS_forecastMod + use LIS_timeMgrMod + + implicit none +! !ARGUMENTS: + integer :: n + integer :: kk + integer :: findex + character(len=*), intent(in) :: era5cdsdir + integer, intent(in) :: yr,mo,da,hr + character(len=*), intent(out) :: instfilename, accfilename, lmlfilename + character(len=*), intent(out) :: prevaccfilename + +! !DESCRIPTION: +! This subroutine puts together ERA5 file names for +! daily netcdf files +! +! The arguments are: +! \begin{description} +! \item[era5cdsdir] +! Name of the ERA5 directory +! \item[yr] +! year +! \item[mo] +! month +! \item[da] +! day of month +! \item[instfilename] +! name of the timestamped ERA5 instantaneous file +! \item[accfilename] +! name of the timestamped ERA5 accumulation file +! \item[lmlfilename] +! name of the timestamped ERA5 lowest model level file +! \item[prevaccfilename] +! name of the previous month timestamped ERA5 accumulation file +! \end{description} +! +!EOP + + character*4 :: cyear + character*2 :: cmonth + character*2 :: cdate + character*2 :: chour + integer :: mn, ss + integer :: iyr,imo,ida,ihr,imn,iss,ts,idoy + real :: igmt + real*8 :: itime + character :: ciyr*4, cimo*2, cida*2 + + mn = 0 + ss = 0 + + if(LIS_rc%forecastMode.eq.0) then !hindcast run + write(unit=cyear, fmt='(i4.4)') yr + write(unit=cmonth,fmt='(i2.2)') mo + write(unit=cdate,fmt='(i2.2)') da + write(unit=chour,fmt='(i2.2)') hr + + instfilename = trim(era5cdsdir)//'/'//cyear//'/era5_forcing_'//cyear//cmonth//'_instant.nc' + lmlfilename = trim(era5cdsdir)//'/'//cyear//'/era5_forcing_'//cyear//cmonth//'_ml_instant.nc' + accfilename = trim(era5cdsdir)//'/'//cyear//'/era5_forcing_'//cyear//cmonth//'_accum.nc' + ! accum file starts from 7z on the 1st of current month and ends on 6z on + ! the 1st of next month, and need previous months + !== roll back for last day of last month + iyr=yr; imo=mo; ida=da + ihr=hr; imn=0; iss=0 + ts = -24*60*60*da + call LIS_tick(itime,idoy,igmt,iyr,imo,ida,ihr,imn,iss,real(ts)) + write(ciyr, '(i4.4)') iyr + write(cimo, '(i2.2)') imo + write(cida, '(i2.2)') ida + if ( imo /= mo ) then + prevaccfilename = trim(era5cdsdir)//'/'//ciyr//'/era5_forcing_'//ciyr//cimo//'_accum.nc' + else + prevaccfilename = "none" + endif + + else !forecast mode + !sample yr, mo, da + + call LIS_sample_forecastDate(n, kk, findex, yr,mo,da) + + write(unit=cyear, fmt='(i4.4)') yr + write(unit=cmonth,fmt='(i2.2)') mo + write(unit=cdate,fmt='(i2.2)') da + + instfilename = trim(era5cdsdir)//'/'//cyear//'/era5_forcing_'//cyear//cmonth//'_instant.nc' + accfilename = trim(era5cdsdir)//'/'//cyear//'/era5_forcing_'//cyear//cmonth//'_accum.nc' + lmlfilename = trim(era5cdsdir)//'/'//cyear//'/era5_forcing_'//cyear//cmonth//'_ml_instant.nc' + ! accum file starts from 7z on the 1st of current month and ends on 6z on + ! the 1st of next month, and need previous months + !== roll back for last day of last month + iyr=yr; imo=mo; ida=da + ihr=hr; imn=0; iss=0 + ts = -24*60*60*da + call LIS_tick(itime,idoy,igmt,iyr,imo,ida,ihr,imn,iss,real(ts)) + write(ciyr, '(i4.4)') iyr + write(cimo, '(i2.2)') imo + write(cida, '(i2.2)') ida + if ( imo /= mo ) then + prevaccfilename = trim(era5cdsdir)//'/'//ciyr//'/era5_forcing_'//ciyr//cimo//'_accum.nc' + else + prevaccfilename = "none" + endif + + endif +end subroutine era5cdsfiles + diff --git a/lis/metforcing/era5cds/read_era5cds.F90 b/lis/metforcing/era5cds/read_era5cds.F90 new file mode 100644 index 0000000000..46e122b836 --- /dev/null +++ b/lis/metforcing/era5cds/read_era5cds.F90 @@ -0,0 +1,1086 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LIS_misc.h" +!BOP +! +! !ROUTINE: read_era5cds +! \label{read_era5cds} +! +! !REVISION HISTORY: +! 23 dec 2019: Sujay Kumar, initial code +! 15 apr 2025: Hiroko Beaudoing, adopted ERA5 routines for the public CDS +! data format +! 23 apr 2026: Hiroko Beaudoing, added handling of near real-time data +! +! !INTERFACE: +subroutine read_era5cds(n, kk, order, year, month, day, hour, read_flag, & + findex,& + instfile, accfile, lmlfile, prevaccfile, ferror) +! !USES: + use LIS_coreMod, only : LIS_rc + use LIS_logMod + use LIS_FORC_AttributesMod + use LIS_metforcingMod, only : LIS_forc + use era5cds_forcingMod, only : era5cds_struc +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: kk + integer, intent(in) :: order ! lower(1) or upper(2) time interval bdry + integer, intent(in) :: year + integer, intent(in) :: month + integer, intent(in) :: day + integer, intent(in) :: hour + logical, intent(in) :: read_flag + integer, intent(in) :: findex + character(len=*), intent(in) :: instfile, accfile, lmlfile, prevaccfile + integer, intent(out) :: ferror + +! +! !DESCRIPTION: +! Mothly ERA5 data is read in at the beginning of simulation and/or month based +! on read_flag flag. For the given time, assigns appropriate forcing fields, +! transforms into 9 LIS forcing parameters and interpolates to the LIS domain. +! +! ERA5 model output variables used to force LIS are provided in 3 files: +! 1) inst, Instantaneous values, available every hour. +! 2) accum, Time integrated values (accumulations) over the past hour. +! 3) lmlinst, Lowest Model Level Instantaneous values, available every hour. +! +! ERA5 FORCING VARIABLES: +! 1. T2M inst, Atmospheric temperature (K) (2m) +! 2. D2M inst, Dewpoint temperature (K) (2m) +! 3. SSRD accum, surface solar radiation downwards [J m**-2] +! 4. STRD accum, surface thermal radiation downwards [J m**-2] +! 5. U10M inst, zonal wind [m/s] (10m) +! 6. V10M inst, meridional wind [m/s] (10m) +! 7. SP inst, surface pressure [Pa] +! 8. LSP accum, large scale precipitation [m] +! 9. CP accum, convective precipitation [m] +! ERA5 LOWEST MODEL LEVEL VARIABLES: +! 1. T inst, air temperature (K) +! 2. Q inst, specific_humidity (kg kg**-1) +! 5. U inst, zonal wind [m/s] +! 6. V inst, meridional wind [m/s] +! +! also inst file includes: +! Z: Geopotential (m**2 s**-2) +! +! NOTE 1: be aware that ECMWF outputs large-scale and convective precipitation +! separately. For total precipitation, need to sum the two fields, +! LSP+CP=TP. +! NOTE 2: SW & LW flux accumulations from time1 are interpolated via new zterp. +! NOTE 3: accum files starts at 7z on the 1st day of month and ends at 6z +! on the 1st day of next month. +! NOTE 4: Current month includes data up to 5 days ago. If the simulation end +! time is past avaiable data, forcing fields are assigned to the last +! available data with [WARN] message. LSM fields will be invalid. + +! The arguments are: +! \begin{description} +! \item[order] +! flag indicating which data to be read (order=1, initial read, +! order=2, assign the next 1 hourly instance ) +! \item[n] +! index of the nest +! \item[name] +! name of the 1 hour ERA5 analysis file +! \item[tscount] +! time step count +! \item[ferror] +! return error code (1 indicates success) +! \end{description} +! +! The routines invoked are: +! \begin{description} +! \item[era5grid\_2\_lisgrid](\ref{era5grid_2_lisgrid}) \newline +! reverse south-north direction and shift 180 deg east-west +! \item[interp\_era5cds\_var](\ref{interp_era5cds_var}) \newline +! spatially interpolate the forcing data using interpolation of choice +! \item[assign\_processed\_era5cdsf](\ref{assign_processed_era5cdsf}) \newline +! assigns interpolated forcing data to the module data structure +! \end{description} +!EOP + + integer :: ftn + integer :: tmpId, qId, uwindId, vwindId, lwdId, psId, rainfId, crainfId + integer :: swdId, timeId + integer :: c,r,k,l,i,j,ll + integer :: tindex,atindex + integer :: mo,rec_size, prev_rec_size + integer :: start_time_index + logical :: file_exists + real :: missingValue + integer :: dimid, dim_len + integer :: cyear, cmonth + character(len=8) :: date + character(len=10) :: time + character(len=5) :: zone + integer, dimension(8) :: current + + real, allocatable :: tair(:,:,:) + real, allocatable :: qair(:,:,:) + real, allocatable :: swd(:,:,:) + real, allocatable :: lwd(:,:,:) + real, allocatable :: uwind(:,:,:) + real, allocatable :: vwind(:,:,:) + real, allocatable :: ps(:,:,:) + real, allocatable :: rainf(:,:,:) + real, allocatable :: crainf(:,:,:) + real, allocatable :: data4d(:,:,:,:) + real, allocatable :: var1(:,:,:) + + real, allocatable, dimension(:,:) :: datalis + + integer, parameter :: n_last_steps = 7 + integer :: days(12) + data days /31,28,31,30,31,30,31,31,30,31,30,31/ + character(len=*), parameter :: dim_name = "valid_time" + real, parameter :: epsln = 0.622 !Constants' values taken from + real, parameter :: A = 2.53E8, B=5.42E3 !Roger&Yau, A Short Course + !in Cloud Physics, pp.12-17 + ! A in kPa, B in K. + real :: p_kPa + + external :: era5grid_2_lisgrid + external :: interp_era5cds_var + external :: assign_processed_era5cdsf + +! __________________________________________________________________________ + + ferror = 0 ! 1 success + tindex = (day - 1)*24 + hour + 1 + atindex = tindex - n_last_steps + + if(read_flag) then +#if (defined USE_NETCDF4) + + if((mod(year,4) .eq. 0 .and. mod(year, 100).ne.0) &!leap year + .or.(mod(year,400) .eq.0)) then + days(2) = 29 + else + days(2) = 28 + endif + + mo = LIS_rc%lnc(n)*LIS_rc%lnr(n) + rec_size = days(month)*24 + era5cds_struc(n)%tdimsize = rec_size + call date_and_time(date, time, zone, current) + cyear = current(1) + cmonth = current(2) + ! present month has incomplete data record + if (cmonth == month .and. cyear == year) then + call LIS_verify(nf90_open(path=trim(instfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + ! Get the ID for the dimension named "var_time" + call LIS_verify(nf90_inq_dimid(ftn, dim_name, dimid), & + 'nf90_inq_dimid failed in read_era5cds') + ! Get the length (size) of that dimension + call LIS_verify(nf90_inquire_dimension(ftn, dimid, len = dim_len), & + 'nf90_inquire_dimension failed in read_era5cds') + call LIS_verify(nf90_close(ftn), & + 'failed to close file in read_era5cds inst in lml') + rec_size = dim_len + write(LIS_logunit,*)'[INFO] in current month, reset rec_size =', & + rec_size + era5cds_struc(n)%tdimsize = rec_size + endif + + allocate(datalis(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold)) + !=== If using the lowest model level forcing (inst and mlinst): + !=== "t","q","u","v", and "sp " ==== + if(era5cds_struc(n)%uselml.eq.1) then + inquire(file=lmlfile,exist=file_exists) + if(file_exists) then + era5cds_struc(n)%ps = LIS_rc%udef + era5cds_struc(n)%tair = LIS_rc%udef + era5cds_struc(n)%qair = LIS_rc%udef + era5cds_struc(n)%uwind = LIS_rc%udef + era5cds_struc(n)%vwind = LIS_rc%udef + write(LIS_logunit,*)'[INFO] Reading ERA5 file (bookend,', & + order,' -',trim(instfile), ')' + + call LIS_verify(nf90_open(path=trim(instfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + + allocate(ps(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + + call LIS_verify(nf90_inq_varid(ftn,'sp',psId), & + 'nf90_inq_varid failed for psurf in read_era5cds') + call LIS_verify(nf90_get_var(ftn,psId, ps),& + 'nf90_get_var failed for ps in read_era5cds') + call LIS_verify(nf90_get_att(ftn,psId, "_FillValue", & + missingValue), & + 'nf90_get_att failed in read_era5cds for sp FillValue') + call LIS_verify(nf90_close(ftn), & + 'failed to close file in read_era5cds inst in lml') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,ps(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%ps(:,l)) + enddo + deallocate(ps) + + write(LIS_logunit,*)'[INFO] Reading ERA5 file (bookend,', & + order,' -',trim(lmlfile), ')' + call LIS_verify(nf90_open(path=trim(lmlfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + + ! model level data are in 4-dimentional + allocate(data4d(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,1,rec_size)) + call LIS_verify(nf90_inq_varid(ftn,'t',tmpId), & + 'nf90_inq_varid failed for t in read_era5cds') + call LIS_verify(nf90_get_var(ftn,tmpId, data4d),& + 'nf90_get_var failed for t in read_era5cds') + call LIS_verify(nf90_get_att(ftn,tmpId, "_FillValue", & + missingValue), & + 'nf90_get_att failed in read_era5cds for t FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,data4d(:,:,1,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%tair(:,l)) + enddo + + call LIS_verify(nf90_inq_varid(ftn,'q',qId), & + 'nf90_inq_varid failed for q in read_era5cds') + call LIS_verify(nf90_get_var(ftn,qId, data4d),& + 'nf90_get_var failed for q in read_era5cds') + call LIS_verify(nf90_get_att(ftn,qId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for q FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,data4d(:,:,1,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%qair(:,l)) + enddo + + call LIS_verify(nf90_inq_varid(ftn,'u',uwindId), & + 'nf90_inq_varid failed for u in read_era5cds') + call LIS_verify(nf90_get_var(ftn,uwindId, data4d),& + 'nf90_get_var failed for u in read_era5cds') + call LIS_verify(nf90_get_att(ftn,uwindId, "_FillValue", & + missingValue), & + 'nf90_get_att failed in read_era5cds for u FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,data4d(:,:,1,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%uwind(:,l)) + enddo + + call LIS_verify(nf90_inq_varid(ftn,'v',vwindId), & + 'nf90_inq_varid failed for v in read_era5cds') + call LIS_verify(nf90_get_var(ftn,vwindId, data4d),& + 'nf90_get_var failed for v in read_era5cds') + call LIS_verify(nf90_get_att(ftn,vwindId, "_FillValue", & + missingValue), & + 'nf90_get_att failed in read_era5cds for v FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,data4d(:,:,1,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%vwind(:,l)) + enddo + deallocate(data4d) + + call LIS_verify(nf90_close(ftn), & + 'failed to close file in read_era5cds lml') + ferror = 1 + else + write(LIS_logunit,*) '[ERR] ',trim(lmlfile)//' does not exist' + call LIS_endrun() + endif + + else !=== If using 2m and 10m forcing (inst): + !=== "t2m","d2m","u10","v10","sp " ==== + + inquire(file=instfile,exist=file_exists) + if(file_exists) then + era5cds_struc(n)%ps = LIS_rc%udef + era5cds_struc(n)%tair = LIS_rc%udef + era5cds_struc(n)%qair = LIS_rc%udef + era5cds_struc(n)%uwind = LIS_rc%udef + era5cds_struc(n)%vwind = LIS_rc%udef + + write(LIS_logunit,*)'[INFO] Reading ERA5 2m10m (bookend,', & + order,' -',trim(instfile), ')' + + call LIS_verify(nf90_open(path=trim(instfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + + allocate(ps(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + + call LIS_verify(nf90_inq_varid(ftn,'sp',psId), & + 'nf90_inq_varid failed for psurf in read_era5cds') + call LIS_verify(nf90_get_var(ftn,psId, ps),& + 'nf90_get_var failed for ps in read_era5cds') + call LIS_verify(nf90_get_att(ftn,psId, "_FillValue", & + missingValue), & + 'nf90_get_att failed in read_era5cds for sp FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,ps(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%ps(:,l)) + enddo + deallocate(ps) + + allocate(qair(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold, & + rec_size)) + call LIS_verify(nf90_inq_varid(ftn,'d2m',qId), & + 'nf90_inq_varid failed for d2m in read_era5cds') + call LIS_verify(nf90_get_var(ftn,qId, qair),& + 'nf90_get_var failed for d2m in read_era5cds') + call LIS_verify(nf90_get_att(ftn,qId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for d2m FillValue') + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,qair(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%qair(:,l)) + enddo + deallocate(qair) + !----------------------------------------------------------------- + ! Calculate specific humidity from Dew point Temp. & Sfc Pressure. + ! Approximate: q~epsilon*e/p, p~p_sfc, e=e_s(Td), e_s=A*exp(-B/T) + !----------------------------------------------------------------- + do l=1,rec_size + do r=1,LIS_rc%lnr(n) + do c=1,LIS_rc%lnc(n) + k=c+(r-1)*LIS_rc%lnc(n) + if(era5cds_struc(n)%ps(k,l).ne.LIS_rc%udef.and.& + era5cds_struc(n)%qair(k,l).ne.LIS_rc%udef) then + p_kPa = era5cds_struc(n)%ps(k,l) / 1000.0 + era5cds_struc(n)%qair(k,l) = & + epsln*(A * EXP(-B/era5cds_struc(n)%qair(k,l))) / p_kPa + endif + enddo + enddo + enddo + + allocate(tair(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LIS_verify(nf90_inq_varid(ftn,'t2m',tmpId), & + 'nf90_inq_varid failed for t2m in read_era5cds') + call LIS_verify(nf90_get_var(ftn,tmpId, tair),& + 'nf90_get_var failed for t2m in read_era5cds') + call LIS_verify(nf90_get_att(ftn,tmpId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for t2m FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,tair(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%tair(:,l)) + enddo + deallocate(tair) + + allocate(uwind(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold, & + rec_size)) + call LIS_verify(nf90_inq_varid(ftn,'u10',uwindId), & + 'nf90_inq_varid failed for u10 in read_era5cds') + call LIS_verify(nf90_get_var(ftn,uwindId, uwind),& + 'nf90_get_var failed for u10 in read_era5cds') + call LIS_verify(nf90_get_att(ftn,uwindId, "_FillValue", & + missingValue), & + 'nf90_get_att failed in read_era5cds for u10 FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,uwind(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%uwind(:,l)) + enddo + deallocate(uwind) + + allocate(vwind(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold, & + rec_size)) + call LIS_verify(nf90_inq_varid(ftn,'v10',vwindId), & + 'nf90_inq_varid failed for v10 in read_era5cds') + call LIS_verify(nf90_get_var(ftn,vwindId, vwind),& + 'nf90_get_var failed for v10 in read_era5cds') + call LIS_verify(nf90_get_att(ftn,vwindId, "_FillValue", & + missingValue), & + 'nf90_get_att failed in read_era5cds for v10 FillValue') + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,vwind(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%vwind(:,l)) + enddo + deallocate(vwind) + + call LIS_verify(nf90_close(ftn), & + 'failed to close file in read_era5cds inst') + ferror = 1 + else + write(LIS_logunit,*) '[ERR] ',trim(instfile)//' does not exist' + call LIS_endrun() + endif + + endif !=== era5cds_struc(n)%uselml + + !=== Read precipitation and radiation fields in accum file + !=== "ssrd","strd","lsp","cp" ==== + inquire(file=accfile,exist=file_exists) + if(file_exists) then + era5cds_struc(n)%rainf = LIS_rc%udef + era5cds_struc(n)%crainf = LIS_rc%udef + era5cds_struc(n)%swd = LIS_rc%udef + era5cds_struc(n)%lwd = LIS_rc%udef + + write(LIS_logunit,*)'[INFO] Reading ERA5 file (bookend,', & + order,' -',trim(accfile), ')' + call LIS_verify(nf90_open(path=trim(accfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + + allocate(rainf(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LIS_verify(nf90_inq_varid(ftn,'lsp',rainfId), & + 'nf90_inq_varid failed for lsp in read_era5cds') + call LIS_verify(nf90_get_var(ftn,rainfId, rainf),& + 'nf90_get_var failed for lsp in read_era5cds') + call LIS_verify(nf90_get_att(ftn,rainfId, "_FillValue", & + missingValue), & + 'nf90_get_att failed in read_era5cds for lsp FillValue') + + ! large scale rainfall: lsp [m] -> mm/s + do l=1,rec_size + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(rainf(i,j,l).ne.missingValue) then + rainf(i,j,l) = rainf(i,j,l)*1000.0/(1.0*60*60) + endif + enddo + enddo + enddo + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,rainf(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .true.,& + era5cds_struc(n)%rainf(:,l)) + enddo + deallocate(rainf) + + allocate(crainf(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold, & + rec_size)) + call LIS_verify(nf90_inq_varid(ftn,'cp',crainfId), & + 'nf90_inq_varid failed for cp in read_era5cds') + call LIS_verify(nf90_get_var(ftn,crainfId, crainf),& + 'nf90_get_var failed for cp in read_era5cds') + call LIS_verify(nf90_get_att(ftn,crainfId, "_FillValue", & + missingValue), & + 'nf90_get_att failed in read_era5cds for cp FillValue') + + ! convective rainfall: cp [m] -> mm/s + do l=1,rec_size + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(crainf(i,j,l).ne.missingValue) then + crainf(i,j,l) = crainf(i,j,l)*1000.0/(1.0*60*60) + endif + enddo + enddo + enddo + + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,crainf(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .true.,& + era5cds_struc(n)%crainf(:,l)) + enddo + deallocate(crainf) + + !----------------------------------------------------------------- + ! rain = lsp + cp [mm/s] + !----------------------------------------------------------------- + do l=1,rec_size + do r=1,LIS_rc%lnr(n) + do c=1,LIS_rc%lnc(n) + k=c+(r-1)*LIS_rc%lnc(n) + if(era5cds_struc(n)%rainf(k,l).ne.LIS_rc%udef.and.& + era5cds_struc(n)%crainf(k,l).ne.LIS_rc%udef) then + era5cds_struc(n)%rainf(k,l) = & + era5cds_struc(n)%rainf(k,l) + & + era5cds_struc(n)%crainf(k,l) + endif + enddo + enddo + enddo + + allocate(swd(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LIS_verify(nf90_inq_varid(ftn,'ssrd',swdId), & + 'nf90_inq_varid failed for ssrd in read_era5cds') + call LIS_verify(nf90_get_var(ftn,swdId, swd),& + 'nf90_get_var failed for swd in read_era5cds') + call LIS_verify(nf90_get_att(ftn,swdId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for swd FillValue') + + ! swd - convert accumulated field to rate + do l=1,rec_size + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(swd(i,j,l).ne.missingValue) then + swd(i,j,l) = swd(i,j,l)/(1.0*60*60) + endif + enddo + enddo + enddo + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,swd(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%swd(:,l)) + enddo + deallocate(swd) + + allocate(lwd(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold,rec_size)) + call LIS_verify(nf90_inq_varid(ftn,'strd',lwdId), & + 'nf90_inq_varid failed for strd in read_era5cds') + call LIS_verify(nf90_get_var(ftn,lwdId, lwd),& + 'nf90_get_var failed for lwd in read_era5cds') + call LIS_verify(nf90_get_att(ftn,lwdId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for lwd FillValue') + + ! lwd - convert accumulated field to rate + do l=1,rec_size + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(lwd(i,j,l).ne.missingValue) then + lwd(i,j,l) = lwd(i,j,l)/(1.0*60*60) + endif + enddo + enddo + enddo + do l=1,rec_size + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,lwd(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%lwd(:,l)) + enddo + deallocate(lwd) + + call LIS_verify(nf90_close(ftn), & + 'failed to close avg file in read_era5cds') + ferror = 1 + else + write(LIS_logunit,*) '[ERR] ',trim(accfile)//' does not exist' + call LIS_endrun() + + endif + !=== Read previous precipitation and radiation fields in accum file + !=== populate prev* arrays + inquire(file=prevaccfile,exist=file_exists) + if(file_exists) then + ! backward fill: read in last 7 time steps of prevaccfile + write(LIS_logunit,*)'[INFO] Reading prev ERA5 file (bookend,', & + order,' -',trim(prevaccfile), ')' + call LIS_verify(nf90_open(path=trim(prevaccfile), mode=NF90_NOWRITE, & + ncid=ftn), 'nf90_open failed in read_era5cds') + + call LIS_verify(nf90_inq_dimid(ftn,'valid_time',timeId), & + 'nf90_inq_dimid failed for timeId in read_era5cds') + call LIS_verify(nf90_inquire_dimension(ftn,timeId,len=prev_rec_size), & + 'nf90_inquire_dimension failed for timeId in read_era5cds') + start_time_index = prev_rec_size - n_last_steps + 1 + + allocate(var1(era5cds_struc(n)%ncold,era5cds_struc(n)%nrold, & + n_last_steps)) + call LIS_verify(nf90_inq_varid(ftn,'lsp',rainfId), & + 'nf90_inq_varid failed for lsp in read_era5cds') + call LIS_verify(nf90_get_var(ftn,rainfId, var1, & + start=(/1,1,start_time_index/), & + count=(/era5cds_struc(n)%ncold,era5cds_struc(n)%nrold, & + n_last_steps/)),& + 'nf90_get_var failed for lsp in read_era5cds') + call LIS_verify(nf90_get_att(ftn,rainfId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for lsp FillValue') + + ! lsp [m] -> mm/s + do l=1,n_last_steps + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(var1(i,j,l).ne.missingValue) then + var1(i,j,l) = var1(i,j,l)*1000.0/(1.0*60*60) + endif + enddo + enddo + enddo + + do l=1,n_last_steps + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,var1(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.true.,& + era5cds_struc(n)%prev_rainf(:,l)) + enddo + + call LIS_verify(nf90_inq_varid(ftn,'cp',crainfId), & + 'nf90_inq_varid failed for cp in read_era5cds') + call LIS_verify(nf90_get_var(ftn,crainfId, var1, & + start=(/1,1,start_time_index/), & + count=(/era5cds_struc(n)%ncold,era5cds_struc(n)%nrold, & + n_last_steps/)),& + 'nf90_get_var failed for cp in read_era5cds') + + ! cp [m] -> mm/s + do l=1,n_last_steps + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(var1(i,j,l).ne.missingValue) then + var1(i,j,l) = var1(i,j,l)*1000.0/(1.0*60*60) + endif + enddo + enddo + enddo + + do l=1,n_last_steps + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,var1(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue,.true.,& + era5cds_struc(n)%prev_crainf(:,l)) + enddo + + !----------------------------------------------------------------- + ! rain = lsp + cp [mm/s] + !----------------------------------------------------------------- + do l=1,n_last_steps + do r=1,LIS_rc%lnr(n) + do c=1,LIS_rc%lnc(n) + k=c+(r-1)*LIS_rc%lnc(n) + if(era5cds_struc(n)%prev_rainf(k,l).ne.LIS_rc%udef.and.& + era5cds_struc(n)%prev_crainf(k,l).ne.LIS_rc%udef) then + era5cds_struc(n)%prev_rainf(k,l) = & + era5cds_struc(n)%prev_rainf(k,l) + & + era5cds_struc(n)%prev_crainf(k,l) + endif + enddo + enddo + enddo + + call LIS_verify(nf90_inq_varid(ftn,'ssrd',swdId), & + 'nf90_inq_varid failed for ssrd in read_era5cds') + call LIS_verify(nf90_get_var(ftn,swdId, var1, & + start=(/1,1,start_time_index/), & + count=(/era5cds_struc(n)%ncold,era5cds_struc(n)%nrold, & + n_last_steps/)),& + 'nf90_get_var failed for swd in read_era5cds') + call LIS_verify(nf90_get_att(ftn,swdId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for swd FillValue') + + ! swd - convert accumulated field to rate + do l=1,n_last_steps + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(var1(i,j,l).ne.missingValue) then + var1(i,j,l) = var1(i,j,l)/(1.0*60*60) + endif + enddo + enddo + enddo + do l=1,n_last_steps + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,var1(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%prev_swd(:,l)) + enddo + + call LIS_verify(nf90_inq_varid(ftn,'strd',lwdId), & + 'nf90_inq_varid failed for strd in read_era5cds') + call LIS_verify(nf90_get_var(ftn,lwdId, var1, & + start=(/1,1,start_time_index/), & + count=(/era5cds_struc(n)%ncold,era5cds_struc(n)%nrold, & + n_last_steps/)),& + 'nf90_get_var failed for lwd in read_era5cds') + call LIS_verify(nf90_get_att(ftn,lwdId, "_FillValue", missingValue), & + 'nf90_get_att failed in read_era5cds for lwd FillValue') + + ! lwd - convert accumulated field to rate + do l=1,n_last_steps + do j=1,era5cds_struc(n)%nrold + do i=1,era5cds_struc(n)%ncold + if(var1(i,j,l).ne.missingValue) then + var1(i,j,l) = var1(i,j,l)/(1.0*60*60) + endif + enddo + enddo + enddo + do l=1,n_last_steps + call era5grid_2_lisgrid(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold,var1(:,:,l),datalis) + call interp_era5cds_var(n,findex,month,datalis,missingValue, & + .false.,& + era5cds_struc(n)%prev_lwd(:,l)) + enddo + deallocate(var1) + + call LIS_verify(nf90_close(ftn), & + 'failed to close avg file in read_era5cds') + else ! prevaccfile exists + ! at start of ERA5 in Jan 1940, no data for 0-6z + write(LIS_logunit,*) '[INFO] ',trim(prevaccfile)//' not available!' + do l = 1, n_last_steps + ll = rec_size - n_last_steps + l + era5cds_struc(n)%prev_rainf(:,l) = LIS_rc%udef + era5cds_struc(n)%prev_crainf(:,l) = LIS_rc%udef + era5cds_struc(n)%prev_swd(:,l) = LIS_rc%udef + era5cds_struc(n)%prev_lwd(:,l) = LIS_rc%udef + enddo + endif + + deallocate(datalis) +#endif + endif !if(read_flag) + + ! make sure current time is within valid data + if (era5cds_struc(n)%tdimsize < tindex) then + write(LIS_logunit,*) & + '[WARN] Current month passed available data, filling with last data' + tindex = era5cds_struc(n)%tdimsize + atindex = era5cds_struc(n)%tdimsize - n_last_steps + endif !(rec_size < tindex) + + call assign_processed_era5cdsf(n,kk,order,1,era5cds_struc(n)%tair(:,tindex)) + call assign_processed_era5cdsf(n,kk,order,2,era5cds_struc(n)%qair(:,tindex)) + call assign_processed_era5cdsf(n,kk,order,5,era5cds_struc(n)%uwind(:,tindex)) + call assign_processed_era5cdsf(n,kk,order,6,era5cds_struc(n)%vwind(:,tindex)) + call assign_processed_era5cdsf(n,kk,order,7,era5cds_struc(n)%ps(:,tindex)) + if ( atindex .le. 0 ) then + atindex = atindex + n_last_steps + call assign_processed_era5cdsf(n,kk,order,3, & + era5cds_struc(n)%prev_swd(:,atindex)) + call assign_processed_era5cdsf(n,kk,order,4, & + era5cds_struc(n)%prev_lwd(:,atindex)) + call assign_processed_era5cdsf(n,kk,order,8, & + era5cds_struc(n)%prev_rainf(:,atindex)) + call assign_processed_era5cdsf(n,kk,order,9, & + era5cds_struc(n)%prev_crainf(:,atindex)) + else + call assign_processed_era5cdsf(n,kk,order,3, & + era5cds_struc(n)%swd(:,atindex)) + call assign_processed_era5cdsf(n,kk,order,4, & + era5cds_struc(n)%lwd(:,atindex)) + call assign_processed_era5cdsf(n,kk,order,8, & + era5cds_struc(n)%rainf(:,atindex)) + call assign_processed_era5cdsf(n,kk,order,9, & + era5cds_struc(n)%crainf(:,atindex)) + endif + if ( .not. read_flag .and. ferror == 0 ) ferror = 1 + +end subroutine read_era5cds + +!BOP +! +! !ROUTINE: interp_era5cds_var +! \label{interp_era5cds_var} +! +! !INTERFACE: +subroutine interp_era5cds_var(n,findex, month, input_var, missingValue, & + pcp_flag, output_var) + +! !USES: + use LIS_coreMod + use LIS_logMod + use LIS_spatialDownscalingMod + use era5cds_forcingMod, only : era5cds_struc +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + implicit none + +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: findex + integer, intent(in) :: month + real, intent(in) :: input_var(era5cds_struc(n)%ncold, & + era5cds_struc(n)%nrold) + real, intent(in) :: missingValue + logical, intent(in) :: pcp_flag + real, intent(out) :: output_var(LIS_rc%lnc(n)*LIS_rc%lnr(n)) + +! +! !DESCRIPTION: +! This subroutine spatially interpolates a ERA5CDS field +! to the LIS running domain +! +!EOP + + integer :: c,r,k,iret + real :: f(era5cds_struc(n)%ncold*era5cds_struc(n)%nrold) + logical*1 :: lb(era5cds_struc(n)%mi) + logical*1 :: lo(LIS_rc%lnc(n)*LIS_rc%lnr(n)) + integer :: input_size + integer :: input_nc, input_nr + integer :: count1,nrec + + external :: conserv_interp + external :: bilinear_interp + external :: neighbor_interp + external :: upscaleByAveraging + +! _____________________________________________________________ + + input_size = era5cds_struc(n)%mi + output_var = LIS_rc%udef + lo = .true. + input_nc = era5cds_struc(n)%ncold + input_nr = era5cds_struc(n)%nrold + +!----------------------------------------------------------------------- +! Apply corrections +!----------------------------------------------------------------------- + + lb = .false. + do r=1,era5cds_struc(n)%nrold + do c=1,era5cds_struc(n)%ncold + k= c+(r-1)*era5cds_struc(n)%ncold + f(k) = input_var(c,r) + if(f(k) .ne. missingValue) then + lb(k) = .true. + endif + enddo + enddo +!----------------------------------------------------------------------- +! Apply downscaling +!----------------------------------------------------------------------- + + if(pcp_flag.and.LIS_rc%pcp_downscale(findex).ne.0) then + call LIS_generatePcpClimoRatioField(n,findex,"ERA5CDS",& + month, & + input_size, & + f, & + lb) + endif + + if(pcp_flag.and.& + trim(era5cds_struc(n)%met_interp).eq."budget-bilinear") then + + call conserv_interp(LIS_rc%gridDesc(n,:),lb,f,lo,& + output_var, & + era5cds_struc(n)%mi,LIS_rc%lnc(n)*LIS_rc%lnr(n),& + LIS_domain(n)%lat, LIS_domain(n)%lon,& + era5cds_struc(n)%w112,era5cds_struc(n)%w122,& + era5cds_struc(n)%w212,era5cds_struc(n)%w222,& + era5cds_struc(n)%n112,era5cds_struc(n)%n122,& + era5cds_struc(n)%n212,era5cds_struc(n)%n222,& + LIS_rc%udef, iret) + + elseif(trim(era5cds_struc(n)%met_interp).eq."bilinear".or.& + trim(era5cds_struc(n)%met_interp).eq."budget-bilinear") then + + call bilinear_interp(LIS_rc%gridDesc(n,:),lb,f,lo,& + output_var, & + era5cds_struc(n)%mi,LIS_rc%lnc(n)*LIS_rc%lnr(n), & + LIS_domain(n)%lat, LIS_domain(n)%lon,& + era5cds_struc(n)%w111,era5cds_struc(n)%w121,& + era5cds_struc(n)%w211,era5cds_struc(n)%w221,& + era5cds_struc(n)%n111,era5cds_struc(n)%n121,& + era5cds_struc(n)%n211,era5cds_struc(n)%n221,& + LIS_rc%udef, iret) + + elseif(trim(era5cds_struc(n)%met_interp).eq."neighbor") then + call neighbor_interp(LIS_rc%gridDesc(n,:),lb,f,lo,& + output_var,era5cds_struc(n)%mi,& + LIS_rc%lnc(n)*LIS_rc%lnr(n),& + LIS_domain(n)%lat, LIS_domain(n)%lon,& + era5cds_struc(n)%n113,LIS_rc%udef,iret) + + elseif(trim(era5cds_struc(n)%met_interp).eq."average") then + call upscaleByAveraging(era5cds_struc(n)%mi, & + LIS_rc%lnc(n)*LIS_rc%lnr(n), LIS_rc%udef, & + era5cds_struc(n)%n111, lb, f, lo, output_var) + + elseif(trim(era5cds_struc(n)%met_interp).eq."none") then + ! at 0.25 degree, no interpolation needed. input_var(input_nc*input_rc) + ! If input grid matches the output grid size: + if( era5cds_struc(n)%subset_nc == input_nc .and. & + era5cds_struc(n)%subset_nr == input_nr ) then + output_var = f + ! Otherwise, subset for same domain resolutions: + else + count1 = 0 + do r = 1, era5cds_struc(n)%subset_nr + do c = 1, era5cds_struc(n)%subset_nc + count1 = count1 + 1 + nrec = (era5cds_struc(n)%lat_line(c,r)-1)*input_nc + & + era5cds_struc(n)%lon_line(c,r) + output_var(count1) = f(nrec) + enddo ! columns + enddo ! rows + endif + + else + write(LIS_logunit,*) '[ERR] Spatial interpolation option '//& + trim(LIS_rc%met_interp(findex))//& + ' not supported for ERA5CDS' + call LIS_endrun() + endif + + if( pcp_flag.and.LIS_rc%pcp_downscale(findex).ne.0 ) then + + call LIS_pcpClimoDownscaling(n, findex, month,& + LIS_rc%lnc(n)*LIS_rc%lnr(n), output_var(:), lo) + + endif + +end subroutine interp_era5cds_var +!BOP +! +! !ROUTINE: era5grid_2_lisgrid +! +! !DESCRIPTION: +! Changes grid_data from ECMWF data convention to GLDAS convention +! +! ERA5: North-to-South around Greenwich Meridian +! Global grid. Data are written as flat binary from "upper left to +! lower right" starting at 0.25-degree grid point center coordinates: +! 0.125E,89.875N and going to 0.125W,89.875S. +! +! LIS: South-to-North around Date Line +! Full global grid. Starts at the southernmost latitude and date line, +! going east and then north. +! +! !REVISION HISTORY: +! 10 Apr 2002: Urszula Jambor; Code adapted from +! ecmwfgrid_2_grid2catgrid, by R. Reichle +! 21 Apr 2025: Hiroko Beaudoing; modified berggrid_2_gldasgrid.F90 +! to work for odd ny. +! +! !INTERFACE: +subroutine era5grid_2_lisgrid( nx, ny, grid_data, out_data ) + !EOP + use LIS_logMod, only: LIS_endrun, LIS_logunit + implicit none + + integer, intent(in) :: nx, ny + real, intent(in), dimension(nx,ny) :: grid_data + real, intent(out), dimension(nx*ny):: out_data + real, dimension(nx,ny):: lis_data + + integer :: i, j, m, n, c + real :: tmp, tmp_data1(nx) + + ! ------------------------------------------------------------------ + ! some checks + + if ((nx /= 1440) .or. (ny /= 720)) then + write (LIS_logunit,*) 'era5grid_2_gldasgrid(): This routine has only been' + write (LIS_logunit,*) 'checked for nx=1440 and ny=720. STOPPING' + call LIS_endrun() + end if + if ((mod(nx,2) /= 0)) then + write (LIS_logunit,*) 'era5grid_2_gldasgrid(): This routine can only work' + write (LIS_logunit,*) 'for even nx. STOPPING.' + call LIS_endrun() + end if + + !------------------------------------------------------------------- + + do j=1,ny + + ! swap latitude bands (North-to-South becomes South-to-North) + n = ny-j+1 + tmp_data1 = grid_data(:,j) + + do i=1,nx/2 + + ! shift longitudes (wrapping around Greenwhich Meridian becomes + ! wrapping around Date Line) + m = i + nx/2 + tmp = tmp_data1(i) + tmp_data1(i) = tmp_data1(m) + tmp_data1(m) = tmp + + end do + + lis_data(:,n) = tmp_data1 + + end do + + ! return output in 1D + c=0 + do j=1,ny + do i=1,nx + c=c+1 + out_data(c) = lis_data(i,j) + end do + end do + +end subroutine era5grid_2_lisgrid + +!BOP +! +! !ROUTINE: assign_processed_era5cdsf +! \label{assign_processed_era5cdsf} +! +! !INTERFACE: +subroutine assign_processed_era5cdsf(n,kk,order,var_index,era5forc) +! !USES: + use LIS_coreMod + use era5cds_forcingMod, only : era5cds_struc +! +! !DESCRIPTION: +! This routine assigns the interpolated ERA5 forcing data +! to the module data structures to be used later for +! time interpolation +! +!EOP + implicit none + + integer :: n + integer :: kk + integer :: order + integer :: var_index + real :: era5forc(LIS_rc%lnc(n)*LIS_rc%lnr(n)) + + integer :: c,r + + do r=1,LIS_rc%lnr(n) + do c=1,LIS_rc%lnc(n) + if(LIS_domain(n)%gindex(c,r).ne.-1) then + if(order.eq.1) then + era5cds_struc(n)%metdata1(kk,var_index,& + LIS_domain(n)%gindex(c,r)) = & + era5forc(c+(r-1)*LIS_rc%lnc(n)) + elseif(order.eq.2) then + era5cds_struc(n)%metdata2(kk,var_index,& + LIS_domain(n)%gindex(c,r)) = & + era5forc(c+(r-1)*LIS_rc%lnc(n)) + endif + endif + enddo + enddo +end subroutine assign_processed_era5cdsf diff --git a/lis/metforcing/era5cds/read_era5cds_elev.F90 b/lis/metforcing/era5cds/read_era5cds_elev.F90 new file mode 100644 index 0000000000..0ae6200287 --- /dev/null +++ b/lis/metforcing/era5cds/read_era5cds_elev.F90 @@ -0,0 +1,97 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LIS_misc.h" +!BOP +! !ROUTINE: read_era5cds_elev +! \label{read_era5cds_elev} +! +! !REVISION HISTORY: +! +! 23 dec 2019; Sujay Kumar; Initial Specificaton +! 04 Mar 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: +subroutine read_era5cds_elev(n,findex) +! !USES: + use LIS_coreMod + use LIS_metforcingMod + use LIS_logMod + use era5cds_forcingMod + use LIS_fileIOMod +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: findex +! !DESCRIPTION: +! +! Opens, reads, and interpolates ERA5CDS model elevation to the LIS +! grid. The data will be used to perform any topographical +! adjustments to the forcing. +! +! The arguments are: +! \begin{description} +! \item[n] +! index of the nest +! \end{description} +! +! The routines invoked are: +! \begin{description} +! \item[ij\_to\_latlon](\ref{ij_to_latlon}) \newline +! computes the lat lon values in LIS grid projection +! \end{description} +!EOP + logical :: file_exists + integer :: nid,elevId + integer :: c,r + real :: elev(LIS_rc%gnc(n),LIS_rc%gnr(n)) + real :: elev_subset(LIS_rc%lnc(n),LIS_rc%lnr(n)) + +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + + inquire(file=LIS_rc%paramfile(n), exist=file_exists) + if(file_exists) then + + write(LIS_logunit,*) " Reading ERA5CDS elevation data ... " + + call LIS_verify(nf90_open(path=LIS_rc%paramfile(n),& + mode=NF90_NOWRITE,ncid=nid),& + 'nf90_open failed in read_era5cds_elev') + call LIS_verify(nf90_inq_varid(nid,"ELEV_ERA5CDS",elevId),& + 'nf90_inq_varid failed in read_era5cds_elev') + call LIS_verify(nf90_get_var(nid,elevId,elev),& + 'nf90_get_var failed in read_era5cds_elev') + call LIS_verify(nf90_close(nid)) + + elev_subset(:,:) = elev(LIS_ews_halo_ind(n,LIS_localPet+1): & + LIS_ewe_halo_ind(n,LIS_localPet+1), & + LIS_nss_halo_ind(n,LIS_localPet+1): & + LIS_nse_halo_ind(n,LIS_localPet+1)) + + do r=1,LIS_rc%lnr(n) + do c=1,LIS_rc%lnc(n) + if(LIS_domain(n)%gindex(c,r).ne.-1) then + LIS_forc(n,findex)%modelelev(LIS_domain(n)%gindex(c,r)) =& + elev_subset(c,r) + endif + enddo + enddo + endif + + write(LIS_logunit,*) & + "[INFO] Finish reading ERA5CDS elevation data" +#endif + +end subroutine read_era5cds_elev + diff --git a/lis/metforcing/era5cds/readcrd_era5cds.F90 b/lis/metforcing/era5cds/readcrd_era5cds.F90 new file mode 100644 index 0000000000..7a92b2411e --- /dev/null +++ b/lis/metforcing/era5cds/readcrd_era5cds.F90 @@ -0,0 +1,73 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! +! !ROUTINE: readcrd_era5cds +! \label{readcrd_era5cds} +! +! !REVISION HISTORY: +! 23 Dec 2019: Sujay Kumar, initial code +! 04 Mar 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: +subroutine readcrd_era5cds() +! !USES: + use ESMF + use LIS_coreMod, only : LIS_rc, LIS_config + use LIS_logMod + use era5cds_forcingMod, only : era5cds_struc +! +! !DESCRIPTION: +! +! This routine reads the options specific to ERA5 forcing +! from the LIS configuration file. +! +!EOP + implicit none + + integer :: n,rc + + call ESMF_ConfigFindLabel(LIS_config,"ERA5CDS forcing directory:",rc=rc) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,era5cds_struc(n)%era5cdsdir,& + rc=rc) + call LIS_verify(rc,& + 'ERA5CDS forcing directory: not defined') + enddo + + call ESMF_ConfigFindLabel(LIS_config,"ERA5CDS surface altitude file:",rc=rc) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,era5cds_struc(n)%era5cdsalt_file,& + rc=rc) + call LIS_verify(rc,& + 'ERA5CDS surface altitude file: not defined') + enddo + + call ESMF_ConfigFindLabel(LIS_config, & + "ERA5CDS use lowest model level forcing:",rc=rc) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,era5cds_struc(n)%uselml,rc=rc) + call LIS_verify(rc,& + 'ERA5CDS use lowest model level forcing: not defined') + enddo + + do n=1,LIS_rc%nnest + write(LIS_logunit,*) '[INFO] Using ERA5CDS forcing' + write(LIS_logunit,*) '[INFO] ERA5CDS forcing directory: ',& + trim(era5cds_struc(n)%era5cdsdir) + write(LIS_logunit,*) '[INFO] ERA5CDS use lowest model level forcing: ',& + era5cds_struc(n)%uselml + + era5cds_struc(n)%era5cdstime1 = 3000.0 + era5cds_struc(n)%era5cdstime2 = 0.0 + + enddo +end subroutine readcrd_era5cds diff --git a/lis/metforcing/era5cds/reset_era5cds.F90 b/lis/metforcing/era5cds/reset_era5cds.F90 new file mode 100644 index 0000000000..ab53ba7cc6 --- /dev/null +++ b/lis/metforcing/era5cds/reset_era5cds.F90 @@ -0,0 +1,41 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! !MODULE: reset_era5cds +! \label{reset_era5cds} +! +! !REVISION HISTORY: +! 23 dec 2019: Sujay Kumar, initial code +! 16 apr 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: +subroutine reset_era5cds +! !USES: + use LIS_coreMod, only : LIS_rc + use LIS_timeMgrMod, only : LIS_date2time + use era5cds_forcingMod +! +! !DESCRIPTION: +! Routine to cleanup allocated structures for era5cds forcing. +! +!EOP + implicit none + integer :: n + + do n=1,LIS_rc%nnest + era5cds_struc(n)%startFlag = .true. + era5cds_struc(n)%dayFlag = .true. + era5cds_struc(n)%era5cdstime1 = 3000.0 + era5cds_struc(n)%era5cdstime2 = 0.0 + era5cds_struc(n)%ringtime = 0.0 + era5cds_struc(n)%reset_flag = .true. + enddo +end subroutine reset_era5cds diff --git a/lis/metforcing/era5cds/timeinterp_era5cds.F90 b/lis/metforcing/era5cds/timeinterp_era5cds.F90 new file mode 100644 index 0000000000..3a3e3345e4 --- /dev/null +++ b/lis/metforcing/era5cds/timeinterp_era5cds.F90 @@ -0,0 +1,265 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2024 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! !ROUTINE: timeinterp_era5cds +! \label{timeinterp_era5cds} +! +! !REVISION HISTORY: +! 23 dec 2019: Sujay Kumar, initial code +! 16 Apr 2025: Hiroko Beudoing, adopted ERA5 routines for the public CDS +! data format +! +! !INTERFACE: + +subroutine timeinterp_era5cds(n,findex) + +! !USES: + use ESMF + use LIS_coreMod + use LIS_metforcingMod + use LIS_FORC_AttributesMod + use LIS_constantsMod + use LIS_timeMgrMod + use LIS_logMod + use era5cds_forcingMod + use LIS_forecastMod + use LIS_ran2_gasdev + use LIS_histDataMod + + implicit none + +! !ARGUMENTS: + integer, intent(in):: n + integer, intent(in):: findex +! +! !DESCRIPTION: +! Temporally interpolates the forcing data to the current model +! timestep. Downward shortwave radiation is interpolated using a +! zenith-angled based approach. Precipitation and downward longwave +! are not temporally interpolated, and the 1 hourly value is used. +! All other variables, including Downward shortwave radiation are +! linearly interpolated between the 1 hourly blocks. +! +! The routines invoked are: +! \begin{description} +! \item[LIS\_time2date](\ref{LIS_time2date}) \newline +! converts the time to a date format +! \item[zterp](\ref{zterp}) \newline +! zenith-angle based interpolation +! \end{description} +!EOP + integer :: k,kk,zdoy + real :: zw1, zw2 + real :: czm, cze, czb + real :: wt1, wt2 + real :: gmt1, gmt2 + integer :: t,index1 + integer :: bdoy,byr,bmo,bda,bhr,bmn + real*8 :: btime + integer :: status + integer :: mfactor,m + type(ESMF_Field) :: tmpField,q2Field,uField,vField,swdField,lwdField + type(ESMF_Field) :: psurfField,pcpField,cpcpField + real,pointer :: tmp(:),q2(:),uwind(:),vwind(:) + real,pointer :: swd(:),lwd(:),psurf(:),pcp(:),cpcp(:) + + external :: zterp + + btime=era5cds_struc(n)%era5cdstime1 + call LIS_time2date(btime,bdoy,gmt1,byr,bmo,bda,bhr,bmn) + btime=era5cds_struc(n)%era5cdstime2 + call LIS_time2date(btime,bdoy,gmt2,byr,bmo,bda,bhr,bmn) + +!=== Interpolate Data in time + wt1=(era5cds_struc(n)%era5cdstime2-LIS_rc%time)/ & + (era5cds_struc(n)%era5cdstime2-era5cds_struc(n)%era5cdstime1) + wt2=1.0-wt1 + + call ESMF_StateGet(LIS_FORC_Base_State(n,findex),LIS_FORC_Tair%varname(1), & + tmpField,& + rc=status) + call LIS_verify(status, 'Error: Enable Tair in the forcing variables list') + + call ESMF_StateGet(LIS_FORC_Base_State(n,findex),LIS_FORC_Qair%varname(1), & + q2Field,& + rc=status) + call LIS_verify(status, 'Error: Enable Qair in the forcing variables list') + + call ESMF_StateGet(LIS_FORC_Base_State(n,findex),LIS_FORC_SWdown%varname(1),& + swdField,& + rc=status) + call LIS_verify(status, 'Error: Enable SWdown in the forcing variables list') + + call ESMF_StateGet(LIS_FORC_Base_State(n,findex),LIS_FORC_LWdown%varname(1),& + lwdField,& + rc=status) + call LIS_verify(status, 'Error: Enable LWdown in the forcing variables list') + + call ESMF_StateGet(LIS_FORC_Base_State(n,findex),LIS_FORC_Wind_E%varname(1),& + uField,& + rc=status) + call LIS_verify(status, 'Error: Enable Wind_E in the forcing variables list') + + call ESMF_StateGet(LIS_FORC_Base_State(n,findex),LIS_FORC_Wind_N%varname(1),& + vField,& + rc=status) + call LIS_verify(status, 'Error: Enable Wind_N in the forcing variables list') + + call ESMF_StateGet(LIS_FORC_Base_State(n,findex),LIS_FORC_Psurf%varname(1), & + psurfField,& + rc=status) + call LIS_verify(status, 'Error: Enable Psurf in the forcing variables list') + + call ESMF_StateGet(LIS_FORC_Base_State(n,findex),LIS_FORC_Rainf%varname(1), & + pcpField,& + rc=status) + call LIS_verify(status, 'Error: Enable Rainf in the forcing variables list') + + call ESMF_StateGet(LIS_FORC_Base_State(n,findex),LIS_FORC_CRainf%varname(1),& + cpcpField,& + rc=status) + call LIS_verify(status, 'Error: Enable CRainf in the forcing variables list') + + call ESMF_FieldGet(tmpField,localDE=0,farrayPtr=tmp,rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(q2Field,localDE=0,farrayPtr=q2,rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(swdField,localDE=0,farrayPtr=swd,rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(lwdField,localDE=0,farrayPtr=lwd,rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(uField,localDE=0,farrayPtr=uwind,rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(vField,localDE=0,farrayPtr=vwind,rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(psurfField,localDE=0,farrayPtr=psurf,rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(pcpField,localDE=0,farrayPtr=pcp,rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(cpcpField,localDE=0,farrayPtr=cpcp,rc=status) + call LIS_verify(status) + + mfactor = LIS_rc%nensem(n)/era5cds_struc(n)%nIter + +!----------------------------------------------------------------------- +! SW down +!----------------------------------------------------------------------- + do k=1,LIS_rc%ntiles(n)/mfactor + do m=1,mfactor + t = m + (k-1)*mfactor + index1 = LIS_domain(n)%tile(t)%index + zdoy=LIS_rc%doy + call zterp(0,LIS_domain(n)%grid(index1)%lat,& + LIS_domain(n)%grid(index1)%lon,& + gmt1,gmt2,LIS_rc%gmt,zdoy,zw1,zw2,czb,cze,czm,LIS_rc) + + kk = LIS_get_iteration_index(n, k, index1, mfactor) + + if (era5cds_struc(n)%metdata1(kk,3,index1).ne.LIS_rc%udef) then + + swd(t) = zw1 * era5cds_struc(n)%metdata1(kk,3,index1) + + if (swd(t).gt.LIS_CONST_SOLAR) then + write(unit=LIS_logunit,fmt=*) & + '[WARN] sw radiation too high in ERA5!!' + write(unit=LIS_logunit,fmt=*)'[WARN] it is',swd(t), zw1 + write(unit=LIS_logunit,fmt=*)'[WARN] at ',t, & + LIS_domain(n)%grid(index1)%lat, & + LIS_domain(n)%grid(index1)%lon + write(unit=LIS_logunit,fmt=*)'[WARN] era5cdsdata1=',& + era5cds_struc(n)%metdata1(kk,3,index1) + write(unit=LIS_logunit,fmt=*)'[WARN] era5cdsdata2=',& + era5cds_struc(n)%metdata2(kk,3,index1) + swd(t) = LIS_CONST_SOLAR + write(unit=LIS_logunit,fmt=*)'[WARN] forcing set to ',swd(t) + endif + endif + + if ((swd(t).ne.LIS_rc%udef).and.(swd(t).lt.0)) then + ! tiny negative swd appear in the high latitudes sometimes + ! WARN in stead of ERR and endrun -- 2/2/2026 HKB + swd(t) = 0.0 + write(LIS_logunit,*) & + '[WARN] timeinterp_era5cds -- forcing not udef but lt0,' + write(LIS_logunit,*)'[WARN] timeinterp_era5cds -- ', & + t,swd(t),era5cds_struc(n)%metdata2(kk,3,index1), & + ' (',LIS_localPet,')' + !call LIS_endrun + endif + + enddo + enddo + +!----------------------------------------------------------------------- +! precip variable - constant rate over the hour +!----------------------------------------------------------------------- + do k=1,LIS_rc%ntiles(n)/mfactor + do m=1,mfactor + t = m + (k-1)*mfactor + index1 = LIS_domain(n)%tile(t)%index + kk = LIS_get_iteration_index(n, k, index1, mfactor) + pcp(t) = era5cds_struc(n)%metdata1(kk,8,index1) + cpcp(t) = era5cds_struc(n)%metdata1(kk,9,index1) + + if ( pcp(t) < 0 ) then + pcp(t) = 0 + endif + if ( cpcp(t) < 0 ) then + cpcp(t) = 0 + endif + enddo + enddo + +!----------------------------------------------------------------------- +! LW down +!----------------------------------------------------------------------- + do k=1,LIS_rc%ntiles(n)/mfactor + do m=1,mfactor + t = m + (k-1)*mfactor + index1 = LIS_domain(n)%tile(t)%index + kk = LIS_get_iteration_index(n, k, index1, mfactor) + ! Longwave down shouldn't be interpolated between ERA5CDS times + lwd(t) = era5cds_struc(n)%metdata1(kk,4,index1) + enddo + enddo +!----------------------------------------------------------------------- +! Linearly interpolate everything else +!----------------------------------------------------------------------- + + do k=1,LIS_rc%ntiles(n)/mfactor + do m=1,mfactor + t = m + (k-1)*mfactor + index1 = LIS_domain(n)%tile(t)%index + kk = LIS_get_iteration_index(n, k, index1, mfactor) + + tmp(t) = era5cds_struc(n)%metdata1(kk,1,index1)*wt1 + & + era5cds_struc(n)%metdata2(kk,1,index1)*wt2 + q2(t) = era5cds_struc(n)%metdata1(kk,2,index1)*wt1 + & + era5cds_struc(n)%metdata2(kk,2,index1)*wt2 + uwind(t) = era5cds_struc(n)%metdata1(kk,5,index1)*wt1+& + era5cds_struc(n)%metdata2(kk,5,index1)*wt2 + vwind(t) = era5cds_struc(n)%metdata1(kk,6,index1)*wt1+& + era5cds_struc(n)%metdata2(kk,6,index1)*wt2 + psurf(t) = era5cds_struc(n)%metdata1(kk,7,index1)*wt1 + & + era5cds_struc(n)%metdata2(kk,7,index1)*wt2 + + enddo + enddo + +end subroutine timeinterp_era5cds + diff --git a/lis/plugins/LIS_metforcing_pluginMod.F90 b/lis/plugins/LIS_metforcing_pluginMod.F90 index bd1643cf54..ee8b843573 100644 --- a/lis/plugins/LIS_metforcing_pluginMod.F90 +++ b/lis/plugins/LIS_metforcing_pluginMod.F90 @@ -144,6 +144,10 @@ subroutine LIS_metforcing_plugin use era5_forcingMod #endif +#if ( defined MF_ERA5CDS ) + use era5cds_forcingMod +#endif + #if ( defined MF_GSWP1 ) use gswp1_forcingMod #endif @@ -406,6 +410,13 @@ subroutine LIS_metforcing_plugin external reset_era5 #endif +#if ( defined MF_ERA5CDS ) + external get_era5cds + external timeinterp_era5cds + external finalize_era5cds + external reset_era5cds +#endif + #if ( defined MF_GSWP1 ) external get_gswp1 external timeinterp_gswp1 @@ -835,6 +846,16 @@ subroutine LIS_metforcing_plugin call registerfinalmetforc(trim(LIS_era5Id)//char(0),finalize_era5) #endif +#if ( defined MF_ERA5CDS) +! - ERA5 Reanalysis from Climate Data Store Forcing: + call registerinitmetforc(trim(LIS_era5cdsId)//char(0),init_ERA5CDS) + call registerretrievemetforc(trim(LIS_era5cdsId)//char(0),get_era5cds) + call registertimeinterpmetforc(trim(LIS_era5cdsId)//char(0), & + timeinterp_era5cds) + call registerresetmetforc(trim(LIS_era5cdsId)//char(0),reset_era5cds) + call registerfinalmetforc(trim(LIS_era5cdsId)//char(0),finalize_era5cds) +#endif + #if ( defined MF_GSWP1 ) ! - GWSP1 Forcing: call registerinitmetforc(trim(LIS_gswp1Id)//char(0),init_GSWP1) diff --git a/lis/plugins/LIS_pluginIndices.F90 b/lis/plugins/LIS_pluginIndices.F90 index 0ca1602f44..c64b8ccdc2 100755 --- a/lis/plugins/LIS_pluginIndices.F90 +++ b/lis/plugins/LIS_pluginIndices.F90 @@ -205,6 +205,7 @@ module LIS_pluginIndices "AWRAL processed forcing" character*50, public, parameter :: LIS_mrmsId = "MRMS" character*50, public, parameter :: LIS_era5Id = "ERA5" + character*50, public, parameter :: LIS_era5cdsId = "ERA5CDS" character*50, public, parameter :: LIS_plumber2Id = "PLUMBER2" character*50, public, parameter :: LIS_gddpId = "GDDP" character*50, public, parameter :: LIS_galwemId = &