module mo_imp_sol

! update O3 production and loss calculations for new HCs (lke/pgmh 5/2/2007)

      use chem_mods, only : clscnt4

      implicit none

      private
      public :: imp_slv_inti
      public :: imp_slv_lt_inti
      public :: imp_sol

      save

      integer, parameter ::  inst = 1, avrg = 2
      real, parameter    ::  rel_err      = 1.e-3
      real, parameter    ::  high_rel_err = 1.e-4
      real, parameter    ::  small_val    = 1.e-40
!-----------------------------------------------------------------------      
!   	    newton-raphson iteration limits
!-----------------------------------------------------------------------      
      integer, parameter :: cut_limit = 5
      integer, parameter :: vec_len   = 64

      integer :: ox_ndx
      integer :: o1d_ndx = -1
      integer :: oh_ndx, ho2_ndx, ch3o2_ndx, po2_ndx, ch3co3_ndx
      integer :: c2h5o2_ndx, isopo2_ndx, macro2_ndx, mco3_ndx, c3h7o2_ndx
      integer :: ro2_ndx, xo2_ndx, no_ndx, no2_ndx, no3_ndx, n2o5_ndx
      integer :: c2h4_ndx, c3h6_ndx, isop_ndx, mvk_ndx, c10h16_ndx
      integer :: ox_p1_ndx, ox_p2_ndx, ox_p3_ndx, ox_p4_ndx, ox_p5_ndx
      integer :: ox_p6_ndx, ox_p7_ndx, ox_p8_ndx, ox_p9_ndx, ox_p10_ndx
      integer :: ox_p11_ndx
      integer :: ox_l1_ndx, ox_l2_ndx, ox_l3_ndx, ox_l4_ndx, ox_l5_ndx
      integer :: ox_l6_ndx, ox_l7_ndx, ox_l8_ndx, ox_l9_ndx, usr4_ndx
      integer :: usr16_ndx, usr17_ndx
!o3_pl+++
      integer :: tolo2_ndx, terpo2_ndx, alko2_ndx, eneo2_ndx, eo2_ndx, meko2_ndx
      integer :: ox_p12_ndx,ox_p13_ndx,ox_p14_ndx,ox_p15_ndx,ox_p16_ndx, ox_p17_ndx
!o3_pl---
      integer :: lt_cnt
      logical :: do_ox_pl = .true.

      type hst_pl
         integer  ::  cnt(2)
         logical  ::  do_hst(2)
      end type hst_pl

      real, private                      ::   epsilon(clscnt4)
      type(hst_pl), private, allocatable ::   imp_hst_prod(:)
      type(hst_pl), private, allocatable ::   imp_hst_loss(:)
      logical, private, allocatable      ::   factor(:)

      contains

      subroutine imp_slv_inti
!-----------------------------------------------------------------------      
!	... initialize the implict solver
!-----------------------------------------------------------------------      

      use chem_mods,    only : clscnt4, implicit
      use mo_grid,      only : pcnstm1
      use mo_histout,   only : hfile, moz_file_cnt
      use mo_chem_utls, only : get_spc_ndx, get_rxt_ndx

      implicit none

!-----------------------------------------------------------------------      
!	... local variables
!-----------------------------------------------------------------------      
      integer :: m, astat, file, timetype
      integer :: il, iu
!o3_pl
      integer :: wrk(27)
      real    :: eps(pcnstm1)

      allocate( factor(implicit%iter_max),stat=astat )
      if( astat /= 0 ) then
         write(*,*) 'imp_slv_inti: failed to allocate factor array; error = ',astat
         call endrun
      end if

      factor(:) = .true.
      eps(:)    = rel_err
      ox_ndx = get_spc_ndx( 'OX' )
      if( ox_ndx < 1 ) then
         ox_ndx  = get_spc_ndx( 'O3' )
         o1d_ndx = get_spc_ndx( 'O1D' )
      end if
      if( ox_ndx > 0 ) then
         eps(ox_ndx) = high_rel_err
      end if
      m = get_spc_ndx( 'NO' )
      if( m > 0 ) then
         eps(m) = high_rel_err
      end if
      m = get_spc_ndx( 'NO2' )
      if( m > 0 ) then
         eps(m) = high_rel_err
      end if
      m = get_spc_ndx( 'NO3' )
      if( m > 0 ) then
         eps(m) = high_rel_err
      end if
      m = get_spc_ndx( 'HNO3' )
      if( m > 0 ) then
         eps(m) = high_rel_err
      end if
      m = get_spc_ndx( 'HO2NO2' )
      if( m > 0 ) then
         eps(m) = high_rel_err
      end if
      m = get_spc_ndx( 'N2O5' )
      if( m > 0 ) then
         eps(m) = high_rel_err
      end if
      m = get_spc_ndx( 'OH' )
      if( m > 0 ) then
         eps(m) = high_rel_err
      end if
      m = get_spc_ndx( 'HO2' )
      if( m > 0 ) then
         eps(m) = high_rel_err
      end if
      do m = 1,max(1,clscnt4)
         epsilon(m) = eps(implicit%clsmap(m))
      end do

      if( ox_ndx > 0 ) then
         ox_p1_ndx = get_rxt_ndx( 'ox_p1' )
         ox_p2_ndx = get_rxt_ndx( 'ox_p2' )
         ox_p3_ndx = get_rxt_ndx( 'ox_p3' )
         ox_p4_ndx = get_rxt_ndx( 'ox_p4' )
         ox_p5_ndx = get_rxt_ndx( 'ox_p5' )
         ox_p6_ndx = get_rxt_ndx( 'ox_p6' )
         ox_p7_ndx = get_rxt_ndx( 'ox_p7' )
         ox_p8_ndx = get_rxt_ndx( 'ox_p8' )
         ox_p9_ndx = get_rxt_ndx( 'ox_p9' )
         ox_p10_ndx = get_rxt_ndx( 'ox_p10' )
         ox_p11_ndx = get_rxt_ndx( 'ox_p11' )
