!$Id: mknoprod.F90 157 2010-01-18 14:21:59Z acosce $
!! =========================================================================
!! INCA - INteraction with Chemistry and Aerosols
!! 
!! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE)
!! 	     Unite mixte CEA-CNRS-UVSQ
!! 
!! Contributors to this INCA subroutine:
!! 
!! Didier Hauglustaine, LSCE, hauglustaine@cea.fr
!! Line Jourdain, SA 
!! 
!! Anne Cozic, LSCE, anne.cozic@cea.fr
!! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr
!!
!! This software is a computer program whose purpose is to simulate the 
!! atmospheric gas phase and aerosol composition. The model is designed to be
!! used within a transport model or a general circulation model. This version 
!! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts
!! for emissions, transport (resolved and sub-grid scale), photochemical 
!! transformations, and scavenging (dry deposition and washout) of chemical 
!! species and aerosols interactively in the GCM. Several versions of the INCA
!! model are currently used depending on the envisaged applications with the 
!! chemistry-climate model. 
!! 
!! This software is governed by the CeCILL  license under French law and
!! abiding by the rules of distribution of free software.  You can  use, 
!! modify and/ or redistribute the software under the terms of the CeCILL
!! license as circulated by CEA, CNRS and INRIA at the following URL
!! "http://www.cecill.info". 
!! 
!! As a counterpart to the access to the source code and  rights to copy,
!! modify and redistribute granted by the license, users are provided only
!! with a limited warranty  and the software's author,  the holder of the
!! economic rights,  and the successive licensors  have only  limited
!! liability. 
!! 
!! In this respect, the user's attention is drawn to the risks associated
!! with loading,  using,  modifying and/or developing or reproducing the
!! software by the user in light of its specific status of free software,
!! that may mean  that it is complicated to manipulate,  and  that  also
!! therefore means  that it is reserved for developers  and  experienced
!! professionals having in-depth computer knowledge. Users are therefore
!! encouraged to load and test the software's suitability as regards their
!! requirements in conditions enabling the security of their systems and/or 
!! data to be ensured and,  more generally, to use and operate it in the 
!! same conditions as regards security. 
!! 
!! The fact that you are presently reading this means that you have had
!! knowledge of the CeCILL license and that you accept its terms.
!! =========================================================================

