
!------------------------------------------------------------------------!
!  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/emis/emis/LTNG_DEFN.F,v 1.6 2012/01/27 19:57:47 sjr Exp $

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

!C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      MODULE LTNG_DEFN

!C----------------------------------------------------------------------
!C Function: production of NO from lightning

!C Revision History:
!C     1 Jul 2010: Jeff Young
!C     1 Jan 2011: Rob Pinder: added support for online calculation of 
!C                 lightning NO from convective precip and NLDN data
!C    10 Mar 2011: Jeff Young: code revisions
!C    11 Apr 2011: S.Roselle: replaced I/O API include files with UTILIO_DEFN
!C    11 May 11 D.Wong: incorporated twoway model implementation
!C     6 May 13 D.Wong: replaced species name RC with RCA in the twoway model
!C    10 Jun 13 D.Wong: modified the code to work in the twoway model mode
!C    24 Sep 13 D.Wong: modified the code to allow finer met data time step
!C                      rather than fixed with one hour data
!C    Aug 15 D.Wong:    Replaced MYPE with IO_PE_INCLUSIVE for parallel I/O
!C                        implementation
!C    April 16 D.Kang: Updating the Lightning NO calculation with NLDN
!C         hourly data directly and using NLDN-RC log linear relationship,
!C         completely changing the data input process   
!C----------------------------------------------------------------------
      USE M3UTILIO

      IMPLICIT NONE

      LOGICAL, SAVE :: LTNG_NO          ! flag for lightning NO emissions
      LOGICAL, SAVE :: DIAG3D           ! flag for diagnotic files 2D or both
      INTEGER, SAVE :: LTNG_MAP         ! map to GC NO

      PUBLIC LTNG_NO, LTNG_MAP,  LTNG_INIT, GET_LTNG
      PRIVATE

!C lightning emis species name
      CHARACTER( 16 ), PARAMETER :: LTSPC = 'NO'
!C lightning log linear regression parameters with RC 
!      REAL, PARAMETER :: Intercept = -1.6
!      REAL, PARAMETER :: Slope  = 1.3
      REAL            :: SQUAREKM 
      INTEGER, PARAMETER :: SPLITCOL  = 320 
      REAL, PARAMETER :: MOLSNCG = 350
      REAL, PARAMETER :: MOLSNIC = 350

      LOGICAL, SAVE :: LPARAM     ! flag to use lightning NO parameters file
      LOGICAL, SAVE :: NLDNSTRIKE   ! true: last sync step this output tstep
      LOGICAL, SAVE :: LTNGDIAG   ! flag to turn on lightning NO diagnostics
      LOGICAL, SAVE :: LINEAR   ! flag to turn on lightning NO diagnostics
      INTEGER, SAVE :: LTLYRS     ! no. of emissions layers
      CHARACTER( 16 ), SAVE :: RC_NAME       ! RC name: old is RC and new is RCA

!C file name for off line lightning NO input; set to "InLine" for inline production
      CHARACTER( 16 ), SAVE :: LNAME
!C file name for lightning NO diagnostic output
      CHARACTER( 16 ), SAVE :: DNAME = 'LTNGOUT'
!C file name for lightning NO diagnostic output
      CHARACTER( 16 ), SAVE :: D3NAME = 'LTNGOUT3D'
!C file name for inline lightning parameters, (only used if LNAME = "InLine")
      CHARACTER( 16 ), SAVE :: ANAME
      CHARACTER( 16 ), SAVE :: ONAME
!C file name for 2dmet file for inline lightning
      CHARACTER( 16 ), SAVE :: MNAME
!C file name for 3Dmet file for inline lightning
      CHARACTER( 16 ), SAVE :: M3NAME
!C file name for ICCG parameters, (only used if LNAME = "InLine")
!C      CHARACTER( 80 ), SAVE :: iccgFile 

!C allocate these if LNAME = 'InLine'
      REAL,    ALLOCATABLE, SAVE :: COLUMN_LTNG_NO( :,: ) ! column total NO
      REAL,    ALLOCATABLE, SAVE :: VDEMIS_LT( :,:,: )   ! lightning emis
      REAL,    ALLOCATABLE, SAVE :: LTNG_PRSFC    ( :,:,:) ! surface pressure
      REAL,    ALLOCATABLE, SAVE :: LTNG_RC       ( :,:,: ) ! convective rainfall
