WRFPOST.F

Go to the documentation of this file.

 
!---------------------------------------------------------------------
!---------------------------------------------------------------------
      PROGRAM wrfpost
 
!
!
!============================================================================================================
!
!     This is an MPI code. All array indexing is with respect to the global indices. Loop indices 
!     look as follows for N MPI tasks.
!
!
!
!  Original                                            New
!  Index                                               Index
!
!   JM ----------------------------------------------- JEND
! JM-1 -                                             - JEND_M
! JM-2 -               MPI TASK N-1                  - JEND_M2
!      -                                             -
!      -                                             -
!      ----------------------------------------------- JSTA, JSTA_M, JSTA_M2
!      ----------------------------------------------- JEND, JEND_M, JEND_M2
!      -                                             -
!      -               MPI TASK N-2                  -
!      -                                             -
!      -                                             -
!      ----------------------------------------------- JSTA, JSTA_M, JSTA_M2
!
!                           .
!                           .
!                           .
!
!      ----------------------------------------------- JEND, JEND_M, JEND_M2
!      -                                             -
!      -               MPI TASK 1                    -
!      -                                             -
!      -                                             -
!      ----------------------------------------------- JSTA, JSTA_M, JSTA_M2
!      ----------------------------------------------- JEND, JEND_M, JEND_M2
!      -                                             - 
!      -               MPI TASK 0                    - 
!    3 -                                             - JSTA_M2
!    2 -                                             - JSTA_M
!    1 ----------------------------------------------- JSTA
!
!     1                                              IM               
!
!
!     Jim Tuccillo
!     Jan 2000
!
!     README - Jim Tuccillo Feb 2001
! 
!     Many common blocks have been replaced by modules to support Fortran
!     "allocate" commands. Many of the 3-D arrays are now allocated to be the
!     exact size required based on the number of MPI tasks. The dimensioning will be 
!        x ( im,jsta_2l:jend_2u,lm)
!     Most 2-D arrays continue to be dimensioned (im,jm). This is fine but please be aware 
!     that the EXCH routine for arrays dimensioned (im,jm) is different than arrays dimensioned
!     (im,jsta_2l:jend_2u). Also, be careful about passing any arrays dimensioned
!     (im,jst_2l:jend_2u,lm). See examples in the code as to the correct calling sequence and
!     EXCH routine to use.
!
!
!     ASYNCHRONOUS I/O HAS BEEN ADDED. THE LAST MPI TASK DOES THE I/O. IF THERE IS
!     ONLY ONE MPI TASK THN TASK ) DOES THE I/O.
!     THE CODE HAS GOTTEN A LITTLE KLUDGY. BASICLY, IM, IMX and IMOUT MUST BE EQUAL
!     AND REPRESENT THE VALUE USED IN THE MODEL. THE SAME HOLDS FOR JM, JMX and JMOUT.
!
!     Jim Tuccillo June 2001
!
!
!===========================================================================================
!
      use netcdf
#if defined(BUILD_WITH_NEMSIO)
      use nemsio_module, only: nemsio_getheadvar, nemsio_gfile, nemsio_init, nemsio_open, &
                               nemsio_getfilehead,nemsio_close