SUBROUTINE MKNOPROD(oro,  & 
   lat,  &
   lon,  &
   area, &
   pmid, &
   zmid, &
   temp, &
   ctop, &
   cbot, &
   nx,   &
   ny)
  !----------------------------------------------------------------------
  !	... NOx from lightning
  ! Line Jourdain, SA, 2001. 
  ! Modified, Didier Hauglustaine, IPSL,  08-2001.
  !--------------------------------------------------------
  USE LIGHTNING
  USE CHEM_CONS
  USE CHEM_MODS
  USE SPECIES_NAMES
  USE INCA_DIM

  IMPLICIT NONE
  
  !--------------------------------------------------------
  !       ... Dummy arguments
  !--------------------------------------------------------
  INTEGER, INTENT(in) :: ctop(PLON) !convective clouds top        
  INTEGER, INTENT(in) :: cbot(PLON) !convective clouds bot
  INTEGER, INTENT(in) :: nx, ny
  REAL, INTENT(in)    :: oro(PLON)         !land-sea mask
  REAL, INTENT(in)    :: area(PLON)        !Grid area (km2)
  REAL, INTENT(in)    :: lat(PLON)         !Not used
  REAL, INTENT(in)    :: pmid(PLON,PLEV)   !Pressure (Pa)
  REAL, INTENT(in)    :: temp(PLON,PLEV)   !Temperature
  REAL, INTENT(in)    :: lon(PLON)         !Not used
  REAL, INTENT(in)    :: zmid(PLON,PLEV)   !Altitude (m)
  
  !--------------------------------------------------------
  !       ... Local variables
  !--------------------------------------------------------
  INTEGER  :: i, k, l, lp
  INTEGER  :: lmax, lmin
  
  REAL, PARAMETER :: secpyr = dayspy * 8.64e4
  REAL, PARAMETER :: ap=0.021
  REAL, PARAMETER :: bp=-0.648
  REAL, PARAMETER :: cp=7.493
  REAL, PARAMETER :: dp=-36.54
  REAL, PARAMETER :: ep=63.09
  REAL, PARAMETER :: zero = 1.e-20
  
  REAL     :: dx, dy
  REAL     :: calibration
  REAL     :: source
  REAL     :: ziso0(PLON)
  REAL     :: fic(PLON), fcg(PLON)
  REAL     :: z(PLON)
  REAL     :: a(PLON), b(PLON)
  REAL     :: globoprod_no(PLON)
  REAL     :: xprod_no(PLON,PLEV)
  REAL     :: no_col(PLON)
  REAL     :: factor(PLON)
  REAL     :: scal(PLON)
  REAL     :: altp(PLON,zdim)
  REAL     :: atc(PLON,zdim-1), atm(PLON,zdim-1), amc(PLON,zdim-1)
  REAL     :: btc(PLON,zdim-1), btm(PLON,zdim-1), bmc(PLON,zdim-1)
  REAL     :: coeff_tcp(PLON,PLEV), coeff_mcp(PLON,PLEV), coeff_tmp(PLON,PLEV)
  REAL     :: somcoeff1(PLON),somcoeff2(PLON),somcoeff3(PLON)
  REAL     :: zinterf(PLON,PLEVP)
  REAL     :: deltaz(PLON,PLEV)
  REAL     :: zmidkm(PLON,PLEV)
  
  zmidkm = zmid*1.e-3

  !     ... Zero all variables for this time-step
  ztop=0.
  flash=0.
  ziso0=0.
  dcold=0.
  pcg=0. 
  fcg=0.
  fic=0.
  globoprod_no=0.
  xprod_no = 0.
  factor=0. 
  scal = 1.0

  !     ... Ensure cloud presence
  DO i=1,PLON
    IF (pmid(i,ctop(i)) < pmid(i,cbot(i))) THEN
        ztop(i)=zmidkm(i,ctop(i))
    ENDIF
  ENDDO

  !     ... Flashes per sec (Price and Rind, 92)
  DO i=1,PLON
    IF (ztop(i) > 4.) THEN
        flash(i)=3.44e-5*ztop(i)**4.90 *oro(i) + 6.40e-4*ztop(i)**1.73 *(1.-oro(i))
    ENDIF
  ENDDO
  
  !     ... Two calibrations applied: 
  !     1. To account for cloud height averaging (see Price and Rind, 1994)
  !     2. To ensure 5 TgN in global and annual mean.
  
  dx = 2.*pi/FLOAT(nx)*r2d
  dy = pi/FLOAT(ny)*r2d
  
  calibration = 9.7241e-1*EXP(4.8203e-2*dx*dy)
  flash(:) = 2. * calibration * flash(:)   !KE

  !     ... Find 0 celcius isotherm
  DO i=1,PLON
    DO l=1,PLEV-1
      IF ( pmid(i,l) > 2.e4 ) THEN
          IF ( (temp(i,l) > 273.15) .AND. (temp(i,l+1) < 273.15) ) THEN
              a(i)=(temp(i,l+1)-temp(i,l))/(zmidkm(i,l+1)-zmidkm(i,l))
              b(i)=temp(i,l)-a(i)*zmidkm(i,l)
              ziso0(i)=(273.15-b(i))/a(i)
          ENDIF
      ENDIF
    ENDDO
  ENDDO
  
  !     ... Intracloud and cloud-to-ground lightning

  DO i=1,PLON

    IF (ziso0(i) < ztop(i)) THEN

        dcold(i)=ztop(i)-ziso0(i)
        
        IF ( (dcold(i) > 5.5) .AND. (dcold(i) < 14.) ) THEN
            z(i)=(ap*dcold(i)**4+bp*dcold(i)**3+cp*dcold(i)**2+dp*dcold(i)+ep)
            pcg(i)=1./(z(i)+1.)
        ELSE IF (dcold(i) <= 5.5) THEN
            pcg(i)=0.
        ELSE IF (dcold(i) >= 14.) THEN
            pcg(i)=2.e-2
        ENDIF
        
        fcg(i) = pcg(i)*flash(i)
        fic(i) = (1.-pcg(i))*flash(i)

    ENDIF
    
  ENDDO

