diff --git a/src/gfs_bufr.fd/CMakeLists.txt b/src/gfs_bufr.fd/CMakeLists.txt index 7df490d0..217a14d1 100644 --- a/src/gfs_bufr.fd/CMakeLists.txt +++ b/src/gfs_bufr.fd/CMakeLists.txt @@ -9,14 +9,14 @@ list(APPEND fortran_src meteorg.f mstadb.f newsig1.f - read_nemsio.f - #read_netcdf.f - read_netcdf_p.f + read_netcdf.f + #read_netcdf_p.f rsearch.f svp.f 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/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/gfsbufr.f b/src/gfs_bufr.fd/gfsbufr.f index c1bb1d35..5d31625c 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,19 @@ 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 +64,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 +76,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 +91,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 +143,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 +164,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 +190,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/meteorg.f b/src/gfs_bufr.fd/meteorg.f index d4fa338d..6e1f81b4 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,9 @@ 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 +! 2024-08-23 Bo Cui Replace sigio_module with the simplified module modpr_module +! ! ! USAGE: CALL PROGRAM meteorg ! INPUT: @@ -46,6 +49,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,18 +65,14 @@ subroutine meteorg(npoint,rlat,rlon,istat,cstat,elevstn, ! !$$$ use netcdf - use nemsio_module - use sigio_module + use modpr_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 :: 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 +83,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 +101,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 +120,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 +128,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 +141,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 +161,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 +188,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 +207,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 +218,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 +229,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 +243,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 +257,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 +271,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 +285,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 +299,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 +355,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 +412,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 +434,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 +445,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 +457,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 +469,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 +478,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 +525,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 +535,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 +556,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 +582,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 +600,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 +609,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 +639,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 +917,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 +1115,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/modpr_module.f90 b/src/gfs_bufr.fd/modpr_module.f90 new file mode 100644 index 00000000..ef9f9290 --- /dev/null +++ b/src/gfs_bufr.fd/modpr_module.f90 @@ -0,0 +1,457 @@ +!------------------------------------------------------------------------------- +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 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 +! +! 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) 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