!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!

C RCS file, release, date & time of last delta, author, state, [and locker]
C $Header: /project/yoj/arc/CCTM/src/init/yamo/opGSI.F,v 1.4 2011/10/21 16:11:24 yoj Exp $ 

C what(1) key, module and SID; SCCS file; date and time of last delta:
C %W% %P% %G% %U%

C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE OPGSI (ccount, CGRID, JDATE, JTIME, TSTEP )

C-----------------------------------------------------------------------
C Function:
C   Create the IO/API netCDF header and open the output CONC file

C Revision history:
C   Jeff - Dec 00 - split out from initscen.F
C                 - move CGRID_MAP into f90 module
C   Jeff - Feb 01 - assumed shape arrays
C   30 Mar 01 J.Young: dyn alloc - Use HGRD_DEFN
C
C    3 Sep 01 David Wong
C     -- let PE 0 open CTM_GSI_1 as new and later on let the rest open
C        it for read and write
C     -- put an explicit barrier before opening a new netCDF file to avoid
C        NCOPEN error
C    7 May 03 J.Young: open and close conc file in processor 0; spin-wait to
C                      deal with nfs network latency for mpich cluster
C   28 Aug 03 J.Young: following Zion Wang at CERT, remove the spin-wait and
C                      simplify opening and closing CONC file (see initscen)
C   31 Jan 05 J.Young: dyn alloc - establish both horizontal & vertical
C                      domain specifications in one module
C   30 May 05 J.Young: optional save derived vert. vel. to conc file
C   21 Jan 06 J.Young: add subset non-reactives capability with NR_CONC.EXT
C                      and subset layers saved to conc file; no longer save
C                      RHOJ to conc file
C   22 Aug 06 J.Young: enable capturing ICs as 0 timestep
C   11 Jun 08 J.Young: clean up
C   30 Mar 10 David Wong: removed unnecessary barrier
C   16 Feb 11 S.Roselle: replaced I/O API include files with UTILIO_DEFN;
C                      removed deprecated TRIMLEN
C-----------------------------------------------------------------------

      USE GRID_CONF             ! horizontal & vertical domain specifications
      USE WVEL_DEFN             ! derived vertical velocity component
      USE STD_CONC              ! standard CONC
      USE UTILIO_DEFN

      USE SE_MODULES         ! stenex
!     USE SUBST_UTIL_MODULE     ! stenex

      IMPLICIT NONE

      INCLUDE SUBST_FILES_ID    ! file name parameters

C Arguments:
      character *(*) ccount

      REAL, POINTER :: CGRID( :,:,:,: )  ! for initial CONC
      INTEGER      JDATE        ! starting date (YYYYDDD)
      INTEGER      JTIME        ! starting time (HHMMSS)
      INTEGER      TSTEP        ! output timestep (HHMMSS)

C Local Variables:

      REAL, ALLOCATABLE :: DBUFF( :,:,: )
      INTEGER      ALLOCSTAT
      INTEGER      STATUS               !  ENVINT status

      CHARACTER( 16 ) :: PNAME = 'OPGSI'
      CHARACTER( 96 ) :: XMSG = ' '
      CHARACTER( 28 ) :: SSTR = ' species saved to GSI file:'

C environment variable for no. of layers from bottom to save on CONC file
      CHARACTER( 16 ) :: NLAYS_CONC = 'NLAYS_CONC'
      character(16) CTMGSI

C  environment variable description
      CHARACTER( 80 ) :: VARDESC

!JDE      INTEGER      LOGDEV

      INTEGER      K, KD, L, SPC, V   ! loop counters
      INTEGER      STRT, FINI         ! loop counters
      INTEGER      INDX
      INTEGER      INO2,VNO2 ! pointer to NO2 information

C-----------------------------------------------------------------------

! JDE add
      CALL SUBST_BARRIER
! end JDE add
! JDE      LOGDEV = INIT3()

C Set output file characteristics based on COORD.EXT and open it
      write(logdev,*)'top of opGSI'
      write(logdev,*)'JDE NLAYS = ',NLAYS
      call flush(logdev)
      CTMGSI='CTM_GSI_'//trim(ccount)

      FTYPE3D = GRDDED3
      SDATE3D = JDATE
      STIME3D = JTIME
      TSTEP3D = TSTEP
!      NVARS3D = N_CSPCS
      !IF ( W_VEL ) NVARS3D = NVARS3D + 1   ! for W_YAMO
      NVARS3D = 1
      NCOLS3D = GL_NCOLS
      NROWS3D = GL_NROWS
!     NLAYS3D = C_NLAYS !JDE 5/2020, use NLAYS from GRID_CONF so NLAYS doesn't
!     rely on conc file layers, which could be 1
      NLAYS3D = NLAYS
      NTHIK3D = 1
      GDTYP3D = GDTYP_GD
      P_ALP3D = P_ALP_GD
      P_BET3D = P_BET_GD 
      P_GAM3D = P_GAM_GD
      XORIG3D = XORIG_GD
      YORIG3D = YORIG_GD
      XCENT3D = XCENT_GD
      YCENT3D = YCENT_GD
      XCELL3D = XCELL_GD
      YCELL3D = YCELL_GD
      VGTYP3D = VGTYP_GD
      VGTOP3D = VGTOP_GD