#endif
      use ctlblk_mod,    only: filenameaer, me, num_procs, num_servers, mpi_comm_comp, datestr,      &
              mpi_comm_inter, filename, ioform, grib, idat, filenameflux, filenamed3d, gdsdegr,      &
              spldef, modelname, ihrst, lsmdef,vtimeunits, tprec, pthresh, datahandle, im, jm, lm,   &
              lp1, lm1, im_jm, isf_surface_physics, nsoil, spl, lsmp1, global, imp_physics,          &
              ista, iend, ista_m, iend_m, ista_2l, iend_2u,                                          &
              jsta, jend, jsta_m, jend_m, jsta_2l, jend_2u, novegtype, icount_calmict, npset, datapd,&
              lsm, fld_info, etafld2_tim, eta2p_tim, mdl2sigma_tim, cldrad_tim, miscln_tim,          &
              mdl2agl_tim, mdl2std_tim, mdl2thandpv_tim, calrad_wcloud_tim,nasa_on,gccpp_on,         &
              fixed_tim, time_output, imin, surfce2_tim, komax, ivegsrc, d3d_on, gocart_on,rdaod,    &
              readxml_tim, spval, fullmodelname, submodelname, hyb_sigp, filenameflat, aqf_on,numx,  &
              run_ifi_tim, slrutah_on, d2d_chem, gtg_on, method_blsn
      use grib2_module,   only: gribit2,num_pset,nrecout,first_grbtbl,grib_info_finalize
      use upp_ifi_mod, only: write_ifi_debug_files
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      implicit none
!
#if defined(BUILD_WITH_NEMSIO)
      type(nemsio_gfile) :: nfile,ffile,rfile
#endif
      include "mpif.h"
!
!     DECLARE VARIABLES.
!     
!     SET HEADER WRITER FLAGS TO TRUE.
!
!temporary vars
!
      real(kind=8) :: time_initpost=0.,initpost_tim=0.,btim,bbtim
      real            rinc(5), untcnvt
      integer      :: status=0,iostatusd3d=0,iostatusflux=0
      integer i,j,iii,l,k,ierr,nrec,ist,lusig,idrt,ncid3d,ncid2d,varid
      integer      :: prntsec,iim,jjm,llm,ioutcount,itmp,iret,iunit,        &
                      iunitd3d,iyear,imn,iday,lcntrl,ieof
      integer      :: iostatusaer
      logical      :: popascal
!
      integer      :: kpo,kth,kpv
      real,dimension(komax) :: po,th,pv
      namelist/nampgb/kpo,po,kth,th,kpv,pv,filenameaer,d3d_on,gocart_on,gccpp_on, nasa_on,gtg_on,method_blsn,popascal &
                     ,hyb_sigp,rdaod,d2d_chem, aqf_on,slrutah_on, vtimeunits,numx,write_ifi_debug_files
      integer      :: itag_ierr
      namelist/model_inputs/filename,ioform,grib,datestr,modelname,submodelname &
                     ,filenameflux,filenameflat
 
      character startdate*19,sysdepinfo*80,iowrfname*3,post_fname*255
      character cgar*1,cdum*4,line*10
!
!------------------------------------------------------------------------------
!     START HERE
!
      call start()
!
!     INITIALIZE MPI
      
      CALL setup_servers(me,                         &
     &                   num_procs,                  &
     &                   num_servers,                &
     &                   mpi_comm_comp,              &
     &                   mpi_comm_inter)
!
!     ME IS THE RANK
!     NUM_PROCS IS THE NUMBER OF TASKS DOING POSTING
!     NUM_SERVERS IS ONE IF THERE ARE MORE THAN ONE TOTAL MPI TASKS, OTHERWISE ZERO
!     MPI_COMM_COMP IS THE INTRACOMMUNICATOR
!     MPI_COMM_INTER IS THE INTERCOMMUNICATOR FOR COMMUNCATION BETWEEN TASK 0 OF THE
!        TASKS DOING THE POSTING AND THE I/O SERVER
!
!
!     IF WE HAVE MORE THAN 1 MPI TASK THEN WE WILL FIRE UP THE IO SERVER
!     THE LAST TASK ( IN THE CONTEXT OF MPI_COMM_WORLD ) IS THE I/O SERVER
!
      if (me == 0) CALL w3tagb('nems     ',0000,0000,0000,'np23   ')
 
      if ( me >= num_procs ) then
!
         call server
!
      else
        spval = 9.9e10
!
!**************************************************************************
!KaYee: Read itag in Fortran Namelist format
!Set default 
       submodelname='NONE'
