From e778e915203fedd8c390d1e20103bd0c189719c4 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Wed, 15 May 2024 23:02:23 -0600 Subject: [PATCH 01/21] Revision to station sampler --- gridcomps/History/MAPL_HistoryGridComp.F90 | 4 +- .../Sampler/MAPL_StationSamplerMod.F90 | 606 ++++++++++++------ 2 files changed, 418 insertions(+), 192 deletions(-) diff --git a/gridcomps/History/MAPL_HistoryGridComp.F90 b/gridcomps/History/MAPL_HistoryGridComp.F90 index 135b982f1a67..50b17a8d9c0d 100644 --- a/gridcomps/History/MAPL_HistoryGridComp.F90 +++ b/gridcomps/History/MAPL_HistoryGridComp.F90 @@ -2429,8 +2429,8 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) list(n)%mask_sampler = MaskSamplerGeosat(cfg,string,clock,_RC) call list(n)%mask_sampler%initialize(items=list(n)%items,bundle=list(n)%bundle,timeinfo=list(n)%timeInfo,vdata=list(n)%vdata,_RC) elseif (list(n)%sampler_spec == 'station') then - list(n)%station_sampler = StationSampler (trim(list(n)%stationIdFile), nskip_line=list(n)%stationSkipLine, _RC) - call list(n)%station_sampler%add_metadata_route_handle(list(n)%bundle,list(n)%timeInfo,vdata=list(n)%vdata,_RC) + list(n)%station_sampler = StationSampler (list(n)%bundle, trim(list(n)%stationIdFile), nskip_line=list(n)%stationSkipLine, _RC) + call list(n)%station_sampler%add_metadata_route_handle(items=list(n)%items,bundle=list(n)%bundle,timeinfo=list(n)%timeInfo,vdata=list(n)%vdata,_RC) else global_attributes = list(n)%global_atts%define_collection_attributes(_RC) if (index(trim(list(n)%output_grid_label), 'SwathGrid') > 0) then diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index 5257f94c375f..7e73b8a74d71 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -9,6 +9,8 @@ module StationSamplerMod use MAPL_BaseMod use MAPL_CommsMod use MAPL_LocstreamRegridderMod + use MPI, only : MPI_INTEGER, MPI_REAL, MPI_REAL8 + use, intrinsic :: iso_fortran_env, only: INT64 use, intrinsic :: iso_fortran_env, only: REAL32 use, intrinsic :: iso_fortran_env, only: REAL64 use, intrinsic :: iso_c_binding, only: C_NULL_CHAR @@ -19,8 +21,11 @@ module StationSamplerMod type :: StationSampler private type(LocStreamFactory) :: LSF - type(ESMF_LocStream) :: esmf_ls + type(ESMF_LocStream) :: LS_rt + type(ESMF_LocStream) :: LS_chunk + type(ESMF_LocStream) :: LS_ds type(LocstreamRegridder) :: regridder + type(ESMF_RouteHandle) :: RH integer :: nstation integer, allocatable :: station_id(:) character(len=ESMF_MAXSTR), allocatable :: station_name(:) @@ -29,7 +34,7 @@ module StationSamplerMod real(kind=REAL64), allocatable :: lats(:) real(kind=REAL64), allocatable :: elevs(:) type(ESMF_FieldBundle) :: bundle - type(FileMetadata) :: fmd + type(FileMetadata) :: metadata type(NetCDF4_FileFormatter) :: formatter type(VerticalData) :: vdata type(TimeData) :: time_info @@ -42,6 +47,7 @@ module StationSamplerMod procedure :: append_file procedure :: get_file_start_time procedure :: compute_time_for_current + procedure :: create_variable => create_metadata_variable end type StationSampler interface StationSampler @@ -50,14 +56,28 @@ module StationSamplerMod contains - function new_StationSampler_readfile (filename,nskip_line, rc) result(sampler) + function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(sampler) use pflogger, only : Logger, logging implicit none type(StationSampler) :: sampler + type(ESMF_FieldBundle), intent(in) :: bundle character(len=*), intent(in) :: filename - integer, optional, intent(in) :: nskip_line + integer, optional, intent(in) :: nskip_line integer, optional, intent(out) :: rc + type(ESMF_VM) :: vm + integer :: mypet, petcount, mpic + type(ESMF_grid) :: grid + integer, allocatable :: sendcount(:), displs(:) + integer :: recvcount + integer :: is, ie, ierr + integer :: M, N, ip + integer :: arr(1) + integer :: nx, nx2, nx_sum + + real(REAL64), allocatable :: lons_chunk(:) + real(REAL64), allocatable :: lats_chunk(:) + integer :: unit, ios, nstation, status integer :: i, j, k, ncount logical :: con1, con2 @@ -73,138 +93,202 @@ function new_StationSampler_readfile (filename,nskip_line, rc) result(sampler) ! ["name_short lat lon elev name_full"] ! - open(newunit=unit, file=trim(filename), form='formatted', & - access='sequential', status='old', _IOSTAT) - ios=0 - nstation=0 - nskip=0 - if (present(nskip_line)) then - nskip=nskip_line - end if - if (nskip>0) then - do i=1, nskip - read(unit, *) - end do - end if - read(unit, '(a100)', IOSTAT=ios) line - call count_substring(line, ',', ncount, _RC) - con1= (ncount>=2 .AND. ncount<=4).OR.(ncount==0) - _ASSERT(con1, 'string sequence in Aeronet file not supported') - if (ncount==0) then - seq='AFFFA' - elseif (ncount==2) then - seq='AFF' - elseif (ncount==3) then - seq='AFFF' - elseif (ncount==4) then - CH1=line(1:1) - con1= (CH1>='a'.AND.CH1<='z').OR.(CH1>='A'.AND.CH1<='Z') - con2= CH1>='0'.AND.CH1<='9' - if (con1) then - seq='AIFFF' - else - if (con2) then - seq='IAFFF' + if ( MAPL_AM_I_ROOT() ) then + open(newunit=unit, file=trim(filename), form='formatted', & + access='sequential', status='old', _IOSTAT) + ios=0 + nstation=0 + nskip=0 + if (present(nskip_line)) then + nskip=nskip_line + end if + if (nskip>0) then + do i=1, nskip + read(unit, *) + end do + end if + read(unit, '(a100)', IOSTAT=ios) line + call count_substring(line, ',', ncount, _RC) + con1= (ncount>=2 .AND. ncount<=4).OR.(ncount==0) + _ASSERT(con1, 'string sequence in Aeronet file not supported') + if (ncount==0) then + seq='AFFFA' + elseif (ncount==2) then + seq='AFF' + elseif (ncount==3) then + seq='AFFF' + elseif (ncount==4) then + CH1=line(1:1) + con1= (CH1>='a'.AND.CH1<='z').OR.(CH1>='A'.AND.CH1<='Z') + con2= CH1>='0'.AND.CH1<='9' + if (con1) then + seq='AIFFF' else - _ASSERT(.false., 'string sequence in Aeronet file not supported') + if (con2) then + seq='IAFFF' + else + _ASSERT(.false., 'string sequence in Aeronet file not supported') + end if end if end if - end if - rewind(unit) - if (nskip>0) then - do i=1, nskip - read(unit, *) + rewind(unit) + if (nskip>0) then + do i=1, nskip + read(unit, *) + end do + end if + ios=0 + do while (ios==0) + read(unit, '(a100)', IOSTAT=ios) line + if (ios==0) nstation=nstation+1 end do - end if - ios=0 - do while (ios==0) - read(unit, '(a100)', IOSTAT=ios) line - if (ios==0) nstation=nstation+1 - end do - sampler%nstation=nstation - allocate(sampler%station_id(nstation), _STAT) - allocate(sampler%station_name(nstation), _STAT) - allocate(sampler%station_fullname(nstation), _STAT) - allocate(sampler%lons(nstation), _STAT) - allocate(sampler%lats(nstation), _STAT) - allocate(sampler%elevs(nstation), _STAT) - - rewind(unit) - if (nskip>0) then - do i=1, nskip - read(unit, *) + sampler%nstation=nstation + allocate(sampler%station_id(nstation), _STAT) + allocate(sampler%station_name(nstation), _STAT) + allocate(sampler%station_fullname(nstation), _STAT) + allocate(sampler%lons(nstation), _STAT) + allocate(sampler%lats(nstation), _STAT) + allocate(sampler%elevs(nstation), _STAT) + + rewind(unit) + if (nskip>0) then + do i=1, nskip + read(unit, *) + end do + end if + do i=1, nstation + if(seq=='IAFFF') then + read(unit, *) & + sampler%station_id(i), & + sampler%station_name(i), & + sampler%lons(i), & + sampler%lats(i) + elseif(seq=='AIFFF') then + read(unit, *) & + sampler%station_name(i), & + sampler%station_id(i), & + sampler%lons(i), & + sampler%lats(i) + elseif(trim(seq)=='AFF' .OR. trim(seq)=='AFFF') then + !!write(6,*) 'i=', i + line='' + read(unit, '(a100)') line + !!write(6,*) 'line=', trim(line) + call CSV_read_line_with_CH_I_R(line, & + sampler%station_name(i), & + sampler%lons(i), & + sampler%lats(i), _RC) + sampler%station_id(i)=i + elseif(trim(seq)=='AFFFA') then + ! NOAA GHCNd + ! Ex: 'CHM00054511 39.9330 116.2830 55.0 BEIJING GSN 54511' + read(unit, *) & + sampler%station_name(i), & + sampler%lats(i), & + sampler%lons(i) + sampler%station_id(i)=i + backspace(unit) + read(unit, '(a100)', IOSTAT=ios) line + j=index(line, '.', BACK=.true.) + line2=line(j+1:) + k=len(line2) + line='' + do j=1, k + CH1=line2(j:j) + con1= (CH1>='a'.AND.CH1<='z').OR.(CH1>='A'.AND.CH1<='Z') + if (con1) exit + enddo + read(line2(j:k), '(a100)') line + line2=trim(line) + k=len(line2) + line='' + do j=1, k + CH1=line2(j:j) + con1= (CH1>='0' .AND. CH1<='9') + if (con1) exit + enddo + if (j>k) j=k + sampler%station_fullname(i) = trim(line2(1:j-1)) + end if end do + close(unit) + lgr => logging%get_logger('HISTORY.sampler') + call lgr%debug('%a %i8', 'nstation=', nstation) + call lgr%debug('%a %a %a', 'sampler%station_name(1:2) : ', & + trim(sampler%station_name(1)), trim(sampler%station_name(2))) + call lgr%debug('%a %f8.2 %f8.2', 'sampler%lons(1:2) : ',& + sampler%lons(1),sampler%lons(2)) + call lgr%debug('%a %f8.2 %f8.2', 'sampler%lats(1:2) : ',& + sampler%lats(1),sampler%lats(2)) + else + nstation=0 + allocate(sampler%station_id(nstation), _STAT) + allocate(sampler%station_name(nstation), _STAT) + allocate(sampler%station_fullname(nstation), _STAT) + allocate(sampler%lons(nstation), _STAT) + allocate(sampler%lats(nstation), _STAT) + allocate(sampler%elevs(nstation), _STAT) end if - do i=1, nstation - if(seq=='IAFFF') then - read(unit, *) & - sampler%station_id(i), & - sampler%station_name(i), & - sampler%lons(i), & - sampler%lats(i) - elseif(seq=='AIFFF') then - read(unit, *) & - sampler%station_name(i), & - sampler%station_id(i), & - sampler%lons(i), & - sampler%lats(i) - elseif(trim(seq)=='AFF' .OR. trim(seq)=='AFFF') then - !!write(6,*) 'i=', i - line='' - read(unit, '(a100)') line - !!write(6,*) 'line=', trim(line) - call CSV_read_line_with_CH_I_R(line, & - sampler%station_name(i), & - sampler%lons(i), & - sampler%lats(i), _RC) - sampler%station_id(i)=i - elseif(trim(seq)=='AFFFA') then - ! NOAA GHCNd - ! Ex: 'CHM00054511 39.9330 116.2830 55.0 BEIJING GSN 54511' - read(unit, *) & - sampler%station_name(i), & - sampler%lats(i), & - sampler%lons(i) - sampler%station_id(i)=i - backspace(unit) - read(unit, '(a100)', IOSTAT=ios) line - j=index(line, '.', BACK=.true.) - line2=line(j+1:) - k=len(line2) - line='' - do j=1, k - CH1=line2(j:j) - con1= (CH1>='a'.AND.CH1<='z').OR.(CH1>='A'.AND.CH1<='Z') - if (con1) exit - enddo - read(line2(j:k), '(a100)') line - line2=trim(line) - k=len(line2) - line='' - do j=1, k - CH1=line2(j:j) - con1= (CH1>='0' .AND. CH1<='9') - if (con1) exit - enddo - if (j>k) j=k - sampler%station_fullname(i) = trim(line2(1:j-1)) - end if - end do - close(unit) - lgr => logging%get_logger('HISTORY.sampler') - call lgr%debug('%a %i8', 'nstation=', nstation) - call lgr%debug('%a %a %a', 'sampler%station_name(1:2) : ', & - trim(sampler%station_name(1)), trim(sampler%station_name(2))) - call lgr%debug('%a %f8.2 %f8.2', 'sampler%lons(1:2) : ',& - sampler%lons(1),sampler%lons(2)) - call lgr%debug('%a %f8.2 %f8.2', 'sampler%lats(1:2) : ',& - sampler%lats(1),sampler%lats(2)) - - !__ 2. create LocStreamFactory, then esmf_ls including route_handle + + + !__ 2. create LocStreamFactory, then LS_rt including route_handle ! - sampler%LSF = LocStreamFactory(sampler%lons, sampler%lats, _RC) - sampler%esmf_ls = sampler%LSF%create_locstream(_RC) + ! grid_A: LS_rt : on root + ! grid_B: LS_chunk : uniform on cores + ! grid_C: LS_ds : bg=CS + ! + call ESMF_VMGetCurrent(vm,_RC) + call ESMF_VMGet(vm, mpiCommunicator=mpic, petCount=petCount, localPet=mypet, _RC) + call MAPL_CommsBcast(vm, DATA=sampler%nstation, N=1, ROOT=MAPL_Root, _RC) + + nx_sum = sampler%nstation + 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) ) + + 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( sampler%lons, sendcount, & + displs, MPI_REAL8, lons_chunk, & + recvcount, MPI_REAL8, 0, mpic, ierr) + + call MPI_Scatterv( sampler%lats, sendcount, & + displs, MPI_REAL8, lats_chunk, & + recvcount, MPI_REAL8, 0, mpic, ierr) + + ! -- root + sampler%LSF = LocStreamFactory(sampler%lons, sampler%lats, _RC) + sampler%LS_rt = sampler%LSF%create_locstream(_RC) + + ! -- chunk + sampler%LSF = LocStreamFactory(lons_chunk,lats_chunk,_RC) + sampler%LS_chunk = sampler%LSF%create_locstream_on_proc(_RC) + + ! -- distributed + call ESMF_FieldBundleGet(bundle,grid=grid,_RC) + sampler%LS_ds = sampler%LSF%create_locstream_on_proc(grid=grid,_RC) + ! ! init ofile sampler%ofile='' @@ -214,10 +298,11 @@ function new_StationSampler_readfile (filename,nskip_line, rc) result(sampler) end function new_StationSampler_readfile - subroutine add_metadata_route_handle (this,bundle,timeInfo,vdata,rc) + subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) class(StationSampler), intent(inout) :: this - type(ESMF_FieldBundle), intent(in) :: bundle - type(TimeData), intent(inout) :: timeInfo + type(GriddedIOitemVector), optional, intent(inout) :: items + type(ESMF_FieldBundle), optional, intent(inout) :: bundle + type(TimeData), optional, intent(inout) :: timeInfo type(VerticalData), optional, intent(inout) :: vdata integer, optional, intent(out) :: rc @@ -237,6 +322,9 @@ subroutine add_metadata_route_handle (this,bundle,timeInfo,vdata,rc) character(len=ESMF_MAXSTR), allocatable :: fieldNameList(:) character(len=ESMF_MAXSTR) :: var_name, long_name, units, vdims + type(ESMF_Field) :: src_field, chunk_field + real(REAL32), pointer :: pt1(:), pt2(:) + !__ 1. metadata add_dimension, ! add_variable for time, latlon, station ! @@ -247,81 +335,123 @@ subroutine add_metadata_route_handle (this,bundle,timeInfo,vdata,rc) else this%vdata = VerticalData(_RC) end if - call this%vdata%append_vertical_metadata(this%fmd,this%bundle,_RC) ! specify lev in fmd + call this%vdata%append_vertical_metadata(this%metadata,this%bundle,_RC) ! specify lev in fmd do_vertical_regrid = (this%vdata%regrid_type /= VERTICAL_METHOD_NONE) if (this%vdata%regrid_type == VERTICAL_METHOD_ETA2LEV) then call this%vdata%get_interpolating_variable(this%bundle,_RC) endif - call timeInfo%add_time_to_metadata(this%fmd,_RC) ! specify time in fmd + call timeInfo%add_time_to_metadata(this%metadata,_RC) ! specify time in fmd this%time_info = timeInfo - call this%fmd%add_dimension('station_index',nstation) + call this%metadata%add_dimension('station_index',nstation) v = Variable(type=pFIO_REAL32, dimensions='station_index') call v%add_attribute('long_name','longitude') call v%add_attribute('unit','degree_east') - call this%fmd%add_variable('longitude',v) + call this%metadata%add_variable('longitude',v) v = Variable(type=pFIO_REAL32, dimensions='station_index') call v%add_attribute('long_name','latitude') call v%add_attribute('unit','degree_north') - call this%fmd%add_variable('latitude',v) + call this%metadata%add_variable('latitude',v) v = Variable(type=pFIO_INT32, dimensions='station_index') - call this%fmd%add_variable('station_id',v) + call this%metadata%add_variable('station_id',v) v = Variable(type=pFIO_STRING, dimensions='station_index') call v%add_attribute('long_name','station name') - call this%fmd%add_variable('station_name',v) + call this%metadata%add_variable('station_name',v) - !__ 2. filemetadata: extract field from bundle, add_variable + !__ 2. filemetadata: + ! create varible with names in metadata; see create_metadata_variable ! - call ESMF_FieldBundleGet(bundle, fieldCount=fieldCount, _RC) - allocate (fieldNameList(fieldCount), _STAT) - call ESMF_FieldBundleGet(bundle, fieldNameList=fieldNameList, _RC) - do i=1, fieldCount - var_name=trim(fieldNameList(i)) - call ESMF_FieldBundleGet(bundle,var_name,field=field,_RC) - call ESMF_FieldGet(field,rank=field_rank,_RC) - call ESMF_AttributeGet(field,name="LONG_NAME",isPresent=is_present,_RC) - if ( is_present ) then - call ESMF_AttributeGet(field, NAME="LONG_NAME",VALUE=long_name, _RC) - else - long_name = var_name - endif - call ESMF_AttributeGet(field,name="UNITS",isPresent=is_present,_RC) - if ( is_present ) then - call ESMF_AttributeGet(field, NAME="UNITS",VALUE=units, _RC) - else - units = 'unknown' - endif - if (field_rank==2) then - vdims = "station_index,time" - v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=[nstation,1]) - else if (field_rank==3) then - vdims = "lev,station_index,time" - call ESMF_FieldGet(field,ungriddedLBound=lb,ungriddedUBound=ub,_RC) - v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=[ub(1)-lb(1)+1,1,1]) + 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 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/)) - call this%fmd%add_variable(trim(var_name),v,_RC) - end do - deallocate (fieldNameList) + call iter%next() + enddo - !__ 3. locstream route handle + !__ 3. route handle CS --> LS_ds ! call ESMF_FieldBundleGet(bundle,grid=grid,_RC) - this%regridder = LocStreamRegridder(grid,this%esmf_ls,_RC) + this%regridder = LocStreamRegridder(grid,this%LS_ds,_RC) + + !__ 4. route handle LS_ds --> LS_chunk + ! + 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,this%RH,_RC) + call ESMF_FieldDestroy(src_field,noGarbage=.true.,_RC) + call ESMF_FieldDestroy(chunk_field,noGarbage=.true.,_RC) + + _RETURN(_SUCCESS) + end subroutine add_metadata_route_handle + + + subroutine create_metadata_variable(this,vname,rc) + class(HistoryTrajectory), intent(inout) :: this + character(len=*), intent(in) :: vname + integer, optional, intent(out) :: rc + + type(ESMF_Field) :: field + type(variable) :: v + logical :: is_present + integer :: field_rank, status + character(len=ESMF_MAXSTR) :: var_name,long_name,units,vdims + 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_AttributeGet(field,name="LONG_NAME",isPresent=is_present,_RC) + if ( is_present ) then + call ESMF_AttributeGet (FIELD, NAME="LONG_NAME",VALUE=long_name, _RC) + else + long_name = var_name + endif + call ESMF_AttributeGet(field,name="UNITS",isPresent=is_present,_RC) + if ( is_present ) then + call ESMF_AttributeGet (FIELD, NAME="UNITS",VALUE=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 + if (field_rank==2) then + vdims = "station_index,time" + v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=[nstation,1]) + else if (field_rank==3) then + vdims = "lev,station_index,time" + call ESMF_FieldGet(field,ungriddedLBound=lb,ungriddedUBound=ub,_RC) + v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=[ub(1)-lb(1)+1,1,1]) + end if + + 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/)) + call this%metadata%add_variable(trim(var_name),v,_RC) _RETURN(_SUCCESS) - end subroutine add_metadata_route_handle + end subroutine create_metadata_variable + subroutine append_file(this,current_time,rc) @@ -342,8 +472,25 @@ subroutine append_file(this,current_time,rc) integer :: i, rank integer :: nx, nz + type(ESMF_VM) :: vm + integer :: mypet, petcount, mpic + 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 + + this%obs_written=this%obs_written+1 + !__ 1. put_var: time variable ! rtimes = this%compute_time_for_current(current_time,_RC) ! rtimes: seconds since opening file @@ -355,26 +502,75 @@ subroutine append_file(this,current_time,rc) start=[this%obs_written],count=[1],_RC) end if - !__ 2. put_var: ungridded_dim from src to dst [regrid] + + !__ 2. regrid + put_var: ungridded_dim from src to dst [regrid] ! - call ESMF_FieldBundleGet(this%bundle, fieldCount=fieldCount, _RC) - allocate (fieldNameList(fieldCount), _STAT) - call ESMF_FieldBundleGet(this%bundle, fieldNameList=fieldNameList, _RC) - do i=1, fieldCount - xname=trim(fieldNameList(i)) - call ESMF_FieldBundleGet(this%bundle,xname,field=src_field,_RC) + + lm = this%vdata%lm + field_2d_rt = ESMF_FieldCreate (this%LS_rt, name='field_2d_rt', typekind=ESMF_TYPEKIND_R4, _RC) + field_3d_rt = ESMF_FieldCreate (this%LS_rt, name='field_3d_rt', typekind=ESMF_TYPEKIND_R4, & + gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) + + field_2d_chunk = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) + 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 + ! redist + nx_sum = this%nstation + 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 ) + + + 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 call ESMF_FieldGet(src_field,farrayptr=p_src_2d,_RC) - dst_field = ESMF_FieldCreate(this%esmf_ls,name=xname, & + dst_field = ESMF_FieldCreate(this%LS_ds,name=xname, & typekind=ESMF_TYPEKIND_R4,_RC) call ESMF_FieldGet(dst_field,farrayptr=p_dst_2d,_RC) call this%regridder%regrid(p_src_2d,p_dst_2d,_RC) + + call ESMF_FieldGet(field_2d_chunk, localDE=0, farrayPtr=p_chunk_2d, _RC ) + call ESMF_FieldRedist(dst_field, field_2d_chunk, this%RH, _RC ) + call MPI_gatherv ( p_chunk_2d, nsend, MPI_REAL, & + p_rt_2d, recvcount, displs, MPI_REAL,& + iroot, mpic, ierr ) + if (mapl_am_i_root()) then - call this%formatter%put_var(xname,p_dst_2d,& + call this%formatter%put_var(xname,p_rt_2d,& start=[1,this%obs_written],count=[this%nstation,1],_RC) end if + call ESMF_FieldDestroy(dst_field,nogarbage=.true.) + else if (rank==3) then call ESMF_FieldGet(src_field,farrayptr=p_src_3d,_RC) call ESMF_FieldGet(src_field,ungriddedLBound=lb,ungriddedUBound=ub,_RC) @@ -382,10 +578,15 @@ subroutine append_file(this,current_time,rc) lb(1)=1 ub(1)=this%vdata%lm end if - dst_field = ESMF_FieldCreate(this%esmf_ls,name=xname,& + dst_field = ESMF_FieldCreate(this%LS_ds,name=xname,& typekind=ESMF_TYPEKIND_R4,ungriddedLBound=lb,ungriddedUBound=ub,_RC) call ESMF_FieldGet(dst_field,farrayptr=p_dst_3d,_RC) call this%regridder%regrid(p_src_3d,p_dst_3d,_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 (mapl_am_i_root()) then nx=size(p_dst_3d,1); nz=size(p_dst_3d,2); allocate(arr(nz, nx), _STAT) arr=reshape(p_dst_3d,[nz,nx],order=[2,1]) @@ -399,7 +600,31 @@ subroutine append_file(this,current_time,rc) _FAIL('grid2LS regridder: rank > 3 not implemented') end if end do - deallocate (fieldNameList) + + + + call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) + call ESMF_FieldGet(src_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 ) + + + + + + call ESMF_FieldBundleGet(this%bundle, fieldCount=fieldCount, _RC) + allocate (fieldNameList(fieldCount), _STAT) + call ESMF_FieldBundleGet(this%bundle, fieldNameList=fieldNameList, _RC) + do i=1, fieldCount + xname=trim(fieldNameList(i)) + + + deallocate (fieldNameList) _RETURN(_SUCCESS) end subroutine append_file @@ -413,14 +638,14 @@ subroutine create_file_handle(this,filename,rc) this%ofile = trim(filename) v = this%time_info%define_time_variable(_RC) - call this%fmd%modify_variable('time',v,_RC) + call this%metadata%modify_variable('time',v,_RC) this%obs_written = 0 if (.not. mapl_am_I_root()) then _RETURN(_SUCCESS) end if call this%formatter%create(trim(filename),_RC) - call this%formatter%write(this%fmd,_RC) + call this%formatter%write(this%metadata,_RC) call this%formatter%put_var('longitude',this%lons,_RC) call this%formatter%put_var('latitude',this%lats,_RC) call this%formatter%put_var('station_id',this%station_id,_RC) @@ -491,7 +716,7 @@ subroutine get_file_start_time(this,start_time,time_units,rc) integer lastspace,since_pos integer year,month,day,hour,min,sec - var => this%fmd%get_variable('time',_RC) + var => this%metadata%get_variable('time',_RC) attr => var%get_attribute('units') ptimeUnits => attr%get_value() select type(pTimeUnits) @@ -570,6 +795,7 @@ subroutine get_file_start_time(this,start_time,time_units,rc) _RETURN(_SUCCESS) end subroutine get_file_start_time + ! TODO: delete and use system utilities when available Subroutine count_substring (str, t, ncount, rc) character (len=*), intent(in) :: str From a19b9126e7be7f12121d11944de3d597f611f426 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Fri, 17 May 2024 08:27:11 -0600 Subject: [PATCH 02/21] update --- .../Sampler/MAPL_StationSamplerMod.F90 | 110 +++++++++++++++--- 1 file changed, 92 insertions(+), 18 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index 7e73b8a74d71..bdcc72031f68 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -464,7 +464,10 @@ subroutine append_file(this,current_time,rc) integer :: ub(1), lb(1) type(ESMF_Field) :: src_field,dst_field real(kind=REAL32), pointer :: p_src_3d(:,:,:),p_src_2d(:,:) - real(kind=REAL32), pointer :: p_dst_3d(:,:),p_dst_2d(:) + real(kind=REAL32), pointer :: p_ds_3d(:,:),p_ds_2d(:) + real(kind=REAL32), pointer :: p_chunk_3d(:,:),p_chunk_2d(:) + real(kind=REAL32), pointer :: p_rt_3d(:,:),p_rt_2d(:) + real(kind=REAL32), allocatable :: arr(:,:) character(len=ESMF_MAXSTR), allocatable :: fieldNameList(:) character(len=ESMF_MAXSTR) :: xname @@ -494,9 +497,6 @@ subroutine append_file(this,current_time,rc) !__ 1. put_var: time variable ! rtimes = this%compute_time_for_current(current_time,_RC) ! rtimes: seconds since opening file - if (this%vdata%regrid_type==VERTICAL_METHOD_ETA2LEV) then - call this%vdata%setup_eta_to_pressure(_RC) - end if if (mapl_am_i_root()) then call this%formatter%put_var('time',rtimes(1:1),& start=[this%obs_written],count=[1],_RC) @@ -507,14 +507,16 @@ subroutine append_file(this,current_time,rc) ! lm = this%vdata%lm - field_2d_rt = ESMF_FieldCreate (this%LS_rt, name='field_2d_rt', typekind=ESMF_TYPEKIND_R4, _RC) - field_3d_rt = ESMF_FieldCreate (this%LS_rt, name='field_3d_rt', typekind=ESMF_TYPEKIND_R4, & + field_ds_2d = ESMF_FieldCreate (this%LS_rt, name='field_2d_rt', typekind=ESMF_TYPEKIND_R4, _RC) + field_ds_3d = ESMF_FieldCreate (this%LS_rt, name='field_3d_rt', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) - field_2d_chunk = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) - field_3d_chunk = ESMF_FieldCreate (this%LS_chunk, name='field_3d_chunk', typekind=ESMF_TYPEKIND_R4, & + field_chunk_2d = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) + field_chunk_3d = 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 ! redist nx_sum = this%nstation @@ -542,24 +544,43 @@ subroutine append_file(this,current_time,rc) allocate (recvcount_v, source = recvcount * lm ) allocate (displs_v, source = displs * lm ) + if (mapl_am_i_root()) then + allocate ( p_rt_2d(nx_sum) ) + else + allocate ( p_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 + + ! + ! reuse + ! the pointers p_ds_2d, p_ds_3d, in field_ds_2d, field_ds_3d + ! the pointers p_chunk_2d, p_chunk_3d in field_chunk_2d, field_chunk_3d + ! gather to p_rt_2d, p_rt_3d + ! 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 call ESMF_FieldGet(src_field,farrayptr=p_src_2d,_RC) - dst_field = ESMF_FieldCreate(this%LS_ds,name=xname, & - typekind=ESMF_TYPEKIND_R4,_RC) - call ESMF_FieldGet(dst_field,farrayptr=p_dst_2d,_RC) - call this%regridder%regrid(p_src_2d,p_dst_2d,_RC) + call ESMF_FieldGet(field_ds_2d,farrayptr=p_ds_2d,_RC) + call ESMF_FieldGet(field_chunk_2d,farrayPtr=p_chunk_2d,_RC) + call ESMF_FieldGet(field_rt_2d,farrayPtr=p_rt_2d,_RC) - call ESMF_FieldGet(field_2d_chunk, localDE=0, farrayPtr=p_chunk_2d, _RC ) - call ESMF_FieldRedist(dst_field, field_2d_chunk, this%RH, _RC ) + call this%regridder%regrid(p_src_2d,p_ds_2d,_RC) + call ESMF_FieldRedist(field_ds_2d, field_chunk_2d, this%RH, _RC ) call MPI_gatherv ( p_chunk_2d, nsend, MPI_REAL, & p_rt_2d, recvcount, displs, MPI_REAL,& iroot, mpic, ierr ) @@ -569,10 +590,61 @@ subroutine append_file(this,current_time,rc) start=[1,this%obs_written],count=[this%nstation,1],_RC) end if - call ESMF_FieldDestroy(dst_field,nogarbage=.true.) - else if (rank==3) then + ! -- regrid + ! -- LS ds->chunk + ! -- chunk->rt + ! 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 + + + call ESMF_FieldGet(src_field,farrayptr=p_src_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 ) + 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 ) + 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) + + + call ESMF_FieldGet(src_field,ungriddedLBound=lb,ungriddedUBound=ub,_RC) if (this%vdata%lm/=(ub(1)-lb(1)+1)) then lb(1)=1 @@ -596,6 +668,8 @@ subroutine append_file(this,current_time,rc) deallocate(arr) end if call ESMF_FieldDestroy(dst_field,nogarbage=.true.) + + else _FAIL('grid2LS regridder: rank > 3 not implemented') end if From 5b8a27ec16ce20cf794fa0e7249198560f5803f7 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Fri, 17 May 2024 11:15:54 -0600 Subject: [PATCH 03/21] . --- .../Sampler/MAPL_StationSamplerMod.F90 | 176 ++++++++---------- 1 file changed, 81 insertions(+), 95 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index bdcc72031f68..55ebe8935c88 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -489,7 +489,9 @@ subroutine append_file(this,current_time,rc) integer :: recvcount integer :: is, ie, ierr integer :: M, N, ip - + + type(GriddedIOitemVectorIterator) :: iter + type(GriddedIOitem), pointer :: item this%obs_written=this%obs_written+1 @@ -591,106 +593,90 @@ subroutine append_file(this,current_time,rc) end if else if (rank==3) then - ! -- regrid - ! -- LS ds->chunk - ! -- chunk->rt - ! - 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 - - - call ESMF_FieldGet(src_field,farrayptr=p_src_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 ) - 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 ) - 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) - - - - 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=xname,& - typekind=ESMF_TYPEKIND_R4,ungriddedLBound=lb,ungriddedUBound=ub,_RC) - call ESMF_FieldGet(dst_field,farrayptr=p_dst_3d,_RC) - call this%regridder%regrid(p_src_3d,p_dst_3d,_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 (mapl_am_i_root()) then - nx=size(p_dst_3d,1); nz=size(p_dst_3d,2); allocate(arr(nz, nx), _STAT) - arr=reshape(p_dst_3d,[nz,nx],order=[2,1]) - call this%formatter%put_var(xname,arr,& - start=[1,1,this%obs_written],count=[nz,nx,1],_RC) - !note: lev,station,time - deallocate(arr) - end if - call ESMF_FieldDestroy(dst_field,nogarbage=.true.) - +!! ! -- regrid +!! ! -- LS ds->chunk +!! ! -- chunk->rt +!! ! +!! 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 +!! +!! +!! call ESMF_FieldGet(src_field,farrayptr=p_src_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 ) +!! 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 ) +!! 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) +!! +!! +!! +!! 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=xname,& +!! typekind=ESMF_TYPEKIND_R4,ungriddedLBound=lb,ungriddedUBound=ub,_RC) +!! call ESMF_FieldGet(dst_field,farrayptr=p_dst_3d,_RC) +!! call this%regridder%regrid(p_src_3d,p_dst_3d,_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 (mapl_am_i_root()) then +!! nx=size(p_dst_3d,1); nz=size(p_dst_3d,2); allocate(arr(nz, nx), _STAT) +!! arr=reshape(p_dst_3d,[nz,nx],order=[2,1]) +!! call this%formatter%put_var(xname,arr,& +!! start=[1,1,this%obs_written],count=[nz,nx,1],_RC) +!! !note: lev,station,time +!! deallocate(arr) +!! end if +!! call ESMF_FieldDestroy(dst_field,nogarbage=.true.) +!! else _FAIL('grid2LS regridder: rank > 3 not implemented') end if end do - - - call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) - call ESMF_FieldGet(src_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 ) - - - - - call ESMF_FieldBundleGet(this%bundle, fieldCount=fieldCount, _RC) allocate (fieldNameList(fieldCount), _STAT) call ESMF_FieldBundleGet(this%bundle, fieldNameList=fieldNameList, _RC) From 871ded0ed90b4944a64421d6af356e1926434e3a Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Mon, 20 May 2024 14:27:54 -0600 Subject: [PATCH 04/21] verion for 2d regrid and gatherV --- .../Sampler/MAPL_StationSamplerMod.F90 | 278 +++++++++--------- 1 file changed, 147 insertions(+), 131 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index bdcc72031f68..a8ba53d720ee 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -4,6 +4,8 @@ module StationSamplerMod use MAPL_ErrorHandlingMod use LocStreamFactoryMod use pFIO + use MAPL_GriddedIOItemMod + use MAPL_GriddedIOItemVectorMod use MAPL_TimeDataMod use MAPL_VerticalDataMod use MAPL_BaseMod @@ -26,10 +28,14 @@ module StationSamplerMod type(ESMF_LocStream) :: LS_ds type(LocstreamRegridder) :: regridder type(ESMF_RouteHandle) :: RH + type(GriddedIOitemVector) :: items + logical :: do_vertical_regrid + integer :: nstation integer, allocatable :: station_id(:) character(len=ESMF_MAXSTR), allocatable :: station_name(:) character(len=ESMF_MAXSTR), allocatable :: station_fullname(:) + character(len=ESMF_MAXSTR) :: index_name_x real(kind=REAL64), allocatable :: lons(:) real(kind=REAL64), allocatable :: lats(:) real(kind=REAL64), allocatable :: elevs(:) @@ -40,6 +46,8 @@ module StationSamplerMod type(TimeData) :: time_info character(LEN=ESMF_MAXPATHLEN) :: ofile integer :: obs_written + + contains procedure :: add_metadata_route_handle procedure :: create_file_handle @@ -93,6 +101,7 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s ! ["name_short lat lon elev name_full"] ! + lgr => logging%get_logger('HISTORY.sampler') if ( MAPL_AM_I_ROOT() ) then open(newunit=unit, file=trim(filename), form='formatted', & access='sequential', status='old', _IOSTAT) @@ -213,16 +222,17 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s end if end do close(unit) - lgr => logging%get_logger('HISTORY.sampler') - call lgr%debug('%a %i8', 'nstation=', nstation) - call lgr%debug('%a %a %a', 'sampler%station_name(1:2) : ', & - trim(sampler%station_name(1)), trim(sampler%station_name(2))) - call lgr%debug('%a %f8.2 %f8.2', 'sampler%lons(1:2) : ',& - sampler%lons(1),sampler%lons(2)) - call lgr%debug('%a %f8.2 %f8.2', 'sampler%lats(1:2) : ',& - sampler%lats(1),sampler%lats(2)) + + write(6,*) 'nstation=', nstation + write(6,*) 'sampler%station_name(1:2) : ', & + trim(sampler%station_name(1)), trim(sampler%station_name(2)) + write(6,*) 'sampler%lons(1:2) : ',& + sampler%lons(1:2) + write(6,*) 'sampler%lats(1:2) : ',& + sampler%lats(1:2) else nstation=0 + sampler%nstation = 0 allocate(sampler%station_id(nstation), _STAT) allocate(sampler%station_name(nstation), _STAT) allocate(sampler%station_fullname(nstation), _STAT) @@ -230,8 +240,19 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s allocate(sampler%lats(nstation), _STAT) allocate(sampler%elevs(nstation), _STAT) end if + sampler%index_name_x = 'station_index' +! call lgr%debug('%a %i8', 'nstation=', nstation) +! call lgr%debug('%a %a %a', 'sampler%station_name(1:2) : ', & +! trim(sampler%station_name(1)), trim(sampler%station_name(2))) +! call lgr%debug('%a %f8.2 %f8.2', 'sampler%lons(1:2) : ',& +! sampler%lons(1),sampler%lons(2)) +! call lgr%debug('%a %f8.2 %f8.2', 'sampler%lats(1:2) : ',& +! sampler%lats(1),sampler%lats(2)) + + + !__ 2. create LocStreamFactory, then LS_rt including route_handle ! ! grid_A: LS_rt : on root @@ -307,6 +328,8 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) integer, optional, intent(out) :: rc type(variable) :: v + type(GriddedIOitemVectorIterator) :: iter + type(GriddedIOitem), pointer :: item type(ESMF_Grid) :: grid type(ESMF_Field) :: field integer :: fieldCount @@ -325,16 +348,21 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) type(ESMF_Field) :: src_field, chunk_field real(REAL32), pointer :: pt1(:), pt2(:) + !__ 1. metadata add_dimension, ! add_variable for time, latlon, station ! - this%bundle = bundle - nstation = this%nstation + 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 + nstation = this%nstation + + call this%vdata%append_vertical_metadata(this%metadata,this%bundle,_RC) ! specify lev in fmd do_vertical_regrid = (this%vdata%regrid_type /= VERTICAL_METHOD_NONE) if (this%vdata%regrid_type == VERTICAL_METHOD_ETA2LEV) then @@ -369,7 +397,7 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() -!! print*, 'list item%xname', trim(item%xname) + 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 @@ -402,7 +430,7 @@ end subroutine add_metadata_route_handle subroutine create_metadata_variable(this,vname,rc) - class(HistoryTrajectory), intent(inout) :: this + class(StationSampler), intent(inout) :: this character(len=*), intent(in) :: vname integer, optional, intent(out) :: rc @@ -411,8 +439,10 @@ subroutine create_metadata_variable(this,vname,rc) logical :: is_present integer :: field_rank, status character(len=ESMF_MAXSTR) :: var_name,long_name,units,vdims + integer :: rank,lb(1),ub(1) 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_AttributeGet(field,name="LONG_NAME",isPresent=is_present,_RC) @@ -435,7 +465,7 @@ subroutine create_metadata_variable(this,vname,rc) if (field_rank==2) then vdims = "station_index,time" - v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=[nstation,1]) + v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=[this%nstation,1]) else if (field_rank==3) then vdims = "lev,station_index,time" call ESMF_FieldGet(field,ungriddedLBound=lb,ungriddedUBound=ub,_RC) @@ -462,37 +492,42 @@ subroutine append_file(this,current_time,rc) integer :: status integer :: fieldCount integer :: ub(1), lb(1) + type(GriddedIOitemVectorIterator) :: iter + type(GriddedIOitem), pointer :: item type(ESMF_Field) :: src_field,dst_field real(kind=REAL32), pointer :: p_src_3d(:,:,:),p_src_2d(:,:) real(kind=REAL32), pointer :: p_ds_3d(:,:),p_ds_2d(:) real(kind=REAL32), pointer :: p_chunk_3d(:,:),p_chunk_2d(:) real(kind=REAL32), pointer :: p_rt_3d(:,:),p_rt_2d(:) - + real(kind=REAL32), allocatable :: p_new_lev(:,:,:) + real(kind=REAL32), allocatable :: arr(:,:) character(len=ESMF_MAXSTR), allocatable :: fieldNameList(:) character(len=ESMF_MAXSTR) :: xname real(kind=ESMF_KIND_R8), allocatable :: rtimes(:) - integer :: i, rank - integer :: nx, nz + + integer :: rank + integer :: i, nz, lm type(ESMF_VM) :: vm - integer :: mypet, petcount, mpic - integer :: nx, nx_sum - integer :: n0 - integer :: arr(1) + integer :: mypet, petcount, mpic, iroot + integer :: n0, nx, nx2 + integer :: na, nb, nx_sum, nsend, nsend_v + + type(ESMF_Field) :: field_ds_2d, field_ds_3d + type(ESMF_Field) :: field_chunk_2d, field_chunk_3d + type(ESMF_Field) :: field_rt_2d, field_rt_3d + 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, allocatable :: recvcount(:), sendcount(:), displs(:) + integer, allocatable :: recvcount_v(:), displs_v(:) integer :: is, ie, ierr integer :: M, N, ip - this%obs_written=this%obs_written+1 - !__ 1. put_var: time variable ! @@ -507,15 +542,17 @@ subroutine append_file(this,current_time,rc) ! lm = this%vdata%lm - field_ds_2d = ESMF_FieldCreate (this%LS_rt, name='field_2d_rt', typekind=ESMF_TYPEKIND_R4, _RC) - field_ds_3d = ESMF_FieldCreate (this%LS_rt, name='field_3d_rt', typekind=ESMF_TYPEKIND_R4, & + field_ds_2d = ESMF_FieldCreate (this%LS_ds, name='field_2d_ds', typekind=ESMF_TYPEKIND_R4, _RC) + field_ds_3d = ESMF_FieldCreate (this%LS_ds, name='field_3d_ds', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) field_chunk_2d = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) field_chunk_3d = ESMF_FieldCreate (this%LS_chunk, name='field_3d_chunk', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) - + field_rt_2d = ESMF_FieldCreate (this%LS_rt, name='field_2d_rt', typekind=ESMF_TYPEKIND_R4, _RC) + field_rt_3d = ESMF_FieldCreate (this%LS_rt, name='field_3d_rt', typekind=ESMF_TYPEKIND_R4, & + gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) ! caution about zero-sized array for MPI ! redist @@ -549,16 +586,6 @@ subroutine append_file(this,current_time,rc) else allocate ( p_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 - ! ! reuse @@ -596,109 +623,98 @@ subroutine append_file(this,current_time,rc) ! -- chunk->rt ! 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 - - - call ESMF_FieldGet(src_field,farrayptr=p_src_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 ) - 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 ) - 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) - - 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=xname,& - typekind=ESMF_TYPEKIND_R4,ungriddedLBound=lb,ungriddedUBound=ub,_RC) - call ESMF_FieldGet(dst_field,farrayptr=p_dst_3d,_RC) - call this%regridder%regrid(p_src_3d,p_dst_3d,_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 (mapl_am_i_root()) then - nx=size(p_dst_3d,1); nz=size(p_dst_3d,2); allocate(arr(nz, nx), _STAT) - arr=reshape(p_dst_3d,[nz,nx],order=[2,1]) - call this%formatter%put_var(xname,arr,& - start=[1,1,this%obs_written],count=[nz,nx,1],_RC) - !note: lev,station,time - deallocate(arr) - end if - call ESMF_FieldDestroy(dst_field,nogarbage=.true.) +! call ESMF_FieldGet(dst_field,farrayptr=p_dst_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 +! +! +! call ESMF_FieldGet(src_field,farrayptr=p_src_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 ) +! 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 ) +! 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) +! +! +! +! 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=xname,& +! typekind=ESMF_TYPEKIND_R4,ungriddedLBound=lb,ungriddedUBound=ub,_RC) +! call ESMF_FieldGet(dst_field,farrayptr=p_dst_3d,_RC) +! call this%regridder%regrid(p_src_3d,p_dst_3d,_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 (mapl_am_i_root()) then +! nx=size(p_dst_3d,1); nz=size(p_dst_3d,2); allocate(arr(nz, nx), _STAT) +! arr=reshape(p_dst_3d,[nz,nx],order=[2,1]) +! call this%formatter%put_var(xname,arr,& +! start=[1,1,this%obs_written],count=[nz,nx,1],_RC) +! !note: lev,station,time +! deallocate(arr) +! end if +! call ESMF_FieldDestroy(dst_field,nogarbage=.true.) else _FAIL('grid2LS regridder: rank > 3 not implemented') end if - end do - - - - call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) - call ESMF_FieldGet(src_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 ) - - - - + endif + end do - call ESMF_FieldBundleGet(this%bundle, fieldCount=fieldCount, _RC) - allocate (fieldNameList(fieldCount), _STAT) - call ESMF_FieldBundleGet(this%bundle, fieldNameList=fieldNameList, _RC) - do i=1, fieldCount - xname=trim(fieldNameList(i)) +!!! call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) +! call ESMF_FieldBundleGet(this%bundle, fieldCount=fieldCount, _RC) +! allocate (fieldNameList(fieldCount), _STAT) +! call ESMF_FieldBundleGet(this%bundle, fieldNameList=fieldNameList, _RC) +! do i=1, fieldCount +! xname=trim(fieldNameList(i)) +! deallocate (fieldNameList) +! end do +! - deallocate (fieldNameList) _RETURN(_SUCCESS) end subroutine append_file From a69347f5f6d317240bda6962014903aaeb0ceafd Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Mon, 20 May 2024 18:06:36 -0600 Subject: [PATCH 05/21] . --- gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index a8ba53d720ee..c7171b07f252 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -507,7 +507,7 @@ subroutine append_file(this,current_time,rc) real(kind=ESMF_KIND_R8), allocatable :: rtimes(:) integer :: rank - integer :: i, nz, lm + integer :: i, j, nz, lm type(ESMF_VM) :: vm integer :: mypet, petcount, mpic, iroot @@ -613,8 +613,14 @@ subroutine append_file(this,current_time,rc) iroot, mpic, ierr ) if (mapl_am_i_root()) then - call this%formatter%put_var(xname,p_rt_2d,& - start=[1,this%obs_written],count=[this%nstation,1],_RC) + do j=1, nx_sum, 500000 + write(6,*) 'p_rt_2d', p_rt_2d(j) + end do + end if + + if (mapl_am_i_root()) then +! call this%formatter%put_var(xname,p_rt_2d,& +! start=[1,this%obs_written],count=[this%nstation,1],_RC) end if else if (rank==3) then From fe6679da3130fae4a101e2534e59f7a35a6ec900 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Tue, 21 May 2024 10:49:37 -0600 Subject: [PATCH 06/21] A version using mpi gatherv --- .../Sampler/MAPL_StationSamplerMod.F90 | 287 +++++++----------- 1 file changed, 118 insertions(+), 169 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index c0417cab225c..cfbb6e11931c 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -30,7 +30,8 @@ module StationSamplerMod type(ESMF_RouteHandle) :: RH type(GriddedIOitemVector) :: items logical :: do_vertical_regrid - + logical :: level_by_level + integer :: nstation integer, allocatable :: station_id(:) character(len=ESMF_MAXSTR), allocatable :: station_name(:) @@ -46,7 +47,6 @@ module StationSamplerMod type(TimeData) :: time_info character(LEN=ESMF_MAXPATHLEN) :: ofile integer :: obs_written - contains procedure :: add_metadata_route_handle @@ -241,16 +241,6 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s allocate(sampler%elevs(nstation), _STAT) end if sampler%index_name_x = 'station_index' - - -! call lgr%debug('%a %i8', 'nstation=', nstation) -! call lgr%debug('%a %a %a', 'sampler%station_name(1:2) : ', & -! trim(sampler%station_name(1)), trim(sampler%station_name(2))) -! call lgr%debug('%a %f8.2 %f8.2', 'sampler%lons(1:2) : ',& -! sampler%lons(1),sampler%lons(2)) -! call lgr%debug('%a %f8.2 %f8.2', 'sampler%lats(1:2) : ',& -! sampler%lats(1),sampler%lats(2)) - !__ 2. create LocStreamFactory, then LS_rt including route_handle @@ -314,7 +304,8 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s ! init ofile sampler%ofile='' sampler%obs_written=0 - + sampler%level_by_level = .true. + _RETURN(_SUCCESS) end function new_StationSampler_readfile @@ -349,8 +340,8 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) real(REAL32), pointer :: pt1(:), pt2(:) - !__ 1. metadata add_dimension, - ! add_variable for time, latlon, station + !__ 1. filemetadata: + ! add_dimension, add_variable for latlon, station ! if(present(bundle)) this%bundle=bundle if(present(items)) this%items=items @@ -362,7 +353,6 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) end if nstation = this%nstation - call this%vdata%append_vertical_metadata(this%metadata,this%bundle,_RC) ! specify lev in fmd do_vertical_regrid = (this%vdata%regrid_type /= VERTICAL_METHOD_NONE) if (this%vdata%regrid_type == VERTICAL_METHOD_ETA2LEV) then @@ -392,12 +382,12 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) !__ 2. filemetadata: - ! create varible with names in metadata; see create_metadata_variable + ! create varible with names in item%xname; see create_metadata_variable ! iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() - print*, 'list item%xname', trim(item%xname) + 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 @@ -431,8 +421,8 @@ end subroutine add_metadata_route_handle subroutine create_metadata_variable(this,vname,rc) class(StationSampler), intent(inout) :: this - character(len=*), intent(in) :: vname - integer, optional, intent(out) :: rc + character(len=*), intent(in) :: vname + integer, optional, intent(out) :: rc type(ESMF_Field) :: field type(variable) :: v @@ -457,11 +447,13 @@ subroutine create_metadata_variable(this,vname,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 + +! -- in future, replace keyword station by index_name_x as in trajectory sampler +! 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 if (field_rank==2) then vdims = "station_index,time" @@ -495,11 +487,14 @@ subroutine append_file(this,current_time,rc) type(GriddedIOitemVectorIterator) :: iter type(GriddedIOitem), pointer :: item type(ESMF_Field) :: src_field,dst_field - real(kind=REAL32), pointer :: p_src_3d(:,:,:),p_src_2d(:,:) - real(kind=REAL32), pointer :: p_ds_3d(:,:),p_ds_2d(:) - real(kind=REAL32), pointer :: p_chunk_3d(:,:),p_chunk_2d(:) - real(kind=REAL32), pointer :: p_rt_3d(:,:),p_rt_2d(:) + real(kind=REAL32), pointer :: p_src_3d(:,:,:),p_src_2d(:,:) ! source + real(kind=REAL32), pointer :: p_dst_3d(:,:) ! destination + real(kind=REAL32), pointer :: p_ds_3d(:,:),p_ds_2d(:) ! distributed LS + real(kind=REAL32), pointer :: p_chunk_3d(:,:),p_chunk_2d(:) ! chunk LS + real(kind=REAL32), pointer :: p_rt_3d(:,:),p_rt_2d(:) ! root LS + real(kind=REAL32), pointer :: p_rt_3d_aux(:,:) real(kind=REAL32), allocatable :: p_new_lev(:,:,:) + real(kind=REAL32), allocatable :: p_dst_t(:,:) real(kind=REAL32), allocatable :: arr(:,:) character(len=ESMF_MAXSTR), allocatable :: fieldNameList(:) @@ -507,7 +502,7 @@ subroutine append_file(this,current_time,rc) real(kind=ESMF_KIND_R8), allocatable :: rtimes(:) integer :: rank - integer :: i, j, nz, lm + integer :: i, j, k, nz, lm type(ESMF_VM) :: vm integer :: mypet, petcount, mpic, iroot @@ -527,13 +522,12 @@ subroutine append_file(this,current_time,rc) integer :: is, ie, ierr integer :: M, N, ip - type(GriddedIOitemVectorIterator) :: iter - type(GriddedIOitem), pointer :: item this%obs_written=this%obs_written+1 !__ 1. put_var: time variable ! + ! rtimes = this%compute_time_for_current(current_time,_RC) ! rtimes: seconds since opening file if (mapl_am_i_root()) then call this%formatter%put_var('time',rtimes(1:1),& @@ -541,22 +535,16 @@ subroutine append_file(this,current_time,rc) end if - !__ 2. regrid + put_var: ungridded_dim from src to dst [regrid] + !__ 2. regrid + put_var: + ! ungridded_dim from src to dst [regrid] ! - lm = this%vdata%lm - field_ds_2d = ESMF_FieldCreate (this%LS_ds, name='field_2d_ds', typekind=ESMF_TYPEKIND_R4, _RC) - field_ds_3d = ESMF_FieldCreate (this%LS_ds, name='field_3d_ds', typekind=ESMF_TYPEKIND_R4, & - gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) - + field_ds_2d = ESMF_FieldCreate (this%LS_ds, name='field_2d_ds', typekind=ESMF_TYPEKIND_R4, _RC) field_chunk_2d = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) - field_chunk_3d = ESMF_FieldCreate (this%LS_chunk, name='field_3d_chunk', typekind=ESMF_TYPEKIND_R4, & - gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) - - field_rt_2d = ESMF_FieldCreate (this%LS_rt, name='field_2d_rt', typekind=ESMF_TYPEKIND_R4, _RC) - field_rt_3d = ESMF_FieldCreate (this%LS_rt, name='field_3d_rt', typekind=ESMF_TYPEKIND_R4, & + dst_field = ESMF_FieldCreate (this%LS_ds, name='dst_field', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) + ! caution about zero-sized array for MPI ! redist nx_sum = this%nstation @@ -579,7 +567,7 @@ subroutine append_file(this,current_time,rc) 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 ) @@ -589,140 +577,101 @@ subroutine append_file(this,current_time,rc) else allocate ( p_rt_2d(1) ) end if - - ! - ! reuse - ! the pointers p_ds_2d, p_ds_3d, in field_ds_2d, field_ds_3d - ! the pointers p_chunk_2d, p_chunk_3d in field_chunk_2d, field_chunk_3d - ! gather to p_rt_2d, p_rt_3d - ! + + ! p_rt_3d (lm, nx) + if (mapl_am_i_root()) then + allocate ( p_rt_3d(lm, nx_sum) ) + allocate ( p_rt_3d_aux(nx_sum, lm) ) + else + allocate ( p_rt_3d(lm, 1) ) + allocate ( p_rt_3d_aux(1,lm) ) + end if + + 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 - call ESMF_FieldGet(src_field,farrayptr=p_src_2d,_RC) - call ESMF_FieldGet(field_ds_2d,farrayptr=p_ds_2d,_RC) - call ESMF_FieldGet(field_chunk_2d,farrayPtr=p_chunk_2d,_RC) - call ESMF_FieldGet(field_rt_2d,farrayPtr=p_rt_2d,_RC) - - call this%regridder%regrid(p_src_2d,p_ds_2d,_RC) - call ESMF_FieldRedist(field_ds_2d, field_chunk_2d, this%RH, _RC ) - call MPI_gatherv ( p_chunk_2d, nsend, MPI_REAL, & - p_rt_2d, recvcount, displs, MPI_REAL,& - iroot, mpic, ierr ) - - if (mapl_am_i_root()) then - do j=1, nx_sum, 500000 - write(6,*) 'p_rt_2d', p_rt_2d(j) - end do - end if - - if (mapl_am_i_root()) then -! call this%formatter%put_var(xname,p_rt_2d,& -! start=[1,this%obs_written],count=[this%nstation,1],_RC) - end if - - else if (rank==3) then - ! -- regrid - ! -- LS ds->chunk - ! -- chunk->rt - ! - call ESMF_FieldGet(src_field,farrayptr=p_src_3d,_RC) - - -! call ESMF_FieldGet(dst_field,farrayptr=p_dst_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,:) + call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) + call ESMF_FieldGet(src_field,rank=rank,_RC) + if (rank==2) then + call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_2d,_RC) + call ESMF_FieldGet(field_ds_2d,localDE=0,farrayptr=p_ds_2d,_RC) + call ESMF_FieldGet(field_chunk_2d,localDE=0,farrayPtr=p_chunk_2d,_RC) + + call this%regridder%regrid(p_src_2d,p_ds_2d,_RC) + call ESMF_FieldRedist(field_ds_2d, field_chunk_2d, this%RH, _RC ) + call MPI_gatherv ( p_chunk_2d, nsend, MPI_REAL, & + p_rt_2d, recvcount, displs, MPI_REAL,& + iroot, mpic, ierr ) + +! if (mapl_am_i_root()) then +! do j=1, nx_sum, 500000 +! write(6,*) 'p_rt_2d', p_rt_2d(j) +! end do ! end if -! end if -! -! -! call ESMF_FieldGet(src_field,farrayptr=p_src_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 ) -! 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 ) -! 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) -! -! -! -! 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=xname,& -! typekind=ESMF_TYPEKIND_R4,ungriddedLBound=lb,ungriddedUBound=ub,_RC) -! call ESMF_FieldGet(dst_field,farrayptr=p_dst_3d,_RC) -! call this%regridder%regrid(p_src_3d,p_dst_3d,_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 (mapl_am_i_root()) then -! nx=size(p_dst_3d,1); nz=size(p_dst_3d,2); allocate(arr(nz, nx), _STAT) -! arr=reshape(p_dst_3d,[nz,nx],order=[2,1]) -! call this%formatter%put_var(xname,arr,& -! start=[1,1,this%obs_written],count=[nz,nx,1],_RC) -! !note: lev,station,time -! deallocate(arr) -! end if -! call ESMF_FieldDestroy(dst_field,nogarbage=.true.) + if (mapl_am_i_root()) then + call this%formatter%put_var(trim(item%xname),p_rt_2d,& + start=[1,this%obs_written],count=[this%nstation,1],_RC) + end if - else - _FAIL('grid2LS regridder: rank > 3 not implemented') - end if - endif - end do - + else if (rank==3) then + ! -- CS-> LS_ds; ds->chunk; gather + ! + call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_3d,_RC) + call ESMF_FieldGet(dst_field,localDE=0,farrayptr=p_dst_3d,_RC) + call this%regridder%regrid(p_src_3d,p_dst_3d,_RC) + + field_ds_3d = ESMF_FieldCreate (this%LS_ds, name='field_3d_ds', typekind=ESMF_TYPEKIND_R4, & + gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) + field_chunk_3d = ESMF_FieldCreate (this%LS_chunk, name='field_3d_chunk', typekind=ESMF_TYPEKIND_R4, & + gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) + call ESMF_FieldGet(field_ds_3d,localDE=0,farrayPtr=p_ds_3d,_RC) + call ESMF_FieldGet(field_chunk_3d,localDE=0,farrayPtr=p_chunk_3d,_RC) + + ! p_ds_3d(lm, nx) + p_ds_3d = reshape(p_dst_3d, shape(p_ds_3d), order=[2,1]) + call ESMF_FieldRedist(field_ds_3d, field_chunk_3d, this%RH, _RC) + + if (this%level_by_level) then + ! p_chunk_3d (lm, nx) + allocate (p_dst_t, source = reshape(p_chunk_3d, [size(p_chunk_3d,2),size(p_chunk_3d,1)], order=[2,1])) + do k = 1, lm + call MPI_gatherv ( p_dst_t(1,k), nsend, MPI_REAL, & + p_rt_3d_aux(1,k), recvcount, displs, MPI_REAL,& + iroot, mpic, ierr ) + end do + p_rt_3d = reshape(p_rt_3d_aux, shape(p_rt_3d), order=[2,1]) + else + call MPI_gatherv ( p_chunk_3d, nsend_v, MPI_REAL, & + p_rt_3d, recvcount_v, displs_v, MPI_REAL,& + iroot, mpic, ierr ) + end if + call ESMF_FieldDestroy(field_ds_3d,noGarbage=.true.,_RC) + call ESMF_FieldDestroy(field_chunk_3d,noGarbage=.true.,_RC) -!!! call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) -! call ESMF_FieldBundleGet(this%bundle, fieldCount=fieldCount, _RC) -! allocate (fieldNameList(fieldCount), _STAT) -! call ESMF_FieldBundleGet(this%bundle, fieldNameList=fieldNameList, _RC) -! do i=1, fieldCount -! xname=trim(fieldNameList(i)) -! deallocate (fieldNameList) -! end do -! +! if (mapl_am_i_root()) then +! do j=1, nx_sum, 500 +! write(6,*) 'p_rt_3d', p_rt_3d(:,j) +! end do +! end if +! if (mapl_am_i_root()) write(6,*) 'regrid + gatherV in 3D' + + if (mapl_am_i_root()) then + nz=size(p_rt_3d,1); nx=size(p_rt_3d,2) + call this%formatter%put_var(trim(item%xname),p_rt_3d,& + start=[1,1,this%obs_written],count=[nz,nx,1],_RC) + !note: lev,station,time + end if + else + _FAIL('grid2LS regridder: rank > 3 not implemented') + end if + endif + + call iter%next() + end do _RETURN(_SUCCESS) end subroutine append_file From 8f86b1e0652b49f19efdd97caf118f6712ac7bb5 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Tue, 21 May 2024 10:52:16 -0600 Subject: [PATCH 07/21] clean --- .../Sampler/MAPL_StationSamplerMod.F90 | 64 +++++++++---------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index cfbb6e11931c..4d5c08b78b6b 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -25,13 +25,13 @@ module StationSamplerMod type(LocStreamFactory) :: LSF type(ESMF_LocStream) :: LS_rt type(ESMF_LocStream) :: LS_chunk - type(ESMF_LocStream) :: LS_ds + type(ESMF_LocStream) :: LS_ds type(LocstreamRegridder) :: regridder type(ESMF_RouteHandle) :: RH type(GriddedIOitemVector) :: items logical :: do_vertical_regrid logical :: level_by_level - + integer :: nstation integer, allocatable :: station_id(:) character(len=ESMF_MAXSTR), allocatable :: station_name(:) @@ -47,7 +47,7 @@ module StationSamplerMod type(TimeData) :: time_info character(LEN=ESMF_MAXPATHLEN) :: ofile integer :: obs_written - + contains procedure :: add_metadata_route_handle procedure :: create_file_handle @@ -82,10 +82,10 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s integer :: M, N, ip integer :: arr(1) integer :: nx, nx2, nx_sum - + real(REAL64), allocatable :: lons_chunk(:) real(REAL64), allocatable :: lats_chunk(:) - + integer :: unit, ios, nstation, status integer :: i, j, k, ncount logical :: con1, con2 @@ -102,7 +102,7 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s ! lgr => logging%get_logger('HISTORY.sampler') - if ( MAPL_AM_I_ROOT() ) then + if ( MAPL_AM_I_ROOT() ) then open(newunit=unit, file=trim(filename), form='formatted', & access='sequential', status='old', _IOSTAT) ios=0 @@ -222,7 +222,7 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s end if end do close(unit) - + write(6,*) 'nstation=', nstation write(6,*) 'sampler%station_name(1:2) : ', & trim(sampler%station_name(1)), trim(sampler%station_name(2)) @@ -232,7 +232,7 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s sampler%lats(1:2) else nstation=0 - sampler%nstation = 0 + sampler%nstation = 0 allocate(sampler%station_id(nstation), _STAT) allocate(sampler%station_name(nstation), _STAT) allocate(sampler%station_fullname(nstation), _STAT) @@ -241,7 +241,7 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s allocate(sampler%elevs(nstation), _STAT) end if sampler%index_name_x = 'station_index' - + !__ 2. create LocStreamFactory, then LS_rt including route_handle ! @@ -252,7 +252,7 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s call ESMF_VMGetCurrent(vm,_RC) call ESMF_VMGet(vm, mpiCommunicator=mpic, petCount=petCount, localPet=mypet, _RC) call MAPL_CommsBcast(vm, DATA=sampler%nstation, N=1, ROOT=MAPL_Root, _RC) - + nx_sum = sampler%nstation ip = mypet ! 0 to M-1 N = nx_sum @@ -299,13 +299,13 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s ! -- distributed call ESMF_FieldBundleGet(bundle,grid=grid,_RC) sampler%LS_ds = sampler%LSF%create_locstream_on_proc(grid=grid,_RC) - + ! ! init ofile sampler%ofile='' sampler%obs_written=0 sampler%level_by_level = .true. - + _RETURN(_SUCCESS) end function new_StationSampler_readfile @@ -339,7 +339,7 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) type(ESMF_Field) :: src_field, chunk_field real(REAL32), pointer :: pt1(:), pt2(:) - + !__ 1. filemetadata: ! add_dimension, add_variable for latlon, station ! @@ -404,7 +404,7 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) this%regridder = LocStreamRegridder(grid,this%LS_ds,_RC) !__ 4. route handle LS_ds --> LS_chunk - ! + ! 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 ) @@ -414,7 +414,7 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) call ESMF_FieldRedistStore(src_field,chunk_field,this%RH,_RC) call ESMF_FieldDestroy(src_field,noGarbage=.true.,_RC) call ESMF_FieldDestroy(chunk_field,noGarbage=.true.,_RC) - + _RETURN(_SUCCESS) end subroutine add_metadata_route_handle @@ -423,7 +423,7 @@ subroutine create_metadata_variable(this,vname,rc) class(StationSampler), intent(inout) :: this character(len=*), intent(in) :: vname integer, optional, intent(out) :: rc - + type(ESMF_Field) :: field type(variable) :: v logical :: is_present @@ -432,7 +432,7 @@ subroutine create_metadata_variable(this,vname,rc) integer :: rank,lb(1),ub(1) 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_AttributeGet(field,name="LONG_NAME",isPresent=is_present,_RC) @@ -448,7 +448,7 @@ subroutine create_metadata_variable(this,vname,rc) units = 'unknown' endif -! -- in future, replace keyword station by index_name_x as in trajectory sampler +! -- in future, replace keyword station by index_name_x as in trajectory sampler ! if (field_rank==2) then ! vdims = this%index_name_x ! else if (field_rank==3) then @@ -494,8 +494,8 @@ subroutine append_file(this,current_time,rc) real(kind=REAL32), pointer :: p_rt_3d(:,:),p_rt_2d(:) ! root LS real(kind=REAL32), pointer :: p_rt_3d_aux(:,:) real(kind=REAL32), allocatable :: p_new_lev(:,:,:) - real(kind=REAL32), allocatable :: p_dst_t(:,:) - + real(kind=REAL32), allocatable :: p_dst_t(:,:) + real(kind=REAL32), allocatable :: arr(:,:) character(len=ESMF_MAXSTR), allocatable :: fieldNameList(:) character(len=ESMF_MAXSTR) :: xname @@ -511,8 +511,8 @@ subroutine append_file(this,current_time,rc) type(ESMF_Field) :: field_ds_2d, field_ds_3d type(ESMF_Field) :: field_chunk_2d, field_chunk_3d - type(ESMF_Field) :: field_rt_2d, field_rt_3d - + type(ESMF_Field) :: field_rt_2d, field_rt_3d + integer :: sec integer, allocatable :: ix(:) ! counter for each obs(k)%nobs_epoch logical :: EX ! file @@ -522,9 +522,9 @@ subroutine append_file(this,current_time,rc) integer :: is, ie, ierr integer :: M, N, ip - + this%obs_written=this%obs_written+1 - + !__ 1. put_var: time variable ! ! @@ -539,14 +539,14 @@ subroutine append_file(this,current_time,rc) ! ungridded_dim from src to dst [regrid] ! lm = this%vdata%lm - field_ds_2d = ESMF_FieldCreate (this%LS_ds, name='field_2d_ds', typekind=ESMF_TYPEKIND_R4, _RC) + field_ds_2d = ESMF_FieldCreate (this%LS_ds, name='field_2d_ds', typekind=ESMF_TYPEKIND_R4, _RC) field_chunk_2d = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) dst_field = ESMF_FieldCreate (this%LS_ds, name='dst_field', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) - - + + ! caution about zero-sized array for MPI - ! redist + ! redist nx_sum = this%nstation call ESMF_VMGetCurrent(vm,_RC) call ESMF_VMGet(vm, mpiCommunicator=mpic, petCount=petCount, localPet=mypet, _RC) @@ -584,7 +584,7 @@ subroutine append_file(this,current_time,rc) allocate ( p_rt_3d_aux(nx_sum, lm) ) else allocate ( p_rt_3d(lm, 1) ) - allocate ( p_rt_3d_aux(1,lm) ) + allocate ( p_rt_3d_aux(1,lm) ) end if @@ -619,7 +619,7 @@ subroutine append_file(this,current_time,rc) else if (rank==3) then ! -- CS-> LS_ds; ds->chunk; gather - ! + ! call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_3d,_RC) call ESMF_FieldGet(dst_field,localDE=0,farrayptr=p_dst_3d,_RC) call this%regridder%regrid(p_src_3d,p_dst_3d,_RC) @@ -628,7 +628,7 @@ subroutine append_file(this,current_time,rc) gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) field_chunk_3d = ESMF_FieldCreate (this%LS_chunk, name='field_3d_chunk', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) - call ESMF_FieldGet(field_ds_3d,localDE=0,farrayPtr=p_ds_3d,_RC) + call ESMF_FieldGet(field_ds_3d,localDE=0,farrayPtr=p_ds_3d,_RC) call ESMF_FieldGet(field_chunk_3d,localDE=0,farrayPtr=p_chunk_3d,_RC) ! p_ds_3d(lm, nx) @@ -843,7 +843,7 @@ subroutine get_file_start_time(this,start_time,time_units,rc) _RETURN(_SUCCESS) end subroutine get_file_start_time - + ! TODO: delete and use system utilities when available Subroutine count_substring (str, t, ncount, rc) character (len=*), intent(in) :: str From 72218e6e8b38ec9d3f1643cf1bde0565306be5a8 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Tue, 21 May 2024 13:31:31 -0600 Subject: [PATCH 08/21] add deallocate(p_dst_t) --- gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index 4d5c08b78b6b..70520a46e4a6 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -387,7 +387,7 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() - print*, 'list item%xname', trim(item%xname) + !!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 @@ -643,6 +643,7 @@ subroutine append_file(this,current_time,rc) p_rt_3d_aux(1,k), recvcount, displs, MPI_REAL,& iroot, mpic, ierr ) end do + deallocate(p_dst_t) p_rt_3d = reshape(p_rt_3d_aux, shape(p_rt_3d), order=[2,1]) else call MPI_gatherv ( p_chunk_3d, nsend_v, MPI_REAL, & From e85a6448eba72e474b8844f1709710064585f5b7 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Tue, 21 May 2024 21:00:32 -0600 Subject: [PATCH 09/21] update --- .../Sampler/MAPL_StationSamplerMod.F90 | 44 +++++++++---------- 1 file changed, 21 insertions(+), 23 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index 70520a46e4a6..096d39d0cf9f 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -431,38 +431,33 @@ subroutine create_metadata_variable(this,vname,rc) character(len=ESMF_MAXSTR) :: var_name,long_name,units,vdims integer :: rank,lb(1),ub(1) integer :: k, ig - + integer, allocatable :: chunksizes(:) call ESMF_FieldBundleGet(this%bundle,vname,field=field,_RC) call ESMF_FieldGet(field,name=var_name,rank=field_rank,_RC) call ESMF_AttributeGet(field,name="LONG_NAME",isPresent=is_present,_RC) + long_name = var_name if ( is_present ) then call ESMF_AttributeGet (FIELD, NAME="LONG_NAME",VALUE=long_name, _RC) - else - long_name = var_name endif call ESMF_AttributeGet(field,name="UNITS",isPresent=is_present,_RC) + units = 'unknown' if ( is_present ) then call ESMF_AttributeGet (FIELD, NAME="UNITS",VALUE=units, _RC) - else - units = 'unknown' endif -! -- in future, replace keyword station by index_name_x as in trajectory sampler -! 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 - - if (field_rank==2) then - vdims = "station_index,time" - v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=[this%nstation,1]) - else if (field_rank==3) then - vdims = "lev,station_index,time" + vdims = "station_index,time" + select case (field_rank) + case(2) + chunksizes = [this%nstation,1] + case(3) + vdims = "lev,"//trim(vdims) call ESMF_FieldGet(field,ungriddedLBound=lb,ungriddedUBound=ub,_RC) - v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=[ub(1)-lb(1)+1,1,1]) - end if + chunksizes = [ub(1)-lb(1)+1,1,1] + case default + _FAIL('unsupported rank') + end select + v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=chunksizes) call v%add_attribute('units',trim(units)) call v%add_attribute('long_name',trim(long_name)) @@ -595,7 +590,8 @@ subroutine append_file(this,current_time,rc) call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) call ESMF_FieldGet(src_field,rank=rank,_RC) - if (rank==2) then + select case (rank) + case(2) call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_2d,_RC) call ESMF_FieldGet(field_ds_2d,localDE=0,farrayptr=p_ds_2d,_RC) call ESMF_FieldGet(field_chunk_2d,localDE=0,farrayPtr=p_chunk_2d,_RC) @@ -617,7 +613,7 @@ subroutine append_file(this,current_time,rc) start=[1,this%obs_written],count=[this%nstation,1],_RC) end if - else if (rank==3) then + case(3) ! -- CS-> LS_ds; ds->chunk; gather ! call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_3d,_RC) @@ -666,9 +662,11 @@ subroutine append_file(this,current_time,rc) start=[1,1,this%obs_written],count=[nz,nx,1],_RC) !note: lev,station,time end if - else + case default _FAIL('grid2LS regridder: rank > 3 not implemented') - end if + end select + else + _FAIL ('ItemType vector not supported') endif call iter%next() From 8a23b743833c2a7884b5ad97b0521bed26503c0b Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Wed, 22 May 2024 13:02:45 -0600 Subject: [PATCH 10/21] test remove mpi barrier --- gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 b/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 index 6201f50e2754..6e4dbb76ae2c 100644 --- a/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 @@ -354,7 +354,6 @@ call ESMF_FieldHaloStore (fieldI4, routehandle=RH_halo, _RC) call ESMF_FieldHalo (fieldI4, routehandle=RH_halo, _RC) - call ESMF_VMBarrier (vm, _RC) k=0 do i=eLB(1), eUB(1) @@ -411,7 +410,7 @@ lons(i) = lons_ptr (ix, jx) lats(i) = lats_ptr (ix, jx) end do - call ESMF_VMBarrier (vm, _RC) + iroot=0 if (mapl_am_i_root()) then @@ -626,7 +625,7 @@ iy = this%index_mask(2,j) p_dst_2d(j) = p_src_2d(ix, iy) end do - call MPI_Barrier(mpic, status) +!! call MPI_Barrier(mpic, status) nsend = nx call MPI_gatherv ( p_dst_2d, nsend, MPI_REAL, & p_dst_2d_full, this%recvcounts, this%displs, MPI_REAL,& @@ -648,7 +647,7 @@ p_dst_3d(m) = p_src_3d(ix, iy, k) end do end do - call MPI_Barrier(mpic, status) +!! call MPI_Barrier(mpic, status) !! write(6,'(2x,a,2x,i5,3x,10f8.1)') 'pet, p_dst_3d(j)', mypet, p_dst_3d(::10) nsend = nx * nz call MPI_gatherv ( p_dst_3d, nsend, MPI_REAL, & From 0c547048109f178003503056eb18a45536f9db8b Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Thu, 23 May 2024 09:16:35 -0600 Subject: [PATCH 11/21] - add timer to put_var for 2D and 3D var in station and mask sampler - removed chunksizes in station sampler --- gridcomps/History/MAPL_HistoryGridComp.F90 | 4 ++-- gridcomps/History/Sampler/MAPL_GeosatMaskMod.F90 | 5 ++++- .../History/Sampler/MAPL_GeosatMaskMod_smod.F90 | 6 ++++++ .../History/Sampler/MAPL_StationSamplerMod.F90 | 13 +++++++++++-- 4 files changed, 23 insertions(+), 5 deletions(-) diff --git a/gridcomps/History/MAPL_HistoryGridComp.F90 b/gridcomps/History/MAPL_HistoryGridComp.F90 index d8a2d6321449..c287dae7b07d 100644 --- a/gridcomps/History/MAPL_HistoryGridComp.F90 +++ b/gridcomps/History/MAPL_HistoryGridComp.F90 @@ -2426,10 +2426,10 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) call list(n)%trajectory%initialize(items=list(n)%items,bundle=list(n)%bundle,timeinfo=list(n)%timeInfo,vdata=list(n)%vdata,_RC) IntState%stampoffset(n) = list(n)%trajectory%epoch_frequency elseif (list(n)%sampler_spec == 'mask') then - list(n)%mask_sampler = MaskSamplerGeosat(cfg,string,clock,_RC) + list(n)%mask_sampler = MaskSamplerGeosat(cfg,string,clock,genstate=GENSTATE,_RC) call list(n)%mask_sampler%initialize(items=list(n)%items,bundle=list(n)%bundle,timeinfo=list(n)%timeInfo,vdata=list(n)%vdata,_RC) elseif (list(n)%sampler_spec == 'station') then - list(n)%station_sampler = StationSampler (list(n)%bundle, trim(list(n)%stationIdFile), nskip_line=list(n)%stationSkipLine, _RC) + list(n)%station_sampler = StationSampler (list(n)%bundle, trim(list(n)%stationIdFile), nskip_line=list(n)%stationSkipLine, genstate=GENSTATE, _RC) call list(n)%station_sampler%add_metadata_route_handle(items=list(n)%items,bundle=list(n)%bundle,timeinfo=list(n)%timeInfo,vdata=list(n)%vdata,_RC) else global_attributes = list(n)%global_atts%define_collection_attributes(_RC) diff --git a/gridcomps/History/Sampler/MAPL_GeosatMaskMod.F90 b/gridcomps/History/Sampler/MAPL_GeosatMaskMod.F90 index 5674a1b2f1ca..c1dc6e1bbd96 100644 --- a/gridcomps/History/Sampler/MAPL_GeosatMaskMod.F90 +++ b/gridcomps/History/Sampler/MAPL_GeosatMaskMod.F90 @@ -20,6 +20,7 @@ module MaskSamplerGeosatMod use MPI use pFIO_FileMetadataMod, only : FileMetadata use pFIO_NetCDF4_FileFormatterMod, only : NetCDF4_FileFormatter + use MAPL_GenericMod, only : MAPL_MetaComp, MAPL_TimerOn, MAPL_TimerOff use, intrinsic :: iso_fortran_env, only: REAL32 use, intrinsic :: iso_fortran_env, only: REAL64 use pflogger, only: Logger, logging @@ -76,6 +77,7 @@ module MaskSamplerGeosatMod real(kind=REAL64), allocatable :: lats(:) integer, allocatable :: recvcounts(:) integer, allocatable :: displs(:) + type(MAPL_MetaComp), pointer :: GENSTATE real(kind=ESMF_KIND_R8), pointer:: obsTime(:) real(kind=ESMF_KIND_R8), allocatable:: t_alongtrack(:) @@ -104,11 +106,12 @@ module MaskSamplerGeosatMod interface - module function MaskSamplerGeosat_from_config(config,string,clock,rc) result(mask) + module function MaskSamplerGeosat_from_config(config,string,clock,GENSTATE,rc) result(mask) type(MaskSamplerGeosat) :: mask type(ESMF_Config), intent(inout) :: config character(len=*), intent(in) :: string type(ESMF_Clock), intent(in) :: clock + type(MAPL_MetaComp), pointer, intent(in), optional :: GENSTATE integer, optional, intent(out) :: rc end function MaskSamplerGeosat_from_config diff --git a/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 b/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 index 6201f50e2754..e72b3783000b 100644 --- a/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 @@ -27,6 +27,8 @@ mask%clock=clock mask%grid_file_name='' + if (present(GENSTATE)) mask%GENSTATE => GENSTATE + call ESMF_ClockGet ( clock, CurrTime=currTime, _RC ) if (mapl_am_I_root()) write(6,*) 'string', string @@ -631,10 +633,12 @@ call MPI_gatherv ( p_dst_2d, nsend, MPI_REAL, & p_dst_2d_full, this%recvcounts, this%displs, MPI_REAL,& iroot, mpic, ierr ) + call MAPL_TimerOn(this%GENSTATE,"put2D") if (mapl_am_i_root()) then call this%formatter%put_var(item%xname,p_dst_2d_full,& start=[1,this%obs_written],count=[this%npt_mask_tot,1],_RC) end if + call MAPL_TimerOff(this%GENSTATE,"put2D") else if (rank==3) then call ESMF_FieldGet(src_field,farrayptr=p_src_3d,_RC) call ESMF_FieldGet(src_field,ungriddedLBound=lb,ungriddedUBound=ub,_RC) @@ -654,6 +658,7 @@ call MPI_gatherv ( p_dst_3d, nsend, MPI_REAL, & p_dst_3d_full, recvcounts_3d, displs_3d, MPI_REAL,& iroot, mpic, ierr ) + call MAPL_TimerOn(this%GENSTATE,"put3D") if (mapl_am_i_root()) then allocate(arr(nz, this%npt_mask_tot), _STAT) arr=reshape(p_dst_3d_full,[nz,this%npt_mask_tot],order=[1,2]) @@ -662,6 +667,7 @@ !note: lev,station,time deallocate(arr, _STAT) end if + call MAPL_TimerOff(this%GENSTATE,"put3D") else _FAIL('grid2LS regridder: rank > 3 not implemented') end if diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index 096d39d0cf9f..d92d54c64851 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -11,6 +11,7 @@ module StationSamplerMod use MAPL_BaseMod use MAPL_CommsMod use MAPL_LocstreamRegridderMod + use MAPL_GenericMod, only : MAPL_MetaComp, MAPL_TimerOn, MAPL_TimerOff use MPI, only : MPI_INTEGER, MPI_REAL, MPI_REAL8 use, intrinsic :: iso_fortran_env, only: INT64 use, intrinsic :: iso_fortran_env, only: REAL32 @@ -31,6 +32,7 @@ module StationSamplerMod type(GriddedIOitemVector) :: items logical :: do_vertical_regrid logical :: level_by_level + type(MAPL_MetaComp), pointer :: GENSTATE integer :: nstation integer, allocatable :: station_id(:) @@ -64,13 +66,14 @@ module StationSamplerMod contains - function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(sampler) + function new_StationSampler_readfile (bundle, filename, nskip_line, GENSTATE, rc) result(sampler) use pflogger, only : Logger, logging implicit none type(StationSampler) :: sampler type(ESMF_FieldBundle), intent(in) :: bundle character(len=*), intent(in) :: filename integer, optional, intent(in) :: nskip_line + type(MAPL_MetaComp), pointer, intent(in), optional :: GENSTATE integer, optional, intent(out) :: rc type(ESMF_VM) :: vm @@ -241,6 +244,7 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, rc) result(s allocate(sampler%elevs(nstation), _STAT) end if sampler%index_name_x = 'station_index' + if (present(GENSTATE)) sampler%GENSTATE => GENSTATE !__ 2. create LocStreamFactory, then LS_rt including route_handle @@ -457,7 +461,8 @@ subroutine create_metadata_variable(this,vname,rc) case default _FAIL('unsupported rank') end select - v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=chunksizes) + ! v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=chunksizes) + v = variable(type=PFIO_REAL32,dimensions=trim(vdims)) call v%add_attribute('units',trim(units)) call v%add_attribute('long_name',trim(long_name)) @@ -608,10 +613,12 @@ subroutine append_file(this,current_time,rc) ! end do ! end if + call MAPL_TimerOn(this%GENSTATE,"put2D") if (mapl_am_i_root()) then call this%formatter%put_var(trim(item%xname),p_rt_2d,& start=[1,this%obs_written],count=[this%nstation,1],_RC) end if + call MAPL_TimerOff(this%GENSTATE,"put2D") case(3) ! -- CS-> LS_ds; ds->chunk; gather @@ -656,12 +663,14 @@ subroutine append_file(this,current_time,rc) ! end if ! if (mapl_am_i_root()) write(6,*) 'regrid + gatherV in 3D' + call MAPL_TimerOn(this%GENSTATE,"put3D") if (mapl_am_i_root()) then nz=size(p_rt_3d,1); nx=size(p_rt_3d,2) call this%formatter%put_var(trim(item%xname),p_rt_3d,& start=[1,1,this%obs_written],count=[nz,nx,1],_RC) !note: lev,station,time end if + call MAPL_TimerOff(this%GENSTATE,"put3D") case default _FAIL('grid2LS regridder: rank > 3 not implemented') end select From d8896c1b60de6e2fce21a41c724563f1cafd82c2 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Thu, 23 May 2024 13:10:09 -0400 Subject: [PATCH 12/21] - remove chunksizes in mask sampler --- gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 b/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 index e72b3783000b..7faee14cd798 100644 --- a/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 @@ -529,11 +529,11 @@ endif if (field_rank==2) then vdims = "mask_index,time" - v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=[this%npt_mask_tot,1]) + v = variable(type=PFIO_REAL32,dimensions=trim(vdims)) else if (field_rank==3) then vdims = "lev,mask_index,time" call ESMF_FieldGet(field,ungriddedLBound=lb,ungriddedUBound=ub,_RC) - v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=[ub(1)-lb(1)+1,1,1]) + v = variable(type=PFIO_REAL32,dimensions=trim(vdims)) end if call v%add_attribute('units', trim(units)) call v%add_attribute('long_name', trim(long_name)) From 38a2cb1eb35bee7d073f6174d5e0a7b3028369d1 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Thu, 23 May 2024 11:35:33 -0600 Subject: [PATCH 13/21] - delete keyword accumulate in mask sampler --- gridcomps/History/Sampler/MAPL_GeosatMaskMod.F90 | 6 +++--- gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 | 7 ++----- 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_GeosatMaskMod.F90 b/gridcomps/History/Sampler/MAPL_GeosatMaskMod.F90 index c1dc6e1bbd96..8fc982729673 100644 --- a/gridcomps/History/Sampler/MAPL_GeosatMaskMod.F90 +++ b/gridcomps/History/Sampler/MAPL_GeosatMaskMod.F90 @@ -94,7 +94,7 @@ module MaskSamplerGeosatMod procedure :: add_metadata procedure :: create_file_handle procedure :: close_file_handle - procedure :: append_file => regrid_accumulate_append_file + procedure :: append_file => regrid_append_file ! procedure :: create_new_bundle procedure :: create_grid => create_Geosat_grid_find_mask procedure :: compute_time_for_current @@ -153,11 +153,11 @@ module subroutine close_file_handle(this,rc) integer, optional, intent(out) :: rc end subroutine close_file_handle - module subroutine regrid_accumulate_append_file(this,current_time,rc) + module subroutine regrid_append_file(this,current_time,rc) class(MaskSamplerGeosat), intent(inout) :: this type(ESMF_Time), intent(inout) :: current_time integer, optional, intent(out) :: rc - end subroutine regrid_accumulate_append_file + end subroutine regrid_append_file module function compute_time_for_current(this,current_time,rc) result(rtime) class(MaskSamplerGeosat), intent(inout) :: this diff --git a/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 b/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 index f4156fc9f0b2..df7917cff814 100644 --- a/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 @@ -326,7 +326,6 @@ obs_lats = lats_ds * MAPL_DEGREES_TO_RADIANS_R8 nx = size ( lons_ds ) allocate ( II(nx), JJ(nx), _STAT ) - call MPI_Barrier(mpic, status) call MAPL_GetHorzIJIndex(nx,II,JJ,lonR8=obs_lons,latR8=obs_lats,grid=grid,_RC) call ESMF_VMBarrier (vm, _RC) @@ -547,7 +546,7 @@ end procedure add_metadata - module procedure regrid_accumulate_append_file + module procedure regrid_append_file ! implicit none integer :: status @@ -627,7 +626,6 @@ iy = this%index_mask(2,j) p_dst_2d(j) = p_src_2d(ix, iy) end do -!! call MPI_Barrier(mpic, status) nsend = nx call MPI_gatherv ( p_dst_2d, nsend, MPI_REAL, & p_dst_2d_full, this%recvcounts, this%displs, MPI_REAL,& @@ -651,7 +649,6 @@ p_dst_3d(m) = p_src_3d(ix, iy, k) end do end do -!! call MPI_Barrier(mpic, status) !! write(6,'(2x,a,2x,i5,3x,10f8.1)') 'pet, p_dst_3d(j)', mypet, p_dst_3d(::10) nsend = nx * nz call MPI_gatherv ( p_dst_3d, nsend, MPI_REAL, & @@ -676,7 +673,7 @@ end do _RETURN(_SUCCESS) - end procedure regrid_accumulate_append_file + end procedure regrid_append_file From 953ff8e7d4761bd339c8c136ac13f9b9b4bd1a7d Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Thu, 23 May 2024 22:36:01 -0600 Subject: [PATCH 14/21] - refine regrid in station sampler --- .../Sampler/MAPL_StationSamplerMod.F90 | 51 ++++++++++++------- 1 file changed, 32 insertions(+), 19 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index d92d54c64851..73d24a9efd29 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -335,7 +335,7 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) integer :: lb(ESMF_MAXDIM) logical :: do_vertical_regrid integer :: status - integer :: i + integer :: i, lm character(len=ESMF_MAXSTR), allocatable :: fieldNameList(:) character(len=ESMF_MAXSTR) :: var_name, long_name, units, vdims @@ -509,9 +509,12 @@ subroutine append_file(this,current_time,rc) integer :: n0, nx, nx2 integer :: na, nb, nx_sum, nsend, nsend_v - type(ESMF_Field) :: field_ds_2d, field_ds_3d - type(ESMF_Field) :: field_chunk_2d, field_chunk_3d - type(ESMF_Field) :: field_rt_2d, field_rt_3d + ! intermediate fields as placeholder + type(ESMF_Field) :: field_ds_2d + type(ESMF_Field) :: field_ds_3d + type(ESMF_Field) :: field_chunk_2d + type(ESMF_Field) :: field_chunk_3d + integer :: sec integer, allocatable :: ix(:) ! counter for each obs(k)%nobs_epoch @@ -538,16 +541,10 @@ subroutine append_file(this,current_time,rc) !__ 2. regrid + put_var: ! ungridded_dim from src to dst [regrid] ! - lm = this%vdata%lm - field_ds_2d = ESMF_FieldCreate (this%LS_ds, name='field_2d_ds', typekind=ESMF_TYPEKIND_R4, _RC) - field_chunk_2d = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) - dst_field = ESMF_FieldCreate (this%LS_ds, name='dst_field', typekind=ESMF_TYPEKIND_R4, & - gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) - - ! caution about zero-sized array for MPI ! redist nx_sum = this%nstation + lm = this%vdata%lm call ESMF_VMGetCurrent(vm,_RC) call ESMF_VMGet(vm, mpiCommunicator=mpic, petCount=petCount, localPet=mypet, _RC) @@ -568,6 +565,7 @@ subroutine append_file(this,current_time,rc) 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 ) @@ -588,11 +586,24 @@ subroutine append_file(this,current_time,rc) end if + !__ Aux. field + ! + field_ds_2d = ESMF_FieldCreate (this%LS_ds, name='field_2d_ds', typekind=ESMF_TYPEKIND_R4, _RC) + field_chunk_2d = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) + dst_field = ESMF_FieldCreate (this%LS_ds, name='dst_field', typekind=ESMF_TYPEKIND_R4, & + gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) + + field_ds_3d = ESMF_FieldCreate (this%LS_ds, name='field_3d_ds', typekind=ESMF_TYPEKIND_R4, & + gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) + field_chunk_3d = ESMF_FieldCreate (this%LS_chunk, name='field_3d_chunk', typekind=ESMF_TYPEKIND_R4, & + gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) + + iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() if (item%itemType == ItemTypeScalar) then - +!! if (mapl_am_i_root()) write(6,*) 'item%xname=', trim(item%xname) call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) call ESMF_FieldGet(src_field,rank=rank,_RC) select case (rank) @@ -600,6 +611,7 @@ subroutine append_file(this,current_time,rc) call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_2d,_RC) call ESMF_FieldGet(field_ds_2d,localDE=0,farrayptr=p_ds_2d,_RC) call ESMF_FieldGet(field_chunk_2d,localDE=0,farrayPtr=p_chunk_2d,_RC) + p_ds_2d = 0.0 call this%regridder%regrid(p_src_2d,p_ds_2d,_RC) call ESMF_FieldRedist(field_ds_2d, field_chunk_2d, this%RH, _RC ) @@ -625,18 +637,15 @@ subroutine append_file(this,current_time,rc) ! call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_3d,_RC) call ESMF_FieldGet(dst_field,localDE=0,farrayptr=p_dst_3d,_RC) + p_dst_3d = 0.0 call this%regridder%regrid(p_src_3d,p_dst_3d,_RC) - field_ds_3d = ESMF_FieldCreate (this%LS_ds, name='field_3d_ds', typekind=ESMF_TYPEKIND_R4, & - gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) - field_chunk_3d = ESMF_FieldCreate (this%LS_chunk, name='field_3d_chunk', typekind=ESMF_TYPEKIND_R4, & - gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) call ESMF_FieldGet(field_ds_3d,localDE=0,farrayPtr=p_ds_3d,_RC) call ESMF_FieldGet(field_chunk_3d,localDE=0,farrayPtr=p_chunk_3d,_RC) ! p_ds_3d(lm, nx) p_ds_3d = reshape(p_dst_3d, shape(p_ds_3d), order=[2,1]) - call ESMF_FieldRedist(field_ds_3d, field_chunk_3d, this%RH, _RC) + call ESMF_FieldRegrid(field_ds_3d, field_chunk_3d, this%RH, _RC) if (this%level_by_level) then ! p_chunk_3d (lm, nx) @@ -653,8 +662,6 @@ subroutine append_file(this,current_time,rc) p_rt_3d, recvcount_v, displs_v, MPI_REAL,& iroot, mpic, ierr ) end if - call ESMF_FieldDestroy(field_ds_3d,noGarbage=.true.,_RC) - call ESMF_FieldDestroy(field_chunk_3d,noGarbage=.true.,_RC) ! if (mapl_am_i_root()) then ! do j=1, nx_sum, 500 @@ -681,6 +688,12 @@ subroutine append_file(this,current_time,rc) call iter%next() end do + call ESMF_FieldDestroy(field_ds_2d, noGarbage=.true., _RC) + call ESMF_FieldDestroy(field_chunk_2d, noGarbage=.true., _RC) + call ESMF_FieldDestroy(field_ds_3d, noGarbage=.true., _RC) + call ESMF_FieldDestroy(field_chunk_3d, noGarbage=.true., _RC) + call ESMF_FieldDestroy(dst_field, noGarbage=.true., _RC) + _RETURN(_SUCCESS) end subroutine append_file From ce6dbd06b87d42fb571abb92363d7189fca9d4db Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Thu, 13 Jun 2024 07:17:14 -0600 Subject: [PATCH 15/21] add notes on redistribute LS_ds to LS_chunk --- gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index 73d24a9efd29..8dc0f98eaf46 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -645,8 +645,10 @@ subroutine append_file(this,current_time,rc) ! p_ds_3d(lm, nx) p_ds_3d = reshape(p_dst_3d, shape(p_ds_3d), order=[2,1]) + ! ... call ESMF_FieldRegrid(field_ds_3d, field_chunk_3d, this%RH, _RC) - + ! redistributed: slower check. + if (this%level_by_level) then ! p_chunk_3d (lm, nx) allocate (p_dst_t, source = reshape(p_chunk_3d, [size(p_chunk_3d,2),size(p_chunk_3d,1)], order=[2,1])) From 1ec7efa276b2bcb561888b0e56d755d0beb069c2 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Thu, 13 Jun 2024 10:22:27 -0600 Subject: [PATCH 16/21] Add a few timers --- gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index 8dc0f98eaf46..b6165e9016fa 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -637,16 +637,24 @@ subroutine append_file(this,current_time,rc) ! call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_3d,_RC) call ESMF_FieldGet(dst_field,localDE=0,farrayptr=p_dst_3d,_RC) + call MAPL_TimerOn(this%GENSTATE,"assign p_dst_3d=0") p_dst_3d = 0.0 - call this%regridder%regrid(p_src_3d,p_dst_3d,_RC) + print*, 'shape(p_dst_3d)', shape(p_dst_3d) + call MAPL_TimerOff(this%GENSTATE,"assign p_dst_3d=0") + call MAPL_TimerOn(this%GENSTATE,"3d_regrid") + call this%regridder%regrid(p_src_3d,p_dst_3d,_RC) + call MAPL_TimerOff(this%GENSTATE,"3d_regrid") + call ESMF_FieldGet(field_ds_3d,localDE=0,farrayPtr=p_ds_3d,_RC) call ESMF_FieldGet(field_chunk_3d,localDE=0,farrayPtr=p_chunk_3d,_RC) + print*, 'shape(p_ds_3d)', shape(p_ds_3d) ! p_ds_3d(lm, nx) p_ds_3d = reshape(p_dst_3d, shape(p_ds_3d), order=[2,1]) ! ... - call ESMF_FieldRegrid(field_ds_3d, field_chunk_3d, this%RH, _RC) + !! call ESMF_FieldRegrid(field_ds_3d, field_chunk_3d, this%RH, _RC) + call ESMF_FieldRedist(field_ds_3d, field_chunk_3d, this%RH, _RC) ! redistributed: slower check. if (this%level_by_level) then From c74f83ee4d155f829705b2469f782f6f24d53f69 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Thu, 13 Jun 2024 23:59:14 -0600 Subject: [PATCH 17/21] use more efficient FA + field for consecutive regridding CS -- LS_ds-- LS_chunk -- LS_root --- .../Sampler/MAPL_StationSamplerMod.F90 | 67 ++++++++++--------- 1 file changed, 36 insertions(+), 31 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index b6165e9016fa..a4702b0a39da 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -415,6 +415,7 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) call ESMF_FieldGet( chunk_field, localDE=0, farrayPtr=pt2, _RC ) pt1=0.0 pt2=0.0 +!! print*, 'shape(pt1 LS_ds)', shape(pt1) call ESMF_FieldRedistStore(src_field,chunk_field,this%RH,_RC) call ESMF_FieldDestroy(src_field,noGarbage=.true.,_RC) call ESMF_FieldDestroy(chunk_field,noGarbage=.true.,_RC) @@ -486,8 +487,10 @@ subroutine append_file(this,current_time,rc) integer :: ub(1), lb(1) type(GriddedIOitemVectorIterator) :: iter type(GriddedIOitem), pointer :: item - type(ESMF_Field) :: src_field,dst_field - real(kind=REAL32), pointer :: p_src_3d(:,:,:),p_src_2d(:,:) ! source + type(ESMF_Grid) :: grid + type(ESMF_Field) :: src_field ! ,dst_field + type(ESMF_Field) :: new_src_field,new_dst_field + real(kind=REAL32), pointer :: p_src_3d(:,:,:),p_src_2d(:,:), qin_3d(:,:,:) ! source real(kind=REAL32), pointer :: p_dst_3d(:,:) ! destination real(kind=REAL32), pointer :: p_ds_3d(:,:),p_ds_2d(:) ! distributed LS real(kind=REAL32), pointer :: p_chunk_3d(:,:),p_chunk_2d(:) ! chunk LS @@ -497,7 +500,7 @@ subroutine append_file(this,current_time,rc) real(kind=REAL32), allocatable :: p_dst_t(:,:) real(kind=REAL32), allocatable :: arr(:,:) - character(len=ESMF_MAXSTR), allocatable :: fieldNameList(:) + character(len=ESMF_MAXSTR), allocatable :: fieldNameList(:) character(len=ESMF_MAXSTR) :: xname real(kind=ESMF_KIND_R8), allocatable :: rtimes(:) @@ -511,11 +514,9 @@ subroutine append_file(this,current_time,rc) ! intermediate fields as placeholder type(ESMF_Field) :: field_ds_2d - type(ESMF_Field) :: field_ds_3d type(ESMF_Field) :: field_chunk_2d type(ESMF_Field) :: field_chunk_3d - integer :: sec integer, allocatable :: ix(:) ! counter for each obs(k)%nobs_epoch logical :: EX ! file @@ -588,16 +589,21 @@ subroutine append_file(this,current_time,rc) !__ Aux. field ! + field_ds_2d = ESMF_FieldCreate (this%LS_ds, name='field_2d_ds', typekind=ESMF_TYPEKIND_R4, _RC) field_chunk_2d = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) - dst_field = ESMF_FieldCreate (this%LS_ds, name='dst_field', typekind=ESMF_TYPEKIND_R4, & - gridToFieldMap=[1],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) + call ESMF_FieldBundleGet(this%bundle,grid=grid,_RC) - field_ds_3d = ESMF_FieldCreate (this%LS_ds, name='field_3d_ds', typekind=ESMF_TYPEKIND_R4, & + new_src_field = ESMF_FieldCreate (grid, name='new_src_field', typekind=ESMF_TYPEKIND_R4, & + gridToFieldMap=[2,3],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) + new_dst_field = ESMF_FieldCreate (this%LS_ds, name='new_dst_field', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) field_chunk_3d = ESMF_FieldCreate (this%LS_chunk, name='field_3d_chunk', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) + call ESMF_FieldGet(new_src_field,localDE=0,farrayPtr=p_src_3d,_RC) + call ESMF_FieldGet(new_dst_field,localDE=0,farrayPtr=p_dst_3d,_RC) + call ESMF_FieldGet(field_chunk_3d,localDE=0,farrayPtr=p_chunk_3d,_RC) iter = this%items%begin() do while (iter /= this%items%end()) @@ -611,7 +617,6 @@ subroutine append_file(this,current_time,rc) call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_2d,_RC) call ESMF_FieldGet(field_ds_2d,localDE=0,farrayptr=p_ds_2d,_RC) call ESMF_FieldGet(field_chunk_2d,localDE=0,farrayPtr=p_chunk_2d,_RC) - p_ds_2d = 0.0 call this%regridder%regrid(p_src_2d,p_ds_2d,_RC) call ESMF_FieldRedist(field_ds_2d, field_chunk_2d, this%RH, _RC ) @@ -635,28 +640,27 @@ subroutine append_file(this,current_time,rc) case(3) ! -- CS-> LS_ds; ds->chunk; gather ! - call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_3d,_RC) - call ESMF_FieldGet(dst_field,localDE=0,farrayptr=p_dst_3d,_RC) - call MAPL_TimerOn(this%GENSTATE,"assign p_dst_3d=0") - p_dst_3d = 0.0 - print*, 'shape(p_dst_3d)', shape(p_dst_3d) - call MAPL_TimerOff(this%GENSTATE,"assign p_dst_3d=0") - - call MAPL_TimerOn(this%GENSTATE,"3d_regrid") - call this%regridder%regrid(p_src_3d,p_dst_3d,_RC) + call ESMF_FieldGet(src_field,localDE=0,farrayptr=qin_3d,_RC) + p_src_3d = reshape(qin_3d,shape(p_src_3d), order = [2,3,1]) + +! print*, 'shape(p_src_3d)', shape(p_src_3d) +! print*, 'shape(p_dst_3d)', shape(p_dst_3d) + + call MAPL_TimerOn(this%GENSTATE,"3d_regrid") +!! call this%regridder%regrid(p_src_3d,p_dst_3d,_RC) + call ESMF_FieldRegrid (new_src_field, new_dst_field, this%regridder%route_handle, _RC) call MAPL_TimerOff(this%GENSTATE,"3d_regrid") - - call ESMF_FieldGet(field_ds_3d,localDE=0,farrayPtr=p_ds_3d,_RC) - call ESMF_FieldGet(field_chunk_3d,localDE=0,farrayPtr=p_chunk_3d,_RC) - print*, 'shape(p_ds_3d)', shape(p_ds_3d) - - ! p_ds_3d(lm, nx) - p_ds_3d = reshape(p_dst_3d, shape(p_ds_3d), order=[2,1]) - ! ... + + + call MAPL_TimerOn(this%GENSTATE,"ESMF_FieldRedist") + ! ... !! call ESMF_FieldRegrid(field_ds_3d, field_chunk_3d, this%RH, _RC) - call ESMF_FieldRedist(field_ds_3d, field_chunk_3d, this%RH, _RC) - ! redistributed: slower check. - + call ESMF_FieldRedist(new_dst_field, field_chunk_3d, this%RH, _RC) + + call MAPL_TimerOff(this%GENSTATE,"ESMF_FieldRedist") + + + call MAPL_TimerOn(this%GENSTATE,"gahterv") if (this%level_by_level) then ! p_chunk_3d (lm, nx) allocate (p_dst_t, source = reshape(p_chunk_3d, [size(p_chunk_3d,2),size(p_chunk_3d,1)], order=[2,1])) @@ -672,6 +676,7 @@ subroutine append_file(this,current_time,rc) p_rt_3d, recvcount_v, displs_v, MPI_REAL,& iroot, mpic, ierr ) end if + call MAPL_TimerOff(this%GENSTATE,"gahterv") ! if (mapl_am_i_root()) then ! do j=1, nx_sum, 500 @@ -700,9 +705,9 @@ subroutine append_file(this,current_time,rc) call ESMF_FieldDestroy(field_ds_2d, noGarbage=.true., _RC) call ESMF_FieldDestroy(field_chunk_2d, noGarbage=.true., _RC) - call ESMF_FieldDestroy(field_ds_3d, noGarbage=.true., _RC) call ESMF_FieldDestroy(field_chunk_3d, noGarbage=.true., _RC) - call ESMF_FieldDestroy(dst_field, noGarbage=.true., _RC) + call ESMF_FieldDestroy(new_dst_field, noGarbage=.true., _RC) + call ESMF_FieldDestroy(new_src_field, noGarbage=.true., _RC) _RETURN(_SUCCESS) end subroutine append_file From 66f3b5834312c7e8d73ca7eca3bff45a983690c3 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Fri, 14 Jun 2024 10:37:55 -0600 Subject: [PATCH 18/21] Add more timing --- .../Sampler/MAPL_StationSamplerMod.F90 | 76 +++++++++---------- 1 file changed, 36 insertions(+), 40 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index a4702b0a39da..a9e0dbee574d 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -566,7 +566,6 @@ subroutine append_file(this,current_time,rc) 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 ) @@ -589,47 +588,45 @@ subroutine append_file(this,current_time,rc) !__ Aux. field ! + call MAPL_TimerOn(this%GENSTATE,"FieldCreate") - field_ds_2d = ESMF_FieldCreate (this%LS_ds, name='field_2d_ds', typekind=ESMF_TYPEKIND_R4, _RC) - field_chunk_2d = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) call ESMF_FieldBundleGet(this%bundle,grid=grid,_RC) - - new_src_field = ESMF_FieldCreate (grid, name='new_src_field', typekind=ESMF_TYPEKIND_R4, & + field_ds_2d = ESMF_FieldCreate (this%LS_ds, name='field_2d_ds', typekind=ESMF_TYPEKIND_R4, _RC) + field_chunk_2d = ESMF_FieldCreate (this%LS_chunk, name='field_2d_chunk', typekind=ESMF_TYPEKIND_R4, _RC) + new_src_field = ESMF_FieldCreate (grid, name='new_src_field', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[2,3],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) - new_dst_field = ESMF_FieldCreate (this%LS_ds, name='new_dst_field', typekind=ESMF_TYPEKIND_R4, & + new_dst_field = ESMF_FieldCreate (this%LS_ds, name='new_dst_field', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) field_chunk_3d = ESMF_FieldCreate (this%LS_chunk, name='field_3d_chunk', typekind=ESMF_TYPEKIND_R4, & gridToFieldMap=[2],ungriddedLBound=[1],ungriddedUBound=[lm],_RC) - call ESMF_FieldGet(new_src_field,localDE=0,farrayPtr=p_src_3d,_RC) - call ESMF_FieldGet(new_dst_field,localDE=0,farrayPtr=p_dst_3d,_RC) - call ESMF_FieldGet(field_chunk_3d,localDE=0,farrayPtr=p_chunk_3d,_RC) + call ESMF_FieldGet(field_ds_2d, localDE=0, farrayptr=p_ds_2d, _RC) + call ESMF_FieldGet(field_chunk_2d,localDE=0, farrayPtr=p_chunk_2d, _RC) + call ESMF_FieldGet(new_src_field, localDE=0, farrayPtr=p_src_3d, _RC) + call ESMF_FieldGet(new_dst_field, localDE=0, farrayPtr=p_dst_3d, _RC) + call ESMF_FieldGet(field_chunk_3d,localDE=0, farrayPtr=p_chunk_3d, _RC) + + call MAPL_TimerOff(this%GENSTATE,"FieldCreate") + + + call MAPL_TimerOn(this%GENSTATE,"Full_regrid") iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() if (item%itemType == ItemTypeScalar) then -!! if (mapl_am_i_root()) write(6,*) 'item%xname=', trim(item%xname) + !! if (mapl_am_i_root()) write(6,*) 'item%xname=', trim(item%xname) call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) call ESMF_FieldGet(src_field,rank=rank,_RC) select case (rank) case(2) call ESMF_FieldGet(src_field,localDE=0,farrayptr=p_src_2d,_RC) - call ESMF_FieldGet(field_ds_2d,localDE=0,farrayptr=p_ds_2d,_RC) - call ESMF_FieldGet(field_chunk_2d,localDE=0,farrayPtr=p_chunk_2d,_RC) - - call this%regridder%regrid(p_src_2d,p_ds_2d,_RC) - call ESMF_FieldRedist(field_ds_2d, field_chunk_2d, this%RH, _RC ) + call ESMF_FieldRegrid (src_field, field_ds_2d, this%regridder%route_handle, _RC) + call ESMF_FieldRedist (field_ds_2d, field_chunk_2d, this%RH, _RC ) call MPI_gatherv ( p_chunk_2d, nsend, MPI_REAL, & p_rt_2d, recvcount, displs, MPI_REAL,& iroot, mpic, ierr ) -! if (mapl_am_i_root()) then -! do j=1, nx_sum, 500000 -! write(6,*) 'p_rt_2d', p_rt_2d(j) -! end do -! end if - call MAPL_TimerOn(this%GENSTATE,"put2D") if (mapl_am_i_root()) then call this%formatter%put_var(trim(item%xname),p_rt_2d,& @@ -641,24 +638,25 @@ subroutine append_file(this,current_time,rc) ! -- CS-> LS_ds; ds->chunk; gather ! call ESMF_FieldGet(src_field,localDE=0,farrayptr=qin_3d,_RC) - p_src_3d = reshape(qin_3d,shape(p_src_3d), order = [2,3,1]) -! print*, 'shape(p_src_3d)', shape(p_src_3d) -! print*, 'shape(p_dst_3d)', shape(p_dst_3d) + !! -- did not improve performance + !!call MAPL_TimerOn(this%GENSTATE,"kloopreshape") + !!do k=1, lm + !! p_src_3d(k,:,:) = qin_3d(:,:,k) + !!end do + !!call MAPL_TimerOff(this%GENSTATE,"kloopreshape") + + call MAPL_TimerOn(this%GENSTATE,"reshape") + p_src_3d = reshape(qin_3d,shape(p_src_3d),order=[2,3,1]) + call MAPL_TimerOff(this%GENSTATE,"reshape") call MAPL_TimerOn(this%GENSTATE,"3d_regrid") -!! call this%regridder%regrid(p_src_3d,p_dst_3d,_RC) call ESMF_FieldRegrid (new_src_field, new_dst_field, this%regridder%route_handle, _RC) call MAPL_TimerOff(this%GENSTATE,"3d_regrid") - - call MAPL_TimerOn(this%GENSTATE,"ESMF_FieldRedist") - ! ... - !! call ESMF_FieldRegrid(field_ds_3d, field_chunk_3d, this%RH, _RC) - call ESMF_FieldRedist(new_dst_field, field_chunk_3d, this%RH, _RC) - - call MAPL_TimerOff(this%GENSTATE,"ESMF_FieldRedist") - + call MAPL_TimerOn(this%GENSTATE,"FieldRedist") + call ESMF_FieldRedist (new_dst_field, field_chunk_3d, this%RH, _RC) + call MAPL_TimerOff(this%GENSTATE,"FieldRedist") call MAPL_TimerOn(this%GENSTATE,"gahterv") if (this%level_by_level) then @@ -678,13 +676,6 @@ subroutine append_file(this,current_time,rc) end if call MAPL_TimerOff(this%GENSTATE,"gahterv") -! if (mapl_am_i_root()) then -! do j=1, nx_sum, 500 -! write(6,*) 'p_rt_3d', p_rt_3d(:,j) -! end do -! end if -! if (mapl_am_i_root()) write(6,*) 'regrid + gatherV in 3D' - call MAPL_TimerOn(this%GENSTATE,"put3D") if (mapl_am_i_root()) then nz=size(p_rt_3d,1); nx=size(p_rt_3d,2) @@ -703,11 +694,16 @@ subroutine append_file(this,current_time,rc) call iter%next() end do + call MAPL_TimerOff(this%GENSTATE,"Full_regrid") + + + call MAPL_TimerOn(this%GENSTATE,"FieldDestroy") call ESMF_FieldDestroy(field_ds_2d, noGarbage=.true., _RC) call ESMF_FieldDestroy(field_chunk_2d, noGarbage=.true., _RC) call ESMF_FieldDestroy(field_chunk_3d, noGarbage=.true., _RC) call ESMF_FieldDestroy(new_dst_field, noGarbage=.true., _RC) call ESMF_FieldDestroy(new_src_field, noGarbage=.true., _RC) + call MAPL_TimerOff(this%GENSTATE,"FieldDestroy") _RETURN(_SUCCESS) end subroutine append_file From 99435930365430dcd0ed9dc0577c1bd7cef78755 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Tue, 18 Jun 2024 07:17:28 -0600 Subject: [PATCH 19/21] add MPI_Barrier --- .../Sampler/MAPL_StationSamplerMod.F90 | 21 +++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index a9e0dbee574d..7b5928eb8be4 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -609,7 +609,7 @@ subroutine append_file(this,current_time,rc) call MAPL_TimerOff(this%GENSTATE,"FieldCreate") - call MAPL_TimerOn(this%GENSTATE,"Full_regrid") +!! call MAPL_TimerOn(this%GENSTATE,"Full_regrid") iter = this%items%begin() do while (iter /= this%items%end()) @@ -654,11 +654,20 @@ subroutine append_file(this,current_time,rc) call ESMF_FieldRegrid (new_src_field, new_dst_field, this%regridder%route_handle, _RC) call MAPL_TimerOff(this%GENSTATE,"3d_regrid") + + + call MPI_Barrier(mpic,ierr) + _VERIFY (ierr) call MAPL_TimerOn(this%GENSTATE,"FieldRedist") call ESMF_FieldRedist (new_dst_field, field_chunk_3d, this%RH, _RC) + call MPI_Barrier(mpic,ierr) + _VERIFY (ierr) call MAPL_TimerOff(this%GENSTATE,"FieldRedist") - call MAPL_TimerOn(this%GENSTATE,"gahterv") + + call MPI_Barrier(mpic,ierr) + _VERIFY (ierr) + call MAPL_TimerOn(this%GENSTATE,"gatherv") if (this%level_by_level) then ! p_chunk_3d (lm, nx) allocate (p_dst_t, source = reshape(p_chunk_3d, [size(p_chunk_3d,2),size(p_chunk_3d,1)], order=[2,1])) @@ -666,6 +675,7 @@ subroutine append_file(this,current_time,rc) call MPI_gatherv ( p_dst_t(1,k), nsend, MPI_REAL, & p_rt_3d_aux(1,k), recvcount, displs, MPI_REAL,& iroot, mpic, ierr ) + _VERIFY (ierr) end do deallocate(p_dst_t) p_rt_3d = reshape(p_rt_3d_aux, shape(p_rt_3d), order=[2,1]) @@ -674,7 +684,10 @@ subroutine append_file(this,current_time,rc) p_rt_3d, recvcount_v, displs_v, MPI_REAL,& iroot, mpic, ierr ) end if - call MAPL_TimerOff(this%GENSTATE,"gahterv") + call MPI_Barrier(mpic,ierr) + _VERIFY (ierr) + call MAPL_TimerOff(this%GENSTATE,"gatherv") + call MAPL_TimerOn(this%GENSTATE,"put3D") if (mapl_am_i_root()) then @@ -694,7 +707,7 @@ subroutine append_file(this,current_time,rc) call iter%next() end do - call MAPL_TimerOff(this%GENSTATE,"Full_regrid") +!! call MAPL_TimerOff(this%GENSTATE,"Full_regrid") call MAPL_TimerOn(this%GENSTATE,"FieldDestroy") From 2b335c5682c143c47199f4f93237aed0dcb98b61 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Wed, 19 Jun 2024 21:52:39 -0400 Subject: [PATCH 20/21] change default for 3d gather to true --- gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 | 3 ++- gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 b/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 index c5f59b10ceec..f19e204d9c04 100644 --- a/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_GeosatMaskMod_smod.F90 @@ -5,13 +5,14 @@ implicit none contains -module function MaskSamplerGeosat_from_config(config,string,clock,rc) result(mask) +module function MaskSamplerGeosat_from_config(config,string,clock,GENSTATE,rc) result(mask) use BinIOMod use pflogger, only : Logger, logging type(MaskSamplerGeosat) :: mask type(ESMF_Config), intent(inout) :: config character(len=*), intent(in) :: string type(ESMF_Clock), intent(in) :: clock + type(MAPL_MetaComp), pointer, intent(in), optional :: GENSTATE integer, optional, intent(out) :: rc type(ESMF_Time) :: currTime diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index 7b5928eb8be4..fb88060e4dda 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -308,7 +308,7 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, GENSTATE, rc ! init ofile sampler%ofile='' sampler%obs_written=0 - sampler%level_by_level = .true. + sampler%level_by_level = .false. ! .true. _RETURN(_SUCCESS) end function new_StationSampler_readfile From 31e69bb8f23924c30742869c84c0d4a1aae5f4fe Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Fri, 21 Jun 2024 15:41:44 -0600 Subject: [PATCH 21/21] code clean up --- .../Sampler/MAPL_StationSamplerMod.F90 | 32 +++---------------- .../History/Sampler/MAPL_TrajectoryMod.F90 | 2 +- 2 files changed, 5 insertions(+), 29 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 index fb88060e4dda..8a94d5ef5665 100644 --- a/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 +++ b/gridcomps/History/Sampler/MAPL_StationSamplerMod.F90 @@ -304,11 +304,10 @@ function new_StationSampler_readfile (bundle, filename, nskip_line, GENSTATE, rc call ESMF_FieldBundleGet(bundle,grid=grid,_RC) sampler%LS_ds = sampler%LSF%create_locstream_on_proc(grid=grid,_RC) - ! ! init ofile sampler%ofile='' sampler%obs_written=0 - sampler%level_by_level = .false. ! .true. + sampler%level_by_level = .false. _RETURN(_SUCCESS) end function new_StationSampler_readfile @@ -391,7 +390,6 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) 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 @@ -415,7 +413,6 @@ subroutine add_metadata_route_handle (this,items,bundle,timeInfo,vdata,rc) call ESMF_FieldGet( chunk_field, localDE=0, farrayPtr=pt2, _RC ) pt1=0.0 pt2=0.0 -!! print*, 'shape(pt1 LS_ds)', shape(pt1) call ESMF_FieldRedistStore(src_field,chunk_field,this%RH,_RC) call ESMF_FieldDestroy(src_field,noGarbage=.true.,_RC) call ESMF_FieldDestroy(chunk_field,noGarbage=.true.,_RC) @@ -462,7 +459,6 @@ subroutine create_metadata_variable(this,vname,rc) case default _FAIL('unsupported rank') end select - ! v = variable(type=PFIO_REAL32,dimensions=trim(vdims),chunksizes=chunksizes) v = variable(type=PFIO_REAL32,dimensions=trim(vdims)) call v%add_attribute('units',trim(units)) @@ -526,12 +522,10 @@ subroutine append_file(this,current_time,rc) integer :: is, ie, ierr integer :: M, N, ip - this%obs_written=this%obs_written+1 !__ 1. put_var: time variable ! - ! rtimes = this%compute_time_for_current(current_time,_RC) ! rtimes: seconds since opening file if (mapl_am_i_root()) then call this%formatter%put_var('time',rtimes(1:1),& @@ -542,8 +536,9 @@ subroutine append_file(this,current_time,rc) !__ 2. regrid + put_var: ! ungridded_dim from src to dst [regrid] ! - ! caution about zero-sized array for MPI - ! redist + ! caution about zero-sized array for MPI + ! redist + ! nx_sum = this%nstation lm = this%vdata%lm call ESMF_VMGetCurrent(vm,_RC) @@ -608,9 +603,6 @@ subroutine append_file(this,current_time,rc) call MAPL_TimerOff(this%GENSTATE,"FieldCreate") - -!! call MAPL_TimerOn(this%GENSTATE,"Full_regrid") - iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() @@ -639,13 +631,6 @@ subroutine append_file(this,current_time,rc) ! call ESMF_FieldGet(src_field,localDE=0,farrayptr=qin_3d,_RC) - !! -- did not improve performance - !!call MAPL_TimerOn(this%GENSTATE,"kloopreshape") - !!do k=1, lm - !! p_src_3d(k,:,:) = qin_3d(:,:,k) - !!end do - !!call MAPL_TimerOff(this%GENSTATE,"kloopreshape") - call MAPL_TimerOn(this%GENSTATE,"reshape") p_src_3d = reshape(qin_3d,shape(p_src_3d),order=[2,3,1]) call MAPL_TimerOff(this%GENSTATE,"reshape") @@ -654,8 +639,6 @@ subroutine append_file(this,current_time,rc) call ESMF_FieldRegrid (new_src_field, new_dst_field, this%regridder%route_handle, _RC) call MAPL_TimerOff(this%GENSTATE,"3d_regrid") - - call MPI_Barrier(mpic,ierr) _VERIFY (ierr) call MAPL_TimerOn(this%GENSTATE,"FieldRedist") @@ -665,8 +648,6 @@ subroutine append_file(this,current_time,rc) call MAPL_TimerOff(this%GENSTATE,"FieldRedist") - call MPI_Barrier(mpic,ierr) - _VERIFY (ierr) call MAPL_TimerOn(this%GENSTATE,"gatherv") if (this%level_by_level) then ! p_chunk_3d (lm, nx) @@ -684,8 +665,6 @@ subroutine append_file(this,current_time,rc) p_rt_3d, recvcount_v, displs_v, MPI_REAL,& iroot, mpic, ierr ) end if - call MPI_Barrier(mpic,ierr) - _VERIFY (ierr) call MAPL_TimerOff(this%GENSTATE,"gatherv") @@ -707,8 +686,6 @@ subroutine append_file(this,current_time,rc) call iter%next() end do -!! call MAPL_TimerOff(this%GENSTATE,"Full_regrid") - call MAPL_TimerOn(this%GENSTATE,"FieldDestroy") call ESMF_FieldDestroy(field_ds_2d, noGarbage=.true., _RC) @@ -935,7 +912,6 @@ subroutine CSV_read_line_with_CH_I_R(line, name, lon, lat, rc) endif _ASSERT (ios==0, 'read error') - k=index(line(j+1:), '.') if (k > 0) then read(line(j+1:), *, iostat=ios) lat diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 index d42d0bcbdf0a..49c2eff96b38 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 @@ -73,7 +73,7 @@ module HistoryTrajectoryMod integer :: obsfile_Ts_index ! for epoch integer :: obsfile_Te_index logical :: active ! case: when no obs. exist - logical :: level_by_level = .true. + 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)