!o3_pl++
         ox_p12_ndx = get_rxt_ndx( 'ox_p12' )
         ox_p13_ndx = get_rxt_ndx( 'ox_p13' )
         ox_p14_ndx = get_rxt_ndx( 'ox_p14' )
         ox_p15_ndx = get_rxt_ndx( 'ox_p15' )
         ox_p16_ndx = get_rxt_ndx( 'ox_p16' )
         ox_p17_ndx = get_rxt_ndx( 'ox_p17' )

         wrk(1:17) = (/ ox_p1_ndx, ox_p2_ndx, ox_p3_ndx, ox_p4_ndx, ox_p5_ndx, &
                        ox_p6_ndx, ox_p7_ndx, ox_p8_ndx, ox_p9_ndx, ox_p10_ndx, ox_p11_ndx, &
                        ox_p12_ndx, ox_p13_ndx, ox_p14_ndx, ox_p15_ndx, ox_p16_ndx, ox_p17_ndx  /)
         if( any( wrk(1:17) < 1 ) ) then
!o3_pl--
            do_ox_pl = .false.
         end if
         if( do_ox_pl ) then
            ox_l1_ndx = get_rxt_ndx( 'ox_l1' )
            ox_l2_ndx = get_rxt_ndx( 'ox_l2' )
            ox_l3_ndx = get_rxt_ndx( 'ox_l3' )
            ox_l4_ndx = get_rxt_ndx( 'ox_l4' )
            ox_l5_ndx = get_rxt_ndx( 'ox_l5' )
            ox_l6_ndx = get_rxt_ndx( 'ox_l6' )
            ox_l7_ndx = get_rxt_ndx( 'ox_l7' )
            ox_l8_ndx = get_rxt_ndx( 'ox_l8' )
            ox_l9_ndx = get_rxt_ndx( 'ox_l9' )
            if( ox_l9_ndx < 1 ) then
               ox_l9_ndx = get_rxt_ndx( 'soa1' )
            end if
            usr4_ndx = get_rxt_ndx( 'usr4' )
            usr16_ndx = get_rxt_ndx( 'usr16' )
            usr17_ndx = get_rxt_ndx( 'usr17' )
            wrk(1:12) = (/ ox_l1_ndx, ox_l2_ndx, ox_l3_ndx, ox_l4_ndx, ox_l5_ndx, &
                           ox_l6_ndx, ox_l7_ndx, ox_l8_ndx, ox_l9_ndx, usr4_ndx, &
                           usr16_ndx, usr17_ndx /)
            if( any( wrk(1:12) < 1 ) ) then
               do_ox_pl = .false.
            end if
         end if

         if( do_ox_pl ) then
            oh_ndx = get_spc_ndx( 'OH' )
            ho2_ndx = get_spc_ndx( 'HO2' )
            ch3o2_ndx = get_spc_ndx( 'CH3O2' )
            po2_ndx = get_spc_ndx( 'PO2' )
            ch3co3_ndx = get_spc_ndx( 'CH3CO3' )
            c2h5o2_ndx = get_spc_ndx( 'C2H5O2' )
            macro2_ndx = get_spc_ndx( 'MACRO2' )
            mco3_ndx = get_spc_ndx( 'MCO3' )
            c3h7o2_ndx = get_spc_ndx( 'C3H7O2' )
            ro2_ndx = get_spc_ndx( 'RO2' )
            xo2_ndx = get_spc_ndx( 'XO2' )
            no_ndx = get_spc_ndx( 'NO' )
            no2_ndx = get_spc_ndx( 'NO2' )
            no3_ndx = get_spc_ndx( 'NO3' )
            n2o5_ndx = get_spc_ndx( 'N2O5' )
            c2h4_ndx = get_spc_ndx( 'C2H4' )
            c3h6_ndx = get_spc_ndx( 'C3H6' )
            isop_ndx = get_spc_ndx( 'ISOP' )
            isopo2_ndx = get_spc_ndx( 'ISOPO2' )
            mvk_ndx = get_spc_ndx( 'MVK' )
            c10h16_ndx = get_spc_ndx( 'C10H16' )
!o3_pl+++
            tolo2_ndx = get_spc_ndx('TOLO2') 
            terpo2_ndx =get_spc_ndx('TERPO2')
            alko2_ndx = get_spc_ndx('ALKO2')
            eneo2_ndx = get_spc_ndx('ENEO2')
            eo2_ndx =   get_spc_ndx('EO2')
            meko2_ndx = get_spc_ndx('MEKO2')
            wrk(1:27) = (/ oh_ndx, ho2_ndx, ch3o2_ndx, po2_ndx, ch3co3_ndx, &
                           c2h5o2_ndx, macro2_ndx, mco3_ndx, c3h7o2_ndx, ro2_ndx, &
                           xo2_ndx, no_ndx, no2_ndx, no3_ndx, n2o5_ndx, &
                           c2h4_ndx, c3h6_ndx, isop_ndx, isopo2_ndx, mvk_ndx, c10h16_ndx, &
                           tolo2_ndx, terpo2_ndx, alko2_ndx, eneo2_ndx, eo2_ndx, meko2_ndx  /)
            if( any( wrk(1:27) < 1 ) ) then