!Set control file name
       filenameflat='postxconfig-NT.txt'
       numx=1
!open namelist
       open(5,file='itag')
       read(5,nml=model_inputs,iostat=itag_ierr,err=888)
888    if (itag_ierr /= 0) then
       print*,'Incorrect namelist variable(s) found in the itag file,stopping.'
       stop
       endif
       if (me == 0) write(6, model_inputs)
       
!       if(MODELNAME == 'NMM')then
!        read(5,1114) VTIMEUNITS
! 1114   format(a4)
!        if (me==0) print*,'VALID TIME UNITS = ', VTIMEUNITS
!       endif
!
 303  format('MODELNAME="',a,'" SUBMODELNAME="',a,'"')
 
        if(me==0) write(*,*)'MODELNAME: ', modelname, submodelname
 
! assume for now that the first date in the stdin file is the start date
        read(datestr,300) iyear,imn,iday,ihrst,imin
        if (me==0) write(*,*) 'in WRFPOST iyear,imn,iday,ihrst,imin',                &
                    iyear,imn,iday,ihrst,imin
 300    format(i4,1x,i2,1x,i2,1x,i2,1x,i2)
 
        idat(1) = imn
        idat(2) = iday
        idat(3) = iyear
        idat(4) = ihrst
        idat(5) = imin
 
 111    format(a256)
 112    format(a19)
 113    format(a20)
 114    format(a8)
 120    format(a5)
 121    format(a4)
 
!KaYee: Read in GFS/FV3 runs in Fortran Namelist Format.
      if(grib=='grib2') then
        gdsdegr = 1.d6
      endif
! 
! set default for kpo, kth, th, kpv, pv     
        kpo = 0
        po  = 0
        kth = 6
        th  = (/310.,320.,350.,450.,550.,650.,(0.,k=kth+1,komax)/) ! isentropic level to output
        kpv = 8
        pv  = (/0.5,-0.5,1.0,-1.0,1.5,-1.5,2.0,-2.0,(0.,k=kpv+1,komax)/)
 
        hyb_sigp    = .true.
        d3d_on      = .false.
        gocart_on   = .false.
        gccpp_on    = .false.
        nasa_on     = .false.
        aqf_on      = .false.
        slrutah_on  = .false.
        gtg_on   = .false.
        method_blsn = .true.
        popascal    = .false.
        filenameaer = ''
        rdaod       = .false.
        d2d_chem     = .false.
        vtimeunits  = ''
 
        read(5,nampgb,iostat=iret,end=119)
 119    continue
        if (me == 0) write(6, nampgb)
       if(mod(num_procs,numx)/=0) then
         if (me==0) then
           print*,'total proces, num_procs=', num_procs 
           print*,'number of subdomain in x direction, numx=', numx 
           print*,'remainder of num_procs/numx = ', mod(num_procs,numx)
           print*,'Warning!!! the remainder of num_procs/numx is not 0, reset numx=1 &
     &             in this run or you adjust numx in the itag file to restart'
         endif
!        stop 9999
         numx=1
         if(me == 0) print*,'Warning!!!  Reset numx as 1, numx=',numx
       endif
       if(numx>num_procs/2) then
         if (me==0) then
           print*,'total proces, num_procs=', num_procs
           print*,'number of subdomain in x direction, numx=', numx
           print*,'Warning!!! numx cannot exceed num_procs/2, reset numx=1 in this run'
           print*,'or you adjust numx in the itag file to restart'
         endif
         numx=1
         if(me == 0) print*,'Warning!!!  Reset numx as 1, numx=',numx
       endif     
        if(me == 0) then
          print*,'komax,iret for nampgb= ',komax,iret 
          print*,'komax,kpo,kth,th,kpv,pv,fileNameAER,nasa_on,popascal= ',komax,kpo        &
     &           ,kth,th(1:kth),kpv,pv(1:kpv),trim(filenameaer),nasa_on,popascal
          print*,'NUM_PROCS=',num_procs
          print*,'numx= ',numx
        endif
 
        IF(trim(ioform) /= 'netcdfpara' .AND. trim(ioform) /= 'netcdf' ) THEN
          numx=1
          if(me == 0) print*,'2D decomposition only supports netcdfpara IO.'
          if(me == 0) print*,'Reset numx= ',numx
        ENDIF
 
        IF(modelname /= 'FV3R' .AND. modelname /= 'GFS') THEN
          numx=1
          if(me == 0) print*,'2D decomposition only supports GFS and FV3R.'
          if(me == 0) print*,'Reset numx= ',numx
        ENDIF
 
