From 3972d46fa6105af2505e69ed858df8eb7e4c6b1c Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Thu, 16 Jan 2025 08:03:11 -0700 Subject: [PATCH] Trajectory sampler: support splitfield in HISTORY.rc (#3329) --- CHANGELOG.md | 1 + base/MAPL_ObsUtil.F90 | 11 ++-- base/Plain_netCDF_Time.F90 | 49 +++++++++++++++ gridcomps/History/MAPL_HistoryGridComp.F90 | 31 +++++----- .../History/Sampler/MAPL_TrajectoryMod.F90 | 1 + .../Sampler/MAPL_TrajectoryMod_smod.F90 | 59 ++++++++++++++++--- 6 files changed, 125 insertions(+), 27 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f0c4478cdc38..a6b309c9ffa0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added loggers when writing or reading weight files - Added new option to AGCM.rc `overwrite_checkpoint` to allow checkpoint files to be overwritten. By default still will not overwrite checkpoints - The trajectory sampler netCDF output variable `location_index_in_iodafile` can be turned off, after we add two control variables: `use_NWP_1_file` and `restore_2_obs_vector` for users. When set to true, the two options will select only one obs file at each Epoch interval, and will rotate the output field index back to the location vector inthe obs file before generating netCDF output. +- Support `splitfield: 1` in HISTORY.rc for trajectory sampler ### Changed diff --git a/base/MAPL_ObsUtil.F90 b/base/MAPL_ObsUtil.F90 index 16df02eebd38..aed4c9adcd59 100644 --- a/base/MAPL_ObsUtil.F90 +++ b/base/MAPL_ObsUtil.F90 @@ -53,7 +53,7 @@ module MAPL_ObsUtilMod character (len=ESMF_MAXSTR) :: var_name_time='' character (len=ESMF_MAXSTR) :: file_name_template='' integer :: ngeoval=0 - integer :: nentry_name=0 + integer :: nfield_name_mx=12 character (len=ESMF_MAXSTR), allocatable :: field_name(:,:) !character (len=ESMF_MAXSTR), allocatable :: field_name(:) end type obs_platform @@ -771,7 +771,7 @@ function copy_platform_nckeys(a, rc) copy_platform_nckeys%var_name_lon = a%var_name_lon copy_platform_nckeys%var_name_lat = a%var_name_lat copy_platform_nckeys%var_name_time = a%var_name_time - copy_platform_nckeys%nentry_name = a%nentry_name + copy_platform_nckeys%nfield_name_mx = a%nfield_name_mx _RETURN(_SUCCESS) end function copy_platform_nckeys @@ -784,7 +784,7 @@ function union_platform(a, b, rc) integer, optional, intent(out) :: rc character (len=ESMF_MAXSTR), allocatable :: field_name_loc(:,:) - integer :: nfield, nentry_name + integer :: nfield, nfield_name_mx integer, allocatable :: tag(:) integer :: i, j, k integer :: status @@ -805,9 +805,9 @@ function union_platform(a, b, rc) enddo union_platform%ngeoval=k nfield=k - nentry_name=union_platform%nentry_name + nfield_name_mx=union_platform%nfield_name_mx if ( allocated (union_platform%field_name) ) deallocate(union_platform%field_name) - allocate(union_platform%field_name(nentry_name, nfield)) + allocate(union_platform%field_name(nfield_name_mx, nfield)) do i=1, a%ngeoval union_platform%field_name(:,i) = a%field_name(:,i) enddo @@ -953,6 +953,7 @@ subroutine fglob(search_name, filename, rc) ! give the last name if (lenmax < slen) then if (MAPL_AM_I_ROOT()) write(6,*) 'pathlen vs filename_max_char_len: ', slen, lenmax _FAIL ('PATHLEN is greater than filename_max_char_len') + STOP 'lenmax < slen' end if if (slen>0) filename(1:slen)=c_filename(1:slen) diff --git a/base/Plain_netCDF_Time.F90 b/base/Plain_netCDF_Time.F90 index b9c163816647..924ba9323978 100644 --- a/base/Plain_netCDF_Time.F90 +++ b/base/Plain_netCDF_Time.F90 @@ -840,6 +840,13 @@ subroutine split_string_by_space (string_in, length_mx, & wc=0 ios=0 string = trim( adjustl(string_in) ) + str_piece(:)='' + i = index (string, mark) + if (i==0) then + nseg=1 + str_piece(1)=string + return + end if do while (ios==0) i = index (string, mark) !!print*, 'index=', i @@ -858,5 +865,47 @@ subroutine split_string_by_space (string_in, length_mx, & end subroutine split_string_by_space + subroutine split_string_by_seperator (string_in, length_mx, seperator_in, & + mxseg, nseg, str_piece, jstatus) + character (len=length_mx), intent (in) :: string_in + integer, intent (in) :: length_mx + character (len=length_mx), intent (in) :: seperator_in + integer, intent (in) :: mxseg + integer, intent (out):: nseg + character (len=length_mx), intent (out):: str_piece(mxseg) + integer, intent (out):: jstatus + character (len=length_mx) :: string_sc, string_oper, string_aux + character (len=1) :: mark, CH + integer :: ios + integer :: wc + integer :: len1, len2 + ! + ! "xxxx; yy; zz; uu, vv," + ! seperator = ";," + ! + + !__ s1. replace seperator by space + ! + string_sc = trim( adjustl(string_in) ) + string_oper = trim( adjustl(seperator_in) ) + len1 = len_trim(string_sc) + len2 = len_trim(string_oper) + string_aux=string_sc + do i = 1, len1 + CH = string_sc(i:i) + do j = 1, len2 + mark = string_oper(j:j) + if (CH==mark) then + string_aux(i:i)=' ' + end if + end do + end do + + !__ s2. split by space + call split_string_by_space (string_aux, length_mx, & + mxseg, nseg, str_piece, jstatus) + + return + end subroutine split_string_by_seperator end module MAPL_scan_pattern_in_file diff --git a/gridcomps/History/MAPL_HistoryGridComp.F90 b/gridcomps/History/MAPL_HistoryGridComp.F90 index bce8c9b46a44..a6e95338cd1e 100644 --- a/gridcomps/History/MAPL_HistoryGridComp.F90 +++ b/gridcomps/History/MAPL_HistoryGridComp.F90 @@ -5464,7 +5464,7 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) integer :: nseg integer :: nseg_ub integer :: nfield, nplatform - integer :: nentry_name + integer :: nfield_name_max logical :: obs_flag integer, allocatable :: map(:) type(Logger), pointer :: lgr @@ -5554,7 +5554,7 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) - ! __ s2.1 scan fields: only determine ngeoval / nentry_name = nword + ! __ s2.1 scan fields: only determine ngeoval / nfield_name_max = nword allocate (str_piece(mxseg)) rewind(unitr) do k=1, nplf @@ -5578,10 +5578,10 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) end if enddo PLFS(k)%ngeoval = ngeoval - PLFS(k)%nentry_name = nseg_ub + nseg_ub = PLFS(k)%nfield_name_mx allocate ( PLFS(k)%field_name (nseg_ub, ngeoval) ) PLFS(k)%field_name = '' - !! print*, 'k, ngeoval, nentry_name', k, ngeoval, nseg_ub + !! print*, 'k, ngeoval, nfield_name_max', k, ngeoval, nseg_ub end do @@ -5679,12 +5679,12 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) call split_string_by_space (line, length_mx, mxseg, & nplatform, str_piece, status) - !! to do: add debug - !write(6,*) 'line for obsplatforms=', trim(line) - !write(6,*) 'split string, nplatform=', nplatform - !write(6,*) 'nplf=', nplf - !write(6,*) 'str_piece=', str_piece(1:nplatform) - + call lgr%debug('%a %a', 'line for obsplatforms=', trim(line)) + call lgr%debug('%a %i6', 'split string, nplatform=', nplatform) + call lgr%debug('%a %i6', 'nplf=', nplf) + !if (mapl_am_i_root()) then + ! write(6,*) ' str_piece=', str_piece(1:nplatform) + !end if ! ! a) union the platform @@ -5700,7 +5700,10 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) end do end do deallocate(str_piece) - !! write(6,*) 'collection n=',n, 'map(:)=', map(:) + !if (mapl_am_i_root()) then + ! write(6,*) 'collection n=',n, 'map(:)=', map(:) + !end if + ! __ write common nc_index,time,lon,lat k=map(1) ! plat form # 1 @@ -5719,10 +5722,10 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) end do nfield = p1%ngeoval - nentry_name = p1%nentry_name + nfield_name_max = p1%nfield_name_mx do j=1, nfield line='' - do i=1, nentry_name + do i=1, nfield_name_max line = trim(line)//' '//trim(p1%field_name(i,j)) enddo if (j==1) then @@ -5741,7 +5744,7 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) write(unitw, '(a)') trim(adjustl(PLFS(k)%file_name_template)) do j=1, PLFS(k)%ngeoval line='' - do i=1, nentry_name + do i=1, nfield_name_max line = trim(line)//' '//trim(adjustl(PLFS(k)%field_name(i,j))) enddo write(unitw, '(a)') trim(adjustl(line)) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 index a8e2050f8e68..7cfccd4b6d40 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 @@ -71,6 +71,7 @@ module HistoryTrajectoryMod integer :: obsfile_Te_index logical :: active ! case: when no obs. exist logical :: level_by_level = .false. + ! ! note ! for MPI_GATHERV of 3D data in procedure :: append_file ! we have choice LEVEL_BY_LEVEL or ALL_AT_ONCE (timing in sec below for extdata) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index 35a7074e1c08..d751474c3483 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -29,22 +29,25 @@ module procedure HistoryTrajectory_from_config use BinIOMod + use MAPL_scan_pattern_in_file 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) :: STR1, line, splitter character(len=ESMF_MAXSTR) :: symd, shms integer :: nline, col integer, allocatable :: ncol(:) character(len=ESMF_MAXSTR), allocatable :: word(:) + character(len=ESMF_MAXSTR), allocatable :: str_piece(:) integer :: nobs, head, jvar logical :: tend integer :: i, j, k, k2, M integer :: count, idx integer :: unitr, unitw + integer :: length_mx, mxseg, nseg type(GriddedIOitem) :: item type(Logger), pointer :: lgr @@ -163,6 +166,9 @@ ! __ s3. retrieve template and geoval, set metadata file_handle lgr => logging%get_logger('HISTORY.sampler') + length_mx = ESMF_MAXSTR + mxseg = 100 + allocate (str_piece(mxseg)) if ( nobs == 0 ) then ! ! treatment-1: @@ -208,11 +214,18 @@ ! 1-item case: file template or one-var ! 2-item : var1 , 'root' case STR1=trim(word(1)) + elseif ( count == 3 ) then + ! the Alias case + the splitField case + ! 3-item : var1 , 'root', var1_alias case + ! 3-item : var1 , 'root', 'TOTEXTTAU470;TOTEXTTAU550;TOTEXTTAU870', + ! 3-item : 'u;v' vector interpolation is not handled + STR1=trim(word(3)) else - ! 3-item : var1 , 'root', var1_alias case STR1=trim(word(3)) + call lgr%debug('%a %i8', 'WARNING: there are more than 3 field_names in platform rcx' ) end if deallocate(word, _STAT) + if ( index(trim(STR1), '-----') == 0 ) then if (head==1 .AND. trim(STR1)/='') then nobs=nobs+1 @@ -221,13 +234,29 @@ head=0 else if (trim(STR1)/='') then - jvar=jvar+1 - idx = index(STR1,";") - if (idx==0) then - traj%obs(nobs)%geoval_xname(jvar) = STR1 + splitter=';,' + call split_string_by_seperator (STR1, length_mx, splitter, mxseg, & + nseg, str_piece, status) + if (count < 3) then + ! case + ! 'var1' + ! 'var1' , 'ROOT' + ! 'u;v' , 'ROOT' + jvar=jvar+1 + if (nseg==1) then + traj%obs(nobs)%geoval_xname(jvar) = STR1 + else + traj%obs(nobs)%geoval_xname(jvar) = trim(str_piece(1)) + traj%obs(nobs)%geoval_yname(jvar) = trim(str_piece(2)) + end if else - traj%obs(nobs)%geoval_xname(jvar) = trim(STR1(1:idx-1)) - traj%obs(nobs)%geoval_yname(jvar) = trim(STR1(idx+1:)) + ! case + ! 'var1' , 'ROOT' , alias + ! 'var1' , 'ROOT' , 'TOTEXTTAU470;TOTEXTTAU550;TOTEXTTAU870,' split_field + do j=1, nseg + jvar=jvar+1 + traj%obs(nobs)%geoval_xname(jvar) = trim(str_piece(j)) + end do end if end if end if @@ -239,6 +268,15 @@ enddo end if + !!if (mapl_am_i_root()) then + !! do k=1, nobs + !! do j=1, traj%obs(k)%ngeoval + !! write(6, '(2x,a,2x,2i10,2x,a)') & + !! 'traj%obs(k)%geoval_xname(j), k, j, xname ', k, j, trim(traj%obs(k)%geoval_xname(j)) + !! end do + !! end do + !!end if + do k=1, traj%nobs_type allocate (traj%obs(k)%metadata, _STAT) @@ -454,6 +492,7 @@ iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() + !!if (mapl_am_I_root()) print*, 'create new bundle, this%items%xname= ', trim(item%xname) 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) @@ -1143,6 +1182,8 @@ iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() + !!if (mapl_am_I_root()) print*, 'regrid: item%xname= ', trim(item%xname) + 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) @@ -1202,6 +1243,8 @@ nx = this%obs(k)%nobs_epoch if (nx>0) then do ig = 1, this%obs(k)%ngeoval + !! print*, 'this%obs(k)%geoval_xname(ig)= ', this%obs(k)%geoval_xname(ig) + 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])