!o3_pl---
               do_ox_pl = .false.
            end if
         end if
      else
         do_ox_pl = .false.
      end if

      if( moz_file_cnt > 0 ) then
         allocate( imp_hst_prod(moz_file_cnt),stat=astat )
         if( astat /= 0 ) then
            write(*,*) 'imp_slv_inti: failed to allocate imp_hst_prod array; error = ',astat
            call endrun
         end if
         if( astat /= 0 ) then
            write(*,*) 'imp_slv_inti: failed to allocate imp_hst_prod array; error = ',astat
            call endrun
         end if
         do file = 1,moz_file_cnt
            imp_hst_prod(file)%do_hst(:) = .false.
            imp_hst_prod(file)%cnt(:)    = 0
         end do
         allocate( imp_hst_loss(moz_file_cnt),stat=astat )
         if( astat /= 0 ) then
            write(*,*) 'imp_slv_inti: failed to allocate imp_hst_loss array; error = ',astat
            call endrun
         end if
         do file = 1,moz_file_cnt
            imp_hst_loss(file)%do_hst(:) = .false.
            imp_hst_loss(file)%cnt(:)    = 0
         end do
!-----------------------------------------------------------------------      
!	... scan for class production to history file(s)
!-----------------------------------------------------------------------      
         do file = 1,moz_file_cnt
            do timetype = inst,avrg
               if( hfile(file)%histout_cnt(14,timetype) > 0 ) then
                  il = hfile(file)%histout_ind(14,timetype)
                  iu = il + hfile(file)%histout_cnt(14,timetype) - 1
                  if( timetype == inst ) then
                     if( any( hfile(file)%inst_map(il:iu)/1000 == 4 ) ) then
                        imp_hst_prod(file)%do_hst(timetype) = .true.
                        do m = il,iu
                           if( hfile(file)%inst_map(m)/1000 == 4 ) then
                              imp_hst_prod(file)%cnt(timetype) = imp_hst_prod(file)%cnt(timetype) + 1
                           end if
                        end do
                        cycle
                     end if
                  else if( timetype == avrg ) then
                     if( any( hfile(file)%timav_map(il:iu)/1000 == 4 ) ) then
                        imp_hst_prod(file)%do_hst(timetype) = .true.
                        do m = il,iu
                           if( hfile(file)%timav_map(m)/1000 == 4 ) then
                              imp_hst_prod(file)%cnt(timetype) = imp_hst_prod(file)%cnt(timetype) + 1
                           end if
                        end do
                        exit
                     end if
                  end if
               end if
            end do
         end do
!-----------------------------------------------------------------------      
!	... scan for class loss to history file(s)
!-----------------------------------------------------------------------      
         do file = 1,moz_file_cnt
            do timetype = inst,avrg
               if( hfile(file)%histout_cnt(15,timetype) > 0 ) then
                  il = hfile(file)%histout_ind(15,timetype)
                  iu = il + hfile(file)%histout_cnt(15,timetype) - 1
                  if( timetype == inst ) then
                     if( any( hfile(file)%inst_map(il:iu)/1000 == 4 ) ) then
                        imp_hst_loss(file)%do_hst(timetype) = .true.
                        do m = il,iu
                           if( hfile(file)%inst_map(m)/1000 == 4 ) then
                              imp_hst_loss(file)%cnt(timetype) = imp_hst_loss(file)%cnt(timetype) + 1
                           end if
                        end do
                        cycle
                     end if
                  else if( timetype == avrg ) then
                     if( any( hfile(file)%timav_map(il:iu)/1000 == 4 ) ) then
                        imp_hst_loss(file)%do_hst(timetype) = .true.
                        do m = il,iu
                           if( hfile(file)%timav_map(m)/1000 == 4 ) then
                              imp_hst_loss(file)%cnt(timetype) = imp_hst_loss(file)%cnt(timetype) + 1
                           end if
                        end do
                        exit
                     end if
                  end if
               end if
            end do
         end do
      end if

      end subroutine imp_slv_inti

      subroutine imp_slv_lt_inti
!-----------------------------------------------------------------------      
!	... initialize the lifetimes for implict solver
!-----------------------------------------------------------------------      

      use mo_lifetime, only : lt_struct
      use mo_histout,  only : moz_file_cnt

!-----------------------------------------------------------------------      
!	... local variables
!-----------------------------------------------------------------------      
      integer :: file

      lt_cnt = 0
      do file = 1,moz_file_cnt
         lt_cnt = max( lt_cnt,sum( lt_struct(4,file)%cnt(:) ) )
      end do

      end subroutine imp_slv_lt_inti

      subroutine imp_sol( base_sol, reaction_rates, het_rates, extfrc, nstep, &
                          delt, hnm, pdel, lat, ip, &
                          plonl, plnplv )
!-----------------------------------------------------------------------
!      	... imp_sol advances the volumetric mixing ratio
!           forward one time step via the fully implicit euler scheme.
!           this source is meant for vector cpus such as those from intel
!-----------------------------------------------------------------------

      use chem_mods,         only : clscnt4, imp_nzcnt, clsze, rxntot, hetcnt, extcnt, implicit
      use mo_histout,        only : outfld, hfile, moz_file_cnt
      use m_tracname,        only : tracnam
      use mo_control,        only : pdiags
      use mo_grid,           only : plev, pcnstm1
      use mo_indprd,         only : indprd
      use mo_imp_prod_loss,  only : imp_prod_loss
      use mo_imp_lin_matrix, only : imp_linmat
      use mo_imp_nln_matrix, only : imp_nlnmat
      use mo_imp_factor,     only : imp_lu_fac
      use mo_imp_solve,      only : imp_lu_slv
      use mo_lifetime,       only : lt_struct

      implicit none