! set up pressure level from POSTGPVARS or DEFAULT
        if(kpo == 0) then
! use default pressure levels
          if(me == 0) then
            print*,'using default pressure levels,spldef=',(spldef(l),l=1,lsmdef)
          endif
          lsm = lsmdef
          do l=1,lsm
            spl(l) = spldef(l)
          end do
        else
! use POSTGPVARS
          if(me == 0) then
            print*,'using pressure levels from POSTGPVARS'
          endif
          lsm = kpo
          if( .not. popascal ) then
            untcnvt = 100.
          else
            untcnvt = 1.
          endif
          if(po(lsm) < po(1))then ! post logic assumes asscending
            do l=1,lsm
              spl(l) = po(lsm-l+1)*untcnvt 
            end do
          else
            do l=1,lsm
              spl(l) = po(l)*untcnvt
            end do
          end if
        end if
        lsmp1 = lsm+1
      
 116    continue
 
! set PTHRESH for different models
        if(modelname == 'NMM')then
          pthresh = 0.000004
        else
          pthresh = 0.000001
        end if  
!Chuang: add dynamical allocation
        if(trim(ioform) == 'netcdf' .OR. trim(ioform) == 'netcdfpara') THEN
         IF(modelname == 'NCAR' .OR. modelname == 'RAPR' .OR. modelname == 'NMM') THEN
          call ext_ncd_ioinit(sysdepinfo,status)
          call ext_ncd_open_for_read( trim(filename), 0, 0, " ",          &
            datahandle, status)
          if ( status /= 0 ) then
            print*,'error opening ',filename, ' Status = ', status ; stop
          endif
          call ext_ncd_get_dom_ti_integer(datahandle                      &
            ,'WEST-EAST_GRID_DIMENSION',iim,1,ioutcount, status )
          im = iim-1
          call ext_ncd_get_dom_ti_integer(datahandle                      &
             ,'SOUTH-NORTH_GRID_DIMENSION',jjm,1,ioutcount, status )
          jm = jjm-1
          call ext_ncd_get_dom_ti_integer(datahandle                      &
            ,'BOTTOM-TOP_GRID_DIMENSION',llm,1,ioutcount, status )
          lm    = llm-1
          lp1   = lm+1
          lm1   = lm-1
          im_jm = im*jm
       
! Read and set global value for surface physics scheme
          call ext_ncd_get_dom_ti_integer(datahandle                      &
            ,'SF_SURFACE_PHYSICS',itmp,1,ioutcount, status )
 
          ! J. Kenyon / 23 Aug 2024: GSL MPAS output erroneously
          ! indicates SF_SURFACE_PHYSICS = 0 (should be 3), so overwrite
          ! here:
          if (modelname == 'RAPR' .and. submodelname == 'MPAS') itmp = 3 
 
          isf_surface_physics = itmp
! set NSOIL to 4 as default for NOAH but change if using other
! SFC scheme
          nsoil = 4
          IF(itmp      == 1) then !thermal diffusion scheme
            nsoil = 5
          ELSE IF(itmp == 3) then ! RUC LSM
            nsoil = 9
          ELSE IF(itmp == 7) then ! Pleim Xu
            nsoil = 2
          END IF
 
          call ext_ncd_ioclose ( datahandle, status )
         ELSE
