From 5e27db3e0a6d74eb684a162645c6f5ae1c7310dd Mon Sep 17 00:00:00 2001 From: "Ruifang.Li" Date: Tue, 19 Sep 2023 20:48:21 +0000 Subject: [PATCH 1/2] Process PM CSV file and convert to PrepBUFR format. --- CMakeLists.txt | 1 + pm/CMakeLists.txt | 19 ++++++ pm/process_pm.f90 | 143 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 163 insertions(+) create mode 100644 pm/CMakeLists.txt create mode 100644 pm/process_pm.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 1d1922a..00cf6cd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -74,3 +74,4 @@ add_subdirectory(update_bc) add_subdirectory(ens_mean_recenter) add_subdirectory(FV3_ensmean_recenter) add_subdirectory(bufrsnd) +add_subdirectory(pm) diff --git a/pm/CMakeLists.txt b/pm/CMakeLists.txt new file mode 100644 index 0000000..82c614f --- /dev/null +++ b/pm/CMakeLists.txt @@ -0,0 +1,19 @@ +list(APPEND src_pm + process_pm.f90) + +add_executable(process_pm.exe ${src_pm}) +target_link_libraries(process_pm.exe PRIVATE ${PROJECT_NAME}::ncio + PRIVATE ${PROJECT_NAME}::esggrid_util + PRIVATE ${PROJECT_NAME}::pesglib + PRIVATE NetCDF::NetCDF_Fortran + PRIVATE gsi::gsi + PRIVATE bufr::bufr_4 + PRIVATE MPI::MPI_Fortran + PRIVATE MPI::MPI_C) + +install( + TARGETS process_pm.exe + EXPORT ${PROJECT_NAME}Exports + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}) diff --git a/pm/process_pm.f90 b/pm/process_pm.f90 new file mode 100644 index 0000000..1df68b3 --- /dev/null +++ b/pm/process_pm.f90 @@ -0,0 +1,143 @@ +program pmbufr + + implicit none + + integer, parameter :: mxmn=35, mxlv=255, maxcnt=5000, pm_limit=-15 + integer :: ireadmg,ireadsb,idate,nmsg,ntb,nsubset,nlvl + integer :: i,j,k,NARG,iret_hdr,iret_ob,interval + integer :: dhr,iret_hdr1,iret_ob1,cnt + integer :: unit_table,unit_out,unit_site,unit_var,nt,pm_len + + character(80):: hdstr='SID XOB YOB DHR TYP T29 SQN PROCN RPT CAT TYPO TSIG' + character(80):: obstr='TPHR QCIND COPOPM' + character(8) :: subset,c_sid + character(10):: analysis_time + character(20):: str,infile,outfile + character(20),allocatable,dimension(:)::sta_id, lat_str, lon_str, elv_str, pm_aqi_str, pm_measured_str, pm_str + + real*8,parameter:: bmiss=10e10_8 + real(8) :: hdr(mxmn),obs(mxmn,mxlv) + real(8) :: rstation_id, lat, lon, elv, pm + + logical :: ifexist + + namelist/setup/analysis_time,infile,outfile,cnt + + ! only c_sid*8 is encoded to bufr file since SID is 64 bit in bufrtable + equivalence(rstation_id,c_sid) + +! read namelist + inquire(file='namelist.pm', EXIST=ifexist ) + if(ifexist) then + open(15, file='namelist.pm') + read(15,setup) + close(15) + write(*,'(a45,a20,3x,2a20,i8)') 'analysis_time, infile, outfile, cnt: ' ,analysis_time, infile, outfile, cnt + write(*,*) + else + write(*,*) "Usage: analysis_time, infile, outfile, pm_count" + stop 123 + + endif + + if (cnt==0) then + write(*,*) 'No observation available, exit' + stop + end if + +! open and read pm2.5 ascii file + unit_var=11 + allocate(sta_id(cnt)) + allocate(lat_str(cnt)) + allocate(lon_str(cnt)) + allocate(elv_str(cnt)) + allocate(pm_aqi_str(cnt)) + allocate(pm_measured_str(cnt)) + allocate(pm_str(cnt)) + + open (unit = unit_var, file = infile,status = 'old') + do i = 1,cnt + read(unit_var,'(7a16)') sta_id(i), lat_str(i), lon_str(i), elv_str(i), pm_aqi_str(i), pm_measured_str(i), pm_str(i) + !write(*,'(a10,7a15)') 'obs ', trim(adjustl(sta_id(i))), trim(adjustl(lat_str(i))), trim(adjustl(lon_str(i))),trim(adjustl(elv_str(i))), & + ! trim(adjustl(pm_aqi_str(i))), trim(adjustl(pm_measured_str(i))), trim(adjustl(pm_str(i))) + !write(*,'(a10,7i15)') 'len ', len(trim(adjustl(sta_id(i)))), len(trim(adjustl(lat_str(i)))), len(trim(adjustl(lon_str(i)))), len(trim(adjustl(elv_str(i)))), & + ! len(trim(adjustl(pm_aqi_str(i)))), len(trim(adjustl(pm_measured_str(i)))), len(trim(adjustl(pm_str(i)))) + end do + + +! open output files and write bufr table to output files + + unit_out=12 + unit_table=13 + open(unit_table,file='pm.bufrtable') + open(unit_out,file=outfile,action='write',form='unformatted') + call openbf(unit_out,'OUT',unit_table) + + !call maxout(20000) ! 20k byte any message size + call datelen(10) + + ! unchaged variables for each profile at current time + subset="ANOWPM" + read(analysis_time,'(I10)') idate ! string to integer + + ! hdstr='SID XOB YOB DHR TYP T29 SQN PROCN RPT CAT TYPO TSIG' + ! obstr='TPHR QCIND COPOPM' + + do j=1, cnt + + hdr=bmiss + obs=bmiss; + + ! write out pm which is not empty as single character ", but has valule like "2.0 + pm_len=len(trim(adjustl(pm_str(j)))) + if (pm_len > 1) then + + str=adjustl(lat_str(j)) + read(str(2:),'(f15.6)') lat + + str=adjustl(lon_str(j)) + read(str(2:),'(f15.6)') lon + + str=adjustl(elv_str(j)) + read(str(2:),'(f15.6)') elv + + str=adjustl(pm_str(j)) + read(str(2:),'(f15.6)') pm + + str=adjustl(sta_id(j)) + c_sid=str(2:9) + !write(*,'(2a15,4f15.5)') 'Obs ', c_sid, lat, lon, elv, pm + hdr(1)=rstation_id; hdr(2)=lon; hdr(3)=lat; hdr(4)=0.0_8; hdr(5)=102; hdr(10)=6 ! Single level report + + + nlvl=1 + obs(1,nlvl)=-1.0_8 + + if (pm > 0.0) then + obs(2,nlvl)=0.0_8 + obs(3,nlvl)=pm*1e-9_8 ! "UG/M3" to KG/M3 + write(*,'(a25,i7,a12,3f10.2,f15.12)') 'goodObs: pm>0',j,c_sid,lat,lon,obs(2,nlvl),obs(3,nlvl) + else if (pm > pm_limit) then + obs(2,nlvl)=0.0_8 + obs(3,nlvl)=0.0_8 + write(*,'(a25,i7,a12,3f10.2,f15.12)') 'limitObs: -15 Date: Thu, 21 Sep 2023 03:38:24 +0000 Subject: [PATCH 2/2] Remove some comments from process_pm.f90 and several libs from CMakeLists.txt --- pm/CMakeLists.txt | 10 ++-------- pm/process_pm.f90 | 11 +++-------- 2 files changed, 5 insertions(+), 16 deletions(-) diff --git a/pm/CMakeLists.txt b/pm/CMakeLists.txt index 82c614f..0d4a639 100644 --- a/pm/CMakeLists.txt +++ b/pm/CMakeLists.txt @@ -2,14 +2,8 @@ list(APPEND src_pm process_pm.f90) add_executable(process_pm.exe ${src_pm}) -target_link_libraries(process_pm.exe PRIVATE ${PROJECT_NAME}::ncio - PRIVATE ${PROJECT_NAME}::esggrid_util - PRIVATE ${PROJECT_NAME}::pesglib - PRIVATE NetCDF::NetCDF_Fortran - PRIVATE gsi::gsi - PRIVATE bufr::bufr_4 - PRIVATE MPI::MPI_Fortran - PRIVATE MPI::MPI_C) +target_link_libraries(process_pm.exe PRIVATE bufr::bufr_4 + PRIVATE MPI::MPI_Fortran) install( TARGETS process_pm.exe diff --git a/pm/process_pm.f90 b/pm/process_pm.f90 index 1df68b3..39a6042 100644 --- a/pm/process_pm.f90 +++ b/pm/process_pm.f90 @@ -58,10 +58,6 @@ program pmbufr open (unit = unit_var, file = infile,status = 'old') do i = 1,cnt read(unit_var,'(7a16)') sta_id(i), lat_str(i), lon_str(i), elv_str(i), pm_aqi_str(i), pm_measured_str(i), pm_str(i) - !write(*,'(a10,7a15)') 'obs ', trim(adjustl(sta_id(i))), trim(adjustl(lat_str(i))), trim(adjustl(lon_str(i))),trim(adjustl(elv_str(i))), & - ! trim(adjustl(pm_aqi_str(i))), trim(adjustl(pm_measured_str(i))), trim(adjustl(pm_str(i))) - !write(*,'(a10,7i15)') 'len ', len(trim(adjustl(sta_id(i)))), len(trim(adjustl(lat_str(i)))), len(trim(adjustl(lon_str(i)))), len(trim(adjustl(elv_str(i)))), & - ! len(trim(adjustl(pm_aqi_str(i)))), len(trim(adjustl(pm_measured_str(i)))), len(trim(adjustl(pm_str(i)))) end do @@ -80,15 +76,13 @@ program pmbufr subset="ANOWPM" read(analysis_time,'(I10)') idate ! string to integer - ! hdstr='SID XOB YOB DHR TYP T29 SQN PROCN RPT CAT TYPO TSIG' - ! obstr='TPHR QCIND COPOPM' do j=1, cnt hdr=bmiss obs=bmiss; - ! write out pm which is not empty as single character ", but has valule like "2.0 + ! write out pm which is not empty as single character ", but has valule e.g. "2.0 pm_len=len(trim(adjustl(pm_str(j)))) if (pm_len > 1) then @@ -106,7 +100,6 @@ program pmbufr str=adjustl(sta_id(j)) c_sid=str(2:9) - !write(*,'(2a15,4f15.5)') 'Obs ', c_sid, lat, lon, elv, pm hdr(1)=rstation_id; hdr(2)=lon; hdr(3)=lat; hdr(4)=0.0_8; hdr(5)=102; hdr(10)=6 ! Single level report @@ -136,6 +129,8 @@ program pmbufr endif enddo + write(*,*) "Process PM is successful" + call closbf(unit_out) close(unit_var)