!     REAL,    ALLOCATABLE, SAVE :: LTNG_CLDT     ( :,: ) ! cloud top
!C allocate these if LPARAM
      REAL,    ALLOCATABLE, SAVE :: STRIKE_FACTOR ( :,: ) ! flashes per unit precip
      REAL,    ALLOCATABLE, SAVE :: LRATIO        ( :,: ) ! ratio derived from NLDN data
      REAL,    ALLOCATABLE, SAVE :: NLDN_STRIKE    ( :,:,: ) ! Hourly NLDN strike data
      REAL,    ALLOCATABLE, SAVE :: OCEAN_MASK    ( :,: ) ! reduce offshore strikes
      REAL,    ALLOCATABLE, SAVE :: SLOPE    ( :,: ) ! reduce offshore strikes
      REAL,    ALLOCATABLE, SAVE :: INTERCEPT    ( :,: ) ! reduce offshore strikes
      REAL,    ALLOCATABLE, SAVE :: SLOPE_lg    ( :,: ) ! reduce offshore strikes
      REAL,    ALLOCATABLE, SAVE :: INTERCEPT_lg    ( :,: ) ! reduce offshore strikes
      REAL,    ALLOCATABLE, SAVE :: ICCG          ( :,: ) ! intercloud strikes per cloud to ground strike
!      REAL,    ALLOCATABLE, SAVE :: MOLSNFLASH    ( :,: ) ! moles N per flash
!      REAL,    ALLOCATABLE, SAVE :: MOLSNIC       ( :,: ) ! moles N per IC flash

      INTEGER, SAVE       :: NCOLS, NROWS
!C    scenario time/date needed for diagnostic output
      INTEGER, SAVE :: LOGDEV
      INTEGER :: JDATE0

!C Vertical coord values
      REAL,    ALLOCATABLE, SAVE :: VGLVSLT( : )

      INTEGER, SAVE :: GXOFF, GYOFF   ! global origin offset from file

      CONTAINS

!C======================================================================
!C Initialize lightning routines
 
         FUNCTION LTNG_INIT ( JDATE, JTIME, TSTEP, NSTEPS) RESULT ( SUCCESS )

         IMPLICIT NONE

!C Includes:
!         INCLUDE SUBST_CONST     ! constants
!         INCLUDE SUBST_FILES_ID  ! file name parameters

!C Arguments:
!        INTEGER, INTENT( IN ) :: JDATE, JTIME, TSTEP
!        INTEGER, INTENT( IN ) :: EMLAYS     ! number of emissions layers
         INTEGER :: JDATE, JTIME, TSTEP, NSTEPS

         LOGICAL SUCCESS

!C External Functions:
         INTEGER, EXTERNAL :: SETUP_LOGDEV
         INTEGER :: LOG_UNIT     ! scenario start time
         INTEGER :: COL, ROW     ! local col and row variable 

!C Local Variables:
         CHARACTER( 16 ), SAVE :: CTM_LTNG_NO = 'CTM_LTNG_NO'  ! env var flag
         CHARACTER( 16 ), SAVE :: LTNG_DIAG3D = 'DIAG3D'  ! env var flag
         CHARACTER( 16 )       :: PNAME = 'LTNG_INIT'
         CHARACTER( 80 )       :: VARDESC   ! env variable description
         CHARACTER( 120 )      :: XMSG = ' '
         real                  :: lat, lon, x, y, val
         INTEGER               :: c, r
         INTEGER               :: gridsFromLand 
         CHARACTER(80)             :: dump

         LOGICAL LTNGPARAM           ! env var to use lightning NO parameters file
         INTEGER LSPCS               ! no. of lightning species
         REAL    CONVEM              ! conversion for emissions rates to Kg/s
         INTEGER i,J, K, L, V, STATUS
         INTEGER               :: cjdate
         CHARACTER( 7 ) :: sjdate
         LOGICAL OK
!-----------------------------------------------------------------------

         SUCCESS = .TRUE.

         LOG_UNIT = INIT3()
!         LOGDEV = SETUP_LOGDEV()

            WRITE( sjdate, '(I7)') JDATE
            READ(sjdate, '(4x, I3)') cjdate
            JDATE0 = JDATE