! use parallel netcdf lib directly to read FV3 output in netCDF
          spval = 9.99e20
          status = nf90_open(trim(filename),ior(nf90_nowrite,nf90_mpiio), &
                   ncid3d,comm=mpi_comm_world,info=mpi_info_null)
          if ( status /= 0 ) then
            print*,'error opening ',filename, ' Status = ', status 
            stop
          endif
          status = nf90_open(trim(filenameflux),ior(nf90_nowrite,nf90_mpiio), &
                   ncid2d,comm=mpi_comm_world,info=mpi_info_null)
          if ( status /= 0 ) then
            print*,'error opening ',filenameflux, ' Status = ', status
            stop
          endif
! read in LSM index and nsoil here
          status=nf90_get_att(ncid2d,nf90_global,'landsfcmdl', isf_surface_physics)
          if(status/=0)then
            print*,'landsfcmdl not found; assigning to 2'
            isf_surface_physics=2 !set LSM physics to 2 for NOAH
          endif
          if(isf_surface_physics<2)then
            isf_surface_physics=2 !set LSM physics to 2 for NOAH
          endif
          status=nf90_get_att(ncid2d,nf90_global,'nsoil', nsoil)
          if(status/=0)then
            print*,'nsoil not found; assigning to 4'
            nsoil=4 !set nsoil to 4 for NOAH
          endif
! read imp_physics
          status=nf90_get_att(ncid2d,nf90_global,'imp_physics',imp_physics)
          if(status/=0)then
            print*,'imp_physics not found; assigning to GFDL 11'
            imp_physics=11
          endif
! get dimesions
          status = nf90_inq_dimid(ncid3d,'grid_xt',varid)
          if ( status /= 0 ) then
           print*,status,varid
           stop 1
          end if
          status = nf90_inquire_dimension(ncid3d,varid,len=im)
          if ( status /= 0 ) then
           print*,status
           stop 1      
          end if   
          status = nf90_inq_dimid(ncid3d,'grid_yt',varid)
          if ( status /= 0 ) then
           print*,status,varid
           stop 1
          end if
          status = nf90_inquire_dimension(ncid3d,varid,len=jm)
          if ( status /= 0 ) then
           print*,status
           stop 1
          end if
          status = nf90_inq_dimid(ncid3d,'pfull',varid)
          if ( status /= 0 ) then
           print*,status,varid
           stop 1
          end if
          status = nf90_inquire_dimension(ncid3d,varid,len=lm)
          if ( status /= 0 ) then
           print*,status
           stop 1
          end if
          lp1   = lm+1
          lm1   = lm-1
          im_jm = im*jm
! set NSOIL to 4 as default for NOAH but change if using other
! SFC scheme
!          NSOIL = 4
         END IF 
 
        ELSE IF(trim(ioform) == 'binary'       .OR.                       &
                trim(ioform) == 'binarympiio' ) THEN
          print*,'WRF Binary format is no longer supported'
          stop 9996
! NEMSIO format
#if defined(BUILD_WITH_NEMSIO)
        ELSE IF(trim(ioform) == 'binarynemsio' .or.                        &
          trim(ioform) == 'binarynemsiompiio' )THEN
      
          spval = 9.99e20
          IF(me == 0)THEN
            call nemsio_init(iret=status)
            call nemsio_open(nfile,trim(filename),'read',iret=status)
            if ( status /= 0 ) then
              print*,'error opening ',filename, ' Status = ', status ; stop
            endif
!---
            call nemsio_getfilehead(nfile,iret=status,nrec=nrec            &
                ,dimx=im,dimy=jm,dimz=lm,nsoil=nsoil)
            if ( status /= 0 ) then
              print*,'error finding model dimensions '; stop
            endif
            call nemsio_getheadvar(nfile,'global',global,iret)
            if (iret /= 0)then 
              print*,"global not found in file-Assigned false"
              global = .false.
            end if
            IF(modelname == 'GFS') global = .true.
