!
      module noz_module
!
! Advance noz by one time step, and partition into no and no2.
!
! Boundary conditions, production and loss for NOZ are defined
! by comp_noz, and referenced by minor_noz. Comp_noz is called
! from a latitude loop in dynamics. After comp_noz, dynamics calls
! minor_noz, which passes this module data to sub minor. Sub
! minor contains 3d mpi calls and its own latitude loops. 
!
! Sub minor_noz also partitions NOZ into NO2 and NO.
!
      use params_module,only: nlevp1
      use addfld_module,only: addfld
      implicit none
!
! Boundary conditions and production and loss terms are allocated
! subdomains by sub alloc_noz (called from allocdata.F).
!
      real,allocatable,dimension(:,:),save   :: noz_ubc ! upper boundary (i,j)
      real,allocatable,dimension(:,:,:),save :: noz_lbc ! lower boundary (i,3,j)
      real,allocatable,dimension(:,:,:),save ::
     |  noz_prod,     ! production of noz (k,i,j)
     |  noz_loss      ! loss of noz       (k,i,j)
!
! Ratios:
      real,dimension(:,:,:),allocatable,save :: ! (nlevp1,lon0:lon1,lat0:lat1)
     |  ratio4,ratio5,rmtru_noz
!
      contains
!-----------------------------------------------------------------------
      subroutine comp_noz(tn,o2,o1,o3,no,no2,n2o,n4s,oh,hox,ne,n2d,o1d,
     |  barm,xnmbar,o2p,lev0,lev1,lon0,lon1,lat)
!
! Advance noz.
!
      use cons_module,only: rmassinv_o3,rmassinv_o2,rmassinv_o1,
     |  rmassinv_n4s,rmassinv_n2d,rmassinv_hox,rmassinv_no,rmass_no,
     |  rmass_no2
      use qrj_module,only: 
     |  pdno2,     ! xjno2
     |  pdnosrb,   ! xjno
     |  pdnoeuv,   ! dnoeuv: photodissociation of no by euv
     |  qnoplya,   ! xjnopn: no+ nighttime ionization by lyman-alpha
     |  qnolya     ! xjnop:  no  ionization by lyman-alpha
      use comp_meta_module,only: ! (nlevp1,lon0:lon1,lat0:lat1)
     |  oh_h,      ! oh/h  (ratio1)
     |  ho2_h,     ! ho2/h (ratio2)
     |  h_hox,     ! h/hox (ratio3)
     |  rmtru_hox  ! "true" mass of hox
      use chemrates_module,only: beta1,beta2,beta3,beta6,beta8,beta9,
     |  beta10,beta11,beta12,beta13,beta14,rk5
      use solgar_module,only: xnozlb ! (nlat)
      implicit none
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn       ,   ! tn (deg K) 
     |  o2       ,   ! o2 (mmr)
     |  o1       ,   ! o1 (mmr)
     |  o3       ,   ! o3 (mmr)
     |  no       ,   ! no (mmr)
     |  no2      ,   ! no2(mmr)
     |  n2o      ,   ! n2o(mmr) from comp_solgar
     |  n4s      ,   ! n4s(mmr)
     |  oh       ,   ! oh (mmr)
     |  hox      ,   ! hox(mmr)
     |  ne       ,   ! ne (cm3)
     |  n2d      ,   ! updated n2d from comp_meta
     |  o1d      ,   ! o(1d) from comp_meta
     |  barm     ,   ! mbar
     |  xnmbar   ,   ! n*mbar
     |  o2p          ! o2+
!
! Local:
      integer :: k,i,i0,i1,nk,nkm1
      real :: test