!C Lightning NO production
!C Set LTNG_NO to Y or T to turn on lightning NO production

         VARDESC = 'Flag for lightning emissions'
         LTNG_NO = ENVYN( CTM_LTNG_NO, VARDESC, .FALSE., STATUS )
         VARDESC = 'Flag for Diagnostic file, Y for 2D only'
         DIAG3D = ENVYN( LTNG_DIAG3D, VARDESC, .FALSE., STATUS )
         IF ( STATUS .NE. 0 ) WRITE( LOGDEV, '(5X, A)' ) VARDESC
         IF ( STATUS .EQ. 1 ) THEN
            XMSG = 'Environment variable improperly formatted'
            CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
            SUCCESS = .FALSE.; RETURN
         ELSE IF ( STATUS .EQ. -1 ) THEN
            XMSG = 'Env. variable set, but empty: Using default:'
            WRITE( LOGDEV, '(5X, A, I9)' ) XMSG, JTIME
         ELSE IF ( STATUS .EQ. -2 ) THEN
            XMSG = 'Environment variable not set ... Using default:'
            WRITE( LOGDEV, '(5X, A, I9)' ) XMSG, JTIME
         END IF

            XMSG = 'Using in-line lightning NO production'
            CALL M3MSG2( XMSG )

!C Is there a file specifying iccg, ocean, and other parameters?
            VARDESC = 'Use parameter file'
            LPARAM = ENVYN( 'LTNGPARAM', VARDESC, .FALSE., STATUS )
            NLDNSTRIKE = ENVYN( 'USE_NLDN', VARDESC, .FALSE., STATUS )
            WRITE(*,*) "NLDNSTRIKE is ", NLDNSTRIKE
            LINEAR = ENVYN( 'LINEAR', VARDESC, .FALSE., STATUS )
            IF ( STATUS .NE. 0 ) WRITE( LOGDEV, '(5X, A)' ) VARDESC
            IF ( STATUS .EQ. 1 ) THEN
               XMSG = 'Environment variable improperly formatted'
               CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
               SUCCESS = .FALSE.; RETURN
            ELSE IF ( STATUS .EQ. -1 ) THEN
               XMSG = 'Env. variable set, but empty: Using default:'
               WRITE( LOGDEV, '(5X, A, I9)' ) XMSG, JTIME
            ELSE IF ( STATUS .EQ. -2 ) THEN
               XMSG = 'Environment variable not set ... Using default:'
               WRITE( LOGDEV, '(5X, A, I9)' ) XMSG, JTIME
            END IF
!C Open lightning parameters file?
               ONAME = PROMPTMFILE(                                &
                         'Enter name for PARAMETER file',         &
                         FSREAD3, 'PARAMETER_FILE', PNAME )
!C Read description of normalized emissions file
             IF(NLDNSTRIKE) then
!              WRITE(*,*) "NLDNSTRIKE is ", NLDNSTRIKE
                 ANAME = PROMPTMFILE(                                &
                       'Enter name for lightning flash file',  &
                          FSREAD3, 'FLASH_FILE', PNAME )
                 if( .not. open3(ANAME,fsread3,pname)) then
                    call m3err(pname, 0, 0,                             &
                  'Could not open Lightnign Flash FILE', .TRUE.)
                  endif
                 IF( .NOT. DIAG3D) THEN
                    IF ( .NOT. DESC3( ANAME ) ) THEN
                      XMSG = 'Could not get description of file "'      &
                         // TRIM( ANAME ) // '"'
                      CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
                      SUCCESS = .FALSE.; RETURN
                     END IF
                 ENDIF
              ELSE
                RC_NAME = 'RC'
              ENDIF                  
             if(DIAG3D) then
                MNAME = PROMPTMFILE(                                   & 
                          'Enter name for gridded met input file', &
                         FSREAD3, 'MET_CRO_2D', PNAME )
                 if( .not. open3(MNAME,fsread3,pname)) then
                     call m3err(pname, 0, 0,                             &
                     'Could not open MET_CRO_2D', .TRUE.)
                 endif

                 M3NAME = PROMPTMFILE(                                   & 
                          'Enter name for gridded met input file', &
                         FSREAD3, 'MET_CRO_3D', PNAME )
                 if( .not. open3(M3NAME,fsread3,pname)) then
                     call m3err(pname, 0, 0,                             &
                     'Could not open MET_CRO_3D', .TRUE.)
                 endif
                 IF ( .NOT. DESC3( M3NAME ) ) THEN
                    XMSG = 'Could not get description of file "'      &
                         // TRIM( M3NAME ) // '"'
                    CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
                    SUCCESS = .FALSE.; RETURN
                 END IF
             ENDIF 

              WRITE(*,*) "MNAME: ", MNAME
              WRITE(*,*) "NLAYS3D ", NLAYS3D 