! global NMMB has i=1 overlap with i=im so post will cut off i=im        
            if(global .and. modelname == 'NMM') im = im-1
 
          END IF
 
          CALL mpi_bcast(im,   1,mpi_integer,0, mpi_comm_comp,status) 
          call mpi_bcast(jm,   1,mpi_integer,0, mpi_comm_comp,status)
          call mpi_bcast(lm,   1,mpi_integer,0, mpi_comm_comp,status)
          call mpi_bcast(nsoil,1,mpi_integer,0, mpi_comm_comp,status)
 
          call mpi_bcast(global,1,mpi_logical,0,mpi_comm_comp,status)
          lp1   = lm+1
          lm1   = lm-1
          im_jm = im*jm
 
! opening GFS flux file
          IF(modelname == 'GFS') THEN
!       iunit=33
            call nemsio_open(ffile,trim(filenameflux),'read',iret=iostatusflux)
            if ( iostatusflux /= 0 ) then
              print*,'error opening ',filenameflux, ' Status = ', iostatusflux
            endif
            iostatusd3d = -1
            iunitd3d    = -1
!
! opening GFS aer file
            call nemsio_open(rfile,trim(filenameaer),'read',iret=iostatusaer)
            if ( iostatusaer /= 0  .and.  me == 0) then
              print*,'error opening AER ',filenameaer, ' Status = ', iostatusaer
            endif
!
!           print*,'iostatusD3D in WRFPOST= ',iostatusD3D
 
          END IF 
#endif
        ELSE
          print*,'UNKNOWN MODEL OUTPUT FORMAT, STOPPING'
          stop 9999
        END IF  
 
 
        CALL mpi_first()
        CALL allocate_all()
     
!
!       INITIALIZE POST COMMON BLOCKS 
!
        lcntrl = 14
        rewind(lcntrl)
 
! EXP. initialize netcdf here instead
        bbtim = mpi_wtime()
        btim = mpi_wtime()
! set default novegtype
        if(modelname == 'GFS')THEN
          novegtype = 13 
          ivegsrc   = 2
        else if(modelname=='NMM' .and. trim(ioform)=='binarynemsio')then
          novegtype = 20
          ivegsrc   = 1
        else if(modelname=='RAPR')then
          novegtype = 20
          ivegsrc   = 1
        else ! USGS
          novegtype = 24
          ivegsrc   = 0
        end if
      
! Reading model output for different models and IO format     
 
        IF(trim(ioform) == 'netcdf' .OR. trim(ioform) == 'netcdfpara') THEN
          IF(modelname == 'NCAR' .OR. modelname == 'RAPR') THEN
            IF (submodelname == 'MPAS') THEN
              CALL initpost_mpas
            ELSE
              CALL initpost
            ENDIF
          ELSE IF (modelname == 'FV3R' .OR. modelname == 'GFS') THEN
! use parallel netcdf library to read output directly
            CALL initpost_netcdf(ncid2d,ncid3d)
          ELSE
            print*,'POST does not have netcdf option for model,',modelname,' STOPPING,'
            stop 9998
          END IF
        ELSE IF(trim(ioform) == 'binarympiio') THEN 
          IF(modelname == 'NCAR' .OR. modelname == 'RAPR' .OR. modelname == 'NMM') THEN
            print*,'WRF BINARY IO FORMAT IS NO LONGER SUPPORTED, STOPPING'
            stop 9996
          ELSE IF(modelname == 'RSM') THEN                            
            print*,'MPI BINARY IO IS NOT YET INSTALLED FOR RSM, STOPPING'
            stop 9997
          ELSE
            print*,'POST does not have mpiio option for this model, STOPPING'
            stop 9998
          END IF