!-----------------------------------------------------------------------
!     	... dummy args
!-----------------------------------------------------------------------
      integer, intent(in) ::   nstep                     ! time step index (zero based)
      integer, intent(in) ::   lat                       ! lat index
      integer, intent(in) ::   ip                        ! longitude tile index
      integer, intent(in) ::   plonl                     ! longitude tile dimension
      integer, intent(in) ::   plnplv                    ! plonl*plev
      real, intent(in)    ::   delt                      ! time step (seconds)
      real, intent(in)    ::   reaction_rates(plnplv,rxntot)
      real, intent(in)    ::   het_rates(plnplv,hetcnt)
      real, intent(in)    ::   extfrc(plnplv,extcnt)
      real, intent(in)    ::   hnm(plnplv)
      real, intent(in)    ::   pdel(plnplv)
      real, intent(inout) ::   base_sol(plnplv,pcnstm1)

!-----------------------------------------------------------------------
!     	... local variables
!-----------------------------------------------------------------------
      type hst_buff
         real, pointer, dimension(:,:) :: buff
      end type hst_buff

      integer ::   nr_iter, &
                   ofl, ofu, &
                   lev, &
                   i, &
                   j, &
                   k, l, &
                   m, &
                   n, &
                   class, file, &
                   cut_cnt, stp_con_cnt
      integer ::   base_ndx
      integer ::   cls_ndx
      integer ::   cls_ndxp
      integer ::   timetype, hndx
      integer ::   astat
      real :: interval_done, dt, dti
      real :: max_delta(max(1,clscnt4))
      real :: sys_jac(plnplv,max(1,imp_nzcnt))
      real :: lin_jac(plnplv,max(1,imp_nzcnt))
      real :: solution(plnplv,max(1,clscnt4))
      real :: forcing(plnplv,max(1,clscnt4))
      real :: iter_invariant(plnplv,max(1,clscnt4))
      real :: prod(plnplv,max(1,clscnt4))
      real :: loss(plnplv,max(1,clscnt4))
      real :: ind_prd(plnplv,max(1,clscnt4))
      real :: wrk(plnplv)
      real :: wrk_buff(plnplv)
      real, allocatable :: base_sol_sav(:,:)
      logical ::   convergence
      logical ::   dump_buff, fill_buff
      logical ::   iter_conv(plnplv)
      logical ::   converged(max(1,clscnt4))
      logical ::   spc_conv(plnplv,max(1,clscnt4))
      character(len=32) :: fldname
      type(hst_buff), allocatable :: prod_buff(:)
      type(hst_buff), allocatable :: loss_buff(:)

      if( moz_file_cnt > 0 ) then
!-----------------------------------------------------------------------      
!        ... allocate history prod, loss buffers
!-----------------------------------------------------------------------      
         allocate( prod_buff(moz_file_cnt), stat=astat )
         if( astat /= 0 ) then
            write(*,*) 'imp_slv_inti: failed to allocate prod_buff; error = ',astat
            call endrun
         end if
         do file = 1,moz_file_cnt
            n = sum( imp_hst_prod(file)%cnt(:) )
            if( n > 0 ) then
               allocate( prod_buff(file)%buff(plnplv,n), stat=astat )
               if( astat /= 0 ) then
                  write(*,*) 'imp_slv: failed to allocate prod buff for file = ',file,'; error = ',astat
                  call endrun
               end if
            else
               nullify( prod_buff(file)%buff )
            end if
         end do
         allocate( loss_buff(moz_file_cnt), stat=astat )
         if( astat /= 0 ) then
            write(*,*) 'imp_slv: failed to allocate loss_buff; error = ',astat
            call endrun
         end if
         do file = 1,moz_file_cnt
            n = sum( imp_hst_loss(file)%cnt(:) )
            if( n > 0 ) then
               allocate( loss_buff(file)%buff(plnplv,n), stat=astat )
               if( astat /= 0 ) then
                  write(*,*) 'imp_slv: failed to allocate loss buff for file = ',file,'; error = ',astat
                  call endrun
               end if
            else
               nullify( loss_buff(file)%buff )
            end if
         end do
      end if

      if( cut_limit > 0 ) then
         allocate( base_sol_sav(plnplv,pcnstm1),stat=astat )
         if( astat /= 0 ) then
            write(*,*) 'imp_slv: failed to allocate base_sol_sav ; error = ',astat
            call endrun
         end if
      end if

      if( implicit%indprd_cnt > 0 .or. extcnt > 0 ) then
!-----------------------------------------------------------------------      
!        ... class independent forcing
!-----------------------------------------------------------------------      
         call indprd( 4, ind_prd, base_sol, extfrc, reaction_rates )
      else
         do m = 1,max(1,clscnt4)
            ind_prd(:,m) = 0.
         end do
      end if
spatial_loop : &
      do ofl = 1,plnplv,vec_len
        ofu = min( plnplv,ofl + vec_len - 1 )
        if( cut_limit > 0 ) then
           do m = 1,pcnstm1
              base_sol_sav(ofl:ofu,m) = base_sol(ofl:ofu,m)
           end do
        end if
!-----------------------------------------------------------------------      
!        ... time step loop
!-----------------------------------------------------------------------      
         dt            = delt
         cut_cnt       = 0
         stp_con_cnt   = 0
         interval_done = 0.
time_step_loop : &
         do
            dti = 1. / dt