!     VGTPUN3D = VGTPUN_GD ! currently, not defined
      DO L = 1, NLAYS3D + 1
         VGLVS3D( L ) = VGLVS_GD( L )
      END DO
      GDNAM3D = GRID_NAME  ! from HGRD_DEFN

      FDESC3D( 1 ) = 'Concentration file output'
      FDESC3D( 2 ) = 'From CMAQ model dyn alloc version CTM'
      FDESC3D( 3 ) = 'Set of variables just NO2'
      FDESC3D( 4 ) = 'For next scenario continuation runs,'
      FDESC3D( 5 ) = 'use the "one-step" CGRID file'
      KD = 5
      !V = CONC_BLEV ! JDE 5/2020
      V = 1
      L = 0
      !DO K = KD + 1, MIN ( C_NLAYS + KD, MXDESC3 ) !JDE change to NLAYS
      DO K = KD + 1, MIN ( NLAYS3D + KD, MXDESC3 )
         L = L + 1
         WRITE( FDESC3D( K ),'( "Layer", I3, " to", I3, " " )' )
     &   V + L - 1, L
      END DO
      IF ( ( KD + 1 + L ) .LT. MXDESC3 ) THEN
         DO K = KD + 1 + L, MXDESC3
            FDESC3D( K ) = ' '
         END DO
      END IF

      WRITE( LOGDEV,* ) ' '
      WRITE( LOGDEV,* ) '      GSI File Header Description:'
      DO K = 1, KD + L
         WRITE( LOGDEV,* ) '    => ', TRIM( FDESC3D( K ) )
      END DO
      V = 0
      FINI = N_C_GC_SPC
      ino2=-1
      vno2=-1
      strt=1
      DO SPC = STRT, FINI
         V = V + 1
         if( C_GC_SPC( V ).eq.'NO2')then
           INO2=SPC
           VNO2=V
           write(LOGDEV,*)'opGSI INO2',INO2,'VNO2',VNO2
           WRITE(LOGDEV,*)'JDE check C_GC_SPC( V ) ',TRIM(C_GC_SPC( V )),' VNAME3D( SPC ) ',TRIM(VNAME3D( SPC ))
           call flush(6)
           exit
         endif
      END DO
      if (INO2<1)then
            XMSG = 'Could not find NO2'
            CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
      endif
      spc=ino2
      v=vno2
      ! JDE changed SPC to 1 to fix file header
      !VTYPE3D( SPC ) = M3REAL
      VTYPE3D( 1 ) = M3REAL
      !WRITE(LOGDEV,*)'JDE0 VNAME3D: ',trim(VNAME3D( SPC ))
      !VNAME3D( SPC ) = C_GC_SPC( V )
      !WRITE(LOGDEV,*)'JDE VNAME3D array: ',VNAME3D(:)
      VNAME3D( 1 ) = C_GC_SPC( V )
      !WRITE(LOGDEV,*)'JDE0 C_GC_SPC: ', trim(C_GC_SPC( V ))
      !WRITE(LOGDEV,*)'JDE0 VNAME3D: ',trim(VNAME3D( SPC ))
      !UNITS3D( SPC ) = 'ppmV'
      UNITS3D( 1 ) = 'ppmV'
      !VDESC3D( SPC ) = 'Variable ' // VNAME3D( SPC )
      !VDESC3D( 1 ) = 'Variable ' // VNAME3D( SPC )
      VDESC3D( 1 ) = 'Variable ' // VNAME3D( 1 )
#if 0
! JDE if 0 this block does not get compiled
      V = 0
      STRT = 1
      FINI = N_C_GC_SPC
      FINI = 1
!     do NO2 only
      DO SPC = STRT, FINI
         V = V + 1
         VTYPE3D( SPC ) = M3REAL
         VNAME3D( SPC ) = C_GC_SPC( V )
         UNITS3D( SPC ) = 'ppmV'
         VDESC3D( SPC ) = 'Variable ' // VNAME3D( SPC )
      END DO
#endif