#if defined(BUILD_WITH_NEMSIO)
        ELSE IF(trim(ioform) == 'binarynemsio') THEN 
          IF(modelname == 'NMM') THEN
            CALL initpost_nems(nrec,nfile)
          ELSE
            print*,'POST does not have nemsio option for model,',modelname,' STOPPING,'
            stop 9998
 
          END IF
       
        ELSE IF(trim(ioform) == 'binarynemsiompiio')THEN
          IF(modelname == 'GFS') THEN
! close nemsio file for serial read
            call nemsio_close(nfile,iret=status)
            call nemsio_close(ffile,iret=status)
            call nemsio_close(rfile,iret=status)
            CALL initpost_gfs_nems_mpiio(iostatusaer)
          ELSE
            print*,'POST does not have nemsio mpi option for model,',modelname, &
            'STOPPING,'
            stop 9999
 
          END IF 
#endif
        ELSE
          print*,'UNKNOWN MODEL OUTPUT FORMAT, STOPPING'
          stop 9999
        END IF 
        initpost_tim  = initpost_tim +(mpi_wtime() - btim)
        IF(me == 0)THEN
          WRITE(6,*)'WRFPOST:  INITIALIZED POST COMMON BLOCKS'
        ENDIF
!
!       IF GRIB2 read out post aviable fields xml file and post control file
!
        if(grib == "grib2") then
          btim=mpi_wtime()
          call read_xml()
          readxml_tim = readxml_tim + (mpi_wtime() - btim)
        endif
! 
!     LOOP OVER THE OUTPUT GRID(S).  FIELD(S) AND  OUTPUT GRID(S) ARE SPECIFIED
!     IN THE CONTROL FILE.  WE PROCESS ONE GRID AND ITS FIELDS AT A TIME.
!     THAT'S WHAT THIS LOOP DOES.
!     
        icount_calmict = 0
        first_grbtbl   = .true.
        npset          = 0
!10   CONTINUE
!     
!        READ CONTROL FILE DIRECTING WHICH FIELDS ON WHICH
!        LEVELS AND TO WHICH GRID TO INTERPOLATE DATA TO.
!        VARIABLE IEOF/=0 WHEN THERE ARE NO MORE GRIDS TO PROCESS.
!
!                      --------    grib1 processing  ---------------
!                                 ------------------
!        if (grib == "grib1") then !DO NOT REVERT TO GRIB1. GRIB1 NOT SUPPORTED ANYMORE
!          IEOF = 0
!          do while (ieof == 0)
!            CALL READCNTRL(kth,IEOF)
!            IF(ME == 0)THEN
!              WRITE(6,*)'POST:  RETURN FROM READCNTRL.  ', 'IEOF=',IEOF
!            ENDIF
!
!           PROCESS SELECTED FIELDS.  FOR EACH SELECTED FIELD/LEVEL
!           WE GO THROUGH THE FOLLOWING STEPS:
!             (1) COMPUTE FIELD IF NEED BE
!             (2) WRITE FIELD TO OUTPUT FILE IN GRIB.
!
!            if (ieof == 0) then
!              CALL PROCESS(kth,kpv,th(1:kth),pv(1:kpv),iostatusD3D)
!              IF(ME == 0)THEN
!                WRITE(6,*)' '
!                WRITE(6,*)'WRFPOST:  PREPARE TO PROCESS NEXT GRID'
!              ENDIF
!            endif
!
!           PROCESS NEXT GRID.
!
!          enddo
!                      --------    grib2 processing  ---------------
!                                 ------------------
!        elseif (grib == "grib2") then
        if (me==0) write(*,*) ' in WRFPOST OUTFORM= ',grib
        if (me==0) write(*,*) '  GRIB1 IS NOT SUPPORTED ANYMORE'    
        if (grib == "grib2") then
          do while (npset < num_pset)
            npset = npset+1
            if (me==0) write(*,*)' in WRFPOST npset=',npset,' num_pset=',num_pset
            CALL set_outflds(kth,th,kpv,pv)
            if (me==0) write(*,*)' in WRFPOST size datapd',size(datapd) 
            if(allocated(datapd)) deallocate(datapd)