!-----------------------------------------------------------------------      
!        ... transfer from base to class array
!-----------------------------------------------------------------------      
            do k = 1,clscnt4
               j = implicit%clsmap(k)
               m = implicit%permute(k)
               solution(ofl:ofu,m) = base_sol(ofl:ofu,j)
            end do
!-----------------------------------------------------------------------      
!        ... set the iteration invariant part of the function f(y)
!        ... if there is "independent" production put it in the forcing
!-----------------------------------------------------------------------      
            if( implicit%indprd_cnt > 0 .or. extcnt > 0 ) then
               do m = 1,clscnt4
                  iter_invariant(ofl:ofu,m) = dti * solution(ofl:ofu,m) + ind_prd(ofl:ofu,m)
               end do
            else
               do m = 1,clscnt4
                  iter_invariant(ofl:ofu,m) = dti * solution(ofl:ofu,m)
               end do
            end if
!-----------------------------------------------------------------------      
!        ... the linear component
!-----------------------------------------------------------------------      
            if( implicit%lin_rxt_cnt > 0 ) then
               call imp_linmat( ofl, ofu, lin_jac, base_sol, reaction_rates, het_rates )
            end if

!=======================================================================
!        the newton-raphson iteration for f(y) = 0
!=======================================================================
            iter_conv(ofl:ofu)  = .false.
iter_loop : do nr_iter = 1,implicit%iter_max
!-----------------------------------------------------------------------      
!        ... the non-linear component
!-----------------------------------------------------------------------      
               if( factor(nr_iter) ) then
                  call imp_nlnmat( ofl, ofu, sys_jac, base_sol, reaction_rates, &
                                   lin_jac, dti )
!-----------------------------------------------------------------------      
!         ... factor the "system" matrix
!-----------------------------------------------------------------------      
                  call imp_lu_fac( ofl, ofu, sys_jac )
               end if      
!-----------------------------------------------------------------------      
!   	... form f(y)
!-----------------------------------------------------------------------      
               call imp_prod_loss( ofl, ofu, prod, loss, base_sol, &
                                   reaction_rates, het_rates )
               do m = 1,clscnt4
                  forcing(ofl:ofu,m) = solution(ofl:ofu,m)*dti &
                       - (iter_invariant(ofl:ofu,m) + prod(ofl:ofu,m) - loss(ofl:ofu,m))
               end do
!-----------------------------------------------------------------------      
!         ... solve for the mixing ratio at t(n+1)
!-----------------------------------------------------------------------      
               call imp_lu_slv( ofl, ofu, sys_jac, forcing )
               do m = 1,clscnt4
                  where( .not. iter_conv(ofl:ofu) )
                     solution(ofl:ofu,m) = solution(ofl:ofu,m) + forcing(ofl:ofu,m)
                  elsewhere
                     forcing(ofl:ofu,m) = 0.
                  endwhere
               end do
!-----------------------------------------------------------------------      
!    	... convergence measures
!-----------------------------------------------------------------------      
               if( nr_iter > 1 ) then
                  do k = 1,clscnt4
                     m = implicit%permute(k)
                     where( abs(solution(ofl:ofu,m)) > small_val )
                        wrk(ofl:ofu) = abs( forcing(ofl:ofu,m)/solution(ofl:ofu,m) )
                     elsewhere
                        wrk(ofl:ofu) = 0.
                     endwhere
                     max_delta(k) = maxval( wrk(ofl:ofu) )
                  end do
               end if
!-----------------------------------------------------------------------      
!   	... limit iterate
!-----------------------------------------------------------------------      
               do m = 1,clscnt4
                  where( solution(ofl:ofu,m) < 0. )
                     solution(ofl:ofu,m) = 0.
                  endwhere
               end do
!-----------------------------------------------------------------------      
!   	... transfer latest solution back to work array
!-----------------------------------------------------------------------      
               do k = 1,clscnt4
                  j = implicit%clsmap(k)
                  m = implicit%permute(k)
                  base_sol(ofl:ofu,j) = solution(ofl:ofu,m)
               end do
!-----------------------------------------------------------------------      
!    	... check for convergence
!-----------------------------------------------------------------------      
               if( nr_iter > 1 ) then
                  do k = 1,clscnt4
                     m = implicit%permute(k)
                     where( abs( forcing(ofl:ofu,m) ) > small_val )
                        spc_conv(ofl:ofu,k) =  abs(forcing(ofl:ofu,m)) <= epsilon(k)*abs(solution(ofl:ofu,m))
                     elsewhere
                        spc_conv(ofl:ofu,k) =  .true.
                     endwhere
                     converged(k) =  all( spc_conv(ofl:ofu,k) )
                  end do
                  convergence = all( converged(:clscnt4) )
                  if( convergence ) then
                     exit iter_loop
                  end if
                  do k = ofl,ofu
                     if( .not. iter_conv(k) ) then
                        iter_conv(k) = all( spc_conv(k,:) )
                     end if
                  end do
               end if
            end do iter_loop

!-----------------------------------------------------------------------      
!    	... check for newton-raphson convergence
!-----------------------------------------------------------------------      
non_conv:   if( .not. convergence ) then
!-----------------------------------------------------------------------      
!   	... non-convergence
!-----------------------------------------------------------------------      
               if( pdiags%imp_slv ) then
                  write(*,'('' imp_sol: time step '',1p,e21.13,'' failed to converge'')') dt
               end if
               stp_con_cnt = 0
               if( cut_cnt < cut_limit ) then
                  cut_cnt = cut_cnt + 1
                  dt      = .5 * dt
                  do m = 1,pcnstm1
                     base_sol(ofl:ofu,m) = base_sol_sav(ofl:ofu,m)
                  end do
                  cycle time_step_loop
               else
                  write(*,'('' imp_sol: failed to converge @ (lat,lev,dt) = '',2i5,1p,e21.13)') lat,lev,dt
                  do m = 1,clscnt4
                     if( .not. converged(m) ) then
                        write(*,'(1x,a8,1x,1pe10.3)') tracnam(implicit%clsmap(m)), max_delta(m)
                     end if
                  end do
               end if
            end if non_conv