C create header
      write(logdev,*)'before open ',trim(CTMGSI)
      call flush(6)
      IF ( MYPE .EQ. 0 ) THEN   ! open new
         write(logdev,*)'CTM_GSI_1 ',trim(CTMGSI)
         WRITE(LOGDEV,*)'JDE VNAME3D(1): ',trim(VNAME3D( 1 ))
         WRITE(LOGDEV,*)'JDE C_GC_SPC(V): ', trim(C_GC_SPC( V ))
         call flush(logdev)
         !IF ( .NOT. OPEN3( CTMGSI, FSCREA3, PNAME ) ) THEN
         IF ( .NOT. OPEN3( CTMGSI, FSNEW3, PNAME ) ) THEN
            XMSG = 'opGSI Could not open ' // CTMGSI // '. File may already exist.' 
            CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
         END IF
      END IF
      CALL SUBST_BARRIER
      if ( .not. io_pe_inclusive ) then !JDE
        IF ( .NOT. OPEN3( CTMGSI, FSREAD3, PNAME ) ) THEN
           XMSG = 'opGSI Could not open ' // CTMGSI
           CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
        END IF
      end if !JDE if PE0
      write(logdev,*)'opened CMT_GSI_ file'
      call flush(logdev)
      CALL SUBST_BARRIER

C write the initial concentrations as step 0 on the conc file
C (inital data assumed to be in correct output units)

      ALLOCATE ( DBUFF( NCOLS,NROWS,NLAYS ), STAT = ALLOCSTAT )
      !ALLOCATE ( DBUFF( NCOLS,NROWS,C_NLAYS ), STAT = ALLOCSTAT ) !JDE 5/2020 NLAYS
      !JDE ALLOCATE ( DBUFF( MY_NCOLS,MY_NROWS,C_NLAYS ), STAT = ALLOCSTAT )
      IF ( ALLOCSTAT .NE. 0 ) THEN
         XMSG = 'Failure allocating DBUFF'
         CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
      END IF

      WRITE( LOGDEV,* ) ' '
      INDX = CONC_MAP( V )

      DBUFF = CGRID( 1:NCOLS,1:NROWS,1:NLAYS,INDX )
      !DBUFF = CGRID( 1:NCOLS,1:NROWS,CONC_BLEV:CONC_ELEV,INDX ) !JDE 5/2020
      !JDE DBUFF = CGRID( 1:MY_NCOLS,1:MY_NROWS,CONC_BLEV:CONC_ELEV,INDX )
      write(logdev,*)'write NO2'
      write(logdev,*)'JDE2 C_GC_SPC: ', trim(C_GC_SPC( V ))
      write(logdev,*)'JDE2 C_GC_SPC: ', trim(C_GC_SPC( SPC ))
      call flush(logdev)

!      IF ( .NOT. WRITE3( CTMGSI, C_GC_SPC( SPC ),
      IF ( .NOT. WRITE3( CTMGSI, C_GC_SPC( V ),
     &                    JDATE, JTIME, DBUFF ) ) THEN
          XMSG = 'opGSI Could not write ' //
     &           TRIM( C_GC_SPC( SPC ) ) //
     &           ' to ' // CTMGSI
          CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
      END IF

      WRITE( LOGDEV,'( 5X, I4, " (", I3, ") ", A )' )
     &                SPC, V, TRIM( C_GC_SPC( SPC ) )

#if 0
      V = 0


!      DO SPC = 1, N_C_GC_SPC
!     NO2 only
      DO SPC = 1,1 
         V = V + 1
         INDX = CONC_MAP( V )

         DBUFF = CGRID( 1:NCOLS,1:NROWS,1:NLAYS,INDX )
         !DBUFF = CGRID( 1:NCOLS,1:NROWS,CONC_BLEV:CONC_ELEV,INDX ) !JDE 5/2020
         !JDE DBUFF = CGRID( 1:MY_NCOLS,1:MY_NROWS,CONC_BLEV:CONC_ELEV,INDX )
         write(logdev,*)'write NO2 ',C_GC_SPC(SPC),' V ',V, 'INDX',INDX
         call flush(logdev)

         IF ( .NOT. WRITE3( CTMGSI, C_GC_SPC( SPC ),
     &                      JDATE, JTIME, DBUFF ) ) THEN
            XMSG = 'Could not write ' //
     &             TRIM( C_GC_SPC( SPC ) ) //
     &             ' to ' // CTMGSI
            CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
         END IF

         WRITE( LOGDEV,'( 5X, I4, " (", I3, ") ", A )' )
     &                SPC, V, TRIM( C_GC_SPC( SPC ) )

      END DO
#endif
      write(logdev,*)'did write3'
      call flush(logdev)


      DEALLOCATE ( DBUFF )

      WRITE( LOGDEV, '( /5X, 3( A, :, 1X ), I8, ":", I6.6 )' )
     &  'Timestep written to', CTMGSI,
     &  'for date and time', JDATE, JTIME
      WRITE( LOGDEV, '(  5X,  A, 1X, I8, ":", I6.6 )' )
     &  'from timestep on initial data files for date and time',
     &   JDATE, JTIME
      !IF ( MYPE .EQ. 0 ) THEN   ! open new
      IF ( .NOT. CLOSE3( CTMGSI ) ) THEN 
         XMSG = 'opgsi Could not close file "' // TRIM( CTMGSI ) // '"'
         CALL M3EXIT( PNAME, 0, 0, XMSG, XSTAT2 )
      END IF
      !END IF

      RETURN
      END