!     ... Change units to fl/s/m2 DH.

      flash(:) = flash(:)/60./area(:)
      flpcg(:) = fcg(:)/60./area(:)
  
  !     ... Column integrated NO production
  globoprod_no(:) = 6.7e25 * (10.*fcg(:)+fic(:)) /60./area(:)
  
  !     ... Pickering et al. (1998) vertical distribution
  DO i=1,PLON
    xprod_no(i,1:ctop(i))=1.
  ENDDO
  
  zinterf(:,1)=0.
  zinterf(:,2)=2.*zmidkm(:,1)
  DO l=2,PLEV
    zinterf(:,l+1)=zinterf(:,l)+2.*(zmidkm(:,l)-zinterf(:,l))
  ENDDO
  
  DO l=1,PLEV
    deltaz(:,l)=(zinterf(:,l+1)-zinterf(:,l))*1.e3
  ENDDO
  
      
  DO i=1,PLON
    IF (flash(i) > zero) THEN
        
        IF ( ABS(lat(i)) >= 30.) THEN
            scal(i)=ztop(i)/14.5
        ELSE
            scal(i)=ztop(i)/16.
        ENDIF
        
    ENDIF
  ENDDO
  
  DO l=1,zdim
    altp(:,l)=alt(l)*scal(:)
  ENDDO
  
  DO i=1,PLON
    DO l=1,zdim-1

      IF (scal(i) > zero) THEN
          atc(i,l)=(coeff_tc(l+1)-coeff_tc(l))/(altp(i,l+1)-altp(i,l))
          btc(i,l)=coeff_tc(l)-atc(i,l)*altp(i,l)
          atm(i,l)=(coeff_tm(l+1)-coeff_tm(l))/(altp(i,l+1)-altp(i,l))
          btm(i,l)=coeff_tm(l)-atm(i,l)*altp(i,l)
          amc(i,l)=(coeff_mc(l+1)-coeff_mc(l))/(altp(i,l+1)-altp(i,l))
          bmc(i,l)=coeff_mc(l)-amc(i,l)*altp(i,l)
      ENDIF
      
    ENDDO
  ENDDO
  
  coeff_mcp=0.
  coeff_tmp=0.
  coeff_tcp=0.
  
  DO i=1,PLON
    DO lp=1,zdim-1
      DO l=1,ctop(i)
        
        IF (flash(i) > zero) THEN
            
            IF (zmidkm(i,l) <= altp(i,1)) THEN
                coeff_tcp(i,l)=coeff_tc(1)
                coeff_tmp(i,l)=coeff_tm(1)
                coeff_mcp(i,l)=coeff_mc(1)
            ENDIF
            
            IF ((zmidkm(i,l) <= altp(i,lp+1)) .AND. (zmidkm(i,l) >= altp(i,lp))) THEN
                coeff_tcp(i,l)=atc(i,lp)*zmidkm(i,l)+btc(i,lp)
                coeff_tmp(i,l)=atm(i,lp)*zmidkm(i,l)+btm(i,lp)
                coeff_mcp(i,l)=amc(i,lp)*zmidkm(i,l)+bmc(i,lp)
            ENDIF
        ENDIF
        
      ENDDO
    ENDDO
  ENDDO
  
  somcoeff1=0.
  somcoeff2=0.
  somcoeff3=0.
  
  DO i=1,PLON
    DO l=1,ctop(i)
      somcoeff1(i)=somcoeff1(i)+coeff_tcp(i,l)*deltaz(i,l)
      somcoeff2(i)=somcoeff2(i)+coeff_tmp(i,l)*deltaz(i,l)
      somcoeff3(i)=somcoeff3(i)+coeff_mcp(i,l)*deltaz(i,l)
    ENDDO
  ENDDO
  
  DO i=1,PLON
    DO l=1,ctop(i)
      
      IF (somcoeff1(i) > zero) THEN
          coeff_tcp(i,l)=coeff_tcp(i,l)/somcoeff1(i)*1.e2
      ENDIF
      
      IF (somcoeff2(i) > zero) THEN
          coeff_tmp(i,l)=coeff_tmp(i,l)/somcoeff2(i)*1.e2
      ENDIF
      
      IF (somcoeff3(i) > zero) THEN
          coeff_mcp(i,l)=coeff_mcp(i,l)/somcoeff3(i)*1.e2
      ENDIF
      
    ENDDO
  ENDDO
  
  DO i=1,PLON
    DO l=1,ctop(i)
      IF (flash(i) > zero) THEN
          IF ( ABS(lat(i)) <= 30. ) THEN
              
              xprod_no(i,l)=coeff_tcp(i,l)*0.01 * oro(i) + coeff_tmp(i,l)*0.01 * (1.-oro(i))
          ELSE 
              
              xprod_no(i,l)=coeff_mcp(i,l)*0.01
          ENDIF
      ENDIF
    ENDDO
  ENDDO
  
  no_col=0.
  DO i=1,PLON
    lmax=ctop(i)
    DO l=1,lmax
      no_col(i) = no_col(i) + xprod_no(i,l) * deltaz(i,l) * 1.e6
    ENDDO
  ENDDO
  
  prod_light=0.
  factor=0.
  DO i=1,PLON
    lmax=ctop(i)
    IF (no_col(i) /= 0.) THEN
        factor(i) = globoprod_no(i) / no_col(i)
    ENDIF
    DO l=1,lmax
      prod_light(i,l) = xprod_no(i,l) * factor(i)
    ENDDO
  ENDDO
  
  source = 0.
  DO i=1,PLON
    DO l=1,PLEV
      source = source + prod_light(i,l)*area(i)*deltaz(i,l)*1.e6
    ENDDO
  ENDDO
  source = source /6.02e23 *14. *secpyr *1.e-12

!Integrated lightning NOx production (should be in kg/m2/s since prod_light is in molec/cm3/s)
  DO i=1,PLON
    DO l=1,PLEV
      prod_light_col(i) = prod_light_col(i) + prod_light(i,l)*1.e6*deltaz(i,l)/6.02e23*14.*1.e-3
    ENDDO
  ENDDO
  
END SUBROUTINE MKNOPROD