From af16f2b84ab809d471ce03b7d8c9402964a7095c Mon Sep 17 00:00:00 2001 From: Qing Liu Date: Fri, 8 Aug 2025 12:52:07 -0400 Subject: [PATCH 01/16] add get_GEOSs2s_v3() met forcing --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 321 ++++++++++++++++++++++++- 1 file changed, 317 insertions(+), 4 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 607dbcf0..8a7b5f23 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -320,9 +320,14 @@ subroutine get_forcing( date_time, force_dtstep, met_path, met_tag, & unlimited_Qair = .true. unlimited_LWdown = .true. + elseif (index(met_tag(1:9), 'GEOSs2sv3')/=0) then + + call get_GEOSs2s_v3( date_time_tmp, met_path, met_tag, N_catd, tile_coord, & + MET_HINTERP, met_force_obs_tile_new, nodata_forcing) + elseif (index(met_tag(1:7), 'GEOSs2s')/=0) then - call get_GEOSs2s( date_time_tmp, met_path, met_tag, N_catd, tile_coord, & + call get_GEOSs2s_v2( date_time_tmp, met_path, met_tag, N_catd, tile_coord, & MET_HINTERP, met_force_obs_tile_new, nodata_forcing, PAR_available) else ! assume forcing from GEOS5 GCM ("DAS" or "MERRA") output @@ -2643,7 +2648,7 @@ end subroutine get_Viviana_OK_precip ! ************************************************************************ - subroutine get_GEOSs2s(date_time, met_path, met_tag, N_catd, tile_coord, & + subroutine get_GEOSs2s_v2(date_time, met_path, met_tag, N_catd, tile_coord, & met_hinterp, met_force_new, nodata_forcing, PAR_available) ! read forcing derived from GEOS S2S output and map to tile space @@ -2758,7 +2763,7 @@ subroutine get_GEOSs2s(date_time, met_path, met_tag, N_catd, tile_coord, & logical :: FCST = .false. logical :: AODAS = .false. - character(len=*), parameter :: Iam = 'get_GEOSs2s' + character(len=*), parameter :: Iam = 'get_GEOSs2s_v2' character(len=400) :: err_msg ! -------------------------------------------------------------------- @@ -3046,7 +3051,315 @@ subroutine get_GEOSs2s(date_time, met_path, met_tag, N_catd, tile_coord, & deallocate(force_array) deallocate(GEOSgcm_name) - end subroutine get_GEOSs2s + end subroutine get_GEOSs2s_v2 + + ! ************************************************************************ + + subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & + met_hinterp, met_force_new, nodata_forcing) + + ! read forcing derived from GEOS S2Sv3 output and map to tile space + ! (using nearest neighbor or bilinear interpolation) + ! + ! forcing derived through post-processing of daily average output from S2S + ! hindcasts/forecasts ("FCST") or the "AODAS" used for initialization + ! (see doc/README.MetForcing_and_BCS.md) + ! + ! implementation follows LDASsa subroutines get_Princeton_netcdf() and get_GEOS(), + ! fzeng, 24 Jun 2019 + ! + ! jkolassa,jmpark,reichle, 10 May - 14 June 2021: + ! modified for GEOSldas; added AODAS; met_tag encodes ID of ensemble + ! member and initial month/day + ! modified to compute precipitation components from total precipitation + ! and air temperature (pre-processing of precipitation components was faulty) + ! qliu, Aug. 2025 + ! modified for S2Sv3 fcst data, remove AODAS option, remove PAR_available option + + use netcdf + implicit none + include 'mpif.h' + + type(date_time_type), intent(in) :: date_time + + character(*), intent(in) :: met_path + character(*), 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(:), intent(inout) :: met_force_new + + real, intent(out) :: nodata_forcing + + ! Hindcast grid and netcdf parameters + + integer, parameter :: GEOSgcm_grid_N_lon = 720 + integer, parameter :: GEOSgcm_grid_N_lat = 361 + real, parameter :: GEOSgcm_grid_dlon = 0.5 + real, parameter :: GEOSgcm_grid_dlat = 0.5 + real, parameter :: GEOSgcm_grid_ll_lon = -180. - GEOSgcm_grid_dlon/2. + real, parameter :: GEOSgcm_grid_ll_lat = -90. - GEOSgcm_grid_dlat/2. + + real, parameter :: nodata_GEOSgcm = -9999. + + integer, parameter :: dt_GEOSgcm_in_hours_FCST = 3 + + integer, parameter :: N_GEOSgcm_vars_FCST = 10 + + ! variable names in "FCST" forcing files match those in MERRA-2(SWGDN, LWGAB) + ! and S2Sv3 (other variables) file specs + + character(40), dimension(N_GEOSgcm_vars_FCST) :: GEOSgcm_name_FCST = & + (/ & + 'PLS ', & ! 1 - flux, largescale_liquid_precipitation, kg m-2 s-1 + 'PCU ', & ! 2 - flux, convective_liquid_precipitation, kg m-2 s-1 + 'SNO ', & ! 3 - flux, snowfall_precipitation, kg m-2 s-1 + 'LWGAB ', & ! 4 - flux, surface_absorbed_longwave_radiation, W m-2 + 'SWGDN ', & ! 5 - flux, surface_incoming_shortwave_flux, W m-2 + 'PS ', & ! 6 - state, surface_pressure, Pa + 'QA ', & ! 7 - state, surface_specific_humidity, kg kg-1 + 'TA ', & ! 8 - state, surface_air_temperature, K + 'SPEED ', & ! 9 - state, surface_wind_speed, m s-1 + 'HLML ' & ! 10 - state, surface_layer_height, m + /) + + ! local variables + + character(40), dimension(:), allocatable :: GEOSgcm_name + + integer :: N_GEOSgcm_vars, dt_GEOSgcm_in_hours, N_lon_tmp, N_lat_tmp + + real :: fnbr(2,2) + + integer, pointer :: i1(:), i2(:), j1(:), j2(:) + real, pointer :: x1(:), x2(:), y1(:), y2(:) + + real, dimension(:,:), allocatable :: force_array + + integer, dimension(3) :: iicount, iistart + integer :: k, hours_in_month, GEOSgcm_var, nciv_data + integer :: fid, rc, nv_id, status + + real :: tol, this_lon, this_lat, Tair_tmp, Totprec_tmp + + character( 8) :: init_YYYYMMDD + character( 4) :: YYYY, ens_num + character( 2) :: MM, DD + character(300) :: fname + + character( 10) :: lat_str = 'latitude' + character( 10) :: lon_str = 'longitude' + + logical :: FCST = .false. + + character(len=*), parameter :: Iam = 'get_GEOSs2s_v3' + character(len=400) :: err_msg + + ! -------------------------------------------------------------------- + ! + ! preparations + + ! parse met_tag + + ! 12345678901234567890123456789 + ! GEOSs2sv2FCST__ensX__YYYYMMDD + ! GEOSs2sv2AODAS + + if (index(met_tag(10:13), 'FCST')/=0) then + + FCST = .true. + + dt_GEOSgcm_in_hours = dt_GEOSgcm_in_hours_FCST + N_GEOSgcm_vars = N_GEOSgcm_vars_FCST + + allocate(GEOSgcm_name(N_GEOSgcm_vars)) + + GEOSgcm_name = GEOSgcm_name_FCST + + ! parse ensemble and initialization month from met_tag + + ens_num = trim(met_tag(16:19)) + init_YYYYMMDD = trim(met_tag(22:29)) + + else + + err_msg = "unknown met_tag" + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + + end if + + nodata_forcing = nodata_GEOSgcm + + 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 + + ! set lon index + + iistart(1) = 1 + iicount(1) = GEOSgcm_grid_N_lon + + ! set lat index + iistart(2) = 1 + iicount(2) = GEOSgcm_grid_N_lat + + ! get time index + + if ( (date_time%min/=0) .or. (date_time%sec/=0) .or. & + (mod(date_time%hour, dt_GEOSgcm_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 + + iistart(3) = hours_in_month / dt_GEOSgcm_in_hours + 1 + iicount(3) = 1 + + ! ---------------------------------------------------------------- + ! + ! open input file + + if (FCST) then + fname = trim(met_path) // '/' // init_YYYYMMDD // '/' // ens_num // '/G +EOSS2S3.' // YYYY // MM // '.nc4' + endif + + call GEOS_openfile(FileOpenedHash,fname,fid,tile_coord,met_hinterp,rc,lat_str,lon_str) + if (rc<0) then + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error opening file') + endif + + ! assign indices from met forcing interpolation + + i1=>local_info%i1 + i2=>local_info%i2 + j1=>local_info%j1 + j2=>local_info%j2 + x1=>local_info%x1 + x2=>local_info%x2 + y1=>local_info%y1 + y2=>local_info%y2 + + ! allocate force_array + + allocate(force_array(N_catd,N_GEOSgcm_vars)) + force_array = nodata_forcing + + ! loop over variables + + do GEOSgcm_var = 1,N_GEOSgcm_vars + + ! init shared memory + N_lon_tmp = -1 + N_lat_tmp = -1 + if (associated(ptrShForce)) then + N_lon_tmp = size(ptrShForce,1) + N_lat_tmp = size(ptrShForce,2) + endif + if( (N_lon_tmp /= GEOSgcm_grid_N_lon) .or. & + (N_lat_tmp /= GEOSgcm_grid_N_lat) ) then + call MAPL_SyncSharedMemory(rc=status) + VERIFY_(status) + if (associated(ptrShForce)) then + call MAPL_DeallocNodeArray(ptrShForce,rc=status) + VERIFY_(status) + endif + call MAPL_AllocateShared(ptrShForce,(/GEOSgcm_grid_N_lon,GEOSgcm_grid_N_lat/),Tr +ansRoot= .true.,rc=status) + VERIFY_(status) + end if + + call MAPL_SyncSharedMemory(rc=status) + VERIFY_(status) + + ! read variable from netcdf file + if (MAPL_AmNodeRoot .or. (.not. MAPL_ShmInitialized)) then + rc= NF90_INQ_VARID( fid, trim(GEOSgcm_name(GEOSgcm_var)), nv_id) + _ASSERT( rc == nf90_noerr, "nf90 error") + rc= NF90_GET_VAR( fid, nv_id, ptrShForce, start=iistart,count=iicount) + end if + + call MAPL_SyncSharedMemory(rc=status) + + ! map variable array to force array using chosen met interpolation method + + select case (MET_HINTERP) + + case(0) ! nearest neighbor interpolation + + do k = 1, N_catd + force_array(k, GEOSgcm_var) = ptrShForce(i1(k), j1(k)) + 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) = ptrShForce(i1(k),j1(k)) + fnbr(1,2) = ptrShForce(i1(k),j2(k)) + fnbr(2,1) = ptrShForce(i2(k),j1(k)) + fnbr(2,2) = ptrShForce(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 + + case default + + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unsupported MET_HINTERP option') + + end select + end do ! GEOSgcm_var + + ! ---------------------------------------------------------------- + + ! convert variables and units of force_array to match met_force_type, + ! put into structure + + ! from GEOSgcm files: + ! + ! FCST/Hindcast + ! + ! force_array(:, 1) = PLS kg/m2/s ([gauge-corr]_largescale_rainfall + ! force_array(:, 2) = PCU kg/m2/s ([gauge-corr]_convective_rainfall + ! force_array(:, 3) = SNO kg/m2/s ([gauge-corr]_snowfall + ! force_array(:, 2) = LWGAB W/m2 (surface_absorbed_longwave_radiation) + ! force_array(:, 3) = SWGDN W/m2 (surface_incoming_shortwave_flux) + ! force_array(:, 4) = PS Pa (surface_pressure) + ! force_array(:, 5) = QA kg/kg (surface_specific_humidity) + ! force_array(:, 6) = TA K (surface_air_temperature) + ! force_array(:, 7) = SPEED m/s (surface_wind_speed) + ! force_array(:, 8) = HLML m (surface_layer_height) + + if (FCST) then + met_force_new%Rainf = force_array(:, 1) + force_array(:, 2) + met_force_new%Rainf_C = 0. + met_force_new%Snowf = force_array(:, 3) + met_force_new%LWdown = force_array(:, 4) + met_force_new%SWdown = 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) + met_force_new%RefH = force_array(:,10) + end if + + deallocate(force_array) + deallocate(GEOSgcm_name) + + end subroutine get_GEOSs2s_v3 + ! ************************************************************************* From 3d7d160a1c1be81ffd80d04d9d514e5c2ddafd21 Mon Sep 17 00:00:00 2001 From: Qing Liu Date: Fri, 8 Aug 2025 15:49:59 -0400 Subject: [PATCH 02/16] update get_GEOSs2s_v3 --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 8a7b5f23..b09e31e3 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -3164,8 +3164,8 @@ subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & ! parse met_tag ! 12345678901234567890123456789 - ! GEOSs2sv2FCST__ensX__YYYYMMDD - ! GEOSs2sv2AODAS + ! GEOSs2sv3FCST__ensX__YYYYMMDD + ! GEOSs2sv3AODAS if (index(met_tag(10:13), 'FCST')/=0) then @@ -3228,8 +3228,7 @@ subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & ! open input file if (FCST) then - fname = trim(met_path) // '/' // init_YYYYMMDD // '/' // ens_num // '/G -EOSS2S3.' // YYYY // MM // '.nc4' + fname = trim(met_path) // '/' // init_YYYYMMDD // '/' // ens_num // '/GEOSS2S3.' // YYYY // MM // '.nc4' endif call GEOS_openfile(FileOpenedHash,fname,fid,tile_coord,met_hinterp,rc,lat_str,lon_str) @@ -3272,8 +3271,7 @@ subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & call MAPL_DeallocNodeArray(ptrShForce,rc=status) VERIFY_(status) endif - call MAPL_AllocateShared(ptrShForce,(/GEOSgcm_grid_N_lon,GEOSgcm_grid_N_lat/),Tr -ansRoot= .true.,rc=status) + call MAPL_AllocateShared(ptrShForce,(/GEOSgcm_grid_N_lon,GEOSgcm_grid_N_lat/),TransRoot= .true.,rc=status) VERIFY_(status) end if From 59f0317ca6f5574483dfe1298e30eb5df56fb842 Mon Sep 17 00:00:00 2001 From: Qing Liu Date: Sun, 10 Aug 2025 12:51:49 -0400 Subject: [PATCH 03/16] update get_s2s_v3() --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 52 ++++++++++++++++++-------- 1 file changed, 36 insertions(+), 16 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index b09e31e3..8b7f341f 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -3102,7 +3102,7 @@ subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & real, parameter :: GEOSgcm_grid_ll_lon = -180. - GEOSgcm_grid_dlon/2. real, parameter :: GEOSgcm_grid_ll_lat = -90. - GEOSgcm_grid_dlat/2. - real, parameter :: nodata_GEOSgcm = -9999. + real, parameter :: nodata_GEOSgcm = 1.e15 integer, parameter :: dt_GEOSgcm_in_hours_FCST = 3 @@ -3138,10 +3138,12 @@ subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & real, dimension(:,:), allocatable :: force_array - integer, dimension(3) :: iicount, iistart + integer, dimension(3) :: iicount, iistart, iistart_tmp integer :: k, hours_in_month, GEOSgcm_var, nciv_data - integer :: fid, rc, nv_id, status - + integer :: fid, rc, nv_id, status, tdimid + integer, parameter :: odd_months(7) = [1, 3, 5, 7, 8, 10, 12] + integer :: ttot_in_month, time_len + real :: tol, this_lon, this_lat, Tair_tmp, Totprec_tmp character( 8) :: init_YYYYMMDD @@ -3149,8 +3151,8 @@ subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & character( 2) :: MM, DD character(300) :: fname - character( 10) :: lat_str = 'latitude' - character( 10) :: lon_str = 'longitude' + character( 3) :: lat_str = 'lat' + character( 3) :: lon_str = 'lon' logical :: FCST = .false. @@ -3219,10 +3221,21 @@ subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & endif hours_in_month = (date_time%day-1)*24 + date_time%hour + + if (any(date_time%month == odd_months)) then + ttot_in_month = 31 * 24 / dt_GEOSgcm_in_hours + elseif (date_time%month == 2) then + if (is_leap_year(date_time%year)) then + ttot_in_month = 29 * 24 / dt_GEOSgcm_in_hours + else + ttot_in_month = 28 * 24 / dt_GEOSgcm_in_hours + endif + else + ttot_in_month = 30 * 24 / dt_GEOSgcm_in_hours + endif iistart(3) = hours_in_month / dt_GEOSgcm_in_hours + 1 iicount(3) = 1 - ! ---------------------------------------------------------------- ! ! open input file @@ -3282,7 +3295,15 @@ subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & if (MAPL_AmNodeRoot .or. (.not. MAPL_ShmInitialized)) then rc= NF90_INQ_VARID( fid, trim(GEOSgcm_name(GEOSgcm_var)), nv_id) _ASSERT( rc == nf90_noerr, "nf90 error") - rc= NF90_GET_VAR( fid, nv_id, ptrShForce, start=iistart,count=iicount) + rc= nf90_inq_dimid(fid,"time",tdimid) + _ASSERT( rc == nf90_noerr, "nf90 error") + rc= nf90_Inquire_Dimension(fid, tdimid, len=time_len) + _ASSERT( rc == nf90_noerr, "nf90 error") + + iistart_tmp = iistart + iistart_tmp(3) = iistart(3) - (ttot_in_month-time_len) + ! some monthly file contains only partial temporal record + rc= NF90_GET_VAR( fid, nv_id, ptrShForce, start=iistart_tmp,count=iicount) end if call MAPL_SyncSharedMemory(rc=status) @@ -3319,7 +3340,6 @@ subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & end select end do ! GEOSgcm_var - ! ---------------------------------------------------------------- ! convert variables and units of force_array to match met_force_type, @@ -3332,13 +3352,13 @@ subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & ! force_array(:, 1) = PLS kg/m2/s ([gauge-corr]_largescale_rainfall ! force_array(:, 2) = PCU kg/m2/s ([gauge-corr]_convective_rainfall ! force_array(:, 3) = SNO kg/m2/s ([gauge-corr]_snowfall - ! force_array(:, 2) = LWGAB W/m2 (surface_absorbed_longwave_radiation) - ! force_array(:, 3) = SWGDN W/m2 (surface_incoming_shortwave_flux) - ! force_array(:, 4) = PS Pa (surface_pressure) - ! force_array(:, 5) = QA kg/kg (surface_specific_humidity) - ! force_array(:, 6) = TA K (surface_air_temperature) - ! force_array(:, 7) = SPEED m/s (surface_wind_speed) - ! force_array(:, 8) = HLML m (surface_layer_height) + ! force_array(:, 4) = LWGAB W/m2 (surface_absorbed_longwave_radiation) + ! force_array(:, 5) = SWGDN W/m2 (surface_incoming_shortwave_flux) + ! force_array(:, 6) = PS Pa (surface_pressure) + ! force_array(:, 7) = QA kg/kg (surface_specific_humidity) + ! force_array(:, 8) = TA K (surface_air_temperature) + ! force_array(:, 9) = SPEED m/s (surface_wind_speed) + ! force_array(:,10) = HLML m (surface_layer_height) if (FCST) then met_force_new%Rainf = force_array(:, 1) + force_array(:, 2) From ff95e4eb728bed44c14a3bba56757d1638f0ba7d Mon Sep 17 00:00:00 2001 From: Qing Liu Date: Wed, 13 Aug 2025 15:44:36 -0400 Subject: [PATCH 04/16] use get_GEOS() for S2Sv3 --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 520 ++++++++----------------- 1 file changed, 154 insertions(+), 366 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 8b7f341f..56a06588 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -38,7 +38,8 @@ module LDAS_ForceMod datetime_lt_refdatetime, & datetime_le_refdatetime, & is_leap_year, & - get_dofyr_pentad + get_dofyr_pentad, & + days_in_month use LDAS_ExceptionsMod, ONLY: & ldas_abort, & @@ -320,11 +321,6 @@ subroutine get_forcing( date_time, force_dtstep, met_path, met_tag, & unlimited_Qair = .true. unlimited_LWdown = .true. - elseif (index(met_tag(1:9), 'GEOSs2sv3')/=0) then - - call get_GEOSs2s_v3( date_time_tmp, met_path, met_tag, N_catd, tile_coord, & - MET_HINTERP, met_force_obs_tile_new, nodata_forcing) - elseif (index(met_tag(1:7), 'GEOSs2s')/=0) then call get_GEOSs2s_v2( date_time_tmp, met_path, met_tag, N_catd, tile_coord, & @@ -3053,332 +3049,6 @@ subroutine get_GEOSs2s_v2(date_time, met_path, met_tag, N_catd, tile_coord, & end subroutine get_GEOSs2s_v2 - ! ************************************************************************ - - subroutine get_GEOSs2s_v3(date_time, met_path, met_tag, N_catd, tile_coord, & - met_hinterp, met_force_new, nodata_forcing) - - ! read forcing derived from GEOS S2Sv3 output and map to tile space - ! (using nearest neighbor or bilinear interpolation) - ! - ! forcing derived through post-processing of daily average output from S2S - ! hindcasts/forecasts ("FCST") or the "AODAS" used for initialization - ! (see doc/README.MetForcing_and_BCS.md) - ! - ! implementation follows LDASsa subroutines get_Princeton_netcdf() and get_GEOS(), - ! fzeng, 24 Jun 2019 - ! - ! jkolassa,jmpark,reichle, 10 May - 14 June 2021: - ! modified for GEOSldas; added AODAS; met_tag encodes ID of ensemble - ! member and initial month/day - ! modified to compute precipitation components from total precipitation - ! and air temperature (pre-processing of precipitation components was faulty) - ! qliu, Aug. 2025 - ! modified for S2Sv3 fcst data, remove AODAS option, remove PAR_available option - - use netcdf - implicit none - include 'mpif.h' - - type(date_time_type), intent(in) :: date_time - - character(*), intent(in) :: met_path - character(*), 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(:), intent(inout) :: met_force_new - - real, intent(out) :: nodata_forcing - - ! Hindcast grid and netcdf parameters - - integer, parameter :: GEOSgcm_grid_N_lon = 720 - integer, parameter :: GEOSgcm_grid_N_lat = 361 - real, parameter :: GEOSgcm_grid_dlon = 0.5 - real, parameter :: GEOSgcm_grid_dlat = 0.5 - real, parameter :: GEOSgcm_grid_ll_lon = -180. - GEOSgcm_grid_dlon/2. - real, parameter :: GEOSgcm_grid_ll_lat = -90. - GEOSgcm_grid_dlat/2. - - real, parameter :: nodata_GEOSgcm = 1.e15 - - integer, parameter :: dt_GEOSgcm_in_hours_FCST = 3 - - integer, parameter :: N_GEOSgcm_vars_FCST = 10 - - ! variable names in "FCST" forcing files match those in MERRA-2(SWGDN, LWGAB) - ! and S2Sv3 (other variables) file specs - - character(40), dimension(N_GEOSgcm_vars_FCST) :: GEOSgcm_name_FCST = & - (/ & - 'PLS ', & ! 1 - flux, largescale_liquid_precipitation, kg m-2 s-1 - 'PCU ', & ! 2 - flux, convective_liquid_precipitation, kg m-2 s-1 - 'SNO ', & ! 3 - flux, snowfall_precipitation, kg m-2 s-1 - 'LWGAB ', & ! 4 - flux, surface_absorbed_longwave_radiation, W m-2 - 'SWGDN ', & ! 5 - flux, surface_incoming_shortwave_flux, W m-2 - 'PS ', & ! 6 - state, surface_pressure, Pa - 'QA ', & ! 7 - state, surface_specific_humidity, kg kg-1 - 'TA ', & ! 8 - state, surface_air_temperature, K - 'SPEED ', & ! 9 - state, surface_wind_speed, m s-1 - 'HLML ' & ! 10 - state, surface_layer_height, m - /) - - ! local variables - - character(40), dimension(:), allocatable :: GEOSgcm_name - - integer :: N_GEOSgcm_vars, dt_GEOSgcm_in_hours, N_lon_tmp, N_lat_tmp - - real :: fnbr(2,2) - - integer, pointer :: i1(:), i2(:), j1(:), j2(:) - real, pointer :: x1(:), x2(:), y1(:), y2(:) - - real, dimension(:,:), allocatable :: force_array - - integer, dimension(3) :: iicount, iistart, iistart_tmp - integer :: k, hours_in_month, GEOSgcm_var, nciv_data - integer :: fid, rc, nv_id, status, tdimid - integer, parameter :: odd_months(7) = [1, 3, 5, 7, 8, 10, 12] - integer :: ttot_in_month, time_len - - real :: tol, this_lon, this_lat, Tair_tmp, Totprec_tmp - - character( 8) :: init_YYYYMMDD - character( 4) :: YYYY, ens_num - character( 2) :: MM, DD - character(300) :: fname - - character( 3) :: lat_str = 'lat' - character( 3) :: lon_str = 'lon' - - logical :: FCST = .false. - - character(len=*), parameter :: Iam = 'get_GEOSs2s_v3' - character(len=400) :: err_msg - - ! -------------------------------------------------------------------- - ! - ! preparations - - ! parse met_tag - - ! 12345678901234567890123456789 - ! GEOSs2sv3FCST__ensX__YYYYMMDD - ! GEOSs2sv3AODAS - - if (index(met_tag(10:13), 'FCST')/=0) then - - FCST = .true. - - dt_GEOSgcm_in_hours = dt_GEOSgcm_in_hours_FCST - N_GEOSgcm_vars = N_GEOSgcm_vars_FCST - - allocate(GEOSgcm_name(N_GEOSgcm_vars)) - - GEOSgcm_name = GEOSgcm_name_FCST - - ! parse ensemble and initialization month from met_tag - - ens_num = trim(met_tag(16:19)) - init_YYYYMMDD = trim(met_tag(22:29)) - - else - - err_msg = "unknown met_tag" - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - - end if - - nodata_forcing = nodata_GEOSgcm - - 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 - - ! set lon index - - iistart(1) = 1 - iicount(1) = GEOSgcm_grid_N_lon - - ! set lat index - iistart(2) = 1 - iicount(2) = GEOSgcm_grid_N_lat - - ! get time index - - if ( (date_time%min/=0) .or. (date_time%sec/=0) .or. & - (mod(date_time%hour, dt_GEOSgcm_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 - - if (any(date_time%month == odd_months)) then - ttot_in_month = 31 * 24 / dt_GEOSgcm_in_hours - elseif (date_time%month == 2) then - if (is_leap_year(date_time%year)) then - ttot_in_month = 29 * 24 / dt_GEOSgcm_in_hours - else - ttot_in_month = 28 * 24 / dt_GEOSgcm_in_hours - endif - else - ttot_in_month = 30 * 24 / dt_GEOSgcm_in_hours - endif - - iistart(3) = hours_in_month / dt_GEOSgcm_in_hours + 1 - iicount(3) = 1 - ! ---------------------------------------------------------------- - ! - ! open input file - - if (FCST) then - fname = trim(met_path) // '/' // init_YYYYMMDD // '/' // ens_num // '/GEOSS2S3.' // YYYY // MM // '.nc4' - endif - - call GEOS_openfile(FileOpenedHash,fname,fid,tile_coord,met_hinterp,rc,lat_str,lon_str) - if (rc<0) then - call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error opening file') - endif - - ! assign indices from met forcing interpolation - - i1=>local_info%i1 - i2=>local_info%i2 - j1=>local_info%j1 - j2=>local_info%j2 - x1=>local_info%x1 - x2=>local_info%x2 - y1=>local_info%y1 - y2=>local_info%y2 - - ! allocate force_array - - allocate(force_array(N_catd,N_GEOSgcm_vars)) - force_array = nodata_forcing - - ! loop over variables - - do GEOSgcm_var = 1,N_GEOSgcm_vars - - ! init shared memory - N_lon_tmp = -1 - N_lat_tmp = -1 - if (associated(ptrShForce)) then - N_lon_tmp = size(ptrShForce,1) - N_lat_tmp = size(ptrShForce,2) - endif - if( (N_lon_tmp /= GEOSgcm_grid_N_lon) .or. & - (N_lat_tmp /= GEOSgcm_grid_N_lat) ) then - call MAPL_SyncSharedMemory(rc=status) - VERIFY_(status) - if (associated(ptrShForce)) then - call MAPL_DeallocNodeArray(ptrShForce,rc=status) - VERIFY_(status) - endif - call MAPL_AllocateShared(ptrShForce,(/GEOSgcm_grid_N_lon,GEOSgcm_grid_N_lat/),TransRoot= .true.,rc=status) - VERIFY_(status) - end if - - call MAPL_SyncSharedMemory(rc=status) - VERIFY_(status) - - ! read variable from netcdf file - if (MAPL_AmNodeRoot .or. (.not. MAPL_ShmInitialized)) then - rc= NF90_INQ_VARID( fid, trim(GEOSgcm_name(GEOSgcm_var)), nv_id) - _ASSERT( rc == nf90_noerr, "nf90 error") - rc= nf90_inq_dimid(fid,"time",tdimid) - _ASSERT( rc == nf90_noerr, "nf90 error") - rc= nf90_Inquire_Dimension(fid, tdimid, len=time_len) - _ASSERT( rc == nf90_noerr, "nf90 error") - - iistart_tmp = iistart - iistart_tmp(3) = iistart(3) - (ttot_in_month-time_len) - ! some monthly file contains only partial temporal record - rc= NF90_GET_VAR( fid, nv_id, ptrShForce, start=iistart_tmp,count=iicount) - end if - - call MAPL_SyncSharedMemory(rc=status) - - ! map variable array to force array using chosen met interpolation method - - select case (MET_HINTERP) - - case(0) ! nearest neighbor interpolation - - do k = 1, N_catd - force_array(k, GEOSgcm_var) = ptrShForce(i1(k), j1(k)) - 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) = ptrShForce(i1(k),j1(k)) - fnbr(1,2) = ptrShForce(i1(k),j2(k)) - fnbr(2,1) = ptrShForce(i2(k),j1(k)) - fnbr(2,2) = ptrShForce(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 - - case default - - call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unsupported MET_HINTERP option') - - end select - end do ! GEOSgcm_var - ! ---------------------------------------------------------------- - - ! convert variables and units of force_array to match met_force_type, - ! put into structure - - ! from GEOSgcm files: - ! - ! FCST/Hindcast - ! - ! force_array(:, 1) = PLS kg/m2/s ([gauge-corr]_largescale_rainfall - ! force_array(:, 2) = PCU kg/m2/s ([gauge-corr]_convective_rainfall - ! force_array(:, 3) = SNO kg/m2/s ([gauge-corr]_snowfall - ! force_array(:, 4) = LWGAB W/m2 (surface_absorbed_longwave_radiation) - ! force_array(:, 5) = SWGDN W/m2 (surface_incoming_shortwave_flux) - ! force_array(:, 6) = PS Pa (surface_pressure) - ! force_array(:, 7) = QA kg/kg (surface_specific_humidity) - ! force_array(:, 8) = TA K (surface_air_temperature) - ! force_array(:, 9) = SPEED m/s (surface_wind_speed) - ! force_array(:,10) = HLML m (surface_layer_height) - - if (FCST) then - met_force_new%Rainf = force_array(:, 1) + force_array(:, 2) - met_force_new%Rainf_C = 0. - met_force_new%Snowf = force_array(:, 3) - met_force_new%LWdown = force_array(:, 4) - met_force_new%SWdown = 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) - met_force_new%RefH = force_array(:,10) - end if - - deallocate(force_array) - deallocate(GEOSgcm_name) - - end subroutine get_GEOSs2s_v3 - - ! ************************************************************************* subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & @@ -3516,6 +3186,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & character(40), dimension(N_M21C_vars, N_defs_cols) :: M21CCOR_defs character(40), dimension(N_M21C_vars, N_defs_cols) :: M21CCSINT_defs character(40), dimension(N_M21C_vars, N_defs_cols) :: M21CCSCOR_defs + character(40), dimension(N_G5DAS_vars, N_defs_cols) :: S2S3FCST_defs + character(40), dimension(N_G5DAS_vars, N_defs_cols) :: S2S3AODAS_defs character(40), dimension(:,:), allocatable :: GEOSgcm_defs @@ -3548,8 +3220,10 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & logical :: minimize_shift, use_prec_corr, use_bkg, tmp_init - logical :: daily_met_files, daily_precipcorr_files - + logical :: daily_met_files, daily_precipcorr_files + + logical :: is_S2S3_fcst + integer :: nv_id, ierr, icount(3), istart(3), lonid, latid character(len=*), parameter :: Iam = 'get_GEOS' @@ -3586,6 +3260,33 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & G5DAS_defs(11,:)=[character(len=40):: 'QLML ','inst','inst1_2d_lfo_Nx','diag','S'] G5DAS_defs(12,:)=[character(len=40):: 'SPEEDLML','inst','inst1_2d_lfo_Nx','diag','S'] + ! ----------------------------------------------------------------------- + ! + ! define GEOS5 S2Sv3 FCST specs + ! + + S2S3FCST_defs( 1,:)=[character(len=40):: 'SWGDN ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 2,:)=[character(len=40):: 'LWGAB ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 3,:)=[character(len=40):: 'PARDR ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 4,:)=[character(len=40):: 'PARDF ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 5,:)=[character(len=40):: 'PCU ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 6,:)=[character(len=40):: 'PLS ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 7,:)=[character(len=40):: 'SNO ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 8,:)=[character(len=40):: 'PS ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 9,:)=[character(len=40):: 'HLML ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs(10,:)=[character(len=40):: 'TA ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs(11,:)=[character(len=40):: 'QA ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs(12,:)=[character(len=40):: 'SPEED ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + + ! ----------------------------------------------------------------------- + ! + ! define GEOS5 S2Sv3 AODAS specs + ! + + S2S3AODAS_defs = S2S3FCST_defs + S2S3AODAS_defs( 5,:)=[character(len=40):: 'PCUCORR','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3AODAS_defs( 6,:)=[character(len=40):: 'PLSCORR','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3AODAS_defs( 7,:)=[character(len=40):: 'SNOCORR','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] ! ----------------------------------------------------------------------- ! @@ -3988,8 +3689,10 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & met_file_ext = 'nc4' daily_met_files = .false. - + precip_corr_file_ext = 'nc4' + + is_S2S3_fcst = .false. if (met_tag(4:8)=='merra') then ! MERRA @@ -4104,7 +3807,39 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & call parse_MERRA2_met_tag( met_path, met_tag, date_time_bkwd, & met_path_bkwd, prec_path_bkwd, met_tag_bkwd, use_prec_corr ) - + + elseif (met_tag(1:8)=='GEOSS2S3') then ! GEOS S2S v3 + + N_GEOSgcm_vars = N_G5DAS_vars + + PAR_available = .false. ! S2Sv3 doesn't hae PARxx + + single_time_in_file = .false. ! FCST are in monthly,AODAS in daily + + if (met_tag(9:12)=='FCST') then + + is_S2S3_fcst = .true. + + GEOSgcm_defs = S2S3FCST_defs + + elseif (met_tag(9:13)=='AODAS') then + + daily_met_files = .true. + + GEOSgcm_defs = S2S3AODAS_defs + + else + + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown "GEOSS2S3[xxx]" met_tag') + + end if + + met_path_fwd = met_path + met_tag_fwd = met_tag + + met_path_bkwd = met_path + met_tag_bkwd = met_tag + else ! GEOS ADAS (FP, GEOSIT) call parse_G5DAS_met_tag( met_path, met_tag, date_time_inst, & @@ -4169,6 +3904,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! if (init==.true.) make second attempt (j=2) to allow for possibly ! missing "diag_sfc" or "tavg" file at date_time_bkwd (and try reading ! the file at date_time_fwd). + if (PAR_available .or. (GEOSgcm_var<3 .or. GEOSgcm_var>4)) then ! skip PARDR/DF for S2Sv3 do j=1,2 @@ -4221,31 +3957,36 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & single_time_in_file = .not. daily_met_files ! MERRA-2 files are daily files end if - - if ( file_exists) then - - exit ! exit j loop after successfully finding file + + + if (.not.(is_S2S3_fcst .and. tmp_init .and. (j==1)) ) then ! always skip first timestep bkwd for S2Sv3 fcst + + if ( file_exists ) then + + exit ! exit j loop after successfully finding file - elseif ( & - (j==1) .and. & - (tmp_init) .and. & - (trim(GEOSgcm_defs(GEOSgcm_var,2))=='tavg') .and. & - (root_logit) ) then + elseif ( & + (j==1) .and. & + (tmp_init) .and. & + (trim(GEOSgcm_defs(GEOSgcm_var,2))=='tavg') .and. & + (root_logit) ) then - if (.not. MERRA_file_specs) write (logunit,'(400A)') & - 'NOTE: Initialization. Data from tavg file are not used ' // & - 'with lfo inst/tavg forcing, but dummy values must be ' // & - 'read from some file for backward compatibility with ' // & - 'MERRA forcing.' + if (.not. MERRA_file_specs) write (logunit,'(400A)') & + 'NOTE: Initialization. Data from tavg file are not used ' // & + 'with lfo inst/tavg forcing, but dummy values must be ' // & + 'read from some file for backward compatibility with ' // & + 'MERRA forcing.' - write (logunit,*) 'try again with different file...' + write (logunit,*) 'try again with different file...' - else + else - call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error finding met forcing file') + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error finding met forcing file') - end if - + end if ! if(file_exists) + + end if ! if(.not.(is_S2S3_fcst xxx) + end do ! j=1,2 ! open file, extract coord info, prep horizontal interpolation info (if not done already) @@ -4292,9 +4033,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! ---------------------------------------------- ! ! read global gridded field of given variable - call LDAS_GetVar( fid, trim(GEOSgcm_defs(GEOSgcm_var,1)), & - YYYYMMDD, HHMMSS, single_time_in_file, local_info, ptrShForce, rc) + YYYYMMDD, HHMMSS, single_time_in_file, is_S2S3_fcst, local_info, ptrShForce, rc) if (rc<0) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error reading gfio file') endif @@ -4382,9 +4122,9 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & HHMMSS = date_time_tmp%hour*10000+date_time_tmp%min*100 +date_time_tmp%sec ! read global gridded field of given variable - + call LDAS_GetVar( fid, trim(GEOSgcm_defs(GEOSgcm_var,1)), & - YYYYMMDD, HHMMSS, .false., local_info, ptrShForce, rc) + YYYYMMDD, HHMMSS, .false.,.false., local_info, ptrShForce, rc) if (rc<0) then err_msg = 'error reading gfio file' @@ -4443,7 +4183,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & end if ! if (fid>0) end if ! if (minimize_shift) .and. [...] - + end if ! if (PAR_avaialble .or. [...] end do ! do GEOSgcm_var = 1,N_GEOSgcm_vars call FileOpenedHash%free( GEOS_closefile,.false. ) @@ -4483,9 +4223,12 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & met_force_new%SWdown = force_array(:, 1) met_force_new%LWdown = force_array(:, 2) - met_force_new%PARdrct = force_array(:, 3) - met_force_new%PARdffs = force_array(:, 4) - + + if (PAR_available) then + met_force_new%PARdrct = force_array(:, 3) + met_force_new%PARdffs = force_array(:, 4) + end if + met_force_new%Psurf = force_array(:, 8) met_force_new%RefH = force_array(:, 9) @@ -4635,7 +4378,7 @@ end subroutine get_GEOS ! ****************************************************************** - subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, local_info, & + subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, is_S2S3_fcst, local_info, & ptrShForce, rc) ! get LDAS forcing variable @@ -4649,6 +4392,7 @@ subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, local_ integer, intent(in) :: yyyymmdd ! Year-month-day, e.g., 19971003 integer, intent(in) :: hhmmss ! Hour-minute-second, e.g., 120000 logical, intent(in) :: single_time_in_file ! if true, no time index is necessary + logical, intent(in) :: is_S2S3_fcst ! S2S3 fcst data are in monthly file type(local_grid), intent(in) :: local_info !OUTPUT PARAMETERS: real, pointer, intent(inout) :: ptrShForce(:,:) ! Gridded data read for this time @@ -4662,6 +4406,8 @@ subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, local_ type(c_ptr) :: c_address integer :: nv_id,imin, jmin, imax, jmax,ierr integer :: DiffDate + integer :: dt_hours, day,month,year, hour + integer :: tdimid, timeFull, timeLen integer :: status character(*), parameter :: Iam="LDAS_getvar" logical :: isCubed ! forcing on cs grid: true/false @@ -4682,7 +4428,26 @@ subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, local_ iicount(3) = 1 endif - if (.not. single_time_in_file ) then ! determine start index + if (is_S2S3_fcst) then + dt_hours = 3 ! forcing time interval + day = mod (yyyymmdd,100) + month = mod (yyyymmdd - day, 10000)/100 + year = (yyyymmdd - month*100 -day) / 10000 + hour = hhmmss / 10000 + if ( MOD(hour-1,dt_hours) .eq. 0) then + ! S2S3 monthly forecast files store a limited temporal range of records, + ! beginning at the forecast initial date minus 1.5 hours and ending at month's end. + ! The following information combined with the time dimension size calculates the + ! record index corresponding to the current timestamp within this file structure. + timeIndex = int(((day-1)*24 + hour) /dt_hours) + 1 ! the index of current time in a month + timeFull = days_in_month(year, month) * 24 / 3 ! number of time record for the full month + else + print *, 'LDAS_GetVar: S2Sv3 forcing require times fall on 1:30,4:30,...' + rc = -6 + return + end if + + elseif (.not. single_time_in_file ) then ! determine start index call GetBegDateTime ( fid, begDate, begTime, incSecs, rc ) if (rc .NE. 0) then print *, 'LDAS_GetVar: could not determine begin_date/begin_time' @@ -4727,6 +4492,14 @@ subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, local_ call c_f_pointer(c_address,tmpShared,shape=icount) rc= NF90_GET_VAR( fid, nv_id, tmpShared, start=istart,count=icount) else + if (is_S2S3_fcst) then + rc= nf90_inq_dimid(fid,"time",tdimid) + _ASSERT( rc == nf90_noerr, "nf90 error") + rc= nf90_Inquire_Dimension(fid, tdimid, len=timeLen) + _ASSERT( rc == nf90_noerr, "nf90 error") + iistart(3) = timeIndex - (timeFull - timeLen) + end if + rc= NF90_GET_VAR( fid, nv_id, ptrShForce, start=iistart,count=iicount) endif _ASSERT( rc == nf90_noerr, "nf90 error") @@ -5870,6 +5643,8 @@ subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_fi character( 16) :: time_stamp character( 4) :: YYYY, HHMM, day_dir character( 2) :: MM, DD + character( 8) :: init_YYYYMMDD ! S2Sv3 initial date, e.g. "20160101" + character( 4) :: ens_num ! S2Sv3 ensemble member, e.g. "ens1" integer :: tmpind, tmpindend @@ -5894,7 +5669,7 @@ subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_fi if (daily_file) then time_stamp(1:8) = YYYY // MM // DD - + elseif ( & index(met_tag,'GEOSIT') > 0 .or. index(met_tag,'geosit') > 0 .or. & index(met_tag,'M21C' ) > 0 .or. index(met_tag,'m21c' ) > 0 & @@ -5929,12 +5704,25 @@ subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_fi day_dir = 'D' // DD // '/' + elseif (met_tag(1:12) == 'GEOSS2S3FCST') then + ! parse met_tag + + ! 1234567890123456789012345678 + ! GEOSS2S3FCST__ensX__YYYYMMDD + ! GEOSS2S3AODAS + + ens_num = trim(met_tag(15:18)) + init_YYYYMMDD = trim(met_tag(21:28)) + + fname = init_YYYYMMDD // '/' // ens_num // '/GEOSS2S3.' // YYYY // MM // '.nc4' + else ! GEOS FP with experiment-specific file names and MERRA-2, e.g., ! ! f525_p5_fp.inst1_2d_lfo_Nx.20200507_0000z.nc4 ! MERRA2_400.inst1_2d_lfo_Nx.20200507.nc4 + ! fname = trim(met_tag) // '.' // trim(GEOSgcm_defs(3)) // '.' // & trim(time_stamp) // '.' // trim(file_ext) From c7d65d3e3a89bf6bfb8f09a69a0c74b78aa97427 Mon Sep 17 00:00:00 2001 From: Rolf Reichle <54944691+gmao-rreichle@users.noreply.github.com> Date: Fri, 15 Aug 2025 11:46:30 -0400 Subject: [PATCH 05/16] Update CHANGELOG.md --- CHANGELOG.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3f1a7245..aaab0b59 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,11 +15,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed -- Cleaned up ldas_setup. Split it to ldas.py and setup_utils.py, +- Cleaned up ldas_setup. Split out ldas.py and setup_utils.py. +- Added reader for surface meteorological forcing from S2S-3. ### Fixed -- Fixed Restart = 1 when the domain is not global +- Fixed Restart=1 when the domain is not global. ### Removed From 6e28d99114ef93460648389691db720768205ebc Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 20 Aug 2025 13:04:50 -0400 Subject: [PATCH 06/16] revised implementation of S2S3 forcing reader; added documentation (LDAS_Forcing.F90, GEOS_MetforceGridComp.F90) --- .../GEOS_MetforceGridComp.F90 | 10 +- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 294 +++++++++--------- 2 files changed, 145 insertions(+), 159 deletions(-) diff --git a/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 b/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 index 308735c7..2dc99129 100644 --- a/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 +++ b/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90 @@ -608,7 +608,7 @@ subroutine Initialize(gc, import, export, clock, rc) integer :: local_nt, k, NUM_ENSEMBLE, i1, i2, j1, j2 integer :: ForceDtStep type(met_force_type) :: mf_nodata - logical :: MERRA_file_specs, ensemble_forcing + logical :: MERRA_file_specs, S2S3_file_specs, ensemble_forcing logical :: backward_looking_fluxes real, pointer :: TileLats(:) real, pointer :: TileLons(:) @@ -768,6 +768,7 @@ subroutine Initialize(gc, import, export, clock, rc) internal%mf%hinterp, & AEROSOL_DEPOSITION, & MERRA_file_specs, & + S2S3_file_specs, & backward_looking_fluxes, & internal%mf%DataNxt, & .true. & ! init @@ -776,7 +777,7 @@ subroutine Initialize(gc, import, export, clock, rc) if (backward_looking_fluxes) & call LDAS_move_new_force_to_old( & - MERRA_file_specs, AEROSOL_DEPOSITION, & + MERRA_file_specs, S2S3_file_specs, AEROSOL_DEPOSITION, & internal%mf%DataNxt, internal%mf%DataPrv ) ! Turn timer off @@ -851,7 +852,7 @@ subroutine Run(gc, import, export, clock, rc) type(met_force_type), pointer, contiguous :: DataTmp(:)=>null() type(met_force_type) :: mf_nodata - logical :: MERRA_file_specs + logical :: MERRA_file_specs, S2S3_file_specs logical :: backward_looking_fluxes integer :: AEROSOL_DEPOSITION ! Export pointers @@ -995,6 +996,7 @@ subroutine Run(gc, import, export, clock, rc) internal%mf%hinterp, & AEROSOL_DEPOSITION, & MERRA_file_specs, & + S2S3_file_specs, & backward_looking_fluxes, & internal%mf%DataNxt, & .false. & ! init @@ -1003,7 +1005,7 @@ subroutine Run(gc, import, export, clock, rc) if (backward_looking_fluxes) & call LDAS_move_new_force_to_old( & - MERRA_file_specs, AEROSOL_DEPOSITION, & + MERRA_file_specs, S2S3_file_specs, AEROSOL_DEPOSITION, & internal%mf%DataNxt, internal%mf%DataPrv ) ! -compute-average-zenith-angle-over-daylight-part-of-forcing-interval- diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 56a06588..71e48621 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -38,8 +38,7 @@ module LDAS_ForceMod datetime_lt_refdatetime, & datetime_le_refdatetime, & is_leap_year, & - get_dofyr_pentad, & - days_in_month + get_dofyr_pentad use LDAS_ExceptionsMod, ONLY: & ldas_abort, & @@ -103,7 +102,8 @@ module LDAS_ForceMod subroutine get_forcing( date_time, force_dtstep, met_path, met_tag, & N_catd, tile_coord, MET_HINTERP, AEROSOL_DEPOSITION, & - MERRA_file_specs, bkwd_looking_fluxes, met_force_obs_tile_new, & + MERRA_file_specs, S2S3_file_specs, & + bkwd_looking_fluxes, met_force_obs_tile_new, & init ) ! Read and check meteorological forcing data for the domain. @@ -166,7 +166,7 @@ subroutine get_forcing( date_time, force_dtstep, met_path, met_tag, & ! intent out: - logical, intent(out) :: MERRA_file_specs + logical, intent(out) :: MERRA_file_specs, S2S3_file_specs logical, intent(out) :: bkwd_looking_fluxes type(met_force_type), dimension(N_catd), intent(out) :: & @@ -222,6 +222,7 @@ subroutine get_forcing( date_time, force_dtstep, met_path, met_tag, & ! initialize MERRA_file_specs = .false. + S2S3_file_specs = .false. bkwd_looking_fluxes = .false. @@ -334,7 +335,8 @@ subroutine get_forcing( date_time, force_dtstep, met_path, met_tag, & N_catd, tile_coord, MET_HINTERP, AEROSOL_DEPOSITION, & supported_option_MET_HINTERP, & supported_option_AEROSOL_DEPOSITION, & - met_force_obs_tile_new, nodata_forcing, PAR_available, MERRA_file_specs, & + met_force_obs_tile_new, & + nodata_forcing, PAR_available, MERRA_file_specs, S2S3_file_specs, & init ) ! subroutine get_GEOS() provided backward-looking fluxes. @@ -432,7 +434,7 @@ end subroutine get_forcing !************************************************************************************** - subroutine LDAS_move_new_force_to_old( MERRA_file_specs, AEROSOL_DEPOSITION, & + subroutine LDAS_move_new_force_to_old( MERRA_file_specs, S2S3_file_specs, AEROSOL_DEPOSITION, & new_force, old_force ) ! move *flux*-type forcing data from "new" to "old"; @@ -441,7 +443,7 @@ subroutine LDAS_move_new_force_to_old( MERRA_file_specs, AEROSOL_DEPOSITION, & implicit none - logical, intent(in) :: MERRA_file_specs + logical, intent(in) :: MERRA_file_specs, S2S3_file_specs integer, intent(in) :: AEROSOL_DEPOSITION type(met_force_type), dimension(:), intent(inout) :: new_force @@ -465,7 +467,7 @@ subroutine LDAS_move_new_force_to_old( MERRA_file_specs, AEROSOL_DEPOSITION, & ! [moved here from below, reichle, 28 Jan 2021] ! treat Wind as flux when forcing with MERRA - if (MERRA_file_specs) then + if (MERRA_file_specs .or. S2S3_file_specs) then old_force%Wind = new_force%Wind new_force%Wind = nodata_generic endif @@ -3055,7 +3057,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & N_catd, tile_coord, MET_HINTERP, AEROSOL_DEPOSITION, & supported_option_MET_HINTERP, & supported_option_AEROSOL_DEPOSITION, & - met_force_new, nodata_forcing, PAR_available, MERRA_file_specs, & + met_force_new, & + nodata_forcing, PAR_available, MERRA_file_specs, S2S3_file_specs, & init ) ! reichle, 5 March 2008 - adapted from get_GEOSgcm_gfio to work with DAS @@ -3155,7 +3158,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & logical, intent(out) :: PAR_available logical, intent(out) :: MERRA_file_specs ! original MERRA specs, not MERRA-2 - + logical, intent(out) :: S2S3_file_specs + ! optional: logical, intent(in), optional :: init @@ -3168,7 +3172,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & integer, parameter :: N_MERRA_vars = 13 integer, parameter :: N_MERRA2_vars = 12 ! same as for G5DAS (excl Aerosol vars) integer, parameter :: N_Aerosol_vars = 60 ! additional aerosol forcing vars for GOSWIM (w/ MERRA-2 only for now) - integer, parameter :: N_M21C_vars = 12 + integer, parameter :: N_M21C_vars = 12 + integer, parameter :: N_S2S3_vars = 12 integer, parameter :: N_MERRA2plusAerosol_vars = N_MERRA2_vars + N_Aerosol_vars @@ -3186,8 +3191,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & character(40), dimension(N_M21C_vars, N_defs_cols) :: M21CCOR_defs character(40), dimension(N_M21C_vars, N_defs_cols) :: M21CCSINT_defs character(40), dimension(N_M21C_vars, N_defs_cols) :: M21CCSCOR_defs - character(40), dimension(N_G5DAS_vars, N_defs_cols) :: S2S3FCST_defs - character(40), dimension(N_G5DAS_vars, N_defs_cols) :: S2S3AODAS_defs + character(40), dimension(N_S2S3_vars, N_defs_cols) :: S2S3FCST_defs + character(40), dimension(N_S2S3_vars, N_defs_cols) :: S2S3AODAS_defs character(40), dimension(:,:), allocatable :: GEOSgcm_defs @@ -3260,34 +3265,6 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & G5DAS_defs(11,:)=[character(len=40):: 'QLML ','inst','inst1_2d_lfo_Nx','diag','S'] G5DAS_defs(12,:)=[character(len=40):: 'SPEEDLML','inst','inst1_2d_lfo_Nx','diag','S'] - ! ----------------------------------------------------------------------- - ! - ! define GEOS5 S2Sv3 FCST specs - ! - - S2S3FCST_defs( 1,:)=[character(len=40):: 'SWGDN ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 2,:)=[character(len=40):: 'LWGAB ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 3,:)=[character(len=40):: 'PARDR ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 4,:)=[character(len=40):: 'PARDF ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 5,:)=[character(len=40):: 'PCU ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 6,:)=[character(len=40):: 'PLS ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 7,:)=[character(len=40):: 'SNO ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 8,:)=[character(len=40):: 'PS ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 9,:)=[character(len=40):: 'HLML ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs(10,:)=[character(len=40):: 'TA ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs(11,:)=[character(len=40):: 'QA ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs(12,:)=[character(len=40):: 'SPEED ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - - ! ----------------------------------------------------------------------- - ! - ! define GEOS5 S2Sv3 AODAS specs - ! - - S2S3AODAS_defs = S2S3FCST_defs - S2S3AODAS_defs( 5,:)=[character(len=40):: 'PCUCORR','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3AODAS_defs( 6,:)=[character(len=40):: 'PLSCORR','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3AODAS_defs( 7,:)=[character(len=40):: 'SNOCORR','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - ! ----------------------------------------------------------------------- ! ! define coupled land/atm DAS file specs (i.e., use bkg.lfo_* files) @@ -3625,23 +3602,54 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! - use "lfo" files ! reichle, 1 Dec 2009 - ! MERRA - ! collection - - MERRA_defs( 1,:)=[character(len=40):: 'SWGDN ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "rad" - MERRA_defs( 2,:)=[character(len=40):: 'LWGAB ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "rad" - MERRA_defs( 3,:)=[character(len=40):: 'PARDR ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "lnd" - MERRA_defs( 4,:)=[character(len=40):: 'PARDF ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "lnd" - MERRA_defs( 5,:)=[character(len=40):: 'PRECTOT','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "lnd" - MERRA_defs( 6,:)=[character(len=40):: 'PRECCON','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "flx" - MERRA_defs( 7,:)=[character(len=40):: 'PRECSNO','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "lnd" - MERRA_defs( 8,:)=[character(len=40):: 'PS ','tavg','tavg1_2d_lfo_Nx','diag','S'] ! "slv" - MERRA_defs( 9,:)=[character(len=40):: 'HLML ','tavg','tavg1_2d_lfo_Nx','diag','S'] ! "flx" - MERRA_defs(10,:)=[character(len=40):: 'TLML ','tavg','tavg1_2d_lfo_Nx','diag','S'] ! "flx" - MERRA_defs(11,:)=[character(len=40):: 'QLML ','tavg','tavg1_2d_lfo_Nx','diag','S'] ! "flx" - MERRA_defs(12,:)=[character(len=40):: 'ULML ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "flx" - MERRA_defs(13,:)=[character(len=40):: 'VLML ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "flx" + ! MERRA + ! collection + + MERRA_defs( 1,:)=[character(len=40):: 'SWGDN ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "rad" + MERRA_defs( 2,:)=[character(len=40):: 'LWGAB ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "rad" + MERRA_defs( 3,:)=[character(len=40):: 'PARDR ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "lnd" + MERRA_defs( 4,:)=[character(len=40):: 'PARDF ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "lnd" + MERRA_defs( 5,:)=[character(len=40):: 'PRECTOT','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "lnd" + MERRA_defs( 6,:)=[character(len=40):: 'PRECCON','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "flx" + MERRA_defs( 7,:)=[character(len=40):: 'PRECSNO','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "lnd" + MERRA_defs( 8,:)=[character(len=40):: 'PS ','tavg','tavg1_2d_lfo_Nx','diag','S'] ! "slv" + MERRA_defs( 9,:)=[character(len=40):: 'HLML ','tavg','tavg1_2d_lfo_Nx','diag','S'] ! "flx" + MERRA_defs(10,:)=[character(len=40):: 'TLML ','tavg','tavg1_2d_lfo_Nx','diag','S'] ! "flx" + MERRA_defs(11,:)=[character(len=40):: 'QLML ','tavg','tavg1_2d_lfo_Nx','diag','S'] ! "flx" + MERRA_defs(12,:)=[character(len=40):: 'ULML ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "flx" + MERRA_defs(13,:)=[character(len=40):: 'VLML ','tavg','tavg1_2d_lfo_Nx','diag','F'] ! "flx" + + ! ----------------------------------------------------------------------- + ! + ! define GEOS S2S3 FCST specs + ! + ! use *only* 3-hourly "tavg" files b/c instantaneous output is not available + ! + + S2S3FCST_defs( 1,:)=[character(len=40):: 'SWGDN ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 2,:)=[character(len=40):: 'LWGAB ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 3,:)=[character(len=40):: 'dummy ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] ! no PARDR for S2S3 + S2S3FCST_defs( 4,:)=[character(len=40):: 'dummy ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] ! no PARDF for S2S3 + S2S3FCST_defs( 5,:)=[character(len=40):: 'PCU ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 6,:)=[character(len=40):: 'PLS ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 7,:)=[character(len=40):: 'SNO ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 8,:)=[character(len=40):: 'PS ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift + S2S3FCST_defs( 9,:)=[character(len=40):: 'HLML ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift + S2S3FCST_defs( 10,:)=[character(len=40):: 'TA ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift + S2S3FCST_defs( 11,:)=[character(len=40):: 'QA ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift + S2S3FCST_defs( 12,:)=[character(len=40):: 'SPEED ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + + ! ----------------------------------------------------------------------- + ! + ! define GEOS S2S3 AODAS specs; same as FCST except for corrected precip + + S2S3AODAS_defs = S2S3FCST_defs + + S2S3AODAS_defs( 5,1)=[character(len=40):: 'PCUCORR'] + S2S3AODAS_defs( 6,1)=[character(len=40):: 'PLSCORR'] + S2S3AODAS_defs( 7,1)=[character(len=40):: 'SNOCORR'] + ! -------------------------------------------------------------------- ! ! preparations @@ -3659,7 +3667,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & tol = abs(nodata_forcing*nodata_tolfrac_generic) - ! all GEOS forcing datasets provide PAR (so far) + ! most GEOS forcing datasets provide PAR PAR_available = .true. @@ -3685,6 +3693,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! initialize to most likely values, overwrite below as needed MERRA_file_specs = .false. + S2S3_file_specs = .false. met_file_ext = 'nc4' @@ -3810,21 +3819,23 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & elseif (met_tag(1:8)=='GEOSS2S3') then ! GEOS S2S v3 - N_GEOSgcm_vars = N_G5DAS_vars + N_GEOSgcm_vars = N_S2S3_vars + + PAR_available = .false. ! S2S3 does not have PAR - PAR_available = .false. ! S2Sv3 doesn't hae PARxx + S2S3_file_specs = .true. - single_time_in_file = .false. ! FCST are in monthly,AODAS in daily + single_time_in_file = .false. ! FCST: monthly files, AODAS: daily files - if (met_tag(9:12)=='FCST') then + if (met_tag(9:12)=='FCST' ) then - is_S2S3_fcst = .true. + is_S2S3_fcst = .true. - GEOSgcm_defs = S2S3FCST_defs + GEOSgcm_defs = S2S3FCST_defs elseif (met_tag(9:13)=='AODAS') then - daily_met_files = .true. + daily_met_files = .true. GEOSgcm_defs = S2S3AODAS_defs @@ -3834,8 +3845,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & end if - met_path_fwd = met_path - met_tag_fwd = met_tag + met_path_fwd = met_path + met_tag_fwd = met_tag met_path_bkwd = met_path met_tag_bkwd = met_tag @@ -3891,6 +3902,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & do GEOSgcm_var = 1,N_GEOSgcm_vars + if (GEOSgcm_defs(GEOSgcm_var,1)=="dummy") cycle ! skip "dummy" variable (e.g., no PAR for S2S3) + ! open GEOS file (G5DAS or MERRA or MERRA-2) ! ! Initial "tavg1_2d_*_Nx" files may not be available. In this case, @@ -3904,7 +3917,6 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! if (init==.true.) make second attempt (j=2) to allow for possibly ! missing "diag_sfc" or "tavg" file at date_time_bkwd (and try reading ! the file at date_time_fwd). - if (PAR_available .or. (GEOSgcm_var<3 .or. GEOSgcm_var>4)) then ! skip PARDR/DF for S2Sv3 do j=1,2 @@ -3958,34 +3970,31 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & end if - - if (.not.(is_S2S3_fcst .and. tmp_init .and. (j==1)) ) then ! always skip first timestep bkwd for S2Sv3 fcst - if ( file_exists ) then - - exit ! exit j loop after successfully finding file + if ( file_exists .and. (.not. is_S2S3_fcst) ) then ! S2S3 FCST has monthly files - elseif ( & - (j==1) .and. & - (tmp_init) .and. & - (trim(GEOSgcm_defs(GEOSgcm_var,2))=='tavg') .and. & - (root_logit) ) then + exit ! exit j loop after successfully finding file - if (.not. MERRA_file_specs) write (logunit,'(400A)') & - 'NOTE: Initialization. Data from tavg file are not used ' // & - 'with lfo inst/tavg forcing, but dummy values must be ' // & - 'read from some file for backward compatibility with ' // & - 'MERRA forcing.' + elseif ( & + (j==1) .and. & + (tmp_init) .and. & + (trim(GEOSgcm_defs(GEOSgcm_var,2))=='tavg') .and. & + (root_logit) ) then - write (logunit,*) 'try again with different file...' + if ( .not. (MERRA_file_specs .or. S2S3_file_specs) ) & + write (logunit,'(400A)') & + 'NOTE: Initialization. Data from tavg file are not used ' // & + 'with lfo inst/tavg forcing, but dummy values must be ' // & + 'read from some file for backward compatibility with ' // & + 'MERRA forcing.' - else + write (logunit,*) 'try again with different file...' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error finding met forcing file') + else - end if ! if(file_exists) - - end if ! if(.not.(is_S2S3_fcst xxx) + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error finding met forcing file') + + end if end do ! j=1,2 @@ -4034,7 +4043,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! ! read global gridded field of given variable call LDAS_GetVar( fid, trim(GEOSgcm_defs(GEOSgcm_var,1)), & - YYYYMMDD, HHMMSS, single_time_in_file, is_S2S3_fcst, local_info, ptrShForce, rc) + YYYYMMDD, HHMMSS, single_time_in_file, local_info, ptrShForce, rc) if (rc<0) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'error reading gfio file') endif @@ -4124,7 +4133,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! read global gridded field of given variable call LDAS_GetVar( fid, trim(GEOSgcm_defs(GEOSgcm_var,1)), & - YYYYMMDD, HHMMSS, .false.,.false., local_info, ptrShForce, rc) + YYYYMMDD, HHMMSS, .false., local_info, ptrShForce, rc) if (rc<0) then err_msg = 'error reading gfio file' @@ -4183,7 +4192,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & end if ! if (fid>0) end if ! if (minimize_shift) .and. [...] - end if ! if (PAR_avaialble .or. [...] + end do ! do GEOSgcm_var = 1,N_GEOSgcm_vars call FileOpenedHash%free( GEOS_closefile,.false. ) @@ -4197,29 +4206,31 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! from GEOSgcm files: ! - ! G5DAS - ! M2INT MERRA - ! M2COR - ! - ! force_array(:, 1) = SWGDN SWGDN W/m2 (downward shortwave) - ! force_array(:, 2) = LWGAB LWGAB W/m2 ("absorbed" longwave) - ! force_array(:, 3) = PARDR PARDR W/m2 (direct PAR) - ! force_array(:, 4) = PARDF PARDF W/m2 (diffuse PAR) - ! force_array(:, 5) = PRECCU[*] PRECTOT kg/m2/s (*see below*) - ! force_array(:, 6) = PRECLS[*] PRECCON kg/m2/s (*see below*) - ! force_array(:, 7) = PRECSN[*] PRECSNO kg/m2/s (*see below*) - ! force_array(:, 8) = PS PS Pa (surface air pressure) - ! force_array(:, 9) = HLML HLML m (height of lowest model level "LML") - ! force_array(:,10) = TLML TLML K (air temperature at LML) - ! force_array(:,11) = QLML QLML kg/kg (air spec humidity at LML) - ! force_array(:,12) = SPEEDLML ULML m/s (wind speed/U-wind at LML) - ! force_array(:,13) = n/a VLML m/s ( V-wind at LML) + ! G5DAS + ! G5BKG + ! GEOSIT + ! M2* + ! M21C* MERRA S2S3* + ! + ! force_array(:, 1) = SWGDN SWGDN SWGDN W/m2 (downward shortwave) + ! force_array(:, 2) = LWGAB LWGAB LWGAB W/m2 ("absorbed" longwave) + ! force_array(:, 3) = PARDR PARDR n/a W/m2 (direct PAR) + ! force_array(:, 4) = PARDF PARDF n/a W/m2 (diffuse PAR) + ! force_array(:, 5) = PRECCU[*] PRECTOT PCU[*] kg/m2/s (*see below*) + ! force_array(:, 6) = PRECLS[*] PRECCON PLS[*] kg/m2/s (*see below*) + ! force_array(:, 7) = PRECSN[*] PRECSNO SNO[*] kg/m2/s (*see below*) + ! force_array(:, 8) = PS PS PS Pa (surface air pressure) + ! force_array(:, 9) = HLML HLML HLML m (height of lowest model level "LML") + ! force_array(:,10) = TLML TLML TA K (air temperature at LML) + ! force_array(:,11) = QLML QLML QA kg/kg (air spec humidity at LML) + ! force_array(:,12) = SPEEDLML ULML SPEED m/s (wind speed/U-wind at LML) + ! force_array(:,13) = n/a VLML n/a 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) + ! PRECCU, PCU kg/m2/s (convective rain) + ! PRECLS, PLS kg/m2/s (large-scale rain) + ! PRECSNO, SNO kg/m2/s (total snow) met_force_new%SWdown = force_array(:, 1) met_force_new%LWdown = force_array(:, 2) @@ -4253,7 +4264,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & end if - else ! G5DAS file specs + else ! other file specs met_force_new(k)%Wind = force_array(k,12) @@ -4274,6 +4285,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & if (MERRA_file_specs) then + ! deal with MERRA precip components + if (force_array(k,5)>0) then met_force_new(k)%Snowf = force_array(k,7) @@ -4297,7 +4310,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & else - ! G5DAS file specs + ! other file specs met_force_new(k)%Rainf = force_array(k,5)+force_array(k,6) met_force_new(k)%Rainf_C = force_array(k,5) @@ -4378,7 +4391,7 @@ end subroutine get_GEOS ! ****************************************************************** - subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, is_S2S3_fcst, local_info, & + subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, local_info, & ptrShForce, rc) ! get LDAS forcing variable @@ -4392,7 +4405,6 @@ subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, is_S2S integer, intent(in) :: yyyymmdd ! Year-month-day, e.g., 19971003 integer, intent(in) :: hhmmss ! Hour-minute-second, e.g., 120000 logical, intent(in) :: single_time_in_file ! if true, no time index is necessary - logical, intent(in) :: is_S2S3_fcst ! S2S3 fcst data are in monthly file type(local_grid), intent(in) :: local_info !OUTPUT PARAMETERS: real, pointer, intent(inout) :: ptrShForce(:,:) ! Gridded data read for this time @@ -4406,8 +4418,6 @@ subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, is_S2S type(c_ptr) :: c_address integer :: nv_id,imin, jmin, imax, jmax,ierr integer :: DiffDate - integer :: dt_hours, day,month,year, hour - integer :: tdimid, timeFull, timeLen integer :: status character(*), parameter :: Iam="LDAS_getvar" logical :: isCubed ! forcing on cs grid: true/false @@ -4428,26 +4438,7 @@ subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, is_S2S iicount(3) = 1 endif - if (is_S2S3_fcst) then - dt_hours = 3 ! forcing time interval - day = mod (yyyymmdd,100) - month = mod (yyyymmdd - day, 10000)/100 - year = (yyyymmdd - month*100 -day) / 10000 - hour = hhmmss / 10000 - if ( MOD(hour-1,dt_hours) .eq. 0) then - ! S2S3 monthly forecast files store a limited temporal range of records, - ! beginning at the forecast initial date minus 1.5 hours and ending at month's end. - ! The following information combined with the time dimension size calculates the - ! record index corresponding to the current timestamp within this file structure. - timeIndex = int(((day-1)*24 + hour) /dt_hours) + 1 ! the index of current time in a month - timeFull = days_in_month(year, month) * 24 / 3 ! number of time record for the full month - else - print *, 'LDAS_GetVar: S2Sv3 forcing require times fall on 1:30,4:30,...' - rc = -6 - return - end if - - elseif (.not. single_time_in_file ) then ! determine start index + if (.not. single_time_in_file ) then ! determine start index call GetBegDateTime ( fid, begDate, begTime, incSecs, rc ) if (rc .NE. 0) then print *, 'LDAS_GetVar: could not determine begin_date/begin_time' @@ -4481,6 +4472,7 @@ subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, is_S2S iistart(3)=timeIndex istart(4) =timeIndex endif + ! node root read and share call MAPL_SyncSharedMemory(rc=status) @@ -4492,14 +4484,6 @@ subroutine LDAS_GetVar(fid, vname, yyyymmdd, hhmmss, single_time_in_file, is_S2S call c_f_pointer(c_address,tmpShared,shape=icount) rc= NF90_GET_VAR( fid, nv_id, tmpShared, start=istart,count=icount) else - if (is_S2S3_fcst) then - rc= nf90_inq_dimid(fid,"time",tdimid) - _ASSERT( rc == nf90_noerr, "nf90 error") - rc= nf90_Inquire_Dimension(fid, tdimid, len=timeLen) - _ASSERT( rc == nf90_noerr, "nf90 error") - iistart(3) = timeIndex - (timeFull - timeLen) - end if - rc= NF90_GET_VAR( fid, nv_id, ptrShForce, start=iistart,count=iicount) endif _ASSERT( rc == nf90_noerr, "nf90 error") @@ -5643,8 +5627,8 @@ subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_fi character( 16) :: time_stamp character( 4) :: YYYY, HHMM, day_dir character( 2) :: MM, DD - character( 8) :: init_YYYYMMDD ! S2Sv3 initial date, e.g. "20160101" - character( 4) :: ens_num ! S2Sv3 ensemble member, e.g. "ens1" + character( 8) :: S2S3_init_YYYYMMDD ! S2S3 initial date, e.g. "20160101" + character( 4) :: S2S3_ens_num ! S2S3 ensemble member, e.g. "ens1" integer :: tmpind, tmpindend @@ -5705,24 +5689,24 @@ subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_fi day_dir = 'D' // DD // '/' elseif (met_tag(1:12) == 'GEOSS2S3FCST') then + ! parse met_tag ! 1234567890123456789012345678 ! GEOSS2S3FCST__ensX__YYYYMMDD ! GEOSS2S3AODAS - - ens_num = trim(met_tag(15:18)) - init_YYYYMMDD = trim(met_tag(21:28)) - - fname = init_YYYYMMDD // '/' // ens_num // '/GEOSS2S3.' // YYYY // MM // '.nc4' + + S2S3_ens_num = trim(met_tag(15:18)) + S2S3_init_YYYYMMDD = trim(met_tag(21:28)) + + fname = S2S3_init_YYYYMMDD // '/' // S2S3_ens_num // '/GEOSS2S3.' // YYYY // MM // '.nc4' else - ! GEOS FP with experiment-specific file names and MERRA-2, e.g., + ! GEOS FP with experiment-specific file names, MERRA-2, etc, e.g., ! ! f525_p5_fp.inst1_2d_lfo_Nx.20200507_0000z.nc4 ! MERRA2_400.inst1_2d_lfo_Nx.20200507.nc4 - ! fname = trim(met_tag) // '.' // trim(GEOSgcm_defs(3)) // '.' // & trim(time_stamp) // '.' // trim(file_ext) From a25422734f93c29bb01a942af9ff325d203811f0 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 20 Aug 2025 13:23:30 -0400 Subject: [PATCH 07/16] fixed syntax error in previous commit, added documentation (LDAS_Forcing.F90) --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 71e48621..a2d6b7d0 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -324,6 +324,9 @@ subroutine get_forcing( date_time, force_dtstep, met_path, met_tag, & elseif (index(met_tag(1:7), 'GEOSs2s')/=0) then + ! IMPORTANT: met_tag = GEOSs2s* for S2S v2 (note lower case "s2s") + ! met_tag = GEOSS2S3* for S2S v3 (note upper case "S2S") --> handled by get_GEOS() in "else" block + call get_GEOSs2s_v2( date_time_tmp, met_path, met_tag, N_catd, tile_coord, & MET_HINTERP, met_force_obs_tile_new, nodata_forcing, PAR_available) @@ -3644,10 +3647,15 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! define GEOS S2S3 AODAS specs; same as FCST except for corrected precip S2S3AODAS_defs = S2S3FCST_defs - - S2S3AODAS_defs( 5,1)=[character(len=40):: 'PCUCORR'] - S2S3AODAS_defs( 6,1)=[character(len=40):: 'PLSCORR'] - S2S3AODAS_defs( 7,1)=[character(len=40):: 'SNOCORR'] + + ! character(40): + ! 1 2 3 4 + ! 1234567890123456789012345678901234567890 + + + S2S3AODAS_defs( 5,1)=[ 'PCUCORR '] + S2S3AODAS_defs( 6,1)=[ 'PLSCORR '] + S2S3AODAS_defs( 7,1)=[ 'SNOCORR '] ! -------------------------------------------------------------------- From fe7b77109647462688b72f4db6d970dce56e3944 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 20 Aug 2025 13:32:38 -0400 Subject: [PATCH 08/16] another attempt at fixing syntax error from previous commit (LDAS_Forcing.F90) --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 40 +++++++++++++------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index a2d6b7d0..63b763d8 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -3629,18 +3629,18 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! use *only* 3-hourly "tavg" files b/c instantaneous output is not available ! - S2S3FCST_defs( 1,:)=[character(len=40):: 'SWGDN ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 2,:)=[character(len=40):: 'LWGAB ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 3,:)=[character(len=40):: 'dummy ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] ! no PARDR for S2S3 - S2S3FCST_defs( 4,:)=[character(len=40):: 'dummy ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] ! no PARDF for S2S3 - S2S3FCST_defs( 5,:)=[character(len=40):: 'PCU ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 6,:)=[character(len=40):: 'PLS ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 7,:)=[character(len=40):: 'SNO ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] - S2S3FCST_defs( 8,:)=[character(len=40):: 'PS ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift - S2S3FCST_defs( 9,:)=[character(len=40):: 'HLML ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift - S2S3FCST_defs( 10,:)=[character(len=40):: 'TA ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift - S2S3FCST_defs( 11,:)=[character(len=40):: 'QA ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift - S2S3FCST_defs( 12,:)=[character(len=40):: 'SPEED ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 1,:)=[character(len=40):: 'SWGDN ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 2,:)=[character(len=40):: 'LWGAB ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 3,:)=[character(len=40):: 'dummy ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] ! no PARDR for S2S3 + S2S3FCST_defs( 4,:)=[character(len=40):: 'dummy ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] ! no PARDF for S2S3 + S2S3FCST_defs( 5,:)=[character(len=40):: 'PCU ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 6,:)=[character(len=40):: 'PLS ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 7,:)=[character(len=40):: 'SNO ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] + S2S3FCST_defs( 8,:)=[character(len=40):: 'PS ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift + S2S3FCST_defs( 9,:)=[character(len=40):: 'HLML ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift + S2S3FCST_defs(10,:)=[character(len=40):: 'TA ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift + S2S3FCST_defs(11,:)=[character(len=40):: 'QA ','tavg','lfo_tavg_3hr_glo_L720x361','diag','S'] ! note "S" --> minimize_shift + S2S3FCST_defs(12,:)=[character(len=40):: 'SPEED ','tavg','lfo_tavg_3hr_glo_L720x361','diag','F'] ! ----------------------------------------------------------------------- ! @@ -3648,14 +3648,14 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & S2S3AODAS_defs = S2S3FCST_defs - ! character(40): - ! 1 2 3 4 - ! 1234567890123456789012345678901234567890 - - - S2S3AODAS_defs( 5,1)=[ 'PCUCORR '] - S2S3AODAS_defs( 6,1)=[ 'PLSCORR '] - S2S3AODAS_defs( 7,1)=[ 'SNOCORR '] + ! character(40): + ! 1 2 3 4 + ! 1234567890123456789012345678901234567890 + + + S2S3AODAS_defs( 5,1) = 'PCUCORR ' + S2S3AODAS_defs( 6,1) = 'PLSCORR ' + S2S3AODAS_defs( 7,1) = 'SNOCORR ' ! -------------------------------------------------------------------- From 3e7481d38603ca641eab84ab6b55cd282c5bc390 Mon Sep 17 00:00:00 2001 From: Rolf Reichle <54944691+gmao-rreichle@users.noreply.github.com> Date: Wed, 20 Aug 2025 14:16:14 -0400 Subject: [PATCH 09/16] add trim() for string comparison (LDAS_Forcing.F90) --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 63b763d8..85a2b142 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -3910,7 +3910,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & do GEOSgcm_var = 1,N_GEOSgcm_vars - if (GEOSgcm_defs(GEOSgcm_var,1)=="dummy") cycle ! skip "dummy" variable (e.g., no PAR for S2S3) + if (trim(GEOSgcm_defs(GEOSgcm_var,1))=="dummy") cycle ! skip "dummy" variable (e.g., no PAR for S2S3) ! open GEOS file (G5DAS or MERRA or MERRA-2) ! From 98c2ea1f445c4841d138dfc8e005b647663ea7ae Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 20 Aug 2025 17:07:14 -0400 Subject: [PATCH 10/16] fix for special case of S2S3 FCST forcing at initialization time (LDAS_Forcing.F90) --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 56 +++++++++++++++++++------- 1 file changed, 41 insertions(+), 15 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 85a2b142..10df68ae 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -3211,6 +3211,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & character( 3) :: met_file_ext character( 3) :: precip_corr_file_ext + character( 8) :: S2S3_init_YYYYMMDD + integer :: N_GEOSgcm_vars, N_lon_tmp, N_lat_tmp real :: this_lon, this_lat @@ -3231,7 +3233,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & logical :: daily_met_files, daily_precipcorr_files logical :: is_S2S3_fcst - + integer :: nv_id, ierr, icount(3), istart(3), lonid, latid character(len=*), parameter :: Iam = 'get_GEOS' @@ -3710,6 +3712,8 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & precip_corr_file_ext = 'nc4' is_S2S3_fcst = .false. + + S2S3_init_YYYYMMDD = 'xxxxxxxx' ! character(8) if (met_tag(4:8)=='merra') then ! MERRA @@ -3825,27 +3829,29 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & met_path_bkwd, prec_path_bkwd, met_tag_bkwd, use_prec_corr ) - elseif (met_tag(1:8)=='GEOSS2S3') then ! GEOS S2S v3 + elseif (met_tag(1:8)=='GEOSS2S3') then ! GEOS S2S v3 N_GEOSgcm_vars = N_S2S3_vars - PAR_available = .false. ! S2S3 does not have PAR + PAR_available = .false. ! S2S3 does not have PAR - S2S3_file_specs = .true. + S2S3_file_specs = .true. - single_time_in_file = .false. ! FCST: monthly files, AODAS: daily files + single_time_in_file = .false. ! FCST: monthly files, AODAS: daily files if (met_tag(9:12)=='FCST' ) then - is_S2S3_fcst = .true. + is_S2S3_fcst = .true. - GEOSgcm_defs = S2S3FCST_defs + GEOSgcm_defs = S2S3FCST_defs + + S2S3_init_YYYYMMDD = met_tag(21:28) ! S2S3 FCST init YYYYMMDD; character(8) elseif (met_tag(9:13)=='AODAS') then - daily_met_files = .true. + daily_met_files = .true. - GEOSgcm_defs = S2S3AODAS_defs + GEOSgcm_defs = S2S3AODAS_defs else @@ -3912,7 +3918,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & if (trim(GEOSgcm_defs(GEOSgcm_var,1))=="dummy") cycle ! skip "dummy" variable (e.g., no PAR for S2S3) - ! open GEOS file (G5DAS or MERRA or MERRA-2) + ! open GEOS file (G5DAS or MERRA or MERRA-2 or ...) ! ! Initial "tavg1_2d_*_Nx" files may not be available. In this case, ! use first available file. For G5DAS file specs, only "PS" is affected @@ -3927,6 +3933,25 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! the file at date_time_fwd). do j=1,2 + + if (is_S2S3_fcst .and. j==1) then + + ! special S2S3 FCST case: must skip j==1 at S2S3 FCST initialization time because (monthly) file + ! exists but does not contain data for "date_time_bkwd" + + write (YYYYMMDD(1:4),'(i4.4)') date_time_tmp%year + write (YYYYMMDD(5:6),'(i2.2)') date_time_tmp%month + write (YYYYMMDD(7:8),'(i2.2)') date_time_tmp%day + + date_time_tmp%hour = 0 + date_time_tmp%min = 0 + date_time_tmp%sec = 0 + + call augment_date_time( -force_dtstep, date_time_tmp ) ! S2S3 fcst is initialized at S2S3_init_YYYYMMDD minus 3 hours + + if datetime_eq_refdatetime( date_time_tmp, date_time_inst ) cycle ! skip to j==2, i.e., try "date_time_fwd" + + end if ! determine time stamp on file and corresponding met_path, prec_path, & met_tag @@ -3950,6 +3975,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & met_path_tmp = met_path_inst prec_path_tmp = prec_path_inst met_tag_tmp = met_tag_inst + else call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown GEOSgcm_defs(2)') @@ -3979,7 +4005,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & end if - if ( file_exists .and. (.not. is_S2S3_fcst) ) then ! S2S3 FCST has monthly files + if ( file_exists ) then exit ! exit j loop after successfully finding file @@ -3996,7 +4022,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & 'read from some file for backward compatibility with ' // & 'MERRA forcing.' - write (logunit,*) 'try again with different file...' + write (logunit,*) 'try again with different file (or time)...' else @@ -5635,8 +5661,8 @@ subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_fi character( 16) :: time_stamp character( 4) :: YYYY, HHMM, day_dir character( 2) :: MM, DD - character( 8) :: S2S3_init_YYYYMMDD ! S2S3 initial date, e.g. "20160101" - character( 4) :: S2S3_ens_num ! S2S3 ensemble member, e.g. "ens1" + character( 8) :: S2S3_init_YYYYMMDD ! S2S3 fcst initialization YYYYMMDD (fcst start time is YYYYMMDD minus 3 hours) + character( 4) :: S2S3_ens_num ! S2S3 fcst ensemble member, e.g. "ens1" integer :: tmpind, tmpindend @@ -5705,7 +5731,7 @@ subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_fi ! GEOSS2S3AODAS S2S3_ens_num = trim(met_tag(15:18)) - S2S3_init_YYYYMMDD = trim(met_tag(21:28)) + S2S3_init_YYYYMMDD = met_tag(21:28) ! character(8) fname = S2S3_init_YYYYMMDD // '/' // S2S3_ens_num // '/GEOSS2S3.' // YYYY // MM // '.nc4' From 03961691165be343981f86329a3293a437b6efc6 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 20 Aug 2025 17:49:32 -0400 Subject: [PATCH 11/16] fixed syntax error in previous commit (LDAS_Forcing.F90) --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 10df68ae..5c615cfc 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -3939,9 +3939,9 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! special S2S3 FCST case: must skip j==1 at S2S3 FCST initialization time because (monthly) file ! exists but does not contain data for "date_time_bkwd" - write (YYYYMMDD(1:4),'(i4.4)') date_time_tmp%year - write (YYYYMMDD(5:6),'(i2.2)') date_time_tmp%month - write (YYYYMMDD(7:8),'(i2.2)') date_time_tmp%day + write (S2S3_init_YYYYMMDD(1:4),'(i4.4)') date_time_tmp%year + write (S2S3_init_YYYYMMDD(5:6),'(i2.2)') date_time_tmp%month + write (S2S3_init_YYYYMMDD(7:8),'(i2.2)') date_time_tmp%day date_time_tmp%hour = 0 date_time_tmp%min = 0 @@ -3949,7 +3949,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & call augment_date_time( -force_dtstep, date_time_tmp ) ! S2S3 fcst is initialized at S2S3_init_YYYYMMDD minus 3 hours - if datetime_eq_refdatetime( date_time_tmp, date_time_inst ) cycle ! skip to j==2, i.e., try "date_time_fwd" + if (datetime_eq_refdatetime( date_time_tmp, date_time_inst )) cycle ! skip to j==2, i.e., try "date_time_fwd" end if From a4538ce1c2489b34cf7626a8488496c29a88adf7 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 21 Aug 2025 08:36:41 -0400 Subject: [PATCH 12/16] add missing use statement (LDAS_Forcing.F90) --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 5c615cfc..db18b693 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -35,6 +35,7 @@ module LDAS_ForceMod use LDAS_DateTimeMod, ONLY: & date_time_type, & augment_date_time, & + datetime_eq_refdatetime, & datetime_lt_refdatetime, & datetime_le_refdatetime, & is_leap_year, & From 30d175f55ea73bcd4cf7f03728a9aedc5e220772 Mon Sep 17 00:00:00 2001 From: Qing Liu Date: Thu, 21 Aug 2025 10:55:07 -0400 Subject: [PATCH 13/16] fix syntax error --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index db18b693..2cb3d795 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -3940,9 +3940,9 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & ! special S2S3 FCST case: must skip j==1 at S2S3 FCST initialization time because (monthly) file ! exists but does not contain data for "date_time_bkwd" - write (S2S3_init_YYYYMMDD(1:4),'(i4.4)') date_time_tmp%year - write (S2S3_init_YYYYMMDD(5:6),'(i2.2)') date_time_tmp%month - write (S2S3_init_YYYYMMDD(7:8),'(i2.2)') date_time_tmp%day + read (S2S3_init_YYYYMMDD(1:4),'(i4.4)') date_time_tmp%year + read (S2S3_init_YYYYMMDD(5:6),'(i2.2)') date_time_tmp%month + read (S2S3_init_YYYYMMDD(7:8),'(i2.2)') date_time_tmp%day date_time_tmp%hour = 0 date_time_tmp%min = 0 From 2872363cde12c979544d6047719d544e9d927df1 Mon Sep 17 00:00:00 2001 From: Qing Liu Date: Thu, 21 Aug 2025 12:52:24 -0400 Subject: [PATCH 14/16] fixing error that resets single_time_in_file --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 2cb3d795..91616c3f 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -4001,7 +4001,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & daily_met_files, met_path_tmp, met_tag_tmp, & GEOSgcm_defs(GEOSgcm_var,:), met_file_ext) - single_time_in_file = .not. daily_met_files ! MERRA-2 files are daily files + single_time_in_file = .not. (daily_met_files .or. is_S2S3_fcst) ! MERRA-2 files are daily files end if From f4a4bea6bfcfbac2353bee036921691293837bfb Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 25 Aug 2025 13:54:11 -0400 Subject: [PATCH 15/16] accommodate "ens10", "ens11", ... for met forcing from S2S3 FCST (LDAS_Forcing.F90) --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 89 ++++++++++++++++++++++---- 1 file changed, 77 insertions(+), 12 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 91616c3f..76bfc06b 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -3212,6 +3212,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & character( 3) :: met_file_ext character( 3) :: precip_corr_file_ext + character( 5) :: S2S3_ens_num character( 8) :: S2S3_init_YYYYMMDD integer :: N_GEOSgcm_vars, N_lon_tmp, N_lat_tmp @@ -3846,7 +3847,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & GEOSgcm_defs = S2S3FCST_defs - S2S3_init_YYYYMMDD = met_tag(21:28) ! S2S3 FCST init YYYYMMDD; character(8) + call parse_S2S3FCST_met_tag( met_tag, S2S3_ens_num, S2S3_init_YYYYMMDD ) elseif (met_tag(9:13)=='AODAS') then @@ -4001,7 +4002,7 @@ subroutine get_GEOS( date_time, force_dtstep, met_path, met_tag, & daily_met_files, met_path_tmp, met_tag_tmp, & GEOSgcm_defs(GEOSgcm_var,:), met_file_ext) - single_time_in_file = .not. (daily_met_files .or. is_S2S3_fcst) ! MERRA-2 files are daily files + single_time_in_file = .not. (daily_met_files .or. is_S2S3_fcst) ! MERRA-2 files are daily files; S2S3FCST are monthly files end if @@ -5635,6 +5636,77 @@ subroutine parse_G5DAS_met_tag( met_path_in, met_tag_in, date_time, & end subroutine parse_G5DAS_met_tag + ! **************************************************************** + + subroutine parse_S2S3FCST_met_tag( met_tag, S2S3_ens_num, S2S3_init_YYYYMMDD ) + + ! reichle, 25 Aug 2025 + + ! parse GEOSS2S3FCST "met_tag": extract S2S3 ens number and initialization YYYYMMDD + ! + ! met_tag = "GEOSS2S3FCST__ens{XX}__[YYYYMMDD]" + ! + ! where + ! + ! ens{XX} = S2S3 ensemble member ("ens1", "ens2", ..., "ens9", "ens10", ..., "ens15") + ! YYYYMMDD = S2S3 fcst initialization YYYYMMDD (fcst start time is YYYYMMDD minus 3 hours) + ! + ! --------------------------------------------------------------------------- + + implicit none + + character(*), intent(in) :: met_tag + + character(5), intent(out) :: S2S3_ens_num + character(8), intent(out) :: S2S3_init_YYYYMMDD + + ! local variables + + integer :: is + + character(len=len(met_tag)) :: tmpstring + + character(len=*), parameter :: Iam = 'parse_S2S3FCST_met_tag' + character(len=400) :: err_msg + + ! ---------------------------------------------------------- + + err_msg = '' ! initialize error message to blank string + + if (met_tag(1:14) /= 'GEOSS2S3FCST__') err_msg = 'met_tag must start with GEOSS2S3FCST__' + + ! cut off leading 'GEOSS2S3FCST__' + + tmpstring = met_tag(15:len_trim(met_tag)) + + ! split met_tag at double underscores + + is = index(tmpstring, '__') + + if (is>0) then + + S2S3_ens_num = tmpstring(1:is-1) + + S2S3_init_YYYYMMDD = tmpstring(is+2:len_trim(tmpstring)) + + if (is<5 .or. is>6 .or. (len_trim(tmpstring)-is-1/=8)) then + + err_msg = 'ens_num or YYYYMMDD in met_tag does not match expectation' + + end if + + else + + err_msg = 'cannot find second double-underscore in met_tag' + + end if + + ! abort if something went wrong + + if (len_trim(err_code)>0) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + + end subroutine parse_S2S3FCST_met_tag + ! **************************************************************** subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_file, met_path, met_tag, & @@ -5663,7 +5735,7 @@ subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_fi character( 4) :: YYYY, HHMM, day_dir character( 2) :: MM, DD character( 8) :: S2S3_init_YYYYMMDD ! S2S3 fcst initialization YYYYMMDD (fcst start time is YYYYMMDD minus 3 hours) - character( 4) :: S2S3_ens_num ! S2S3 fcst ensemble member, e.g. "ens1" + character( 4) :: S2S3_ens_num ! S2S3 fcst ensemble member, e.g. "ens1", "ens12" integer :: tmpind, tmpindend @@ -5724,17 +5796,10 @@ subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_fi day_dir = 'D' // DD // '/' elseif (met_tag(1:12) == 'GEOSS2S3FCST') then - - ! parse met_tag - - ! 1234567890123456789012345678 - ! GEOSS2S3FCST__ensX__YYYYMMDD - ! GEOSS2S3AODAS - S2S3_ens_num = trim(met_tag(15:18)) - S2S3_init_YYYYMMDD = met_tag(21:28) ! character(8) + call parse_S2S3FCST_met_tag( met_tag, S2S3_ens_num, S2S3_init_YYYYMMDD ) - fname = S2S3_init_YYYYMMDD // '/' // S2S3_ens_num // '/GEOSS2S3.' // YYYY // MM // '.nc4' + fname = S2S3_init_YYYYMMDD // '/' // trim(S2S3_ens_num) // '/GEOSS2S3.' // YYYY // MM // '.nc4' else From 28948433d5f7492d32d3ca517174a4630e7cc3a8 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 25 Aug 2025 14:06:59 -0400 Subject: [PATCH 16/16] fixed syntax errors in previous commit (LDAS_Forcing.F90) --- GEOSmetforce_GridComp/LDAS_Forcing.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GEOSmetforce_GridComp/LDAS_Forcing.F90 b/GEOSmetforce_GridComp/LDAS_Forcing.F90 index 76bfc06b..b0b201da 100644 --- a/GEOSmetforce_GridComp/LDAS_Forcing.F90 +++ b/GEOSmetforce_GridComp/LDAS_Forcing.F90 @@ -5703,7 +5703,7 @@ subroutine parse_S2S3FCST_met_tag( met_tag, S2S3_ens_num, S2S3_init_YYYYMMDD ) ! abort if something went wrong - if (len_trim(err_code)>0) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + if (len_trim(err_msg)>0) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end subroutine parse_S2S3FCST_met_tag @@ -5735,7 +5735,7 @@ subroutine get_GEOS_forcing_filename(fname_full,file_exists, date_time, daily_fi character( 4) :: YYYY, HHMM, day_dir character( 2) :: MM, DD character( 8) :: S2S3_init_YYYYMMDD ! S2S3 fcst initialization YYYYMMDD (fcst start time is YYYYMMDD minus 3 hours) - character( 4) :: S2S3_ens_num ! S2S3 fcst ensemble member, e.g. "ens1", "ens12" + character( 5) :: S2S3_ens_num ! S2S3 fcst ensemble member, e.g. "ens1", "ens12" integer :: tmpind, tmpindend