!              WRITE(*,*) "NLDNSTRIKE = ", NLDNSTRIKE

           NCOLS = NCOLS3D
           NROWS = NROWS3D 
           LTLYRS = NLAYS3D

            ALLOCATE( COLUMN_LTNG_NO( NCOLS,NROWS ),                   &
                     LTNG_PRSFC    ( NCOLS,NROWS,NSTEPS ),                    &
                     VGLVSLT       (0:LTLYRS),                          &
                     VDEMIS_LT( NCOLS,NROWS,LTLYRS ),                   &
                     LTNG_RC       ( NCOLS,NROWS,NSTEPS ), STAT = STATUS )
            ALLOCATE( NLDN_STRIKE( NCOLS,NROWS,NSTEPS ))
!    &                LTNG_CLDT     ( NCOLS,NROWS ),
            IF ( STATUS .NE. 0 ) THEN
               XMSG = 'COLUMN_LTNG_NO, LTNG_PRSFC, LTNG_RC, or LTNG_CLDT'  &
                   // '  memory allocation failed'
               CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
               SUCCESS = .FALSE.; RETURN
            END IF
               WRITE(*,*) "LTLYRS = ", LTLYRS
!              WRITE(*,*) "LPARM = ", LPARAM

            VGLVSLT = VGLVS3D(1:(LTLYRS+1))
            WRITE(*,*) "VGLVS = ", VGLVSLT
            IF ( LPARAM ) THEN
               XMSG = 'Using lightning NO parameters file'
               CALL M3MSG2( XMSG )

!C...  set Lambert projection
      IF( .NOT. SETLAM( SNGL( P_ALP3D ),                               &
           SNGL( P_BET3D ),                                            &
           SNGL( P_GAM3D ),                                            &
           SNGL( XCENT3D ),                                            &
          SNGL( YCENT3D ) ) ) THEN
         CALL M3EXIT( 'CELL_LL', 0, 0,                                 &
             'Lambert projection setup error', 2 )
      ENDIF
      write(*,*) P_BET3D, P_GAM3D, XCENT3D, YCENT3D, P_ALP3D
      SQUAREKM = (XCELL3D * 1.0e-3) * (YCELL3D * 1.0e-3)
!C Check grid definition (initialize if first call)
              ALLOCATE( LRATIO       ( NCOLS,NROWS ),                 &
                        ICCG         ( NCOLS,NROWS ),                 &
                        OCEAN_MASK   ( NCOLS,NROWS ),                 &
                        SLOPE        ( NCOLS,NROWS ),                 &
                        INTERCEPT    ( NCOLS,NROWS ),                 &
                        SLOPE_lg        ( NCOLS,NROWS ),              &
                        INTERCEPT_lg    ( NCOLS,NROWS ),              &
                        STRIKE_FACTOR( NCOLS,NROWS ), STAT = STATUS )
               IF ( STATUS .NE. 0 ) THEN
                  XMSG = 'LRATIO, ICCG, OCEAN_MASK, MOLSNFLASH or STRIKE_FACTOR'  &
                      // '  memory allocation failed'
                  CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
                  SUCCESS = .FALSE.; RETURN
               END IF