!-----------------------------------------------------------------------      
!   	... check for interval done
!-----------------------------------------------------------------------      
            interval_done = interval_done + dt
            if( abs( delt - interval_done ) <= .0001 ) then
               exit time_step_loop
            else
!-----------------------------------------------------------------------      
!   	... transfer latest solution back to base save array
!-----------------------------------------------------------------------      
               if( convergence ) then
                  stp_con_cnt = stp_con_cnt + 1
               end if
               do m = 1,pcnstm1
                  base_sol_sav(ofl:ofu,m) = base_sol(ofl:ofu,m)
               end do
               if( stp_con_cnt >= 2 ) then
                  dt = 2.*dt
                  stp_con_cnt = 0
               end if
               dt = min( dt,delt-interval_done )
               if( pdiags%imp_slv ) then
                  write(*,'('' imp_sol: new time step '',1p,e21.13)') dt
               end if
            end if
         end do time_step_loop
!-----------------------------------------------------------------------      
!   	... check for history prod, loss
!-----------------------------------------------------------------------      
hist_buff_loop : &
         do file = 1,moz_file_cnt
            do timetype = inst,avrg
               fill_buff = imp_hst_prod(file)%do_hst(timetype)
               if( timetype == inst ) then
                  fill_buff = fill_buff .and. hfile(file)%wrhstts
               end if
has_prod_buff: if( fill_buff ) then
                     hndx = 0
                     do n = 1,hfile(file)%histout_cnt(14,timetype)
                        if( timetype == inst ) then
                           class = hfile(file)%inst_map(hfile(file)%histout_ind(14,inst)+n-1)/1000
                        else
                           class = hfile(file)%timav_map(hfile(file)%histout_ind(14,avrg)+n-1)/1000
                        end if
                        if( class == 4 ) then
                           if( timetype == inst ) then
                              cls_ndx = mod( hfile(file)%inst_map(hfile(file)%histout_ind(14,inst)+n-1),1000 )
                           else
                              cls_ndx = mod( hfile(file)%timav_map(hfile(file)%histout_ind(14,avrg)+n-1),1000 )
                           end if
                           hndx = hndx + 1
                           l = implicit%clsmap(cls_ndx)
                           if( l == ox_ndx ) then
!-----------------------------------------------------------------------      
!         ... ozone production (only valid for the troposphere!)
!-----------------------------------------------------------------------      
                              if( do_ox_pl ) then
                                 prod_buff(file)%buff(ofl:ofu,hndx) = &
                                    (reaction_rates(ofl:ofu,ox_p1_ndx)*base_sol(ofl:ofu,ho2_ndx) &
                                     + reaction_rates(ofl:ofu,ox_p2_ndx) *base_sol(ofl:ofu,ch3o2_ndx) &
                                     + reaction_rates(ofl:ofu,ox_p3_ndx) *base_sol(ofl:ofu,po2_ndx) &
                                     + reaction_rates(ofl:ofu,ox_p4_ndx) *base_sol(ofl:ofu,ch3co3_ndx) &
                                     + reaction_rates(ofl:ofu,ox_p5_ndx) *base_sol(ofl:ofu,c2h5o2_ndx) &
!o3_pl+++
!                                    + .88*reaction_rates(ofl:ofu,ox_p6_ndx)*base_sol(ofl:ofu,isopo2_ndx) &
                                     + .92*reaction_rates(ofl:ofu,ox_p6_ndx)*base_sol(ofl:ofu,isopo2_ndx) &
!o3_pl---
                                     + .985*reaction_rates(ofl:ofu,ox_p7_ndx)*base_sol(ofl:ofu,macro2_ndx) &
                                     + reaction_rates(ofl:ofu,ox_p8_ndx)*base_sol(ofl:ofu,mco3_ndx) & 
                                     + reaction_rates(ofl:ofu,ox_p9_ndx)*base_sol(ofl:ofu,c3h7o2_ndx) & 
                                     + reaction_rates(ofl:ofu,ox_p10_ndx)*base_sol(ofl:ofu,ro2_ndx) &
!o3_pl+++
!                                    + reaction_rates(ofl:ofu,ox_p11_ndx)*base_sol(ofl:ofu,xo2_ndx)) * base_sol(ofl:ofu,no_ndx)
                                     + reaction_rates(ofl:ofu,ox_p11_ndx)*base_sol(ofl:ofu,xo2_ndx) &
                                     + .9*reaction_rates(ofl:ofu,ox_p12_ndx)*base_sol(ofl:ofu,tolo2_ndx) &
                                     +    reaction_rates(ofl:ofu,ox_p13_ndx)*base_sol(ofl:ofu,terpo2_ndx)&
                                     + .9*reaction_rates(ofl:ofu,ox_p14_ndx)*base_sol(ofl:ofu,alko2_ndx) &
                                     +    reaction_rates(ofl:ofu,ox_p15_ndx)*base_sol(ofl:ofu,eneo2_ndx) &
                                     +    reaction_rates(ofl:ofu,ox_p16_ndx)*base_sol(ofl:ofu,eo2_ndx) &
                                     +    reaction_rates(ofl:ofu,ox_p17_ndx)*base_sol(ofl:ofu,meko2_ndx)) * base_sol(ofl:ofu,no_ndx)
