From 71468666da0452fc9b13a19f4e093c5e1930c13c Mon Sep 17 00:00:00 2001 From: "Bo.Cui" Date: Mon, 19 Aug 2024 14:30:50 -0400 Subject: [PATCH 1/6] Update bufr sounding codes to handle one forecast at a time --- src/gfs_bufr.fd/CMakeLists.txt | 2 +- src/gfs_bufr.fd/CMakeLists.txt_org | 53 ++ src/gfs_bufr.fd/buff.f | 4 +- src/gfs_bufr.fd/buff.f_org | 91 ++ src/gfs_bufr.fd/gfsbufr.f | 152 ++-- src/gfs_bufr.fd/gfsbufr.f_org | 275 ++++++ src/gfs_bufr.fd/meteorg.f | 389 ++------ src/gfs_bufr.fd/meteorg.f_org | 1326 ++++++++++++++++++++++++++++ 8 files changed, 1906 insertions(+), 386 deletions(-) create mode 100644 src/gfs_bufr.fd/CMakeLists.txt_org create mode 100644 src/gfs_bufr.fd/buff.f_org create mode 100644 src/gfs_bufr.fd/gfsbufr.f_org create mode 100644 src/gfs_bufr.fd/meteorg.f_org diff --git a/src/gfs_bufr.fd/CMakeLists.txt b/src/gfs_bufr.fd/CMakeLists.txt index 7df490d0..09c13e13 100644 --- a/src/gfs_bufr.fd/CMakeLists.txt +++ b/src/gfs_bufr.fd/CMakeLists.txt @@ -10,7 +10,7 @@ list(APPEND fortran_src mstadb.f newsig1.f read_nemsio.f - #read_netcdf.f + read_netcdf.f read_netcdf_p.f rsearch.f svp.f diff --git a/src/gfs_bufr.fd/CMakeLists.txt_org b/src/gfs_bufr.fd/CMakeLists.txt_org new file mode 100644 index 00000000..7df490d0 --- /dev/null +++ b/src/gfs_bufr.fd/CMakeLists.txt_org @@ -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/buff.f b/src/gfs_bufr.fd/buff.f index 82acf036..b5970b2d 100644 --- a/src/gfs_bufr.fd/buff.f +++ b/src/gfs_bufr.fd/buff.f @@ -28,7 +28,9 @@ subroutine buff(nint1,nend1,nint3,nend3,npoint,idate,jdate,levs, do np = 1, npoint C OPEN BUFR OUTPUT FILE. write(fnbufr,fmto) dird(1:lss),istat(np),jdate - print *, ' fnbufr =', fnbufr + if(np==1.or.np==100) then + print *, ' fnbufr =', fnbufr + endif open(unit=19,file=fnbufr,form='unformatted', & status='new', iostat=ios) IF ( ios .ne. 0 ) THEN diff --git a/src/gfs_bufr.fd/buff.f_org b/src/gfs_bufr.fd/buff.f_org new file mode 100644 index 00000000..82acf036 --- /dev/null +++ b/src/gfs_bufr.fd/buff.f_org @@ -0,0 +1,91 @@ + subroutine buff(nint1,nend1,nint3,nend3,npoint,idate,jdate,levs, + & dird,lss,istat,sbset,seqflg,clist,npp,wrkd) + character*512 dird, fnbufr, fmto +!! integer nint, nend, npoint, idate(4), levs, jdate + integer nint1, nend1, nint3, nend3 + integer npoint, idate(4), levs, jdate + real*8 data(6*levs+25), wrkd(1) + integer idtln, nf, nfile, np + integer lss, istat(npoint), ios + CHARACTER*8 SBSET + LOGICAL SEQFLG(4) + CHARACTER*80 CLIST(4) + INTEGER NPP(4) + CHARACTER*8 SEQNAM(4) + FMTO = '(A,".",I6.6,".",I10)' + idtln = 8 + nfile = 20 +C print *, 'inside buff.f nint1,nend1,nint3,nend3,jdate=' +C print *, nint1,nend1,nint3,nend3,jdate + do nf = 0, nend1, nint1 + nfile = nfile + 1 + rewind nfile + enddo + do nf = nend1+nint3, nend3, nint3 + nfile = nfile + 1 + rewind nfile + enddo + do np = 1, npoint +C OPEN BUFR OUTPUT FILE. + write(fnbufr,fmto) dird(1:lss),istat(np),jdate + print *, ' fnbufr =', fnbufr + open(unit=19,file=fnbufr,form='unformatted', + & status='new', iostat=ios) + IF ( ios .ne. 0 ) THEN + WRITE (6,*) ' CANNOT open ', 19 + STOP + END IF + CALL OPENBF ( 19, 'OUT', 1 ) + nfile = 20 + do nf = 0, nend1, nint1 + nfile = nfile + 1 + read(nfile) data + if(np.eq.1) then + print *, ' creating bufr file for np, nfile =', + & np, nfile + endif +CC WRITE DATA MESSAGE TO BUFR OUTPUT FILE. +CC LUNTBL=-9 BECAUSE BUFR TABLE FILE NOT USED HERE. +CC SEQNAM=XXXXXX BECAUSE MNEMONIC SEQUENCE NAMES NOT USED HERE. + CALL BFRIZE ( -9, 19, SBSET, + & idate(4), iDATE(2), + & iDATE(3), iDATE(1), + & 'XXXXXX', SEQFLG, 4, .FALSE., DATA, levs, + & CLIST, NPP, WRKD, IRET ) + IF ( IRET .NE. 0 ) THEN + PRINT *,' BFRIZE FAILED ' + ENDIF +c 300 continue + enddo +C 3hourly output starts here +!! print *, 'buff.f nfile,nend1+nint3,nend3,nint3=' +!! print *, nfile,nend1+nint3,nend3,nint3 + do nf = nend1+nint3, nend3, nint3 + nfile = nfile + 1 + read(nfile) data + if(np.eq.1) then + print *, ' creating bufr file for np, nfile =', + & np, nfile + endif +C print *, 'read2 in fort(nfile) =', nfile +CC WRITE DATA MESSAGE TO BUFR OUTPUT FILE. +CC LUNTBL=-9 BECAUSE BUFR TABLE FILE NOT USED HERE. +CC SEQNAM=XXXXXX BECAUSE MNEMONIC SEQUENCE NAMES NOT USED HERE. + CALL BFRIZE ( -9, 19, SBSET, + & idate(4), iDATE(2), + & iDATE(3), iDATE(1), + & 'XXXXXX', SEQFLG, 4, .FALSE., DATA, levs, + & CLIST, NPP, WRKD, IRET ) + IF ( IRET .NE. 0 ) THEN + PRINT *,' BFRIZE FAILED ' + ENDIF + enddo + CALL BFRIZE ( 0, 19, SBSET, + & IDATE(4), IDATE(2), + & IDATE(3), IDATE(1), + & 'XXXXXX', SEQFLG, 4, .FALSE., DATA, levs, + & CLIST, NPP, WRKD, IRET ) + call closbf(19) + enddo + return + end diff --git a/src/gfs_bufr.fd/gfsbufr.f b/src/gfs_bufr.fd/gfsbufr.f index c1bb1d35..a3f1acbf 100644 --- a/src/gfs_bufr.fd/gfsbufr.f +++ b/src/gfs_bufr.fd/gfsbufr.f @@ -14,6 +14,8 @@ program meteormrf C 17-02-27 GUANG PING LOU: CHANGE MODEL OUTPUT READ-IN TO HOURLY C TO 120 HOURS AND 3 HOURLY TO 180 HOURS. C 19-07-16 GUANG PING LOU: CHANGE FROM NEMSIO TO GRIB2. +C 24-08-08 Bo Cui: UPDATE TO HANDLE ONE FORECAST AT A TIME +C REMOVE NEMSIO INPUT FILES C C C USAGE: @@ -42,25 +44,21 @@ program meteormrf C C$$$ use netcdf - use mpi - use nemsio_module use sigio_module implicit none -!! include 'mpif.h' integer,parameter:: nsta=3000 integer,parameter:: ifile=11 integer,parameter:: levso=64 integer(sigio_intkind):: irets - type(nemsio_gfile) :: gfile integer ncfsig, nsig - integer istat(nsta), idate(4), jdate + integer istat(nsta), idate(4), jdate, nfhour integer :: levs,nstart,nend,nint,nsfc,levsi,im,jm integer :: npoint,np,ist,is,iret,lss,nss,nf,nsk,nfile integer :: ielev integer :: lsfc real :: alat,alon,rla,rlo real :: wrkd(1),dummy - real rlat(nsta), rlon(nsta), elevstn(nsta) + real rlat(nsta), rlon(nsta), elevstn(nsta), fhour integer iidum(nsta),jjdum(nsta) integer nint1, nend1, nint3, nend3, np1 integer landwater(nsta) @@ -68,7 +66,7 @@ program meteormrf character*4 t3 character*4 cstat(nsta) character*32 desc - character*512 dird, fnsig + character*512 dird, fnsig,fngrib,fngrib2 logical f00, makebufr CHARACTER*8 SBSET LOGICAL SEQFLG(4) @@ -80,6 +78,7 @@ program meteormrf integer :: error, ncid, id_var,dimid character(len=10) :: dim_nam character(len=6) :: fformat + character(len=100) :: long_name !added from Cory integer :: iope, ionproc integer, allocatable :: iocomms(:) @@ -94,14 +93,10 @@ program meteormrf C namelist /nammet/ levs, makebufr, dird, & nstart, nend, nint, nend1, nint1, - & nint3, nsfc, f00, fformat, np1 + & nint3, nsfc, f00, fformat, np1, + & fnsig,fngrib,fngrib2 - call mpi_init(ierr) - call mpi_comm_rank(MPI_COMM_WORLD,mrank,ierr) - call mpi_comm_size(MPI_COMM_WORLD,msize,ierr) - if(mrank.eq.0) then - CALL W3TAGB('METEOMRF',1999,0202,0087,'NP23') - endif + CALL W3TAGB('METEOMRF',1999,0202,0087,'NP23') open(5,file='gfsparm') read(5,nammet) write(6,nammet) @@ -150,7 +145,7 @@ program meteormrf enddo endif 98 FORMAT (3I6, 2F9.2) - if (mrank.eq.0.and.makebufr) then + if (makebufr) then REWIND 1 READ (1,100) SBSET 100 FORMAT ( ////// 2X, A8 ) @@ -171,40 +166,23 @@ program meteormrf lss = lss - 1 END DO C - endif - nsig = 11 - nss = nstart + nint - if(f00) nss = nstart -c do nf = nss, nend, nint - ntot = (nend - nss) / nint + 1 - ntask = mrank/(float(msize)/float(ntot)) - nf = ntask * nint + nss - print*,'n0 ntot nint nss mrank msize' - print*, n0,ntot,nint,nss,mrank,msize - print*,'nf, ntask= ', nf, ntask + else ! else of makebufr + +!! nfile - output data file channel, start from fort.21 +!! nf - forecast hour + + nf=nstart if(nf .le. nend1) then nfile = 21 + (nf / nint1) else nfile = 21 + (nend1/nint1) + (nf-nend1)/nint3 endif print*, 'nf,nint,nfile = ',nf,nint,nfile - if(nf.le.nend) then - if(nf.lt.10) then - fnsig = 'sigf0' - write(fnsig(6:6),'(i1)') nf - ncfsig = 6 - elseif(nf.lt.100) then - fnsig = 'sigf' - write(fnsig(5:6),'(i2)') nf - ncfsig = 6 - else - fnsig = 'sigf' - write(fnsig(5:7),'(i3)') nf - ncfsig = 7 - endif - print *, 'Opening file : ',fnsig + print *, 'Opening atmos file : ',trim(fnsig) + print *, 'Opening surface file : ',trim(fngrib) + print *, 'Opening surface file 2 : ',trim(fngrib2) -!! read in either nemsio or NetCDF files +!! read in NetCDF files if (fformat == 'netcdf') then error=nf90_open(trim(fnsig),nf90_nowrite,ncid) error=nf90_inq_dimid(ncid,"grid_xt",dimid) @@ -214,62 +192,52 @@ program meteormrf error=nf90_inq_dimid(ncid,"pfull",dimid) error=nf90_inquire_dimension(ncid,dimid,dim_nam,levsi) error=nf90_close(ncid) - print*,'NetCDF file im,jm,lm= ',im,jm,levs,levsi - - else - call nemsio_init(iret=irets) - print *,'nemsio_init, iret=',irets - call nemsio_open(gfile,trim(fnsig),'read',iret=irets) - if ( irets /= 0 ) then - print*,"fail to open nems atmos file";stop - endif - - call nemsio_getfilehead(gfile,iret=irets - & ,dimx=im,dimy=jm,dimz=levsi) - if( irets /= 0 ) then - print*,'error finding model dimensions '; stop - endif - print*,'nemsio file im,jm,lm= ',im,jm,levsi - call nemsio_close(gfile,iret=irets) - endif - allocate (iocomms(0:ntot)) - if (fformat == 'netcdf') then - print*,'iocomms= ', iocomms - call mpi_comm_split(MPI_COMM_WORLD,ntask,0,iocomms(ntask),ierr) - call mpi_comm_rank(iocomms(ntask), iope, ierr) - call mpi_comm_size(iocomms(ntask), ionproc, ierr) +! print*,'NetCDF file im,jm,lm= ',im,jm,levs,levsi call meteorg(npoint,rlat,rlon,istat,cstat,elevstn, - & nf,nfile,fnsig,jdate,idate, + & nf,nfile,fnsig,fngrib,fngrib2,jdate,idate, & levsi,im,jm,nsfc, & landwater,nend1, nint1, nint3, iidum,jjdum,np1, - & fformat,iocomms(ntask),iope,ionproc) - call mpi_barrier(iocomms(ntask), ierr) - call mpi_comm_free(iocomms(ntask), ierr) - else -!! For nemsio input - call meteorg(npoint,rlat,rlon,istat,cstat,elevstn, - & nf,nfile,fnsig,jdate,idate, - & levs,im,jm,nsfc, - & landwater,nend1, nint1, nint3, iidum,jjdum,np1, - & fformat,iocomms(ntask),iope,ionproc) - endif - endif - call mpi_barrier(mpi_comm_world,ierr) - call mpi_finalize(ierr) - if(mrank.eq.0) then - print *, ' starting to make bufr files' - print *, ' makebufr= ', makebufr - print *, 'nint1,nend1,nint3,nend= ',nint1,nend1,nint3,nend -!! idate = 0 7 1 2019 -!! jdate = 2019070100 + & fformat) + endif ! end of process + + endif ! endif of makebufr if(makebufr) then - nend3 = nend - call buff(nint1,nend1,nint3,nend3, + +! read in NetCDF file header info +! sample of idate and jdate +! idate = 0 7 1 2019 +! jdate = 2019070100 + + if (fformat == 'netcdf') then + error=nf90_open(trim(fnsig),nf90_nowrite,ncid) + error=nf90_inq_varid(ncid, "time", id_var) + error=nf90_get_var(ncid, id_var, nfhour) + error=nf90_get_att(ncid,id_var,"units",long_name) + error=nf90_close(ncid) + endif + + read(long_name(13:16),"(i4)")idate(4) + read(long_name(18:19),"(i2)")idate(2) + read(long_name(21:22),"(i2)")idate(3) + read(long_name(24:25),"(i2)")idate(1) + fhour=float(nfhour) + jdate = idate(4)*1000000 + idate(2)*10000+ + & idate(3)*100 + idate(1) + + print *, ' starting to make bufr files' + print *, ' makebufr= ', makebufr + print *, ' processing forecast hour ', fhour + print *, 'nint1,nend1,nint3,nend= ',nint1,nend1,nint3,nend + print *, 'idate,jdate=',idate,jdate + + nend3 = nend + call buff(nint1,nend1,nint3,nend3, & npoint,idate,jdate,levso, & dird,lss,istat,sbset,seqflg,clist,npp,wrkd) - CALL W3TAGE('METEOMRF') - endif - endif + CALL W3TAGE('METEOMRF') + + endif ! end of makebufr + end diff --git a/src/gfs_bufr.fd/gfsbufr.f_org b/src/gfs_bufr.fd/gfsbufr.f_org new file mode 100644 index 00000000..c1bb1d35 --- /dev/null +++ b/src/gfs_bufr.fd/gfsbufr.f_org @@ -0,0 +1,275 @@ + program meteormrf +C$$$ MAIN PROGRAM DOCUMENTATION BLOCK +C +C MAIN PROGRAM: METEOMRF +C PRGMMR: PAN ORG: NP23 DATE: 1999-07-21 +C +C ABSTRACT: Creates BUFR meteogram files for the AVN and MRF. +C +C PROGRAM HISTORY LOG: +C 99-07-21 Hualu Pan +C 16-09-27 HUIYA CHUANG MODIFY TO READ GFS NEMS OUTPUT ON GRID SPACE +C 16-10-15 HUIYA CHUANG: CONSOLIDATE TO READ FLUX FIELDS IN THIS +C PACKAGE TOO AND THIS SPEEDS UP BFS BUFR BY 3X +C 17-02-27 GUANG PING LOU: CHANGE MODEL OUTPUT READ-IN TO HOURLY +C TO 120 HOURS AND 3 HOURLY TO 180 HOURS. +C 19-07-16 GUANG PING LOU: CHANGE FROM NEMSIO TO GRIB2. +C +C +C USAGE: +C INPUT FILES: +C FTxxF001 - UNITS 11 THRU 49 +C PARM - UNIT 5 (STANDARD READ) +C +C OUTPUT FILES: (INCLUDING SCRATCH FILES) +C FTxxF001 - UNITS 51 THRU 79 +C FTxxF001 - UNIT 6 (STANDARD PRINTFILE) +C +C SUBPROGRAMS CALLED: (LIST ALL CALLED FROM ANYWHERE IN CODES) +C UNIQUE: - ROUTINES THAT ACCOMPANY SOURCE FOR COMPILE +C LIBRARY: +C W3LIB - +C +C EXIT STATES: +C COND = 0 - SUCCESSFUL RUN +C =NNNN - TROUBLE OR SPECIAL FLAG - SPECIFY NATURE +C +C REMARKS: LIST CAVEATS, OTHER HELPFUL HINTS OR INFORMATION +C +C ATTRIBUTES: +C LANGUAGE: INDICATE EXTENSIONS, COMPILER OPTIONS +C MACHINE: IBM SP +C +C$$$ + use netcdf + use mpi + use nemsio_module + use sigio_module + implicit none +!! include 'mpif.h' + integer,parameter:: nsta=3000 + integer,parameter:: ifile=11 + integer,parameter:: levso=64 + integer(sigio_intkind):: irets + type(nemsio_gfile) :: gfile + integer ncfsig, nsig + integer istat(nsta), idate(4), jdate + integer :: levs,nstart,nend,nint,nsfc,levsi,im,jm + integer :: npoint,np,ist,is,iret,lss,nss,nf,nsk,nfile + integer :: ielev + integer :: lsfc + real :: alat,alon,rla,rlo + real :: wrkd(1),dummy + real rlat(nsta), rlon(nsta), elevstn(nsta) + integer iidum(nsta),jjdum(nsta) + integer nint1, nend1, nint3, nend3, np1 + integer landwater(nsta) + character*1 ns, ew + character*4 t3 + character*4 cstat(nsta) + character*32 desc + character*512 dird, fnsig + logical f00, makebufr + CHARACTER*8 SBSET + LOGICAL SEQFLG(4) + CHARACTER*80 CLIST(4) + INTEGER NPP(4) + CHARACTER*8 SEQNAM(4) + integer ierr, mrank, msize,ntask + integer n0, ntot + integer :: error, ncid, id_var,dimid + character(len=10) :: dim_nam + character(len=6) :: fformat + !added from Cory + integer :: iope, ionproc + integer, allocatable :: iocomms(:) +C + DATA SBSET / 'ABCD1234' / +C + DATA SEQFLG / .FALSE., .TRUE., .FALSE., .FALSE. / +C + DATA SEQNAM / 'HEADR', 'PROFILE', 'CLS1' ,'D10M' / +c DATA SEQNAM / 'HEADR', 'PRES TMDB UWND VWND SPFH OMEG', +c & 'CLS1' ,'D10M' / +C + namelist /nammet/ levs, makebufr, dird, + & nstart, nend, nint, nend1, nint1, + & nint3, nsfc, f00, fformat, np1 + + call mpi_init(ierr) + call mpi_comm_rank(MPI_COMM_WORLD,mrank,ierr) + call mpi_comm_size(MPI_COMM_WORLD,msize,ierr) + if(mrank.eq.0) then + CALL W3TAGB('METEOMRF',1999,0202,0087,'NP23') + endif + open(5,file='gfsparm') + read(5,nammet) + write(6,nammet) + npoint = 0 + 99 FORMAT (I6, F6.2,A1, F7.2,A1,1X,A4,1X,I2, A28, I4) + do np = 1, nsta+2 + read(8,99,end=200) IST,ALAT,NS,ALON,EW,T3,lsfc,DESC,IELEV +CC print*," IST,ALAT,NS,ALON,EW,T3,lsfc,DESC,IELEV= " +CC print*, IST,ALAT,NS,ALON,EW,T3,lsfc,DESC,IELEV + if(alat.lt.95.) then + npoint = npoint + 1 + RLA = 9999. + IF (NS .EQ. 'N') RLA = ALAT + IF (NS .EQ. 'S') RLA = -ALAT + RLO = 9999. + IF (EW .EQ. 'E') RLO = ALON + IF (EW .EQ. 'W') RLO = -ALON + rlat(npoint) = rla + rlon(npoint) = rlo + istat(npoint) = ist + cstat(npoint) = T3 + elevstn(npoint) = ielev + + if(lsfc .le. 9) then + landwater(npoint) = 2 !!nearest + else if(lsfc .le. 19) then + landwater(npoint) = 1 !!land + else if(lsfc .ge. 20) then + landwater(npoint) = 0 !!water + endif + endif + enddo + 200 continue + if(npoint.le.0) then + print *, ' station list file is empty, abort program' + call abort + elseif(npoint.gt.nsta) then + print *, ' number of station exceeds nsta, abort program' + call abort + endif +! print*,'npoint= ', npoint +! print*,'np,IST,idum,jdum,rlat(np),rlon(np)= ' + if(np1 == 0) then + do np = 1, npoint + read(7,98) IST, iidum(np), jjdum(np), ALAT, ALON + enddo + endif + 98 FORMAT (3I6, 2F9.2) + if (mrank.eq.0.and.makebufr) then + REWIND 1 + READ (1,100) SBSET + 100 FORMAT ( ////// 2X, A8 ) + PRINT 120, SBSET + 120 FORMAT ( ' SBSET=#', A8, '#' ) + REWIND 1 +C +C READ PARM NAMES AND NUMBER OF PARM NAMES FROM BUFR TABLE. + DO IS = 1,4 + CALL BFRHDR ( 1, SEQNAM(IS), SEQFLG(IS), + X CLIST(IS), NPP(IS), IRET ) + IF ( IRET .NE. 0 ) THEN + PRINT*, ' CALL BFRHDR IRET=', IRET + ENDIF + ENDDO + lss = len ( dird ) + DO WHILE ( dird (lss:lss) .eq. ' ' ) + lss = lss - 1 + END DO +C + endif + nsig = 11 + nss = nstart + nint + if(f00) nss = nstart +c do nf = nss, nend, nint + ntot = (nend - nss) / nint + 1 + ntask = mrank/(float(msize)/float(ntot)) + nf = ntask * nint + nss + print*,'n0 ntot nint nss mrank msize' + print*, n0,ntot,nint,nss,mrank,msize + print*,'nf, ntask= ', nf, ntask + if(nf .le. nend1) then + nfile = 21 + (nf / nint1) + else + nfile = 21 + (nend1/nint1) + (nf-nend1)/nint3 + endif + print*, 'nf,nint,nfile = ',nf,nint,nfile + if(nf.le.nend) then + if(nf.lt.10) then + fnsig = 'sigf0' + write(fnsig(6:6),'(i1)') nf + ncfsig = 6 + elseif(nf.lt.100) then + fnsig = 'sigf' + write(fnsig(5:6),'(i2)') nf + ncfsig = 6 + else + fnsig = 'sigf' + write(fnsig(5:7),'(i3)') nf + ncfsig = 7 + endif + print *, 'Opening file : ',fnsig + +!! read in either nemsio or NetCDF files + if (fformat == 'netcdf') then + error=nf90_open(trim(fnsig),nf90_nowrite,ncid) + error=nf90_inq_dimid(ncid,"grid_xt",dimid) + error=nf90_inquire_dimension(ncid,dimid,dim_nam,im) + error=nf90_inq_dimid(ncid,"grid_yt",dimid) + error=nf90_inquire_dimension(ncid,dimid,dim_nam,jm) + error=nf90_inq_dimid(ncid,"pfull",dimid) + error=nf90_inquire_dimension(ncid,dimid,dim_nam,levsi) + error=nf90_close(ncid) + print*,'NetCDF file im,jm,lm= ',im,jm,levs,levsi + + else + call nemsio_init(iret=irets) + print *,'nemsio_init, iret=',irets + call nemsio_open(gfile,trim(fnsig),'read',iret=irets) + if ( irets /= 0 ) then + print*,"fail to open nems atmos file";stop + endif + + call nemsio_getfilehead(gfile,iret=irets + & ,dimx=im,dimy=jm,dimz=levsi) + if( irets /= 0 ) then + print*,'error finding model dimensions '; stop + endif + print*,'nemsio file im,jm,lm= ',im,jm,levsi + call nemsio_close(gfile,iret=irets) + endif + allocate (iocomms(0:ntot)) + if (fformat == 'netcdf') then + print*,'iocomms= ', iocomms + call mpi_comm_split(MPI_COMM_WORLD,ntask,0,iocomms(ntask),ierr) + call mpi_comm_rank(iocomms(ntask), iope, ierr) + call mpi_comm_size(iocomms(ntask), ionproc, ierr) + + call meteorg(npoint,rlat,rlon,istat,cstat,elevstn, + & nf,nfile,fnsig,jdate,idate, + & levsi,im,jm,nsfc, + & landwater,nend1, nint1, nint3, iidum,jjdum,np1, + & fformat,iocomms(ntask),iope,ionproc) + call mpi_barrier(iocomms(ntask), ierr) + call mpi_comm_free(iocomms(ntask), ierr) + else +!! For nemsio input + call meteorg(npoint,rlat,rlon,istat,cstat,elevstn, + & nf,nfile,fnsig,jdate,idate, + & levs,im,jm,nsfc, + & landwater,nend1, nint1, nint3, iidum,jjdum,np1, + & fformat,iocomms(ntask),iope,ionproc) + endif + endif + call mpi_barrier(mpi_comm_world,ierr) + call mpi_finalize(ierr) + if(mrank.eq.0) then + print *, ' starting to make bufr files' + print *, ' makebufr= ', makebufr + print *, 'nint1,nend1,nint3,nend= ',nint1,nend1,nint3,nend +!! idate = 0 7 1 2019 +!! jdate = 2019070100 + + if(makebufr) then + nend3 = nend + call buff(nint1,nend1,nint3,nend3, + & npoint,idate,jdate,levso, + & dird,lss,istat,sbset,seqflg,clist,npp,wrkd) + CALL W3TAGE('METEOMRF') + endif + endif + end diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f index d4fa338d..852c125d 100644 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -1,8 +1,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, - & nf,nfile,fnsig,jdate,idate, + & nf,nfile,fnsig,fngrib,fngrib2,jdate,idate, & levs,im,jm,kdim, & landwater,nend1,nint1,nint3,iidum,jjdum,np1, - & fformat,iocomms,iope,ionproc) + & fformat) !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . @@ -36,6 +36,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! 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 +! 2024-08-08 Bo Cui UPDATE TO HANDLE ONE FORECAST AT A TIME, REMOVE NEMSIO INPUT FILES ! ! USAGE: CALL PROGRAM meteorg ! INPUT: @@ -46,6 +47,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! elevstn(npoint) - station elevation (m) ! nf - forecast cycle ! fnsig - sigma file name +! fngrib - surface file name +! fngrib2 - surface file name ! idate(4) - date ! levs - input vertical layers ! kdim - sfc file dimension @@ -60,7 +63,6 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! !$$$ use netcdf - use nemsio_module use sigio_module use physcons use mersenne_twister @@ -68,10 +70,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, use funcphys implicit none include 'mpif.h' - type(nemsio_gfile) :: gfile - type(nemsio_gfile) :: ffile - type(nemsio_gfile) :: ffile2 - integer :: nfile,npoint,levs,kdim + integer :: nfile,npoint,levs,kdim,nsoil + character(len=10) :: dim_nam integer :: nfile1 integer :: i,j,im,jm,kk,idum,jdum,idvc,idsl ! idsl Integer(sigio_intkind) semi-lagrangian id @@ -82,7 +82,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, integer :: idate(4),nij,nflx2,np,k,l,nf,nfhour,np1 integer :: idate_nems(7) integer :: iret,jdate,leveta,lm,lp1 - character*512 :: fnsig,fngrib + character*512 :: fnsig,fngrib,fngrib2 !! real*8 :: data(6*levs+25) real*8 :: data2(6*levso+25) real*8 :: rstat1 @@ -100,6 +100,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, real,dimension(im,jm,15) :: dum2d real,dimension(im,jm,levs) :: t3d, q3d, uh, vh,omega3d real,dimension(im,jm,levs) :: delpz +! real,dimension(im,jm,4) :: soilt3d + real, allocatable :: soilt3d(:,:,:) real,dimension(im,jm,levs+1) :: pint, zint real,dimension(npoint,levs) :: gridu,gridv,omega,qnew,zp real,dimension(npoint,levs) :: p1,pd3,ttnew @@ -117,7 +119,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, integer :: n3dfercld,iseedl integer :: istat(npoint) logical :: trace -!! logical, parameter :: debugprint=.true. +! logical, parameter :: debugprint=.true. logical, parameter :: debugprint=.false. character lprecip_accu*3 real, parameter :: ERAD=6.371E6 @@ -125,7 +127,6 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, real :: ap integer :: nf1, fint integer :: nend1, nint1, nint3 - character*150 :: fngrib2 integer recn_dpres,recn_delz,recn_dzdt integer :: jrec equivalence (cstat1,rstat1) @@ -139,7 +140,6 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, character(7) :: zone character(3) :: Zreverse character(20) :: VarName,LayName - integer iocomms,iope,ionproc nij = 12 !! nflx = 6 * levs @@ -160,9 +160,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! enddo if(fformat .eq. "netcdf") then - print*,'iocomms inside meteorg.f=', iocomms - error=nf90_open(trim(fnsig),ior(nf90_nowrite,nf90_mpiio), - & ncid,comm=iocomms, info = mpi_info_null) + error=nf90_open(trim(fnsig),nf90_nowrite,ncid) error=nf90_get_att(ncid,nf90_global,"ak",vdummy) do k = 1, levs+1 vcoord(k,1)=vdummy(levs-k+1) @@ -189,38 +187,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, error=nf90_get_var(ncid, id_var, gdlon) 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 - - call nemsio_open(gfile,trim(fnsig),'read',iret=iret) - call nemsio_getfilehead(gfile,iret=iret - + ,idate=idate_nems(1:7),nfhour=nfhour - + ,idvc=idvc,idsl=idsl,lat=dum1d,lon=dum1d2 - + ,vcoord=vcoordnems) - - do k=1,levs+1 - vcoord(k,1)=vcoordnems(k,1,1) - vcoord(k,2)=vcoordnems(k,2,1) - end do - idate(1)=idate_nems(4) - idate(2)=idate_nems(2) - idate(3)=idate_nems(3) - idate(4)=idate_nems(1) - fhour=float(nfhour) - print *, ' processing forecast hour ', fhour - print *, ' idate =', idate - jdate = idate(4)*1000000 + idate(2)*10000+ - & idate(3)*100 + idate(1) - print *, 'jdate = ', jdate - print *, 'Total number of stations = ', npoint - ap = 0.0 - do j=1,jm - do i=1,im - gdlat(i,j)=dum1d((j-1)*im+i) - gdlon(i,j)=dum1d2((j-1)*im+i) - end do - end do - endif !end read in nemsio hearder +!!end read NetCDF hearder info + endif !end read in netcdf if(debugprint) then do k=1,levs+1 @@ -238,14 +206,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='hgtsfc' Zreverse='yes' - call read_netcdf_p(ncid,im,jm,1,VarName,hgt,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'surface hgt not found' - else - VarName='hgt' - LayName='sfc' - call read_nemsio(gfile,im,jm,1,VarName,LayName,hgt, - & iope,ionproc,iocomms,error) + call read_netcdf(ncid,im,jm,1,VarName,hgt,Zreverse, + & error) if (error /= 0) print*,'surface hgt not found' endif if(debugprint)print*,'sample sfc h= ',hgt(im/5,jm/4) @@ -255,15 +217,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='pressfc' Zreverse='yes' - call read_netcdf_p(ncid,im,jm,1,VarName,pint(:,:,1), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'surface pressure not found' - else - VarName='pres' - LayName='sfc' - call read_nemsio(gfile,im,jm,1,VarName, - & LayName,pint(:,:,1),error) + call read_netcdf(ncid,im,jm,1,VarName,pint(:,:,1), + & Zreverse,error) if (error /= 0) print*,'surface pressure not found' endif if(debugprint)print*,'sample sfc P= ',pint(im/2,jm/4,1), @@ -273,13 +228,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='tmp' Zreverse='yes' - call read_netcdf_p(ncid,im,jm,levs,VarName,t3d,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'temp not found' - else - VarName='tmp' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName,t3d,error) + call read_netcdf(ncid,im,jm,levs,VarName,t3d,Zreverse, + & error) if (error /= 0) print*,'temp not found' endif if(debugprint) then @@ -292,13 +242,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='spfh' Zreverse='yes' - call read_netcdf_p(ncid,im,jm,levs,VarName,q3d,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'spfh not found' - else - VarName='spfh' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName,q3d,error) + call read_netcdf(ncid,im,jm,levs,VarName,q3d,Zreverse, + & error) if (error /= 0) print*,'spfh not found' endif if(debugprint) then @@ -311,13 +256,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='ugrd' Zreverse='yes' - call read_netcdf_p(ncid,im,jm,levs,VarName,uh,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'ugrd not found' - else - VarName='ugrd' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName,uh,error) + call read_netcdf(ncid,im,jm,levs,VarName,uh,Zreverse, + & error) if (error /= 0) print*,'ugrd not found' endif if(debugprint) then @@ -330,13 +270,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='vgrd' Zreverse='yes' - call read_netcdf_p(ncid,im,jm,levs,VarName,vh,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'vgrd not found' - else - VarName='vgrd' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName,vh,error) + call read_netcdf(ncid,im,jm,levs,VarName,vh,Zreverse, + & error) if (error /= 0) print*,'vgrd not found' endif if(debugprint) then @@ -349,14 +284,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='dzdt' Zreverse='yes' - call read_netcdf_p(ncid,im,jm,levs,VarName,omega3d,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'dzdt not found' - else - VarName='dzdt' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName, - & omega3d,error) + call read_netcdf(ncid,im,jm,levs,VarName,omega3d,Zreverse, + & error) if (error /= 0) print*,'dzdt not found' endif if(debugprint) then @@ -369,14 +298,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='dpres' Zreverse='no' - call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'dpres not found' - else - VarName='dpres' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName, - & delpz,error) + call read_netcdf(ncid,im,jm,levs,VarName,delpz,Zreverse, + & error) if (error /= 0) print*,'dpres not found' endif if(debugprint) then @@ -431,13 +354,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='delz' Zreverse='no' - call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'delz not found' - else - VarName='delz' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName,delpz,error) + call read_netcdf(ncid,im,jm,levs,VarName,delpz,Zreverse, + & error) if (error /= 0) print*,'delz not found' endif if(debugprint) then @@ -493,51 +411,18 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, nf1 = nf - nint3 endif if ( nf == 0 ) nf1=0 - if(nf==0) then - fngrib='flxf00' - elseif(nf.lt.10) then - fngrib='flxf0' - write(fngrib(6:6),'(i1)') nf - elseif(nf.lt.100) then - fngrib='flxf' - write(fngrib(5:6),'(i2)') nf - else - fngrib='flxf' - write(fngrib(5:7),'(i3)') nf - endif - if(nf1==0) then - fngrib2='flxf00' - elseif(nf1.lt.10) then - fngrib2='flxf0' - write(fngrib2(6:6),'(i1)') nf1 - elseif(nf1.lt.100) then - fngrib2='flxf' - write(fngrib2(5:6),'(i2)') nf1 - else - fngrib2='flxf' - write(fngrib2(5:7),'(i3)') nf1 - endif if (fformat == 'netcdf') then error=nf90_open(trim(fngrib),nf90_nowrite,ncid) !open T-nint below error=nf90_open(trim(fngrib2),nf90_nowrite,ncid2) if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) - 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) endif ! land water mask if (fformat == 'netcdf') then VarName='land' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,lwmask,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'lwmask not found' - else - VarName='land' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName,lwmask,error) + call read_netcdf(ncid,im,jm,1,VarName,lwmask,Zreverse, + & error) if (error /= 0) print*,'lwmask not found' endif if(debugprint) @@ -548,15 +433,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='tmpsfc' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,1), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'tmpsfc not found' - else - VarName='tmp' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - & dum2d(:,:,1),error) + call read_netcdf(ncid,im,jm,1,VarName,dum2d(:,:,1), + & Zreverse,error) if (error /= 0) print*,'tmpsfc not found' endif if(debugprint) @@ -566,15 +444,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='tmp2m' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,2), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'tmp2m not found' - else - VarName='tmp' - LayName='2 m above gnd' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,2),error) + call read_netcdf(ncid,im,jm,1,VarName,dum2d(:,:,2), + & Zreverse,error) if (error /= 0) print*,'tmp2m not found' endif if(debugprint) @@ -585,15 +456,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='spfh2m' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,3), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'spfh2m not found' - else - VarName='spfh' - LayName='2 m above gnd' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,3),error) + call read_netcdf(ncid,im,jm,1,VarName,dum2d(:,:,3), + & Zreverse,error) if (error /= 0) print*,'spfh2m not found' endif if(debugprint) @@ -604,15 +468,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='ugrd10m' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,4), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'ugrd10m not found' - else - VarName='ugrd' - LayName='10 m above gnd' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,4),error) + call read_netcdf(ncid,im,jm,1,VarName,dum2d(:,:,4), + & Zreverse,error) if (error /= 0) print*,'ugrd10m not found' endif @@ -620,33 +477,45 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='vgrd10m' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,5), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'vgrd10m not found' - else - VarName='vgrd' - LayName='10 m above gnd' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,5),error) + call read_netcdf(ncid,im,jm,1,VarName,dum2d(:,:,5), + & Zreverse,error) if (error /= 0) print*,'vgrd10m not found' endif ! soil T if (fformat == 'netcdf') then - VarName='soilt1' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,6), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'soilt1 not found' - else - VarName='tmp' - LayName='0-10 cm down' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,6),error) - if (error /= 0) print*,'soil T not found' - endif + + error=nf90_inq_varid(ncid,'zsoil',dimid) + + if(error/=0)then !read soil avriables in 2D + + VarName='soilt1' + Zreverse='no' + call read_netcdf(ncid,im,jm,1,VarName,dum2d(:,:,6), + & Zreverse,error) + if (error /= 0) print*,'soilt1 not found' + + else !read soil variables in 3D + + error=nf90_inq_dimid(ncid,"zsoil",dimid) + error=nf90_inquire_dimension(ncid,dimid,dim_nam,nsoil) + + print*,'NetCDF file 3D soil = ',error, nsoil + allocate(soilt3d(im,jm,nsoil)) + + VarName='soilt' + Zreverse='no' + call read_netcdf(ncid,im,jm,nsoil,VarName,soilt3d, + & Zreverse,error) + if (error /= 0) print*,'3D soil temp not found' + dum2d(1:im,1:jm,6)=soilt3d(1:im,1:jm,1) +! print*,'soilt3d(im/2,jm/2,:)= ', soilt3d(im/2,jm/2,:) +! print*,'dum2d(im/2,jm/2,6)= ', dum2d(im/2,jm/2,6) + deallocate(soilt3d) + endif + + endif + if(debugprint) + print*,'sample soil T= ',dum2d(im/2,jm/4,6),dum2d(im/2,jm/3,6), + dum2d(im/2,jm/2,6) @@ -655,15 +524,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='snod' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,7), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'snod not found' - else - VarName='snod' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,7),error) + call read_netcdf(ncid,im,jm,1,VarName,dum2d(:,:,7), + & Zreverse,error) if (error /= 0) print*,'snod not found' endif @@ -672,15 +534,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='lhtfl' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,8), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'lhtfl not found' - else - VarName='lhtfl' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,8),error) + call read_netcdf(ncid,im,jm,1,VarName,dum2d(:,:,8), + & Zreverse,error) if (error /= 0) print*,'lhtfl not found' endif if(debugprint) @@ -700,21 +555,12 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='prate_ave' Zreverse='no' -!! call read_netcdf_p(ncid,im,jm,1,VarName,apcp,Zreverse,error) !current hour - call read_netcdf_p(ncid,im,jm,1,VarName,apcp,Zreverse, - & iope,ionproc,iocomms,error) -!! call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse,error) !earlier hour - call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, - & iope,ionproc,iocomms,error) +!! call read_netcdf(ncid,im,jm,1,VarName,apcp,Zreverse,error) !current hour + call read_netcdf(ncid,im,jm,1,VarName,apcp,Zreverse, + & error) +!! call read_netcdf(ncid2,im,jm,1,VarName,cpcp,Zreverse,error) !earlier hour + call read_netcdf(ncid2,im,jm,1,VarName,cpcp,Zreverse,error) if (error /= 0) print*,'prate_ave not found' - else - VarName='prate_ave' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + apcp,error) - call read_nemsio(ffile2,im,jm,1,VarName,LayName, - + cpcp,error) - if (error /= 0) print*,'prate_ave2 not found' endif if(debugprint) & print*,'sample fhour ,3= ', fhour, @@ -735,19 +581,11 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='cprat_ave' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,apcp,Zreverse, - & iope,ionproc,iocomms,error) - call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, - & iope,ionproc,iocomms,error) + call read_netcdf(ncid,im,jm,1,VarName,apcp,Zreverse, + & error) + call read_netcdf(ncid2,im,jm,1,VarName,cpcp,Zreverse, + & error) if (error /= 0) print*,'cprat_ave not found' - else - VarName='cprat_ave' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + apcp,error) - call read_nemsio(ffile2,im,jm,1,VarName,LayName, - + cpcp,error) - if (error /= 0) print*,'cprat_ave2 not found' endif ap=fhour-fint do j=1,jm @@ -761,15 +599,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='weasd' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,11), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'weasd not found' - else - VarName='weasd' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,11),error) + call read_netcdf(ncid,im,jm,1,VarName,dum2d(:,:,11), + & Zreverse,error) if (error /= 0) print*,'weasd not found' endif @@ -777,48 +608,27 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then VarName='tcdc_avelcl' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName, - & dum2d(:,:,12),Zreverse, - & iope,ionproc,iocomms,error) + call read_netcdf(ncid,im,jm,1,VarName, + & dum2d(:,:,12),Zreverse,error) if (error /= 0) print*,'tcdc_avelcl not found' - else - VarName='tcdc_ave' - LayName='low cld lay' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,12),error) - if (error /= 0) print*,'low cld lay not found' endif ! mid cloud fraction if (fformat == 'netcdf') then VarName='tcdc_avemcl' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName, - & dum2d(:,:,13),Zreverse, - & iope,ionproc,iocomms,error) + call read_netcdf(ncid,im,jm,1,VarName, + & dum2d(:,:,13),Zreverse,error) if (error /= 0) print*,'tcdc_avemcl not found' - else - VarName='tcdc_ave' - LayName='mid cld lay' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,13),error) - if (error /= 0) print*,'mid cld lay not found' endif ! high cloud fraction if (fformat == 'netcdf') then VarName='tcdc_avehcl' Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName, - & dum2d(:,:,14),Zreverse, - & iope,ionproc,iocomms,error) + call read_netcdf(ncid,im,jm,1,VarName, + & dum2d(:,:,14),Zreverse,error) if (error /= 0) print*,'tcdc_avehcl not found' - else - VarName='tcdc_ave' - LayName='high cld lay' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,14),error) - if (error /= 0) print*,'high cld lay not found' endif if(debugprint) @@ -828,9 +638,6 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, if (fformat == 'netcdf') then error=nf90_close(ncid) error=nf90_close(ncid2) - else - call nemsio_close(ffile,iret=error) - call nemsio_close(ffile2,iret=error) endif call date_and_time(date,time,zone,clocking) ! print *,'10reading surface data end= ', clocking @@ -1109,7 +916,6 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! ! prepare buffer data ! - if(iope == 0) then call gfuncphys do np = 1, npoint pi3(np,1)=psn(np)*1000 @@ -1308,7 +1114,6 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, !! write(nfile) data write(nfile) data2 enddo !End loop over stations np - endif call date_and_time(date,time,zone,clocking) ! print *,'13reading write data end= ', clocking print *,'13date, time, zone',date, time, zone diff --git a/src/gfs_bufr.fd/meteorg.f_org b/src/gfs_bufr.fd/meteorg.f_org new file mode 100644 index 00000000..d4fa338d --- /dev/null +++ b/src/gfs_bufr.fd/meteorg.f_org @@ -0,0 +1,1326 @@ + subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, + & nf,nfile,fnsig,jdate,idate, + & levs,im,jm,kdim, + & landwater,nend1,nint1,nint3,iidum,jjdum,np1, + & fformat,iocomms,iope,ionproc) + +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! . . . . +! 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 +! PRESSURE (line 290) +! CORRECT THE TEMPERATURE ADJUSTMENT (line 238) +! 2014-03-27 DANA CARLIS UNIFY CODE WITH GFS FORECAST MODEL PRECIP +! TYPE CALCULATION +! 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 +! AND RECALCULATE HOURLY AND 3 HOURLY OUTPUT DEPENDING +! 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 +! 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: +! npoint - number of points +! rlat(npint) - latitude +! rlon(npoint) - longtitude +! istat(npoint) - station id +! elevstn(npoint) - station elevation (m) +! nf - forecast cycle +! fnsig - sigma file name +! idate(4) - date +! levs - input vertical layers +! kdim - sfc file dimension +! +! OUTPUT: +! nfile - output data file channel +! jdate - date YYYYMMDDHH +! +! ATTRIBUTES: +! LANGUAGE: +! MACHINE: IBM SP +! +!$$$ + use netcdf + use nemsio_module + use sigio_module + use physcons + use mersenne_twister +! use funcphys, only : gfuncphys + use funcphys + implicit none + include 'mpif.h' + type(nemsio_gfile) :: gfile + type(nemsio_gfile) :: ffile + type(nemsio_gfile) :: ffile2 + integer :: nfile,npoint,levs,kdim + integer :: nfile1 + integer :: i,j,im,jm,kk,idum,jdum,idvc,idsl +! idsl Integer(sigio_intkind) semi-lagrangian id +! idvc Integer(sigio_intkind) vertical coordinate id +! (=1 for sigma, =2 for ec-hybrid, =3 for ncep hybrid) + integer,parameter :: nvcoord=2 + integer,parameter :: levso=64 + integer :: idate(4),nij,nflx2,np,k,l,nf,nfhour,np1 + integer :: idate_nems(7) + integer :: iret,jdate,leveta,lm,lp1 + character*512 :: fnsig,fngrib +!! real*8 :: data(6*levs+25) + real*8 :: data2(6*levso+25) + real*8 :: rstat1 + character*8 :: cstat1 + character*4 :: cstat(npoint) + real :: fhour,pp,ppn,qs,qsn,esn,es,psfc,ppi,dtemp,nd + real :: t,q,u,v,td,tlcl,plcl,qw,tw,xlat,xlon + integer,dimension(npoint):: landwater + integer,dimension(im,jm):: lwmask + real,dimension(im,jm):: apcp, cpcp + real,dimension(npoint,2+levs*3):: grids + real,dimension(npoint) :: rlat,rlon,pmsl,ps,psn,elevstn + real,dimension(im*jm) :: dum1d,dum1d2 + real,dimension(im,jm) :: gdlat, hgt, gdlon + real,dimension(im,jm,15) :: dum2d + real,dimension(im,jm,levs) :: t3d, q3d, uh, vh,omega3d + real,dimension(im,jm,levs) :: delpz + real,dimension(im,jm,levs+1) :: pint, zint + real,dimension(npoint,levs) :: gridu,gridv,omega,qnew,zp + real,dimension(npoint,levs) :: p1,pd3,ttnew + real,dimension(npoint,levs) :: z1 + real,dimension(npoint,levs+1) :: pi3 + real :: zp2(2) + real,dimension(kdim,npoint) :: sfc + real,dimension(1,levs+1) :: prsi,phii + real,dimension(1,levs) :: gt0,gq0,prsl,phy_f3d + real :: PREC,TSKIN,SR,randomno(1,2) + real :: DOMR,DOMZR,DOMIP,DOMS + real :: vcoord(levs+1,nvcoord),vdummy(levs+1) + real :: vcoordnems(levs+1,3,2) + real :: rdum + integer :: n3dfercld,iseedl + integer :: istat(npoint) + logical :: trace +!! logical, parameter :: debugprint=.true. + logical, parameter :: debugprint=.false. + character lprecip_accu*3 + real, parameter :: ERAD=6.371E6 + real, parameter :: DTR=3.1415926/180. + real :: ap + integer :: nf1, fint + integer :: nend1, nint1, nint3 + character*150 :: fngrib2 + integer recn_dpres,recn_delz,recn_dzdt + integer :: jrec + equivalence (cstat1,rstat1) + integer iidum(npoint),jjdum(npoint) + integer :: error, ncid, ncid2, id_var,dimid + character(len=100) :: long_name + character(len=6) :: fformat + integer,dimension(8) :: clocking + character(10) :: date + character(12) :: time + character(7) :: zone + character(3) :: Zreverse + character(20) :: VarName,LayName + integer iocomms,iope,ionproc + + nij = 12 +!! nflx = 6 * levs + nflx2 = 6 * levso + recn_dpres = 0 + recn_delz = 0 + recn_dzdt = 0 + jrec = 0 + lprecip_accu='yes' + + idvc=2 + idsl=1 +!read in NetCDF file header info + print*,"fformat= ", fformat +! print*,'meteorg.f, idum,jdum= ' +! do np = 1, npoint +! print*, iidum(np), jjdum(np) +! enddo + + if(fformat .eq. "netcdf") then + print*,'iocomms inside meteorg.f=', iocomms + error=nf90_open(trim(fnsig),ior(nf90_nowrite,nf90_mpiio), + & ncid,comm=iocomms, info = mpi_info_null) + error=nf90_get_att(ncid,nf90_global,"ak",vdummy) + do k = 1, levs+1 + vcoord(k,1)=vdummy(levs-k+1) + enddo + error=nf90_get_att(ncid,nf90_global,"bk",vdummy) + do k = 1, levs+1 + vcoord(k,2)=vdummy(levs-k+1) + enddo + error=nf90_inq_varid(ncid, "time", id_var) + error=nf90_get_var(ncid, id_var, nfhour) + print*, "nfhour:",nfhour + error=nf90_get_att(ncid,id_var,"units",long_name) +!! print*,'time units',' -- ',trim(long_name) + read(long_name(13:16),"(i4)")idate(4) + read(long_name(18:19),"(i2)")idate(2) + read(long_name(21:22),"(i2)")idate(3) + read(long_name(24:25),"(i2)")idate(1) + fhour=float(nfhour) + print*,'date= ', idate + jdate = idate(4)*1000000 + idate(2)*10000+ + & idate(3)*100 + idate(1) + print *, 'jdate = ', jdate + error=nf90_inq_varid(ncid, "lon", id_var) + error=nf90_get_var(ncid, id_var, gdlon) + 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 + + call nemsio_open(gfile,trim(fnsig),'read',iret=iret) + call nemsio_getfilehead(gfile,iret=iret + + ,idate=idate_nems(1:7),nfhour=nfhour + + ,idvc=idvc,idsl=idsl,lat=dum1d,lon=dum1d2 + + ,vcoord=vcoordnems) + + do k=1,levs+1 + vcoord(k,1)=vcoordnems(k,1,1) + vcoord(k,2)=vcoordnems(k,2,1) + end do + idate(1)=idate_nems(4) + idate(2)=idate_nems(2) + idate(3)=idate_nems(3) + idate(4)=idate_nems(1) + fhour=float(nfhour) + print *, ' processing forecast hour ', fhour + print *, ' idate =', idate + jdate = idate(4)*1000000 + idate(2)*10000+ + & idate(3)*100 + idate(1) + print *, 'jdate = ', jdate + print *, 'Total number of stations = ', npoint + ap = 0.0 + do j=1,jm + do i=1,im + gdlat(i,j)=dum1d((j-1)*im+i) + 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 + print*,'vcoord(k,2)= ', k, vcoord(k,2) + end do + print*,'sample lat= ',gdlat(im/5,jm/4) + + ,gdlat(im/5,jm/3),gdlat(im/5,jm/2) + print*,'sample lon= ',gdlon(im/5,jm/4) + + ,gdlon(im/5,jm/3),gdlon(im/5,jm/2) + endif +! topography + if (fformat == 'netcdf') then + VarName='hgtsfc' + Zreverse='yes' + call read_netcdf_p(ncid,im,jm,1,VarName,hgt,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'surface hgt not found' + else + VarName='hgt' + LayName='sfc' + call read_nemsio(gfile,im,jm,1,VarName,LayName,hgt, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'surface hgt not found' + endif + if(debugprint)print*,'sample sfc h= ',hgt(im/5,jm/4) + + ,hgt(im/5,jm/3),hgt(im/5,jm/2) + +! surface pressure (Pa) + if (fformat == 'netcdf') then + VarName='pressfc' + Zreverse='yes' + call read_netcdf_p(ncid,im,jm,1,VarName,pint(:,:,1), + & Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'surface pressure not found' + else + VarName='pres' + LayName='sfc' + call read_nemsio(gfile,im,jm,1,VarName, + & LayName,pint(:,:,1),error) + if (error /= 0) print*,'surface pressure not found' + endif + if(debugprint)print*,'sample sfc P= ',pint(im/2,jm/4,1), + + pint(im/2,jm/3,1),pint(im/2,jm/2,1) + +! temperature using NetCDF + if (fformat == 'netcdf') then + VarName='tmp' + Zreverse='yes' + call read_netcdf_p(ncid,im,jm,levs,VarName,t3d,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'temp not found' + else + VarName='tmp' + LayName='mid layer' + call read_nemsio(gfile,im,jm,levs,VarName,LayName,t3d,error) + if (error /= 0) print*,'temp not found' + endif + if(debugprint) then + print*,'sample T at lev=1 to levs ' + do k = 1, levs + print*,k, t3d(im/2,jm/3,k) + enddo + endif +! specific humidity + if (fformat == 'netcdf') then + VarName='spfh' + Zreverse='yes' + call read_netcdf_p(ncid,im,jm,levs,VarName,q3d,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'spfh not found' + else + VarName='spfh' + LayName='mid layer' + call read_nemsio(gfile,im,jm,levs,VarName,LayName,q3d,error) + if (error /= 0) print*,'spfh not found' + endif + if(debugprint) then + print*,'sample Q at lev=1 to levs ' + do k = 1, levs + print*,k, q3d(im/2,jm/3,k) + enddo + endif +! U wind + if (fformat == 'netcdf') then + VarName='ugrd' + Zreverse='yes' + call read_netcdf_p(ncid,im,jm,levs,VarName,uh,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'ugrd not found' + else + VarName='ugrd' + LayName='mid layer' + call read_nemsio(gfile,im,jm,levs,VarName,LayName,uh,error) + if (error /= 0) print*,'ugrd not found' + endif + if(debugprint) then + print*,'sample U at lev=1 to levs ' + do k = 1, levs + print*,k, uh(im/2,jm/3,k) + enddo + endif +! V wind + if (fformat == 'netcdf') then + VarName='vgrd' + Zreverse='yes' + call read_netcdf_p(ncid,im,jm,levs,VarName,vh,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'vgrd not found' + else + VarName='vgrd' + LayName='mid layer' + call read_nemsio(gfile,im,jm,levs,VarName,LayName,vh,error) + if (error /= 0) print*,'vgrd not found' + endif + if(debugprint) then + print*,'sample V at lev=1 to levs ' + do k = 1, levs + print*,k, vh(im/2,jm/3,k) + enddo + endif +! dzdt !added by Guang Ping Lou for FV3GFS + if (fformat == 'netcdf') then + VarName='dzdt' + Zreverse='yes' + call read_netcdf_p(ncid,im,jm,levs,VarName,omega3d,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'dzdt not found' + else + VarName='dzdt' + LayName='mid layer' + call read_nemsio(gfile,im,jm,levs,VarName,LayName, + & omega3d,error) + if (error /= 0) print*,'dzdt not found' + endif + if(debugprint) then + print*,'sample dzdt at lev=1 to levs ' + do k = 1, levs + print*,k, omega3d(im/2,jm/3,k) + enddo + endif +! dpres !added by Guang Ping Lou for FV3GFS (interface pressure delta) + if (fformat == 'netcdf') then + VarName='dpres' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'dpres not found' + else + VarName='dpres' + LayName='mid layer' + call read_nemsio(gfile,im,jm,levs,VarName,LayName, + & delpz,error) + if (error /= 0) print*,'dpres not found' + endif + if(debugprint) then + print*,'sample delp at lev=1 to levs ' + do k = 1, levs + print*,k, delpz(im/2,jm/3,k) + enddo + endif +! compute interface pressure + if(recn_dpres == -9999) then + do k=2,levs+1 + do j=1,jm + do i=1,im + pint(i,j,k)=vcoord(k,1) + + +vcoord(k,2)*pint(i,j,1) + end do + end do + end do + else +! compute pint using dpres from top down if DZDT is used + if (fformat == 'netcdf') then + do j=1,jm + do i=1,im + pint(i,j,levs+1) = delpz(i,j,1) + end do + end do + do k=levs,2,-1 + kk=levs-k+2 + do j=1,jm + do i=1,im + pint(i,j,k) = pint(i,j,k+1) + delpz(i,j,kk) + end do + end do + end do + else + do k=2,levs+1 + do j=1,jm + do i=1,im + pint(i,j,k) = pint(i,j,k-1) - delpz(i,j,k-1) + end do + end do + end do + endif + if(debugprint) then + print*,'sample interface pressure pint at lev =1 to levs ' + do k = 1, levs+1 + print*,k, pint(im/2,jm/3,k),pint(im/3,jm/8,k) + enddo + endif + endif +! delz !added by Guang Ping Lou for FV3GFS ("height thickness" with unit "meters" bottom up) + if (fformat == 'netcdf') then + VarName='delz' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'delz not found' + else + VarName='delz' + LayName='mid layer' + call read_nemsio(gfile,im,jm,levs,VarName,LayName,delpz,error) + if (error /= 0) print*,'delz not found' + endif + if(debugprint) then + print*,'sample delz at lev=1 to levs ' + do k = 1, levs + print*,k, delpz(im/2,jm/3,k) + enddo + endif + +! compute interface height (meter) + if(recn_delz == -9999) then + print*, 'using calculated height' + else +! compute zint using delz from bot up if DZDT is used + if (fformat == 'netcdf') then + do j=1,jm + do i=1,im + zint(i,j,1) = 0.0 + end do + end do + do k=2,levs+1 + kk=levs-k+1 + do j=1,jm + do i=1,im + zint(i,j,k) = zint(i,j,k-1) - delpz(i,j,kk) + end do + end do + end do + else + do k=2,levs+1 + do j=1,jm + do i=1,im + zint(i,j,k) = zint(i,j,k-1) + delpz(i,j,k-1) + end do + end do + end do + endif + if(debugprint) then + print*,'sample interface height zint at lev =1 to levs ' + do k = 1, levs+1 + print*,k, zint(im/2,jm/3,k),zint(im/3,jm/8,k) + enddo + endif + endif + +! close up this NetCDF file + error=nf90_close(ncid) + +! Now open up NetCDF surface files + if ( nf .le. nend1 ) then + nf1 = nf - nint1 + else + nf1 = nf - nint3 + endif + if ( nf == 0 ) nf1=0 + if(nf==0) then + fngrib='flxf00' + elseif(nf.lt.10) then + fngrib='flxf0' + write(fngrib(6:6),'(i1)') nf + elseif(nf.lt.100) then + fngrib='flxf' + write(fngrib(5:6),'(i2)') nf + else + fngrib='flxf' + write(fngrib(5:7),'(i3)') nf + endif + if(nf1==0) then + fngrib2='flxf00' + elseif(nf1.lt.10) then + fngrib2='flxf0' + write(fngrib2(6:6),'(i1)') nf1 + elseif(nf1.lt.100) then + fngrib2='flxf' + write(fngrib2(5:6),'(i2)') nf1 + else + fngrib2='flxf' + write(fngrib2(5:7),'(i3)') nf1 + endif + if (fformat == 'netcdf') then + error=nf90_open(trim(fngrib),nf90_nowrite,ncid) +!open T-nint below + error=nf90_open(trim(fngrib2),nf90_nowrite,ncid2) + if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) + 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) + endif +! land water mask + if (fformat == 'netcdf') then + VarName='land' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,lwmask,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'lwmask not found' + else + VarName='land' + LayName='sfc' + call read_nemsio(ffile,im,jm,1,VarName,LayName,lwmask,error) + if (error /= 0) print*,'lwmask not found' + endif + if(debugprint) + + print*,'sample land mask= ',lwmask(im/2,jm/4), + + lwmask(im/2,jm/3) + +! surface T + if (fformat == 'netcdf') then + VarName='tmpsfc' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,1), + & Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'tmpsfc not found' + else + VarName='tmp' + LayName='sfc' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + & dum2d(:,:,1),error) + if (error /= 0) print*,'tmpsfc not found' + endif + if(debugprint) + + print*,'sample sfc T= ',dum2d(im/2,jm/4,1),dum2d(im/2,jm/3,1), + + dum2d(im/2,jm/2,1) +! 2m T + if (fformat == 'netcdf') then + VarName='tmp2m' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,2), + & Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'tmp2m not found' + else + VarName='tmp' + LayName='2 m above gnd' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,2),error) + if (error /= 0) print*,'tmp2m not found' + endif + if(debugprint) + + print*,'sample 2m T= ',dum2d(im/2,jm/4,2),dum2d(im/2,jm/3,2), + + dum2d(im/2,jm/2,2) + +! 2m Q + if (fformat == 'netcdf') then + VarName='spfh2m' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,3), + & Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'spfh2m not found' + else + VarName='spfh' + LayName='2 m above gnd' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,3),error) + if (error /= 0) print*,'spfh2m not found' + endif + if(debugprint) + + print*,'sample 2m Q= ',dum2d(im/2,jm/4,3),dum2d(im/2,jm/3,3), + + dum2d(im/2,jm/2,3) + +! U10 + if (fformat == 'netcdf') then + VarName='ugrd10m' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,4), + & Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'ugrd10m not found' + else + VarName='ugrd' + LayName='10 m above gnd' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,4),error) + if (error /= 0) print*,'ugrd10m not found' + endif + +! V10 + if (fformat == 'netcdf') then + VarName='vgrd10m' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,5), + & Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'vgrd10m not found' + else + VarName='vgrd' + LayName='10 m above gnd' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,5),error) + if (error /= 0) print*,'vgrd10m not found' + endif + +! soil T + if (fformat == 'netcdf') then + VarName='soilt1' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,6), + & Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'soilt1 not found' + else + VarName='tmp' + LayName='0-10 cm down' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,6),error) + if (error /= 0) print*,'soil T not found' + endif + if(debugprint) + + print*,'sample soil T= ',dum2d(im/2,jm/4,6),dum2d(im/2,jm/3,6), + + dum2d(im/2,jm/2,6) + +! snow depth + if (fformat == 'netcdf') then + VarName='snod' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,7), + & Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'snod not found' + else + VarName='snod' + LayName='sfc' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,7),error) + if (error /= 0) print*,'snod not found' + endif + +! evaporation +!instantaneous surface latent heat net flux + if (fformat == 'netcdf') then + VarName='lhtfl' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,8), + & Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'lhtfl not found' + else + VarName='lhtfl' + LayName='sfc' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,8),error) + if (error /= 0) print*,'lhtfl not found' + endif + if(debugprint) + + print*,'evaporation latent heat net flux= ', + + dum2d(im/2,jm/4,8),dum2d(im/2,jm/3,8) + if(debugprint) + + print*,'evaporation latent heat net flux stn 000692)= ', + + dum2d(2239,441,8) + +! total precip + if ( nf .le. nend1 ) then + fint = nint1 + else + fint = nint3 + endif +! for accumulated precipitation: + if (fformat == 'netcdf') then + VarName='prate_ave' + Zreverse='no' +!! call read_netcdf_p(ncid,im,jm,1,VarName,apcp,Zreverse,error) !current hour + call read_netcdf_p(ncid,im,jm,1,VarName,apcp,Zreverse, + & iope,ionproc,iocomms,error) +!! call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse,error) !earlier hour + call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'prate_ave not found' + else + VarName='prate_ave' + LayName='sfc' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + apcp,error) + call read_nemsio(ffile2,im,jm,1,VarName,LayName, + + cpcp,error) + if (error /= 0) print*,'prate_ave2 not found' + endif + if(debugprint) + & print*,'sample fhour ,3= ', fhour, + & '1sample precip rate= ',apcp(im/2,jm/3),cpcp(im/2,jm/3) + ap=fhour-fint + do j=1,jm + do i=1,im + dum2d(i,j,9) =(apcp(i,j)*fhour-cpcp(i,j)*ap)*3600.0 + end do + end do + + if(debugprint) + & print*,'sample fhour ,5= ', fhour, + & 'sample total precip= ',dum2d(im/2,jm/4,9), + + dum2d(im/2,jm/3,9),dum2d(im/2,jm/2,9) + +! convective precip + if (fformat == 'netcdf') then + VarName='cprat_ave' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,apcp,Zreverse, + & iope,ionproc,iocomms,error) + call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'cprat_ave not found' + else + VarName='cprat_ave' + LayName='sfc' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + apcp,error) + call read_nemsio(ffile2,im,jm,1,VarName,LayName, + + cpcp,error) + if (error /= 0) print*,'cprat_ave2 not found' + endif + ap=fhour-fint + 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 + +! water equi + if (fformat == 'netcdf') then + VarName='weasd' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,11), + & Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'weasd not found' + else + VarName='weasd' + LayName='sfc' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,11),error) + if (error /= 0) print*,'weasd not found' + endif + +! low cloud fraction + if (fformat == 'netcdf') then + VarName='tcdc_avelcl' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName, + & dum2d(:,:,12),Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'tcdc_avelcl not found' + else + VarName='tcdc_ave' + LayName='low cld lay' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,12),error) + if (error /= 0) print*,'low cld lay not found' + endif + +! mid cloud fraction + if (fformat == 'netcdf') then + VarName='tcdc_avemcl' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName, + & dum2d(:,:,13),Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'tcdc_avemcl not found' + else + VarName='tcdc_ave' + LayName='mid cld lay' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,13),error) + if (error /= 0) print*,'mid cld lay not found' + endif + +! high cloud fraction + if (fformat == 'netcdf') then + VarName='tcdc_avehcl' + Zreverse='no' + call read_netcdf_p(ncid,im,jm,1,VarName, + & dum2d(:,:,14),Zreverse, + & iope,ionproc,iocomms,error) + if (error /= 0) print*,'tcdc_avehcl not found' + else + VarName='tcdc_ave' + LayName='high cld lay' + call read_nemsio(ffile,im,jm,1,VarName,LayName, + + dum2d(:,:,14),error) + if (error /= 0) print*,'high cld lay not found' + endif + + if(debugprint) + + print*,'sample high cloud frac= ',dum2d(im/2,jm/4,14), + + dum2d(im/2,jm/3,14),dum2d(im/2,jm/2,14) + + if (fformat == 'netcdf') then + error=nf90_close(ncid) + error=nf90_close(ncid2) + else + call nemsio_close(ffile,iret=error) + call nemsio_close(ffile2,iret=error) + endif + call date_and_time(date,time,zone,clocking) +! print *,'10reading surface data end= ', clocking + print *,'10date, time, zone',date, time, zone +! +! get the nearest neighbor i,j from the table +! + do np=1, npoint +! use read in predetermined i,j + if (np1==0) then + idum=iidum(np) + jdum=jjdum(np) + + else +! find nearest neighbor + rdum=rlon(np) + if(rdum<0.)rdum=rdum+360. + + do j=1,jm-1 + do i=1,im-1 + if((rdum>=gdlon(i,j) .and. rdum<=gdlon(i+1,j)) .and. + + (rlat(np)<=gdlat(i,j).and.rlat(np)>=gdlat(i,j+1)) ) then + if(landwater(np) == 2)then + idum=i + jdum=j + exit + else if(landwater(np) == lwmask(i,j))then + idum=i + jdum=j !1 + exit + else if(landwater(np) == lwmask(i+1,j))then + idum=i+1 + jdum=j ! 2 + exit + else if(landwater(np) == lwmask(i-1,j))then + idum=i-1 + jdum=j ! 3 + exit + else if(landwater(np) == lwmask(i,j+1))then + idum=i + jdum=j+1 ! 4 + exit + else if(landwater(np) == lwmask(i,j-1))then + idum=i + jdum=j-1 ! 5 + exit + else if(landwater(np) == lwmask(i+1,j-1))then + idum=i+1 + jdum=j-1 ! 6 + exit + else if(landwater(np) == lwmask(i+1,j+1))then + idum=i+1 + jdum=j+1 ! 7 + exit + else if(landwater(np) == lwmask(i-1,j+1))then + idum=i-1 + jdum=j+1 ! 8 + exit + else if(landwater(np) == lwmask(i-1,j-1))then + idum=i-1 + jdum=j-1 ! 9 + exit + else if(landwater(np) == lwmask(i,j+2))then + idum=i + jdum=j+2 ! 10 + exit + else if(landwater(np) == lwmask(i+2,j))then + idum=i+2 + jdum=j !11 + exit + else if(landwater(np) == lwmask(i,j-2))then + idum=i + jdum=j-2 ! 12 + exit + else if(landwater(np) == lwmask(i-2,j))then + idum=i-2 + jdum=j !13 + exit + else if(landwater(np) == lwmask(i-2,j+1))then + idum=i-2 + jdum=j+1 ! 14 + exit + else if(landwater(np) == lwmask(i-1,j+2))then + idum=i-1 + jdum=j+2 !15 + exit + else if(landwater(np) == lwmask(i+1,j+2))then + idum=i+1 + jdum=j+2 !16 + exit + else if(landwater(np) == lwmask(i+2,j+1))then + idum=i+2 + jdum=j+1 !17 + exit + else if(landwater(np) == lwmask(i+2,j-1))then + idum=i+2 + jdum=j-1 !18 + exit + else if(landwater(np) == lwmask(i+1,j-2))then + idum=i+1 + jdum=j-2 !19 + exit + else if(landwater(np) == lwmask(i-1,j-2))then + idum=i-1 + jdum=j-2 !20 + exit + else if(landwater(np) == lwmask(i-2,j-1))then + idum=i-2 + jdum=j-1 !21 + exit + else if(landwater(np) == lwmask(i-2,j-2))then + idum=i-2 + jdum=j-2 !22 + exit + else if(landwater(np) == lwmask(i+2,j-2))then + idum=i+2 + jdum=j-2 !23 + exit + else if(landwater(np) == lwmask(i+2,j+2))then + idum=i+2 + jdum=j+2 !24 + exit + else if(landwater(np) == lwmask(i-2,j+2))then + idum=i-2 + jdum=j+2 !25 + exit + else if(landwater(np) == lwmask(i+3,j))then + idum=i+3 + jdum=j !26 + exit + else if(landwater(np) == lwmask(i-3,j))then + idum=i-3 + jdum=j !27 + exit + else if(landwater(np) == lwmask(i,j+3))then + idum=i + jdum=j+3 !28 + exit + else if(landwater(np) == lwmask(i,j-3))then + idum=i + jdum=j-3 !29 + exit + else +CC print*,'no matching land sea mask np,landwater,i,j,mask= ' +CC print*, np,landwater(np),i,j,lwmask(i,j) +CC print*, ' So it takes i,j ' + idum=i + jdum=j + exit + end if + end if + end do + end do + + idum=max0(min0(idum,im),1) + jdum=max0(min0(jdum,jm),1) + endif !! read in i,j ends here + if (fhour==0.0) then + if(debugprint) then + write(nij,98) np,idum,jdum,rlat(np),rlon(np) + 98 FORMAT (3I6, 2F9.2) + if(elevstn(np)==-999.) elevstn(np)=hgt(idum,jdum) + write(9,99) np,rlat(np),rlon(np),elevstn(np),hgt(idum,jdum) + 99 FORMAT (I6, 4F9.2) + if(np==1 .or.np==100)print*,'nearest neighbor for station ',np + + ,idum,jdum,rlon(np),rlat(np),lwmask(i,j),landwater(np) + endif + endif + + 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) + sfc(12,np)=dum2d(idum,jdum,9) + sfc(11,np)=dum2d(idum,jdum,10) + sfc(10,np)=dum2d(idum,jdum,11) + sfc(27,np)=dum2d(idum,jdum,12) + sfc(26,np)=dum2d(idum,jdum,13) + sfc(25,np)=dum2d(idum,jdum,14) + sfc(34,np)=dum2d(idum,jdum,4) + sfc(35,np)=dum2d(idum,jdum,5) + sfc(30,np)=dum2d(idum,jdum,2) + sfc(31,np)=dum2d(idum,jdum,3) + +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) + + do k=1,levs + 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) + gridv(np,k)=vh(idum,jdum,k) + p1(np,k)=pint(idum,jdum,k+1) + 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 +! !ps in kPa + ps(np) = grids(np,2)/1000. !! surface pressure + enddo + +! +! ----------------- +! Put topo(1),surf press(2),vir temp(3:66),and specifi hum(67:130) in grids +! for each station +!! if(recn_dzdt == 0 ) then !!DZDT + do k = 1, levs + do np = 1, npoint + omega(np,k) = grids(np,2+levs*2+k) + enddo + enddo + if(debugprint) + + print*,'sample (omega) dzdt ', (omega(3,k),k=1,levs) +! +! move surface pressure to the station surface from the model surface +! + do np = 1, npoint +! +! when the station elevation information in the table says missing, +! use the model elevation +! +! print *, "elevstn = ", elevstn(np) + if(elevstn(np)==-999.) elevstn(np) = grids(np,1) + psn(np) = ps(np) + call sigio_modpr(1,1,levs,nvcoord,idvc, + & idsl,vcoord,iret, + & 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) + if(np==11)print*,'pd3= ', pd3(np,1:levs) + enddo +! +!! test removing height adjustments + print*, 'do not do height adjustments' +! +! get sea-level pressure (Pa) and layer geopotential height +! + do k = 1, levs + do np = 1, npoint + ttnew(np,k) = grids(np,k+2) + qnew(np,k) = grids(np,k+levs+2) + enddo + enddo + + do np=1,npoint +!! call gslp(levs,elevstn(np),ps(np)*1000, + call gslp(levs,grids(np,1),ps(np)*1000, + & p1(np,1:levs),ttnew(np,1:levs),qnew(np,1:levs), + & pmsl(np),zp(np,1:levs),zp2(1:2)) + enddo + print *, 'call gslp pmsl= ', (pmsl(np),np=1,20) + if(recn_delz == -9999) then + print*, 'using calculated height ' + else + print*, 'using model height m' + do k = 1, levs + do np=1, npoint + zp(np,k) = z1(np,k) + enddo + enddo + 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) +! +! prepare buffer data +! + if(iope == 0) then + call gfuncphys + do np = 1, npoint + pi3(np,1)=psn(np)*1000 + do k=1,levs + pi3(np,k+1)=pi3(np,k)-pd3(np,k) !layer pressure (Pa) + enddo +!! ==ivalence (cstat1,rstat1) + cstat1=cstat(np) +!! data(1) = ifix(fhour+.2) * 3600 ! FORECAST TIME (SEC) +!! data(2) = istat(np) ! STATION NUMBER +!! data(3) = rstat1 ! STATION CHARACTER ID +!! data(4) = rlat(np) ! LATITUDE (DEG N) +!! data(5) = rlon(np) ! LONGITUDE (DEG E) +!! data(6) = elevstn(np) ! STATION ELEVATION (M) + data2(1) = ifix(fhour+.2) * 3600 ! FORECAST TIME (SEC) + data2(2) = istat(np) ! STATION NUMBER + data2(3) = rstat1 ! STATION CHARACTER ID + data2(4) = rlat(np) ! LATITUDE (DEG N) + data2(5) = rlon(np) ! LONGITUDE (DEG E) + data2(6) = elevstn(np) ! STATION ELEVATION (M) + psfc = 10. * psn(np) ! convert to MB + leveta = 1 + do k = 1, levs + kk= k/2 + 1 +! +! look for the layer above 500 mb for precip type computation +! + if(pi3(np,k).ge.50000.) leveta = k + ppi = pi3(np,k) + t = grids(np,k+2) + q = max(1.e-8,grids(np,2+k+levs)) + u = gridu(np,k) + v = gridv(np,k) +!! data((k-1)*6+7) = p1(np,k) ! PRESSURE (PA) at integer layer +!! data((k-1)*6+8) = t ! TEMPERATURE (K) +!! data((k-1)*6+9) = u ! U WIND (M/S) +!! data((k-1)*6+10) = v ! V WIND (M/S) +!! data((k-1)*6+11) = q ! HUMIDITY (KG/KG) +!! data((k-1)*6+12) = omega(np,k)*100. ! Omega (pa/sec) !changed to dzdt(cm/s) if available + 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+11) = q + data2((kk-1)*6+12) = omega(np,k)*100. + endif +!changed to dzdt(cm/s) if available + enddo +! +! process surface flux file fields +! +!! data(8+nflx) = psfc * 100. ! SURFACE PRESSURE (PA) +!! data(7+nflx) = pmsl(np) + data2(8+nflx2) = psfc * 100. ! SURFACE PRESSURE (PA) + 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(37,np) = data2(7+nflx2) * .01 + 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: +! It was decided that no corrctions were needed due to higher model +! resolution. +! +! if(sfc(10,np)==0.) then +! sfc(30,np) = sfc(30,np) + dtemp +! sfc(5,np) = sfc(5,np) + dtemp +! endif +! if(sfc(10,np).gt.0..and.sfc(5,np).lt.273.16) then +! sfc(5,np) = sfc(5,np) + dtemp +! if(sfc(5,np).gt.273.16) then +! dtemp = sfc(5,np) - 273.16 +! sfc(5,np) = 273.16 +! endif +! sfc(30,np) = sfc(30,np) + dtemp +! endif +! +!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 +! 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 + if(debugprint) + + print*,'evaporation (stn 000692)= ',sfc(17,np) +!! data(9+nflx) = sfc(5,np) ! tsfc (K) +!! data(10+nflx) = sfc(6,np) ! 10cm soil temp (K) +!!! data(11+nflx) = sfc(17,np)/28.3565 ! evaporation (kg/m**2) from (W m-2) +!! data(11+nflx) = sfc(17,np)*0.035265 ! evaporation (kg/m**2) from (W m-2) +!! data(12+nflx) = sfc(12,np) ! total precip (m) +!! data(13+nflx) = sfc(11,np) ! convective precip (m) +!! data(14+nflx) = sfc(10,np) ! water equi. snow (m) +!! data(15+nflx) = sfc(27,np) ! low cloud (%) +!! data(16+nflx) = sfc(26,np) ! mid cloud +!! data(17+nflx) = sfc(25,np) ! high cloud +!! data(18+nflx) = sfc(34,np) ! U10 (m/s) +!! data(19+nflx) = sfc(35,np) ! V10 (m/s) +!! data(20+nflx) = sfc(30,np) ! T2 (K) +!! data(21+nflx) = sfc(31,np) ! Q2 (K) + +!! data(22+nflx) = 0. +!! data(23+nflx) = 0. +!! data(24+nflx) = 0. +!! data(25+nflx) = 0. +!! create 64 level bufr files + data2(9+nflx2) = sfc(5,np) ! tsfc (K) + data2(10+nflx2) = sfc(6,np) ! 10cm soil temp (K) +!! data2(11+nflx2) = sfc(17,np)/28.3565 ! evaporation (kg/m**2) from (W m-2) + data2(11+nflx2) = sfc(17,np)*0.035265 ! evaporation (kg/m**2) from (W m-2) + data2(12+nflx2) = sfc(12,np) ! total precip (m) + data2(13+nflx2) = sfc(11,np) ! convective precip (m) + data2(14+nflx2) = sfc(10,np) ! water equi. snow (m) + data2(15+nflx2) = sfc(27,np) ! low cloud (%) + data2(16+nflx2) = sfc(26,np) ! mid cloud + data2(17+nflx2) = sfc(25,np) ! high cloud + data2(18+nflx2) = sfc(34,np) ! U10 (m/s) + data2(19+nflx2) = sfc(35,np) ! V10 (m/s) + data2(20+nflx2) = sfc(30,np) ! T2 (K) + data2(21+nflx2) = sfc(31,np) ! Q2 (K) + + data2(22+nflx2) = 0. + data2(23+nflx2) = 0. + data2(24+nflx2) = 0. + data2(25+nflx2) = 0. + nd = 0 + trace = .false. + DOMS=0. + DOMR=0. + DOMIP=0. + DOMZR=0. + if(np==1.or.np==2) nd = 1 + if(np==1.or.np==2) trace = .true. + + 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) + t = grids(np,k+2) + q = max(0.,grids(np,2+k+levs)) + u = gridu(np,k) + v = gridv(np,k) + if(q.gt.1.e-6.and.pp.ge.20000.) then + call tdew(td,t,q,pp) + call lcl(tlcl,plcl,t,pp,q) + call mstadb(qw,tw,pp,q,tlcl,plcl) + else + td = t - 30. + tw = t - 30. + endif +! Calpreciptype input variables + gt0(1,k)= t + gq0(1,k) = q + prsl(1,k) = pp + prsi(1,k)=ppi + phii(1,k)=zp(np,k) !height in meters + enddo +! Use GFS routine calpreciptype.f to calculate precip type + xlat=rlat(np) + xlon=rlon(np) + lm=leveta + lp1=leveta+1 +!! PREC=data(12+nflx) + PREC=data2(12+nflx2) + n3dfercld=1 !if =3 then use Ferriers Explicit Precip Type + TSKIN=1. !used in Ferriers Explicit Precip Scheme + SR=1. !used in Ferriers Explicit Precip Scheme + iseedl=jdate + call random_setseed(iseedl) + call random_number(randomno) + call calpreciptype(1,1,1,1,lm,lp1,randomno,xlat,xlon, !input + & gt0,gq0,prsl,prsi,PREC,phii,n3dfercld,TSKIN,SR,phy_f3d, !input + & DOMR,DOMZR,DOMIP,DOMS) ! Output vars + endif +!! data(nflx + 22) = DOMS +!! data(nflx + 23) = DOMIP +!! data(nflx + 24) = DOMZR +!! data(nflx + 25) = DOMR + data2(nflx2 + 22) = DOMS + data2(nflx2 + 23) = DOMIP + data2(nflx2 + 24) = DOMZR + data2(nflx2 + 25) = DOMR + if(np==1.or.np==100) then + print *, ' surface fields for hour', nf, 'np =', np + print *, (data2(l+nflx2),l=1,25) + print *, ' temperature sounding' + print 6101, (data2((k-1)*6+8),k=1,levso) + print *, ' omega sounding' + print *, (data2((k-1)*6+12),k=1,levso) + endif +C print *, 'in meteorg nfile1= ', nfile1 +!! write(nfile) data + write(nfile) data2 + enddo !End loop over stations np + endif + call date_and_time(date,time,zone,clocking) +! print *,'13reading write data end= ', clocking + print *,'13date, time, zone',date, time, zone + print *, 'in meteorg nf,nfile,nfhour= ', nf,nfile,nfhour + print *, 'Finished writing bufr data file' + 6101 format(2x,6f12.3) + 6102 format(2x,6f12.5) + 6103 format(2x,6f12.5) +! + close(unit=nfile) + return + 910 print *, ' error reading surface flux file' + end + +!----------------------------------------------------------------------- From 8aefee28d11a7ec9c04216c7825656ea8adfb7ab Mon Sep 17 00:00:00 2001 From: "Bo.Cui" Date: Mon, 19 Aug 2024 14:42:38 -0400 Subject: [PATCH 2/6] delete files --- src/gfs_bufr.fd/CMakeLists.txt_org | 53 -- src/gfs_bufr.fd/buff.f_org | 91 -- src/gfs_bufr.fd/gfsbufr.f_org | 275 ------ src/gfs_bufr.fd/meteorg.f_org | 1326 ---------------------------- 4 files changed, 1745 deletions(-) delete mode 100644 src/gfs_bufr.fd/CMakeLists.txt_org delete mode 100644 src/gfs_bufr.fd/buff.f_org delete mode 100644 src/gfs_bufr.fd/gfsbufr.f_org delete mode 100644 src/gfs_bufr.fd/meteorg.f_org diff --git a/src/gfs_bufr.fd/CMakeLists.txt_org b/src/gfs_bufr.fd/CMakeLists.txt_org deleted file mode 100644 index 7df490d0..00000000 --- a/src/gfs_bufr.fd/CMakeLists.txt_org +++ /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/buff.f_org b/src/gfs_bufr.fd/buff.f_org deleted file mode 100644 index 82acf036..00000000 --- a/src/gfs_bufr.fd/buff.f_org +++ /dev/null @@ -1,91 +0,0 @@ - subroutine buff(nint1,nend1,nint3,nend3,npoint,idate,jdate,levs, - & dird,lss,istat,sbset,seqflg,clist,npp,wrkd) - character*512 dird, fnbufr, fmto -!! integer nint, nend, npoint, idate(4), levs, jdate - integer nint1, nend1, nint3, nend3 - integer npoint, idate(4), levs, jdate - real*8 data(6*levs+25), wrkd(1) - integer idtln, nf, nfile, np - integer lss, istat(npoint), ios - CHARACTER*8 SBSET - LOGICAL SEQFLG(4) - CHARACTER*80 CLIST(4) - INTEGER NPP(4) - CHARACTER*8 SEQNAM(4) - FMTO = '(A,".",I6.6,".",I10)' - idtln = 8 - nfile = 20 -C print *, 'inside buff.f nint1,nend1,nint3,nend3,jdate=' -C print *, nint1,nend1,nint3,nend3,jdate - do nf = 0, nend1, nint1 - nfile = nfile + 1 - rewind nfile - enddo - do nf = nend1+nint3, nend3, nint3 - nfile = nfile + 1 - rewind nfile - enddo - do np = 1, npoint -C OPEN BUFR OUTPUT FILE. - write(fnbufr,fmto) dird(1:lss),istat(np),jdate - print *, ' fnbufr =', fnbufr - open(unit=19,file=fnbufr,form='unformatted', - & status='new', iostat=ios) - IF ( ios .ne. 0 ) THEN - WRITE (6,*) ' CANNOT open ', 19 - STOP - END IF - CALL OPENBF ( 19, 'OUT', 1 ) - nfile = 20 - do nf = 0, nend1, nint1 - nfile = nfile + 1 - read(nfile) data - if(np.eq.1) then - print *, ' creating bufr file for np, nfile =', - & np, nfile - endif -CC WRITE DATA MESSAGE TO BUFR OUTPUT FILE. -CC LUNTBL=-9 BECAUSE BUFR TABLE FILE NOT USED HERE. -CC SEQNAM=XXXXXX BECAUSE MNEMONIC SEQUENCE NAMES NOT USED HERE. - CALL BFRIZE ( -9, 19, SBSET, - & idate(4), iDATE(2), - & iDATE(3), iDATE(1), - & 'XXXXXX', SEQFLG, 4, .FALSE., DATA, levs, - & CLIST, NPP, WRKD, IRET ) - IF ( IRET .NE. 0 ) THEN - PRINT *,' BFRIZE FAILED ' - ENDIF -c 300 continue - enddo -C 3hourly output starts here -!! print *, 'buff.f nfile,nend1+nint3,nend3,nint3=' -!! print *, nfile,nend1+nint3,nend3,nint3 - do nf = nend1+nint3, nend3, nint3 - nfile = nfile + 1 - read(nfile) data - if(np.eq.1) then - print *, ' creating bufr file for np, nfile =', - & np, nfile - endif -C print *, 'read2 in fort(nfile) =', nfile -CC WRITE DATA MESSAGE TO BUFR OUTPUT FILE. -CC LUNTBL=-9 BECAUSE BUFR TABLE FILE NOT USED HERE. -CC SEQNAM=XXXXXX BECAUSE MNEMONIC SEQUENCE NAMES NOT USED HERE. - CALL BFRIZE ( -9, 19, SBSET, - & idate(4), iDATE(2), - & iDATE(3), iDATE(1), - & 'XXXXXX', SEQFLG, 4, .FALSE., DATA, levs, - & CLIST, NPP, WRKD, IRET ) - IF ( IRET .NE. 0 ) THEN - PRINT *,' BFRIZE FAILED ' - ENDIF - enddo - CALL BFRIZE ( 0, 19, SBSET, - & IDATE(4), IDATE(2), - & IDATE(3), IDATE(1), - & 'XXXXXX', SEQFLG, 4, .FALSE., DATA, levs, - & CLIST, NPP, WRKD, IRET ) - call closbf(19) - enddo - return - end diff --git a/src/gfs_bufr.fd/gfsbufr.f_org b/src/gfs_bufr.fd/gfsbufr.f_org deleted file mode 100644 index c1bb1d35..00000000 --- a/src/gfs_bufr.fd/gfsbufr.f_org +++ /dev/null @@ -1,275 +0,0 @@ - program meteormrf -C$$$ MAIN PROGRAM DOCUMENTATION BLOCK -C -C MAIN PROGRAM: METEOMRF -C PRGMMR: PAN ORG: NP23 DATE: 1999-07-21 -C -C ABSTRACT: Creates BUFR meteogram files for the AVN and MRF. -C -C PROGRAM HISTORY LOG: -C 99-07-21 Hualu Pan -C 16-09-27 HUIYA CHUANG MODIFY TO READ GFS NEMS OUTPUT ON GRID SPACE -C 16-10-15 HUIYA CHUANG: CONSOLIDATE TO READ FLUX FIELDS IN THIS -C PACKAGE TOO AND THIS SPEEDS UP BFS BUFR BY 3X -C 17-02-27 GUANG PING LOU: CHANGE MODEL OUTPUT READ-IN TO HOURLY -C TO 120 HOURS AND 3 HOURLY TO 180 HOURS. -C 19-07-16 GUANG PING LOU: CHANGE FROM NEMSIO TO GRIB2. -C -C -C USAGE: -C INPUT FILES: -C FTxxF001 - UNITS 11 THRU 49 -C PARM - UNIT 5 (STANDARD READ) -C -C OUTPUT FILES: (INCLUDING SCRATCH FILES) -C FTxxF001 - UNITS 51 THRU 79 -C FTxxF001 - UNIT 6 (STANDARD PRINTFILE) -C -C SUBPROGRAMS CALLED: (LIST ALL CALLED FROM ANYWHERE IN CODES) -C UNIQUE: - ROUTINES THAT ACCOMPANY SOURCE FOR COMPILE -C LIBRARY: -C W3LIB - -C -C EXIT STATES: -C COND = 0 - SUCCESSFUL RUN -C =NNNN - TROUBLE OR SPECIAL FLAG - SPECIFY NATURE -C -C REMARKS: LIST CAVEATS, OTHER HELPFUL HINTS OR INFORMATION -C -C ATTRIBUTES: -C LANGUAGE: INDICATE EXTENSIONS, COMPILER OPTIONS -C MACHINE: IBM SP -C -C$$$ - use netcdf - use mpi - use nemsio_module - use sigio_module - implicit none -!! include 'mpif.h' - integer,parameter:: nsta=3000 - integer,parameter:: ifile=11 - integer,parameter:: levso=64 - integer(sigio_intkind):: irets - type(nemsio_gfile) :: gfile - integer ncfsig, nsig - integer istat(nsta), idate(4), jdate - integer :: levs,nstart,nend,nint,nsfc,levsi,im,jm - integer :: npoint,np,ist,is,iret,lss,nss,nf,nsk,nfile - integer :: ielev - integer :: lsfc - real :: alat,alon,rla,rlo - real :: wrkd(1),dummy - real rlat(nsta), rlon(nsta), elevstn(nsta) - integer iidum(nsta),jjdum(nsta) - integer nint1, nend1, nint3, nend3, np1 - integer landwater(nsta) - character*1 ns, ew - character*4 t3 - character*4 cstat(nsta) - character*32 desc - character*512 dird, fnsig - logical f00, makebufr - CHARACTER*8 SBSET - LOGICAL SEQFLG(4) - CHARACTER*80 CLIST(4) - INTEGER NPP(4) - CHARACTER*8 SEQNAM(4) - integer ierr, mrank, msize,ntask - integer n0, ntot - integer :: error, ncid, id_var,dimid - character(len=10) :: dim_nam - character(len=6) :: fformat - !added from Cory - integer :: iope, ionproc - integer, allocatable :: iocomms(:) -C - DATA SBSET / 'ABCD1234' / -C - DATA SEQFLG / .FALSE., .TRUE., .FALSE., .FALSE. / -C - DATA SEQNAM / 'HEADR', 'PROFILE', 'CLS1' ,'D10M' / -c DATA SEQNAM / 'HEADR', 'PRES TMDB UWND VWND SPFH OMEG', -c & 'CLS1' ,'D10M' / -C - namelist /nammet/ levs, makebufr, dird, - & nstart, nend, nint, nend1, nint1, - & nint3, nsfc, f00, fformat, np1 - - call mpi_init(ierr) - call mpi_comm_rank(MPI_COMM_WORLD,mrank,ierr) - call mpi_comm_size(MPI_COMM_WORLD,msize,ierr) - if(mrank.eq.0) then - CALL W3TAGB('METEOMRF',1999,0202,0087,'NP23') - endif - open(5,file='gfsparm') - read(5,nammet) - write(6,nammet) - npoint = 0 - 99 FORMAT (I6, F6.2,A1, F7.2,A1,1X,A4,1X,I2, A28, I4) - do np = 1, nsta+2 - read(8,99,end=200) IST,ALAT,NS,ALON,EW,T3,lsfc,DESC,IELEV -CC print*," IST,ALAT,NS,ALON,EW,T3,lsfc,DESC,IELEV= " -CC print*, IST,ALAT,NS,ALON,EW,T3,lsfc,DESC,IELEV - if(alat.lt.95.) then - npoint = npoint + 1 - RLA = 9999. - IF (NS .EQ. 'N') RLA = ALAT - IF (NS .EQ. 'S') RLA = -ALAT - RLO = 9999. - IF (EW .EQ. 'E') RLO = ALON - IF (EW .EQ. 'W') RLO = -ALON - rlat(npoint) = rla - rlon(npoint) = rlo - istat(npoint) = ist - cstat(npoint) = T3 - elevstn(npoint) = ielev - - if(lsfc .le. 9) then - landwater(npoint) = 2 !!nearest - else if(lsfc .le. 19) then - landwater(npoint) = 1 !!land - else if(lsfc .ge. 20) then - landwater(npoint) = 0 !!water - endif - endif - enddo - 200 continue - if(npoint.le.0) then - print *, ' station list file is empty, abort program' - call abort - elseif(npoint.gt.nsta) then - print *, ' number of station exceeds nsta, abort program' - call abort - endif -! print*,'npoint= ', npoint -! print*,'np,IST,idum,jdum,rlat(np),rlon(np)= ' - if(np1 == 0) then - do np = 1, npoint - read(7,98) IST, iidum(np), jjdum(np), ALAT, ALON - enddo - endif - 98 FORMAT (3I6, 2F9.2) - if (mrank.eq.0.and.makebufr) then - REWIND 1 - READ (1,100) SBSET - 100 FORMAT ( ////// 2X, A8 ) - PRINT 120, SBSET - 120 FORMAT ( ' SBSET=#', A8, '#' ) - REWIND 1 -C -C READ PARM NAMES AND NUMBER OF PARM NAMES FROM BUFR TABLE. - DO IS = 1,4 - CALL BFRHDR ( 1, SEQNAM(IS), SEQFLG(IS), - X CLIST(IS), NPP(IS), IRET ) - IF ( IRET .NE. 0 ) THEN - PRINT*, ' CALL BFRHDR IRET=', IRET - ENDIF - ENDDO - lss = len ( dird ) - DO WHILE ( dird (lss:lss) .eq. ' ' ) - lss = lss - 1 - END DO -C - endif - nsig = 11 - nss = nstart + nint - if(f00) nss = nstart -c do nf = nss, nend, nint - ntot = (nend - nss) / nint + 1 - ntask = mrank/(float(msize)/float(ntot)) - nf = ntask * nint + nss - print*,'n0 ntot nint nss mrank msize' - print*, n0,ntot,nint,nss,mrank,msize - print*,'nf, ntask= ', nf, ntask - if(nf .le. nend1) then - nfile = 21 + (nf / nint1) - else - nfile = 21 + (nend1/nint1) + (nf-nend1)/nint3 - endif - print*, 'nf,nint,nfile = ',nf,nint,nfile - if(nf.le.nend) then - if(nf.lt.10) then - fnsig = 'sigf0' - write(fnsig(6:6),'(i1)') nf - ncfsig = 6 - elseif(nf.lt.100) then - fnsig = 'sigf' - write(fnsig(5:6),'(i2)') nf - ncfsig = 6 - else - fnsig = 'sigf' - write(fnsig(5:7),'(i3)') nf - ncfsig = 7 - endif - print *, 'Opening file : ',fnsig - -!! read in either nemsio or NetCDF files - if (fformat == 'netcdf') then - error=nf90_open(trim(fnsig),nf90_nowrite,ncid) - error=nf90_inq_dimid(ncid,"grid_xt",dimid) - error=nf90_inquire_dimension(ncid,dimid,dim_nam,im) - error=nf90_inq_dimid(ncid,"grid_yt",dimid) - error=nf90_inquire_dimension(ncid,dimid,dim_nam,jm) - error=nf90_inq_dimid(ncid,"pfull",dimid) - error=nf90_inquire_dimension(ncid,dimid,dim_nam,levsi) - error=nf90_close(ncid) - print*,'NetCDF file im,jm,lm= ',im,jm,levs,levsi - - else - call nemsio_init(iret=irets) - print *,'nemsio_init, iret=',irets - call nemsio_open(gfile,trim(fnsig),'read',iret=irets) - if ( irets /= 0 ) then - print*,"fail to open nems atmos file";stop - endif - - call nemsio_getfilehead(gfile,iret=irets - & ,dimx=im,dimy=jm,dimz=levsi) - if( irets /= 0 ) then - print*,'error finding model dimensions '; stop - endif - print*,'nemsio file im,jm,lm= ',im,jm,levsi - call nemsio_close(gfile,iret=irets) - endif - allocate (iocomms(0:ntot)) - if (fformat == 'netcdf') then - print*,'iocomms= ', iocomms - call mpi_comm_split(MPI_COMM_WORLD,ntask,0,iocomms(ntask),ierr) - call mpi_comm_rank(iocomms(ntask), iope, ierr) - call mpi_comm_size(iocomms(ntask), ionproc, ierr) - - call meteorg(npoint,rlat,rlon,istat,cstat,elevstn, - & nf,nfile,fnsig,jdate,idate, - & levsi,im,jm,nsfc, - & landwater,nend1, nint1, nint3, iidum,jjdum,np1, - & fformat,iocomms(ntask),iope,ionproc) - call mpi_barrier(iocomms(ntask), ierr) - call mpi_comm_free(iocomms(ntask), ierr) - else -!! For nemsio input - call meteorg(npoint,rlat,rlon,istat,cstat,elevstn, - & nf,nfile,fnsig,jdate,idate, - & levs,im,jm,nsfc, - & landwater,nend1, nint1, nint3, iidum,jjdum,np1, - & fformat,iocomms(ntask),iope,ionproc) - endif - endif - call mpi_barrier(mpi_comm_world,ierr) - call mpi_finalize(ierr) - if(mrank.eq.0) then - print *, ' starting to make bufr files' - print *, ' makebufr= ', makebufr - print *, 'nint1,nend1,nint3,nend= ',nint1,nend1,nint3,nend -!! idate = 0 7 1 2019 -!! jdate = 2019070100 - - if(makebufr) then - nend3 = nend - call buff(nint1,nend1,nint3,nend3, - & npoint,idate,jdate,levso, - & dird,lss,istat,sbset,seqflg,clist,npp,wrkd) - CALL W3TAGE('METEOMRF') - endif - endif - end diff --git a/src/gfs_bufr.fd/meteorg.f_org b/src/gfs_bufr.fd/meteorg.f_org deleted file mode 100644 index d4fa338d..00000000 --- a/src/gfs_bufr.fd/meteorg.f_org +++ /dev/null @@ -1,1326 +0,0 @@ - subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, - & nf,nfile,fnsig,jdate,idate, - & levs,im,jm,kdim, - & landwater,nend1,nint1,nint3,iidum,jjdum,np1, - & fformat,iocomms,iope,ionproc) - -!$$$ SUBPROGRAM DOCUMENTATION BLOCK -! . . . . -! 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 -! PRESSURE (line 290) -! CORRECT THE TEMPERATURE ADJUSTMENT (line 238) -! 2014-03-27 DANA CARLIS UNIFY CODE WITH GFS FORECAST MODEL PRECIP -! TYPE CALCULATION -! 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 -! AND RECALCULATE HOURLY AND 3 HOURLY OUTPUT DEPENDING -! 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 -! 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: -! npoint - number of points -! rlat(npint) - latitude -! rlon(npoint) - longtitude -! istat(npoint) - station id -! elevstn(npoint) - station elevation (m) -! nf - forecast cycle -! fnsig - sigma file name -! idate(4) - date -! levs - input vertical layers -! kdim - sfc file dimension -! -! OUTPUT: -! nfile - output data file channel -! jdate - date YYYYMMDDHH -! -! ATTRIBUTES: -! LANGUAGE: -! MACHINE: IBM SP -! -!$$$ - use netcdf - use nemsio_module - use sigio_module - use physcons - use mersenne_twister -! use funcphys, only : gfuncphys - use funcphys - implicit none - include 'mpif.h' - type(nemsio_gfile) :: gfile - type(nemsio_gfile) :: ffile - type(nemsio_gfile) :: ffile2 - integer :: nfile,npoint,levs,kdim - integer :: nfile1 - integer :: i,j,im,jm,kk,idum,jdum,idvc,idsl -! idsl Integer(sigio_intkind) semi-lagrangian id -! idvc Integer(sigio_intkind) vertical coordinate id -! (=1 for sigma, =2 for ec-hybrid, =3 for ncep hybrid) - integer,parameter :: nvcoord=2 - integer,parameter :: levso=64 - integer :: idate(4),nij,nflx2,np,k,l,nf,nfhour,np1 - integer :: idate_nems(7) - integer :: iret,jdate,leveta,lm,lp1 - character*512 :: fnsig,fngrib -!! real*8 :: data(6*levs+25) - real*8 :: data2(6*levso+25) - real*8 :: rstat1 - character*8 :: cstat1 - character*4 :: cstat(npoint) - real :: fhour,pp,ppn,qs,qsn,esn,es,psfc,ppi,dtemp,nd - real :: t,q,u,v,td,tlcl,plcl,qw,tw,xlat,xlon - integer,dimension(npoint):: landwater - integer,dimension(im,jm):: lwmask - real,dimension(im,jm):: apcp, cpcp - real,dimension(npoint,2+levs*3):: grids - real,dimension(npoint) :: rlat,rlon,pmsl,ps,psn,elevstn - real,dimension(im*jm) :: dum1d,dum1d2 - real,dimension(im,jm) :: gdlat, hgt, gdlon - real,dimension(im,jm,15) :: dum2d - real,dimension(im,jm,levs) :: t3d, q3d, uh, vh,omega3d - real,dimension(im,jm,levs) :: delpz - real,dimension(im,jm,levs+1) :: pint, zint - real,dimension(npoint,levs) :: gridu,gridv,omega,qnew,zp - real,dimension(npoint,levs) :: p1,pd3,ttnew - real,dimension(npoint,levs) :: z1 - real,dimension(npoint,levs+1) :: pi3 - real :: zp2(2) - real,dimension(kdim,npoint) :: sfc - real,dimension(1,levs+1) :: prsi,phii - real,dimension(1,levs) :: gt0,gq0,prsl,phy_f3d - real :: PREC,TSKIN,SR,randomno(1,2) - real :: DOMR,DOMZR,DOMIP,DOMS - real :: vcoord(levs+1,nvcoord),vdummy(levs+1) - real :: vcoordnems(levs+1,3,2) - real :: rdum - integer :: n3dfercld,iseedl - integer :: istat(npoint) - logical :: trace -!! logical, parameter :: debugprint=.true. - logical, parameter :: debugprint=.false. - character lprecip_accu*3 - real, parameter :: ERAD=6.371E6 - real, parameter :: DTR=3.1415926/180. - real :: ap - integer :: nf1, fint - integer :: nend1, nint1, nint3 - character*150 :: fngrib2 - integer recn_dpres,recn_delz,recn_dzdt - integer :: jrec - equivalence (cstat1,rstat1) - integer iidum(npoint),jjdum(npoint) - integer :: error, ncid, ncid2, id_var,dimid - character(len=100) :: long_name - character(len=6) :: fformat - integer,dimension(8) :: clocking - character(10) :: date - character(12) :: time - character(7) :: zone - character(3) :: Zreverse - character(20) :: VarName,LayName - integer iocomms,iope,ionproc - - nij = 12 -!! nflx = 6 * levs - nflx2 = 6 * levso - recn_dpres = 0 - recn_delz = 0 - recn_dzdt = 0 - jrec = 0 - lprecip_accu='yes' - - idvc=2 - idsl=1 -!read in NetCDF file header info - print*,"fformat= ", fformat -! print*,'meteorg.f, idum,jdum= ' -! do np = 1, npoint -! print*, iidum(np), jjdum(np) -! enddo - - if(fformat .eq. "netcdf") then - print*,'iocomms inside meteorg.f=', iocomms - error=nf90_open(trim(fnsig),ior(nf90_nowrite,nf90_mpiio), - & ncid,comm=iocomms, info = mpi_info_null) - error=nf90_get_att(ncid,nf90_global,"ak",vdummy) - do k = 1, levs+1 - vcoord(k,1)=vdummy(levs-k+1) - enddo - error=nf90_get_att(ncid,nf90_global,"bk",vdummy) - do k = 1, levs+1 - vcoord(k,2)=vdummy(levs-k+1) - enddo - error=nf90_inq_varid(ncid, "time", id_var) - error=nf90_get_var(ncid, id_var, nfhour) - print*, "nfhour:",nfhour - error=nf90_get_att(ncid,id_var,"units",long_name) -!! print*,'time units',' -- ',trim(long_name) - read(long_name(13:16),"(i4)")idate(4) - read(long_name(18:19),"(i2)")idate(2) - read(long_name(21:22),"(i2)")idate(3) - read(long_name(24:25),"(i2)")idate(1) - fhour=float(nfhour) - print*,'date= ', idate - jdate = idate(4)*1000000 + idate(2)*10000+ - & idate(3)*100 + idate(1) - print *, 'jdate = ', jdate - error=nf90_inq_varid(ncid, "lon", id_var) - error=nf90_get_var(ncid, id_var, gdlon) - 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 - - call nemsio_open(gfile,trim(fnsig),'read',iret=iret) - call nemsio_getfilehead(gfile,iret=iret - + ,idate=idate_nems(1:7),nfhour=nfhour - + ,idvc=idvc,idsl=idsl,lat=dum1d,lon=dum1d2 - + ,vcoord=vcoordnems) - - do k=1,levs+1 - vcoord(k,1)=vcoordnems(k,1,1) - vcoord(k,2)=vcoordnems(k,2,1) - end do - idate(1)=idate_nems(4) - idate(2)=idate_nems(2) - idate(3)=idate_nems(3) - idate(4)=idate_nems(1) - fhour=float(nfhour) - print *, ' processing forecast hour ', fhour - print *, ' idate =', idate - jdate = idate(4)*1000000 + idate(2)*10000+ - & idate(3)*100 + idate(1) - print *, 'jdate = ', jdate - print *, 'Total number of stations = ', npoint - ap = 0.0 - do j=1,jm - do i=1,im - gdlat(i,j)=dum1d((j-1)*im+i) - 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 - print*,'vcoord(k,2)= ', k, vcoord(k,2) - end do - print*,'sample lat= ',gdlat(im/5,jm/4) - + ,gdlat(im/5,jm/3),gdlat(im/5,jm/2) - print*,'sample lon= ',gdlon(im/5,jm/4) - + ,gdlon(im/5,jm/3),gdlon(im/5,jm/2) - endif -! topography - if (fformat == 'netcdf') then - VarName='hgtsfc' - Zreverse='yes' - call read_netcdf_p(ncid,im,jm,1,VarName,hgt,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'surface hgt not found' - else - VarName='hgt' - LayName='sfc' - call read_nemsio(gfile,im,jm,1,VarName,LayName,hgt, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'surface hgt not found' - endif - if(debugprint)print*,'sample sfc h= ',hgt(im/5,jm/4) - + ,hgt(im/5,jm/3),hgt(im/5,jm/2) - -! surface pressure (Pa) - if (fformat == 'netcdf') then - VarName='pressfc' - Zreverse='yes' - call read_netcdf_p(ncid,im,jm,1,VarName,pint(:,:,1), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'surface pressure not found' - else - VarName='pres' - LayName='sfc' - call read_nemsio(gfile,im,jm,1,VarName, - & LayName,pint(:,:,1),error) - if (error /= 0) print*,'surface pressure not found' - endif - if(debugprint)print*,'sample sfc P= ',pint(im/2,jm/4,1), - + pint(im/2,jm/3,1),pint(im/2,jm/2,1) - -! temperature using NetCDF - if (fformat == 'netcdf') then - VarName='tmp' - Zreverse='yes' - call read_netcdf_p(ncid,im,jm,levs,VarName,t3d,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'temp not found' - else - VarName='tmp' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName,t3d,error) - if (error /= 0) print*,'temp not found' - endif - if(debugprint) then - print*,'sample T at lev=1 to levs ' - do k = 1, levs - print*,k, t3d(im/2,jm/3,k) - enddo - endif -! specific humidity - if (fformat == 'netcdf') then - VarName='spfh' - Zreverse='yes' - call read_netcdf_p(ncid,im,jm,levs,VarName,q3d,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'spfh not found' - else - VarName='spfh' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName,q3d,error) - if (error /= 0) print*,'spfh not found' - endif - if(debugprint) then - print*,'sample Q at lev=1 to levs ' - do k = 1, levs - print*,k, q3d(im/2,jm/3,k) - enddo - endif -! U wind - if (fformat == 'netcdf') then - VarName='ugrd' - Zreverse='yes' - call read_netcdf_p(ncid,im,jm,levs,VarName,uh,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'ugrd not found' - else - VarName='ugrd' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName,uh,error) - if (error /= 0) print*,'ugrd not found' - endif - if(debugprint) then - print*,'sample U at lev=1 to levs ' - do k = 1, levs - print*,k, uh(im/2,jm/3,k) - enddo - endif -! V wind - if (fformat == 'netcdf') then - VarName='vgrd' - Zreverse='yes' - call read_netcdf_p(ncid,im,jm,levs,VarName,vh,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'vgrd not found' - else - VarName='vgrd' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName,vh,error) - if (error /= 0) print*,'vgrd not found' - endif - if(debugprint) then - print*,'sample V at lev=1 to levs ' - do k = 1, levs - print*,k, vh(im/2,jm/3,k) - enddo - endif -! dzdt !added by Guang Ping Lou for FV3GFS - if (fformat == 'netcdf') then - VarName='dzdt' - Zreverse='yes' - call read_netcdf_p(ncid,im,jm,levs,VarName,omega3d,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'dzdt not found' - else - VarName='dzdt' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName, - & omega3d,error) - if (error /= 0) print*,'dzdt not found' - endif - if(debugprint) then - print*,'sample dzdt at lev=1 to levs ' - do k = 1, levs - print*,k, omega3d(im/2,jm/3,k) - enddo - endif -! dpres !added by Guang Ping Lou for FV3GFS (interface pressure delta) - if (fformat == 'netcdf') then - VarName='dpres' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'dpres not found' - else - VarName='dpres' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName, - & delpz,error) - if (error /= 0) print*,'dpres not found' - endif - if(debugprint) then - print*,'sample delp at lev=1 to levs ' - do k = 1, levs - print*,k, delpz(im/2,jm/3,k) - enddo - endif -! compute interface pressure - if(recn_dpres == -9999) then - do k=2,levs+1 - do j=1,jm - do i=1,im - pint(i,j,k)=vcoord(k,1) - + +vcoord(k,2)*pint(i,j,1) - end do - end do - end do - else -! compute pint using dpres from top down if DZDT is used - if (fformat == 'netcdf') then - do j=1,jm - do i=1,im - pint(i,j,levs+1) = delpz(i,j,1) - end do - end do - do k=levs,2,-1 - kk=levs-k+2 - do j=1,jm - do i=1,im - pint(i,j,k) = pint(i,j,k+1) + delpz(i,j,kk) - end do - end do - end do - else - do k=2,levs+1 - do j=1,jm - do i=1,im - pint(i,j,k) = pint(i,j,k-1) - delpz(i,j,k-1) - end do - end do - end do - endif - if(debugprint) then - print*,'sample interface pressure pint at lev =1 to levs ' - do k = 1, levs+1 - print*,k, pint(im/2,jm/3,k),pint(im/3,jm/8,k) - enddo - endif - endif -! delz !added by Guang Ping Lou for FV3GFS ("height thickness" with unit "meters" bottom up) - if (fformat == 'netcdf') then - VarName='delz' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,levs,VarName,delpz,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'delz not found' - else - VarName='delz' - LayName='mid layer' - call read_nemsio(gfile,im,jm,levs,VarName,LayName,delpz,error) - if (error /= 0) print*,'delz not found' - endif - if(debugprint) then - print*,'sample delz at lev=1 to levs ' - do k = 1, levs - print*,k, delpz(im/2,jm/3,k) - enddo - endif - -! compute interface height (meter) - if(recn_delz == -9999) then - print*, 'using calculated height' - else -! compute zint using delz from bot up if DZDT is used - if (fformat == 'netcdf') then - do j=1,jm - do i=1,im - zint(i,j,1) = 0.0 - end do - end do - do k=2,levs+1 - kk=levs-k+1 - do j=1,jm - do i=1,im - zint(i,j,k) = zint(i,j,k-1) - delpz(i,j,kk) - end do - end do - end do - else - do k=2,levs+1 - do j=1,jm - do i=1,im - zint(i,j,k) = zint(i,j,k-1) + delpz(i,j,k-1) - end do - end do - end do - endif - if(debugprint) then - print*,'sample interface height zint at lev =1 to levs ' - do k = 1, levs+1 - print*,k, zint(im/2,jm/3,k),zint(im/3,jm/8,k) - enddo - endif - endif - -! close up this NetCDF file - error=nf90_close(ncid) - -! Now open up NetCDF surface files - if ( nf .le. nend1 ) then - nf1 = nf - nint1 - else - nf1 = nf - nint3 - endif - if ( nf == 0 ) nf1=0 - if(nf==0) then - fngrib='flxf00' - elseif(nf.lt.10) then - fngrib='flxf0' - write(fngrib(6:6),'(i1)') nf - elseif(nf.lt.100) then - fngrib='flxf' - write(fngrib(5:6),'(i2)') nf - else - fngrib='flxf' - write(fngrib(5:7),'(i3)') nf - endif - if(nf1==0) then - fngrib2='flxf00' - elseif(nf1.lt.10) then - fngrib2='flxf0' - write(fngrib2(6:6),'(i1)') nf1 - elseif(nf1.lt.100) then - fngrib2='flxf' - write(fngrib2(5:6),'(i2)') nf1 - else - fngrib2='flxf' - write(fngrib2(5:7),'(i3)') nf1 - endif - if (fformat == 'netcdf') then - error=nf90_open(trim(fngrib),nf90_nowrite,ncid) -!open T-nint below - error=nf90_open(trim(fngrib2),nf90_nowrite,ncid2) - if(error /= 0)print*,'file not open',trim(fngrib), trim(fngrib2) - 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) - endif -! land water mask - if (fformat == 'netcdf') then - VarName='land' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,lwmask,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'lwmask not found' - else - VarName='land' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName,lwmask,error) - if (error /= 0) print*,'lwmask not found' - endif - if(debugprint) - + print*,'sample land mask= ',lwmask(im/2,jm/4), - + lwmask(im/2,jm/3) - -! surface T - if (fformat == 'netcdf') then - VarName='tmpsfc' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,1), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'tmpsfc not found' - else - VarName='tmp' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - & dum2d(:,:,1),error) - if (error /= 0) print*,'tmpsfc not found' - endif - if(debugprint) - + print*,'sample sfc T= ',dum2d(im/2,jm/4,1),dum2d(im/2,jm/3,1), - + dum2d(im/2,jm/2,1) -! 2m T - if (fformat == 'netcdf') then - VarName='tmp2m' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,2), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'tmp2m not found' - else - VarName='tmp' - LayName='2 m above gnd' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,2),error) - if (error /= 0) print*,'tmp2m not found' - endif - if(debugprint) - + print*,'sample 2m T= ',dum2d(im/2,jm/4,2),dum2d(im/2,jm/3,2), - + dum2d(im/2,jm/2,2) - -! 2m Q - if (fformat == 'netcdf') then - VarName='spfh2m' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,3), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'spfh2m not found' - else - VarName='spfh' - LayName='2 m above gnd' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,3),error) - if (error /= 0) print*,'spfh2m not found' - endif - if(debugprint) - + print*,'sample 2m Q= ',dum2d(im/2,jm/4,3),dum2d(im/2,jm/3,3), - + dum2d(im/2,jm/2,3) - -! U10 - if (fformat == 'netcdf') then - VarName='ugrd10m' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,4), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'ugrd10m not found' - else - VarName='ugrd' - LayName='10 m above gnd' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,4),error) - if (error /= 0) print*,'ugrd10m not found' - endif - -! V10 - if (fformat == 'netcdf') then - VarName='vgrd10m' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,5), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'vgrd10m not found' - else - VarName='vgrd' - LayName='10 m above gnd' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,5),error) - if (error /= 0) print*,'vgrd10m not found' - endif - -! soil T - if (fformat == 'netcdf') then - VarName='soilt1' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,6), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'soilt1 not found' - else - VarName='tmp' - LayName='0-10 cm down' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,6),error) - if (error /= 0) print*,'soil T not found' - endif - if(debugprint) - + print*,'sample soil T= ',dum2d(im/2,jm/4,6),dum2d(im/2,jm/3,6), - + dum2d(im/2,jm/2,6) - -! snow depth - if (fformat == 'netcdf') then - VarName='snod' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,7), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'snod not found' - else - VarName='snod' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,7),error) - if (error /= 0) print*,'snod not found' - endif - -! evaporation -!instantaneous surface latent heat net flux - if (fformat == 'netcdf') then - VarName='lhtfl' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,8), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'lhtfl not found' - else - VarName='lhtfl' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,8),error) - if (error /= 0) print*,'lhtfl not found' - endif - if(debugprint) - + print*,'evaporation latent heat net flux= ', - + dum2d(im/2,jm/4,8),dum2d(im/2,jm/3,8) - if(debugprint) - + print*,'evaporation latent heat net flux stn 000692)= ', - + dum2d(2239,441,8) - -! total precip - if ( nf .le. nend1 ) then - fint = nint1 - else - fint = nint3 - endif -! for accumulated precipitation: - if (fformat == 'netcdf') then - VarName='prate_ave' - Zreverse='no' -!! call read_netcdf_p(ncid,im,jm,1,VarName,apcp,Zreverse,error) !current hour - call read_netcdf_p(ncid,im,jm,1,VarName,apcp,Zreverse, - & iope,ionproc,iocomms,error) -!! call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse,error) !earlier hour - call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'prate_ave not found' - else - VarName='prate_ave' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + apcp,error) - call read_nemsio(ffile2,im,jm,1,VarName,LayName, - + cpcp,error) - if (error /= 0) print*,'prate_ave2 not found' - endif - if(debugprint) - & print*,'sample fhour ,3= ', fhour, - & '1sample precip rate= ',apcp(im/2,jm/3),cpcp(im/2,jm/3) - ap=fhour-fint - do j=1,jm - do i=1,im - dum2d(i,j,9) =(apcp(i,j)*fhour-cpcp(i,j)*ap)*3600.0 - end do - end do - - if(debugprint) - & print*,'sample fhour ,5= ', fhour, - & 'sample total precip= ',dum2d(im/2,jm/4,9), - + dum2d(im/2,jm/3,9),dum2d(im/2,jm/2,9) - -! convective precip - if (fformat == 'netcdf') then - VarName='cprat_ave' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,apcp,Zreverse, - & iope,ionproc,iocomms,error) - call read_netcdf_p(ncid2,im,jm,1,VarName,cpcp,Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'cprat_ave not found' - else - VarName='cprat_ave' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + apcp,error) - call read_nemsio(ffile2,im,jm,1,VarName,LayName, - + cpcp,error) - if (error /= 0) print*,'cprat_ave2 not found' - endif - ap=fhour-fint - 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 - -! water equi - if (fformat == 'netcdf') then - VarName='weasd' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName,dum2d(:,:,11), - & Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'weasd not found' - else - VarName='weasd' - LayName='sfc' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,11),error) - if (error /= 0) print*,'weasd not found' - endif - -! low cloud fraction - if (fformat == 'netcdf') then - VarName='tcdc_avelcl' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName, - & dum2d(:,:,12),Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'tcdc_avelcl not found' - else - VarName='tcdc_ave' - LayName='low cld lay' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,12),error) - if (error /= 0) print*,'low cld lay not found' - endif - -! mid cloud fraction - if (fformat == 'netcdf') then - VarName='tcdc_avemcl' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName, - & dum2d(:,:,13),Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'tcdc_avemcl not found' - else - VarName='tcdc_ave' - LayName='mid cld lay' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,13),error) - if (error /= 0) print*,'mid cld lay not found' - endif - -! high cloud fraction - if (fformat == 'netcdf') then - VarName='tcdc_avehcl' - Zreverse='no' - call read_netcdf_p(ncid,im,jm,1,VarName, - & dum2d(:,:,14),Zreverse, - & iope,ionproc,iocomms,error) - if (error /= 0) print*,'tcdc_avehcl not found' - else - VarName='tcdc_ave' - LayName='high cld lay' - call read_nemsio(ffile,im,jm,1,VarName,LayName, - + dum2d(:,:,14),error) - if (error /= 0) print*,'high cld lay not found' - endif - - if(debugprint) - + print*,'sample high cloud frac= ',dum2d(im/2,jm/4,14), - + dum2d(im/2,jm/3,14),dum2d(im/2,jm/2,14) - - if (fformat == 'netcdf') then - error=nf90_close(ncid) - error=nf90_close(ncid2) - else - call nemsio_close(ffile,iret=error) - call nemsio_close(ffile2,iret=error) - endif - call date_and_time(date,time,zone,clocking) -! print *,'10reading surface data end= ', clocking - print *,'10date, time, zone',date, time, zone -! -! get the nearest neighbor i,j from the table -! - do np=1, npoint -! use read in predetermined i,j - if (np1==0) then - idum=iidum(np) - jdum=jjdum(np) - - else -! find nearest neighbor - rdum=rlon(np) - if(rdum<0.)rdum=rdum+360. - - do j=1,jm-1 - do i=1,im-1 - if((rdum>=gdlon(i,j) .and. rdum<=gdlon(i+1,j)) .and. - + (rlat(np)<=gdlat(i,j).and.rlat(np)>=gdlat(i,j+1)) ) then - if(landwater(np) == 2)then - idum=i - jdum=j - exit - else if(landwater(np) == lwmask(i,j))then - idum=i - jdum=j !1 - exit - else if(landwater(np) == lwmask(i+1,j))then - idum=i+1 - jdum=j ! 2 - exit - else if(landwater(np) == lwmask(i-1,j))then - idum=i-1 - jdum=j ! 3 - exit - else if(landwater(np) == lwmask(i,j+1))then - idum=i - jdum=j+1 ! 4 - exit - else if(landwater(np) == lwmask(i,j-1))then - idum=i - jdum=j-1 ! 5 - exit - else if(landwater(np) == lwmask(i+1,j-1))then - idum=i+1 - jdum=j-1 ! 6 - exit - else if(landwater(np) == lwmask(i+1,j+1))then - idum=i+1 - jdum=j+1 ! 7 - exit - else if(landwater(np) == lwmask(i-1,j+1))then - idum=i-1 - jdum=j+1 ! 8 - exit - else if(landwater(np) == lwmask(i-1,j-1))then - idum=i-1 - jdum=j-1 ! 9 - exit - else if(landwater(np) == lwmask(i,j+2))then - idum=i - jdum=j+2 ! 10 - exit - else if(landwater(np) == lwmask(i+2,j))then - idum=i+2 - jdum=j !11 - exit - else if(landwater(np) == lwmask(i,j-2))then - idum=i - jdum=j-2 ! 12 - exit - else if(landwater(np) == lwmask(i-2,j))then - idum=i-2 - jdum=j !13 - exit - else if(landwater(np) == lwmask(i-2,j+1))then - idum=i-2 - jdum=j+1 ! 14 - exit - else if(landwater(np) == lwmask(i-1,j+2))then - idum=i-1 - jdum=j+2 !15 - exit - else if(landwater(np) == lwmask(i+1,j+2))then - idum=i+1 - jdum=j+2 !16 - exit - else if(landwater(np) == lwmask(i+2,j+1))then - idum=i+2 - jdum=j+1 !17 - exit - else if(landwater(np) == lwmask(i+2,j-1))then - idum=i+2 - jdum=j-1 !18 - exit - else if(landwater(np) == lwmask(i+1,j-2))then - idum=i+1 - jdum=j-2 !19 - exit - else if(landwater(np) == lwmask(i-1,j-2))then - idum=i-1 - jdum=j-2 !20 - exit - else if(landwater(np) == lwmask(i-2,j-1))then - idum=i-2 - jdum=j-1 !21 - exit - else if(landwater(np) == lwmask(i-2,j-2))then - idum=i-2 - jdum=j-2 !22 - exit - else if(landwater(np) == lwmask(i+2,j-2))then - idum=i+2 - jdum=j-2 !23 - exit - else if(landwater(np) == lwmask(i+2,j+2))then - idum=i+2 - jdum=j+2 !24 - exit - else if(landwater(np) == lwmask(i-2,j+2))then - idum=i-2 - jdum=j+2 !25 - exit - else if(landwater(np) == lwmask(i+3,j))then - idum=i+3 - jdum=j !26 - exit - else if(landwater(np) == lwmask(i-3,j))then - idum=i-3 - jdum=j !27 - exit - else if(landwater(np) == lwmask(i,j+3))then - idum=i - jdum=j+3 !28 - exit - else if(landwater(np) == lwmask(i,j-3))then - idum=i - jdum=j-3 !29 - exit - else -CC print*,'no matching land sea mask np,landwater,i,j,mask= ' -CC print*, np,landwater(np),i,j,lwmask(i,j) -CC print*, ' So it takes i,j ' - idum=i - jdum=j - exit - end if - end if - end do - end do - - idum=max0(min0(idum,im),1) - jdum=max0(min0(jdum,jm),1) - endif !! read in i,j ends here - if (fhour==0.0) then - if(debugprint) then - write(nij,98) np,idum,jdum,rlat(np),rlon(np) - 98 FORMAT (3I6, 2F9.2) - if(elevstn(np)==-999.) elevstn(np)=hgt(idum,jdum) - write(9,99) np,rlat(np),rlon(np),elevstn(np),hgt(idum,jdum) - 99 FORMAT (I6, 4F9.2) - if(np==1 .or.np==100)print*,'nearest neighbor for station ',np - + ,idum,jdum,rlon(np),rlat(np),lwmask(i,j),landwater(np) - endif - endif - - 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) - sfc(12,np)=dum2d(idum,jdum,9) - sfc(11,np)=dum2d(idum,jdum,10) - sfc(10,np)=dum2d(idum,jdum,11) - sfc(27,np)=dum2d(idum,jdum,12) - sfc(26,np)=dum2d(idum,jdum,13) - sfc(25,np)=dum2d(idum,jdum,14) - sfc(34,np)=dum2d(idum,jdum,4) - sfc(35,np)=dum2d(idum,jdum,5) - sfc(30,np)=dum2d(idum,jdum,2) - sfc(31,np)=dum2d(idum,jdum,3) - -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) - - do k=1,levs - 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) - gridv(np,k)=vh(idum,jdum,k) - p1(np,k)=pint(idum,jdum,k+1) - 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 -! !ps in kPa - ps(np) = grids(np,2)/1000. !! surface pressure - enddo - -! -! ----------------- -! Put topo(1),surf press(2),vir temp(3:66),and specifi hum(67:130) in grids -! for each station -!! if(recn_dzdt == 0 ) then !!DZDT - do k = 1, levs - do np = 1, npoint - omega(np,k) = grids(np,2+levs*2+k) - enddo - enddo - if(debugprint) - + print*,'sample (omega) dzdt ', (omega(3,k),k=1,levs) -! -! move surface pressure to the station surface from the model surface -! - do np = 1, npoint -! -! when the station elevation information in the table says missing, -! use the model elevation -! -! print *, "elevstn = ", elevstn(np) - if(elevstn(np)==-999.) elevstn(np) = grids(np,1) - psn(np) = ps(np) - call sigio_modpr(1,1,levs,nvcoord,idvc, - & idsl,vcoord,iret, - & 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) - if(np==11)print*,'pd3= ', pd3(np,1:levs) - enddo -! -!! test removing height adjustments - print*, 'do not do height adjustments' -! -! get sea-level pressure (Pa) and layer geopotential height -! - do k = 1, levs - do np = 1, npoint - ttnew(np,k) = grids(np,k+2) - qnew(np,k) = grids(np,k+levs+2) - enddo - enddo - - do np=1,npoint -!! call gslp(levs,elevstn(np),ps(np)*1000, - call gslp(levs,grids(np,1),ps(np)*1000, - & p1(np,1:levs),ttnew(np,1:levs),qnew(np,1:levs), - & pmsl(np),zp(np,1:levs),zp2(1:2)) - enddo - print *, 'call gslp pmsl= ', (pmsl(np),np=1,20) - if(recn_delz == -9999) then - print*, 'using calculated height ' - else - print*, 'using model height m' - do k = 1, levs - do np=1, npoint - zp(np,k) = z1(np,k) - enddo - enddo - 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) -! -! prepare buffer data -! - if(iope == 0) then - call gfuncphys - do np = 1, npoint - pi3(np,1)=psn(np)*1000 - do k=1,levs - pi3(np,k+1)=pi3(np,k)-pd3(np,k) !layer pressure (Pa) - enddo -!! ==ivalence (cstat1,rstat1) - cstat1=cstat(np) -!! data(1) = ifix(fhour+.2) * 3600 ! FORECAST TIME (SEC) -!! data(2) = istat(np) ! STATION NUMBER -!! data(3) = rstat1 ! STATION CHARACTER ID -!! data(4) = rlat(np) ! LATITUDE (DEG N) -!! data(5) = rlon(np) ! LONGITUDE (DEG E) -!! data(6) = elevstn(np) ! STATION ELEVATION (M) - data2(1) = ifix(fhour+.2) * 3600 ! FORECAST TIME (SEC) - data2(2) = istat(np) ! STATION NUMBER - data2(3) = rstat1 ! STATION CHARACTER ID - data2(4) = rlat(np) ! LATITUDE (DEG N) - data2(5) = rlon(np) ! LONGITUDE (DEG E) - data2(6) = elevstn(np) ! STATION ELEVATION (M) - psfc = 10. * psn(np) ! convert to MB - leveta = 1 - do k = 1, levs - kk= k/2 + 1 -! -! look for the layer above 500 mb for precip type computation -! - if(pi3(np,k).ge.50000.) leveta = k - ppi = pi3(np,k) - t = grids(np,k+2) - q = max(1.e-8,grids(np,2+k+levs)) - u = gridu(np,k) - v = gridv(np,k) -!! data((k-1)*6+7) = p1(np,k) ! PRESSURE (PA) at integer layer -!! data((k-1)*6+8) = t ! TEMPERATURE (K) -!! data((k-1)*6+9) = u ! U WIND (M/S) -!! data((k-1)*6+10) = v ! V WIND (M/S) -!! data((k-1)*6+11) = q ! HUMIDITY (KG/KG) -!! data((k-1)*6+12) = omega(np,k)*100. ! Omega (pa/sec) !changed to dzdt(cm/s) if available - 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+11) = q - data2((kk-1)*6+12) = omega(np,k)*100. - endif -!changed to dzdt(cm/s) if available - enddo -! -! process surface flux file fields -! -!! data(8+nflx) = psfc * 100. ! SURFACE PRESSURE (PA) -!! data(7+nflx) = pmsl(np) - data2(8+nflx2) = psfc * 100. ! SURFACE PRESSURE (PA) - 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(37,np) = data2(7+nflx2) * .01 - 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: -! It was decided that no corrctions were needed due to higher model -! resolution. -! -! if(sfc(10,np)==0.) then -! sfc(30,np) = sfc(30,np) + dtemp -! sfc(5,np) = sfc(5,np) + dtemp -! endif -! if(sfc(10,np).gt.0..and.sfc(5,np).lt.273.16) then -! sfc(5,np) = sfc(5,np) + dtemp -! if(sfc(5,np).gt.273.16) then -! dtemp = sfc(5,np) - 273.16 -! sfc(5,np) = 273.16 -! endif -! sfc(30,np) = sfc(30,np) + dtemp -! endif -! -!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 -! 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 - if(debugprint) - + print*,'evaporation (stn 000692)= ',sfc(17,np) -!! data(9+nflx) = sfc(5,np) ! tsfc (K) -!! data(10+nflx) = sfc(6,np) ! 10cm soil temp (K) -!!! data(11+nflx) = sfc(17,np)/28.3565 ! evaporation (kg/m**2) from (W m-2) -!! data(11+nflx) = sfc(17,np)*0.035265 ! evaporation (kg/m**2) from (W m-2) -!! data(12+nflx) = sfc(12,np) ! total precip (m) -!! data(13+nflx) = sfc(11,np) ! convective precip (m) -!! data(14+nflx) = sfc(10,np) ! water equi. snow (m) -!! data(15+nflx) = sfc(27,np) ! low cloud (%) -!! data(16+nflx) = sfc(26,np) ! mid cloud -!! data(17+nflx) = sfc(25,np) ! high cloud -!! data(18+nflx) = sfc(34,np) ! U10 (m/s) -!! data(19+nflx) = sfc(35,np) ! V10 (m/s) -!! data(20+nflx) = sfc(30,np) ! T2 (K) -!! data(21+nflx) = sfc(31,np) ! Q2 (K) - -!! data(22+nflx) = 0. -!! data(23+nflx) = 0. -!! data(24+nflx) = 0. -!! data(25+nflx) = 0. -!! create 64 level bufr files - data2(9+nflx2) = sfc(5,np) ! tsfc (K) - data2(10+nflx2) = sfc(6,np) ! 10cm soil temp (K) -!! data2(11+nflx2) = sfc(17,np)/28.3565 ! evaporation (kg/m**2) from (W m-2) - data2(11+nflx2) = sfc(17,np)*0.035265 ! evaporation (kg/m**2) from (W m-2) - data2(12+nflx2) = sfc(12,np) ! total precip (m) - data2(13+nflx2) = sfc(11,np) ! convective precip (m) - data2(14+nflx2) = sfc(10,np) ! water equi. snow (m) - data2(15+nflx2) = sfc(27,np) ! low cloud (%) - data2(16+nflx2) = sfc(26,np) ! mid cloud - data2(17+nflx2) = sfc(25,np) ! high cloud - data2(18+nflx2) = sfc(34,np) ! U10 (m/s) - data2(19+nflx2) = sfc(35,np) ! V10 (m/s) - data2(20+nflx2) = sfc(30,np) ! T2 (K) - data2(21+nflx2) = sfc(31,np) ! Q2 (K) - - data2(22+nflx2) = 0. - data2(23+nflx2) = 0. - data2(24+nflx2) = 0. - data2(25+nflx2) = 0. - nd = 0 - trace = .false. - DOMS=0. - DOMR=0. - DOMIP=0. - DOMZR=0. - if(np==1.or.np==2) nd = 1 - if(np==1.or.np==2) trace = .true. - - 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) - t = grids(np,k+2) - q = max(0.,grids(np,2+k+levs)) - u = gridu(np,k) - v = gridv(np,k) - if(q.gt.1.e-6.and.pp.ge.20000.) then - call tdew(td,t,q,pp) - call lcl(tlcl,plcl,t,pp,q) - call mstadb(qw,tw,pp,q,tlcl,plcl) - else - td = t - 30. - tw = t - 30. - endif -! Calpreciptype input variables - gt0(1,k)= t - gq0(1,k) = q - prsl(1,k) = pp - prsi(1,k)=ppi - phii(1,k)=zp(np,k) !height in meters - enddo -! Use GFS routine calpreciptype.f to calculate precip type - xlat=rlat(np) - xlon=rlon(np) - lm=leveta - lp1=leveta+1 -!! PREC=data(12+nflx) - PREC=data2(12+nflx2) - n3dfercld=1 !if =3 then use Ferriers Explicit Precip Type - TSKIN=1. !used in Ferriers Explicit Precip Scheme - SR=1. !used in Ferriers Explicit Precip Scheme - iseedl=jdate - call random_setseed(iseedl) - call random_number(randomno) - call calpreciptype(1,1,1,1,lm,lp1,randomno,xlat,xlon, !input - & gt0,gq0,prsl,prsi,PREC,phii,n3dfercld,TSKIN,SR,phy_f3d, !input - & DOMR,DOMZR,DOMIP,DOMS) ! Output vars - endif -!! data(nflx + 22) = DOMS -!! data(nflx + 23) = DOMIP -!! data(nflx + 24) = DOMZR -!! data(nflx + 25) = DOMR - data2(nflx2 + 22) = DOMS - data2(nflx2 + 23) = DOMIP - data2(nflx2 + 24) = DOMZR - data2(nflx2 + 25) = DOMR - if(np==1.or.np==100) then - print *, ' surface fields for hour', nf, 'np =', np - print *, (data2(l+nflx2),l=1,25) - print *, ' temperature sounding' - print 6101, (data2((k-1)*6+8),k=1,levso) - print *, ' omega sounding' - print *, (data2((k-1)*6+12),k=1,levso) - endif -C print *, 'in meteorg nfile1= ', nfile1 -!! write(nfile) data - write(nfile) data2 - enddo !End loop over stations np - endif - call date_and_time(date,time,zone,clocking) -! print *,'13reading write data end= ', clocking - print *,'13date, time, zone',date, time, zone - print *, 'in meteorg nf,nfile,nfhour= ', nf,nfile,nfhour - print *, 'Finished writing bufr data file' - 6101 format(2x,6f12.3) - 6102 format(2x,6f12.5) - 6103 format(2x,6f12.5) -! - close(unit=nfile) - return - 910 print *, ' error reading surface flux file' - end - -!----------------------------------------------------------------------- From 979192f5ffe0308d6645f26914d26f2eb4d65a47 Mon Sep 17 00:00:00 2001 From: "Bo.Cui" Date: Wed, 21 Aug 2024 09:39:55 -0400 Subject: [PATCH 3/6] comment the line with code not used --- src/gfs_bufr.fd/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gfs_bufr.fd/CMakeLists.txt b/src/gfs_bufr.fd/CMakeLists.txt index 09c13e13..5328092c 100644 --- a/src/gfs_bufr.fd/CMakeLists.txt +++ b/src/gfs_bufr.fd/CMakeLists.txt @@ -11,7 +11,7 @@ list(APPEND fortran_src newsig1.f read_nemsio.f read_netcdf.f - read_netcdf_p.f + #read_netcdf_p.f rsearch.f svp.f tdew.f From ee6bf3b7e0b317f2004af7921a5d985fc89e501c Mon Sep 17 00:00:00 2001 From: "Bo.Cui" Date: Fri, 23 Aug 2024 15:00:50 -0400 Subject: [PATCH 4/6] Add module modpr_module.f90 --- src/gfs_bufr.fd/CMakeLists.txt | 4 +- src/gfs_bufr.fd/gfsbufr.f | 2 - src/gfs_bufr.fd/meteorg.f | 4 +- src/gfs_bufr.fd/modpr_module.f90 | 456 +++++++++++++++++++++++++++++++ src/gfs_bufr.fd/modstuff1.f | 3 +- 5 files changed, 462 insertions(+), 7 deletions(-) create mode 100644 src/gfs_bufr.fd/modpr_module.f90 diff --git a/src/gfs_bufr.fd/CMakeLists.txt b/src/gfs_bufr.fd/CMakeLists.txt index 5328092c..217a14d1 100644 --- a/src/gfs_bufr.fd/CMakeLists.txt +++ b/src/gfs_bufr.fd/CMakeLists.txt @@ -9,7 +9,6 @@ list(APPEND fortran_src meteorg.f mstadb.f newsig1.f - read_nemsio.f read_netcdf.f #read_netcdf_p.f rsearch.f @@ -17,6 +16,7 @@ list(APPEND fortran_src tdew.f terp3.f vintg.f + modpr_module.f90 ) list(APPEND fortran_src_free @@ -40,10 +40,8 @@ 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) diff --git a/src/gfs_bufr.fd/gfsbufr.f b/src/gfs_bufr.fd/gfsbufr.f index a3f1acbf..5d31625c 100644 --- a/src/gfs_bufr.fd/gfsbufr.f +++ b/src/gfs_bufr.fd/gfsbufr.f @@ -44,12 +44,10 @@ program meteormrf C C$$$ use netcdf - use sigio_module implicit none integer,parameter:: nsta=3000 integer,parameter:: ifile=11 integer,parameter:: levso=64 - integer(sigio_intkind):: irets integer ncfsig, nsig integer istat(nsta), idate(4), jdate, nfhour integer :: levs,nstart,nend,nint,nsfc,levsi,im,jm diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f index 852c125d..83ed46a0 100644 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -37,6 +37,8 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! adjustment ! 2023-03-28 Bo Cui Fix compilation error with "-check all" for gfs_bufrsnd ! 2024-08-08 Bo Cui UPDATE TO HANDLE ONE FORECAST AT A TIME, REMOVE NEMSIO INPUT FILES +! 2024-08-23 Bo Cui Replace sigio_module with the simplified module modpr_module +! ! ! USAGE: CALL PROGRAM meteorg ! INPUT: @@ -63,7 +65,7 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! !$$$ use netcdf - use sigio_module + use modpr_module use physcons use mersenne_twister ! use funcphys, only : gfuncphys diff --git a/src/gfs_bufr.fd/modpr_module.f90 b/src/gfs_bufr.fd/modpr_module.f90 new file mode 100644 index 00000000..70ffbec9 --- /dev/null +++ b/src/gfs_bufr.fd/modpr_module.f90 @@ -0,0 +1,456 @@ +!------------------------------------------------------------------------------- +module modpr_module +!$$$ Module Documentation Block +! +! Original module sourced from: /apps/ops/prod/libs/build/v1.1.0/pkg/sigio-v2.3.2/src/sigio_module.f +! +! Modifications: +! - Copied the original sigio_module.f to the folder gfs_bufr.fd and renamed it to modpr_module.f90 +! - Retained only the sigio_modpr subroutine and its related code +! - Removed all other subroutines and associated code that were present in sigio_module.f +! +! Reason for modification: +! - To simplify the module by including only the necessary subroutines +! +! Program History Log: +! 2024-08-23 Bo Cui: Simplified the module by retaining only sigio_modpr subroutine +! +! Below are Original Module Documentation Block from sigio_module.f +! +! Module: sigio_module API for global spectral sigma file I/O +! Prgmmr: iredell Org: w/nx23 date: 1999-01-18 +! +! Abstract: This module provides an Application Program Interface +! for performing I/O on the sigma restart file of the global spectral model. +! Functions include opening, reading, writing, and closing as well as +! allocating and deallocating data buffers used in the transfers. +! The I/O performed here is sequential. +! The transfers are limited to header records or data records. +! +! Program History Log: +! 1999-01-18 Mark Iredell +! 2013-10-14 Fanglin Yang: Added dynamics restart fields (ixga etc) +! and restructureed physics restart fields (ixgr etc). +! 2018-05-11 Mark Iredell: Added error check for NEMSIO file. +! +! Public Variables: +! sigio_lhead1 Integer parameter length of first header record (=32) +! sigio_charkind Integer parameter kind or length of passed characters (=8) +! sigio_intkind Integer parameter kind or length of passed integers (=4) +! sigio_realkind Integer parameter kind or length of passed reals (=4) +! sigio_dblekind Integer parameter kind or length of passed longreals (=8) +! sigio_realfill Real(sigio_realkind) parameter fill value (=-9999.) +! sigio_dblefill Real(sigio_dblekind) parameter fill value (=-9999.) +! +! Public Defined Types: +! sigio_head Sigma file header information +! clabsig Character(sigio_lhead1) ON85 label +! (obsolescent) +! fhour Real(sigio_realkind) forecast hour +! idate Integer(sigio_intkind)(4) initial date +! (hour, month, day, 4-digit year) +! si Real(sigio_realkind)(101) sigma interfaces +! (obsolescent) +! sl Real(sigio_realkind)(100) sigma levels +! (obsolescent) +! ak Real(sigio_realkind)(101) hybrid interface a +! (obsolescent) +! bk Real(sigio_realkind)(101) hybrid interface b +! (obsolescent) +! jcap Integer(sigio_intkind) spectral truncation +! levs Integer(sigio_intkind) number of levels +! itrun Integer(sigio_intkind) truncation flag +! (=1 for triangular) +! iorder Integer(sigio_intkind) coefficient order flag +! (=2 for ibm order) +! irealf Integer(sigio_intkind) floating point flag +! (=1 for 4-byte ieee, =2 for 8-byte ieee) +! igen Integer(sigio_intkind) model generating flag +! latf Integer(sigio_intkind) latitudes in dynamics +! (=(jcap+1)*3/2) +! lonf Integer(sigio_intkind) longitudes in dynamics +! (>=(jcap+1)*3 appropriate for fft) +! latb Integer(sigio_intkind) latitudes in physics +! lonb Integer(sigio_intkind) longitudes in physics +! latr Integer(sigio_intkind) latitudes in radiation +! lonr Integer(sigio_intkind) longitudes in radiation +! ntrac Integer(sigio_intkind) number of tracers +! icen2 Integer(sigio_intkind) subcenter id +! iens Integer(sigio_intkind)(2) ensemble ids +! idpp Integer(sigio_intkind) processing id +! idsl Integer(sigio_intkind) semi-lagrangian id +! idvc Integer(sigio_intkind) vertical coordinate id +! (=1 for sigma, =2 for ec-hybrid, =3 for ncep hybrid) +! idvm Integer(sigio_intkind) mass variable id +! idvt Integer(sigio_intkind) tracer variable id +! idrun Integer(sigio_intkind) run id +! idusr Integer(sigio_intkind) user-defined id +! pdryini Real(sigio_realkind) global mean dry air pressure (kPa) +! (obsolescent) +! ncldt Integer(sigio_intkind) number of cloud types +! ixgr Integer(sigio_intkind) extra fileds for physics. +! ixgr=00000000 no extra fields +! ixgr=0000000a zhao micro, a=1: zhao1, two 3d, one 2d, and nxss=0 +! a=2: zhao2, four 3d, three 2d, and nxss=0 +! a=3: zhao2, four 3d, three 2d, and nxss=1 +! ixgr=000000b0 ferrier micro, b=1: three 3d, one 2d, and nxss=0 +! ferrier micro, b=2: three 3d, one 2d, and nxss=1 +! ixgr=00000c00 c=1, pdf cld, three 3d +! ixga Integer(sigio_intkind) extra fileds for dynamics. +! ixga=00000000 no extra fields +! ixga=0000000a zflxtvd micro, ntrac 3d, zero 2d +! ixga=000000b0 (reserved for) joe-sela semi-lag gfs +! ivs Integer(sigio_intkind) version number +! nvcoord Integer(sigio_intkind) number of vcoord profiles +! The following variables should be allocated with sigio_alhead: +! vcoord Real(sigio_realkind)((levs+1),nvcoord) vcoord profiles +! cfvars Character(8)(5+ntrac) field variable names +! The following variables should not be modified by the user: +! nxgr Integer(sigio_intkind) number of extra physics grid fields +! nxss Integer(sigio_intkind) number of extra scalars +! nxga Integer(sigio_intkind) number of extra dynamics grid fields +! nhead Integer(sigio_intkind) number of header records +! ndata Integer(sigio_intkind) number of data records +! lhead Integer(sigio_intkind)(nhead) header record lengths +! ldata Integer(sigio_intkind)(ndata) data record lengths +! +! sigio_data Sigma file data fields +! hs Real(sigio_realkind)(:) pointer to spectral +! coefficients of surface height in m +! ps Real(sigio_realkind)(:) pointer to spectral +! coefficients of log of surface pressure over 1 kPa +! t Real(sigio_realkind)(:,:) pointer to spectral +! coefficients of virtual temperature by level in K +! d Real(sigio_realkind)(:,:) pointer to spectral +! coefficients of divergence by level in 1/second +! z Real(sigio_realkind)(:,:) pointer to spectral +! coefficients of vorticity by level in 1/second +! q Real(sigio_realkind)(:,:,:) pointer to spectral +! coefficients of tracers by level and tracer number +! in specific densities +! xgr Real(sigio_realkind)(:,:,:) pointer to extra grid fields +! by longitude, latitude and number of extra physics grid fields +! xss Real(sigio_realkind)(:) pointer to scalar array +! xga Real(sigio_realkind)(:,:,:) pointer to extra dynamics grid fields +! by longitude, latitude and number of extra grid fields +! +! sigio_dbta Sigma file longreal data fields +! hs Real(sigio_dblekind)(:) pointer to spectral +! coefficients of surface height in m +! ps Real(sigio_dblekind)(:) pointer to spectral +! coefficients of log of surface pressure over 1 kPa +! t Real(sigio_dblekind)(:,:) pointer to spectral +! coefficients of virtual temperature by level in K +! d Real(sigio_dblekind)(:,:) pointer to spectral +! coefficients of divergence by level in 1/second +! z Real(sigio_dblekind)(:,:) pointer to spectral +! coefficients of vorticity by level in 1/second +! q Real(sigio_dblekind)(:,:,:) pointer to spectral +! coefficients of tracers by level and tracer number +! in specific densities +! xgr Real(sigio_dblekind)(:,:,:) pointer to extra physics grid fields +! by longitude, latitude and number of extra grid fields +! xss Real(sigio_dblekind)(:) pointer to scalar array +! xga Real(sigio_dblekind)(:,:,:) pointer to extra dynamics grid fields +! by longitude, latitude and number of extra grid fields +! +! Public Subprograms: +! sigio_modpr Compute model pressures +! im Integer(sigio_intkind) input number of points +! ix Integer(sigio_intkind) input first dimension +! km Integer(sigio_intkind) input number of levels +! nvcoord Integer(sigio_intkind) input number of vertical coords +! idvc Integer(sigio_intkind) input vertical coordinate id +! (1 for sigma and 2 for hybrid) +! idsl Integer(sigio_intkind) input type of sigma structure +! (1 for phillips or 2 for mean) +! vcoord Real(sigio_realkind)(km+1,nvcoord) input vertical coords +! for idvc=1, nvcoord=1: sigma interface +! for idvc=2, nvcoord=2: hybrid interface a and b +! iret Integer(sigio_intkind) output return code +! ps Real(sigio_realkind)(ix) input optional surface pressure (Pa) +! tv Real(sigio_realkind)(ix,km) input optional virtual temperature (K) +! pd Real(sigio_realkind)(ix,km) output optional delta pressure (Pa) +! pm Real(sigio_realkind)(ix,km) output optional layer pressure (Pa) +! +! +! Remarks: +! (1) The sigma file format follows: +! For ivs=198410: +! ON85 label (32 bytes) +! Header information record containing +! real forecast hour, initial date, sigma interfaces, sigma levels, +! padding to allow for 100 levels, and finally 44 identifier words +! containing JCAP, LEVS, NTRAC, IREALF, etc. (250 4-byte words) +! (word size in the remaining records depends on the value of IREALF) +! Orography (NC words, where NC=(JCAP+1)*(JCAP+2)) +! Log surface pressure (NC words) +! Temperature (LEVS records of NC words) +! Divergence & Vorticity interleaved (2*LEVS records of NC words) +! Tracers (LEVS*NTRAC records of NC words) +! Extra grid fields (NXGR records of LONB*LATB words) +! For ivs=200509: +! Label containing +! 'GFS ','SIG ',ivs,nhead,ndata,reserved(3) (8 4-byte words) +! Header records +! lhead(nhead),ldata(ndata) (nhead+ndata 4-byte words) +! fhour, idate(4), jcap, levs, itrun, iorder, irealf, igen, +! latf, lonf, latb, lonb, latr, lonr, ntrac, nvcoord, +! icen2, iens(2), idpp, idsl, idvc, idvm, idvt, idrun, idusr, +! pdryini, ncldt, ixgr, ixga,reserved(17) (50 4-byte words) +! vcoord((levs+1)*nvcoord 4-byte words) +! cfvars(5+ntrac 8-byte character words) +! Data records (word size depends on irealf) +! orography (nc words, where nc=(jcap+1)*(jcap+2)) +! log surface pressure (nc words) +! temperature (levs records of nc words) +! divergence (levs records of nc words) +! vorticity (levs records of nc words) +! tracers (levs*ntrac records of nc words) +! scalars (nxss words) +! extra physics grid fields (nxgr records of lonb*latb words) +! extra scalars (nxss words) +! extra dynamics grid fields (nxga records of lonf*latf words) +! +! (2) Possible return codes: +! 0 Successful call +! -1 Open or close I/O error +! -2 Header record I/O error (possible EOF) +! -3 Allocation or deallocation error +! -4 Data record I/O error +! -5 Insufficient data dimensions allocated +! -6 Attempted to read a NEMSIO file +! +! Examples: +! (1) Read the entire sigma file 'sigf24' and +! print out the global mean temperature profile. +! +! use sigio_module +! type(sigio_head):: head +! type(sigio_data):: data +! call sigio_srohdc(11,'sigf24',head,data,iret) +! print '(f8.2)',data%t(1,head%levs:1:-1)/sqrt(2.) +! end +! +! Attributes: +! Language: Fortran 90 +! +!$$$ + implicit none + private +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! Public Variables + integer,parameter,public:: sigio_lhead1=32 + integer,parameter,public:: sigio_intkind=4,sigio_realkind=4,sigio_dblekind=8 + integer,parameter,public:: sigio_charkind=8 + real(sigio_intkind),parameter,public:: sigio_intfill=-9999_sigio_intkind + real(sigio_realkind),parameter,public:: sigio_realfill=-9999._sigio_realkind + real(sigio_dblekind),parameter,public:: sigio_dblefill=-9999._sigio_dblekind +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! Public Types + type,public:: sigio_head + character(sigio_lhead1):: clabsig=' ' + real(sigio_realkind):: fhour=sigio_realfill + integer(sigio_intkind):: idate(4)=sigio_intfill + real(sigio_realkind):: si(101)=sigio_realfill + real(sigio_realkind):: sl(100)=sigio_realfill + real(sigio_realkind):: ak(101)=sigio_realfill + real(sigio_realkind):: bk(101)=sigio_realfill + integer(sigio_intkind):: jcap=sigio_intfill + integer(sigio_intkind):: levs=sigio_intfill + integer(sigio_intkind):: itrun=sigio_intfill + integer(sigio_intkind):: iorder=sigio_intfill + integer(sigio_intkind):: irealf=sigio_intfill + integer(sigio_intkind):: igen=sigio_intfill + integer(sigio_intkind):: latf=sigio_intfill + integer(sigio_intkind):: lonf=sigio_intfill + integer(sigio_intkind):: latb=sigio_intfill + integer(sigio_intkind):: lonb=sigio_intfill + integer(sigio_intkind):: latr=sigio_intfill + integer(sigio_intkind):: lonr=sigio_intfill + integer(sigio_intkind):: ntrac=sigio_intfill + integer(sigio_intkind):: icen2=sigio_intfill + integer(sigio_intkind):: iens(2)=sigio_intfill + integer(sigio_intkind):: idpp=sigio_intfill + integer(sigio_intkind):: idsl=sigio_intfill + integer(sigio_intkind):: idvc=sigio_intfill + integer(sigio_intkind):: idvm=sigio_intfill + integer(sigio_intkind):: idvt=sigio_intfill + integer(sigio_intkind):: idrun=sigio_intfill + integer(sigio_intkind):: idusr=sigio_intfill + real(sigio_realkind):: pdryini=sigio_realfill + integer(sigio_intkind):: ncldt=sigio_intfill + integer(sigio_intkind):: ixgr=sigio_intfill + integer(sigio_intkind):: ixga=sigio_intfill + integer(sigio_intkind):: ivs=sigio_intfill + integer(sigio_intkind):: nvcoord=sigio_intfill + real(sigio_realkind),allocatable:: vcoord(:,:) + character(sigio_charkind),allocatable:: cfvars(:) + integer(sigio_intkind):: nxgr=sigio_intfill + integer(sigio_intkind):: nxss=sigio_intfill + integer(sigio_intkind):: nxga=sigio_intfill + integer(sigio_intkind):: nhead=sigio_intfill + integer(sigio_intkind):: ndata=sigio_intfill + integer(sigio_intkind),allocatable:: lhead(:) + integer(sigio_intkind),allocatable:: ldata(:) + real(sigio_realkind), allocatable :: cpi(:), ri(:) +! real(sigio_realkind):: cpi(100)=sigio_realfill +! real(sigio_realkind):: ri(100)=sigio_realfill + end type + type,public:: sigio_data + real(sigio_realkind),pointer:: hs(:)=>null() + real(sigio_realkind),pointer:: ps(:)=>null() + real(sigio_realkind),pointer:: t(:,:)=>null() + real(sigio_realkind),pointer:: d(:,:)=>null() + real(sigio_realkind),pointer:: z(:,:)=>null() + real(sigio_realkind),pointer:: q(:,:,:)=>null() + real(sigio_realkind),pointer:: xgr(:,:,:)=>null() + real(sigio_realkind),pointer:: xss(:)=>null() + real(sigio_realkind),pointer:: xga(:,:,:)=>null() + end type + type,public:: sigio_dbta + real(sigio_dblekind),pointer:: hs(:)=>null() + real(sigio_dblekind),pointer:: ps(:)=>null() + real(sigio_dblekind),pointer:: t(:,:)=>null() + real(sigio_dblekind),pointer:: d(:,:)=>null() + real(sigio_dblekind),pointer:: z(:,:)=>null() + real(sigio_dblekind),pointer:: q(:,:,:)=>null() + real(sigio_dblekind),pointer:: xgr(:,:,:)=>null() + real(sigio_dblekind),pointer:: xss(:)=>null() + real(sigio_dblekind),pointer:: xga(:,:,:)=>null() + end type +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! Public Subprograms + public sigio_modpr +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! Private Variables +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! Private Types + type sigio_head2 + sequence + real(sigio_realkind):: fhour + integer(sigio_intkind):: idate(4) + real(sigio_realkind):: sisl(2*100+1) + real(sigio_realkind):: ext(44) + end type +contains +!------------------------------------------------------------------------------- + subroutine sigio_modpr(im,ix,km,nvcoord,idvc,idsl,vcoord,iret,& + ps,t,pd,dpddps,dpddt,pm,dpmdps,dpmdt) + implicit none + integer,intent(in):: im,ix,km,nvcoord,idvc,idsl + real,intent(in):: vcoord(km+1,nvcoord) + integer,intent(out):: iret + real,intent(in),optional:: ps(ix),t(ix,km) + real,intent(out),optional:: pd(ix,km),pm(ix,km) + real,intent(out),optional:: dpddps(ix,km),dpddt(ix,km) + real,intent(out),optional:: dpmdps(ix,km),dpmdt(ix,km) + real(sigio_dblekind),parameter:: rocp=287.05/1004.6,rocpr=1/rocp + real(sigio_dblekind),parameter:: t00=300. + integer id1,id2 + real(sigio_dblekind) pid(im),dpiddps(im),dpiddt(im),tid(im),pidk(im) + real(sigio_dblekind) piu,dpiudps,dpiudt,tiu,piuk + real(sigio_dblekind) pmm,dpmdpid,dpmdpiu + real(sigio_dblekind) pmk + integer i,k +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + if((idvc.eq.0.or.idvc.eq.1).and.nvcoord.eq.1.and.present(ps)) then + id1=11 + elseif(idvc.eq.2.and.nvcoord.eq.2.and.present(ps)) then + id1=22 + elseif(idvc.eq.3.and.nvcoord.eq.3.and.all(vcoord(:,3).eq.0).and.present(ps)) then + id1=22 + elseif(idvc.eq.3.and.nvcoord.eq.2.and.present(ps).and.present(t)) then + id1=32 + elseif(idvc.eq.3.and.nvcoord.eq.3.and.present(ps).and.present(t)) then + id1=33 + else + id1=0 + endif + if(idsl.eq.0.or.idsl.eq.1) then + id2=1 + elseif(idsl.eq.2) then + id2=2 + else + id2=0 + endif + iret=0 +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + if(id1.gt.0.and.id2.gt.0) then +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) + do i=1,im + pid(i)=ps(i) + dpiddps(i)=1 + dpiddt(i)=0 + tid(i)=0 + if(id2.eq.1) pidk(i)=pid(i)**rocp + enddo +!$OMP END PARALLEL DO + +!!$OMP PARALLEL DO DEFAULT(SHARED) & +!!$OMP& PRIVATE(i,k,piu,dpiudps,dpiudt,tiu,piuk,pmk,pmm,dpmdpid,dpmdpiu) & +!!$OMP& PRIVATE(pid,dpiddps,dpiddt,tid,pidk) + + do k=1,km +!$OMP PARALLEL DO DEFAULT(SHARED) & +!$OMP& PRIVATE(i,piu,dpiudps,dpiudt,tiu,piuk,pmk,pmm,dpmdpid,dpmdpiu) + do i=1,im + select case(id1) + case(11) + piu=vcoord(k+1,1)*ps(i) + dpiudps=vcoord(k+1,1) + dpiudt=0 + case(22) + piu=vcoord(k+1,1)+vcoord(k+1,2)*ps(i) + dpiudps=vcoord(k+1,2) + dpiudt=0 + case(32) + tiu=(t(i,k)+t(i,min(k+1,km)))/2 + piu=vcoord(k+1,2)*ps(i)+vcoord(k+1,1)*(tiu/t00)**rocpr + dpiudps=vcoord(k+1,2) + dpiudt=vcoord(k+1,1)*(tiu/t00)**rocpr*rocpr/tiu + if(k.lt.km) dpiudt=dpiudt/2 + case(33) + tiu=(t(i,k)+t(i,min(k+1,km)))/2 + piu=vcoord(k+1,1)+vcoord(k+1,2)*ps(i)+vcoord(k+1,3)*(tiu/t00)**rocpr + dpiudps=vcoord(k+1,2) + dpiudt=vcoord(k+1,3)*(tiu/t00)**rocpr*rocpr/tiu + if(k.lt.km) dpiudt=dpiudt/2 + end select + if(present(pd)) pd(i,k)=pid(i)-piu + if(present(dpddps)) dpddps(i,k)=dpiddps(i)-dpiudps + if(present(dpddt)) dpddt(i,k)=dpiddt(i)-dpiudt +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + select case(id2) + case(1) + piuk=piu**rocp + pmk=(pid(i)*pidk(i)-piu*piuk)/((rocp+1)*(pid(i)-piu)) + pmm=pmk**rocpr + dpmdpid=rocpr*pmm/(pid(i)-piu)*(pidk(i)/pmk-1) + dpmdpiu=rocpr*pmm/(pid(i)-piu)*(1-piuk/pmk) + case(2) + pmm=(pid(i)+piu)/2 + dpmdpid=0.5 + dpmdpiu=0.5 + end select + if(present(pm)) pm(i,k)=pmm + if(present(dpmdps)) dpmdps(i,k)=dpmdpid*dpiddps(i)+dpmdpiu*dpiudps + if(present(dpmdt)) dpmdt(i,k)=dpmdpid*dpiddt(i)+dpmdpiu*dpiudt +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + pid(i)=piu + dpiddps(i)=dpiudps + dpiddt(i)=dpiudt + tid(i)=tiu + if(id2.eq.1) pidk(i)=piuk + enddo +!$OMP END PARALLEL DO + enddo +!!$OMP END PARALLEL DO + else + if(id1.le.0) iret=iret+1 + if(id2.le.0) iret=iret+2 + endif +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + end subroutine +!------------------------------------------------------------------------------- +end module diff --git a/src/gfs_bufr.fd/modstuff1.f b/src/gfs_bufr.fd/modstuff1.f index 95d41383..7bc075a5 100644 --- a/src/gfs_bufr.fd/modstuff1.f +++ b/src/gfs_bufr.fd/modstuff1.f @@ -11,6 +11,7 @@ subroutine modstuff(km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,& ! ! Program history log: ! 1999-10-18 Mark Iredell +! 2024-08-23 Bo Cui Replace sigio_module with the simplified module modpr_module ! ! Usage: call modstuff(km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,& ! pd,pi,pm,aps,apm,os,om,px,py) @@ -41,7 +42,7 @@ subroutine modstuff(km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,& ! Language: Fortran 90 ! !$$$ - use sigio_module + use modpr_module implicit none integer,intent(in):: km,idvc,idsl,nvcoord real,intent(in):: vcoord(km+1,nvcoord) From 15bf330b084d8991c28ab260141306718d033f1e Mon Sep 17 00:00:00 2001 From: "Bo.Cui" Date: Fri, 23 Aug 2024 15:12:23 -0400 Subject: [PATCH 5/6] Add module modpr_module.f90 --- src/gfs_bufr.fd/modpr_module.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/gfs_bufr.fd/modpr_module.f90 b/src/gfs_bufr.fd/modpr_module.f90 index 70ffbec9..ef9f9290 100644 --- a/src/gfs_bufr.fd/modpr_module.f90 +++ b/src/gfs_bufr.fd/modpr_module.f90 @@ -10,7 +10,8 @@ module modpr_module ! - Removed all other subroutines and associated code that were present in sigio_module.f ! ! Reason for modification: -! - To simplify the module by including only the necessary subroutines +! - To remove the dependency on sigio +! - To simplify the module by including only the necessary one ! ! Program History Log: ! 2024-08-23 Bo Cui: Simplified the module by retaining only sigio_modpr subroutine From 3004dff63ebddcdfa24708dcc24620c1fe9d2d87 Mon Sep 17 00:00:00 2001 From: "Bo.Cui" Date: Fri, 23 Aug 2024 16:51:03 -0400 Subject: [PATCH 6/6] remove "mpif.h" from meteorg.f and read_netcdf.f --- src/gfs_bufr.fd/meteorg.f | 1 - src/gfs_bufr.fd/read_netcdf.f | 1 - 2 files changed, 2 deletions(-) diff --git a/src/gfs_bufr.fd/meteorg.f b/src/gfs_bufr.fd/meteorg.f index 83ed46a0..6e1f81b4 100644 --- a/src/gfs_bufr.fd/meteorg.f +++ b/src/gfs_bufr.fd/meteorg.f @@ -71,7 +71,6 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! use funcphys, only : gfuncphys use funcphys implicit none - include 'mpif.h' integer :: nfile,npoint,levs,kdim,nsoil character(len=10) :: dim_nam integer :: nfile1 diff --git a/src/gfs_bufr.fd/read_netcdf.f b/src/gfs_bufr.fd/read_netcdf.f index a024323b..05af922b 100644 --- a/src/gfs_bufr.fd/read_netcdf.f +++ b/src/gfs_bufr.fd/read_netcdf.f @@ -5,7 +5,6 @@ subroutine read_netcdf(ncid,im,jm,levs, use netcdf implicit none - include 'mpif.h' character(len=20),intent(in) :: VarName character(len=3),intent(in) :: Zreverse integer,intent(in) :: ncid,im,jm,levs