!               OCEAN_MASK = 1.0
             DO i = 1, NSTEPS
               IF(DIAG3D) THEN
                  IF ( .NOT. READ3( MNAME, 'PRSFC', 1,   &
                             JDATE, JTIME, LTNG_PRSFC(:,:,i) ) ) THEN
                    XMSG = 'Could not read PRSFC from ' // TRIM( MNAME )
                    CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT2 )     
                   END IF
                ENDIF
              IF(NLDNSTRIKE) THEN 
                 IF ( .NOT. READ3( ANAME, "LNT", 1,               &
!                 IF ( .NOT. READ3( ANAME, "NLDNstrk", 1,               &
                                  JDATE, JTIME, NLDN_STRIKE(:,:,i) ) ) THEN
                   XMSG = 'Could not read LNT from' // TRIM( ANAME )
                 END IF
              ELSE
                 IF ( .NOT. READ3( MNAME, RC_NAME, 1,   &
                               JDATE, JTIME, LTNG_RC(:,:,i) ) ) THEN
                   XMSG = 'Could not read RC from ' // TRIM( MNAME )
                   CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT2 )     
               END IF
              ENDIF
               CALL NEXTIME( JDATE, JTIME, TSTEP)
            ENDDO
              if( .not. open3(ONAME,fsread3,pname)) then
               call m3err(pname, jdate, jtime,                               &
              'Could not open Ocean Mask FILE', .TRUE.)
              endif
               IF ( .NOT. READ3( ONAME, "OCNMASK", 1,                &
                                  jdate, jtime, OCEAN_MASK ) ) THEN
                  XMSG = 'Could not read OCNMASK from' // TRIM( ONAME ) // '"'
               END IF
               IF ( .NOT. READ3( ONAME, "SLOPE", 1,                  &
                                  jdate, jtime, SLOPE ) ) THEN
                  XMSG = 'Could not read SLOPE from' // TRIM( ONAME ) // '"'
               END IF
!               IF ( .NOT. READ3( ONAME, "INTERCEPT", 1,              &
!                                  jdate, jtime, INTERCEPT ) ) THEN
!                  XMSG = 'CCuld not read INTERCEPT from' // TRIM( ONAME ) // '"'
!               END IF
               IF ( .NOT. READ3( ONAME, "SLOPE_lg", 1,               &
                                  jdate, jtime, SLOPE_lg ) ) THEN
                  XMSG = 'Could not read SLOPE_lg from' // TRIM( ONAME ) // '"'
               END IF
!               IF ( .NOT. READ3( ONAME, "INTERCEPT_lg",             &
!                                  jdate, jtime, INTERCEPT_lg ) ) THEN
!                  XMSG = 'Could not read INTERCEPT_lg from' // TRIM( ONAME ) // '"'
!               END IF

!               IF ( .NOT. XTRACT3( ONAME, "INTERCEPT_lg", 1, 1,                &
!                                   1, 299, 1, 459,    &
!                                   0, 0, INTERCEPT_lg ) ) Then
!                  XMSG = 'Could not interpolate OPEN from ' // ONAME
!                  Call M3EXIT( PNAME, 0, 0, XMSG, XSTAT1 )
!               End If

               IF(cjdate > 90 .and. cjdate < 274) THEN 
                  IF ( .NOT. READ3( ONAME, "ICCG_SUM", 1,          &
                                  jdate, jtime, ICCG ) ) THEN
                    XMSG = 'Could not read ICCG from' // TRIM( ONAME ) // '"'
                  END IF
                ELSE
                  IF ( .NOT. READ3( ONAME, "ICCG_WIN", 1,          &
                                  jdate, jtime, ICCG ) ) THEN
                    XMSG = 'Could not read ICCG from' // TRIM( ONAME ) // '"'
                  END IF
                ENDIF
              
            END IF   ! LPARAM

!C Lightning Diagnostic file?
            LTNGDIAG = ENVYN( 'LTNGDIAG',                     &
                   'Write lightning NO production diagnostic file?',  &
                   .FALSE., STATUS )
           
            IF ( LTNGDIAG ) THEN
!C (all but variables-table and horizontal domain in description is borrowed from MNAME)
               SDATE3D = JDATE0
               NVARS3D = 1 
               VNAME3D( 1 ) = LTSPC
               VTYPE3D( 1 ) = M3REAL
               UNITS3D( 1 ) = 'moles/s'

               FDESC3D = ' '   ! array assignment
               FDESC3D( 1 ) = 'Gridded lightning NO production from CMAQ'
               FDESC3D( 2 ) = '/from/ ' // PNAME
               FDESC3D( 3 ) = '/Version/ CMAQ'
            IF(DIAG3D) THEN
                VDESC3D( 1 ) = 'hourly 3D NO produced from inline lightning'
                IF ( .NOT. OPEN3( D3NAME, FSUNKN3, PNAME ) ) THEN
                     XMSG = 'Could not open ' // TRIM(  D3NAME )
                     CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
                END IF
           ENDIF
               NLAYS3D = 1
               VDESC3D( 1 ) = 'hourly average NO produced from NLDN lightning'
!C Open output file (mol/s)
                  IF ( .NOT. OPEN3( DNAME, FSNEW3, PNAME ) ) THEN
                     XMSG = 'Could not open ' // TRIM( DNAME )
                     CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
                  END IF

            END IF   ! LTNGDIAG

         END FUNCTION LTNG_INIT

!C======================================================================
!C Get NO produced from lightning in VDEMIS_LT

         SUBROUTINE GET_LTNG ( JDATE, JTIME, TSTEP )

         IMPLICIT NONE

!     need to change this to set this path in the Makefile
!        include "/usr/local/intel/mkl/include/mkl_vml.f77"
!         INCLUDE SUBST_CONST     ! constants
!         INCLUDE SUBST_FILES_ID  ! file name parameters

         INTEGER JDATE, JTIME, TSTEP
         INTEGER INDX

         REAL,   PARAMETER :: CONVPA = 1.0E-2  ! convert Pa to hPa
!        REAL,   PARAMETER :: WK = 8.0         ! shape parameter for weibull distribution
!        REAL,   PARAMETER :: WLAMBDA = 700.0  ! scale parameter for weibull distribution
         REAL,   PARAMETER :: WMU = 350.0      ! mean
         REAL,   PARAMETER :: WSIGMA = 200.0   ! standard deviation
         REAL,   PARAMETER :: W2MU = 600.0     ! mean
         REAL,   PARAMETER :: W2SIGMA = 50.0   ! standard deviation
         REAL,   PARAMETER :: SQRT2 = 1.414213562731
         REAL,   PARAMETER :: SFACTOR1 = 0.95  ! the scaling factor for the wider distribution (WMU350)
         REAL,   PARAMETER :: SFACTOR2 = 0.12  ! the scaling factor for the wider distribution (WMU350

         INTEGER COL, ROW, LAY ! iterator variables

         REAL    PCALC    ! pressure level for NO vertical distribution (hPa)
         REAL    BOTTOM   ! pressure at bottom of grid cell (hPa)
         REAL    TOP      ! pressure at top of grid cell (hPa)        
         REAL    BOTTOM_FRAC, TOP_FRAC ! their difference is the fraction of lightning NO in this grid cell
         REAL    BOTTOM_FRAC2, TOP_FRAC2
         REAL    SUM_FRAC ! stores the sum of vertical fractions to re-normalize the column
         REAL    WEIGHT ! used to normalize emissions to total amount
         REAL    inErfB, inErfT !  nputs to error funciton calculation
         REAL    outErfB, outErfT ! outputs from error funciton calculation
         REAL :: LTEMIS( LTLYRS )
         REAL :: LTEMIS2( LTLYRS )
         REAL    XCELLR, YCELLR   ! cell spacing ratio to 36Km
         REAL    FLASH_FAC        ! lightning flashes factor

         LOGICAL, SAVE :: LASTTIC   ! true: last sync step this output tstep
         REAL    DIVFAC   ! averaging factor for diagnostic file

         CHARACTER( 16 ) :: MNAME
         CHARACTER( 16 ) :: PNAME = 'GET_LTNG'
         CHARACTER( 120 ) :: XMSG = ' '

         INTEGER               :: SPC
         INTEGER               :: CDATE, CTIME
         INTEGER               :: cjdate 
         CHARACTER( 7 ) :: sjdate 

         LOGICAL, SAVE :: FIRSTIME = .TRUE.

!-----------------------------------------------------------------------

!         IF ( LNAME .EQ. "InLine" ) THEN
!C case of inline lightning NO production


            WRITE( LOGDEV,'(5X, A)' ) 'InLine LNOx production'
            WRITE( sjdate, '(I7)') JDATE
            READ(sjdate, '(4x, I3)') cjdate 
!C Open met file

!C Get domain window info for met_cro_2d file

!           IF ( .NOT. INTERPX( MNAME, "RC", PNAME,  ! convective rain in cm

!C Iterate over each grid cell and distribute lightning NO vertically
            COLUMN_LTNG_NO( :,: ) = 0.0
            IF ( LPARAM ) THEN
                INDX = JTIME/10000 + 1
               DO ROW = 1, NROWS 
                  DO COL = 1, NCOLS
                    IF ( NLDNSTRIKE ) THEN
                     COLUMN_LTNG_NO( COL,ROW ) =                &
                        (NLDN_STRIKE(COL,ROW, INDX)             &
                         *SQUAREKM                              &
                         * OCEAN_MASK( COL,ROW )                &
                         *(MOLSNCG + MOLSNIC * ICCG( COL,ROW )))        &
                       / ( 60.0 * 60.0 )         ! get time units right:emisions are in the unit of moles/s
                   ELSE 
                     IF(LTNG_RC(COL, ROW, INDX) > 0) THEN
!                        SCL_FACTOR = 0.8   !min(0.2, EXP(INTERCEPT_lg(COL, ROW)))
                        
                        IF( .NOT. LINEAR .AND. OCEAN_MASK(COL,ROW).GT. 0.2)  THEN
                            INTERCEPT_LG(COL,ROW) = 0.0
                            COLUMN_LTNG_NO( COL,ROW ) =               &
                            ( EXP( SLOPE_LG( COL,ROW ) * LOG( LTNG_RC(COL,ROW, INDX))  &
                            + INTERCEPT_LG( COL,ROW) )                &
                            * SQUAREKM                                    &
                            * OCEAN_MASK( COL,ROW )            &
                            * ( MOLSNCG + (MOLSNIC * ICCG( COL,ROW ) )))  &
                            / ( 60.0 * 60.0 )         ! get time units right:emisions are in the unit of moles/s
!                        WRITE(*,*) "Intercept_lt = ", INTERCEPT_LG( COL,ROW) 
                        ELSE 
                           COLUMN_LTNG_NO( COL,ROW ) =                    &
                           ( ( SLOPE( COL,ROW ) * LTNG_RC( COL,ROW, INDX ) + INTERCEPT( COL,ROW ) )  &
                           * SQUAREKM  &
                           * OCEAN_MASK( COL,ROW )  &
                           * ( MOLSNCG + ( MOLSNIC * ICCG( COL,ROW ) ) ) )  &
                           / ( 60.0 * 60.0 )      
                           IF ( COLUMN_LTNG_NO( COL,ROW ) .LT. 0 ) COLUMN_LTNG_NO( COL,ROW ) = 0.0
                      END IF
                    ELSE
                       COLUMN_LTNG_NO( COL,ROW ) = 0.0
                    ENDIF
                  ENDIF

                  END DO
               END DO

            END IF

           VDEMIS_LT = 0.0   ! array assignment

            DO ROW = 1, NROWS
               DO COL = 1, NCOLS

!C check to see if there are lightning strikes for this grid cell
!C only calculate lightning for cloud top greater than 6500 meters
!                 write( logdev,* ) col, row, column_ltng_no( col,row ),
!    &                              ltng_cldt( col,row )
!                 IF ( COLUMN_LTNG_NO( COL,ROW ) .GT. 0.0  .AND.
!    &                 LTNG_CLDT( COL,ROW ) .GT. 6500.0 )  THEN

                  IF ( COLUMN_LTNG_NO( COL,ROW ) .LE. 0.0 ) CYCLE
                  SUM_FRAC = 0.0
                  LTEMIS = 0.0   ! array assignment

                  DO LAY = 1, LTLYRS 
!                    WRITE(*,*) "LAY = ", LAY, " VGLVSLT = ", VGLVSLT( LAY-1 )
!C Get pressures: Use SIGMA values and surface pres.
!p=sigma*(psfc-ptop)+ptop
                     BOTTOM = ( VGLVSLT( LAY-1 )                     &
                        * ( LTNG_PRSFC( COL,ROW,INDX ) - VGTOP3D )  &
                        + VGTOP3D ) * CONVPA
!                           write( logdev,* ) "bottom: ", bottom
                     TOP    = ( VGLVSLT( LAY )                       &
                        * ( LTNG_PRSFC( COL,ROW,INDX ) - VGTOP3D )  &
                           + VGTOP3D ) * CONVPA
!                           write( logdev,* ) "top: ", top

!C Find the bottom and top of each layer, and calculate the fraction
!C of the column emissions for that layer
!C Use normal distribution, mean = wmu, standard deviation = wsigma
                     inErfB = ( BOTTOM - WMU ) / ( WSIGMA * SQRT2 )
                     inErfT = ( TOP - WMU ) / ( WSIGMA * SQRT2 )
                     outErfB = ERF( inErfB )
                     outErfT = ERF( inErfT )
                     BOTTOM_FRAC = 0.5 * ( 1.0 + outErfB )
                     TOP_FRAC    = 0.5 * ( 1.0 + outErfT )

!C Find the bottom and top of each layer, and calculate the fraction
!C of the column emissions for that layer
!C use normal distribution, mean = wmu, standard deviation = wsigma
                     inErfB = ( BOTTOM - W2MU ) / ( W2SIGMA * SQRT2 )
                     inErfT = ( TOP - W2MU ) / ( W2SIGMA * SQRT2 )
                     outErfB = ERF( inErfB )
                     outErfT = ERF( inErfT )
                     BOTTOM_FRAC2 = 0.5 * ( 1.0 + outErfB )
                     TOP_FRAC2    = 0.5 * ( 1.0 + outErfT )

!C Add weighted contribution to this level
                     WEIGHT = ( BOTTOM_FRAC - TOP_FRAC ) * SFACTOR1         &
                           + ( BOTTOM_FRAC2 - TOP_FRAC2 ) * SFACTOR2

                     LTEMIS( LAY ) = WEIGHT * COLUMN_LTNG_NO( COL,ROW )

!C Sum weights in order to normalize to 1
                     SUM_FRAC = SUM_FRAC + WEIGHT
!                    write( logdev,* ) col, row, lay,
!     &                                bottom_frac, top_frac,
!     &                                bottom_frac2, top_frac2,
!     sum_frac,
!     &                                column_ltng_no( col,row )

!C If emissions are less than 0, generate an error message in the log
                     IF ( LTEMIS( LAY ) .LT. 0.0 ) THEN
                        WRITE(*,*) "VGTOP = ",VGTOP3D, " LSF = ",LTNG_PRSFC( COL,ROW,INDX) 
                        WRITE(*,*) "Bf = ", BOTTOM_FRAC, " Tf = ", TOP_FRAC
                        WRITE(*,*) "B = ", BOTTOM, " T = ", TOP
                        WRITE(*,*) "W = ", WEIGHT, " CLNO = ", COLUMN_LTNG_NO( COL,ROW )
                        WRITE(*,*) "C = ", COL, " R = ", ROW, " L = ",LAY, "LTEMIS = ", LTEMIS(LAY)
                        WRITE( LOGDEV,* ) LTEMIS( LAY ),          &
                                COLUMN_LTNG_NO( COL,ROW ),        &
                                         BOTTOM_FRAC, TOP_FRAC
                        XMSG = '*** Ltng NO emis is less than zero'
                        CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT2 )
                     END IF

!                    write( logdev,* ) "VDEMIS_LT ", VDEMIS_LT(
!                    1,COL,ROW,LAY )

                  END DO   ! end layers loop

                  DO LAY = 1, LTLYRS 
!C Re-normalize, in some cases area under the error function is not 1
                     VDEMIS_LT( COL,ROW,LAY ) = LTEMIS( LAY ) / SUM_FRAC
!                    write( logdev,* ) col, row, lay,
!     &                                vdemis_lt( col,row,lay )
                  END DO        ! layers renormalized

               END DO   ! columns
            END DO   ! rows

                 
!C Write lightning NO to the diagnostic file
            WRITE(*,*) "seg here???"
            WRITE(*,*) "LTNGD = ", LTNGDIAG 

            IF ( LTNGDIAG ) THEN
                  IF ( .NOT. WRITE3( DNAME, LTSPC, JDATE, JTIME, COLUMN_LTNG_NO ) )  THEN
                     XMSG = 'Could not write to ' // TRIM( DNAME )
                     CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT2 )
                  ELSE
                     WRITE( LOGDEV,94040 )               &
                         'Timestep written to', TRIM( DNAME ),   &
                         'for date and time', JDATE, JTIME
                  END IF
              IF ( .NOT. WRITE3( D3NAME, LTSPC, JDATE, JTIME, VDEMIS_LT ) )  THEN
                     XMSG = 'Could not write to ' // TRIM( D3NAME )
                     CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT2 )
                  ELSE
                     WRITE( LOGDEV,94040 )                     &
                      'Timestep written to', TRIM(D3NAME ),    &
                         'for date and time', JDATE, JTIME
                  END IF


            END IF  ! diagnostics turned on

         RETURN

!C------------------  Format  Statements   ------------------------------

94040 FORMAT( /5X, 3( A, :, 1X ), I8, ":", I6.6 )
94042 FORMAT( /5X, A, 1X, I8, ":", I6.6, 1X, 1PE13.5 )
        
         END SUBROUTINE GET_LTNG
      
      END MODULE LTNG_DEFN

