module clsm_ensdrv_force_routines ! collection of *forcing* subroutines for enkf_driver ! (originally these routines were in clsm_ensdrv_drv_routines.F90) ! reichle, 13 Aug 2008 use clsm_ensdrv_glob_param, ONLY: & logit, & logunit, & nodata_generic, & nodata_tolfrac_generic, & nodata_tol_generic use MAPL_ConstantsMod, ONLY: & stefan_boltzmann => MAPL_STFBOL, & Tzero => MAPL_TICE use MAPL_SatVaporMod, ONLY: & MAPL_EQsat use driver_types, ONLY: & met_force_type, & assignment (=) use tile_coord_types, ONLY: & tile_coord_type use date_time_util, ONLY: & date_time_type, & augment_date_time, & datetime_lt_refdatetime, & datetime_le_refdatetime, & is_leap_year, & get_dofyr_pentad use ldas_exceptions, ONLY: & ldas_abort, & LDAS_GENERIC_ERROR USE update_model_paras, ONLY : & upd_params, SYN_PRCP, LAT_SHIFT, LON_SHIFT implicit none include 'netcdf.inc' ! everything is private by default unless made public private public :: get_forcing public :: repair_forcing !!public :: ignore_SWNET_for_snow !! !!logical :: ignore_SWNET_for_snow ! fixes sibalb bug in MERRA snow albedo ! ! MERRA had a snow albedo bug in sibalb. Earlier versions of LDASsa ! included a variable "ignore_SWnet_for_snow" which was meant to address ! this bug but never worked as intended in the LDASsa MPI parallel ! implementation because of a bug in LDASsa. ! - reichle, 2 Apr 2015 real, parameter :: DEFAULT_REFH = 10. ! m real, parameter :: SWDN_MAX = 1360. ! W/m2 real, parameter :: LWDN_EMISS_MIN = 0.5 ! min effective emissivity for LWdown real, parameter :: LWDN_EMISS_MAX = 1.0 ! max effective emissivity for LWdown character(10), private :: tmpstring10 character(40), private :: tmpstring40 contains ! ******************************************************************** subroutine get_forcing( date_time, force_dtstep, met_path, met_tag, & N_catd, tile_coord, met_hinterp, & met_force_obs_tile_new, move_met_force_obs_new_to_old, & init, alb_from_SWnet ) ! Read and check meteorological forcing data for the domain. ! ! time convention: ! - forcing states (such as Tair) are snapshots at date_time ! - forcing fluxes (such as SWdn) are time avg over *subsequent* forcing ! interval (date_time:date_time+force_dtstep) ! ! The above time convention is heritage from older versions of the ! off-line driver and creates problems with "operational" forcing ! data from GEOS5. For "operational" integrations, the forward-looking ! forcing fluxes are not available for "met_force_obs_tile_new". ! ! As a work-around, the output parameter "move_met_force_obs_new_to_old" ! is used to treat "operational" forcing data sets accordingly in the ! main program. This work-around replaces an older work-around that ! was less efficient. ! ! When LDASsa is integrated within the coupled GEOS5 DAS, initial (time-avg) ! "tavg1_2d_*_Nx" files are not available. Use optional "init" flag to ! deal with this situation. ! ! reichle, 28 March 2006 ! reichle, 13 March 2008 - added optional "init" flag ! qliu+reichle, 12 Aug 2008 - new field RefH (reference height) in met_force_type ! reichle, 25 Sep 2009 - removed unneeded inputs ! reichle, 23 Feb 2016 - new and more efficient work-around to make GEOS-5 ! forcing work with LDASsa time convention for forcing data implicit none type(date_time_type), intent(in) :: date_time integer, intent(in) :: force_dtstep character(200), intent(in) :: met_path character(80), intent(in) :: met_tag integer, intent(in) :: N_catd, met_hinterp type(tile_coord_type), dimension(:), pointer :: tile_coord ! input type(met_force_type), dimension(N_catd), intent(out) :: & met_force_obs_tile_new type(met_force_type), intent(out) :: move_met_force_obs_new_to_old logical, intent(in), optional :: init, alb_from_SWnet ! local variables real :: nodata_forcing type(date_time_type) :: date_time_tmp logical :: MERRA_file_specs, GEOS_forcing character(len=*), parameter :: Iam = 'get_forcing' character(len=400) :: err_msg ! -------------------------------------------------------------- ! ! shift forcing date if so indicated by met_tag (for twin experiments, ! see function shift_forcing_date for details) - reichle, 6 Apr 2007 date_time_tmp = shift_forcing_date(met_tag, date_time) ! set reference height to default value (if appropriate, will be overwriten ! within specific subroutine) met_force_obs_tile_new%RefH = DEFAULT_REFH ! set SWnet, PARdrct, PARdffs to nodata_generic ! (Note that nodata_forcing is set to the native no-data-value ! in the individual get_*() subroutines and used to communicate with ! check_forcing_nodata. AFTER the call to check_forcing_nodata all forcing ! fields EXCEPT SWnet must NOT be no-data values, and SWnet must be ! nodata_generic if unavailable.) ! ! reichle+qliu, 8 Oct 2008 ! reichle, 23 Feb 2009 -- same goes for ParDrct, ParDffs ! reichle, 5 Mar 2009 -- deleted ParDrct, ParDffs after testing found no impact ! reichle, 22 Jul 2010 -- fixed treatment of SWnet nodata values ! reichle, 20 Dec 2011 -- reinstated PARdrct and PARdffs for MERRA-Land file specs ! --------------------------------------------------------------------------------- ! ! initialize MERRA_file_specs = .false. GEOS_forcing = .false. move_met_force_obs_new_to_old = 0. met_force_obs_tile_new%SWnet = nodata_generic met_force_obs_tile_new%PARdrct = nodata_generic met_force_obs_tile_new%PARdffs = nodata_generic ! --------------------------------------------------------------------------------- ! ! get forcing in tile space if (index(met_tag, 'Berg_netcdf')/=0) then call get_Berg_netcdf( date_time_tmp, met_path, N_catd, tile_coord, & met_force_obs_tile_new, nodata_forcing) ! check for nodata values and unphysical values ! (check here, not outside "if" block, because of GEOSgcm case) call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, & met_force_obs_tile_new ) call repair_forcing( N_catd, met_force_obs_tile_new, & echo=.true., tile_coord=tile_coord ) elseif (index(met_tag, 'GLDAS_2x2_5_netcdf')/=0) then call get_GLDAS_2x2_5_netcdf(date_time_tmp, met_path, N_catd, tile_coord, & met_force_obs_tile_new, nodata_forcing) ! check for nodata values and unphysical values ! (check here, not outside "if" block, because of GEOSgcm case) call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, & met_force_obs_tile_new ) call repair_forcing( N_catd, met_force_obs_tile_new, & echo=.true., tile_coord=tile_coord ) elseif (index(met_tag, 'Viviana_OK')/=0) then ! vmaggion & reichle, 17 July 2008 ! ! use 2x2.5 deg GLDAS for all forcing fields except precip call get_GLDAS_2x2_5_netcdf(date_time_tmp, met_path, N_catd, tile_coord, & met_force_obs_tile_new, nodata_forcing) if (index(met_tag, 'Viviana_OK_nopert')/=0) then call get_Viviana_OK_precip(10, date_time_tmp, met_path, met_tag, & N_catd, tile_coord, met_force_obs_tile_new) end if ! check for nodata values and unphysical values ! (check here, not outside "if" block, because of GEOSgcm case) call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, & met_force_obs_tile_new ) call repair_forcing( N_catd, met_force_obs_tile_new, & echo=.true., tile_coord=tile_coord ) elseif (index(met_tag, 'GSWP2_1x1_netcdf')/=0) then call get_GSWP2_1x1_netcdf( date_time_tmp, met_path, N_catd, tile_coord, & met_force_obs_tile_new, nodata_forcing) ! check for nodata values and unphysical values ! (check here, not outside "if" block, because of GEOSgcm case) call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, & met_force_obs_tile_new ) call repair_forcing( N_catd, met_force_obs_tile_new, & echo=.true., tile_coord=tile_coord ) elseif (index(met_tag, 'RedArk_ASCII')/=0) then call get_RedArk_ASCII( date_time_tmp, met_path, N_catd, tile_coord, & met_force_obs_tile_new, nodata_forcing) ! check for nodata values and unphysical values ! (check here, not outside "if" block, because of GEOSgcm case) call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, & met_force_obs_tile_new ) call repair_forcing( N_catd, met_force_obs_tile_new, & echo=.true., tile_coord=tile_coord ) elseif (index(met_tag, 'RedArk_GOLD')/=0) then call get_RedArk_GOLD( date_time_tmp, met_path, N_catd, tile_coord, & met_force_obs_tile_new, nodata_forcing) ! check for nodata values and unphysical values ! (check here, not outside "if" block, because of GEOSgcm case) call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, & met_force_obs_tile_new ) call repair_forcing( N_catd, met_force_obs_tile_new, & echo=.true., tile_coord=tile_coord ) elseif (index(met_tag, 'RedArk_Princeton')/=0) then call get_RedArk_Princeton( date_time_tmp, met_path, N_catd, tile_coord, & met_force_obs_tile_new, nodata_forcing) ! check for nodata values and unphysical values ! (check here, not outside "if" block, because of GEOSgcm case) call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, & met_force_obs_tile_new ) call repair_forcing( N_catd, met_force_obs_tile_new, & echo=.true., tile_coord=tile_coord ) elseif (index(met_tag, 'Princeton_netcdf')/=0) then ! tyamada+reichle, 17 Jul 2007 call get_Princeton_netcdf( date_time_tmp, met_path, N_catd, tile_coord, & met_force_obs_tile_new, nodata_forcing) ! check for nodata values and unphysical values ! (check here, not outside "if" block, because of GEOSgcm case) call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, & met_force_obs_tile_new ) call repair_forcing( N_catd, met_force_obs_tile_new, & echo=.true., tile_coord=tile_coord ) elseif (index(met_tag, 'conus_0.5d_netcdf')/=0) then ! sarith+reichle, 17 Jul 2007 call get_conus_netcdf( date_time_tmp, met_path, N_catd, tile_coord, & met_force_obs_tile_new, nodata_forcing) ! check for nodata values and unphysical values ! (check here, not outside "if" block, because of GEOSgcm case) call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, & met_force_obs_tile_new ) call repair_forcing( N_catd, met_force_obs_tile_new, & echo=.true., tile_coord=tile_coord ) else ! assume forcing from GEOS5 GCM ("DAS" or "MERRA") output if (logit) write (logunit,*) 'get_forcing(): assuming GEOS-5 forcing data set' GEOS_forcing = .true. ! note "met_tag" in call to get_GEOSgcm_gfio (interface differs ! from other get_* subroutines) !call get_GEOSgcm_gfio( date_time_tmp, met_path, met_tag, & ! N_catd, tile_coord, & ! met_force_obs_tile_new, nodata_forcing) call get_GEOS( date_time_tmp, force_dtstep, & met_path, met_tag, N_catd, tile_coord, met_hinterp, & met_force_obs_tile_new, nodata_forcing, MERRA_file_specs, init ) ! check for nodata values and unphysical values ! (check here, not outside "if" block, because of GEOSgcm case) call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, & met_force_obs_tile_new ) ! call repair_forcing with switch "unlimited_Qair=.true." for GEOS5 forcing ! (default is to limit Qair so that it does not exceed Qair_sat) ! reichle+qliu, 8 Oct 2008 ! ! likewise for "unlimited_LWdown=.true." ! reichle, 11 Feb 2009 call repair_forcing( N_catd, met_force_obs_tile_new, & echo=.true., tile_coord=tile_coord, & unlimited_Qair=.true., unlimited_LWdown=.true. ) ! Subroutine get_GEOS() reads forcing fluxes from "previous" ! interval, not from "subsequent" interval, because in operational ! applications the "subsequent" fluxes for "met_force_new" are not ! available. The following lines restore consistency with the ! time convention stated above. Note that only "old" fluxes ! are needed in subroutine interpolate_to_timestep(), and ! "met_force_obs_tile_new" is set to nodata for forcing fluxes. move_met_force_obs_new_to_old%Rainf_C = 1. move_met_force_obs_new_to_old%Rainf = 1. move_met_force_obs_new_to_old%Snowf = 1. move_met_force_obs_new_to_old%LWdown = 1. move_met_force_obs_new_to_old%SWdown = 1. move_met_force_obs_new_to_old%SWnet = 1. move_met_force_obs_new_to_old%PARdrct = 1. move_met_force_obs_new_to_old%PARdffs = 1. if (MERRA_file_specs) move_met_force_obs_new_to_old%Wind = 1. end if ! make sure SWnet is generally available if needed to back out albedo ! (only works for GEOS forcing, e.g., MERRA, FP, FP-IT) ! ! NOTE: need to check here because GEOS forcing is the default ! forcing data set in the "if... elseif... elseif... else..." ! statement above, that is, it is only asserted here whether ! the requested forcing data are GEOS. if (present(alb_from_SWnet)) then if ( (.not. GEOS_forcing) .and. (alb_from_SWnet) ) then ! stop if per nml inputs the albedo should be backed ! out from SWnet but forcing data is not GEOS err_msg = 'requested nml input [alb_from_SWnet=.true.] ' // & '*only* works with GEOS forcing' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if end if ! stop if anything other than nearest-neighbor met forcing interpolation was ! requested for a non-GEOS5 forcing dataset if ((.not. GEOS_forcing) .and. (met_hinterp>0)) then err_msg = 'for non-GEOS forcing, only ' // & 'nearest-neighbor interpolation is available' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if end subroutine get_forcing ! **************************************************************** subroutine get_Berg_netcdf(date_time, met_path, N_catd, tile_coord, & met_force_new, nodata_forcing ) ! read Berg_NetCDF files and extract forcings in tile space ! (uses nearest neighbor interpolation) ! reichle, 25 May 2005 ! reichle, 23 Feb 2009 -- revised treatment of Rainf_C implicit none type(date_time_type), intent(in) :: date_time character(200), intent(in) :: met_path integer, intent(in) :: N_catd type(tile_coord_type), dimension(:), pointer :: tile_coord ! input type(met_force_type), dimension(N_catd), intent(inout) :: met_force_new real, intent(out) :: nodata_forcing ! Berg grid and netcdf parameters integer, parameter :: berg_grid_N_lon = 720 integer, parameter :: berg_grid_N_lat = 360 real, parameter :: berg_grid_ll_lon = -180. real, parameter :: berg_grid_ll_lat = -90. real, parameter :: berg_grid_dlon = .5 real, parameter :: berg_grid_dlat = .5 integer, parameter :: N_berg_compressed = 67420 ! Berg forcing time step in hours integer, parameter :: dt_berg_in_hours = 6 integer, parameter :: nciv_land_i = 3 integer, parameter :: nciv_land_j = 4 integer, parameter :: nciv_data = 7 integer, parameter :: N_berg_vars = 8 real, parameter :: nodata_berg = 1.e20 character(40), dimension(N_berg_vars), parameter :: berg_dir = (/ & 'PREC-CONV/', & 'PREC-TOTL/', & 'PRES-SRF/ ', & 'RAD-LW/ ', & 'RAD-SW/ ', & 'TEMP-AIR/ ', & 'TEMP-DEW/ ', & 'WIND/ ' /) character(40), dimension(N_berg_vars), parameter :: berg_name = (/ & 'RainfSnowf_C_ecmwf', & 'RainfSnowf_ecmwf ', & 'PSurf_ecmwf ', & 'LWdown_ecmwf ', & 'SWdown_ecmwf ', & 'Tair_ecmwf ', & 'Tdew_ecmwf ', & 'Wind_ecmwf ' /) ! local variables !!real, parameter :: Tzero = 273.16 real :: tol real, dimension(berg_grid_N_lon,berg_grid_N_lat) :: tmp_grid integer, dimension(N_berg_compressed) :: land_i_berg, land_j_berg integer, dimension(N_catd) :: i_ind, j_ind real, dimension(N_berg_compressed) :: tmp_vec real, dimension(N_catd,N_berg_vars) :: force_array integer, dimension(2) :: start, icount integer :: k, n, hours_in_month, berg_var, ierr, ncid real :: this_lon, this_lat, dt_berg_in_seconds character(4) :: YYYY character(2) :: MM character(300) :: fname character(len=*), parameter :: Iam = 'get_Berg_netcdf' character(len=400) :: err_msg ! -------------------------------------------------------------------- dt_berg_in_seconds = real(3600*dt_berg_in_hours) nodata_forcing = nodata_berg tol = abs(nodata_forcing*nodata_tolfrac_generic) ! assemble year and month strings write (YYYY,'(i4.4)') date_time%year write (MM, '(i2.2)') date_time%month ! find out which data are needed ! compressed space dimension (always read global vector) start(1) = 1 icount(1) = N_berg_compressed ! time dimension (first entry in Berg_NetCDF file is at 0Z) if ( (date_time%min/=0) .or. (date_time%sec/=0) .or. & (mod(date_time%hour,dt_berg_in_hours)/=0) ) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'timing error') end if hours_in_month = (date_time%day-1)*24 + date_time%hour start(2) = hours_in_month / dt_berg_in_hours + 1 icount(2) = 1 ! ---------------------------------------------- ! ! compute indices for nearest neighbor interpolation from Berg grid ! to tile space ! ! (NOTE: this should at some point be replaced with a regridding ! subroutine that interpolates from the ! native forcing grid to the GCM atmospheric grid that is used ! to cut catchments into tiles - then "standard" grid2tile ! using tile_coord%atm_i and tile_coord%atm_j applies. ! reichle, 26 May 2005) do k=1,N_catd ! ll_lon and ll_lat refer to lower left corner of grid cell ! (as opposed to the grid point in the center of the grid cell) this_lon = tile_coord(k)%com_lon this_lat = tile_coord(k)%com_lat i_ind(k) = ceiling( (this_lon - berg_grid_ll_lon)/berg_grid_dlon ) j_ind(k) = ceiling( (this_lat - berg_grid_ll_lat)/berg_grid_dlat ) ! NOTE: For a "date line on center" grid and (180-dlon/2) < lon < 180 ! we now have i_ind=(grid%N_lon+1) ! This would need to be fixed. if (i_ind(k)>berg_grid_N_lon) i_ind(k)=1 end do ! ------------------------------------------------------ ! ! read compression parameters (same for all data variables and time steps) berg_var = 1 fname = trim(met_path) // trim(berg_dir(berg_var)) // '/' // YYYY & // '/' // trim(berg_name(berg_var)) // '.' // YYYY // MM // '.nc' if (logit) write (logunit,*) 'get netcdf compression params from ' // trim(fname) ierr = NF_OPEN(fname,NF_NOWRITE,ncid) if (ierr/=0) then err_msg = 'error opening netcdf file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if ierr = NF_GET_VARA_INT(ncid, nciv_land_i, start, icount, land_i_berg) ierr = NF_GET_VARA_INT(ncid, nciv_land_j, start, icount, land_j_berg) ierr = NF_CLOSE(ncid) ! ------------------------------------------------------ ! ! get forcing data do berg_var = 1,N_berg_vars ! open file, read compressed data, and put on global grid fname = trim(met_path) // trim(berg_dir(berg_var)) // '/' // YYYY & // '/' // trim(berg_name(berg_var)) // '.' // YYYY // MM // '.nc' if (logit) write (logunit,*) 'opening ' // trim(fname) ierr = NF_OPEN(fname,NF_NOWRITE,ncid) if (ierr/=0) then err_msg = 'error opening netcdf file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if ierr = NF_GET_VARA_REAL(ncid, nciv_data, start, icount, tmp_vec ) ierr = NF_CLOSE(ncid) tmp_grid = nodata_forcing do n=1,N_berg_compressed tmp_grid(land_i_berg(n), land_j_berg(n) ) = tmp_vec(n) end do ! interpolate to tile space ! (NOTE: This should at some point be replaced with a regridding ! subroutine that interpolates from the ! native forcing grid to the GCM atmospheric grid that is used ! to cut catchments into tiles - then "standard" grid2tile ! using tile_coord%atm_i and tile_coord%atm_j applies. ! reichle, 26 May 2005) do k=1,N_catd force_array(k,berg_var) = tmp_grid(i_ind(k), j_ind(k)) end do end do ! convert variables and units of force_array to match met_force_type, ! put into structure ! from Berg files: ! ! force_array(:,1) = PREC-CONV = Rainf_C+Snowf_C kg/m2 (6h total) ??? ! force_array(:,2) = PREC-TOTL = RainfSnowf kg/m2 (6h total) ! force_array(:,3) = PRES-SRF = PSurf Pa ! force_array(:,4) = RAD-LW = LWdown W/m2 ! force_array(:,5) = RAD-SW = SWdown W/m2 ! force_array(:,6) = TEMP-AIR = Tair K ! force_array(:,7) = TEMP-DEW = Tdew K ! force_array(:,8) = WIND = Wind m/s met_force_new%Psurf = force_array(:,3) met_force_new%LWdown = force_array(:,4) met_force_new%SWdown = force_array(:,5) met_force_new%Tair = force_array(:,6) met_force_new%Wind = force_array(:,8) ! get specific humidity from dew point temperature do k=1,N_catd if ( abs(force_array(k,3)-nodata_berg) DIFFERENT FROM GSWP2 FORMAT!!! <------ ! ! At 0Z for first day of month: ! - for fluxes read first entry of that month ! - for states read first entry of that month ! At 6Z for first day of month: ! - for fluxes read second entry of that month ! - for states read second entry of that month ! and so on... select case (this_var) case (1,2,3,4,5) ! "fluxes" date_time_tmp = date_time case (6,7,8,9) ! "states" date_time_tmp = date_time !!call augment_date_time(-int_dt_RedArk_GOLD_in_seconds,date_time_tmp) case default call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error') end select ! assemble year and month strings write (YYYY,'(i4.4)') date_time_tmp%year write (MM, '(i2.2)') date_time_tmp%month write (DD, '(i2.2)') date_time_tmp%day write (HHMM,'(i4.4)') date_time_tmp%hour*100 + date_time_tmp%min ! assemble file name, open file fname = trim(met_path) // '/' // trim(RedArk_GOLD_name(this_var)) & // '/Y' // YYYY // '/M' // MM // '/' & // trim(RedArk_GOLD_name(this_var)) // '_RedArk_' // & YYYY // MM // DD // '_' // HHMM if (logit) write (logunit,*) 'opening ' // trim(fname) open(10,file=fname,form='formatted',action='read') ! read compressed data, and put into tile space do i=1,RedArk_GOLD_N_cells read (10,*) tmp_vec(i) end do close(10,status='keep') do k=1,N_catd force_array(k,this_var) = tmp_vec( ind_gold2tile(k) ) end do end do ! convert variables and units of force_array to match met_force_type, ! put into structure ! from RedArk GOLD files: ! ! force_array(:, 1) = SWdown W/m2 ! force_array(:, 2) = LWdown W/m2 ! force_array(:, 3) = Rainf_C kg/m2/s ! force_array(:, 4) = Rainf kg/m2/s ! force_array(:, 5) = Snowf kg/m2/s ! force_array(:, 6) = PSurf Pa ! force_array(:, 7) = Qair kg/kg ! force_array(:, 8) = Tair K ! force_array(:, 9) = Wind m/s met_force_new%SWdown = force_array(:,1) met_force_new%LWdown = force_array(:,2) met_force_new%Rainf_C = force_array(:,3) met_force_new%Rainf = force_array(:,4) met_force_new%Snowf = force_array(:,5) met_force_new%Psurf = force_array(:,6) met_force_new%Qair = force_array(:,7) met_force_new%Tair = force_array(:,8) met_force_new%Wind = force_array(:,9) end subroutine get_RedArk_GOLD ! *************************************************************************** subroutine get_RedArk_Princeton(date_time, met_path, N_catd, tile_coord, & met_force_new, nodata_forcing ) ! read RedArk Princeton files and map to forcings in tile space ! (uses pre-computed nearest neighbor interpolation) ! reichle, 6 Apr 2007 implicit none type(date_time_type), intent(in) :: date_time character(200), intent(in) :: met_path integer, intent(in) :: N_catd type(tile_coord_type), dimension(:), pointer :: tile_coord ! input type(met_force_type), dimension(N_catd), intent(inout) :: met_force_new real, intent(out) :: nodata_forcing ! interpolation parameters integer, parameter :: RedArk_Princeton_N_cells = 70 character(80), parameter :: RedArk_Princeton_coord_file = & 'Princeton_forcing_cells_2_RedArkOSSE_tiles.dat' ! Princeton_REDARK forcing time step in hours integer, parameter :: dt_RedArk_Princeton_in_hours = 3 integer, parameter :: N_RedArk_Princeton_vars = 7 real, parameter :: nodata_RedArk_Princeton = -9999. ! ????? character(40), dimension(N_RedArk_Princeton_vars), parameter :: & RedArk_Princeton_name = (/ & 'dswrf', & ! 1 - flux 'dlwrf', & ! 2 - flux 'prcp ', & ! 3 - flux 'pres ', & ! 4 - state 'shum ', & ! 5 - state 'tas ', & ! 6 - state 'wind ' /) ! 7 - state ! local variables real :: tol integer, dimension(N_catd) :: ind_princeton2tile real, dimension(RedArk_Princeton_N_cells) :: tmp_vec real, dimension(N_catd,N_RedArk_Princeton_vars) :: force_array type(date_time_type) :: date_time_tmp integer :: int_dt_RedArk_Princeton_in_secs, i, k, this_var, tmp_tile_id character(4) :: YYYY, HHMM character(2) :: MM, DD character(300) :: fname character(len=*), parameter :: Iam = 'get_RedArk_Princeton' character(len=400) :: err_msg ! -------------------------------------------------------------------- int_dt_RedArk_Princeton_in_secs = 3600*dt_RedArk_Princeton_in_hours nodata_forcing = nodata_RedArk_Princeton tol = abs(nodata_forcing*nodata_tolfrac_generic) ! ---------------------------------------------- ! ! load indices for nearest neighbor interpolation from Princeton RedArk vec ! to tile space fname = trim(met_path) // '/' // trim(RedArk_Princeton_coord_file) open (10, file=fname, form='formatted', action='read') do i=1,N_catd read (10,*) tmp_tile_id, ind_princeton2tile(i) if (tmp_tile_id /= tile_coord(i)%tile_id) then err_msg = 'RedArk Princeton2tile error' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if end do ! ------------------------------------------------------ ! ! get forcing data do this_var = 1,N_RedArk_Princeton_vars ! time dimension ! ! *** STATEMENTS BELOW STILL NEED TO BE VERIFIED!!!! - reichle, 6 Apr 2007 **** ! ! First entry in Princeton_RedArk file is at 0Z for states, ! with fluxes for 0Z-3Z ! ------> DIFFERENT FROM GSWP2 FORMAT!!! <------ ! ! At 0Z for first day of month: ! - for fluxes read first entry of that month ! - for states read first entry of that month ! At 3Z for first day of month: ! - for fluxes read second entry of that month ! - for states read second entry of that month ! and so on... select case (this_var) case (1,2,3) ! "fluxes" date_time_tmp = date_time case (4,5,6,7) ! "states" date_time_tmp = date_time !!call augment_date_time(-int_dt_RedArk_Princeton_in_secs,date_time_tmp) case default call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error') end select ! assemble year and month strings write (YYYY,'(i4.4)') date_time_tmp%year write (MM, '(i2.2)') date_time_tmp%month write (DD, '(i2.2)') date_time_tmp%day write (HHMM,'(i4.4)') date_time_tmp%hour*100 + date_time_tmp%min ! assemble file name, open file fname = trim(met_path) // '/' // trim(RedArk_Princeton_name(this_var)) & // '/Y' // YYYY // '/M' // MM // '/' & // trim(RedArk_Princeton_name(this_var)) // '_RedArk_' // & YYYY // MM // DD // '_' // HHMM if (logit) write (logunit,*) 'opening ' // trim(fname) open(10,file=fname,form='formatted',action='read') ! read compressed data, and put into tile space do i=1,RedArk_Princeton_N_cells read (10,*) tmp_vec(i) end do close(10,status='keep') do k=1,N_catd force_array(k,this_var) = tmp_vec( ind_princeton2tile(k) ) end do end do ! convert variables and units of force_array to match met_force_type, ! put into structure ! from RedArk Princeton files: ! ! force_array(:, 1) = SWdown W/m2 ! force_array(:, 2) = LWdown W/m2 ! force_array(:, 3) = RainfSnowf kg/m2/s ! force_array(:, 4) = PSurf Pa ! force_array(:, 5) = Qair kg/kg ! force_array(:, 6) = Tair K ! force_array(:, 7) = Wind m/s met_force_new%SWdown = force_array(:,1) met_force_new%LWdown = force_array(:,2) met_force_new%Psurf = force_array(:,4) met_force_new%Qair = force_array(:,5) met_force_new%Tair = force_array(:,6) met_force_new%Wind = force_array(:,7) do k=1,N_catd ! rain and snow: ! partition total precip into rain and snow according to Tair ! set convective precip to zero met_force_new(k)%Rainf_C = 0. if ( met_force_new(k)%Tair < Tzero ) then met_force_new(k)%Rainf = 0. met_force_new(k)%Snowf = force_array(k,3) else met_force_new(k)%Rainf = force_array(k,3) met_force_new(k)%Snowf = 0. end if end do end subroutine get_RedArk_Princeton ! **************************************************************** subroutine get_Princeton_netcdf(date_time, met_path, N_catd, tile_coord, & met_force_new, nodata_forcing ) ! read Princeton_NetCDF files and extract forcing in tile space ! (uses nearest neighbor interpolation) ! tyamada+reichle, 19 Jul 2007 implicit none type(date_time_type), intent(in) :: date_time character(200), intent(in) :: met_path integer, intent(in) :: N_catd type(tile_coord_type), dimension(:), pointer :: tile_coord ! input type(met_force_type), dimension(N_catd), intent(inout) :: met_force_new real, intent(out) :: nodata_forcing ! Princeton_netcdf grid and netcdf parameters integer, parameter :: Princeton_grid_N_lon = 360 integer, parameter :: Princeton_grid_N_lat = 180 real, parameter :: Princeton_grid_ll_lon = 0. real, parameter :: Princeton_grid_ll_lat = -90. real, parameter :: Princeton_grid_dlon = 1. real, parameter :: Princeton_grid_dlat = 1. ! Princeton_netcdf forcing time step in hours integer, parameter :: dt_Princeton_in_hours = 3 integer, parameter :: nciv_data = 5 ! (1=lon, 2=lat, 3=z, 4=time, 5=data) integer, parameter :: N_Princeton_vars = 7 real, parameter :: nodata_Princeton = 2.e20 character(40), dimension(N_Princeton_vars), parameter :: Princeton_name = & (/ & 'dswrf', & ! 1 - flux 'dlwrf', & ! 2 - flux 'prcp ', & ! 3 - flux 'pres ', & ! 4 - state 'shum ', & ! 5 - state 'tas ', & ! 6 - state 'wind ' & ! 7 - state /) ! local variables integer, dimension(N_catd) :: i_ind, j_ind real, dimension(Princeton_grid_N_lon, Princeton_grid_N_lat) :: tmp_grid real, dimension(N_catd, N_Princeton_vars) :: force_array integer, dimension(4) :: start, icount integer :: k, hours_in_year, Princeton_var, ierr, ncid real :: tol, this_lon, this_lat, dt_Princeton_in_seconds character( 4) :: YYYY, HHMM character( 2) :: MM, DD character(300) :: fname character(len=*), parameter :: Iam = 'get_Princeton_netcdf' character(len=400) :: err_msg ! ---------------------------------------------------------------- dt_Princeton_in_seconds = real(3600*dt_Princeton_in_hours) nodata_forcing = nodata_Princeton tol = abs(nodata_forcing*nodata_tolfrac_generic) ! assemble year and month strings write (YYYY, '(i4.4)') date_time%year write (MM, '(i2.2)') date_time%month write (DD, '(i2.2)') date_time%day write (HHMM, '(i4.4)') date_time%hour*100 + date_time%min ! set lon index start(1) = 1 icount(1) = 360 ! set lat index start(2) = 1 icount(2) = 180 ! set z index start(3) = 1 icount(3) = 1 ! get time index if ( (date_time%min/=0) .or. (date_time%sec/=0) .or. & (mod(date_time%hour, dt_Princeton_in_hours)/=0) ) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'timing ERROR!!') endif hours_in_year = (date_time%dofyr-1)*24 + date_time%hour start(4) = hours_in_year / dt_Princeton_in_hours + 1 icount(4) = 1 ! ---------------------------------------------------------------- do k=1,N_catd this_lon = tile_coord(k)%com_lon this_lat = tile_coord(k)%com_lat ! change lon units for compatibility with Princeton netcdf ! grid which starts at the Greenwich Meridian and goes ! eastward (lon=0:360) if (this_lon<0.) this_lon = this_lon + 360. i_ind(k) = ceiling((this_lon - Princeton_grid_ll_lon)/Princeton_grid_dlon) j_ind(k) = ceiling((this_lat - Princeton_grid_ll_lat)/Princeton_grid_dlat) ! not sure this is quite right -- reichle, 24 Feb 2009 i_ind(k) = mod(i_ind(k), Princeton_grid_N_lon) ! fixed per suggestion from Greg Walker, - reichle, 26 Aug 2015 if(i_ind(k) < 1) i_ind(k) = i_ind(k) + Princeton_grid_N_lon enddo ! ---------------------------------------------------------------- ! ! open input file do Princeton_var = 1,N_Princeton_vars fname = trim(met_path) // '/' // trim(Princeton_name(Princeton_var)) & // '_3hourly_' // YYYY // '-' // YYYY // '.nc' if (logit) write(logunit,*) 'opening' // trim(fname) ierr = NF_OPEN(fname, NF_NOWRITE, ncid) if (ierr/=0) then err_msg = 'error opening netcdf file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if ierr = NF_GET_VARA_REAL(ncid, nciv_data, start, icount, tmp_grid) ierr = NF_CLOSE(ncid) do k = 1, N_catd force_array(k, Princeton_var) = tmp_grid(i_ind(k), j_ind(k)) enddo enddo ! ---------------------------------------------------------------- met_force_new%SWdown = force_array(:, 1) met_force_new%LWdown = force_array(:, 2) met_force_new%Psurf = force_array(:, 4) met_force_new%Qair = force_array(:, 5) met_force_new%Tair = force_array(:, 6) met_force_new%Wind = force_array(:, 7) do k=1, N_catd met_force_new(k)%Rainf_C = 0. if ( met_force_new(k)%Tair < Tzero ) then met_force_new(k)%Rainf = 0. met_force_new(k)%Snowf = force_array(k,3) else met_force_new(k)%Rainf = force_array(k,3) met_force_new(k)%Snowf = 0. endif enddo end subroutine get_Princeton_netcdf ! **************************************************************** subroutine get_conus_netcdf(date_time, met_path, N_catd, tile_coord, & met_force_new, nodata_forcing ) ! read conus_NetCDF files and extract forcing in tile space ! (uses nearest neighbor interpolation) ! sarith+reichle, 19 Jul 2007 implicit none type(date_time_type), intent(in) :: date_time character(200), intent(in) :: met_path integer, intent(in) :: N_catd type(tile_coord_type), dimension(:), pointer :: tile_coord ! input type(met_force_type), dimension(N_catd), intent(inout) :: met_force_new real, intent(out) :: nodata_forcing ! conus_netcdf grid and netcdf parameters integer, parameter :: conus_grid_N_lon = 115 integer, parameter :: conus_grid_N_lat = 48 real, parameter :: conus_grid_ll_lon = -125. real, parameter :: conus_grid_ll_lat = 25. real, parameter :: conus_grid_dlon = 0.5 real, parameter :: conus_grid_dlat = 0.5 ! conus_netcdf forcing time step in hours integer, parameter :: dt_conus_in_hours = 1 integer, parameter :: N_conus_vars = 7 real, parameter :: nodata_conus = 1.e20 character(40), dimension(N_conus_vars), parameter :: conus_name = & (/ & 'fsds ', & ! 1 - flux 'flds ', & ! 2 - flux 'precs', & ! 3 - flux 'tbot ', & ! 4 - state 'wind ', & ! 5 - state 'psrf ', & ! 6 - state 'qbot ' & ! 7 - state /) integer :: nciv_data ! local variables integer, dimension(N_catd) :: i_ind, j_ind real, dimension(conus_grid_N_lon, conus_grid_N_lat) :: tmp_grid real, dimension(N_catd, N_conus_vars) :: force_array integer, dimension(3) :: start, icount integer :: k, hours_in_month, conus_var, ierr, ncid real :: tol, this_lon, this_lat, dt_conus_in_seconds character( 4) :: YYYY, HHMM character( 2) :: MM, DD character(300) :: fname character(len=*), parameter :: Iam = 'get_conus_netcdf' character(len=400) :: err_msg ! ---------------------------------------------------------------- dt_conus_in_seconds = real(3600*dt_conus_in_hours) nodata_forcing = nodata_conus tol = abs(nodata_forcing*nodata_tolfrac_generic) ! assemble year and month strings write (YYYY, '(i4.4)') date_time%year write (MM, '(i2.2)') date_time%month write (DD, '(i2.2)') date_time%day write (HHMM, '(i4.4)') date_time%hour*100 + date_time%min ! set lon index start(1) = 1 icount(1) = conus_grid_N_lon ! set lat index start(2) = 1 icount(2) = conus_grid_N_lat ! set z index !start(3) = 1 !icount(3) = 1 ! get time index if ( (date_time%min/=0) .or. (date_time%sec/=0) .or. & (mod(date_time%hour, dt_conus_in_hours)/=0) ) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'timing ERROR!!') endif hours_in_month = (date_time%day-1)*24 + date_time%hour start(3) = hours_in_month / dt_conus_in_hours + 1 icount(3) = 1 ! ---------------------------------------------------------------- do k=1,N_catd this_lon = tile_coord(k)%com_lon this_lat = tile_coord(k)%com_lat i_ind(k) = ceiling((this_lon - conus_grid_ll_lon)/conus_grid_dlon) j_ind(k) = ceiling((this_lat - conus_grid_ll_lat)/conus_grid_dlat) enddo ! ---------------------------------------------------------------- ! ! open input file fname = trim(met_path) // '/' // YYYY//'-'//MM//'.nc' if (logit) write (logunit,*) 'opening' // trim(fname) ierr = NF_OPEN(fname, NF_NOWRITE, ncid) if (ierr/=0) then err_msg = 'error opening conus_netcdf file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if do conus_var = 1,N_conus_vars ! nciv_data = order of fields in netcdf file: ! ! 5=lat, 6=lon, 7=fsds, 8=flds, 9=precs, 10=tbot, 11=wind, 12=psrf, 13=qbot nciv_data = 6 + conus_var ierr = NF_GET_VARA_REAL(ncid, nciv_data , start, icount, tmp_grid) do k = 1, N_catd force_array(k, conus_var) = tmp_grid(i_ind(k), j_ind(k)) end do enddo ierr = NF_CLOSE(ncid) ! ---------------------------------------------------------------- met_force_new%SWdown = force_array(:, 1) met_force_new%LWdown = force_array(:, 2) met_force_new%Psurf = force_array(:, 6) met_force_new%Qair = force_array(:, 7) met_force_new%Tair = force_array(:, 4) met_force_new%Wind = force_array(:, 5) do k=1, N_catd met_force_new(k)%Rainf_C = 0. if ( met_force_new(k)%Tair < Tzero ) then met_force_new(k)%Rainf = 0. met_force_new(k)%Snowf = force_array(k,3) else met_force_new(k)%Rainf = force_array(k,3) met_force_new(k)%Snowf = 0. endif enddo end subroutine get_conus_netcdf ! **************************************************************** subroutine get_GLDAS_2x2_5_netcdf(date_time, met_path, N_catd, tile_coord, & met_force_new, nodata_forcing ) ! read GLDAS_NetCDF files and extract forcings in tile space ! (uses nearest neighbor interpolation) ! reichle, 30 Jun 2005 implicit none type(date_time_type), intent(in) :: date_time character(200), intent(in) :: met_path integer, intent(in) :: N_catd type(tile_coord_type), dimension(:), pointer :: tile_coord ! input type(met_force_type), dimension(N_catd), intent(inout) :: met_force_new real, intent(out) :: nodata_forcing ! GLDAS grid and netcdf parameters integer, parameter :: gldas_grid_N_lon = 144 integer, parameter :: gldas_grid_N_lat = 76 real, parameter :: gldas_grid_ll_lon = -181.25 real, parameter :: gldas_grid_ll_lat = -61. real, parameter :: gldas_grid_dlon = 2.5 real, parameter :: gldas_grid_dlat = 2. integer, parameter :: N_gldas_compressed = 3023 ! GLDAS forcing time step in hours integer, parameter :: dt_gldas_in_hours = 3 integer, parameter :: nciv_land_i = 3 integer, parameter :: nciv_land_j = 4 integer, parameter :: nciv_data = 7 integer, parameter :: N_gldas_vars = 10 real, parameter :: nodata_gldas = 1.e20 character(40), dimension(N_gldas_vars), parameter :: gldas_name = (/ & 'SWdown_gldas ', & ! 1 'LWdown_gldas ', & ! 2 'Snowf_gldas ', & ! 3 'Rainf_gldas ', & ! 4 'Tair_gldas ', & ! 5 'Qair_gldas ', & ! 6 'Wind_E_gldas ', & ! 7 'Wind_N_gldas ', & ! 8 'PSurf_gldas ', & ! 9 'RainfSnowf_C_gldas' /) ! 10 ??? ! local variables real :: tol real, dimension(gldas_grid_N_lon,gldas_grid_N_lat) :: tmp_grid integer, dimension(N_gldas_compressed) :: land_i_gldas, land_j_gldas integer, dimension(N_catd) :: i_ind, j_ind real, dimension(N_gldas_compressed) :: tmp_vec real, dimension(N_catd,N_gldas_vars) :: force_array integer, dimension(2) :: start, icount integer :: k, n, hours_in_month, gldas_var, ierr, ncid real :: this_lon, this_lat, dt_gldas_in_seconds character(4) :: YYYY character(2) :: MM character(300) :: fname character(len=*), parameter :: Iam = 'get_GLDAS_netcdf' character(len=400) :: err_msg ! -------------------------------------------------------------------- dt_gldas_in_seconds = real(3600*dt_gldas_in_hours) nodata_forcing = nodata_gldas tol = abs(nodata_forcing*nodata_tolfrac_generic) ! assemble year and month strings write (YYYY,'(i4.4)') date_time%year write (MM, '(i2.2)') date_time%month ! find out which data are needed ! compressed space dimension (always read global vector) start(1) = 1 icount(1) = N_gldas_compressed ! time dimension (first entry in GLDAS_NetCDF file is at 0Z) if ( (date_time%min/=0) .or. (date_time%sec/=0) .or. & (mod(date_time%hour,dt_gldas_in_hours)/=0) ) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'timing ERROR!!') end if hours_in_month = (date_time%day-1)*24 + date_time%hour start(2) = hours_in_month / dt_gldas_in_hours + 1 icount(2) = 1 ! ---------------------------------------------- ! ! compute indices for nearest neighbor interpolation from GLDAS grid ! to tile space ! ! (NOTE: this should at some point be replaced with a regridding ! subroutine that interpolates from the ! native forcing grid to the GCM atmospheric grid that is used ! to cut catchments into tiles - then "standard" grid2tile ! using tile_coord%atm_i and tile_coord%atm_j applies. ! reichle, 26 May 2005) do k=1,N_catd ! ll_lon and ll_lat refer to lower left corner of grid cell ! (as opposed to the grid point in the center of the grid cell) this_lon = tile_coord(k)%com_lon this_lat = tile_coord(k)%com_lat i_ind(k) = ceiling( (this_lon - gldas_grid_ll_lon)/gldas_grid_dlon ) j_ind(k) = ceiling( (this_lat - gldas_grid_ll_lat)/gldas_grid_dlat ) ! NOTE: For a "date line on center" grid and (180-dlon/2) < lon < 180 ! we now have i_ind=(grid%N_lon+1) ! This needs to be fixed. if (i_ind(k)>gldas_grid_N_lon) i_ind(k)=1 end do ! ------------------------------------------------------ ! ! read compression parameters (same for all data variables and time steps) gldas_var = 1 fname = trim(met_path) // trim(gldas_name(gldas_var)) // '/' // YYYY & // '/' // trim(gldas_name(gldas_var)) // '.' // YYYY // MM // '.nc' if (logit) write (logunit,*) 'get netcdf compression params from ' // trim(fname) ierr = NF_OPEN(fname,NF_NOWRITE,ncid) if (ierr/=0) then err_msg = 'error opening netcdf file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if ierr = NF_GET_VARA_INT(ncid, nciv_land_i, start, icount, land_i_gldas) ierr = NF_GET_VARA_INT(ncid, nciv_land_j, start, icount, land_j_gldas) ierr = NF_CLOSE(ncid) ! ------------------------------------------------------ ! ! get forcing data do gldas_var = 1,N_gldas_vars ! open file, read compressed data, and put on global grid fname = trim(met_path) // trim(gldas_name(gldas_var)) // '/' // YYYY & // '/' // trim(gldas_name(gldas_var)) // '.' // YYYY // MM // & '.nc' if (logit) write (logunit,*) 'opening ' // trim(fname) ierr = NF_OPEN(fname,NF_NOWRITE,ncid) if (ierr/=0) then err_msg = 'error opening netcdf file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if ierr = NF_GET_VARA_REAL(ncid, nciv_data, start, icount, tmp_vec ) ierr = NF_CLOSE(ncid) tmp_grid = nodata_forcing do n=1,N_gldas_compressed tmp_grid(land_i_gldas(n), land_j_gldas(n) ) = tmp_vec(n) end do ! interpolate to tile space ! (NOTE: This should at some point be replaced with a regridding ! subroutine that interpolates from the ! native forcing grid to the GCM atmospheric grid that is used ! to cut catchments into tiles - then "standard" grid2tile ! using tile_coord%atm_i and tile_coord%atm_j applies. ! reichle, 26 May 2005) do k=1,N_catd force_array(k,gldas_var) = tmp_grid(i_ind(k), j_ind(k)) end do end do ! convert variables and units of force_array to match met_force_type, ! put into structure ! from GLDAS files: ! ! force_array(:, 1) = SWdown W/m2 ! force_array(:, 2) = LWdown W/m2 ! force_array(:, 3) = Snowf kg/m2 (3h total) ! force_array(:, 4) = Rainf kg/m2 (3h total) ! force_array(:, 5) = Tair K ! force_array(:, 6) = Qair kg/kg ! force_array(:, 7) = Wind_E m/s ! force_array(:, 8) = Wind_N m/s ! force_array(:, 9) = PSurf Pa ! force_array(:,10) = RainfSnowf_C kg/m2 (3h total) ??? met_force_new%SWdown = force_array(:,1) met_force_new%LWdown = force_array(:,2) met_force_new%Tair = force_array(:,5) met_force_new%Qair = force_array(:,6) met_force_new%Psurf = force_array(:,9) ! get wind speed from wind vector do k=1,N_catd if ( abs(force_array(k,7)-nodata_gldas)GEOSgcm_grid_N_lon) i1(k)=1 ! for lat/lon shift ! ----------------- this_lon = tile_coord(k)%com_lon this_lat = tile_coord(k)%com_lat this_lon = this_lon + LON_SHIFT this_lat = this_lat + LAT_SHIFT i2(k) = ceiling((this_lon - GEOSgcm_grid_ll_lon)/GEOSgcm_grid_dlon) j2(k) = ceiling((this_lat - GEOSgcm_grid_ll_lat)/GEOSgcm_grid_dlat) ! NOTE: For a "date line on center" grid and (180-dlon/2) < lon < 180 ! we now have i1=(grid%N_lon+1) ! This needs to be fixed as follows: if (i2(k)>GEOSgcm_grid_N_lon) i2(k)=1 end do case (1) ! bilinear interpolation ! compute indices of nearest neighbors needed for bilinear ! interpolation from GEOSgcm grid to tile space do k=1,N_catd ! ll_lon and ll_lat refer to lower left corner of grid cell ! (as opposed to the grid point in the center of the grid cell) ! pchakrab: For bilinear interpolation, for each tile, we need: ! x1, x2, y1, y2 (defining the co-ords of four neighbors) and ! i1, i2, j1, j2 (defining the indices of four neighbors) ! find nearest neighbor forcing grid cell ("1") ! com of kth tile this_lon = tile_coord(k)%com_lon this_lat = tile_coord(k)%com_lat icur = ceiling((this_lon - GEOSgcm_grid_ll_lon)/GEOSgcm_grid_dlon) jcur = ceiling((this_lat - GEOSgcm_grid_ll_lat)/GEOSgcm_grid_dlat) ! wrap-around if (icur>GEOSgcm_grid_N_lon) icur = 1 if (jcur>GEOSgcm_grid_N_lat) then err_msg = "encountered tile near the poles" call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if xcur = real(icur-1)*GEOSgcm_grid_dlon - 180.0 ycur = real(jcur-1)*GEOSgcm_grid_dlat - 90.0 i1(k) = icur j1(k) = jcur x1(k) = xcur ! lon of grid cell center y1(k) = ycur ! lat of grid cell center ! find forcing grid cell ("2") diagonally across from icur, jcur tmp_lon = this_lon + 0.5*GEOSgcm_grid_dlon tmp_lat = this_lat + 0.5*GEOSgcm_grid_dlat inew = ceiling((tmp_lon - GEOSgcm_grid_ll_lon)/GEOSgcm_grid_dlon) jnew = ceiling((tmp_lat - GEOSgcm_grid_ll_lat)/GEOSgcm_grid_dlat) if (inew==icur) inew = inew - 1 if (jnew==jcur) jnew = jnew - 1 xnew = real(inew-1)*GEOSgcm_grid_dlon - 180.0 ynew = real(jnew-1)*GEOSgcm_grid_dlat - 90.0 ! wrap-around if (inew==0) inew = GEOSgcm_grid_N_lon if (inew>GEOSgcm_grid_N_lon) inew = 1 if ((jnew==0) .or. (jnew>GEOSgcm_grid_N_lat)) then err_msg = "encountered tile near the poles" call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if i2(k) = inew j2(k) = jnew x2(k) = xnew ! lon of grid cell center y2(k) = ynew ! lat of grid cell center end do case default err_msg = "unknown horizontal interpolation method" call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end select end if ! ---------------------------------------------- ! ! read global gridded field of given variable if ( (use_prec_corr) .and. (GEOSgcm_defs(GEOSgcm_var,1)(1:4)=='PREC') ) then ierr = NF_INQ_VARID( fid, GEOSgcm_defs(GEOSgcm_var,1), nv_id) ierr = NF_GET_VARA_REAL( fid, nv_id, istart, icount, tmp_grid) ierr = NF_CLOSE(fid) else call GFIO_GetVar( fid, trim(GEOSgcm_defs(GEOSgcm_var,1)), & YYYYMMDD, HHMMSS, GEOSgcm_grid_N_lon, GEOSgcm_grid_N_lat, & 0, 1, tmp_grid, rc ) if (rc<0) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error reading gfio file') end if call GFIO_Close ( fid, rc ) end if ! interpolate to tile space select case (met_hinterp) case (0) ! nearest-neighbor interpolation do k=1,N_catd if(GEOSgcm_defs(GEOSgcm_var,1)(1:4)=='PREC') then force_array(k,GEOSgcm_var) = tmp_grid(i2(k), j2(k)) else force_array(k,GEOSgcm_var) = tmp_grid(i1(k), j1(k)) endif end do case (1) ! bilinear interpolation do k=1,N_catd this_lon = tile_coord(k)%com_lon this_lat = tile_coord(k)%com_lat fnbr(1,1) = tmp_grid(i1(k),j1(k)) fnbr(1,2) = tmp_grid(i1(k),j2(k)) fnbr(2,1) = tmp_grid(i2(k),j1(k)) fnbr(2,2) = tmp_grid(i2(k),j2(k)) !DEC$ FORCEINLINE force_array(k,GEOSgcm_var) = BilinearInterpolation(this_lon, this_lat, & x1(k), x2(k), y1(k), y2(k), fnbr, nodata_forcing, tol) end do end select ! ---------------------------------------------- ! ! For non-flux "S" forcing variables that are read from MERRA "tavg" files ! *optionally* read forward-looking 'tavg' file (if available) and interpolate ! in time. ! ! Doing so minimizes the time shift between the "true" (but unavailable) MERRA ! instantaneous forcing values and their off-line time interpolated equivalent ! -- at the expense of a dampened diurnal cycle. ! ! reichle+qliu, 8 Oct 2008 minimize_shift = .true. ! minimize_shift should only affect "MERRA" forcing - reichle, 27 Feb 2012 if ((minimize_shift) .and. & (trim(GEOSgcm_defs(GEOSgcm_var,2))=='tavg') .and. & (trim(GEOSgcm_defs(GEOSgcm_var,5))=='S') ) then date_time_tmp = date_time_tavg_fwd ! open file call get_GEOS_gfio_openfile( date_time_tmp, daily_met_files, & met_path_tmp, met_tag_tmp, & GEOSgcm_defs(GEOSgcm_var,:), met_file_ext, fid) if (fid>0) then YYYYMMDD = date_time_tmp%year*10000+date_time_tmp%month*100+date_time_tmp%day HHMMSS = date_time_tmp%hour*10000+date_time_tmp%min*100 +date_time_tmp%sec ! read global gridded field of given variable call GFIO_GetVar( fid, trim(GEOSgcm_defs(GEOSgcm_var,1)), & YYYYMMDD, HHMMSS, GEOSgcm_grid_N_lon, GEOSgcm_grid_N_lat, & 0, 1, tmp_grid, rc ) if (rc<0) then err_msg = 'error reading gfio file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if call GFIO_Close ( fid, rc ) ! interpolate to tile space and in time select case (met_hinterp) case (0) ! nearest-neighbor interpolation do k=1,N_catd if ( abs(force_array(k,GEOSgcm_var) -nodata_GEOSgcm)0) end if ! if (minimize_shift) .and. [...] end do deallocate(tmp_grid) deallocate(GEOSgcm_defs) ! -------------------------------------------------------------------- ! convert variables and units of force_array to match met_force_type, ! put into structure ! from GEOSgcm files: ! ! G5DAS ! M2INT MERRA ! M2COR ! ! force_array(:, 1) = SWGDN SWGDN W/m2 (downward shortwave) ! force_array(:, 2) = SWLAND SWLAND W/m2 (net shortwave) ! force_array(:, 3) = LWGAB LWGAB W/m2 ("absorbed" longwave) ! force_array(:, 4) = PARDR PARDR W/m2 (direct PAR) ! force_array(:, 5) = PARDF PARDF W/m2 (diffuse PAR) ! force_array(:, 6) = PRECCU[*] PRECTOT kg/m2/s (*see below*) ! force_array(:, 7) = PRECLS[*] PRECCON kg/m2/s (*see below*) ! force_array(:, 8) = PRECSN[*] PRECSNO kg/m2/s (*see below*) ! force_array(:, 9) = PS PS Pa (surface air pressure) ! force_array(:,10) = HLML HLML m (height of lowest model level "LML") ! force_array(:,11) = TLML TLML K (air temperature at LML) ! force_array(:,12) = QLML QLML kg/kg (air spec humidity at LML) ! force_array(:,13) = SPEEDLML ULML m/s (wind speed/U-wind at LML) ! force_array(:,14) = n/a VLML m/s ( V-wind at LML) ! ! PRECTOT kg/m2/s (total rain+snow) = PRECCU+PRECLS+PRECSNO ! PRECCON kg/m2/s (convective rain+snow) ! PRECCU kg/m2/s (convective rain) ! PRECLS kg/m2/s (large-scale rain) ! PRECSNO kg/m2/s (total snow) met_force_new%SWdown = force_array(:, 1) met_force_new%SWnet = force_array(:, 2) met_force_new%LWdown = force_array(:, 3) met_force_new%PARdrct = force_array(:, 4) met_force_new%PARdffs = force_array(:, 5) met_force_new%Psurf = force_array(:, 9) met_force_new%RefH = force_array(:,10) met_force_new%Tair = force_array(:,11) met_force_new%Qair = force_array(:,12) do k=1,N_catd ! get wind speed if (MERRA_file_specs) then if ( abs(force_array(k,13)-nodata_GEOSgcm)0) then met_force_new(k)%Snowf = force_array(k,8) ! total_rain = total_precip - total_snow met_force_new(k)%Rainf = force_array(k,6) - force_array(k,8) ! conv_rain = (conv_precip/total_precip) * total_rain met_force_new(k)%Rainf_C = & force_array(k,7)/force_array(k,6)*met_force_new(k)%Rainf else met_force_new(k)%Rainf = 0. met_force_new(k)%Rainf_C = 0. met_force_new(k)%Snowf = 0. end if else ! G5DAS file specs met_force_new(k)%Rainf = force_array(k,6)+force_array(k,7) met_force_new(k)%Rainf_C = force_array(k,6) met_force_new(k)%Snowf = force_array(k,8) end if end if end do deallocate(force_array) end subroutine get_GEOS ! **************************************************************** function BilinearInterpolation(x,y,x1,x2,y1,y2,fnbr,UNDEF,tol) result(fxy) ! pchakrab: function to compute (bilinear) interpolated ! value f(x,y). If we know the value at 4 points ! f11=f(x1,y1), f12=f(x1,y2), f21=f(x2,y1) and f22=f(x2,y2), ! the interpolated value, fxy, is computed as ! ! f11(x2-x)(y2-y) + f21(x-x1)(y2-y) + f12(x2-x)(y-y1) + f22(x-x1)(y-y1) ! --------------------------------------------------------------------- ! (x2-x1)(y2-y1) ! ! NOTE 1: UNDEF is the nodata value (1e15) ! NOTE 2: If all the neighbouring f values are UNDEF, fxy = UNDEF ! else, the UNDEF value at a corner is replaced by an ! average of non-UNDEF values before fxy is computed real, intent(in) :: x, y, x1, x2, y1, y2, fnbr(2,2) real, intent(in) :: UNDEF, tol real :: fxy ! output ! local real :: floc(2,2), fsum, dx1, dx2, dy1, dy2 logical :: NoData(2,2) integer :: i, j, numNoData floc = fnbr ! check for nodata values ! and compute the sum of defined f values NoData = .false. numNodata = 0 fsum = 0.0 do j=1,2 do i=1,2 if (abs(floc(i,j)-UNDEF)0) then where (NoData) floc = fsum/real(4-numNoData) end if dx1 = x-x1 dx2 = x2-x dy1 = y-y1 dy2 = y2-y fxy = ( & floc(1,1)*dx2*dy2 + floc(2,1)*dx1*dy2 + & floc(1,2)*dx2*dy1 + floc(2,2)*dx1*dy1 & ) / ((x2-x1)*(y2-y1)) end if end function BilinearInterpolation ! **************************************************************** subroutine parse_MERRA_met_tag( met_path_in, met_tag_in, date_time, & met_path_default, met_path_prec, met_tag_out, use_prec_corr ) ! reichle, 1 Dec 2009 ! parse MERRA "met_tag", extract MERRA stream, assemble data paths ! ! Convention for driver_inputs*nml file or command line arguments: ! ! met_path = "/*/merra_land/" ! ! eg, on discover: /gpfsm/dnb51/projects/p15/iau/merra_land/ ! on discover: /discover/nobackup/qliu/merra_land/ ! on land01: /merra_land/ ! ! ! met_tag = "d5_merra_[STREAM]__[GCM-TAG]{__prec[PREC]}" ! ! where {__prec[PREC]} is optional and where ! ! STREAM = 'jan79', 'jan89', 'jan98', or 'cross' ! GCM-TAG = 'GEOSdas-2_1_4', ... ! PREC = 'CMAPvS', 'GPCPv1.1', ... ! ! examples: ! ! STREAM = 'jan89' : use only Stream 2 MERRA data ! STREAM = 'cross' : integrate across more than one stream ! GCM-TAG = 'GEOSdas-2_1_4' : tag of standard MERRA (may differ for replays) ! PREC : identifier for precip corrections to MERRA forcing ! ! --------------------------------------------------------------------------- implicit none character(200), intent(in) :: met_path_in character( 80), intent(in) :: met_tag_in type(date_time_type), intent(in) :: date_time character(200), intent(out) :: met_path_default, met_path_prec character( 80), intent(out) :: met_tag_out logical, intent(out) :: use_prec_corr ! local variables integer :: is, ie type(date_time_type) :: dt1, dt2 character( 5) :: stream character(80) :: gcm_tag, prec_tag character(len=*), parameter :: Iam = 'parse_MERRA_met_tag' character(len=400) :: err_msg ! ---------------------------------------------------------- ! define intervals that determine which MERRA stream is used ! in integrations that "cross" multiple streams ! Stream 1 ('jan79') --> 16 Dec 1978 - 31 Dec 1989 ! Stream 2 ('jan89') --> 1 Jan 1990 - 31 Dec 1998 ! Stream 3 ('jan98') --> 1 Jan 1999 - present ! dates before dt1 use Stream 1 dt1%year = 1993 dt1%month = 1 dt1%day = 1 dt1%hour = 0 dt1%min = 0 dt1%sec = 0 ! otherwise, dates before dt2 use Stream 2 dt2%year = 2001 dt2%month = 1 dt2%day = 1 dt2%hour = 0 dt2%min = 0 dt2%sec = 0 ! ---------------------------------------------------- ! initialize met_tag_out = repeat(' ', len(met_tag_out)) ! define which stream to use if (met_tag_in(10:14)=='cross') then if (datetime_lt_refdatetime( date_time, dt1 )) then stream = 'jan79' elseif (datetime_lt_refdatetime( date_time, dt2 )) then stream = 'jan89' else stream = 'jan98' end if met_tag_out = met_tag_in(1:9) // stream else met_tag_out = met_tag_in(1:14) end if ! ----------------------------------------------------- ! ! identify GCM tag and which precip corrections to use, ! assemble met_path accordingly ! ! met_tag = "d5_merra_[STREAM]__[GCM-TAG]{__prec[PREC]}" ! ! where {__prec[PREC]} is optional is = index( met_tag_in, '__') ie = index( met_tag_in, '__', .true.) if (is/=ie) then ! using precip corrections gcm_tag = met_tag_in(is+2:ie-1) met_path_default = trim(met_path_in) // '/MERRA_land_forcing/' // & trim(gcm_tag) // '/' prec_tag = met_tag_in(ie+6:len(met_tag_in)) met_path_prec = trim(met_path_in) // '/precip_corr_' // trim(prec_tag) // & '/' // trim(gcm_tag) // '/' use_prec_corr = .true. else ! not using precip corrections gcm_tag = met_tag_in(is+2:len(met_tag_in)) met_path_default = trim(met_path_in) // '/MERRA_land_forcing/' // & trim(gcm_tag) // '/' prec_tag = repeat(' ', len(prec_tag)) met_path_prec = met_path_default use_prec_corr = .false. ! check if prec_tag was accidentally appended with a single underscore if (len_trim(met_tag_in)>29) then err_msg = 'questionable met_tag_in, not enough double underscores' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if end if ! use "ignore_SWNET_for_snow" to fix snow albedo bug in MERRA forcing: ! ! For "MERRA-Land" do *not* use MERRA SWNET for snow-covered surface ! because of the bug in the call to sibalb() in MERRA. This is communicated ! to subroutine propagate_cat() via "ignore_SWNET_for_snow=.true." ! In all other cases, "ignore_SWNET_for_snow=.false." ! !!if (use_prec_corr) then !! !! ! MERRA-Land !! !! ignore_SWNET_for_snow = .true. !! !! if (logit) write (logunit,*) 'ignore_SWNET_for_snow = ', ignore_SWNET_for_snow !! !!else !! !! ! MERRA replay !! !! ignore_SWNET_for_snow = .false. !! !!end if ! ! The above fix did not work in MPI because subroutine get_forcing() is ! only called by the master process. All other processes are unaware of ! any changes to "ignore_SWnet_for_snow" from its uninitialized value ! because an MPI broadcast was missing. ! As of April 2015, "ignore_SWnet_for_snow" is no longer meaningful. ! - reichle, 2 Apr 2015 end subroutine parse_MERRA_met_tag ! **************************************************************** ! **************************************************************** subroutine parse_MERRA2_met_tag( met_path_in, met_tag_in, date_time, & met_path_default, met_path_prec, met_tag_out, use_prec_corr ) ! reichle, 27 Jul 2015 ! parse MERRA2 "met_tag", extract MERRA stream, assemble data paths ! ! met_tag = "M2[xxx]_[STREAM]{__prec[PREC]}" ! ! where {__prec[PREC]} is optional and where ! ! [xxx] = 'GCM' or 'COR' ! STREAM = '100', '200', '300', '400', or 'cross' ! PREC = 'CMAPvS', 'GPCPv1.1', ... ! ! examples: ! ! STREAM = '200' : use only Stream 2 MERRA data ! STREAM = 'cross' : integrate across more than one stream ! PREC : identifier for corrected precip data ! ! --------------------------------------------------------------------------- implicit none character(200), intent(in) :: met_path_in character( 80), intent(in) :: met_tag_in type(date_time_type), intent(in) :: date_time character(200), intent(out) :: met_path_default, met_path_prec character( 80), intent(out) :: met_tag_out logical, intent(out) :: use_prec_corr ! local variables integer :: is type(date_time_type) :: dt1, dt2, dt3 character(10) :: stream character(80) :: prec_tag character(len=*), parameter :: Iam = 'parse_MERRA2_met_tag' character(len=400) :: err_msg ! ---------------------------------------------------------- ! define intervals that determine which MERRA-2 stream is used ! in integrations that "cross" multiple streams ! ! 1/1/1980 - 12/31/1991: MERRA2_100 (Stream 1) ! 1/1/1992 - 12/31/2000: MERRA2_200 (Stream 2) ! 1/1/2001 - 12/31/2010: MERRA2_300 (Stream 3) ! 1/1/2011 - present: MERRA2_400 (Stream 4) ! dates before dt1 use Stream 1 dt1%year = 1992 dt1%month = 1 dt1%day = 1 dt1%hour = 0 dt1%min = 0 dt1%sec = 0 ! otherwise, dates before dt2 use Stream 2 dt2%year = 2001 dt2%month = 1 dt2%day = 1 dt2%hour = 0 dt2%min = 0 dt2%sec = 0 ! otherwise, dates before dt3 use Stream 3 dt3%year = 2011 dt3%month = 1 dt3%day = 1 dt3%hour = 0 dt3%min = 0 dt3%sec = 0 ! ---------------------------------------------------- ! initialize met_tag_out = repeat(' ', len(met_tag_out)) stream = repeat(' ', len(stream )) ! define which stream to use if (met_tag_in(7:11)=='cross') then if (datetime_lt_refdatetime( date_time, dt1 )) then stream = 'MERRA2_100' elseif (datetime_lt_refdatetime( date_time, dt2 )) then stream = 'MERRA2_200' elseif (datetime_lt_refdatetime( date_time, dt3 )) then stream = 'MERRA2_300' else stream = 'MERRA2_400' end if met_tag_out = trim(stream) else met_tag_out = 'MERRA2_' // met_tag_in(7:9) end if met_path_default = trim(met_path_in) // '/' ! ----------------------------------------------------- ! ! identify which precip corrections to use, ! assemble met_path accordingly ! ! met_tag = "M2[xxx]_[STREAM]{__prec[PREC]}" ! ! where {__prec[PREC]} is optional is = index( met_tag_in, '__') if (is>0) then ! using precip corrections prec_tag = met_tag_in(is+6:len(met_tag_in)) met_path_prec = trim(met_path_in) // '/precip_corr_' // trim(prec_tag) // '/' use_prec_corr = .true. else ! not using precip corrections prec_tag = repeat(' ', len(prec_tag)) met_path_prec = met_path_default use_prec_corr = .false. ! check if prec_tag was accidentally appended with a single underscore if (len_trim(met_tag_in)>11) then err_msg = 'questionable met_tag_in, not enough double underscores' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if end if end subroutine parse_MERRA2_met_tag ! **************************************************************** subroutine parse_G5DAS_met_tag( met_path_in, met_tag_in, date_time, & met_path_default, met_path_prec, met_tag_out, use_prec_corr, & use_Predictor ) ! parse G5DAS "met_tag" ! ! Convention for driver_inputs*nml file or command line arguments: ! ! met_tag = "[G5DAS-NAME]{__prec[PREC]}{__Nx+-}" ! ! where {__prec[PREC]} and {__Nx+-} are optional and where ! ! G5DAS-NAME : name of standard G5DAS forcing (must not contain "__"!) ! e.g., "e5110_fp", "d591_rpit1", "d591_fpit", ... ! for cross-stream forcing, use "cross_d5124_RPFPIT" or "cross_FP" ! PREC : identifier for precip corrections to G5DAS forcing (eg., 'GPCPv1.1') ! ! If {__Nx+-} is present, set flag for use forcing files from the DAS/GCM Predictor ! segment. (Default is to use forcing from Corrector segment.) ! ! reichle, 3 Jun 2013 ! reichle, 20 Sep 2013: restructured, added "use_Predictor" ! reichle, 23 Sep 2013: added "cross-stream" capability ! reichle, 5 Dec 2014: updated "cross-stream" dates ! reichle, 14 Sep 2015: revised FP "cross-stream" dates ! reichle, 28 Dec 2016: added 5.12.4 RPIT/FPIT streams ! reichle, 19 Jan 2017: added FP transition from e5131 to f516 ! ! --------------------------------------------------------------------------- implicit none character(200), intent(in) :: met_path_in character( 80), intent(in) :: met_tag_in type(date_time_type), intent(in) :: date_time character(200), intent(out) :: met_path_default, met_path_prec character( 80), intent(out) :: met_tag_out logical, intent(out) :: use_prec_corr logical, intent(out) :: use_Predictor ! local variables integer :: is, ie, ii, N_opt_tag character(80) :: prec_tag, stream character(80), dimension(2) :: tmp_tag type(date_time_type) :: dt_end_d591_rpit1 type(date_time_type) :: dt_end_d591_rpit2 type(date_time_type) :: dt_end_d591_rpit3 type(date_time_type) :: dt_end_d591_fpit type(date_time_type) :: dt_end_d5124_rpit1 type(date_time_type) :: dt_end_d5124_rpit2 type(date_time_type) :: dt_end_d5124_rpit3 type(date_time_type) :: dt_end_e5110_fp type(date_time_type) :: dt_end_e5130_fp type(date_time_type) :: dt_end_e5131_fp type(date_time_type) :: dt_end_f516_fp character(len=*), parameter :: Iam = 'parse_G5DAS_met_tag' character(len=400) :: err_msg ! ---------------------------------------------------------- ! ! define transition times between RP-IT, FP-IT or FP streams ! if "cross-stream" forcing is requested ! | stream start | stream end (as of 5 Dec 2014) ! ---------------------------------------- ! d591_rpit1 | 1 Jan 2000 | 1 Jun 2007 ! d591_rpit2 | 1 Jun 2006 | 30 Dec 2011 ! d591_rpit3 | 1 Jan 2011 | 6 May 2013 ! d591_fpit | 30 Dec 2012 | (present) dt_end_d591_rpit1%year = 2007 dt_end_d591_rpit1%month = 6 dt_end_d591_rpit1%day = 1 dt_end_d591_rpit1%hour = 0 dt_end_d591_rpit1%min = 0 dt_end_d591_rpit1%sec = 0 dt_end_d591_rpit2%year = 2011 dt_end_d591_rpit2%month = 12 dt_end_d591_rpit2%day = 1 dt_end_d591_rpit2%hour = 0 dt_end_d591_rpit2%min = 0 dt_end_d591_rpit2%sec = 0 dt_end_d591_rpit3%year = 2013 dt_end_d591_rpit3%month = 5 dt_end_d591_rpit3%day = 1 dt_end_d591_rpit3%hour = 0 dt_end_d591_rpit3%min = 0 dt_end_d591_rpit3%sec = 0 dt_end_d591_fpit%year = 9999 dt_end_d591_fpit%month = 1 dt_end_d591_fpit%day = 1 dt_end_d591_fpit%hour = 0 dt_end_d591_fpit%min = 0 dt_end_d591_fpit%sec = 0 ! | stream start | stream end (as of 28 Dec 2016) ! ---------------------------------------- ! d5124_rpit1 | 1 Jan 2000 | 1 Jan 2004 ! d5124_rpit2 | 1 Jan 2004 | 1 Jan 2012 ! d5124_rpit3 | 1 Jan 2012 | (present) dt_end_d5124_rpit1%year = 2004 dt_end_d5124_rpit1%month = 1 dt_end_d5124_rpit1%day = 1 dt_end_d5124_rpit1%hour = 0 dt_end_d5124_rpit1%min = 0 dt_end_d5124_rpit1%sec = 0 dt_end_d5124_rpit2%year = 2012 dt_end_d5124_rpit2%month = 1 dt_end_d5124_rpit2%day = 1 dt_end_d5124_rpit2%hour = 0 dt_end_d5124_rpit2%min = 0 dt_end_d5124_rpit2%sec = 0 dt_end_d5124_rpit3%year = 9999 dt_end_d5124_rpit3%month = 1 dt_end_d5124_rpit3%day = 1 dt_end_d5124_rpit3%hour = 0 dt_end_d5124_rpit3%min = 0 dt_end_d5124_rpit3%sec = 0 ! --------------------------------- ! ! FP streams ! ! Stream start/end in terms of availability ! in the GEOS-5 archive (approximately): ! ! | stream start | stream end ! ---------------------------------------- ! e5110_fp | 11 Jun 2013 | 19 Aug 2014 ! e5130_fp | 1 Aug 2014 | 4 May 2015 ! e5131_fp | 1 May 2015 | 24 Jan 2017 ! e5131_fp | 24 Jan 2015 | (present) ! ! Official stream transition times (as defined ! by GMAO ops group) are: ! ! FP e5110 --> e5130 : 20 Aug 2014, 6z ADAS analysis ! FP e5130 --> e5131 : 1 May 2015, 6z ADAS analysis ! FP e5131 --> f516 : 24 Jan 2017, 6z ADAS analysis ! ! Official stream transition times refer to the definition ! of the official FP files with generic file names on the ! NCCS data portal. ! ! Note that "lfo" files for the 6z analysis start from 4z, ! that is, the exact transition time is slightly different ! from the LDAS perspective. ! ! Define LDASsa "cross" streams dates for 0z of the ! first day that contains only "new" data for ! compatibility with the SMAP L4 Ops system and because of the ! asymmetric availability of the data on the NCCS data portal. ! ! - reichle, 14 Sep 2015 ! ! Revised Aug 2014 cross-over date to Aug 20 at 0z because ! e5110 data are not available for Aug 20 ! - reichle+qliu, 29 Jan 2016 dt_end_e5110_fp%year = 2014 dt_end_e5110_fp%month = 8 dt_end_e5110_fp%day = 20 dt_end_e5110_fp%hour = 0 dt_end_e5110_fp%min = 0 dt_end_e5110_fp%sec = 0 dt_end_e5130_fp%year = 2015 dt_end_e5130_fp%month = 5 dt_end_e5130_fp%day = 2 dt_end_e5130_fp%hour = 0 dt_end_e5130_fp%min = 0 dt_end_e5130_fp%sec = 0 dt_end_e5131_fp%year = 2017 dt_end_e5131_fp%month = 1 dt_end_e5131_fp%day = 24 dt_end_e5131_fp%hour = 3 dt_end_e5131_fp%min = 0 dt_end_e5131_fp%sec = 0 dt_end_f516_fp%year = 9999 dt_end_f516_fp%month = 1 dt_end_f516_fp%day = 1 dt_end_f516_fp%hour = 0 dt_end_f516_fp%min = 0 dt_end_f516_fp%sec = 0 ! ---------------------------------------------------- ! initialize met_tag_out = repeat(' ', len(met_tag_out)) stream = repeat(' ', len(stream )) use_prec_corr = .false. use_Predictor = .false. ! ----------------------------------------------------- ! ! identify "GCM" tag and whether to use precip corrections and/or ! forcing from the DAS/GCM Predictor segment (default: Corrector segment) ! ! met_tag = "[G5DAS-NAME]{__prec[PREC]}{__Nx+-}" ! ! where {__prec[PREC]} and {__Nx+-} are optional is = index( met_tag_in, '__') ie = index( met_tag_in, '__', .true.) ! determine how many optional tag segments are present if (is==0) then N_opt_tag = 0 met_tag_out = met_tag_in elseif (is==ie) then N_opt_tag = 1 met_tag_out = met_tag_in(1:(is-1)) tmp_tag(1) = met_tag_in((is+2):len(met_tag_in)) else ! make sure there are at most two optional tag segments ! "ii" should be the index for the next "__" after "is" ii = index( met_tag_in(is+2:len(met_tag_in)), '__') if (is+1+ii==ie) then N_opt_tag = 2 met_tag_out = met_tag_in(1:(is-1)) tmp_tag(1) = met_tag_in((is+2):ie-1) tmp_tag(2) = met_tag_in((ie+2):len(met_tag_in)) else err_msg = 'invalid met_tag_in, too many double underscores' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if end if ! resolve cross-stream if requested if (met_tag_out(1:17)=='cross_d591_RPFPIT') then if (datetime_lt_refdatetime( date_time, dt_end_d591_rpit1 )) then stream = 'd591_rpit1_jan00' ! use d591 RP-IT stream 1 elseif (datetime_lt_refdatetime( date_time, dt_end_d591_rpit2 )) then stream = 'd591_rpit2_jun06' ! use d591 RP-IT stream 2 elseif (datetime_lt_refdatetime( date_time, dt_end_d591_rpit3 )) then stream = 'd591_rpit3_jan11' ! use d591 RP-IT stream 3 else stream = 'd591_fpit' ! use d591 FP-IT end if elseif (met_tag_out(1:18)=='cross_d5124_RPFPIT') then if (datetime_lt_refdatetime( date_time, dt_end_d5124_rpit1 )) then stream = 'd5124_rpit_jan00' ! use d5124 RP-IT stream 1 elseif (datetime_lt_refdatetime( date_time, dt_end_d5124_rpit2 )) then stream = 'd5124_rpit_jan04' ! use d5124 RP-IT stream 2 else stream = 'd5124_rpit_jan12' ! use d5124 RP-IT stream 3 end if elseif (met_tag_out(1:8)=='cross_FP') then if (datetime_lt_refdatetime( date_time, dt_end_e5110_fp )) then stream = 'e5110_fp' ! use GEOS-5.11.0 output elseif (datetime_lt_refdatetime( date_time, dt_end_e5130_fp )) then stream = 'e5130_fp' ! use GEOS-5.13.0 output elseif (datetime_le_refdatetime( date_time, dt_end_e5131_fp )) then ! Note "less-than-or-equal" (_le_) above stream = 'e5131_fp' ! use GEOS-5.13.1 output else stream = 'f516_fp' ! use GEOS-5.16 output end if else stream = met_tag_out end if met_tag_out = stream ! interpret optional tag segments do ii=1,N_opt_tag if (tmp_tag(ii)(1:4)=='prec') then use_prec_corr = .true. prec_tag = tmp_tag(ii)(5:len(tmp_tag(ii))) elseif (tmp_tag(ii)(1:4)=='Nx+-') then use_Predictor = .true. else err_msg = 'invalid met_tag_in, unknown optional tag segment' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if end do ! get path to met forcing (except precip) met_path_default = trim(met_path_in) // '/' ! get path to precip if (use_prec_corr) then ! using precip corrections met_path_prec = trim(met_path_in) // '/precip_corr_' // trim(prec_tag) // '/' else ! *not* using precip corrections met_path_prec = met_path_default end if ! Double-check if optional tag segments were somehow missed, e.g., ! because they were accidentally appended with single underscores. ! Assumes that "prec" and "Nx+-" never appear in GEOS-5 product names. ! Does NOT protect against spelling errors, e.g., "__perc" or "__Nx-+". ! - reichle, 24 Nov 2015 if ((.not. use_prec_corr) .and. (index( met_tag_in, 'prec')>0)) then err_msg = 'questionable met_tag_in: includes "prec" but use_prec_corr=.false.' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if if ((.not. use_Predictor) .and. (index( met_tag_in, 'Nx+-')>0)) then err_msg = 'questionable met_tag_in: includes "Nx+-" but use_Predictor=.false.' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if end subroutine parse_G5DAS_met_tag ! **************************************************************** subroutine get_GEOS_gfio_openfile( date_time, daily_file, met_path, met_tag, & GEOSgcm_defs, file_ext, fid, inquire_only ) ! reichle, 27 Jul 2015 - added "daily_file" option ! (because MERRA-2 provides aggregated daily files that contain 24 hourly fields) implicit none type(date_time_type), intent(in) :: date_time logical, intent(in) :: daily_file character(200), intent(in) :: met_path character( 80), intent(in) :: met_tag character( 40), dimension(5), intent(in) :: GEOSgcm_defs character( 3), intent(in) :: file_ext integer, intent(out) :: fid logical, optional, intent(in) :: inquire_only ! local variables character(300) :: fname, fname_full, fname_full_tmp character( 14) :: time_stamp character( 4) :: YYYY, HHMM character( 2) :: MM, DD integer :: i, rc logical :: inquire_only_tmp, file_exists character(len=*), parameter :: Iam = 'get_GEOS_gfio_openfile' character :: err_msg ! ------------------------------------------------------------- ! ! initialize inquire_only_tmp = .false. if (present(inquire_only)) inquire_only_tmp = inquire_only fid = -9999 ! assemble date/time strings write (YYYY,'(i4.4)') date_time%year write (MM, '(i2.2)') date_time%month write (DD, '(i2.2)') date_time%day write (HHMM,'(i4.4)') date_time%hour*100+date_time%min ! deal with absence of "minutes" in "bkg.sfc" file name ! (replace minutes with blanks, then trim, see below) if (trim(GEOSgcm_defs(3))=='bkg.sfc') HHMM = HHMM(1:2) // ' ' ! assemble file name time_stamp = repeat(' ', len(time_stamp)) if (daily_file) then time_stamp(1:8) = YYYY // MM // DD else time_stamp = YYYY // MM // DD // '_' // trim(HHMM) // 'z' end if fname = trim(met_tag) // '.' // trim(GEOSgcm_defs(3)) // '.' // & trim(time_stamp) // '.' // file_ext ! ---------------------------------------------- ! ! open file ! ! Try reading the files directly inside directory "met_path/" first (because in ! coupled DAS mode met_path=workdir, and the files are simply sitting there). ! If this fails, try reading the files in "met_path/met_tag/*/Yyyyy/Mmm/" ! as in the archived directory structure. do i=1,2 if (i==1) then fname_full = trim(met_path) // '/' // trim(fname) fname_full_tmp = fname_full ! remember for error log below elseif (i==2) then fname_full = trim(met_path) // '/' // trim(met_tag) // '/' // & trim(GEOSgcm_defs(4)) // '/Y' // YYYY // '/M' // MM // '/' // trim(fname) else call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error') end if inquire(file=fname_full, exist=file_exists) if (file_exists) then fid = 0 if (.not. inquire_only_tmp) then call Gfio_Open ( fname_full, 1, fid, rc ) if (rc/=0) then fid = -9999 else if (logit) write (logunit,'(400A)') & 'read ' // trim(GEOSgcm_defs(1)) // ' from ' // trim(fname_full) end if end if exit ! if file exists exit loop through i end if end do ! create error log if all attempts to open a file failed if ((fid<0) .and. (.not. inquire_only_tmp)) then if (logit) then write (logunit,*) & 'get_GEOS_gfio_openfile(): Unsuccessfully tried to open files:' write (logunit,'(400A)') trim(fname_full_tmp) write (logunit,'(400A)') trim(fname_full) end if end if end subroutine get_GEOS_gfio_openfile ! **************************************************************** subroutine get_GEOS_netcdf_openfile( date_time, met_path, met_tag, & GEOSgcm_defs, file_ext, ncid, ierr ) implicit none type(date_time_type), intent(in) :: date_time character(200), intent(in) :: met_path character( 80), intent(in) :: met_tag character( 40), dimension(5), intent(in) :: GEOSgcm_defs character( 3), intent(in) :: file_ext integer, intent(out) :: ncid, ierr ! local variables character(100) :: fname character(200) :: fdir character(300) :: fname_full, fname_full_tmp character( 4) :: YYYY, HHMM character( 2) :: MM, DD ! ------------------------------------------------------------- ! ! assemble date/time strings write (YYYY,'(i4.4)') date_time%year write (MM, '(i2.2)') date_time%month write (DD, '(i2.2)') date_time%day write (HHMM,'(i4.4)') date_time%hour*100+date_time%min ! assemble file name fname = trim(met_tag) // '.' // trim(GEOSgcm_defs(3)) // & '.' // YYYY // MM // DD // '_' // trim(HHMM) // 'z.' // & trim(file_ext) ! assemble dir name without "/Mmm" (month) dir fdir = trim(met_path) // '/' // trim(met_tag) // '/' // & trim(GEOSgcm_defs(4)) // '/' // 'Y' // YYYY // '/' ! ----------------------------------------------------------------------- ! try opening file with "/Mmm" (month) dir ! (standard for corrected G5DAS precip) fname_full = trim(fdir) // 'M' // MM // '/' // trim(fname) ierr = NF_OPEN(fname_full, NF_NOWRITE, ncid) if (ierr/=0) then fname_full_tmp = fname_full ! remember for error log below ! try again *without* "/Mmm" (month) dir ! (for corrected MERRA precip) fname_full = trim(fdir) // trim(fname) ierr = NF_OPEN(fname_full, NF_NOWRITE, ncid) if (ierr/=0) then if (logit) then write (logunit,*) & 'get_GEOS_netcdf_openfile(): Unsuccessfully tried to open files:' write (logunit,'(400A)') trim(fname_full_tmp) write (logunit,'(400A)') trim(fname_full) end if else if (logit) write (logunit,'(400A)') & 'read ' // trim(GEOSgcm_defs(1)) // ' from ' // trim(fname_full) end if else if (logit) write (logunit,'(400A)') & 'read ' // trim(GEOSgcm_defs(1)) // ' from ' // trim(fname_full) end if end subroutine get_GEOS_netcdf_openfile ! **************************************************************** subroutine get_GSWP2_1x1_netcdf(date_time, met_path, N_catd, tile_coord, & met_force_new, nodata_forcing ) ! read GSWP2_NetCDF files and extract forcings in tile space ! (uses nearest neighbor interpolation) ! reichle, 28 Jul 2005 implicit none type(date_time_type), intent(in) :: date_time character(200), intent(in) :: met_path integer, intent(in) :: N_catd type(tile_coord_type), dimension(:), pointer :: tile_coord ! input type(met_force_type), dimension(N_catd), intent(inout) :: met_force_new real, intent(out) :: nodata_forcing ! GEOS5-DAS grid and netcdf parameters integer, parameter :: gswp2_grid_N_lon = 360 integer, parameter :: gswp2_grid_N_lat = 150 real, parameter :: gswp2_grid_ll_lon = -180. real, parameter :: gswp2_grid_ll_lat = -60. real, parameter :: gswp2_grid_dlon = 1. real, parameter :: gswp2_grid_dlat = 1. integer, parameter :: N_gswp2_compressed = 15238 ! GSWP2 forcing time step in hours integer, parameter :: dt_gswp2_in_hours = 3 integer, parameter :: nciv_land = 3 integer, parameter :: nciv_data = 6 integer, parameter :: N_gswp2_vars = 9 real, parameter :: nodata_gswp2 = 1.e20 character(40), dimension(N_gswp2_vars), parameter :: gswp2_name = (/ & 'SWdown_srb ', & ! 1 - flux 'LWdown_srb ', & ! 2 - flux 'Rainf_C_gswp ', & ! 3 - flux 'Rainf_gswp ', & ! 4 - flux 'Snowf_gswp ', & ! 5 - flux 'PSurf_ecor ', & ! 6 - state 'Qair_cru ', & ! 7 - state 'Tair_cru ', & ! 8 - state 'Wind_ncep ' /) ! 9 - state ! local variables real :: tol real, dimension(gswp2_grid_N_lon,gswp2_grid_N_lat) :: tmp_grid integer, dimension(N_gswp2_compressed) :: land_gswp2 integer, dimension(N_gswp2_compressed) :: land_i_gswp2, land_j_gswp2 integer, dimension(N_catd) :: i_ind, j_ind real, dimension(N_gswp2_compressed) :: tmp_vec real, dimension(N_catd,N_gswp2_vars) :: force_array integer, dimension(2) :: start, icount integer :: k, n, hours_in_month, gswp2_var, ierr, ncid integer :: int_dt_gswp2_in_seconds real :: this_lon, this_lat type(date_time_type) :: date_time_tmp character(4) :: YYYY character(2) :: MM character(300) :: fname character(len=*), parameter :: Iam = 'get_GSWP2_1x1_netcdf' character(len=400) :: err_msg ! -------------------------------------------------------------------- int_dt_gswp2_in_seconds = 3600*dt_gswp2_in_hours nodata_forcing = nodata_gswp2 tol = abs(nodata_forcing*nodata_tolfrac_generic) ! ---------------------------------------------- ! ! compute indices for nearest neighbor interpolation from GSWP2 grid ! to tile space ! ! (NOTE: this should at some point be replaced with a regridding ! subroutine that interpolates from the ! native forcing grid to the GCM atmospheric grid that is used ! to cut catchments into tiles - then "standard" grid2tile ! using tile_coord%atm_i and tile_coord%atm_j applies. ! reichle, 26 May 2005) do k=1,N_catd ! ll_lon and ll_lat refer to lower left corner of grid cell ! (as opposed to the grid point in the center of the grid cell) this_lon = tile_coord(k)%com_lon this_lat = tile_coord(k)%com_lat ! i_ind, j_ind count eastward and northward ! (note that lat/lon coordinates in GSWP2 netcdf files ! are eastward and southward) i_ind(k) = ceiling( (this_lon - gswp2_grid_ll_lon)/gswp2_grid_dlon ) j_ind(k) = ceiling( (this_lat - gswp2_grid_ll_lat)/gswp2_grid_dlat ) end do ! ------------------------------------------------------ ! ! space dimension is same for all variables start(1) = 1 icount(1) = N_gswp2_compressed ! check for possible error with time if ( (date_time%min/=0) .or. (date_time%sec/=0) .or. & (mod(date_time%hour,dt_gswp2_in_hours)/=0) ) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'timing ERROR!!') end if ! ------------------------------------------------------ ! ! get forcing data do gswp2_var = 1,N_gswp2_vars ! time dimension ! ! First entry in GSWP2_NetCDF file is at 3Z, with fluxes for 0Z-3Z ! ! At 0Z for first day of month: ! - for fluxes read first entry of that month ! - for states read last entry of preceding month ! At 3Z for first day of month: ! - for fluxes read second entry of that month ! - for states read first entry of that month ! and so on... select case (gswp2_var) case (1,2,3,4,5) ! "fluxes" date_time_tmp = date_time case (6,7,8,9) ! "states" date_time_tmp = date_time call augment_date_time( -int_dt_gswp2_in_seconds, date_time_tmp ) case default call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error') end select hours_in_month = (date_time_tmp%day-1)*24 + date_time_tmp%hour start(2) = hours_in_month / dt_gswp2_in_hours + 1 icount(2) = 1 ! assemble year and month strings write (YYYY,'(i4.4)') date_time_tmp%year write (MM, '(i2.2)') date_time_tmp%month ! assemble file name, open file fname = trim(met_path) // trim(gswp2_name(gswp2_var)) // '/' & // '/' // trim(gswp2_name(gswp2_var)) // YYYY // MM // '.nc' if (logit) write (logunit,*) 'opening ' // trim(fname) ierr = NF_OPEN(fname,NF_NOWRITE,ncid) if (ierr/=0) then err_msg = 'error opening netcdf file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if ! ------------------------------------------------------ ! ! read compression parameters (same for all data variables and time steps) if (gswp2_var == 1) then if (logit) write (logunit,*) 'get netcdf compression params from ' // trim(fname) if (ierr/=0) then err_msg = 'error opening netcdf file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if ierr = NF_GET_VARA_INT(ncid, nciv_land, start, icount, land_gswp2) ! land_i_gswp2, land_j_gswp2 count eastward and southward ! (note that lat/lon coordinates in GSWP2 netcdf files ! are eastward and southward) do k=1,N_gswp2_compressed land_j_gswp2(k) = (land_gswp2(k)-1)/gswp2_grid_N_lon + 1 land_i_gswp2(k) = & land_gswp2(k) - (land_j_gswp2(k)-1)*gswp2_grid_N_lon end do end if ! ------------------------------------------------------ ! ! read compressed data, and put on global grid ierr = NF_GET_VARA_REAL(ncid, nciv_data, start, icount, tmp_vec ) ierr = NF_CLOSE(ncid) tmp_grid = nodata_forcing do n=1,N_gswp2_compressed tmp_grid(land_i_gswp2(n), land_j_gswp2(n) ) = tmp_vec(n) end do ! flip tmp_grid ! (land_j_gswp2 counts southward, whereas j_ind counts northward) tmp_grid = tmp_grid(:,gswp2_grid_N_lat:1:-1) !!do k=1,gswp2_grid_N_lon !! do n=1,gswp2_grid_N_lat !! !! write (999,*) k, n, tmp_grid(k,n) !! !! end do !!end do !!write (logunit,*) 'debug stop here' !!stop ! interpolate to tile space ! (NOTE: This should at some point be replaced with a regridding ! subroutine that interpolates from the ! native forcing grid to the GCM atmospheric grid that is used ! to cut catchments into tiles - then "standard" grid2tile ! using tile_coord%atm_i and tile_coord%atm_j applies. ! reichle, 26 May 2005) do k=1,N_catd force_array(k,gswp2_var) = tmp_grid(i_ind(k), j_ind(k)) end do end do ! convert variables and units of force_array to match met_force_type, ! put into structure ! from GSWP2 files: ! ! force_array(:, 1) = SWdown_srb W/m2 ! force_array(:, 2) = LWdown_srb W/m2 ! force_array(:, 3) = Rainf_C_gswp kg/m2/s ! force_array(:, 4) = Rainf_gswp kg/m2/s ! force_array(:, 5) = Snowf_gswp kg/m2/s ! force_array(:, 6) = PSurf_ecor Pa ! force_array(:, 7) = Qair_cru kg/kg ! force_array(:, 8) = Tair_cru K ! force_array(:, 9) = Wind_ncep m/s met_force_new%SWdown = force_array(:,1) met_force_new%LWdown = force_array(:,2) met_force_new%Rainf_C = force_array(:,3) met_force_new%Rainf = force_array(:,4) met_force_new%Snowf = force_array(:,5) met_force_new%Psurf = force_array(:,6) met_force_new%Qair = force_array(:,7) met_force_new%Tair = force_array(:,8) met_force_new%Wind = force_array(:,9) end subroutine get_GSWP2_1x1_netcdf ! **************************************************************** subroutine check_forcing_nodata( N_catd, tile_coord, nodata_forcing, met_force ) ! check input forcing for no-data-values and unphysical values ! ! If no-data-value is encountered, use value from "next" catchment, ! where "next" is next in line (not necessarily next in distance!). ! if that does not work, give up ! ! reset unphysical values as best as possible ! ! reichle, 13 May 2003 ! reichle, 13 Jun 2005 implicit none integer, intent(in) :: N_catd type(tile_coord_type), dimension(:), pointer :: tile_coord ! input real, intent(in) :: nodata_forcing type(met_force_type), dimension(N_catd), intent(inout) :: met_force ! local variables integer :: i real :: tol ! ------------------------------------------------------------ call check_forcing_nodata_2(N_catd,tile_coord,nodata_forcing,met_force%Tair ) call check_forcing_nodata_2(N_catd,tile_coord,nodata_forcing,met_force%Qair ) call check_forcing_nodata_2(N_catd,tile_coord,nodata_forcing,met_force%Psurf ) call check_forcing_nodata_2(N_catd,tile_coord,nodata_forcing,met_force%Rainf_C) call check_forcing_nodata_2(N_catd,tile_coord,nodata_forcing,met_force%Rainf ) call check_forcing_nodata_2(N_catd,tile_coord,nodata_forcing,met_force%Snowf ) call check_forcing_nodata_2(N_catd,tile_coord,nodata_forcing,met_force%LWdown ) call check_forcing_nodata_2(N_catd,tile_coord,nodata_forcing,met_force%SWdown ) call check_forcing_nodata_2(N_catd,tile_coord,nodata_forcing,met_force%Wind ) ! do NOT call check_forcing_nodata_2() for "RefH" and "SWnet" (these are typically ! from GCM or DAS files and should not have any problems to begin with) ! reichle+qliu, 8 Oct 2008 ! for SWnet change nodata_forcing to nodata_generic -- reichle, 22 Jul 2010 tol = abs(nodata_forcing*nodata_tolfrac_generic) do i=1,N_catd if ( abs(met_force(i)%SWnet-nodata_forcing) < tol ) then met_force(i)%SWnet = nodata_generic end if end do end subroutine check_forcing_nodata ! ***************************************************************** subroutine check_forcing_nodata_2( N_catd, tile_coord, nodata_forcing, force_vec ) ! helper subroutine for check_forcing_nodata() ! ! If no-data-value is encountered, use value from "next" catchment, ! where "next" is next in line (not necessarily next in distance!). ! if that does not work, give up ! ! reichle, 13 Jun 2005 implicit none integer, intent(in) :: N_catd type(tile_coord_type), dimension(:), pointer :: tile_coord ! input real, intent(in) :: nodata_forcing real, dimension(N_catd), intent(inout) :: force_vec ! local variables real :: tol integer :: i, i_next, N_black ! set the following logical to .true. to generate a blacklist ! for a given forcing data set (it will end up in the file ! "fort.9999") logical, parameter :: create_blacklist = .false. character(len=*), parameter :: Iam = 'check_forcing_nodata_2' character(len=400) :: err_msg ! ------------------------------------------------------------ N_black = 0 tol = abs(nodata_forcing*nodata_tolfrac_generic) do i=1,N_catd ! no-data-value checks i_next = min(i+1,N_catd) if (abs(force_vec(i)-nodata_forcing)tol) then if (logit) write (logunit,*) 'forcing has no-data-value in tile ID = ', & tile_coord(i)%tile_id force_vec(i)=force_vec(i_next) else write (tmpstring10,*) tile_coord(i)%tile_id write (tmpstring40,*) tile_coord(i_next)%tile_id err_msg = 'forcing has no-data-value in tile ID = ' // & trim(tmpstring10) // ' and ' // trim(tmpstring40) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if end if end if end do if (create_blacklist) then write (logunit,*) '---------------------------------------------------------------' write (logunit,*) ' found N_black = ',N_black, ' tiles that should be blacklisted' err_msg = 'blacklist now in file fort.9999' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if end subroutine check_forcing_nodata_2 ! ****************************************************************** subroutine repair_forcing( N_catd, met_force, echo, tile_coord, fieldname, & unlimited_Qair, unlimited_LWdown ) ! check forcing for unphysical values, reset, print optional warning ! ! if optional input "fieldname" is present, only the corresponding ! forcing fields will be repaired ! ! reichle, 30 Mar 2004 ! reichle+qliu, 8 Oct 2008 - added optional input "unlimited_Qair" ! reichle, 11 Feb 2009 - added optional input "unlimited_LWdown" ! reichle, 2 Dec 2009 - eliminated "trim(field)==" to avoid memory ! leak due to Absoft 9 bug ! reichle, 24 Dec 2013 - limited number of warnings to N_tile_warn_max ! STILL TO DO: fix format specification for warning ! statements to avoid line breaks within ! a given statement implicit none integer, intent(in) :: N_catd type(met_force_type), dimension(N_catd), intent(inout) :: met_force logical, intent(in), optional :: echo, unlimited_Qair, unlimited_LWdown type(tile_coord_type), dimension(:), pointer, optional :: tile_coord ! in character(40), intent(in), optional :: fieldname ! local variables ! turn warnings off if warnings have been printed for N_warn_max tiles integer, parameter :: N_tile_warn_max = 20 ! min/max values for allowable range of forcing fields real, parameter :: min_Tair = 190. ! [K] real, parameter :: max_Tair = 340. ! [K] real, parameter :: max_PSurf = 115000. ! [Pa] ! slack parameters beyond which warnings are printed (if requested) ! specific humidity is often a tiny bit above saturated specific humidity ! (but sometimes much larger...) ! ! ALWAYS limit Qair <= Qair_sat ! ! *warning* ONLY for Qair <= tmpfac*Qair_sat ! ! tmpfac=1.02 corresponds approximately to relative humidity above 1.02 real, parameter :: tmpfac_warn_Qair = 1.2 ! 1.02 real, parameter :: tmpadd_warn_SWnet = 1.e-2 ! [W/m2] real, parameter :: tmpadd_warn_PAR = tmpadd_warn_SWnet real, parameter :: tmp_warn_Prec = 3.e-10 ! [m/s] (1.e-10m/s ~ 3mm/year) logical :: warn, unlimited_Qair_tmp, unlimited_LWdown_tmp, problem_tile integer :: i, kk real :: tmp_LWdown, min_LWdown, max_LWdown, Qair_sat, tmp_maxPar character(10) :: SWDN_MAX_string character(10) :: min_Tair_string, max_Tair_string, max_Psurf_string character(13) :: tmpstr13a, tmpstr13b character(16) :: tile_id_str, tmpstr16 character(40) :: field character(len=*), parameter :: Iam = 'repair_forcing' character(len=400) :: err_msg ! ------------------------------------------------------------ if (present(echo)) then if (present(tile_coord)) then warn = echo else call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'inconsistent optional args') end if else warn = .false. end if if (.not. logit) warn = .false. if (present(unlimited_Qair)) then unlimited_Qair_tmp = unlimited_Qair else unlimited_Qair_tmp = .false. end if if (present(unlimited_LWdown)) then unlimited_LWdown_tmp = unlimited_LWdown else unlimited_LWdown_tmp = .false. end if ! -------------------------------- if (present(fieldname)) then field = fieldname else field = 'all' end if ! -------------------------------- if (warn) then write (SWDN_MAX_string, '(f10.1)') SWDN_MAX write (min_Tair_string, '(f10.1)') min_Tair write (max_Tair_string, '(f10.1)') max_Tair write (max_PSurf_string,'(f10.1)') max_PSurf end if ! -------------------------------- kk = 0 ! counter for number of problematic tiles do i=1,N_catd problem_tile = .false. if (warn) write (tile_id_str,'(i16)') tile_coord(i)%tile_id ! convert integer to string ! precip and snowfall must be non-negative ! fix Rainf first, otherwise cannot use Rainf to constrain Rainf_C ! reichle, 11 Aug 2010 if (field(1:3)=='all' .or. field(1:7)=='Rainf ' ) then if ((warn) .and. (met_force(i)%Rainf < -tmp_warn_Prec)) then write (tmpstr13a,'(e13.5)') met_force(i)%Rainf ! convert real to string write (logunit,'(200A)') 'repair_forcing: Rainf < 0. in tile ID ' // & tile_id_str // ': met_force(i)%Rainf = ' // tmpstr13a problem_tile=.true. end if met_force(i)%Rainf = max( 0., met_force(i)%Rainf) end if if (field(1:3)=='all' .or. field(1:7)=='Rainf_C') then if ((warn) .and. (met_force(i)%Rainf_C < -tmp_warn_Prec)) then write (tmpstr13a,'(e13.5)') met_force(i)%Rainf_C ! convert real to string write (logunit,'(200A)') 'repair_forcing: Rainf_C < 0. in tile ID ' //& tile_id_str // ': met_force(i)%Rainf_C = ' // tmpstr13a problem_tile=.true. end if met_force(i)%Rainf_C = max( 0., met_force(i)%Rainf_C) ! make sure convective precip is less than total precip if ((warn) .and. (met_force(i)%Rainf+tmp_warn_Prec < met_force(i)%Rainf_C)) then write (tmpstr13a,'(e13.5)') met_force(i)%Rainf ! convert real to string write (tmpstr13b,'(e13.5)') met_force(i)%Rainf_C ! convert real to string write (logunit,'(200A)') 'repair_forcing: Rainf < Rainf_C in tile ID ' // & tile_id_str // ': met_force(i)%Rainf = ' // tmpstr13a // & ', met_force(i)%Rainf_C = ' // tmpstr13b problem_tile=.true. end if met_force(i)%Rainf_C = min( met_force(i)%Rainf, met_force(i)%Rainf_C) end if if (field(1:3)=='all' .or. field(1:7)=='Snowf ' ) then if ((warn) .and. (met_force(i)%Snowf < -tmp_warn_Prec)) then write (tmpstr13a,'(e13.5)') met_force(i)%Snowf ! convert real to string write (logunit,'(200A)') 'repair_forcing: Snowf < 0. in tile ID ' //& tile_id_str // ': met_force(i)%Snowf = ' // tmpstr13a problem_tile=.true. end if met_force(i)%Snowf = max( 0., met_force(i)%Snowf) end if ! -------------------------------- if (field(1:3)=='all' .or. field(1:7)=='Tair ') then ! temperatures must not be too low or too high ! NOTE: "warn" is turned on when repair_forcing is called first ! time after the forcing has been read from files if ((warn) .and. (met_force(i)%Tair < 190.)) then write (tmpstr13a,'(e13.5)') met_force(i)%Tair ! convert real to string write (logunit,'(200A)') & 'repair_forcing: Tair < '//min_Tair_string//' in tile ID ' // & tile_id_str // ': met_force(i)%Tair = ' // tmpstr13a call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'Tair too low') end if if ((warn) .and. (met_force(i)%Tair > 340.)) then write (tmpstr13a,'(e13.5)') met_force(i)%Tair ! convert real to string write (logunit,'(200A)') & 'repair_forcing: Tair > '//max_Tair_string//' in tile ID ' // & tile_id_str // ': met_force(i)%Tair = ' // tmpstr13a call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'Tair too high') end if end if ! -------------------------------- if (field(1:3)=='all' .or. field(1:7)=='Psurf ') then ! surface air pressure must not be too high if ((warn) .and. (met_force(i)%Psurf > max_Psurf)) then write (tmpstr13a,'(e13.5)') met_force(i)%PSurf ! convert real to string write (logunit,'(200A)') & 'repair_forcing: Psurf > '//max_PSurf_string//' in tile ID ' // & tile_id_str // ': met_force(i)%PSurf = ' // tmpstr13a call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'Psurf too high') end if end if ! -------------------------------- ! specific humidity must not be negative if (field(1:3)=='all' .or. field(1:7)=='Qair ') then if ((warn) .and. (met_force(i)%Qair < 0.)) then write (tmpstr13a,'(e13.5)') met_force(i)%Qair ! convert real to string write (logunit,'(200A)') 'repair_forcing: Qair < 0. in tile ID ' // & tile_id_str // ': met_force(i)%Qair = ' // tmpstr13a problem_tile=.true. end if met_force(i)%Qair = max( 0., met_force(i)%Qair ) ! specific humidity should not exceed saturated specific humidity if (.not. unlimited_Qair_tmp) then Qair_sat = MAPL_EQsat(met_force(i)%Tair, met_force(i)%Psurf ) if ((warn) .and. (met_force(i)%Qair > tmpfac_warn_Qair*Qair_sat)) then write (tmpstr13a,'(e13.5)') met_force(i)%Qair ! convert real to string write (tmpstr13b,'(e13.5)') met_force(i)%Qair/Qair_sat write (logunit,'(200A)') 'repair_forcing: Qair > Qair_sat in tile ID ' // & tile_id_str // ': met_force(i)%Qair = ' // tmpstr13a // & ', met_force(i)%Qair/Qair_sat = ' // tmpstr13b problem_tile=.true. end if met_force(i)%Qair = min(Qair_sat, met_force(i)%Qair) end if end if ! -------------------------------- ! wind must be positive ! (zero wind creates problem in turbulence calculations; ! warn only if wind is negative) if (field(1:3)=='all' .or. field(1:7)=='Wind ') then if ((warn) .and. (met_force(i)%Wind < 0.)) then write (tmpstr13a,'(e13.5)') met_force(i)%Wind ! convert real to string write (logunit,'(200A)') 'repair_forcing: Wind < 0. in tile ID ' //& tile_id_str // ': met_force(i)%Wind = ' // tmpstr13a problem_tile=.true. end if met_force(i)%Wind = max( 0.0001, met_force(i)%Wind) end if ! -------------------------------- if (field(1:3)=='all' .or. field(1:7)=='LWdown ') then if (.not. unlimited_LWdown_tmp) then ! make sure radiation is between min and max tmp_LWdown = met_force(i)%Tair*met_force(i)%Tair tmp_LWdown = stefan_boltzmann*tmp_LWdown*tmp_LWdown min_LWdown = LWDN_EMISS_MIN*tmp_LWdown max_LWdown = LWDN_EMISS_MAX*tmp_LWdown if ((warn) .and. (met_force(i)%LWdown < min_LWdown)) then write (tmpstr13a,'(e13.5)') met_force(i)%LWdown ! convert real to string write (tmpstr13b,'(e13.5)') min_LWdown write (logunit,'(200A)') 'repair_forcing: LWdown < min_LWdown in tile ID ' // & tile_id_str // ': met_force(i)%LWdown = ' // tmpstr13a // & ', min_LWdown = ' // tmpstr13b problem_tile=.true. end if met_force(i)%LWdown = max( min_LWdown, met_force(i)%LWdown) if ((warn) .and. (met_force(i)%LWdown > max_LWdown)) then write (tmpstr13a,'(e13.5)') met_force(i)%LWdown ! convert real to string write (tmpstr13b,'(e13.5)') max_LWdown write (logunit,'(200A)') 'repair_forcing: LWdown > max_LWdown in tile ID ' // & tile_id_str // ': met_force(i)%LWdown = ' // tmpstr13a // & ', max_LWdown = ' // tmpstr13b problem_tile=.true. end if met_force(i)%LWdown = min( max_LWdown, met_force(i)%LWdown) end if end if ! ----------------------------------- if (field(1:3)=='all' .or. field(1:7)=='SWdown ') then if ((warn) .and. (met_force(i)%SWdown < 0. )) then write (tmpstr13a,'(e13.5)') met_force(i)%SWdown ! convert real to string write (logunit,'(200A)') 'repair_forcing: SWdown < 0. in tile ID ' //& tile_id_str // ': met_force(i)%SWdown = ' // tmpstr13a problem_tile=.true. end if met_force(i)%SWdown = max( 0., met_force(i)%SWdown) ! shortwave must be less than solar constant ! (see also subroutine interpolate_to_timestep() for interpolated ! shortwave forcing) if ((warn) .and. (met_force(i)%SWdown > SWDN_MAX)) then write (tmpstr13a,'(e13.5)') met_force(i)%SWdown ! convert real to string write (logunit,'(200A)') 'repair_forcing: SWdown > ' // SWDN_MAX_string // & ' in tile ID ' // tile_id_str // ': met_force(i)%SWdown = ' // tmpstr13a problem_tile=.true. end if met_force(i)%SWdown = min( SWDN_MAX, met_force(i)%SWdown) end if ! ----------------------------------- ! SWnet is no-data-value for most forcing data sets (except MERRA, G5DAS) if(abs(met_force(i)%SWnet-nodata_generic)>nodata_tol_generic) then if (field(1:3)=='all' .or. field(1:7)=='SWnet ') then if ((warn) .and. (met_force(i)%SWnet < 0. )) then write (tmpstr13a,'(e13.5)') met_force(i)%SWnet ! convert real to string write (logunit,'(200A)') 'repair_forcing: SWnet < 0. in tile ID ' //& tile_id_str // ': met_force(i)%SWnet = ' // tmpstr13a problem_tile=.true. end if met_force(i)%SWnet = max( 0., met_force(i)%SWnet) ! net solar radiation must be less than solar incoming radiation if ( (warn) .and. & (met_force(i)%SWnet > met_force(i)%SWdown+tmpadd_warn_SWnet)) then write (tmpstr13a,'(e13.5)') met_force(i)%SWnet ! convert real to string write (tmpstr13b,'(e13.5)') met_force(i)%SWdown ! convert real to string write (logunit,'(200A)') 'repair_forcing: SWnet > SWdown in tile ID ' // & tile_id_str // ': met_force(i)%SWnet = ' // tmpstr13a // & ', met_force(i)%SWdown = ' // tmpstr13b problem_tile=.true. end if met_force(i)%SWnet = min( met_force(i)%SWnet,met_force(i)%SWdown) end if end if ! ----------------------------------- ! PARdffs is no-data-value for most forcing data sets (except MERRA, G5DAS) if(abs(met_force(i)%PARdffs-nodata_generic)>nodata_tol_generic) then if (field(1:3)=='all' .or. field(1:7)=='PARdffs') then if ((warn) .and. (met_force(i)%PARdffs < 0. )) then write (tmpstr13a,'(e13.5)') met_force(i)%PARdffs ! convert real to string write (logunit,'(200A)') 'repair_forcing: PARdffs < 0. in tile ID ' //& tile_id_str // ': met_force(i)%PARdffs = ' // tmpstr13a problem_tile=.true. end if met_force(i)%PARdffs = max( 0., met_force(i)%PARdffs) ! upper threshold for PARdffs = 0.6 of solar incoming radiation ! (after analysis if May and Dec cases from MERRA Scout) ! - reichle, 24 Feb 2009) ! updated to threshold of 0.8 after brief analysis of hourly MERRA data ! for Jul 2003 (courtesy of Greg Walker) ! - reichle, 20 Dec 2011 tmp_maxPar = 0.8*met_force(i)%SWdown if ((warn) .and. (met_force(i)%PARdffs > tmp_maxPar+tmpadd_warn_PAR )) then write (tmpstr13a,'(e13.5)') tmp_maxPar ! convert real to string write (tmpstr13b,'(e13.5)') met_force(i)%PARdffs ! convert real to string write (logunit,'(200A)') 'repair_forcing: PARdffs > ' // tmpstr13a // & ' in tile ID ' // tile_id_str // ': met_force(i)%PARdffs = ' // & tmpstr13b problem_tile=.true. end if met_force(i)%PARdffs = min( met_force(i)%PARdffs,tmp_maxPar) end if end if ! ----------------------------------- ! PARdrct is no-data-value for most forcing data sets (except MERRA, G5DAS) ! ! MUST "repair" PARdffs *before* PARdrct if(abs(met_force(i)%PARdrct-nodata_generic)>nodata_tol_generic) then if (field(1:3)=='all' .or. field(1:7)=='PARdrct') then if ((warn) .and. (met_force(i)%PARdrct < 0. )) then write (tmpstr13a,'(e13.5)') met_force(i)%PARdrct ! convert real to string write (logunit,'(200A)') 'repair_forcing: PARdrct < 0. in tile ID ' //& tile_id_str // ': met_force(i)%PARdrct = ' // tmpstr13a problem_tile=.true. end if met_force(i)%PARdrct = max( 0., met_force(i)%PARdrct) ! upper threshold for PARdrct = SWdown - PARdffs tmp_maxPar = met_force(i)%SWdown - met_force(i)%PARdffs if ((warn) .and. (met_force(i)%PARdrct > tmp_maxPar+tmpadd_warn_PAR )) then write (tmpstr13a,'(e13.5)') tmp_maxPar ! convert real to string write (tmpstr13b,'(e13.5)') met_force(i)%PARdrct ! convert real to string write (logunit,'(200A)') 'repair_forcing: PARdrct > ' // tmpstr13a // & ' in tile ID ' // tile_id_str // ': met_force(i)%PARdrct = ' // & tmpstr13b problem_tile=.true. end if met_force(i)%PARdrct = min( met_force(i)%PARdrct,tmp_maxPar) end if end if ! ------------------------------------------------------ ! ! count problematic tiles if (problem_tile) kk=kk+1 ! turn off warnings if number of problem tiles gets too large if ((warn) .and. (kk>N_tile_warn_max)) then warn = .false. ! turn off warnings for the remainder of the loop through tiles write (tmpstr16,'(i16)') kk ! convert integer to string write (logunit,'(200A)') & 'repair_forcing: turning OFF warnings after detecting ' // & trim(tmpstr16) // ' tiles with problematic forcing' end if ! ------------------------------------------------------ end do end subroutine repair_forcing ! ****************************************************************** type(date_time_type) function shift_forcing_date( met_tag, date_time ) ! shift date_time by years or days, useful for twin experiments ! ! examples: ! ! met_tag = "RedArk_ASCII_shift-1year" uses 1980 forcing for 1981 and so on ! met_tag = "RedArk_ASCII_shift+3day" uses Jan 4 forcing for Jan 1 and so on ! ! reichle, 6 Apr 2007 implicit none character(80) :: met_tag type(date_time_type) :: date_time, date_time_tmp integer :: is, ie, shift character(300) :: tmpstring300 ! initialize date_time_tmp = date_time ! check whether met_tag asks for shift is = index( met_tag, 'shift' ) if (is>0) then ! make sure this is only used for pre-specified forcing data sets ! for now permit use with all "RedArk" data sets (reichle, 6 Apr 2007) if (index(met_tag,'RedArk')==0) then tmpstring300 = 'shift_forcing_date(): Are you sure? ' // & 'If so, edit source code and recompile.' if (logit) write (logunit,*) tmpstring300 write(6,*) tmpstring300 write(0,*) tmpstring300 stop end if ie = index( met_tag, 'year') if (ie>0) then read (met_tag(is+5:ie-1),'(i1)') shift date_time_tmp%year = date_time%year + shift ! deal with leap year issues if ( is_leap_year(date_time_tmp%year) .or. & is_leap_year(date_time%year) ) then if (date_time%month==2 .and. date_time%day==29) & date_time_tmp%day = 28 call get_dofyr_pentad( date_time_tmp ) end if end if ie = index( met_tag, 'day') if (ie>0) then read (met_tag(is+5:ie-1),'(i1)') shift call augment_date_time( 86400*shift, date_time_tmp ) end if end if shift_forcing_date = date_time_tmp end function shift_forcing_date ! ****************************************************************** end module CLSM_ensdrv_force_routines ! ------------------------------------------------------------------- #if 0 program ut_parse_G5DAS_met_tag implicit none integer, parameter :: N_met_tag_in=6 ! "in" character(200) :: met_path_in character( 80) :: met_tag_in ! "out" character(200) :: met_path_default, met_path_prec character( 80) :: met_tag_out logical :: use_prec_corr logical :: use_Predictor ! other integer :: ii character( 80), dimension(N_met_tag_in) :: met_tag_in_vec met_path_in = 'mymetpathin' met_tag_in_vec = (/ & 'gcmexpname' , & 'gcmexpname__Nx+-' , & 'gcmexpname__precCORRPREC', & 'gcmexpname__precCORRPREC__Nx+-' , & 'gcmexpname__Nx+-__precCORRPREC' , & 'gcmexpname__qrecCORRPREC__Nx+-' /) do ii=1,N_met_tag_in met_tag_in = met_tag_in_vec(ii) write (*,*) '----------------------------------------------' write (*,*) 'ii = ', ii write (*,*) 'met_tag_in = ', trim(met_tag_in) call parse_G5DAS_met_tag( met_path_in, met_tag_in, & met_path_default, met_path_prec, met_tag_out, use_prec_corr, use_Predictor ) write (*,*) 'met_path_default = ', trim(met_path_default) write (*,*) 'met_path_prec = ', trim(met_path_prec) write (*,*) 'met_tag_out = ', trim(met_tag_out) write (*,*) 'use_prec_corr = ', use_prec_corr write (*,*) 'use_Predictor = ', use_Predictor end do end program ut_parse_G5DAS_met_tag #endif #if 0 program test_shift_forcing_date use clsm_ensdrv_functions use date_time_util type(date_time_type) :: date_time, date_time_tmp integer :: dtstep, iter, i character(80) :: met_tag ! ------- met_tag = 'RedArk_OSSE_shift+3day' met_tag = 'Princeton_shift+3day' dtstep = 21600 iter = 120 date_time%year = 1992 ! 4-digit year date_time%month = 2 ! month in year date_time%day = 1 ! day in month date_time%hour = 3 ! hour of day date_time%min = 0 ! minute of hour date_time%sec = 0 ! seconds of minute date_time%pentad = -9999 ! pentad of year date_time%dofyr = -9999 ! day of year do i=1,iter call augment_date_time(dtstep, date_time) date_time_tmp = shift_forcing_date(met_tag, date_time) write (*,'(16i5)') date_time, date_time_tmp end do end program test_shift_forcing_date #endif ! *********** EOF **************************************************