!o3_pl---
                              end if
                           else
                              j = implicit%permute(cls_ndx)
                              prod_buff(file)%buff(ofl:ofu,hndx) = prod(ofl:ofu,j) + ind_prd(ofl:ofu,j)
                           end if
                        end if
                     end do
                  end if has_prod_buff
                  fill_buff = imp_hst_loss(file)%do_hst(timetype)
                  if( timetype == inst ) then
                     fill_buff = fill_buff .and. hfile(file)%wrhstts
               end if
has_loss_buff: if( fill_buff ) then
                     hndx = 0
                     do n = 1,hfile(file)%histout_cnt(15,timetype)
                        if( timetype == inst ) then
                           class = hfile(file)%inst_map(hfile(file)%histout_ind(15,inst)+n-1)/1000
                        else
                           class = hfile(file)%timav_map(hfile(file)%histout_ind(15,avrg)+n-1)/1000
                        end if
                        if( class == 4 ) then
                           if( timetype == inst ) then
                              cls_ndx = mod( hfile(file)%inst_map(hfile(file)%histout_ind(15,inst)+n-1),1000 )
                           else
                              cls_ndx = mod( hfile(file)%timav_map(hfile(file)%histout_ind(15,avrg)+n-1),1000 )
                           end if
                           hndx = hndx + 1
                           l = implicit%clsmap(cls_ndx)
                           if( l == ox_ndx ) then
                              if( do_ox_pl ) then
!-----------------------------------------------------------------------      
!         ... ozone destruction (only valid for the troposphere!)
!             also include ox loss from no2+oh, n2o5+aerosol, no3+aerosol
!-----------------------------------------------------------------------      
                                 if( o1d_ndx < 1 ) then
                                    loss_buff(file)%buff(ofl:ofu,hndx) =  reaction_rates(ofl:ofu,ox_l1_ndx)
                                 else
                                    loss_buff(file)%buff(ofl:ofu,hndx) =  reaction_rates(ofl:ofu,ox_l1_ndx) &
                                                                    * base_sol(ofl:ofu,o1d_ndx)/base_sol(ofl:ofu,ox_ndx)
                                 end if
                                 loss_buff(file)%buff(ofl:ofu,hndx) =  loss_buff(file)%buff(ofl:ofu,hndx) &
                                   + reaction_rates(ofl:ofu,ox_l2_ndx) *base_sol(ofl:ofu,oh_ndx) & 
                                   + reaction_rates(ofl:ofu,ox_l3_ndx) *base_sol(ofl:ofu,ho2_ndx) & 
                                   + reaction_rates(ofl:ofu,ox_l6_ndx) *base_sol(ofl:ofu,c2h4_ndx) &
                                   + reaction_rates(ofl:ofu,ox_l4_ndx) *base_sol(ofl:ofu,c3h6_ndx) &
                                   + .9*reaction_rates(ofl:ofu,ox_l5_ndx) *base_sol(ofl:ofu,isop_ndx) &
                                   + .8*(reaction_rates(ofl:ofu,ox_l7_ndx)*base_sol(ofl:ofu,mvk_ndx) &
                                         + reaction_rates(ofl:ofu,ox_l8_ndx)*base_sol(ofl:ofu,macro2_ndx)) &
                                   + .235*reaction_rates(ofl:ofu,ox_l9_ndx)*base_sol(ofl:ofu,c10h16_ndx) &
                                   + (reaction_rates(ofl:ofu,usr4_ndx) * base_sol(ofl:ofu,no2_ndx) * base_sol(ofl:ofu,oh_ndx) &
                                      + 3. * reaction_rates(ofl:ofu,usr16_ndx) * base_sol(ofl:ofu,n2o5_ndx) &
                                      + 2. * reaction_rates(ofl:ofu,usr17_ndx) * base_sol(ofl:ofu,no3_ndx)) &
                                     /max( base_sol(ofl:ofu,ox_ndx),1.e-20 )
                                 loss_buff(file)%buff(ofl:ofu,hndx) = loss_buff(file)%buff(ofl:ofu,hndx) * base_sol(ofl:ofu,l)
                              end if
                           else
                              j = implicit%permute(cls_ndx)
                              loss_buff(file)%buff(ofl:ofu,hndx) = loss(ofl:ofu,j)
                           end if
                        end if
                     end do
               end if has_loss_buff
            end do
         end do hist_buff_loop
      end do spatial_loop

