#include "MAPL_ErrLog.h" #include "unused_dummy.H" submodule (HistoryTrajectoryMod) HistoryTrajectory_implement use ESMF use MAPL_ErrorHandlingMod use MAPL_KeywordEnforcerMod use LocStreamFactoryMod use MAPL_LocstreamRegridderMod use MAPL_FileMetadataUtilsMod use pFIO use MAPL_GriddedIOItemMod use MAPL_GriddedIOItemVectorMod use MAPL_TimeDataMod use MAPL_VerticalDataMod use MAPL_BaseMod use MAPL_CommsMod use MAPL_SortMod use MAPL_NetCDF use MAPL_StringTemplate use Plain_netCDF_Time use MAPL_ObsUtilMod use MPI, only : MPI_INTEGER, MPI_REAL, MPI_REAL8 use, intrinsic :: iso_fortran_env, only: REAL32 use, intrinsic :: iso_fortran_env, only: REAL64 use, intrinsic :: iso_fortran_env, only: INT64 implicit none contains module procedure HistoryTrajectory_from_config use BinIOMod use pflogger, only : Logger, logging type(ESMF_Time) :: currTime type(ESMF_TimeInterval) :: epoch_frequency type(ESMF_TimeInterval) :: obs_time_span integer :: time_integer, second integer :: status character(len=ESMF_MAXSTR) :: STR1, line character(len=ESMF_MAXSTR) :: symd, shms integer :: nline, col integer, allocatable :: ncol(:) character(len=ESMF_MAXSTR), allocatable :: word(:) integer :: nobs, head, jvar logical :: tend integer :: i, j, k, k2, M integer :: count, idx integer :: unitr, unitw type(GriddedIOitem) :: item type(Logger), pointer :: lgr traj%clock=clock if (present(GENSTATE)) traj%GENSTATE => GENSTATE call ESMF_ClockGet ( clock, CurrTime=currTime, _RC ) call ESMF_ConfigGetAttribute(config, value=time_integer, label=trim(string)//'Epoch:', default=0, _RC) _ASSERT(time_integer /= 0, 'Epoch value in config wrong') second = hms_2_s(time_integer) call ESMF_TimeIntervalSet(epoch_frequency, s=second, _RC) traj%Epoch = time_integer traj%RingTime = currTime traj%epoch_frequency = epoch_frequency traj%alarm = ESMF_AlarmCreate( clock=clock, RingInterval=epoch_frequency, & RingTime=traj%RingTime, sticky=.false., _RC ) call ESMF_ConfigGetAttribute(config, value=traj%index_name_x, default="", & label=trim(string) // 'index_name_x:', _RC) call ESMF_ConfigGetAttribute(config, value=traj%var_name_lon_full, default="", & label=trim(string) // 'var_name_lon:', _RC) call ESMF_ConfigGetAttribute(config, value=traj%var_name_lat_full, default="", & label=trim(string) // 'var_name_lat:', _RC) call ESMF_ConfigGetAttribute(config, value=traj%var_name_time_full, default="", & label=trim(string) // 'var_name_time:', _RC) call ESMF_ConfigGetAttribute(config, value=STR1, default="", & label=trim(string) // 'obs_file_begin:', _RC) if (trim(STR1)=='') then traj%obsfile_start_time = currTime call ESMF_TimeGet(currTime, timestring=STR1, _RC) if (mapl_am_I_root()) then write(6,105) 'obs_file_begin missing, default = currTime :', trim(STR1) endif else call ESMF_TimeSet(traj%obsfile_start_time, STR1, _RC) if (mapl_am_I_root()) then write(6,105) 'obs_file_begin provided: ', trim(STR1) end if end if call ESMF_ConfigGetAttribute(config, value=STR1, default="", & label=trim(string) // 'obs_file_end:', _RC) if (trim(STR1)=='') then call ESMF_TimeIntervalSet(obs_time_span, d=14, _RC) traj%obsfile_end_time = traj%obsfile_start_time + obs_time_span call ESMF_TimeGet(traj%obsfile_end_time, timestring=STR1, _RC) if (mapl_am_I_root()) then write(6,105) 'obs_file_end missing, default = begin+14D:', trim(STR1) endif else call ESMF_TimeSet(traj%obsfile_end_time, STR1, _RC) if (mapl_am_I_root()) then write(6,105) 'obs_file_end provided:', trim(STR1) end if end if call ESMF_ConfigGetAttribute(config, value=STR1, default="", & label=trim(string) // 'obs_file_interval:', _RC) _ASSERT(STR1/='', 'fatal error: obs_file_interval not provided in RC file') if (mapl_am_I_root()) write(6,105) 'obs_file_interval:', trim(STR1) if (mapl_am_I_root()) write(6,106) 'Epoch (second) :', second i= index( trim(STR1), ' ' ) if (i>0) then symd=STR1(1:i-1) shms=STR1(i+1:) else symd='' shms=trim(STR1) endif call convert_twostring_2_esmfinterval (symd, shms, traj%obsfile_interval, _RC) traj%active = .true. ! __ s1. overall print call ESMF_ConfigGetDim(config, nline, col, label=trim(string)//'obs_files:', rc=rc) _ASSERT(rc==0 .AND. nline > 0, 'obs_files not found') !! write(6,*) 'nline, col', nline, col allocate(ncol(1:nline), _STAT) call ESMF_ConfigFindLabel( config, trim(string)//'obs_files:', _RC ) do i = 1, nline call ESMF_ConfigNextLine(config, _RC) ncol(i) = ESMF_ConfigGetLen(config, _RC) !! write(6,*) 'line', i, 'ncol(i)', ncol(i) enddo ! __ s2. find nobs && distinguish design with vs wo '------' nobs=0 call ESMF_ConfigFindLabel( config, trim(string)//'obs_files:', _RC) do i=1, nline call ESMF_ConfigNextLine( config, tableEnd=tend, _RC) call ESMF_ConfigGetAttribute( config, STR1, _RC) if ( index(trim(STR1), '-----') > 0 ) nobs=nobs+1 enddo ! __ s3. retrieve template and geoval, set metadata file_handle lgr => logging%get_logger('HISTORY.sampler') if ( nobs == 0 ) then ! ! treatment-1: ! _FAIL('this setting in HISTORY.rc obs_files: is not supported, stop') traj%nobs_type = nline ! here .rc format cannot have empty spaces allocate (traj%obs(nline), _STAT) call ESMF_ConfigFindLabel( config, trim(string)//'obs_files:', _RC) do i=1, nline call ESMF_ConfigNextLine( config, tableEnd=tend, _RC) call ESMF_ConfigGetAttribute( config, traj%obs(i)%input_template, _RC) traj%obs(i)%export_all_geoval = .true. enddo else ! !-- selectively output geovals ! treatment-2: ! traj%nobs_type = nobs allocate (traj%obs(nobs), _STAT) ! nobs=0 ! reuse counter head=1 jvar=0 ! ! count '------' in history.rc as special markers for ngeoval ! call ESMF_ConfigFindLabel(config, trim(string)//'obs_files:', _RC) do i=1, nline call ESMF_ConfigNextLine(config, tableEnd=tend, _RC) M = ncol(i) _ASSERT(M>=1, '# of columns should be >= 1') allocate (word(M), _STAT) count=0 do col=1, M call ESMF_ConfigGetAttribute(config, STR1, _RC) if (trim(STR1)/=',') then count=count+1 word(count) = extract_unquoted_item(STR1) end if enddo if (count ==1 .or. count==2) then ! 1-item case: file template or one-var ! 2-item : var1 , 'root' case STR1=trim(word(1)) else ! 3-item : var1 , 'root', var1_alias case STR1=trim(word(3)) end if deallocate(word, _STAT) if ( index(trim(STR1), '-----') == 0 ) then if (head==1 .AND. trim(STR1)/='') then nobs=nobs+1 traj%obs(nobs)%input_template = trim(STR1) traj%obs(nobs)%export_all_geoval = .false. head=0 else if (trim(STR1)/='') then jvar=jvar+1 idx = index(STR1,";") if (idx==0) then traj%obs(nobs)%geoval_xname(jvar) = STR1 else traj%obs(nobs)%geoval_xname(jvar) = trim(STR1(1:idx-1)) traj%obs(nobs)%geoval_yname(jvar) = trim(STR1(idx+1:)) end if end if end if else traj%obs(nobs)%ngeoval=jvar head=1 jvar=0 endif enddo end if do k=1, traj%nobs_type allocate (traj%obs(k)%metadata, _STAT) if (mapl_am_i_root()) then allocate (traj%obs(k)%file_handle, _STAT) end if end do call lgr%debug('%a %i8', 'nobs_type=', traj%nobs_type) do i=1, traj%nobs_type call lgr%debug('%a %i4 %a %a', 'obs(', i, ') input_template =', & trim(traj%obs(i)%input_template)) k=index(traj%obs(i)%input_template , '/', back=.true.) j=index(traj%obs(i)%input_template(k+1:), '%') if (j>0) then ! normal case: geos_atmosphere/aircraft.%y4%m2%d2T%h2%n2%S2Z.nc4 traj%obs(i)%name = traj%obs(i)%input_template(k+1:k+j-1) else ! different case: Y%y4/M%m2/.../this.nc or ./this k2=index(traj%obs(i)%input_template(k+1:), '.') if (k2>0) then traj%obs(i)%name = traj%obs(i)%input_template(k+1:k+k2) else traj%obs(i)%name = trim(traj%obs(i)%input_template(k+1:)) end if end if end do _RETURN(_SUCCESS) 105 format (1x,a,2x,a) 106 format (1x,a,2x,i8) end procedure HistoryTrajectory_from_config ! !-- integrate both initialize and reinitialize ! module procedure initialize_ integer :: status type(ESMF_Grid) :: grid type(variable) :: v type(GriddedIOitemVectorIterator) :: iter type(GriddedIOitem), pointer :: item type(ESMF_Time) :: currTime integer :: k ! if (mapl_am_i_root()) write(6,'(2x,a,10(2x,L5))') & ! 'traj initialize_ : present(reinitialize), reinitialize =', & ! present(reinitialize), reinitialize if (.not. present(reinitialize)) then if(present(bundle)) this%bundle=bundle if(present(items)) this%items=items if(present(timeInfo)) this%time_info=timeInfo if (present(vdata)) then this%vdata=vdata else this%vdata=VerticalData(_RC) end if !if (mapl_am_i_root()) write(6,'(2x,a,10(2x,L5))') & ! 'traj initialize_ : initialize : not present ' else if (reinitialize) then do k=1, this%nobs_type allocate (this%obs(k)%metadata, _STAT) if (mapl_am_i_root()) then allocate (this%obs(k)%file_handle, _STAT) end if end do !if (mapl_am_i_root()) write(6,'(2x,a,10(2x,L5))') & ! 'traj initialize_ : initialize : TRUE' end if end if do k=1, this%nobs_type call this%vdata%append_vertical_metadata(this%obs(k)%metadata,this%bundle,_RC) end do this%do_vertical_regrid = (this%vdata%regrid_type /= VERTICAL_METHOD_NONE) if (this%vdata%regrid_type == VERTICAL_METHOD_ETA2LEV) call this%vdata%get_interpolating_variable(this%bundle,_RC) call ESMF_ClockGet (this%clock, CurrTime=currTime, _RC) call get_obsfile_Tbracket_from_epoch(currTime, & this%obsfile_start_time, this%obsfile_end_time, & this%obsfile_interval, this%epoch_frequency, & this%obsfile_Ts_index, this%obsfile_Te_index, _RC) if (this%obsfile_Te_index < 0) then if (mapl_am_I_root()) then write(6,*) "model start time is earlier than obsfile_start_time" write(6,*) "solution: adjust obsfile_start_time and Epoch in rc file" end if _FAIL("obs file not found at init time") endif call this%create_grid(_RC) call ESMF_FieldBundleGet(this%bundle,grid=grid,_RC) this%regridder = LocStreamRegridder(grid,this%LS_ds,_RC) this%output_bundle = this%create_new_bundle(_RC) this%acc_bundle = this%create_new_bundle(_RC) do k=1, this%nobs_type call this%obs(k)%metadata%add_dimension(this%index_name_x, this%obs(k)%nobs_epoch) if (this%time_info%integer_time) then v = Variable(type=PFIO_INT32,dimensions=this%index_name_x) else v = Variable(type=PFIO_REAL64,dimensions=this%index_name_x) end if call v%add_attribute('units', this%datetime_units) call v%add_attribute('long_name', 'dateTime') call this%obs(k)%metadata%add_variable(this%var_name_time,v) v = Variable(type=PFIO_INT32,dimensions=this%index_name_x) call v%add_attribute('units', '1') call v%add_attribute('long_name', 'Location index in corresponding IODA file') call this%obs(k)%metadata%add_variable(this%location_index_name,v) v = variable(type=PFIO_REAL64,dimensions=this%index_name_x) call v%add_attribute('units','degrees_east') call v%add_attribute('long_name','longitude') call this%obs(k)%metadata%add_variable(this%var_name_lon,v) v = variable(type=PFIO_REAL64,dimensions=this%index_name_x) call v%add_attribute('units','degrees_north') call v%add_attribute('long_name','latitude') call this%obs(k)%metadata%add_variable(this%var_name_lat,v) end do ! push varible names down to each obs(k); see create_metadata_variable iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() !! print*, 'list item%xname', trim(item%xname) if (item%itemType == ItemTypeScalar) then call this%create_variable(item%xname,_RC) else if (item%itemType == ItemTypeVector) then call this%create_variable(item%xname,_RC) call this%create_variable(item%yname,_RC) end if call iter%next() enddo _RETURN(_SUCCESS) end procedure initialize_ module procedure create_metadata_variable type(ESMF_Field) :: field type(variable) :: v logical :: is_present integer :: field_rank, status character(len=ESMF_MAXSTR) :: var_name,long_name,units,vdims type(ESMF_Info) :: infoh integer :: k, ig call ESMF_FieldBundleGet(this%bundle,vname,field=field,_RC) call ESMF_FieldGet(field,name=var_name,rank=field_rank,_RC) call ESMF_InfoGetFromHost(field,infoh,_RC) is_present = ESMF_InfoIsPresent(infoh,"LONG_NAME",_RC) if ( is_present ) then call ESMF_InfoGet(infoh,"LONG_NAME",long_name,_RC) else long_name = var_name endif is_present = ESMF_InfoIsPresent(infoh,"UNITS",_RC) if ( is_present ) then call ESMF_InfoGet(infoh,"UNITS",units,_RC) else units = 'unknown' endif if (field_rank==2) then vdims = this%index_name_x else if (field_rank==3) then vdims = trim(this%index_name_x)//",lev" end if v = variable(type=PFIO_REAL32,dimensions=trim(vdims)) call v%add_attribute('units',trim(units)) call v%add_attribute('long_name',trim(long_name)) call v%add_attribute('missing_value',MAPL_UNDEF) call v%add_attribute('_FillValue',MAPL_UNDEF) call v%add_attribute('valid_range',(/-MAPL_UNDEF,MAPL_UNDEF/)) do k = 1, this%nobs_type do ig = 1, this%obs(k)%ngeoval if (trim(var_name) == trim(this%obs(k)%geoval_xname(ig))) then call this%obs(k)%metadata%add_variable(trim(var_name),v,_RC) !! if (mapl_am_i_root()) write(6, '(2x,a,/,10(2x,a))') & !! 'Traj: create_metadata_variable: vname, var_name, this%obs(k)%geoval_xname(ig)', & !! trim(vname), trim(var_name), trim(this%obs(k)%geoval_xname(ig)) endif enddo enddo _RETURN(_SUCCESS) end procedure create_metadata_variable module procedure create_new_bundle type(GriddedIOitemVectorIterator) :: iter type(GriddedIOitem), pointer :: item type(ESMF_Field) :: src_field,dst_field integer :: rank,lb(1),ub(1) integer :: status new_bundle = ESMF_FieldBundleCreate(_RC) iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() if (item%itemType == ItemTypeScalar) then call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) call ESMF_FieldGet(src_field,rank=rank,_RC) if (rank==2) then dst_field = ESMF_FieldCreate(this%LS_ds,name=trim(item%xname), & typekind=ESMF_TYPEKIND_R4,_RC) else if (rank==3) then call ESMF_FieldGet(src_field,ungriddedLBound=lb,ungriddedUBound=ub,_RC) if (this%vdata%lm/=(ub(1)-lb(1)+1)) then lb(1)=1 ub(1)=this%vdata%lm end if dst_field = ESMF_FieldCreate(this%LS_ds,name=trim(item%xname), & typekind=ESMF_TYPEKIND_R4,ungriddedLBound=lb,ungriddedUBound=ub,_RC) end if call MAPL_FieldBundleAdd(new_bundle,dst_field,_RC) else if (item%itemType == ItemTypeVector) then !! _FAIL("ItemTypeVector not yet supported") end if call iter%next() enddo _RETURN(_SUCCESS) end procedure create_new_bundle module procedure create_file_handle use pflogger, only : Logger, logging integer :: status integer :: k character(len=ESMF_MAXSTR) :: filename type(Logger), pointer :: lgr if (.NOT. this%active) then _RETURN(ESMF_SUCCESS) endif if (this%nobs_epoch_sum==0) then rc=0 return endif lgr => logging%get_logger('HISTORY.sampler') do k=1, this%nobs_type call this%obs(k)%metadata%modify_dimension(this%index_name_x, this%obs(k)%nobs_epoch) enddo if (mapl_am_I_root()) then do k=1, this%nobs_type if (this%obs(k)%nobs_epoch > 0) then filename=trim(this%obs(k)%name)//trim(filename_suffix) call lgr%debug('%a %a', & "Sampling to new file : ",trim(filename)) call this%obs(k)%file_handle%create(trim(filename),_RC) call this%obs(k)%file_handle%write(this%obs(k)%metadata,_RC) end if enddo end if _RETURN(_SUCCESS) end procedure create_file_handle module procedure close_file_handle integer :: status integer :: k if (.NOT. this%active) then _RETURN(ESMF_SUCCESS) endif if (this%nobs_epoch_sum==0) then rc=0 return endif if (mapl_am_I_root()) then do k=1, this%nobs_type if (this%obs(k)%nobs_epoch > 0) then call this%obs(k)%file_handle%close(_RC) end if end do end if _RETURN(_SUCCESS) end procedure close_file_handle module procedure create_grid use pflogger, only: Logger, logging character(len=ESMF_MAXSTR) :: filename integer(ESMF_KIND_I4) :: num_times integer :: len integer :: len_full integer :: status type(Logger), pointer :: lgr character(len=ESMF_MAXSTR) :: grp_name character(len=ESMF_MAXSTR) :: timeunits_file character :: new_char(ESMF_MAXSTR) real(REAL64), allocatable :: lons_full(:), lats_full(:) real(REAL64), allocatable :: times_R8_full(:) real(REAL64) :: t_shift integer, allocatable :: obstype_id_full(:) integer, allocatable :: location_index_ioda_full(:) integer, allocatable :: IA_full(:) real(ESMF_KIND_R8), pointer :: ptAT(:) type(ESMF_routehandle) :: RH type(ESMF_Time) :: timeset(2) type(ESMF_Time) :: current_time type(ESMF_Time) :: time0 type(ESMF_TimeInterval) :: dt type(ESMF_Grid) :: grid type(ESMF_VM) :: vm integer :: mypet, petcount, mpic integer :: i, j, k, L, ii, jj integer :: fid_s, fid_e integer(kind=ESMF_KIND_I8) :: j0, j1 integer(kind=ESMF_KIND_I8) :: jt1, jt2 integer(kind=ESMF_KIND_I8) :: nstart, nend real(kind=ESMF_KIND_R8) :: jx0, jx1 integer :: nx, nx_sum integer :: n0 integer :: arr(1) integer :: sec integer, allocatable :: ix(:) ! counter for each obs(k)%nobs_epoch integer :: nx2 logical :: EX ! file logical :: zero_obs integer, allocatable :: sendcount(:), displs(:) integer :: recvcount integer :: is, ie, ierr integer :: M, N, ip real(REAL64), allocatable :: lons_chunk(:) real(REAL64), allocatable :: lats_chunk(:) real(REAL64), allocatable :: times_R8_chunk(:) lgr => logging%get_logger('HISTORY.sampler') call ESMF_VMGetCurrent(vm,_RC) call ESMF_VMGet(vm, mpiCommunicator=mpic, petCount=petCount, localPet=mypet, _RC) if (this%index_name_x == '') then ! !-- non IODA case / non netCDF ! _FAIL('non-IODA format is not implemented here') end if ! !-- IODA case ! i=index(this%var_name_lon_full, '/') if (i==0) then grp_name = '' call lgr%debug('%a', 'grp_name not found') else grp_name = this%var_name_lon_full(1:i-1) end if this%var_name_lon = this%var_name_lon_full(i+1:) i=index(this%var_name_lat_full, '/') this%var_name_lat = this%var_name_lat_full(i+1:) i=index(this%var_name_time_full, '/') this%var_name_time= this%var_name_time_full(i+1:) this%location_index_name = 'location_index_in_iodafile' call lgr%debug('%a', 'grp_name,this%index_name_x,this%var_name_lon,this%var_name_lat,this%var_name_time') call lgr%debug('%a %a %a %a %a', & trim(grp_name),trim(this%index_name_x),trim(this%var_name_lon),& trim(this%var_name_lat),trim(this%var_name_time)) L=0 fid_s=this%obsfile_Ts_index fid_e=this%obsfile_Te_index call lgr%debug('%a %i10 %i10', & 'fid_s, fid_e', fid_s, fid_e) arr(1)=0 ! len_full if (mapl_am_I_root()) then len = 0 do k=1, this%nobs_type j = max (fid_s, L) do while (j<=fid_e) filename = get_filename_from_template_use_index( & this%obsfile_start_time, this%obsfile_interval, & j, this%obs(k)%input_template, EX, _RC) if (EX) then call lgr%debug('%a %i10', 'exist: filename fid j :', j) call lgr%debug('%a %a', 'exist: true filename :', trim(filename)) call get_ncfile_dimension(filename, tdim=num_times, key_time=this%index_name_x, _RC) len = len + num_times else call lgr%debug('%a %i10', 'non-exist: filename fid j :', j) call lgr%debug('%a %a', 'non-exist: missing filename:', trim(filename)) end if j=j+1 enddo enddo arr(1)=len if (len>0) then allocate(lons_full(len),lats_full(len),_STAT) allocate(times_R8_full(len),_STAT) allocate(obstype_id_full(len),_STAT) allocate(location_index_ioda_full(len),_STAT) allocate(IA_full(len),_STAT) call lgr%debug('%a %i12', 'nobs from input file:', len) len = 0 ii = 0 do k=1, this%nobs_type j = max (fid_s, L) do while (j<=fid_e) filename = get_filename_from_template_use_index( & this%obsfile_start_time, this%obsfile_interval, & j, this%obs(k)%input_template, EX, _RC) if (EX) then ii = ii + 1 call get_ncfile_dimension(trim(filename), tdim=num_times, key_time=this%index_name_x, _RC) call get_v1d_netcdf_R8 (filename, this%var_name_lon, lons_full(len+1:), num_times, group_name=grp_name) call get_v1d_netcdf_R8 (filename, this%var_name_lat, lats_full(len+1:), num_times, group_name=grp_name) call get_v1d_netcdf_R8 (filename, this%var_name_time, times_R8_full(len+1:), num_times, group_name=grp_name) call get_attribute_from_group (filename, grp_name, this%var_name_time, "units", timeunits_file) if (ii == 1) then this%datetime_units = trim(timeunits_file) call lgr%debug('%a %a', 'datetime_units from 1st file:', trim(timeunits_file)) end if call diff_two_timeunits (this%datetime_units, timeunits_file, t_shift, _RC) times_R8_full(len+1:len+num_times) = times_R8_full(len+1:len+num_times) + t_shift obstype_id_full(len+1:len+num_times) = k do jj = 1, num_times location_index_ioda_full(len+jj) = jj end do len = len + num_times end if j=j+1 enddo enddo end if end if call ESMF_VMAllFullReduce(vm, sendData=arr, recvData=nx_sum, & count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) if (nx_sum == 0) then allocate(this%lons(0),this%lats(0),_STAT) allocate(this%times_R8(0),_STAT) allocate(this%obstype_id(0),_STAT) allocate(this%location_index_ioda(0),_STAT) this%epoch_index(1:2) = 0 this%nobs_epoch = 0 this%nobs_epoch_sum = 0 ! ! empty shell to keep regridding and destroy_RH_LS to work ! this%locstream_factory = LocStreamFactory(this%lons,this%lats,_RC) this%LS_rt = this%locstream_factory%create_locstream(_RC) call ESMF_FieldBundleGet(this%bundle,grid=grid,_RC) this%LS_ds = this%locstream_factory%create_locstream(grid=grid,_RC) this%fieldB = ESMF_FieldCreate (this%LS_ds, name='B_time', typekind=ESMF_TYPEKIND_R8, _RC) call ESMF_FieldGet( this%fieldB, localDE=0, farrayPtr=this%obsTime) this%obsTime= -1.d0 rc = 0 return end if call MAPL_CommsBcast(vm, this%datetime_units, N=ESMF_MAXSTR, ROOT=MAPL_Root, _RC) if (mapl_am_I_root()) then call sort_index (times_R8_full, IA_full, _RC) location_index_ioda_full(:) = IA_full(:) ! NVHPC dies with NVFORTRAN-S-0155-Could not resolve generic procedure sort_multi_arrays_by_time call sort_four_arrays_by_time(lons_full, lats_full, times_R8_full, obstype_id_full, _RC) call ESMF_ClockGet(this%clock,currTime=current_time,_RC) timeset(1) = current_time timeset(2) = current_time + this%epoch_frequency call time_esmf_2_nc_int (timeset(1), this%datetime_units, j0, _RC) sec = hms_2_s(this%Epoch) j1 = j0 + int(sec, kind=ESMF_KIND_I8) jx0 = real ( j0, kind=ESMF_KIND_R8) jx1 = real ( j1, kind=ESMF_KIND_R8) nstart=1; nend=size(times_R8_full) call bisect( times_R8_full, jx0, jt1, n_LB=int(nstart, ESMF_KIND_I8), n_UB=int(nend, ESMF_KIND_I8), rc=rc) call bisect( times_R8_full, jx1, jt2, n_LB=int(nstart, ESMF_KIND_I8), n_UB=int(nend, ESMF_KIND_I8), rc=rc) call lgr%debug ('%a %i20 %i20', 'nstart, nend', nstart, nend) call lgr%debug ('%a %f20.1 %f20.1', 'j0[currT] j1[T+Epoch] w.r.t. timeunit ', jx0, jx1) call lgr%debug ('%a %f20.1 %f20.1', 'x0[times(1)] xn[times(N)] w.r.t. timeunit ', & times_R8_full(1), times_R8_full(nend)) call lgr%debug ('%a %i20 %i20', 'jt1, jt2 [final intercepted position]', jt1, jt2) ! if (jt1==jt2) then ! _FAIL('Epoch Time is too small, empty grid is generated, increase Epoch') ! endif !-- shift the zero item to index 1 zero_obs = .false. if (jt1/=jt2) then zero_obs = .false. if (jt1==0) jt1=1 else ! at most one obs point exist, set it .true. zero_obs = .true. !! if (jt1==0) jt1=1 end if ! !-- exclude the out-of-range case ! if ( zero_obs ) then allocate(this%lons(0),this%lats(0),_STAT) allocate(this%times_R8(0),_STAT) allocate(this%obstype_id(0),_STAT) allocate(this%location_index_ioda(0),_STAT) this%epoch_index(1:2)=0 this%nobs_epoch = 0 nx=0 arr(1)=nx else !! doulbe check ! (x1, x2] design in bisect this%epoch_index(1)= jt1 + 1 !! ! (x1, x2] design in bisect !! if (jt1==0) then !! this%epoch_index(1)= 1 !! else !! this%epoch_index(1)= jt1 !! endif _ASSERT(jt2<=len, 'bisect index for this%epoch_index(2) failed') if (jt2==0) then this%epoch_index(2)= 1 else this%epoch_index(2)= jt2 endif nx= this%epoch_index(2) - this%epoch_index(1) + 1 this%nobs_epoch = nx allocate(this%lons(nx),this%lats(nx),_STAT) allocate(this%times_R8(nx),_STAT) allocate(this%obstype_id(nx),_STAT) allocate(this%location_index_ioda(nx),_STAT) j=this%epoch_index(1) do i=1, nx this%lons(i) = lons_full(j) this%lats(i) = lats_full(j) this%times_R8(i) = times_R8_full(j) this%obstype_id(i) = obstype_id_full(j) this%location_index_ioda(i) = location_index_ioda_full(j) j=j+1 enddo arr(1)=nx do k=1, this%nobs_type this%obs(k)%nobs_epoch = 0 enddo do j = this%epoch_index(1), this%epoch_index(2) k = obstype_id_full(j) this%obs(k)%nobs_epoch = this%obs(k)%nobs_epoch + 1 enddo do k=1, this%nobs_type nx2 = this%obs(k)%nobs_epoch allocate (this%obs(k)%lons(nx2), _STAT) allocate (this%obs(k)%lats(nx2), _STAT) allocate (this%obs(k)%times_R8(nx2), _STAT) allocate (this%obs(k)%location_index_ioda(nx2), _STAT) enddo allocate(ix(this%nobs_type), _STAT) ix(:)=0 j=this%epoch_index(1) do i=1, nx k = obstype_id_full(j) ix(k) = ix(k) + 1 this%obs(k)%lons(ix(k)) = lons_full(j) this%obs(k)%lats(ix(k)) = lats_full(j) this%obs(k)%times_R8(ix(k)) = times_R8_full(j) this%obs(k)%location_index_ioda(ix(k)) = location_index_ioda_full(j) !if (mod(k,10**8)==1) then ! write(6,*) 'this%obs(k)%times_R8(ix(k))', this%obs(k)%times_R8(ix(k)) !endif j=j+1 enddo deallocate(ix, _STAT) deallocate(lons_full, lats_full, times_R8_full, obstype_id_full, location_index_ioda_full, _STAT) call lgr%debug('%a %i12 %i12 %i12', & 'epoch_index(1:2), nx', this%epoch_index(1), & this%epoch_index(2), this%nobs_epoch) do k=1, this%nobs_type call lgr%debug('%a %i4 %a %i12', & 'obs(', k, ')%nobs_epoch', this%obs(k)%nobs_epoch ) enddo end if else allocate(this%lons(0),this%lats(0),_STAT) allocate(this%times_R8(0),_STAT) allocate(this%obstype_id(0),_STAT) allocate(this%location_index_ioda(0),_STAT) this%epoch_index(1:2)=0 this%nobs_epoch = 0 nx=0 arr(1)=nx endif call ESMF_VMAllFullReduce(vm, sendData=arr, recvData=nx_sum, & count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) this%nobs_epoch_sum = nx_sum call lgr%debug('%a %i20', 'nobservation points=', nx_sum) ! !__ s1. distrubute data chunk for the locstream points : mpi_scatterV !__ s2. create LS on parallel processors ! caution about zero-sized array for MPI ! ip = mypet ! 0 to M-1 N = nx_sum M = petCount recvcount = int(ip+1, INT64) * int(N, INT64) / int(M, INT64) - & int(ip , INT64) * int(N, INT64) / int(M, INT64) call lgr%debug('%a %i12 %i12', 'ip, recvcount', ip, recvcount) allocate ( sendcount (petCount) ) allocate ( displs (petCount) ) do ip=0, M-1 sendcount(ip+1) = int(ip+1, INT64) * int(N, INT64) / int(M, INT64) - & int(ip , INT64) * int(N, INT64) / int(M, INT64) end do displs(1)=0 do i = 2, petCount displs(i) = displs(i-1) + sendcount(i-1) end do allocate ( lons_chunk (recvcount) ) allocate ( lats_chunk (recvcount) ) allocate ( times_R8_chunk (recvcount) ) arr(1) = recvcount call ESMF_VMAllFullReduce(vm, sendData=arr, recvData=nx2, & count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) _ASSERT( nx2 == nx_sum, 'Erorr in recvcount' ) call MPI_Scatterv( this%lons, sendcount, & displs, MPI_REAL8, lons_chunk, & recvcount, MPI_REAL8, 0, mpic, ierr) _VERIFY(ierr) call MPI_Scatterv( this%lats, sendcount, & displs, MPI_REAL8, lats_chunk, & recvcount, MPI_REAL8, 0, mpic, ierr) _VERIFY(ierr) call MPI_Scatterv( this%times_R8, sendcount, & displs, MPI_REAL8, times_R8_chunk, & recvcount, MPI_REAL8, 0, mpic, ierr) _VERIFY(ierr) ! -- root this%locstream_factory = LocStreamFactory(this%lons,this%lats,_RC) this%LS_rt = this%locstream_factory%create_locstream(_RC) ! -- proc this%locstream_factory = LocStreamFactory(lons_chunk,lats_chunk,_RC) this%LS_chunk = this%locstream_factory%create_locstream_on_proc(_RC) call ESMF_FieldBundleGet(this%bundle,grid=grid,_RC) this%LS_ds = this%locstream_factory%create_locstream_on_proc(grid=grid,_RC) this%fieldA = ESMF_FieldCreate (this%LS_chunk, name='A_time', typekind=ESMF_TYPEKIND_R8, _RC) this%fieldB = ESMF_FieldCreate (this%LS_ds, name='B_time', typekind=ESMF_TYPEKIND_R8, _RC) call ESMF_FieldGet( this%fieldA, localDE=0, farrayPtr=ptAT) call ESMF_FieldGet( this%fieldB, localDE=0, farrayPtr=this%obsTime) ptAT(:) = times_R8_chunk(:) this%obsTime= -1.d0 call ESMF_FieldRedistStore (this%fieldA, this%fieldB, RH, _RC) call ESMF_FieldRedist (this%fieldA, this%fieldB, RH, _RC) call ESMF_FieldRedistRelease(RH, noGarbage=.true., _RC) call ESMF_FieldDestroy(this%fieldA,nogarbage=.true.,_RC) ! defer destroy fieldB at regen_grid step ! _RETURN(_SUCCESS) end procedure create_grid module procedure append_file use pflogger, only: Logger, logging type(GriddedIOitemVectorIterator) :: iter type(GriddedIOitem), pointer :: item type(ESMF_RouteHandle) :: RH type(Logger), pointer :: lgr type(ESMF_Field) :: src_field, dst_field type(ESMF_Field) :: acc_field type(ESMF_Field) :: acc_field_2d_rt, acc_field_3d_rt real(REAL32), pointer :: p_acc_3d(:,:),p_acc_2d(:) real(REAL32), pointer :: p_acc_rt_2d(:) real(REAL32), pointer :: p_src(:,:),p_dst(:,:), p_dst_t(:,:) ! _t: transpose real(REAL32), pointer :: p_dst_rt(:,:), p_acc_rt_3d(:,:) real(REAL32), pointer :: pt1(:), pt2(:) type(ESMF_Field) :: acc_field_2d_chunk, acc_field_3d_chunk, chunk_field real(REAL32), pointer :: p_acc_chunk_3d(:,:),p_acc_chunk_2d(:) integer :: is, ie, nx integer :: lm integer :: rank integer :: status integer :: j, k, ig integer, allocatable :: ix(:) type(ESMF_VM) :: vm integer :: mypet, petcount, mpic, iroot integer :: na, nb, nx_sum, nsend integer, allocatable :: RecvCount(:), displs(:) integer :: i, ierr integer :: nsend_v integer, allocatable :: recvcount_v(:), displs_v(:) integer :: ip, M, N if (.NOT. this%active) then _RETURN(ESMF_SUCCESS) endif if (this%nobs_epoch_sum==0) then rc=0 return endif lgr => logging%get_logger('HISTORY.sampler') is=1 do k = 1, this%nobs_type !-- limit nx < 2**32 (integer*4) nx=this%obs(k)%nobs_epoch if (nx >0) then if (mapl_am_i_root()) then call this%obs(k)%file_handle%put_var(this%var_name_time, real(this%obs(k)%times_R8), & start=[is], count=[nx], _RC) call this%obs(k)%file_handle%put_var(this%var_name_lon, this%obs(k)%lons, & start=[is], count=[nx], _RC) call this%obs(k)%file_handle%put_var(this%var_name_lat, this%obs(k)%lats, & start=[is], count=[nx], _RC) call this%obs(k)%file_handle%put_var(this%location_index_name, this%obs(k)%location_index_ioda, & start=[is], count=[nx], _RC) end if end if enddo ! get RH from 2d field src_field = ESMF_FieldCreate(this%LS_ds,typekind=ESMF_TYPEKIND_R4,gridToFieldMap=[1],_RC) chunk_field = ESMF_FieldCreate(this%LS_chunk,typekind=ESMF_TYPEKIND_R4,gridToFieldMap=[1],_RC) call ESMF_FieldGet( src_field, localDE=0, farrayPtr=pt1, _RC ) call ESMF_FieldGet( chunk_field, localDE=0, farrayPtr=pt2, _RC ) pt1=0.0 pt2=0.0 call ESMF_FieldRedistStore(src_field,chunk_field,RH,_RC) call ESMF_FieldDestroy(src_field,noGarbage=.true.,_RC) call ESMF_FieldDestroy(chunk_field,noGarbage=.true.,_RC) ! redist and put_var lm = this%vdata%lm acc_field_2d_rt = ESMF_FieldCreate (this%LS_rt, name='field_2d_rt', typekind=ESMF_TYPEKIND_R4, _RC) acc_field_3d_rt = ESMF_FieldCreate (this%LS_rt, name='field_3d_rt', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) acc_field_2d_chunk = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) acc_field_3d_chunk = ESMF_FieldCreate (this%LS_chunk, name='field_3d_chunk', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) ! ! caution about zero-sized array for MPI ! nx_sum = this%nobs_epoch_sum call ESMF_VMGetCurrent(vm,_RC) call ESMF_VMGet(vm, mpiCommunicator=mpic, petCount=petCount, localPet=mypet, _RC) iroot = 0 ip = mypet N = nx_sum M = petCount nsend = int(ip+1, INT64) * int(N, INT64) / int(M, INT64) - & int(ip , INT64) * int(N, INT64) / int(M, INT64) allocate ( recvcount (petCount) ) allocate ( displs (petCount) ) do ip=0, M-1 recvcount(ip+1) = int(ip+1, INT64) * int(N, INT64) / int(M, INT64) - & int(ip , INT64) * int(N, INT64) / int(M, INT64) end do displs(1)=0 do i = 2, petCount displs(i) = displs(i-1) + recvcount(i-1) end do nsend_v = nsend * lm ! vertical allocate (recvcount_v, source = recvcount * lm ) allocate (displs_v, source = displs * lm ) if (mapl_am_i_root()) then allocate ( p_acc_rt_2d(nx_sum) ) else allocate ( p_acc_rt_2d(1) ) end if ! ! p_dst (lm, nx) if (mapl_am_i_root()) then allocate ( p_acc_rt_3d(nx_sum,lm) ) allocate ( p_dst_rt(lm, nx_sum) ) else allocate ( p_acc_rt_3d(1,lm) ) allocate ( p_dst_rt(lm, 1) ) end if iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() if (item%itemType == ItemTypeScalar) then call ESMF_FieldBundleGet(this%acc_bundle,trim(item%xname),field=acc_field,_RC) call ESMF_FieldGet(acc_field,rank=rank,_RC) if (rank==1) then call ESMF_FieldGet( acc_field, localDE=0, farrayPtr=p_acc_2d, _RC ) call ESMF_FieldGet( acc_field_2d_chunk, localDE=0, farrayPtr=p_acc_chunk_2d, _RC ) call ESMF_FieldRedist( acc_field, acc_field_2d_chunk, RH, _RC ) call MPI_gatherv ( p_acc_chunk_2d, nsend, MPI_REAL, & p_acc_rt_2d, recvcount, displs, MPI_REAL,& iroot, mpic, ierr ) _VERIFY(ierr) if (mapl_am_i_root()) then ! !-- pack fields to obs(k)%p2d and put_var ! is=1 ie=this%epoch_index(2)-this%epoch_index(1)+1 do k=1, this%nobs_type nx = this%obs(k)%nobs_epoch allocate (this%obs(k)%p2d(nx), _STAT) enddo allocate(ix(this%nobs_type), _STAT) ix(:)=0 do j=is, ie k = this%obstype_id(j) ix(k) = ix(k) + 1 this%obs(k)%p2d(ix(k)) = p_acc_rt_2d(j) enddo do k=1, this%nobs_type if (ix(k) /= this%obs(k)%nobs_epoch) then print*, 'obs_', k, ' : ix(k) /= this%obs(k)%nobs_epoch' print*, 'obs_', k, ' : this%obs(k)%nobs_epoch, ix(k) =', this%obs(k)%nobs_epoch, ix(k) _FAIL('test ix(k) failed') endif enddo deallocate(ix, _STAT) do k=1, this%nobs_type is = 1 nx = this%obs(k)%nobs_epoch if (nx>0) then do ig = 1, this%obs(k)%ngeoval if (trim(item%xname) == trim(this%obs(k)%geoval_xname(ig))) then call this%obs(k)%file_handle%put_var(trim(item%xname), this%obs(k)%p2d(1:nx), & start=[is],count=[nx]) end if enddo endif enddo do k=1, this%nobs_type deallocate (this%obs(k)%p2d, _STAT) enddo end if else if (rank==2) then call ESMF_FieldGet( acc_field, localDE=0, farrayPtr=p_acc_3d, _RC) dst_field=ESMF_FieldCreate(this%LS_chunk,typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) src_field=ESMF_FieldCreate(this%LS_ds,typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) call ESMF_FieldGet(src_field,localDE=0,farrayPtr=p_src,_RC) call ESMF_FieldGet(dst_field,localDE=0,farrayPtr=p_dst,_RC) p_src= reshape(p_acc_3d,shape(p_src), order=[2,1]) call ESMF_FieldRegrid(src_field,dst_field,RH,_RC) if (this%level_by_level) then ! p_dst (lm, nx) allocate ( p_dst_t, source = reshape ( p_dst, [size(p_dst,2),size(p_dst,1)], order=[2,1] ) ) do k = 1, lm call MPI_gatherv ( p_dst_t(1,k), nsend, MPI_REAL, & p_acc_rt_3d(1,k), recvcount, displs, MPI_REAL,& iroot, mpic, ierr ) _VERIFY(ierr) end do deallocate (p_dst_t) else call MPI_gatherv ( p_dst, nsend_v, MPI_REAL, & p_dst_rt, recvcount_v, displs_v, MPI_REAL,& iroot, mpic, ierr ) _VERIFY(ierr) p_acc_rt_3d = reshape ( p_dst_rt, shape(p_acc_rt_3d), order=[2,1] ) end if call ESMF_FieldDestroy(dst_field,noGarbage=.true.,_RC) call ESMF_FieldDestroy(src_field,noGarbage=.true.,_RC) if (mapl_am_i_root()) then ! !-- pack fields to obs(k)%p3d and put_var ! is=1 ie=this%epoch_index(2)-this%epoch_index(1)+1 do k=1, this%nobs_type nx = this%obs(k)%nobs_epoch allocate (this%obs(k)%p3d(nx, size(p_acc_rt_3d,2)), _STAT) enddo allocate(ix(this%nobs_type), _STAT) ix(:)=0 do j=is, ie k = this%obstype_id(j) ix(k) = ix(k) + 1 this%obs(k)%p3d(ix(k),:) = p_acc_rt_3d(j,:) enddo deallocate(ix, _STAT) do k=1, this%nobs_type is = 1 nx = this%obs(k)%nobs_epoch if (nx>0) then do ig = 1, this%obs(k)%ngeoval if (trim(item%xname) == trim(this%obs(k)%geoval_xname(ig))) then call this%obs(k)%file_handle%put_var(trim(item%xname), this%obs(k)%p3d(:,:), & start=[is,1],count=[nx,size(p_acc_rt_3d,2)]) end if end do endif enddo do k=1, this%nobs_type deallocate (this%obs(k)%p3d, _STAT) enddo end if endif else if (item%itemType == ItemTypeVector) then _FAIL("ItemTypeVector not yet supported") end if call iter%next() enddo call ESMF_FieldDestroy(acc_field_2d_chunk, noGarbage=.true., _RC) call ESMF_FieldDestroy(acc_field_3d_chunk, noGarbage=.true., _RC) call ESMF_FieldRedistRelease(RH, noGarbage=.true., _RC) _RETURN(_SUCCESS) end procedure append_file module procedure regrid_accumulate_on_xsubset integer :: x_subset(2) type(ESMF_Time) :: timeset(2) type(ESMF_Time) :: current_time type(ESMF_TimeInterval) :: dur type(GriddedIOitemVectorIterator) :: iter type(GriddedIOitem), pointer :: item type(ESMF_Field) :: src_field,dst_field,acc_field integer :: rank real(REAL32), allocatable :: p_new_lev(:,:,:) real(REAL32), pointer :: p_src_3d(:,:,:),p_src_2d(:,:) real(REAL32), pointer :: p_dst_3d(:,:),p_dst_2d(:) real(REAL32), pointer :: p_acc_3d(:,:),p_acc_2d(:) type(ESMF_VM) :: vm integer :: mypet, petcount integer :: is, ie, nx_sum integer :: status integer :: arr(1) if (.NOT. this%active) then _RETURN(ESMF_SUCCESS) endif if (this%nobs_epoch_sum==0) then _RETURN(ESMF_SUCCESS) endif if (this%nobs_epoch_sum==0) then rc=0 return endif call ESMF_ClockGet(this%clock,currTime=current_time,_RC) call ESMF_ClockGet(this%clock,timeStep=dur, _RC ) timeset(1) = current_time - dur timeset(2) = current_time call this%get_x_subset(timeset, x_subset, _RC) is=x_subset(1) ie=x_subset(2) !! write(6,'(2x,a,4i10)') 'in regrid_accumulate is, ie=', is, ie ! ! __ I designed a method to return from regridding if no valid points exist ! in reality for 29 ioda platforms and dt > 20 sec, we donot need this ! !!arr(1)=1 !!if (.NOT. (is > 0 .AND. is <= ie )) arr(1)=0 !!call ESMF_VMGetCurrent(vm,_RC) !!call ESMF_VMGet(vm, localPet=mypet, petCount=petCount, _RC) !!call ESMF_VMAllFullReduce(vm, sendData=arr, recvData=nx_sum, & !! count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) !!if ( nx_sum == 0 ) then !! write(6, '(2x,a,2x,3i10)') 'invalid points, mypet, is, ie =', mypet, is, ie !! ! no valid points to regrid !! _RETURN(ESMF_SUCCESS) !!else !! write(6, '(2x,a,2x,3i10)') ' valid points, mypet, is, ie =', mypet, is, ie !!end if if (this%vdata%regrid_type==VERTICAL_METHOD_ETA2LEV) then call this%vdata%setup_eta_to_pressure(_RC) endif iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() if (item%itemType == ItemTypeScalar) then call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) call ESMF_FieldBundleGet(this%output_bundle,trim(item%xname),field=dst_field,_RC) call ESMF_FieldBundleGet(this%acc_bundle,trim(item%xname),field=acc_field,_RC) call ESMF_FieldGet(src_field,rank=rank,_RC) if (rank==2) then call ESMF_FieldGet(src_field,farrayptr=p_src_2d,_RC) call ESMF_FieldGet(dst_field,farrayptr=p_dst_2d,_RC) call ESMF_FieldGet(acc_field,farrayptr=p_acc_2d,_RC) !! print*, 'size(src,dst,acc)', size(p_src_2d), size(p_dst_2d), size(p_acc_2d) call this%regridder%regrid(p_src_2d,p_dst_2d,_RC) if (is > 0 .AND. is <= ie ) then p_acc_2d(is:ie) = p_dst_2d(is:ie) endif !!if (is>0) write(6,'(a)') 'regrid_accu: p_dst_2d' !!if (is>0) write(6,'(10f7.1)') p_dst_2d else if (rank==3) then call ESMF_FieldGet(src_field,farrayptr=p_src_3d,_RC) call ESMF_FieldGet(dst_field,farrayptr=p_dst_3d,_RC) call ESMF_FieldGet(acc_field,farrayptr=p_acc_3d,_RC) if (this%vdata%regrid_type==VERTICAL_METHOD_ETA2LEV) then allocate(p_new_lev(size(p_src_3d,1),size(p_src_3d,2),this%vdata%lm),_STAT) call this%vdata%regrid_eta_to_pressure(p_src_3d,p_new_lev,_RC) call this%regridder%regrid(p_new_lev,p_dst_3d,_RC) if (is > 0 .AND. is <= ie ) then p_acc_3d(is:ie,:) = p_dst_3d(is:ie,:) end if else call this%regridder%regrid(p_src_3d,p_dst_3d,_RC) if (is > 0 .AND. is <= ie ) then p_acc_3d(is:ie,:) = p_dst_3d(is:ie,:) end if end if end if else if (item%itemType == ItemTypeVector) then _FAIL("ItemTypeVector not yet supported") end if call iter%next() enddo _RETURN(ESMF_SUCCESS) end procedure regrid_accumulate_on_xsubset module procedure destroy_rh_regen_LS integer :: status integer :: numVars, i, k character(len=ESMF_MAXSTR), allocatable :: names(:) type(ESMF_Field) :: field type(ESMF_Time) :: currTime if (.NOT. this%active) then _RETURN(ESMF_SUCCESS) endif call ESMF_FieldDestroy(this%fieldB,nogarbage=.true.,_RC) call this%locstream_factory%destroy_locstream(this%LS_rt, _RC) call this%locstream_factory%destroy_locstream(this%LS_ds, _RC) call this%regridder%destroy(_RC) deallocate (this%lons, this%lats, & this%times_R8, this%obstype_id, this%location_index_ioda, _STAT) do k=1, this%nobs_type deallocate (this%obs(k)%metadata, _STAT) if (mapl_am_i_root()) then deallocate (this%obs(k)%file_handle, _STAT) end if end do if (mapl_am_i_root()) then do k=1, this%nobs_type if (allocated (this%obs(k)%lons)) then deallocate (this%obs(k)%lons) end if if (allocated (this%obs(k)%lats)) then deallocate (this%obs(k)%lats) end if if (allocated (this%obs(k)%times_R8)) then deallocate (this%obs(k)%times_R8) end if if (allocated (this%obs(k)%location_index_ioda)) then deallocate (this%obs(k)%location_index_ioda) end if if (allocated(this%obs(k)%p2d)) then deallocate (this%obs(k)%p2d) endif if (allocated(this%obs(k)%p3d)) then deallocate (this%obs(k)%p3d) endif end do end if call ESMF_FieldBundleGet(this%acc_bundle,fieldCount=numVars,_RC) allocate(names(numVars), _STAT) call ESMF_FieldBundleGet(this%acc_bundle,fieldNameList=names,_RC) do i=1,numVars call ESMF_FieldBundleGet(this%acc_bundle,trim(names(i)),field=field,_RC) call ESMF_FieldDestroy(field,noGarbage=.true., _RC) enddo call ESMF_FieldBundleDestroy(this%acc_bundle,noGarbage=.true.,_RC) deallocate(names, _STAT) call ESMF_FieldBundleGet(this%output_bundle,fieldCount=numVars,_RC) allocate(names(numVars), _STAT) call ESMF_FieldBundleGet(this%output_bundle,fieldNameList=names,_RC) do i=1,numVars call ESMF_FieldBundleGet(this%output_bundle,trim(names(i)),field=field,_RC) call ESMF_FieldDestroy(field,noGarbage=.true., _RC) enddo call ESMF_FieldBundleDestroy(this%output_bundle,noGarbage=.true.,_RC) deallocate(names, _STAT) call ESMF_ClockGet ( this%clock, CurrTime=currTime, _RC ) if (currTime > this%obsfile_end_time) then this%active = .false. _RETURN(ESMF_SUCCESS) end if this%epoch_index(1:2)=0 call this%initialize(reinitialize=.true., _RC) _RETURN(ESMF_SUCCESS) end procedure destroy_rh_regen_LS module procedure get_x_subset type (ESMF_Time) :: T1, T2 real (ESMF_KIND_R8) :: rT1, rT2 integer(ESMF_KIND_I8) :: i1, i2 integer(ESMF_KIND_I8) :: index1, index2, lb, ub integer :: jlo, jhi integer :: status T1= interval(1) T2= interval(2) call time_esmf_2_nc_int (T1, this%datetime_units, i1, _RC) call time_esmf_2_nc_int (T2, this%datetime_units, i2, _RC) rT1=real(i1, kind=ESMF_KIND_R8) rT2=real(i2, kind=ESMF_KIND_R8) jlo = 1 !! !! I choose UB = N+1 not N, because my sub. bisect find n: Y(n)<x<=Y(n+1) jhi= size(this%obstime) + 1 !write(6,*) 'size(this%obstime)=', size(this%obstime) if (jhi==1) then x_subset(1:2)=0 _RETURN(_SUCCESS) endif lb=int(jlo, ESMF_KIND_I8) ub=int(jhi, ESMF_KIND_I8) call bisect( this%obstime, rT1, index1, n_LB=lb, n_UB=ub, rc=rc) call bisect( this%obstime, rT2, index2, n_LB=lb, n_UB=ub, rc=rc) ! (x1, x2] design in bisect ! simple version x_subset(1) = index1+1 x_subset(2) = index2 !! write(6,'(2x,a,2x,2i10)') 'mod vers. get_x_subset, LB,UB=', x_subset(1:2) _RETURN(_SUCCESS) end procedure get_x_subset function extract_unquoted_item(string_list) result(item) character(:), allocatable :: item character(*), intent(in) :: string_list integer :: i integer :: j character(1) :: QUOTE = "'" i = index(string_list( 1:), QUOTE) j = index(string_list(i+1:), QUOTE)+i if( i.ne.0 ) then item = adjustl( string_list(i+1:j-1) ) else item = adjustl( string_list) endif end function extract_unquoted_item end submodule HistoryTrajectory_implement