From d717332de27acd8cf323810dde4f1fcfb054d995 Mon Sep 17 00:00:00 2001 From: "Bo.Cui" Date: Fri, 5 May 2023 10:52:19 -0400 Subject: [PATCH 1/9] Update gfs_bufr.fd and tocsbufr.fd --- src/gfs_bufr.fd/CMakeLists.txt | 53 ---------- src/gfs_bufr.fd/bfrhdr.f | 0 src/gfs_bufr.fd/bfrize.f | 0 src/gfs_bufr.fd/buff.f | 0 src/gfs_bufr.fd/calpreciptype.f | 2 +- src/gfs_bufr.fd/calwxt_gfs_baldwin.f | 0 src/gfs_bufr.fd/calwxt_gfs_ramer.f | 0 src/gfs_bufr.fd/funcphys.f | 0 src/gfs_bufr.fd/gfsbufr.f | 0 src/gfs_bufr.fd/gslp.f | 0 src/gfs_bufr.fd/lcl.f | 0 src/gfs_bufr.fd/machine.f | 0 src/gfs_bufr.fd/makefile_module | 79 ++++++++++++++ src/gfs_bufr.fd/meteorg.f | 151 ++++++++++++++------------- src/gfs_bufr.fd/modstuff1.f | 0 src/gfs_bufr.fd/mstadb.f | 0 src/gfs_bufr.fd/newsig1.f | 0 src/gfs_bufr.fd/physcons.f | 0 src/gfs_bufr.fd/rsearch.f | 0 src/gfs_bufr.fd/svp.f | 0 src/gfs_bufr.fd/tdew.f | 0 src/gfs_bufr.fd/terp3.f | 0 src/gfs_bufr.fd/vintg.f | 0 src/tocsbufr.fd/CMakeLists.txt | 19 ---- src/tocsbufr.fd/makefile_module | 82 +++++++++++++++ src/tocsbufr.fd/tocsbufr.f | 0 26 files changed, 238 insertions(+), 148 deletions(-) delete mode 100644 src/gfs_bufr.fd/CMakeLists.txt mode change 100644 => 100755 src/gfs_bufr.fd/bfrhdr.f mode change 100644 => 100755 src/gfs_bufr.fd/bfrize.f mode change 100644 => 100755 src/gfs_bufr.fd/buff.f mode change 100644 => 100755 src/gfs_bufr.fd/calwxt_gfs_baldwin.f mode change 100644 => 100755 src/gfs_bufr.fd/calwxt_gfs_ramer.f mode change 100644 => 100755 src/gfs_bufr.fd/funcphys.f mode change 100644 => 100755 src/gfs_bufr.fd/gfsbufr.f mode change 100644 => 100755 src/gfs_bufr.fd/gslp.f mode change 100644 => 100755 src/gfs_bufr.fd/lcl.f mode change 100644 => 100755 src/gfs_bufr.fd/machine.f create mode 100755 src/gfs_bufr.fd/makefile_module mode change 100644 => 100755 src/gfs_bufr.fd/meteorg.f mode change 100644 => 100755 src/gfs_bufr.fd/modstuff1.f mode change 100644 => 100755 src/gfs_bufr.fd/mstadb.f mode change 100644 => 100755 src/gfs_bufr.fd/newsig1.f mode change 100644 => 100755 src/gfs_bufr.fd/physcons.f mode change 100644 => 100755 src/gfs_bufr.fd/rsearch.f mode change 100644 => 100755 src/gfs_bufr.fd/svp.f mode change 100644 => 100755 src/gfs_bufr.fd/tdew.f mode change 100644 => 100755 src/gfs_bufr.fd/terp3.f mode change 100644 => 100755 src/gfs_bufr.fd/vintg.f delete mode 100644 src/tocsbufr.fd/CMakeLists.txt create mode 100755 src/tocsbufr.fd/makefile_module mode change 100644 => 100755 src/tocsbufr.fd/tocsbufr.f diff --git a/src/gfs_bufr.fd/CMakeLists.txt b/src/gfs_bufr.fd/CMakeLists.txt deleted file mode 100644 index 7df490d0..00000000 --- a/src/gfs_bufr.fd/CMakeLists.txt +++ /dev/null @@ -1,53 +0,0 @@ -list(APPEND fortran_src - bfrhdr.f - bfrize.f - buff.f - #calwxt_gfs_baldwin.f - #calwxt_gfs_ramer.f - gfsbufr.f - lcl.f - meteorg.f - mstadb.f - newsig1.f - read_nemsio.f - #read_netcdf.f - read_netcdf_p.f - rsearch.f - svp.f - tdew.f - terp3.f - vintg.f -) - -list(APPEND fortran_src_free - calpreciptype.f - funcphys.f - gslp.f - machine.f - modstuff1.f - physcons.f -) - -if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -convert big_endian -fp-model source") - set_source_files_properties(${fortran_src_free} PROPERTIES COMPILE_FLAGS "-free") -elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fconvert=big-endian") - set_source_files_properties(${fortran_src_free} PROPERTIES COMPILE_FLAGS "-ffree-form") -endif() - -set(exe_name gfs_bufr.x) -add_executable(${exe_name} ${fortran_src} ${fortran_src_free}) -target_link_libraries(${exe_name} PRIVATE NetCDF::NetCDF_Fortran - bacio::bacio_4 - sigio::sigio - sp::sp_4 - w3emc::w3emc_4 - nemsio::nemsio - bufr::bufr_4) - -if(OpenMP_Fortran_FOUND) - target_link_libraries(${exe_name} PRIVATE OpenMP::OpenMP_Fortran) -endif() - -install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/gfs_bufr.fd/bfrhdr.f b/src/gfs_bufr.fd/bfrhdr.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/bfrize.f b/src/gfs_bufr.fd/bfrize.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/buff.f b/src/gfs_bufr.fd/buff.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/calpreciptype.f b/src/gfs_bufr.fd/calpreciptype.f index 23072313..bca669c4 100644 --- a/src/gfs_bufr.fd/calpreciptype.f +++ b/src/gfs_bufr.fd/calpreciptype.f @@ -77,7 +77,7 @@ SUBROUTINE CALPRECIPTYPE(kdt,nrcm,im,ix,lm,lp1,randomno, & ALLOCATE ( RH(LM),TD8(LM),TWET8(LM) ) ! Create look up table - call gfuncphys +! call gfuncphys time_vert = 0. time_ncep = 0. diff --git a/src/gfs_bufr.fd/calwxt_gfs_baldwin.f b/src/gfs_bufr.fd/calwxt_gfs_baldwin.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/calwxt_gfs_ramer.f b/src/gfs_bufr.fd/calwxt_gfs_ramer.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/funcphys.f b/src/gfs_bufr.fd/funcphys.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/gfsbufr.f b/src/gfs_bufr.fd/gfsbufr.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/gslp.f b/src/gfs_bufr.fd/gslp.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/lcl.f b/src/gfs_bufr.fd/lcl.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/machine.f b/src/gfs_bufr.fd/machine.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/makefile_module b/src/gfs_bufr.fd/makefile_module new file mode 100755 index 00000000..d9d5374a --- /dev/null +++ b/src/gfs_bufr.fd/makefile_module @@ -0,0 +1,79 @@ +##################################################################################### +# gfs_bufr using module compile standard +# # 11/08/2019 guang.ping.lou@noaa.gov: Create NetCDF version +# ##################################################################################### +# set -eux +# + +FC = $(myFC) $(myFCFLAGS) +CPP = $(myCPP) $(myCPPFLAGS) + +FFLAGS = -I$(NETCDF_INCLUDES) \ + -I$(NEMSIO_INC) \ + -I$(SIGIO_INC) \ + -I$(W3EMC_INC4) + +LIBS = -L$(NETCDF_LIBRARIES) -lnetcdff -lnetcdf \ + -L$(HDF5_LIBRARIES) -lhdf5_hl -lhdf5 -lz \ + $(NEMSIO_LIB) \ + $(W3EMC_LIB4) \ + $(W3NCO_LIB4) \ + $(BUFR_LIB4) \ + $(BACIO_LIB4) \ + $(SP_LIB4) \ + $(SIGIO_LIB) + +SRCM = gfsbufr.f +OBJS = physcons.o funcphys.o meteorg.o bfrhdr.o newsig1.o terp3.o\ + bfrize.o vintg.o buff.o rsearch.o \ + svp.o calpreciptype.o lcl.o mstadb.o tdew.o\ + machine.o gslp.o modstuff1.o read_nemsio.o read_netcdf_p.o + +CMD = ../../exec/gfs_bufr + +$(CMD): $(SRCM) $(OBJS) + $(FC) $(FFLAGS) $(SRCM) $(OBJS) $(LIBS) -o $(CMD) + +machine.o: machine.f + $(FC) $(FFLAGS) -free -c machine.f +physcons.o: physcons.f machine.o + $(FC) $(FFLAGS) -free -c physcons.f +funcphys.o: funcphys.f physcons.o + $(FC) $(FFLAGS) -free -c funcphys.f +gslp.o: gslp.f + $(FC) $(FFLAGS) -free -c gslp.f +modstuff1.o: modstuff1.f + $(FC) $(INC) $(FFLAGS) -free -c modstuff1.f +meteorg.o: meteorg.f physcons.o funcphys.o + $(FC) $(INC) $(FFLAGS) -c meteorg.f +read_netcdf_p.o: read_netcdf_p.f + $(FC) $(INC) $(FFLAGS) -c read_netcdf_p.f +read_nemsio.o: read_nemsio.f + $(FC) $(INC) $(FFLAGS) -c read_nemsio.f +bfrhdr.o: bfrhdr.f + $(FC) $(FFLAGS) -c bfrhdr.f +newsig1.o: newsig1.f + $(FC) $(FFLAGS) -c newsig1.f +terp3.o: terp3.f + $(FC) $(FFLAGS) -c terp3.f +bfrize.o: bfrize.f + $(FC) $(FFLAGS) -c bfrize.f +vintg.o: vintg.f + $(FC) $(FFLAGS) -c vintg.f +buff.o: buff.f + $(FC) $(FFLAGS) -c buff.f +rsearch.o: rsearch.f + $(FC) $(FFLAGS) -c rsearch.f +svp.o: svp.f + $(FC) $(FFLAGS) -c svp.f +calpreciptype.o: calpreciptype.f physcons.o funcphys.o + $(FC) $(FFLAGS) -FR -c calpreciptype.f +lcl.o: lcl.f + $(FC) $(FFLAGS) -c lcl.f +mstadb.o: mstadb.f + $(FC) $(FFLAGS) -c mstadb.f +tdew.o: tdew.f + $(FC) $(FFLAGS) -c tdew.f + +clean: + /bin/rm -f $(OBJS) *.mod gfs_bufr diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f old mode 100644 new mode 100755 index 6b7c2c7d..d777d673 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -6,15 +6,15 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . -! SUBPROGRAM: meteorg +! SUBPROGRAM: meteorg ! PRGMMR: HUALU PAN ORG: W/NMC23 DATE: 1999-07-21 ! ! ABSTRACT: Creates BUFR meteogram files for the AVN and MRF. ! ! PROGRAM HISTORY LOG: -! 1999-07-21 HUALU PAN -! 2007-02-02 FANGLIN YANG EXPAND FOR HYBRID COORDINATES USING SIGIO -! 2009-07-24 FANGLIN YANG CHANGE OUTPUT PRESSURE TO INTEGER-LAYER +! 1999-07-21 HUALU PAN +! 2007-02-02 FANGLIN YANG EXPAND FOR HYBRID COORDINATES USING SIGIO +! 2009-07-24 FANGLIN YANG CHANGE OUTPUT PRESSURE TO INTEGER-LAYER ! PRESSURE (line 290) ! CORRECT THE TEMPERATURE ADJUSTMENT (line 238) ! 2014-03-27 DANA CARLIS UNIFY CODE WITH GFS FORECAST MODEL PRECIP @@ -22,19 +22,20 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! 2016-09-27 HUIYA CHUANG MODIFY TO READ GFS NEMS OUTPUT ON GRID SPACE ! 2017-02-27 GUANG PING LOU CHANGE OUTPUT PRECIPITATION TO HOURLY AMOUNT ! TO 120 HOURS AND 3 HOURLY TO 180 HOURS. -! 2018-02-01 GUANG PING LOU INGEST FV3GFS NEMSIO ACCUMULATED PRECIPITATION +! 2018-02-01 GUANG PING LOU INGEST FV3GFS NEMSIO ACCUMULATED PRECIPITATION ! AND RECALCULATE HOURLY AND 3 HOURLY OUTPUT DEPENDING -! ON LOGICAL VALUE OF precip_accu. +! ON LOGICAL VALUE OF precip_accu. ! 2018-02-08 GUANG PING LOU ADDED READING IN AND USING DZDT AS VERTICAL VELOCITY ! 2018-02-16 GUANG PING LOU ADDED READING IN AND USING MODEL DELP AND DELZ ! 2018-02-21 GUANG PING LOU THIS VERSION IS BACKWARD COMPATIBLE TO GFS MODEL ! 2018-03-27 GUANG PING LOU CHANGE STATION ELEVATION CORRECTION LAPSE RATE FROM 0.01 TO 0.0065 -! 2018-03-28 GUANG PING LOU GENERALIZE TIME INTERVAL +! 2018-03-28 GUANG PING LOU GENERALIZE TIME INTERVAL ! 2019-07-08 GUANG PING LOU ADDED STATION CHARACTER IDS ! 2019-10-08 GUANG PING LOU MODIFY TO READ IN NetCDF FILES. RETAIN NEMSIO ! RELATED CALLS AND CLEAN UP THE CODE. ! 2020-04-24 GUANG PING LOU Clean up code and remove station height ! adjustment +! 2023-03-28 Bo Cui Fix compilation error with "-check all" for gfs_bufrsnd ! ! USAGE: CALL PROGRAM meteorg ! INPUT: @@ -43,28 +44,29 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! rlon(npoint) - longtitude ! istat(npoint) - station id ! elevstn(npoint) - station elevation (m) -! nf - forecast cycle -! fnsig - sigma file name -! idate(4) - date +! nf - forecast cycle +! fnsig - sigma file name +! idate(4) - date ! levs - input vertical layers -! kdim - sfc file dimension +! kdim - sfc file dimension ! -! OUTPUT: -! nfile - output data file channel -! jdate - date YYYYMMDDHH +! OUTPUT: +! nfile - output data file channel +! jdate - date YYYYMMDDHH ! ! ATTRIBUTES: -! LANGUAGE: +! LANGUAGE: ! MACHINE: IBM SP ! !$$$ use netcdf use nemsio_module - use sigio_module + use sigio_module use physcons use mersenne_twister +! use funcphys, only : gfuncphys use funcphys - implicit none + implicit none include 'mpif.h' type(nemsio_gfile) :: gfile type(nemsio_gfile) :: ffile @@ -93,7 +95,6 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, real,dimension(im,jm):: apcp, cpcp real,dimension(npoint,2+levs*3):: grids real,dimension(npoint) :: rlat,rlon,pmsl,ps,psn,elevstn - real,dimension(1) :: psone real,dimension(im*jm) :: dum1d,dum1d2 real,dimension(im,jm) :: gdlat, hgt, gdlon real,dimension(im,jm,15) :: dum2d @@ -189,7 +190,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, error=nf90_inq_varid(ncid, "lat", id_var) error=nf90_get_var(ncid, id_var, gdlat) !!end read NetCDF hearder info, read nemsio below if necessary - else + else call nemsio_open(gfile,trim(fnsig),'read',iret=iret) call nemsio_getfilehead(gfile,iret=iret @@ -219,14 +220,14 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, gdlon(i,j)=dum1d2((j-1)*im+i) end do end do - - endif !end read in nemsio hearder + + endif !end read in nemsio hearder if(debugprint) then - do k=1,levs+1 + do k=1,levs+1 print*,'vcoord(k,1)= ', k, vcoord(k,1) end do - do k=1,levs+1 + do k=1,levs+1 print*,'vcoord(k,2)= ', k, vcoord(k,2) end do print*,'sample lat= ',gdlat(im/5,jm/4) @@ -241,7 +242,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,1,VarName,hgt,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'surface hgt not found' - else + else VarName='hgt' LayName='sfc' call read_nemsio(gfile,im,jm,1,VarName,LayName,hgt, @@ -259,7 +260,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'surface pressure not found' - else + else VarName='pres' LayName='sfc' call read_nemsio(gfile,im,jm,1,VarName, @@ -276,7 +277,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,t3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'temp not found' - else + else VarName='tmp' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,t3d,error) @@ -295,7 +296,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,q3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'spfh not found' - else + else VarName='spfh' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,q3d,error) @@ -314,7 +315,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,uh,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'ugrd not found' - else + else VarName='ugrd' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,uh,error) @@ -333,7 +334,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,vh,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'vgrd not found' - else + else VarName='vgrd' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,vh,error) @@ -352,7 +353,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,omega3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'dzdt not found' - else + else VarName='dzdt' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName, @@ -372,7 +373,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'dpres not found' - else + else VarName='dpres' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName, @@ -391,7 +392,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, do j=1,jm do i=1,im pint(i,j,k)=vcoord(k,1) - + +vcoord(k,2)*pint(i,j,1) + + +vcoord(k,2)*pint(i,j,1) end do end do end do @@ -434,7 +435,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'delz not found' - else + else VarName='delz' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,delpz,error) @@ -522,7 +523,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !open T-nint below error=nf90_open(trim(fngrib2),nf90_nowrite,ncid2) if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) - else + else call nemsio_open(ffile,trim(fngrib),'read',iret=error) call nemsio_open(ffile2,trim(fngrib2),'read',iret=error) if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) @@ -534,7 +535,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,1,VarName,lwmask,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'lwmask not found' - else + else VarName='land' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName,lwmask,error) @@ -552,7 +553,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tmpsfc not found' - else + else VarName='tmp' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -570,7 +571,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tmp2m not found' - else + else VarName='tmp' LayName='2 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -589,7 +590,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'spfh2m not found' - else + else VarName='spfh' LayName='2 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -608,7 +609,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'ugrd10m not found' - else + else VarName='ugrd' LayName='10 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -624,7 +625,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'vgrd10m not found' - else + else VarName='vgrd' LayName='10 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -640,7 +641,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'soilt1 not found' - else + else VarName='tmp' LayName='0-10 cm down' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -659,7 +660,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'snod not found' - else + else VarName='snod' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -676,7 +677,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'lhtfl not found' - else + else VarName='lhtfl' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -707,7 +708,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'prate_ave not found' - else + else VarName='prate_ave' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -740,7 +741,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'cprat_ave not found' - else + else VarName='cprat_ave' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -753,7 +754,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, do j=1,jm do i=1,im dum2d(i,j,10)=(apcp(i,j)*fhour-cpcp(i,j)*ap)*3600.0 - & + & end do end do @@ -765,7 +766,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'weasd not found' - else + else VarName='weasd' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -781,7 +782,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,12),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avelcl not found' - else + else VarName='tcdc_ave' LayName='low cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -797,7 +798,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,13),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avemcl not found' - else + else VarName='tcdc_ave' LayName='mid cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -813,7 +814,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,14),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avehcl not found' - else + else VarName='tcdc_ave' LayName='high cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -845,7 +846,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, jdum=jjdum(np) else -! find nearest neighbor +! find nearest neighbor rdum=rlon(np) if(rdum<0.)rdum=rdum+360. @@ -863,7 +864,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, exit else if(landwater(np) == lwmask(i+1,j))then idum=i+1 - jdum=j ! 2 + jdum=j ! 2 exit else if(landwater(np) == lwmask(i-1,j))then idum=i-1 @@ -1002,7 +1003,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, grids(np,1)=hgt(idum,jdum) grids(np,2)=pint(idum,jdum,1) - + sfc(5,np)=dum2d(idum,jdum,1) sfc(6,np)=dum2d(idum,jdum,6) sfc(17,np)=dum2d(idum,jdum,8) @@ -1019,10 +1020,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, CC There may be cases where convective precip is greater than total precip CC due to rounding and interpolation errors, correct it here -G.P. Lou: - if(sfc(11,np) .gt. sfc(12,np)) sfc(11,np)=sfc(12,np) + if(sfc(11,np) .gt. sfc(12,np)) sfc(11,np)=sfc(12,np) do k=1,levs - grids(np,k+2)=t3d(idum,jdum,k) + grids(np,k+2)=t3d(idum,jdum,k) grids(np,k+2+levs)=q3d(idum,jdum,k) grids(np,k+2+2*levs)=omega3d(idum,jdum,k) gridu(np,k)=uh(idum,jdum,k) @@ -1031,10 +1032,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, z1(np,k)=zint(idum,jdum,k+1) !! p1(np,k)=0.5*(pint(idum,jdum,k)+pint(idum,jdum,k+1)) !! z1(np,k)=0.5*(zint(idum,jdum,k)+zint(idum,jdum,k+1)) - + end do - end do - + end do + print*,'finish finding nearest neighbor for each station' do np = 1, npoint @@ -1049,7 +1050,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !! if(recn_dzdt == 0 ) then !!DZDT do k = 1, levs do np = 1, npoint - omega(np,k) = grids(np,2+levs*2+k) + omega(np,k) = grids(np,2+levs*2+k) enddo enddo if(debugprint) @@ -1065,10 +1066,9 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! print *, "elevstn = ", elevstn(np) if(elevstn(np)==-999.) elevstn(np) = grids(np,1) psn(np) = ps(np) - psone = ps(np) call sigio_modpr(1,1,levs,nvcoord,idvc, & idsl,vcoord,iret, - & ps=psone*1000,pd=pd3(np,1:levs)) + & ps=psn(np)*1000,pd=pd3(np,1:levs)) grids(np,2) = log(psn(np)) if(np==11)print*,'station H,grud H,psn,ps,new pm', & elevstn(np),grids(np,1),psn(np),ps(np) @@ -1094,7 +1094,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & pmsl(np),zp(np,1:levs),zp2(1:2)) enddo print *, 'call gslp pmsl= ', (pmsl(np),np=1,20) - if(recn_delz == -9999) then + if(recn_delz == -9999) then print*, 'using calculated height ' else print*, 'using model height m' @@ -1106,11 +1106,12 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, endif print*,'finish computing MSLP' print*,'finish computing zp ', (zp(11,k),k=1,levs) - print*,'finish computing zp2(11-12) ', zp2(11),zp2(12) + print*,'finish computing zp2(1-2) ', zp2(1),zp2(2) ! ! prepare buffer data ! if(iope == 0) then + call gfuncphys do np = 1, npoint pi3(np,1)=psn(np)*1000 do k=1,levs @@ -1138,7 +1139,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! look for the layer above 500 mb for precip type computation ! if(pi3(np,k).ge.50000.) leveta = k - ppi = pi3(np,k) + ppi = pi3(np,k) t = grids(np,k+2) q = max(1.e-8,grids(np,2+k+levs)) u = gridu(np,k) @@ -1152,8 +1153,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (mod(k,2)>0) then data2((kk-1)*6+7) = p1(np,k) data2((kk-1)*6+8) = t - data2((kk-1)*6+9) = u - data2((kk-1)*6+10) = v + data2((kk-1)*6+9) = u + data2((kk-1)*6+10) = v data2((kk-1)*6+11) = q data2((kk-1)*6+12) = omega(np,k)*100. endif @@ -1163,16 +1164,16 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! process surface flux file fields ! !! data(8+nflx) = psfc * 100. ! SURFACE PRESSURE (PA) -!! data(7+nflx) = pmsl(np) +!! data(7+nflx) = pmsl(np) data2(8+nflx2) = psfc * 100. ! SURFACE PRESSURE (PA) - data2(7+nflx2) = pmsl(np) + data2(7+nflx2) = pmsl(np) !! dtemp = .0065 * (grids(np,1) - elevstn(np)) !! dtemp = .0100 * (grids(np,1) - elevstn(np)) !! sfc(37,np) = data(6+nflx) * .01 !! sfc(37,np) = data(7+nflx) * .01 -!! sfc(39,np) = zp2(2) !500 hPa height +!! sfc(39,np) = zp2(2) !500 hPa height sfc(37,np) = data2(7+nflx2) * .01 - sfc(39,np) = zp2(2) !500 hPa height + sfc(39,np) = zp2(2) !500 hPa height ! ! do height correction if there is no snow or if the temp is less than 0 ! G.P.LOU: @@ -1194,10 +1195,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! !G.P. Lou 20200501: !convert instantaneous surface latent heat net flux to surface -!evapolation 1 W m-2 = 0.0864 MJ m-2 day-1 +!evapolation 1 W m-2 = 0.0864 MJ m-2 day-1 ! and 1 mm day-1 = 2.45 MJ m-2 day-1 ! equivament to 0.0864/2.54 = 0.035265 -! equivament to 2.54/0.0864 = 28.3565 +! equivament to 2.54/0.0864 = 28.3565 if(debugprint) + print*,'evaporation (stn 000692)= ',sfc(17,np) !! data(9+nflx) = sfc(5,np) ! tsfc (K) @@ -1250,8 +1251,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if(sfc(12,np).gt.0.) then !check for precip then calc precip type do k = 1, leveta+1 - pp = p1(np,k) - ppi = pi3(np,k) + pp = p1(np,k) + ppi = pi3(np,k) t = grids(np,k+2) q = max(0.,grids(np,2+k+levs)) u = gridu(np,k) @@ -1269,7 +1270,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, gq0(1,k) = q prsl(1,k) = pp prsi(1,k)=ppi - phii(1,k)=zp(np,k) !height in meters + phii(1,k)=zp(np,k) !height in meters enddo ! Use GFS routine calpreciptype.f to calculate precip type xlat=rlat(np) diff --git a/src/gfs_bufr.fd/modstuff1.f b/src/gfs_bufr.fd/modstuff1.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/mstadb.f b/src/gfs_bufr.fd/mstadb.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/newsig1.f b/src/gfs_bufr.fd/newsig1.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/physcons.f b/src/gfs_bufr.fd/physcons.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/rsearch.f b/src/gfs_bufr.fd/rsearch.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/svp.f b/src/gfs_bufr.fd/svp.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/tdew.f b/src/gfs_bufr.fd/tdew.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/terp3.f b/src/gfs_bufr.fd/terp3.f old mode 100644 new mode 100755 diff --git a/src/gfs_bufr.fd/vintg.f b/src/gfs_bufr.fd/vintg.f old mode 100644 new mode 100755 diff --git a/src/tocsbufr.fd/CMakeLists.txt b/src/tocsbufr.fd/CMakeLists.txt deleted file mode 100644 index 5f55927b..00000000 --- a/src/tocsbufr.fd/CMakeLists.txt +++ /dev/null @@ -1,19 +0,0 @@ -list(APPEND fortran_src - tocsbufr.f -) - -if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -convert big_endian -fp-model source") -elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fconvert=big-endian") -endif() - -set(exe_name tocsbufr.x) -add_executable(${exe_name} ${fortran_src}) -target_link_libraries(${exe_name} PRIVATE bacio::bacio_4 - sigio::sigio - sp::sp_4 - w3emc::w3emc_4 - bufr::bufr_4) - -install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/tocsbufr.fd/makefile_module b/src/tocsbufr.fd/makefile_module new file mode 100755 index 00000000..b6310ff7 --- /dev/null +++ b/src/tocsbufr.fd/makefile_module @@ -0,0 +1,82 @@ +SHELL=/bin/sh +# +# This makefile was produced by /usr/bin/fmgen at 11:21:07 AM on 10/28/94 +# If it is invoked by the command line +# make -f makefile +# it will compile the fortran modules indicated by SRCS into the object +# modules indicated by OBJS and produce an executable named a.out. +# +# If it is invoked by the command line +# make -f makefile a.out.prof +# it will compile the fortran modules indicated by SRCS into the object +# modules indicated by OBJS and produce an executable which profiles +# named a.out.prof. +# +# To remove all the objects but leave the executables use the command line +# make -f makefile clean +# +# To remove everything but the source files use the command line +# make -f makefile clobber +# +# To remove the source files created by /usr/bin/fmgen and this makefile +# use the command line +# make -f makefile void +# +# The parameters SRCS and OBJS should not need to be changed. If, however, +# you need to add a new module add the name of the source module to the +# SRCS parameter and add the name of the resulting object file to the OBJS +# parameter. The new modules are not limited to fortran, but may be C, YACC, +# LEX, or CAL. An explicit rule will need to be added for PASCAL modules. +# +SRCS= tocsbufr.f + +OBJS= tocsbufr.o + +# Tunable parameters +# +# FC Name of the fortran compiling system to use +# LDFLAGS Flags to the loader +# LIBS List of libraries +# CMD Name of the executable +# PROFLIB Library needed for profiling +# +FC = $(myFC) +LDFLAGS = $(myFCFLAGS) +LIBS = $(W3EMC_LIB4) \ + $(W3NCO_LIB4) \ + $(BUFR_LIB4) \ + $(BACIO_LIB4) \ + $(SP_LIB4) \ + $(SIGIO_LIB4) +CMD = ../../exec/tocsbufr +PROFLIB = -lprof + +# To perform the default compilation, use the first line +# To compile with flowtracing turned on, use the second line +# To compile giving profile additonal information, use the third line +# WARNING: SIMULTANEOUSLY PROFILING AND FLOWTRACING IS NOT RECOMMENDED +FFLAGS = $(FFLAGSM) +#FFLAGS = -F +#FFLAGS = -Wf"-ez" + +# Lines from here on down should not need to be changed. They are the +# actual rules which make uses to build a.out. +# +all: $(CMD) + +$(CMD): $(OBJS) + $(FC) $(LDFLAGS) -o $(@) $(OBJS) $(LIBS) + +# Make the profiled version of the command and call it a.out.prof +# +$(CMD).prof: $(OBJS) + $(FC) $(LDFLAGS) -o $(@) $(OBJS) $(PROFLIB) $(LIBS) + +clean: + -rm -f $(OBJS) + +clobber: clean + -rm -f $(CMD) $(CMD).prof + +void: clobber + -rm -f $(SRCS) makefile diff --git a/src/tocsbufr.fd/tocsbufr.f b/src/tocsbufr.fd/tocsbufr.f old mode 100644 new mode 100755 From 337fc0242ac7566e19dfd19791d17615233e8df9 Mon Sep 17 00:00:00 2001 From: "Bo.Cui" Date: Fri, 5 May 2023 13:55:57 -0400 Subject: [PATCH 2/9] Revert "Update gfs_bufr.fd and tocsbufr.fd" This reverts commit d717332de27acd8cf323810dde4f1fcfb054d995. --- src/gfs_bufr.fd/CMakeLists.txt | 53 ++++++++++ src/gfs_bufr.fd/bfrhdr.f | 0 src/gfs_bufr.fd/bfrize.f | 0 src/gfs_bufr.fd/buff.f | 0 src/gfs_bufr.fd/calpreciptype.f | 2 +- src/gfs_bufr.fd/calwxt_gfs_baldwin.f | 0 src/gfs_bufr.fd/calwxt_gfs_ramer.f | 0 src/gfs_bufr.fd/funcphys.f | 0 src/gfs_bufr.fd/gfsbufr.f | 0 src/gfs_bufr.fd/gslp.f | 0 src/gfs_bufr.fd/lcl.f | 0 src/gfs_bufr.fd/machine.f | 0 src/gfs_bufr.fd/makefile_module | 79 -------------- src/gfs_bufr.fd/meteorg.f | 151 +++++++++++++-------------- src/gfs_bufr.fd/modstuff1.f | 0 src/gfs_bufr.fd/mstadb.f | 0 src/gfs_bufr.fd/newsig1.f | 0 src/gfs_bufr.fd/physcons.f | 0 src/gfs_bufr.fd/rsearch.f | 0 src/gfs_bufr.fd/svp.f | 0 src/gfs_bufr.fd/tdew.f | 0 src/gfs_bufr.fd/terp3.f | 0 src/gfs_bufr.fd/vintg.f | 0 src/tocsbufr.fd/CMakeLists.txt | 19 ++++ src/tocsbufr.fd/makefile_module | 82 --------------- src/tocsbufr.fd/tocsbufr.f | 0 26 files changed, 148 insertions(+), 238 deletions(-) create mode 100644 src/gfs_bufr.fd/CMakeLists.txt mode change 100755 => 100644 src/gfs_bufr.fd/bfrhdr.f mode change 100755 => 100644 src/gfs_bufr.fd/bfrize.f mode change 100755 => 100644 src/gfs_bufr.fd/buff.f mode change 100755 => 100644 src/gfs_bufr.fd/calwxt_gfs_baldwin.f mode change 100755 => 100644 src/gfs_bufr.fd/calwxt_gfs_ramer.f mode change 100755 => 100644 src/gfs_bufr.fd/funcphys.f mode change 100755 => 100644 src/gfs_bufr.fd/gfsbufr.f mode change 100755 => 100644 src/gfs_bufr.fd/gslp.f mode change 100755 => 100644 src/gfs_bufr.fd/lcl.f mode change 100755 => 100644 src/gfs_bufr.fd/machine.f delete mode 100755 src/gfs_bufr.fd/makefile_module mode change 100755 => 100644 src/gfs_bufr.fd/meteorg.f mode change 100755 => 100644 src/gfs_bufr.fd/modstuff1.f mode change 100755 => 100644 src/gfs_bufr.fd/mstadb.f mode change 100755 => 100644 src/gfs_bufr.fd/newsig1.f mode change 100755 => 100644 src/gfs_bufr.fd/physcons.f mode change 100755 => 100644 src/gfs_bufr.fd/rsearch.f mode change 100755 => 100644 src/gfs_bufr.fd/svp.f mode change 100755 => 100644 src/gfs_bufr.fd/tdew.f mode change 100755 => 100644 src/gfs_bufr.fd/terp3.f mode change 100755 => 100644 src/gfs_bufr.fd/vintg.f create mode 100644 src/tocsbufr.fd/CMakeLists.txt delete mode 100755 src/tocsbufr.fd/makefile_module mode change 100755 => 100644 src/tocsbufr.fd/tocsbufr.f diff --git a/src/gfs_bufr.fd/CMakeLists.txt b/src/gfs_bufr.fd/CMakeLists.txt new file mode 100644 index 00000000..7df490d0 --- /dev/null +++ b/src/gfs_bufr.fd/CMakeLists.txt @@ -0,0 +1,53 @@ +list(APPEND fortran_src + bfrhdr.f + bfrize.f + buff.f + #calwxt_gfs_baldwin.f + #calwxt_gfs_ramer.f + gfsbufr.f + lcl.f + meteorg.f + mstadb.f + newsig1.f + read_nemsio.f + #read_netcdf.f + read_netcdf_p.f + rsearch.f + svp.f + tdew.f + terp3.f + vintg.f +) + +list(APPEND fortran_src_free + calpreciptype.f + funcphys.f + gslp.f + machine.f + modstuff1.f + physcons.f +) + +if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -convert big_endian -fp-model source") + set_source_files_properties(${fortran_src_free} PROPERTIES COMPILE_FLAGS "-free") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fconvert=big-endian") + set_source_files_properties(${fortran_src_free} PROPERTIES COMPILE_FLAGS "-ffree-form") +endif() + +set(exe_name gfs_bufr.x) +add_executable(${exe_name} ${fortran_src} ${fortran_src_free}) +target_link_libraries(${exe_name} PRIVATE NetCDF::NetCDF_Fortran + bacio::bacio_4 + sigio::sigio + sp::sp_4 + w3emc::w3emc_4 + nemsio::nemsio + bufr::bufr_4) + +if(OpenMP_Fortran_FOUND) + target_link_libraries(${exe_name} PRIVATE OpenMP::OpenMP_Fortran) +endif() + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/gfs_bufr.fd/bfrhdr.f b/src/gfs_bufr.fd/bfrhdr.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/bfrize.f b/src/gfs_bufr.fd/bfrize.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/buff.f b/src/gfs_bufr.fd/buff.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/calpreciptype.f b/src/gfs_bufr.fd/calpreciptype.f index bca669c4..23072313 100644 --- a/src/gfs_bufr.fd/calpreciptype.f +++ b/src/gfs_bufr.fd/calpreciptype.f @@ -77,7 +77,7 @@ SUBROUTINE CALPRECIPTYPE(kdt,nrcm,im,ix,lm,lp1,randomno, & ALLOCATE ( RH(LM),TD8(LM),TWET8(LM) ) ! Create look up table -! call gfuncphys + call gfuncphys time_vert = 0. time_ncep = 0. diff --git a/src/gfs_bufr.fd/calwxt_gfs_baldwin.f b/src/gfs_bufr.fd/calwxt_gfs_baldwin.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/calwxt_gfs_ramer.f b/src/gfs_bufr.fd/calwxt_gfs_ramer.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/funcphys.f b/src/gfs_bufr.fd/funcphys.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/gfsbufr.f b/src/gfs_bufr.fd/gfsbufr.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/gslp.f b/src/gfs_bufr.fd/gslp.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/lcl.f b/src/gfs_bufr.fd/lcl.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/machine.f b/src/gfs_bufr.fd/machine.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/makefile_module b/src/gfs_bufr.fd/makefile_module deleted file mode 100755 index d9d5374a..00000000 --- a/src/gfs_bufr.fd/makefile_module +++ /dev/null @@ -1,79 +0,0 @@ -##################################################################################### -# gfs_bufr using module compile standard -# # 11/08/2019 guang.ping.lou@noaa.gov: Create NetCDF version -# ##################################################################################### -# set -eux -# - -FC = $(myFC) $(myFCFLAGS) -CPP = $(myCPP) $(myCPPFLAGS) - -FFLAGS = -I$(NETCDF_INCLUDES) \ - -I$(NEMSIO_INC) \ - -I$(SIGIO_INC) \ - -I$(W3EMC_INC4) - -LIBS = -L$(NETCDF_LIBRARIES) -lnetcdff -lnetcdf \ - -L$(HDF5_LIBRARIES) -lhdf5_hl -lhdf5 -lz \ - $(NEMSIO_LIB) \ - $(W3EMC_LIB4) \ - $(W3NCO_LIB4) \ - $(BUFR_LIB4) \ - $(BACIO_LIB4) \ - $(SP_LIB4) \ - $(SIGIO_LIB) - -SRCM = gfsbufr.f -OBJS = physcons.o funcphys.o meteorg.o bfrhdr.o newsig1.o terp3.o\ - bfrize.o vintg.o buff.o rsearch.o \ - svp.o calpreciptype.o lcl.o mstadb.o tdew.o\ - machine.o gslp.o modstuff1.o read_nemsio.o read_netcdf_p.o - -CMD = ../../exec/gfs_bufr - -$(CMD): $(SRCM) $(OBJS) - $(FC) $(FFLAGS) $(SRCM) $(OBJS) $(LIBS) -o $(CMD) - -machine.o: machine.f - $(FC) $(FFLAGS) -free -c machine.f -physcons.o: physcons.f machine.o - $(FC) $(FFLAGS) -free -c physcons.f -funcphys.o: funcphys.f physcons.o - $(FC) $(FFLAGS) -free -c funcphys.f -gslp.o: gslp.f - $(FC) $(FFLAGS) -free -c gslp.f -modstuff1.o: modstuff1.f - $(FC) $(INC) $(FFLAGS) -free -c modstuff1.f -meteorg.o: meteorg.f physcons.o funcphys.o - $(FC) $(INC) $(FFLAGS) -c meteorg.f -read_netcdf_p.o: read_netcdf_p.f - $(FC) $(INC) $(FFLAGS) -c read_netcdf_p.f -read_nemsio.o: read_nemsio.f - $(FC) $(INC) $(FFLAGS) -c read_nemsio.f -bfrhdr.o: bfrhdr.f - $(FC) $(FFLAGS) -c bfrhdr.f -newsig1.o: newsig1.f - $(FC) $(FFLAGS) -c newsig1.f -terp3.o: terp3.f - $(FC) $(FFLAGS) -c terp3.f -bfrize.o: bfrize.f - $(FC) $(FFLAGS) -c bfrize.f -vintg.o: vintg.f - $(FC) $(FFLAGS) -c vintg.f -buff.o: buff.f - $(FC) $(FFLAGS) -c buff.f -rsearch.o: rsearch.f - $(FC) $(FFLAGS) -c rsearch.f -svp.o: svp.f - $(FC) $(FFLAGS) -c svp.f -calpreciptype.o: calpreciptype.f physcons.o funcphys.o - $(FC) $(FFLAGS) -FR -c calpreciptype.f -lcl.o: lcl.f - $(FC) $(FFLAGS) -c lcl.f -mstadb.o: mstadb.f - $(FC) $(FFLAGS) -c mstadb.f -tdew.o: tdew.f - $(FC) $(FFLAGS) -c tdew.f - -clean: - /bin/rm -f $(OBJS) *.mod gfs_bufr diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f old mode 100755 new mode 100644 index d777d673..6b7c2c7d --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -6,15 +6,15 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . -! SUBPROGRAM: meteorg +! SUBPROGRAM: meteorg ! PRGMMR: HUALU PAN ORG: W/NMC23 DATE: 1999-07-21 ! ! ABSTRACT: Creates BUFR meteogram files for the AVN and MRF. ! ! PROGRAM HISTORY LOG: -! 1999-07-21 HUALU PAN -! 2007-02-02 FANGLIN YANG EXPAND FOR HYBRID COORDINATES USING SIGIO -! 2009-07-24 FANGLIN YANG CHANGE OUTPUT PRESSURE TO INTEGER-LAYER +! 1999-07-21 HUALU PAN +! 2007-02-02 FANGLIN YANG EXPAND FOR HYBRID COORDINATES USING SIGIO +! 2009-07-24 FANGLIN YANG CHANGE OUTPUT PRESSURE TO INTEGER-LAYER ! PRESSURE (line 290) ! CORRECT THE TEMPERATURE ADJUSTMENT (line 238) ! 2014-03-27 DANA CARLIS UNIFY CODE WITH GFS FORECAST MODEL PRECIP @@ -22,20 +22,19 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! 2016-09-27 HUIYA CHUANG MODIFY TO READ GFS NEMS OUTPUT ON GRID SPACE ! 2017-02-27 GUANG PING LOU CHANGE OUTPUT PRECIPITATION TO HOURLY AMOUNT ! TO 120 HOURS AND 3 HOURLY TO 180 HOURS. -! 2018-02-01 GUANG PING LOU INGEST FV3GFS NEMSIO ACCUMULATED PRECIPITATION +! 2018-02-01 GUANG PING LOU INGEST FV3GFS NEMSIO ACCUMULATED PRECIPITATION ! AND RECALCULATE HOURLY AND 3 HOURLY OUTPUT DEPENDING -! ON LOGICAL VALUE OF precip_accu. +! ON LOGICAL VALUE OF precip_accu. ! 2018-02-08 GUANG PING LOU ADDED READING IN AND USING DZDT AS VERTICAL VELOCITY ! 2018-02-16 GUANG PING LOU ADDED READING IN AND USING MODEL DELP AND DELZ ! 2018-02-21 GUANG PING LOU THIS VERSION IS BACKWARD COMPATIBLE TO GFS MODEL ! 2018-03-27 GUANG PING LOU CHANGE STATION ELEVATION CORRECTION LAPSE RATE FROM 0.01 TO 0.0065 -! 2018-03-28 GUANG PING LOU GENERALIZE TIME INTERVAL +! 2018-03-28 GUANG PING LOU GENERALIZE TIME INTERVAL ! 2019-07-08 GUANG PING LOU ADDED STATION CHARACTER IDS ! 2019-10-08 GUANG PING LOU MODIFY TO READ IN NetCDF FILES. RETAIN NEMSIO ! RELATED CALLS AND CLEAN UP THE CODE. ! 2020-04-24 GUANG PING LOU Clean up code and remove station height ! adjustment -! 2023-03-28 Bo Cui Fix compilation error with "-check all" for gfs_bufrsnd ! ! USAGE: CALL PROGRAM meteorg ! INPUT: @@ -44,29 +43,28 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! rlon(npoint) - longtitude ! istat(npoint) - station id ! elevstn(npoint) - station elevation (m) -! nf - forecast cycle -! fnsig - sigma file name -! idate(4) - date +! nf - forecast cycle +! fnsig - sigma file name +! idate(4) - date ! levs - input vertical layers -! kdim - sfc file dimension +! kdim - sfc file dimension ! -! OUTPUT: -! nfile - output data file channel -! jdate - date YYYYMMDDHH +! OUTPUT: +! nfile - output data file channel +! jdate - date YYYYMMDDHH ! ! ATTRIBUTES: -! LANGUAGE: +! LANGUAGE: ! MACHINE: IBM SP ! !$$$ use netcdf use nemsio_module - use sigio_module + use sigio_module use physcons use mersenne_twister -! use funcphys, only : gfuncphys use funcphys - implicit none + implicit none include 'mpif.h' type(nemsio_gfile) :: gfile type(nemsio_gfile) :: ffile @@ -95,6 +93,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, real,dimension(im,jm):: apcp, cpcp real,dimension(npoint,2+levs*3):: grids real,dimension(npoint) :: rlat,rlon,pmsl,ps,psn,elevstn + real,dimension(1) :: psone real,dimension(im*jm) :: dum1d,dum1d2 real,dimension(im,jm) :: gdlat, hgt, gdlon real,dimension(im,jm,15) :: dum2d @@ -190,7 +189,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, error=nf90_inq_varid(ncid, "lat", id_var) error=nf90_get_var(ncid, id_var, gdlat) !!end read NetCDF hearder info, read nemsio below if necessary - else + else call nemsio_open(gfile,trim(fnsig),'read',iret=iret) call nemsio_getfilehead(gfile,iret=iret @@ -220,14 +219,14 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, gdlon(i,j)=dum1d2((j-1)*im+i) end do end do - - endif !end read in nemsio hearder + + endif !end read in nemsio hearder if(debugprint) then - do k=1,levs+1 + do k=1,levs+1 print*,'vcoord(k,1)= ', k, vcoord(k,1) end do - do k=1,levs+1 + do k=1,levs+1 print*,'vcoord(k,2)= ', k, vcoord(k,2) end do print*,'sample lat= ',gdlat(im/5,jm/4) @@ -242,7 +241,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,1,VarName,hgt,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'surface hgt not found' - else + else VarName='hgt' LayName='sfc' call read_nemsio(gfile,im,jm,1,VarName,LayName,hgt, @@ -260,7 +259,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'surface pressure not found' - else + else VarName='pres' LayName='sfc' call read_nemsio(gfile,im,jm,1,VarName, @@ -277,7 +276,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,t3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'temp not found' - else + else VarName='tmp' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,t3d,error) @@ -296,7 +295,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,q3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'spfh not found' - else + else VarName='spfh' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,q3d,error) @@ -315,7 +314,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,uh,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'ugrd not found' - else + else VarName='ugrd' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,uh,error) @@ -334,7 +333,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,vh,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'vgrd not found' - else + else VarName='vgrd' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,vh,error) @@ -353,7 +352,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,omega3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'dzdt not found' - else + else VarName='dzdt' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName, @@ -373,7 +372,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'dpres not found' - else + else VarName='dpres' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName, @@ -392,7 +391,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, do j=1,jm do i=1,im pint(i,j,k)=vcoord(k,1) - + +vcoord(k,2)*pint(i,j,1) + + +vcoord(k,2)*pint(i,j,1) end do end do end do @@ -435,7 +434,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'delz not found' - else + else VarName='delz' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,delpz,error) @@ -523,7 +522,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !open T-nint below error=nf90_open(trim(fngrib2),nf90_nowrite,ncid2) if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) - else + else call nemsio_open(ffile,trim(fngrib),'read',iret=error) call nemsio_open(ffile2,trim(fngrib2),'read',iret=error) if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) @@ -535,7 +534,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,1,VarName,lwmask,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'lwmask not found' - else + else VarName='land' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName,lwmask,error) @@ -553,7 +552,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tmpsfc not found' - else + else VarName='tmp' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -571,7 +570,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tmp2m not found' - else + else VarName='tmp' LayName='2 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -590,7 +589,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'spfh2m not found' - else + else VarName='spfh' LayName='2 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -609,7 +608,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'ugrd10m not found' - else + else VarName='ugrd' LayName='10 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -625,7 +624,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'vgrd10m not found' - else + else VarName='vgrd' LayName='10 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -641,7 +640,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'soilt1 not found' - else + else VarName='tmp' LayName='0-10 cm down' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -660,7 +659,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'snod not found' - else + else VarName='snod' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -677,7 +676,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'lhtfl not found' - else + else VarName='lhtfl' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -708,7 +707,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'prate_ave not found' - else + else VarName='prate_ave' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -741,7 +740,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'cprat_ave not found' - else + else VarName='cprat_ave' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -754,7 +753,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, do j=1,jm do i=1,im dum2d(i,j,10)=(apcp(i,j)*fhour-cpcp(i,j)*ap)*3600.0 - & + & end do end do @@ -766,7 +765,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'weasd not found' - else + else VarName='weasd' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -782,7 +781,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,12),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avelcl not found' - else + else VarName='tcdc_ave' LayName='low cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -798,7 +797,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,13),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avemcl not found' - else + else VarName='tcdc_ave' LayName='mid cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -814,7 +813,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,14),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avehcl not found' - else + else VarName='tcdc_ave' LayName='high cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -846,7 +845,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, jdum=jjdum(np) else -! find nearest neighbor +! find nearest neighbor rdum=rlon(np) if(rdum<0.)rdum=rdum+360. @@ -864,7 +863,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, exit else if(landwater(np) == lwmask(i+1,j))then idum=i+1 - jdum=j ! 2 + jdum=j ! 2 exit else if(landwater(np) == lwmask(i-1,j))then idum=i-1 @@ -1003,7 +1002,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, grids(np,1)=hgt(idum,jdum) grids(np,2)=pint(idum,jdum,1) - + sfc(5,np)=dum2d(idum,jdum,1) sfc(6,np)=dum2d(idum,jdum,6) sfc(17,np)=dum2d(idum,jdum,8) @@ -1020,10 +1019,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, CC There may be cases where convective precip is greater than total precip CC due to rounding and interpolation errors, correct it here -G.P. Lou: - if(sfc(11,np) .gt. sfc(12,np)) sfc(11,np)=sfc(12,np) + if(sfc(11,np) .gt. sfc(12,np)) sfc(11,np)=sfc(12,np) do k=1,levs - grids(np,k+2)=t3d(idum,jdum,k) + grids(np,k+2)=t3d(idum,jdum,k) grids(np,k+2+levs)=q3d(idum,jdum,k) grids(np,k+2+2*levs)=omega3d(idum,jdum,k) gridu(np,k)=uh(idum,jdum,k) @@ -1032,10 +1031,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, z1(np,k)=zint(idum,jdum,k+1) !! p1(np,k)=0.5*(pint(idum,jdum,k)+pint(idum,jdum,k+1)) !! z1(np,k)=0.5*(zint(idum,jdum,k)+zint(idum,jdum,k+1)) - + end do - end do - + end do + print*,'finish finding nearest neighbor for each station' do np = 1, npoint @@ -1050,7 +1049,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !! if(recn_dzdt == 0 ) then !!DZDT do k = 1, levs do np = 1, npoint - omega(np,k) = grids(np,2+levs*2+k) + omega(np,k) = grids(np,2+levs*2+k) enddo enddo if(debugprint) @@ -1066,9 +1065,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! print *, "elevstn = ", elevstn(np) if(elevstn(np)==-999.) elevstn(np) = grids(np,1) psn(np) = ps(np) + psone = ps(np) call sigio_modpr(1,1,levs,nvcoord,idvc, & idsl,vcoord,iret, - & ps=psn(np)*1000,pd=pd3(np,1:levs)) + & ps=psone*1000,pd=pd3(np,1:levs)) grids(np,2) = log(psn(np)) if(np==11)print*,'station H,grud H,psn,ps,new pm', & elevstn(np),grids(np,1),psn(np),ps(np) @@ -1094,7 +1094,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & pmsl(np),zp(np,1:levs),zp2(1:2)) enddo print *, 'call gslp pmsl= ', (pmsl(np),np=1,20) - if(recn_delz == -9999) then + if(recn_delz == -9999) then print*, 'using calculated height ' else print*, 'using model height m' @@ -1106,12 +1106,11 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, endif print*,'finish computing MSLP' print*,'finish computing zp ', (zp(11,k),k=1,levs) - print*,'finish computing zp2(1-2) ', zp2(1),zp2(2) + print*,'finish computing zp2(11-12) ', zp2(11),zp2(12) ! ! prepare buffer data ! if(iope == 0) then - call gfuncphys do np = 1, npoint pi3(np,1)=psn(np)*1000 do k=1,levs @@ -1139,7 +1138,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! look for the layer above 500 mb for precip type computation ! if(pi3(np,k).ge.50000.) leveta = k - ppi = pi3(np,k) + ppi = pi3(np,k) t = grids(np,k+2) q = max(1.e-8,grids(np,2+k+levs)) u = gridu(np,k) @@ -1153,8 +1152,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (mod(k,2)>0) then data2((kk-1)*6+7) = p1(np,k) data2((kk-1)*6+8) = t - data2((kk-1)*6+9) = u - data2((kk-1)*6+10) = v + data2((kk-1)*6+9) = u + data2((kk-1)*6+10) = v data2((kk-1)*6+11) = q data2((kk-1)*6+12) = omega(np,k)*100. endif @@ -1164,16 +1163,16 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! process surface flux file fields ! !! data(8+nflx) = psfc * 100. ! SURFACE PRESSURE (PA) -!! data(7+nflx) = pmsl(np) +!! data(7+nflx) = pmsl(np) data2(8+nflx2) = psfc * 100. ! SURFACE PRESSURE (PA) - data2(7+nflx2) = pmsl(np) + data2(7+nflx2) = pmsl(np) !! dtemp = .0065 * (grids(np,1) - elevstn(np)) !! dtemp = .0100 * (grids(np,1) - elevstn(np)) !! sfc(37,np) = data(6+nflx) * .01 !! sfc(37,np) = data(7+nflx) * .01 -!! sfc(39,np) = zp2(2) !500 hPa height +!! sfc(39,np) = zp2(2) !500 hPa height sfc(37,np) = data2(7+nflx2) * .01 - sfc(39,np) = zp2(2) !500 hPa height + sfc(39,np) = zp2(2) !500 hPa height ! ! do height correction if there is no snow or if the temp is less than 0 ! G.P.LOU: @@ -1195,10 +1194,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! !G.P. Lou 20200501: !convert instantaneous surface latent heat net flux to surface -!evapolation 1 W m-2 = 0.0864 MJ m-2 day-1 +!evapolation 1 W m-2 = 0.0864 MJ m-2 day-1 ! and 1 mm day-1 = 2.45 MJ m-2 day-1 ! equivament to 0.0864/2.54 = 0.035265 -! equivament to 2.54/0.0864 = 28.3565 +! equivament to 2.54/0.0864 = 28.3565 if(debugprint) + print*,'evaporation (stn 000692)= ',sfc(17,np) !! data(9+nflx) = sfc(5,np) ! tsfc (K) @@ -1251,8 +1250,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if(sfc(12,np).gt.0.) then !check for precip then calc precip type do k = 1, leveta+1 - pp = p1(np,k) - ppi = pi3(np,k) + pp = p1(np,k) + ppi = pi3(np,k) t = grids(np,k+2) q = max(0.,grids(np,2+k+levs)) u = gridu(np,k) @@ -1270,7 +1269,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, gq0(1,k) = q prsl(1,k) = pp prsi(1,k)=ppi - phii(1,k)=zp(np,k) !height in meters + phii(1,k)=zp(np,k) !height in meters enddo ! Use GFS routine calpreciptype.f to calculate precip type xlat=rlat(np) diff --git a/src/gfs_bufr.fd/modstuff1.f b/src/gfs_bufr.fd/modstuff1.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/mstadb.f b/src/gfs_bufr.fd/mstadb.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/newsig1.f b/src/gfs_bufr.fd/newsig1.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/physcons.f b/src/gfs_bufr.fd/physcons.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/rsearch.f b/src/gfs_bufr.fd/rsearch.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/svp.f b/src/gfs_bufr.fd/svp.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/tdew.f b/src/gfs_bufr.fd/tdew.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/terp3.f b/src/gfs_bufr.fd/terp3.f old mode 100755 new mode 100644 diff --git a/src/gfs_bufr.fd/vintg.f b/src/gfs_bufr.fd/vintg.f old mode 100755 new mode 100644 diff --git a/src/tocsbufr.fd/CMakeLists.txt b/src/tocsbufr.fd/CMakeLists.txt new file mode 100644 index 00000000..5f55927b --- /dev/null +++ b/src/tocsbufr.fd/CMakeLists.txt @@ -0,0 +1,19 @@ +list(APPEND fortran_src + tocsbufr.f +) + +if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -convert big_endian -fp-model source") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fconvert=big-endian") +endif() + +set(exe_name tocsbufr.x) +add_executable(${exe_name} ${fortran_src}) +target_link_libraries(${exe_name} PRIVATE bacio::bacio_4 + sigio::sigio + sp::sp_4 + w3emc::w3emc_4 + bufr::bufr_4) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/tocsbufr.fd/makefile_module b/src/tocsbufr.fd/makefile_module deleted file mode 100755 index b6310ff7..00000000 --- a/src/tocsbufr.fd/makefile_module +++ /dev/null @@ -1,82 +0,0 @@ -SHELL=/bin/sh -# -# This makefile was produced by /usr/bin/fmgen at 11:21:07 AM on 10/28/94 -# If it is invoked by the command line -# make -f makefile -# it will compile the fortran modules indicated by SRCS into the object -# modules indicated by OBJS and produce an executable named a.out. -# -# If it is invoked by the command line -# make -f makefile a.out.prof -# it will compile the fortran modules indicated by SRCS into the object -# modules indicated by OBJS and produce an executable which profiles -# named a.out.prof. -# -# To remove all the objects but leave the executables use the command line -# make -f makefile clean -# -# To remove everything but the source files use the command line -# make -f makefile clobber -# -# To remove the source files created by /usr/bin/fmgen and this makefile -# use the command line -# make -f makefile void -# -# The parameters SRCS and OBJS should not need to be changed. If, however, -# you need to add a new module add the name of the source module to the -# SRCS parameter and add the name of the resulting object file to the OBJS -# parameter. The new modules are not limited to fortran, but may be C, YACC, -# LEX, or CAL. An explicit rule will need to be added for PASCAL modules. -# -SRCS= tocsbufr.f - -OBJS= tocsbufr.o - -# Tunable parameters -# -# FC Name of the fortran compiling system to use -# LDFLAGS Flags to the loader -# LIBS List of libraries -# CMD Name of the executable -# PROFLIB Library needed for profiling -# -FC = $(myFC) -LDFLAGS = $(myFCFLAGS) -LIBS = $(W3EMC_LIB4) \ - $(W3NCO_LIB4) \ - $(BUFR_LIB4) \ - $(BACIO_LIB4) \ - $(SP_LIB4) \ - $(SIGIO_LIB4) -CMD = ../../exec/tocsbufr -PROFLIB = -lprof - -# To perform the default compilation, use the first line -# To compile with flowtracing turned on, use the second line -# To compile giving profile additonal information, use the third line -# WARNING: SIMULTANEOUSLY PROFILING AND FLOWTRACING IS NOT RECOMMENDED -FFLAGS = $(FFLAGSM) -#FFLAGS = -F -#FFLAGS = -Wf"-ez" - -# Lines from here on down should not need to be changed. They are the -# actual rules which make uses to build a.out. -# -all: $(CMD) - -$(CMD): $(OBJS) - $(FC) $(LDFLAGS) -o $(@) $(OBJS) $(LIBS) - -# Make the profiled version of the command and call it a.out.prof -# -$(CMD).prof: $(OBJS) - $(FC) $(LDFLAGS) -o $(@) $(OBJS) $(PROFLIB) $(LIBS) - -clean: - -rm -f $(OBJS) - -clobber: clean - -rm -f $(CMD) $(CMD).prof - -void: clobber - -rm -f $(SRCS) makefile diff --git a/src/tocsbufr.fd/tocsbufr.f b/src/tocsbufr.fd/tocsbufr.f old mode 100755 new mode 100644 From f1a65dcb719d875647d1e09900b58b658df29b46 Mon Sep 17 00:00:00 2001 From: "Bo.Cui" Date: Fri, 5 May 2023 14:03:02 -0400 Subject: [PATCH 3/9] update gfs_bufr the lastest version with bugzilla fix --- src/gfs_bufr.fd/calpreciptype.f | 2 +- src/gfs_bufr.fd/meteorg.f | 151 ++++++++++++++++---------------- 2 files changed, 77 insertions(+), 76 deletions(-) diff --git a/src/gfs_bufr.fd/calpreciptype.f b/src/gfs_bufr.fd/calpreciptype.f index 23072313..bca669c4 100644 --- a/src/gfs_bufr.fd/calpreciptype.f +++ b/src/gfs_bufr.fd/calpreciptype.f @@ -77,7 +77,7 @@ SUBROUTINE CALPRECIPTYPE(kdt,nrcm,im,ix,lm,lp1,randomno, & ALLOCATE ( RH(LM),TD8(LM),TWET8(LM) ) ! Create look up table - call gfuncphys +! call gfuncphys time_vert = 0. time_ncep = 0. diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f index 6b7c2c7d..d777d673 100644 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -6,15 +6,15 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . -! SUBPROGRAM: meteorg +! SUBPROGRAM: meteorg ! PRGMMR: HUALU PAN ORG: W/NMC23 DATE: 1999-07-21 ! ! ABSTRACT: Creates BUFR meteogram files for the AVN and MRF. ! ! PROGRAM HISTORY LOG: -! 1999-07-21 HUALU PAN -! 2007-02-02 FANGLIN YANG EXPAND FOR HYBRID COORDINATES USING SIGIO -! 2009-07-24 FANGLIN YANG CHANGE OUTPUT PRESSURE TO INTEGER-LAYER +! 1999-07-21 HUALU PAN +! 2007-02-02 FANGLIN YANG EXPAND FOR HYBRID COORDINATES USING SIGIO +! 2009-07-24 FANGLIN YANG CHANGE OUTPUT PRESSURE TO INTEGER-LAYER ! PRESSURE (line 290) ! CORRECT THE TEMPERATURE ADJUSTMENT (line 238) ! 2014-03-27 DANA CARLIS UNIFY CODE WITH GFS FORECAST MODEL PRECIP @@ -22,19 +22,20 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! 2016-09-27 HUIYA CHUANG MODIFY TO READ GFS NEMS OUTPUT ON GRID SPACE ! 2017-02-27 GUANG PING LOU CHANGE OUTPUT PRECIPITATION TO HOURLY AMOUNT ! TO 120 HOURS AND 3 HOURLY TO 180 HOURS. -! 2018-02-01 GUANG PING LOU INGEST FV3GFS NEMSIO ACCUMULATED PRECIPITATION +! 2018-02-01 GUANG PING LOU INGEST FV3GFS NEMSIO ACCUMULATED PRECIPITATION ! AND RECALCULATE HOURLY AND 3 HOURLY OUTPUT DEPENDING -! ON LOGICAL VALUE OF precip_accu. +! ON LOGICAL VALUE OF precip_accu. ! 2018-02-08 GUANG PING LOU ADDED READING IN AND USING DZDT AS VERTICAL VELOCITY ! 2018-02-16 GUANG PING LOU ADDED READING IN AND USING MODEL DELP AND DELZ ! 2018-02-21 GUANG PING LOU THIS VERSION IS BACKWARD COMPATIBLE TO GFS MODEL ! 2018-03-27 GUANG PING LOU CHANGE STATION ELEVATION CORRECTION LAPSE RATE FROM 0.01 TO 0.0065 -! 2018-03-28 GUANG PING LOU GENERALIZE TIME INTERVAL +! 2018-03-28 GUANG PING LOU GENERALIZE TIME INTERVAL ! 2019-07-08 GUANG PING LOU ADDED STATION CHARACTER IDS ! 2019-10-08 GUANG PING LOU MODIFY TO READ IN NetCDF FILES. RETAIN NEMSIO ! RELATED CALLS AND CLEAN UP THE CODE. ! 2020-04-24 GUANG PING LOU Clean up code and remove station height ! adjustment +! 2023-03-28 Bo Cui Fix compilation error with "-check all" for gfs_bufrsnd ! ! USAGE: CALL PROGRAM meteorg ! INPUT: @@ -43,28 +44,29 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! rlon(npoint) - longtitude ! istat(npoint) - station id ! elevstn(npoint) - station elevation (m) -! nf - forecast cycle -! fnsig - sigma file name -! idate(4) - date +! nf - forecast cycle +! fnsig - sigma file name +! idate(4) - date ! levs - input vertical layers -! kdim - sfc file dimension +! kdim - sfc file dimension ! -! OUTPUT: -! nfile - output data file channel -! jdate - date YYYYMMDDHH +! OUTPUT: +! nfile - output data file channel +! jdate - date YYYYMMDDHH ! ! ATTRIBUTES: -! LANGUAGE: +! LANGUAGE: ! MACHINE: IBM SP ! !$$$ use netcdf use nemsio_module - use sigio_module + use sigio_module use physcons use mersenne_twister +! use funcphys, only : gfuncphys use funcphys - implicit none + implicit none include 'mpif.h' type(nemsio_gfile) :: gfile type(nemsio_gfile) :: ffile @@ -93,7 +95,6 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, real,dimension(im,jm):: apcp, cpcp real,dimension(npoint,2+levs*3):: grids real,dimension(npoint) :: rlat,rlon,pmsl,ps,psn,elevstn - real,dimension(1) :: psone real,dimension(im*jm) :: dum1d,dum1d2 real,dimension(im,jm) :: gdlat, hgt, gdlon real,dimension(im,jm,15) :: dum2d @@ -189,7 +190,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, error=nf90_inq_varid(ncid, "lat", id_var) error=nf90_get_var(ncid, id_var, gdlat) !!end read NetCDF hearder info, read nemsio below if necessary - else + else call nemsio_open(gfile,trim(fnsig),'read',iret=iret) call nemsio_getfilehead(gfile,iret=iret @@ -219,14 +220,14 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, gdlon(i,j)=dum1d2((j-1)*im+i) end do end do - - endif !end read in nemsio hearder + + endif !end read in nemsio hearder if(debugprint) then - do k=1,levs+1 + do k=1,levs+1 print*,'vcoord(k,1)= ', k, vcoord(k,1) end do - do k=1,levs+1 + do k=1,levs+1 print*,'vcoord(k,2)= ', k, vcoord(k,2) end do print*,'sample lat= ',gdlat(im/5,jm/4) @@ -241,7 +242,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,1,VarName,hgt,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'surface hgt not found' - else + else VarName='hgt' LayName='sfc' call read_nemsio(gfile,im,jm,1,VarName,LayName,hgt, @@ -259,7 +260,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'surface pressure not found' - else + else VarName='pres' LayName='sfc' call read_nemsio(gfile,im,jm,1,VarName, @@ -276,7 +277,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,t3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'temp not found' - else + else VarName='tmp' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,t3d,error) @@ -295,7 +296,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,q3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'spfh not found' - else + else VarName='spfh' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,q3d,error) @@ -314,7 +315,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,uh,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'ugrd not found' - else + else VarName='ugrd' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,uh,error) @@ -333,7 +334,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,vh,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'vgrd not found' - else + else VarName='vgrd' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,vh,error) @@ -352,7 +353,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,omega3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'dzdt not found' - else + else VarName='dzdt' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName, @@ -372,7 +373,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'dpres not found' - else + else VarName='dpres' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName, @@ -391,7 +392,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, do j=1,jm do i=1,im pint(i,j,k)=vcoord(k,1) - + +vcoord(k,2)*pint(i,j,1) + + +vcoord(k,2)*pint(i,j,1) end do end do end do @@ -434,7 +435,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'delz not found' - else + else VarName='delz' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,delpz,error) @@ -522,7 +523,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !open T-nint below error=nf90_open(trim(fngrib2),nf90_nowrite,ncid2) if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) - else + else call nemsio_open(ffile,trim(fngrib),'read',iret=error) call nemsio_open(ffile2,trim(fngrib2),'read',iret=error) if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) @@ -534,7 +535,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,1,VarName,lwmask,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'lwmask not found' - else + else VarName='land' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName,lwmask,error) @@ -552,7 +553,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tmpsfc not found' - else + else VarName='tmp' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -570,7 +571,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tmp2m not found' - else + else VarName='tmp' LayName='2 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -589,7 +590,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'spfh2m not found' - else + else VarName='spfh' LayName='2 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -608,7 +609,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'ugrd10m not found' - else + else VarName='ugrd' LayName='10 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -624,7 +625,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'vgrd10m not found' - else + else VarName='vgrd' LayName='10 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -640,7 +641,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'soilt1 not found' - else + else VarName='tmp' LayName='0-10 cm down' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -659,7 +660,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'snod not found' - else + else VarName='snod' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -676,7 +677,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'lhtfl not found' - else + else VarName='lhtfl' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -707,7 +708,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'prate_ave not found' - else + else VarName='prate_ave' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -740,7 +741,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'cprat_ave not found' - else + else VarName='cprat_ave' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -753,7 +754,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, do j=1,jm do i=1,im dum2d(i,j,10)=(apcp(i,j)*fhour-cpcp(i,j)*ap)*3600.0 - & + & end do end do @@ -765,7 +766,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'weasd not found' - else + else VarName='weasd' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -781,7 +782,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,12),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avelcl not found' - else + else VarName='tcdc_ave' LayName='low cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -797,7 +798,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,13),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avemcl not found' - else + else VarName='tcdc_ave' LayName='mid cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -813,7 +814,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,14),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avehcl not found' - else + else VarName='tcdc_ave' LayName='high cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -845,7 +846,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, jdum=jjdum(np) else -! find nearest neighbor +! find nearest neighbor rdum=rlon(np) if(rdum<0.)rdum=rdum+360. @@ -863,7 +864,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, exit else if(landwater(np) == lwmask(i+1,j))then idum=i+1 - jdum=j ! 2 + jdum=j ! 2 exit else if(landwater(np) == lwmask(i-1,j))then idum=i-1 @@ -1002,7 +1003,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, grids(np,1)=hgt(idum,jdum) grids(np,2)=pint(idum,jdum,1) - + sfc(5,np)=dum2d(idum,jdum,1) sfc(6,np)=dum2d(idum,jdum,6) sfc(17,np)=dum2d(idum,jdum,8) @@ -1019,10 +1020,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, CC There may be cases where convective precip is greater than total precip CC due to rounding and interpolation errors, correct it here -G.P. Lou: - if(sfc(11,np) .gt. sfc(12,np)) sfc(11,np)=sfc(12,np) + if(sfc(11,np) .gt. sfc(12,np)) sfc(11,np)=sfc(12,np) do k=1,levs - grids(np,k+2)=t3d(idum,jdum,k) + grids(np,k+2)=t3d(idum,jdum,k) grids(np,k+2+levs)=q3d(idum,jdum,k) grids(np,k+2+2*levs)=omega3d(idum,jdum,k) gridu(np,k)=uh(idum,jdum,k) @@ -1031,10 +1032,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, z1(np,k)=zint(idum,jdum,k+1) !! p1(np,k)=0.5*(pint(idum,jdum,k)+pint(idum,jdum,k+1)) !! z1(np,k)=0.5*(zint(idum,jdum,k)+zint(idum,jdum,k+1)) - + end do - end do - + end do + print*,'finish finding nearest neighbor for each station' do np = 1, npoint @@ -1049,7 +1050,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !! if(recn_dzdt == 0 ) then !!DZDT do k = 1, levs do np = 1, npoint - omega(np,k) = grids(np,2+levs*2+k) + omega(np,k) = grids(np,2+levs*2+k) enddo enddo if(debugprint) @@ -1065,10 +1066,9 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! print *, "elevstn = ", elevstn(np) if(elevstn(np)==-999.) elevstn(np) = grids(np,1) psn(np) = ps(np) - psone = ps(np) call sigio_modpr(1,1,levs,nvcoord,idvc, & idsl,vcoord,iret, - & ps=psone*1000,pd=pd3(np,1:levs)) + & ps=psn(np)*1000,pd=pd3(np,1:levs)) grids(np,2) = log(psn(np)) if(np==11)print*,'station H,grud H,psn,ps,new pm', & elevstn(np),grids(np,1),psn(np),ps(np) @@ -1094,7 +1094,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & pmsl(np),zp(np,1:levs),zp2(1:2)) enddo print *, 'call gslp pmsl= ', (pmsl(np),np=1,20) - if(recn_delz == -9999) then + if(recn_delz == -9999) then print*, 'using calculated height ' else print*, 'using model height m' @@ -1106,11 +1106,12 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, endif print*,'finish computing MSLP' print*,'finish computing zp ', (zp(11,k),k=1,levs) - print*,'finish computing zp2(11-12) ', zp2(11),zp2(12) + print*,'finish computing zp2(1-2) ', zp2(1),zp2(2) ! ! prepare buffer data ! if(iope == 0) then + call gfuncphys do np = 1, npoint pi3(np,1)=psn(np)*1000 do k=1,levs @@ -1138,7 +1139,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! look for the layer above 500 mb for precip type computation ! if(pi3(np,k).ge.50000.) leveta = k - ppi = pi3(np,k) + ppi = pi3(np,k) t = grids(np,k+2) q = max(1.e-8,grids(np,2+k+levs)) u = gridu(np,k) @@ -1152,8 +1153,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (mod(k,2)>0) then data2((kk-1)*6+7) = p1(np,k) data2((kk-1)*6+8) = t - data2((kk-1)*6+9) = u - data2((kk-1)*6+10) = v + data2((kk-1)*6+9) = u + data2((kk-1)*6+10) = v data2((kk-1)*6+11) = q data2((kk-1)*6+12) = omega(np,k)*100. endif @@ -1163,16 +1164,16 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! process surface flux file fields ! !! data(8+nflx) = psfc * 100. ! SURFACE PRESSURE (PA) -!! data(7+nflx) = pmsl(np) +!! data(7+nflx) = pmsl(np) data2(8+nflx2) = psfc * 100. ! SURFACE PRESSURE (PA) - data2(7+nflx2) = pmsl(np) + data2(7+nflx2) = pmsl(np) !! dtemp = .0065 * (grids(np,1) - elevstn(np)) !! dtemp = .0100 * (grids(np,1) - elevstn(np)) !! sfc(37,np) = data(6+nflx) * .01 !! sfc(37,np) = data(7+nflx) * .01 -!! sfc(39,np) = zp2(2) !500 hPa height +!! sfc(39,np) = zp2(2) !500 hPa height sfc(37,np) = data2(7+nflx2) * .01 - sfc(39,np) = zp2(2) !500 hPa height + sfc(39,np) = zp2(2) !500 hPa height ! ! do height correction if there is no snow or if the temp is less than 0 ! G.P.LOU: @@ -1194,10 +1195,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! !G.P. Lou 20200501: !convert instantaneous surface latent heat net flux to surface -!evapolation 1 W m-2 = 0.0864 MJ m-2 day-1 +!evapolation 1 W m-2 = 0.0864 MJ m-2 day-1 ! and 1 mm day-1 = 2.45 MJ m-2 day-1 ! equivament to 0.0864/2.54 = 0.035265 -! equivament to 2.54/0.0864 = 28.3565 +! equivament to 2.54/0.0864 = 28.3565 if(debugprint) + print*,'evaporation (stn 000692)= ',sfc(17,np) !! data(9+nflx) = sfc(5,np) ! tsfc (K) @@ -1250,8 +1251,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if(sfc(12,np).gt.0.) then !check for precip then calc precip type do k = 1, leveta+1 - pp = p1(np,k) - ppi = pi3(np,k) + pp = p1(np,k) + ppi = pi3(np,k) t = grids(np,k+2) q = max(0.,grids(np,2+k+levs)) u = gridu(np,k) @@ -1269,7 +1270,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, gq0(1,k) = q prsl(1,k) = pp prsi(1,k)=ppi - phii(1,k)=zp(np,k) !height in meters + phii(1,k)=zp(np,k) !height in meters enddo ! Use GFS routine calpreciptype.f to calculate precip type xlat=rlat(np) From ae6fc823c54d2aa44bf208cb2ed8731a94f5790d Mon Sep 17 00:00:00 2001 From: BoCui-NOAA <53531984+BoCui-NOAA@users.noreply.github.com> Date: Fri, 5 May 2023 15:16:15 -0400 Subject: [PATCH 4/9] Update meteorg.f clean out extra whitespace --- src/gfs_bufr.fd/meteorg.f | 44 +++++++++++++++++++-------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f index d777d673..dd08e815 100644 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -6,13 +6,13 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . -! SUBPROGRAM: meteorg +! SUBPROGRAM: meteorg ! PRGMMR: HUALU PAN ORG: W/NMC23 DATE: 1999-07-21 ! ! ABSTRACT: Creates BUFR meteogram files for the AVN and MRF. ! ! PROGRAM HISTORY LOG: -! 1999-07-21 HUALU PAN +! 1999-07-21 HUALU PAN ! 2007-02-02 FANGLIN YANG EXPAND FOR HYBRID COORDINATES USING SIGIO ! 2009-07-24 FANGLIN YANG CHANGE OUTPUT PRESSURE TO INTEGER-LAYER ! PRESSURE (line 290) @@ -44,15 +44,15 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! rlon(npoint) - longtitude ! istat(npoint) - station id ! elevstn(npoint) - station elevation (m) -! nf - forecast cycle -! fnsig - sigma file name -! idate(4) - date +! nf - forecast cycle +! fnsig - sigma file name +! idate(4) - date ! levs - input vertical layers -! kdim - sfc file dimension +! kdim - sfc file dimension ! -! OUTPUT: +! OUTPUT: ! nfile - output data file channel -! jdate - date YYYYMMDDHH +! jdate - date YYYYMMDDHH ! ! ATTRIBUTES: ! LANGUAGE: @@ -221,10 +221,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, end do end do - endif !end read in nemsio hearder + endif !end read in nemsio hearder if(debugprint) then - do k=1,levs+1 + do k=1,levs+1 print*,'vcoord(k,1)= ', k, vcoord(k,1) end do do k=1,levs+1 @@ -754,7 +754,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, do j=1,jm do i=1,im dum2d(i,j,10)=(apcp(i,j)*fhour-cpcp(i,j)*ap)*3600.0 - & + & end do end do @@ -1003,7 +1003,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, grids(np,1)=hgt(idum,jdum) grids(np,2)=pint(idum,jdum,1) - + sfc(5,np)=dum2d(idum,jdum,1) sfc(6,np)=dum2d(idum,jdum,6) sfc(17,np)=dum2d(idum,jdum,8) @@ -1035,7 +1035,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, end do end do - + print*,'finish finding nearest neighbor for each station' do np = 1, npoint @@ -1106,7 +1106,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, endif print*,'finish computing MSLP' print*,'finish computing zp ', (zp(11,k),k=1,levs) - print*,'finish computing zp2(1-2) ', zp2(1),zp2(2) + print*,'finish computing zp2(1-2) ', zp2(1),zp2(2) ! ! prepare buffer data ! @@ -1139,7 +1139,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! look for the layer above 500 mb for precip type computation ! if(pi3(np,k).ge.50000.) leveta = k - ppi = pi3(np,k) + ppi = pi3(np,k) t = grids(np,k+2) q = max(1.e-8,grids(np,2+k+levs)) u = gridu(np,k) @@ -1153,7 +1153,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (mod(k,2)>0) then data2((kk-1)*6+7) = p1(np,k) data2((kk-1)*6+8) = t - data2((kk-1)*6+9) = u + data2((kk-1)*6+9) = u data2((kk-1)*6+10) = v data2((kk-1)*6+11) = q data2((kk-1)*6+12) = omega(np,k)*100. @@ -1164,16 +1164,16 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! process surface flux file fields ! !! data(8+nflx) = psfc * 100. ! SURFACE PRESSURE (PA) -!! data(7+nflx) = pmsl(np) +!! data(7+nflx) = pmsl(np) data2(8+nflx2) = psfc * 100. ! SURFACE PRESSURE (PA) - data2(7+nflx2) = pmsl(np) + data2(7+nflx2) = pmsl(np) !! dtemp = .0065 * (grids(np,1) - elevstn(np)) !! dtemp = .0100 * (grids(np,1) - elevstn(np)) !! sfc(37,np) = data(6+nflx) * .01 !! sfc(37,np) = data(7+nflx) * .01 -!! sfc(39,np) = zp2(2) !500 hPa height +!! sfc(39,np) = zp2(2) !500 hPa height sfc(37,np) = data2(7+nflx2) * .01 - sfc(39,np) = zp2(2) !500 hPa height + sfc(39,np) = zp2(2) !500 hPa height ! ! do height correction if there is no snow or if the temp is less than 0 ! G.P.LOU: @@ -1251,8 +1251,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if(sfc(12,np).gt.0.) then !check for precip then calc precip type do k = 1, leveta+1 - pp = p1(np,k) - ppi = pi3(np,k) + pp = p1(np,k) + ppi = pi3(np,k) t = grids(np,k+2) q = max(0.,grids(np,2+k+levs)) u = gridu(np,k) From db151742c13856c5e5b15dc4f90e8ff8ef36c70c Mon Sep 17 00:00:00 2001 From: BoCui-NOAA <53531984+BoCui-NOAA@users.noreply.github.com> Date: Fri, 5 May 2023 15:29:30 -0400 Subject: [PATCH 5/9] Update meteorg.f --- src/gfs_bufr.fd/meteorg.f | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f index dd08e815..32959fee 100644 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -12,9 +12,9 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! ABSTRACT: Creates BUFR meteogram files for the AVN and MRF. ! ! PROGRAM HISTORY LOG: -! 1999-07-21 HUALU PAN -! 2007-02-02 FANGLIN YANG EXPAND FOR HYBRID COORDINATES USING SIGIO -! 2009-07-24 FANGLIN YANG CHANGE OUTPUT PRESSURE TO INTEGER-LAYER +! 1999-07-21 HUALU PAN +! 2007-02-02 FANGLIN YANG EXPAND FOR HYBRID COORDINATES USING SIGIO +! 2009-07-24 FANGLIN YANG CHANGE OUTPUT PRESSURE TO INTEGER-LAYER ! PRESSURE (line 290) ! CORRECT THE TEMPERATURE ADJUSTMENT (line 238) ! 2014-03-27 DANA CARLIS UNIFY CODE WITH GFS FORECAST MODEL PRECIP @@ -22,9 +22,9 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! 2016-09-27 HUIYA CHUANG MODIFY TO READ GFS NEMS OUTPUT ON GRID SPACE ! 2017-02-27 GUANG PING LOU CHANGE OUTPUT PRECIPITATION TO HOURLY AMOUNT ! TO 120 HOURS AND 3 HOURLY TO 180 HOURS. -! 2018-02-01 GUANG PING LOU INGEST FV3GFS NEMSIO ACCUMULATED PRECIPITATION +! 2018-02-01 GUANG PING LOU INGEST FV3GFS NEMSIO ACCUMULATED PRECIPITATION ! AND RECALCULATE HOURLY AND 3 HOURLY OUTPUT DEPENDING -! ON LOGICAL VALUE OF precip_accu. +! ON LOGICAL VALUE OF precip_accu. ! 2018-02-08 GUANG PING LOU ADDED READING IN AND USING DZDT AS VERTICAL VELOCITY ! 2018-02-16 GUANG PING LOU ADDED READING IN AND USING MODEL DELP AND DELZ ! 2018-02-21 GUANG PING LOU THIS VERSION IS BACKWARD COMPATIBLE TO GFS MODEL From 334b33cd2d51f86488eca89a5bea34effaea058d Mon Sep 17 00:00:00 2001 From: BoCui-NOAA <53531984+BoCui-NOAA@users.noreply.github.com> Date: Fri, 5 May 2023 15:50:44 -0400 Subject: [PATCH 6/9] Update meteorg.f --- src/gfs_bufr.fd/meteorg.f | 77 +++++++++++++++++++-------------------- 1 file changed, 38 insertions(+), 39 deletions(-) diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f index 32959fee..dc69c835 100644 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -29,7 +29,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! 2018-02-16 GUANG PING LOU ADDED READING IN AND USING MODEL DELP AND DELZ ! 2018-02-21 GUANG PING LOU THIS VERSION IS BACKWARD COMPATIBLE TO GFS MODEL ! 2018-03-27 GUANG PING LOU CHANGE STATION ELEVATION CORRECTION LAPSE RATE FROM 0.01 TO 0.0065 -! 2018-03-28 GUANG PING LOU GENERALIZE TIME INTERVAL +! 2018-03-28 GUANG PING LOU GENERALIZE TIME INTERVAL ! 2019-07-08 GUANG PING LOU ADDED STATION CHARACTER IDS ! 2019-10-08 GUANG PING LOU MODIFY TO READ IN NetCDF FILES. RETAIN NEMSIO ! RELATED CALLS AND CLEAN UP THE CODE. @@ -44,29 +44,29 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! rlon(npoint) - longtitude ! istat(npoint) - station id ! elevstn(npoint) - station elevation (m) -! nf - forecast cycle -! fnsig - sigma file name -! idate(4) - date +! nf - forecast cycle +! fnsig - sigma file name +! idate(4) - date ! levs - input vertical layers -! kdim - sfc file dimension +! kdim - sfc file dimension ! -! OUTPUT: -! nfile - output data file channel -! jdate - date YYYYMMDDHH +! OUTPUT: +! nfile - output data file channel +! jdate - date YYYYMMDDHH ! ! ATTRIBUTES: -! LANGUAGE: +! LANGUAGE: ! MACHINE: IBM SP ! !$$$ use netcdf use nemsio_module - use sigio_module + use sigio_module use physcons use mersenne_twister ! use funcphys, only : gfuncphys use funcphys - implicit none + implicit none include 'mpif.h' type(nemsio_gfile) :: gfile type(nemsio_gfile) :: ffile @@ -190,7 +190,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, error=nf90_inq_varid(ncid, "lat", id_var) error=nf90_get_var(ncid, id_var, gdlat) !!end read NetCDF hearder info, read nemsio below if necessary - else + else call nemsio_open(gfile,trim(fnsig),'read',iret=iret) call nemsio_getfilehead(gfile,iret=iret @@ -220,14 +220,13 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, gdlon(i,j)=dum1d2((j-1)*im+i) end do end do - endif !end read in nemsio hearder if(debugprint) then do k=1,levs+1 print*,'vcoord(k,1)= ', k, vcoord(k,1) end do - do k=1,levs+1 + do k=1,levs+1 print*,'vcoord(k,2)= ', k, vcoord(k,2) end do print*,'sample lat= ',gdlat(im/5,jm/4) @@ -242,7 +241,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,1,VarName,hgt,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'surface hgt not found' - else + else VarName='hgt' LayName='sfc' call read_nemsio(gfile,im,jm,1,VarName,LayName,hgt, @@ -260,7 +259,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'surface pressure not found' - else + else VarName='pres' LayName='sfc' call read_nemsio(gfile,im,jm,1,VarName, @@ -277,7 +276,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,t3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'temp not found' - else + else VarName='tmp' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,t3d,error) @@ -296,7 +295,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,q3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'spfh not found' - else + else VarName='spfh' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,q3d,error) @@ -315,7 +314,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,uh,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'ugrd not found' - else + else VarName='ugrd' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,uh,error) @@ -334,7 +333,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,vh,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'vgrd not found' - else + else VarName='vgrd' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,vh,error) @@ -353,7 +352,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,omega3d,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'dzdt not found' - else + else VarName='dzdt' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName, @@ -373,7 +372,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'dpres not found' - else + else VarName='dpres' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName, @@ -435,7 +434,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'delz not found' - else + else VarName='delz' LayName='mid layer' call read_nemsio(gfile,im,jm,levs,VarName,LayName,delpz,error) @@ -523,7 +522,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !open T-nint below error=nf90_open(trim(fngrib2),nf90_nowrite,ncid2) if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) - else + else call nemsio_open(ffile,trim(fngrib),'read',iret=error) call nemsio_open(ffile2,trim(fngrib2),'read',iret=error) if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) @@ -535,7 +534,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid,im,jm,1,VarName,lwmask,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'lwmask not found' - else + else VarName='land' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName,lwmask,error) @@ -553,7 +552,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tmpsfc not found' - else + else VarName='tmp' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -571,7 +570,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tmp2m not found' - else + else VarName='tmp' LayName='2 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -590,7 +589,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'spfh2m not found' - else + else VarName='spfh' LayName='2 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -609,7 +608,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'ugrd10m not found' - else + else VarName='ugrd' LayName='10 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -625,7 +624,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'vgrd10m not found' - else + else VarName='vgrd' LayName='10 m above gnd' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -641,7 +640,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'soilt1 not found' - else + else VarName='tmp' LayName='0-10 cm down' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -660,7 +659,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'snod not found' - else + else VarName='snod' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -677,7 +676,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'lhtfl not found' - else + else VarName='lhtfl' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -708,7 +707,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'prate_ave not found' - else + else VarName='prate_ave' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -741,7 +740,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'cprat_ave not found' - else + else VarName='cprat_ave' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -766,7 +765,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'weasd not found' - else + else VarName='weasd' LayName='sfc' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -782,7 +781,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,12),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avelcl not found' - else + else VarName='tcdc_ave' LayName='low cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -798,7 +797,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,13),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avemcl not found' - else + else VarName='tcdc_ave' LayName='mid cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, @@ -814,7 +813,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & dum2d(:,:,14),Zreverse, & iope,ionproc,iocomms,error) if (error /= 0) print*,'tcdc_avehcl not found' - else + else VarName='tcdc_ave' LayName='high cld lay' call read_nemsio(ffile,im,jm,1,VarName,LayName, From 748e1546f390ba25e7c2d814e1a02a8f12b2ca85 Mon Sep 17 00:00:00 2001 From: BoCui-NOAA <53531984+BoCui-NOAA@users.noreply.github.com> Date: Fri, 5 May 2023 15:56:38 -0400 Subject: [PATCH 7/9] Update meteorg.f clean out extra whitespace --- src/gfs_bufr.fd/meteorg.f | 44 +++++++++++++++++++-------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f index dc69c835..9c289ca0 100644 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -391,7 +391,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, do j=1,jm do i=1,im pint(i,j,k)=vcoord(k,1) - + +vcoord(k,2)*pint(i,j,1) + + +vcoord(k,2)*pint(i,j,1) end do end do end do @@ -753,7 +753,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, do j=1,jm do i=1,im dum2d(i,j,10)=(apcp(i,j)*fhour-cpcp(i,j)*ap)*3600.0 - & + & end do end do @@ -845,7 +845,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, jdum=jjdum(np) else -! find nearest neighbor +! find nearest neighbor rdum=rlon(np) if(rdum<0.)rdum=rdum+360. @@ -863,7 +863,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, exit else if(landwater(np) == lwmask(i+1,j))then idum=i+1 - jdum=j ! 2 + jdum=j ! 2 exit else if(landwater(np) == lwmask(i-1,j))then idum=i-1 @@ -1002,7 +1002,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, grids(np,1)=hgt(idum,jdum) grids(np,2)=pint(idum,jdum,1) - + sfc(5,np)=dum2d(idum,jdum,1) sfc(6,np)=dum2d(idum,jdum,6) sfc(17,np)=dum2d(idum,jdum,8) @@ -1019,10 +1019,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, CC There may be cases where convective precip is greater than total precip CC due to rounding and interpolation errors, correct it here -G.P. Lou: - if(sfc(11,np) .gt. sfc(12,np)) sfc(11,np)=sfc(12,np) + if(sfc(11,np) .gt. sfc(12,np)) sfc(11,np)=sfc(12,np) do k=1,levs - grids(np,k+2)=t3d(idum,jdum,k) + grids(np,k+2)=t3d(idum,jdum,k) grids(np,k+2+levs)=q3d(idum,jdum,k) grids(np,k+2+2*levs)=omega3d(idum,jdum,k) gridu(np,k)=uh(idum,jdum,k) @@ -1033,7 +1033,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !! z1(np,k)=0.5*(zint(idum,jdum,k)+zint(idum,jdum,k+1)) end do - end do + end do print*,'finish finding nearest neighbor for each station' @@ -1049,7 +1049,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !! if(recn_dzdt == 0 ) then !!DZDT do k = 1, levs do np = 1, npoint - omega(np,k) = grids(np,2+levs*2+k) + omega(np,k) = grids(np,2+levs*2+k) enddo enddo if(debugprint) @@ -1093,7 +1093,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, & pmsl(np),zp(np,1:levs),zp2(1:2)) enddo print *, 'call gslp pmsl= ', (pmsl(np),np=1,20) - if(recn_delz == -9999) then + if(recn_delz == -9999) then print*, 'using calculated height ' else print*, 'using model height m' @@ -1138,7 +1138,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! look for the layer above 500 mb for precip type computation ! if(pi3(np,k).ge.50000.) leveta = k - ppi = pi3(np,k) + ppi = pi3(np,k) t = grids(np,k+2) q = max(1.e-8,grids(np,2+k+levs)) u = gridu(np,k) @@ -1152,8 +1152,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (mod(k,2)>0) then data2((kk-1)*6+7) = p1(np,k) data2((kk-1)*6+8) = t - data2((kk-1)*6+9) = u - data2((kk-1)*6+10) = v + data2((kk-1)*6+9) = u + data2((kk-1)*6+10) = v data2((kk-1)*6+11) = q data2((kk-1)*6+12) = omega(np,k)*100. endif @@ -1163,16 +1163,16 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! process surface flux file fields ! !! data(8+nflx) = psfc * 100. ! SURFACE PRESSURE (PA) -!! data(7+nflx) = pmsl(np) +!! data(7+nflx) = pmsl(np) data2(8+nflx2) = psfc * 100. ! SURFACE PRESSURE (PA) - data2(7+nflx2) = pmsl(np) + data2(7+nflx2) = pmsl(np) !! dtemp = .0065 * (grids(np,1) - elevstn(np)) !! dtemp = .0100 * (grids(np,1) - elevstn(np)) !! sfc(37,np) = data(6+nflx) * .01 !! sfc(37,np) = data(7+nflx) * .01 -!! sfc(39,np) = zp2(2) !500 hPa height +!! sfc(39,np) = zp2(2) !500 hPa height sfc(37,np) = data2(7+nflx2) * .01 - sfc(39,np) = zp2(2) !500 hPa height + sfc(39,np) = zp2(2) !500 hPa height ! ! do height correction if there is no snow or if the temp is less than 0 ! G.P.LOU: @@ -1194,10 +1194,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! !G.P. Lou 20200501: !convert instantaneous surface latent heat net flux to surface -!evapolation 1 W m-2 = 0.0864 MJ m-2 day-1 +!evapolation 1 W m-2 = 0.0864 MJ m-2 day-1 ! and 1 mm day-1 = 2.45 MJ m-2 day-1 ! equivament to 0.0864/2.54 = 0.035265 -! equivament to 2.54/0.0864 = 28.3565 +! equivament to 2.54/0.0864 = 28.3565 if(debugprint) + print*,'evaporation (stn 000692)= ',sfc(17,np) !! data(9+nflx) = sfc(5,np) ! tsfc (K) @@ -1250,8 +1250,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if(sfc(12,np).gt.0.) then !check for precip then calc precip type do k = 1, leveta+1 - pp = p1(np,k) - ppi = pi3(np,k) + pp = p1(np,k) + ppi = pi3(np,k) t = grids(np,k+2) q = max(0.,grids(np,2+k+levs)) u = gridu(np,k) @@ -1269,7 +1269,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, gq0(1,k) = q prsl(1,k) = pp prsi(1,k)=ppi - phii(1,k)=zp(np,k) !height in meters + phii(1,k)=zp(np,k) !height in meters enddo ! Use GFS routine calpreciptype.f to calculate precip type xlat=rlat(np) From a52049e4640d76b4f4e02b5f59479a1a4123fdaf Mon Sep 17 00:00:00 2001 From: BoCui-NOAA <53531984+BoCui-NOAA@users.noreply.github.com> Date: Fri, 5 May 2023 15:58:14 -0400 Subject: [PATCH 8/9] Update meteorg.f clean out extra whitespace --- src/gfs_bufr.fd/meteorg.f | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f index 9c289ca0..35577032 100644 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -1031,10 +1031,10 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, z1(np,k)=zint(idum,jdum,k+1) !! p1(np,k)=0.5*(pint(idum,jdum,k)+pint(idum,jdum,k+1)) !! z1(np,k)=0.5*(zint(idum,jdum,k)+zint(idum,jdum,k+1)) - + end do end do - + print*,'finish finding nearest neighbor for each station' do np = 1, npoint From b98538d994a9ab4b9203b89a6a9f5a0a5ce98e73 Mon Sep 17 00:00:00 2001 From: BoCui-NOAA <53531984+BoCui-NOAA@users.noreply.github.com> Date: Fri, 5 May 2023 16:21:33 -0400 Subject: [PATCH 9/9] Update meteorg.f clean out extra whitespace --- src/gfs_bufr.fd/meteorg.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f index 35577032..337a3c7f 100644 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -1105,7 +1105,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, endif print*,'finish computing MSLP' print*,'finish computing zp ', (zp(11,k),k=1,levs) - print*,'finish computing zp2(1-2) ', zp2(1),zp2(2) + print*,'finish computing zp2(1-2) ', zp2(1),zp2(2) ! ! prepare buffer data !