!Jesse x-decomposition
!           allocate(datapd(im,1:jend-jsta+1,nrecout+100))
            allocate(datapd(1:iend-ista+1,1:jend-jsta+1,nrecout+100))
!$omp parallel do private(i,j,k)
            do k=1,nrecout+100
              do j=1,jend+1-jsta
!Jesse x-decomposition
!               do i=1,im
                do i =1,iend+1-ista
                  datapd(i,j,k) = 0.
                enddo
              enddo
            enddo
            call get_postfilename(post_fname)
            if (me==0) write(*,*)'post_fname=',trim(post_fname)
            if (me==0) write(*,*)'get_postfilename,post_fname=',trim(post_fname), &
                      'npset=',npset, 'num_pset=',num_pset,            &
                      'iSF_SURFACE_PHYSICS=',isf_surface_physics
!     
!           PROCESS SELECTED FIELDS.  FOR EACH SELECTED FIELD/LEVEL
!           WE GO THROUGH THE FOLLOWING STEPS:
!             (1) COMPUTE FIELD IF NEED BE
!             (2) WRITE FIELD TO OUTPUT FILE IN GRIB.
!
            CALL process(kth,kpv,th(1:kth),pv(1:kpv),iostatusd3d)
            IF(me == 0) WRITE(6,*)'WRFPOST:  PREPARE TO PROCESS NEXT GRID'
!
!           write(*,*)'enter gribit2 before mpi_barrier'
            call mpi_barrier(mpi_comm_comp,ierr)
 
!           if(me==0)call w3tage('bf grb2  ')
!           write(*,*)'enter gribit2 after mpi barrier'
            call gribit2(post_fname)
            deallocate(datapd)
            deallocate(fld_info)
!
!           PROCESS NEXT GRID.
!
          enddo
 
        endif
!     
!-------
        call grib_info_finalize()
!
        IF(me == 0) THEN
          WRITE(6,*)' '
          WRITE(6,*)'ALL GRIDS PROCESSED.'
          WRITE(6,*)' '
        ENDIF
!
        call de_allocate
 
!       GO TO 98
 1000   CONTINUE
!exp      call ext_ncd_ioclose ( DataHandle, Status )
!
        IF(me == 0) THEN
         print*, 'INITPOST_tim = ',  initpost_tim
         print*, 'MDLFLD_tim = ',  etafld2_tim
         print*, 'MDL2P_tim =  ',eta2p_tim 
         print*, 'MDL2SIGMA_tim =  ',mdl2sigma_tim 
         print*, 'MDL2AGL_tim =  ',mdl2agl_tim 
         print*, 'SURFCE_tim =  ',surfce2_tim
         print*, 'CLDRAD_tim =  ',cldrad_tim 
         print*, 'MISCLN_tim = ',miscln_tim
         print*, 'MDL2STD_tim =  ',mdl2std_tim    
         print*, 'FIXED_tim = ',fixed_tim
         print*, 'MDL2THANDPV_tim =  ',mdl2thandpv_tim
         print*, 'CALRAD_WCLOUD_tim = ',calrad_wcloud_tim    
         print*, 'RUN_IFI_tim = ',run_ifi_tim
         print*, 'Total time = ',(mpi_wtime() - bbtim)
         print*, 'Time for OUTPUT = ',time_output
         print*, 'Time for READxml = ',readxml_tim
        endif
!     
!       END OF PROGRAM.
!
!
!       MPI_LAST WILL SHUTDOWN THE IO SERVER, IF IT EXISTS
!
        CALL mpi_last
!
!
      end if
!
!
!
!      call summary()
      if (me == 0) CALL w3tage('UNIFIED_POST')
      CALL mpi_finalize(ierr)
 
 
      stop 0
 
      END