!
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
      do i=lon0,lon1
        do k=lev0,lev1-1
          ratio4(k,i,lat) =                                        ! s1
     |      (beta9(k,i,lat)*o3(k,i)*rmassinv_o3+beta10(k,i,lat)*
     |      hox(k,i)/rmtru_hox(k,i,lat)*ho2_h(k,i,lat)*h_hox(k,i,lat))/
     |      (beta11*o1(k,i)*rmassinv_o1+beta12(k,i,lat)*o3(k,i)*
     |      rmassinv_o3+0.5*(pdno2(k,i,lat)+pdno2(k+1,i,lat))/
     |      xnmbar(k,i)+beta13*n4s(k,i)*rmassinv_n4s)

          ratio5(k,i,lat) = 1./(1.+ratio4(k,i,lat))                ! s2
          if (ratio5(k,i,lat) < 1.e-6) ratio5(k,i,lat) = 1.e-6
          test = (1.-ratio5(k,i,lat))/ratio5(k,i,lat)
          if (ratio4(k,i,lat) < test) ratio4(k,i,lat) = test

          rmtru_noz(k,i,lat) = ratio5(k,i,lat)*              ! rmtru1
     |      (rmass_no+ratio4(k,i,lat)*rmass_no2)

          if (ratio4(k,i,lat) < 1.e-6) ratio4(k,i,lat) = 1.e-6
          if (ratio5(k,i,lat) < 1.e-6) ratio5(k,i,lat) = 1.e-6
          if (rmtru_noz(k,i,lat) < 1.e-6) rmtru_noz(k,i,lat) = 1.e-6

        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!