!-----------------------------------------------------------------------      
!    	... check for implicit species lifetime output
!           first check instantaneous then time averaged
!-----------------------------------------------------------------------      
has_lt:     if( moz_file_cnt > 0 .and. lt_cnt > 0 ) then
lifetime_loop : &
               do file = 1,moz_file_cnt
                  if( hfile(file)%wrhstts .and. lt_struct(4,file)%cnt(1) > 0 ) then
                     l = 0
                     do n = 1,hfile(file)%histout_cnt(19,1)
                        hndx  = hfile(file)%histout_ind(19,1)+n-1
                        m     = hfile(file)%inst_map(hndx)
                        if( m/1000 == 4 ) then
                           cls_ndx   = mod( m,1000 )
                           cls_ndxp  = implicit%permute(cls_ndx)
                           base_ndx  = implicit%clsmap(cls_ndx)
                           l = l + 1
                           do lev = 1,plev
                             ofl = (lev - 1)*plonl + 1
                             ofu = ofl + plonl - 1
                             lt_struct(4,file)%lt_mass_inst(1:plonl,lat,ip,l) = &
                               lt_struct(4,file)%lt_mass_inst(1:plonl,lat,ip,l) &
                               + base_sol(ofl:ofu,base_ndx)*pdel(ofl:ofu)
                             lt_struct(4,file)%lt_loss_inst(1:plonl,lat,ip,l) = &
                               lt_struct(4,file)%lt_loss_inst(1:plonl,lat,ip,l) &
                               + pdel(ofl:ofu)*loss(ofl:ofu,cls_ndxp)
                           end do
                        end if
                     end do
                  end if
                  if( lt_struct(4,file)%cnt(2) > 0 ) then
                     l = 0
                     do n = 1,hfile(file)%histout_cnt(19,2)
                        hndx  = hfile(file)%histout_ind(19,2)+n-1
                        m     = hfile(file)%timav_map(hndx)
                        if( m/1000 == 4 ) then
                           cls_ndx   = mod( m,1000 )
                           cls_ndxp  = implicit%permute(cls_ndx)
                           base_ndx  = implicit%clsmap(cls_ndx)
                           l = l + 1
                           do lev = 1,plev
                             ofl = (lev - 1)*plonl + 1
                             ofu = ofl + plonl - 1
                             lt_struct(4,file)%lt_mass_avrg(1:plonl,lat,ip,l) = &
                               lt_struct(4,file)%lt_mass_avrg(1:plonl,lat,ip,l) &
                               + base_sol(ofl:ofu,base_ndx)*pdel(ofl:ofu)
                             lt_struct(4,file)%lt_loss_avrg(1:plonl,lat,ip,l) = &
                               lt_struct(4,file)%lt_loss_avrg(1:plonl,lat,ip,l) &
                               + pdel(ofl:ofu)*loss(ofl:ofu,cls_ndxp)
                           end do
                        end if
                     end do
                  end if
               end do lifetime_loop
            end if has_lt

!-----------------------------------------------------------------------      
!    	... check for implicit species production and loss output
!           first check production; instantaneous then time averaged
!           then  check loss; instantaneous then time averaged
!-----------------------------------------------------------------------      
      do file = 1,moz_file_cnt
         do timetype = inst,avrg
            dump_buff = imp_hst_prod(file)%do_hst(timetype)
            if( timetype == inst ) then
               dump_buff = dump_buff .and. hfile(file)%wrhstts
            end if
            if( dump_buff ) then
               hndx = 0
               do n = 1,hfile(file)%histout_cnt(14,timetype)
                  if( timetype == inst ) then
                     class = hfile(file)%inst_map(hfile(file)%histout_ind(14,timetype)+n-1)/1000
                  else
                     class = hfile(file)%timav_map(hfile(file)%histout_ind(14,timetype)+n-1)/1000
                  end if
                  if( class == 4 ) then
                     if( timetype == inst ) then
                        cls_ndx = mod( hfile(file)%inst_map(hfile(file)%histout_ind(14,inst)+n-1),1000 )
                        fldname = hfile(file)%hist_inst(hfile(file)%histout_ind(14,timetype)+n-1)
                     else
                        cls_ndx = mod( hfile(file)%timav_map(hfile(file)%histout_ind(14,avrg)+n-1),1000 )
                        fldname = hfile(file)%hist_timav(hfile(file)%histout_ind(14,timetype)+n-1)
                     end if
                     hndx        = hndx + 1
                     wrk_buff(:) = prod_buff(file)%buff(:,hndx) * hnm(:)
                     call outfld( fldname, wrk_buff, plonl, ip, lat, file )
                  end if
               end do
            end if
            dump_buff = imp_hst_loss(file)%do_hst(timetype)
            if( timetype == inst ) then
               dump_buff = dump_buff .and. hfile(file)%wrhstts
            end if
            if( dump_buff ) then
               hndx = 0
               do n = 1,hfile(file)%histout_cnt(15,timetype)
                  if( timetype == inst ) then
                     class = hfile(file)%inst_map(hfile(file)%histout_ind(15,timetype)+n-1)/1000
                  else
                     class = hfile(file)%timav_map(hfile(file)%histout_ind(15,timetype)+n-1)/1000
                  end if
                  if( class == 4 ) then
                     if( timetype == inst ) then
                        cls_ndx = mod( hfile(file)%inst_map(hfile(file)%histout_ind(15,inst)+n-1),1000 )
                        fldname = hfile(file)%hist_inst(hfile(file)%histout_ind(15,timetype)+n-1)
                     else
                        cls_ndx = mod( hfile(file)%timav_map(hfile(file)%histout_ind(15,avrg)+n-1),1000 )
                        fldname = hfile(file)%hist_timav(hfile(file)%histout_ind(15,timetype)+n-1)
                     end if
                     hndx        = hndx + 1
                     l           = implicit%clsmap(cls_ndx)
                     wrk_buff(:) = loss_buff(file)%buff(:,hndx) * hnm(:)
                     call outfld( fldname, wrk_buff, plonl, ip, lat, file )
                  end if
               end do
            end if
         end do
      end do

      if( allocated( prod_buff ) ) then
         do file = 1,moz_file_cnt
            if( associated( prod_buff(file)%buff ) ) then
               deallocate( prod_buff(file)%buff )
            end if
         end do
         deallocate( prod_buff )
      end if
      if( allocated( loss_buff ) ) then
         do file = 1,moz_file_cnt
            if( associated( loss_buff(file)%buff ) ) then
               deallocate( loss_buff(file)%buff )
            end if
         end do
         deallocate( loss_buff )
      end if

      end subroutine imp_sol

      end module mo_imp_sol