! Some fields are saved only up to level lev1-1:
!     call addfld('RATIO4'  ,' ',' ',ratio4(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('RATIO5'  ,' ',' ',ratio5(lev0:lev1-1,:,lat),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('XNMBAR'  ,' ',' ',xnmbar(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('BETA1'   ,' ',' ',beta1(lev0:lev1-1,lon0:lon1,lat),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('RMTRU'   ,' ',' ',rmtru_hox(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('OH_H'    ,' ',' ',oh_h(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('H_HOX'   ,' ',' ',h_hox(lev0:lev1-1,:,lat),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('RMTRUNOZ',' ',' ',rmtru_noz(lev0:lev1-1,:,lat),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('MBAR',' ',' ',barm(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      do i=lon0,lon1
        do k=lev0,lev1-1
          noz_prod(k,i,lat) = xnmbar(k,i)**2*(o2(k,i)*rmassinv_o2*  ! s2
     |      (beta1(k,i,lat)*n4s(k,i)*rmassinv_n4s+beta2*n2d(k,i)*
     |      rmassinv_n2d)+beta8*hox(k,i)/rmtru_hox(k,i,lat)*
     |      oh_h(k,i,lat)*h_hox(k,i,lat)*n4s(k,i)*rmassinv_n4s)+2.*
     |      beta14*o1d(k,i)*rmassinv_o1*xnmbar(k,i)*n2o(k,i)*
     |      xnmbar(k,i)/(0.5*(barm(k,i)+barm(k+1,i)))
          noz_prod(k,i,lat) = noz_prod(k,i,lat)*rmtru_noz(k,i,lat)*
     |      rmassinv_no
! 
          noz_loss(k,i,lat) = -(xnmbar(k,i)*(beta3(k,i,lat)*n4s(k,i)* ! s1
     |      rmassinv_n4s+beta6*n2d(k,i)*rmassinv_n2d)+
     |      0.5*(pdnosrb(k,i,lat)+pdnosrb(k+1,i,lat)+
     |           qnolya (k,i,lat)+qnolya (k+1,i,lat)+
     |           qnoplya(k,i,lat)+qnoplya(k+1,i,lat)+
     |           pdnoeuv(k,i,lat)+pdnoeuv(k+1,i,lat))+
     |      rk5*o2p(k,i))*ratio5(k,i,lat)-xnmbar(k,i)*(beta13*n4s(k,i)*
     |      rmassinv_n4s+beta12(k,i,lat)*o3(k,i)*rmassinv_o3)*
     |      ratio4(k,i,lat)*ratio5(k,i,lat)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('NOZ_PROD',' ',' ',noz_prod(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('NOZ_LOSS',' ',' ',noz_loss(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      do i=lon0,lon1 
        noz_lbc(i,1,lat) = 0.    ! t1
        noz_lbc(i,2,lat) = 1.    ! t2
!
! Number density mixing ratio of noz at lower boundary
! (xnozlb is from solgar module solgar.F)
!
        noz_lbc(i,3,lat) = -xnozlb(lat)*rmtru_noz(lev0,i,lat)/
     |    barm(lev0,i)  ! t3
!
! Zero diffusive flux at top:
        noz_ubc(i,lat) = 0.      ! t4
      enddo ! i=lon0,lon1 
!
      end subroutine comp_noz
!-----------------------------------------------------------------------
      subroutine minor_noz(tn,o2,ox,w,difkk,noz,noz_nm,hdnoz,noz_out,
     |  noznm_out,no_out,no2_out,lev0,lev1,lon0,lon1,lat0,lat1)
!
! Call minor to complete noz. This is called from dynamics, once
! per time step after comp_noz (which is called from a latitude loop).
! Inputs are at 3d subdomains. 
!
      use compdat_module,only: wnoz ! (nlevp1)
      use cons_module,only: rmass_no,rmass_no2
#ifdef MPI
      use mpi_module,only: mp_periodic_f3d
#endif
      implicit none
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(in) ::
     |  tn,      ! neutral temperature (deg K)
     |  o2,      ! molecular oxygen (mmr)
     |  ox,      ! oxygen family (mmr)
     |  w,       ! omega (vertical velocity)
     |  difkk,   ! eddy viscosity (from mgw)
     |  noz,     ! noz (mmr)
     |  noz_nm,  ! noz at time n-1
     |  hdnoz    ! horizontal diffusion
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  noz_out,    ! noz output
     |  noznm_out,  ! noz output at time n-1
     |  no_out,     ! no output
     |  no2_out     ! no2 output
!
! Local:
      integer :: i,k,ibndbot,ibndtop,nk,nkm1,i0,i1,lat
      real,parameter :: xynoz = 1.e-3
      real,parameter :: phi_noz(3) = (/0.814, 0.866, 0.926/)
      real,parameter :: alfa_noz = 0. ! thermal diffusion coefficient
      real :: test
      real :: ftmp(lev0:lev1,lon0:lon1,lat0:lat1,2)
      real,parameter :: noz_min = 1.e-10
      integer :: nzero

      ibndbot = 0 ! ilbc in minor
      ibndtop = 0 ! iubc in minor
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1

      call minor(tn,o2,ox,w,difkk,noz,noz_nm,hdnoz,noz_lbc,
     |  noz_ubc,noz_loss,noz_prod,wnoz,rmass_no,phi_noz,alfa_noz,
     |  xynoz,ibndbot,ibndtop,noz_out,noznm_out,
     |  lev0,lev1,lon0,lon1,lat0,lat1,0)

#ifdef MPI
      ftmp(:,:,:,1) = noz_out(:,lon0:lon1,lat0:lat1)
      ftmp(:,:,:,2) = noznm_out(:,lon0:lon1,lat0:lat1)

      call mp_periodic_f3d(ftmp,lev0,lev1,lon0,lon1,lat0,lat1,2)

      noz_out(:,lon0:lon1,lat0:lat1)   = ftmp(:,:,:,1)
      noznm_out(:,lon0:lon1,lat0:lat1) = ftmp(:,:,:,2)
#else
!
! Set periodic points for non-MPI serial run:
      do k=lev0,lev1
        noz_out(k,1:2,:) = noz_out(k,lon1-3:lon1-2,:)
        noz_out(k,lon1-1:lon1,:) = noz_out(k,3:4,:)
        noznm_out(k,1:2,:) = noznm_out(k,lon1-3:lon1-2,:)
        noznm_out(k,lon1-1:lon1,:) = noznm_out(k,3:4,:)
      enddo
#endif
!
! Insure positive mmr:
! 8/6/09: noz_out can go to zero. This does not cause the model
!         to crash, but in some cases can cause tgcmproc_f90 to
!         hang in cplbam (sub contour in plt.F) when contouring 
!         log10 of NOZ, but sometimes it works, properly inserting 
!         missing values and plotting white boxes in the plots. 
!         Saw this at equatorial latitudes around zp -12 in slt 
!         noon longitude slice.
!
      nzero = 0
      do lat=lat0,lat1
        do i=lon0,lon1
          do k=lev0,lev1
            if (noz_out(k,i,lat) <= 0.) then
!             write(6,"('comp_noz found value <= 0: k,i,lat=',3i5,
!    |          ' noz_out(k,i,lat)=',e12.4,' setting it to noz_min=',
!    |          e12.4)") k,i,lat,noz_out(k,i,lat),noz_min
              noz_out(k,i,lat) = noz_min
              nzero = nzero+1
            endif
          enddo
        enddo
      enddo
!
! 8/6/09: This is typically 1-15 values:
!     if (nzero > 0) write(6,"('comp_noz: nzero in noz_out=',i6)")
!    |  nzero

!     do lat=lat0,lat1
!       call addfld('NOZ_MNR' ,' ',' ',noz_out  (:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('NOZM_MNR',' ',' ',noznm_out(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!     enddo
!
! Partition NOZ into NO2 and NO:
!
! Print for debug:
!     do lat=lat0,lat1
!       do i=lon0,lon1
!         do k=lev0,lev1-1
!           test = xynoz*wnoz(k)
!           if (noz_out(k,i,lat) < test) noz_out(k,i,lat) = test
!         enddo
!         write(6,"('comp_noz: lat=',i4,' i=',i4,
!    |      ' noz_out(lev0:lev1-1,i,lat)=',/,(6e12.4))")
!    |      lat,i,noz_out(lev0:lev1-1,i,lat)   
!       enddo
!     enddo

      do lat=lat0,lat1
        do i=lon0,lon1
          do k=lev0,lev1-1 ! changed 12/17/03: lev1 to lev1-1
            test = xynoz*wnoz(k)
            if (noz_out(k,i,lat) < test) noz_out(k,i,lat) = test
            no_out(k,i,lat) = ratio5(k,i,lat)*noz_out(k,i,lat)*
     |        rmass_no/rmtru_noz(k,i,lat)
            no2_out(k,i,lat) = ratio4(k,i,lat)*ratio5(k,i,lat)*
     |        noz_out(k,i,lat)*rmass_no2/rmtru_noz(k,i,lat)
          enddo ! k=lev0,lev1
!
! Set top level using log interpolation:
!
!         write(6,"('minor_noz: i=',i4,' lat=',i4,' no_out(:,i,lat)=',
!    |      /(6e12.4))") i,lat,
!    |      no_out(:,i,lat)

          no_out(lev1,i,lat) = no_out(lev1-1,i,lat)**2/
     |      no_out(lev1-2,i,lat)
          no2_out(lev1,i,lat) = no2_out(lev1-1,i,lat)**2/
     |      no2_out(lev1-2,i,lat)
        enddo ! lon=lon0,lon1

!       call addfld('NOZ_OUT' ,' ',' ',noz_out  (:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('NOZNMOUT',' ',' ',noznm_out(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('NO_OUT'  ,' ',' ',no_out   (:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('NO2_OUT' ,' ',' ',no2_out  (:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

      enddo ! lat=lat0,lat1
!
      end subroutine minor_noz
!-----------------------------------------------------------------------
      subroutine alloc_noz(lon0,lon1,lat0,lat1)
!
! Allocate subdomains (without ghost cells) to module data for boundary
! conditions and production and loss terms. This is called once per run
! from sub allocdata (allocdata.F).
!
! Args:
      integer,intent(in) :: lon0,lon1,lat0,lat1
!
! Local:
      integer :: istat
!
! Allocate subdomains to boundary conditions:
      allocate(noz_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_noz: error allocating',
     |  ' noz_ubc: stat=',i3)") istat
      allocate(noz_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_noz: error allocating',
     |  ' noz_lbc: stat=',i3)") istat
      noz_ubc = 0. ; noz_lbc = 0.
!
! Allocate subdomains to production and loss:
      allocate(noz_prod(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_noz: error allocating',
     |  ' noz_prod: stat=',i3)") istat
      allocate(noz_loss(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_noz: error allocating',
     |  ' noz_loss: stat=',i3)") istat
      noz_prod = 0. ; noz_loss = 0.
!
! Allocate rmtru_noz and ratios:
      allocate(ratio4(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_noz: error allocating',
     |  ' ratio4: stat=',i3)") istat
      allocate(ratio5(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_noz: error allocating',
     |  ' ratio5: stat=',i3)") istat
      allocate(rmtru_noz(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_noz: error allocating',
     |  ' rmtru_noz: stat=',i3)") istat
      ratio4 = 0. ; ratio5 = 0. ; rmtru_noz = 0.
!
!     write(6,"('alloc_noz: allocated module data')")
!
      end subroutine alloc_noz
!-----------------------------------------------------------------------
      end module noz_module
