!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!#################################################################

!------------------------------------------------------------------------------
!                    TM5                                                      !
!------------------------------------------------------------------------------
!BOP
!
! !MODULE:  convection
!
! !DESCRIPTION: methods to init and apply convection
!\\
!\\
! !INTERFACE: 
!
module convection
!
! !USES:
!
  use GO, only : gol, goPr, goErr

  implicit none

  private
!
! !PUBLIC MEMBER FUNCTIONS:
!
  public  ::  Convection_Init, Convection_Done
  public  ::  Convec
!
! !PRIVATE DATA MEMBERS:
!
  character(len=*), parameter  ::  mname = 'convection'
!
! !REVISION HISTORY: 
!   30 Mar 2010 - P. Le Sager - re-instate test on wet dep flag in convec
!
! !REMARKS:
!
!EOP
!------------------------------------------------------------------------------

contains



  ! ================================================================
  
  
  subroutine Convection_Init( status )
  
    use dims       , only : nregions
    use Meteo      , only : Set
    use Meteo      , only : entu_dat, entd_dat, detu_dat, detd_dat
    use Meteo      , only : mfw_dat, omega_dat
#ifndef without_wet_deposition
    use Meteo      , only : cp_dat
#endif

    ! --- in/out --------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Convection_Init'
    
    ! --- local -------------------------------------
    
    integer           ::  region
    
    ! --- begin --------------------------------
    
    ! loop over all (zoom) regions:
    do region = 1, nregions
    
      ! entrements/detrements
      call Set( entu_dat(region), status, used=.true. )
      call Set( entd_dat(region), status, used=.true. )
      call Set( detu_dat(region), status, used=.true. )
      call Set( detd_dat(region), status, used=.true. )

#ifndef without_wet_deposition
      ! convective precip:
      call Set(   cp_dat(region), status, used=.true. )
#endif
    
      ! omega is used when creating updr/downdr from raw data;
      ! computed from mfw :
      call Set(   mfw_dat(region), status, used=.true. )
      call Set( omega_dat(region), status, used=.true. )
    
    end do
    
    ! ok
    status = 0
    
  end subroutine Convection_Init
  
  
  ! ***
  
  
  subroutine Convection_Done( status )
  
    ! --- in/out --------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Convection_Done'

    ! --- begin --------------------------------
    
    ! nothing to be done

    ! ok
    status = 0
    
  end subroutine Convection_Done
  
  
  
  ! ================================================================
  

  !-----------------------------------------------------------------------
  !
  !****   convec          - mix tracers by convection     v 8.5
  !
  !       programmed by           mh      mpi HH          1-oct-1991
  !       modified by             mh      mpi HH         11-jun-1994
  !
  !       purpose
  !       -------
  !       mix tracers vertically by convection
  !
  !       interface
  !       ---------
  !       call convec
  !
  !       method
  !       ------
  !       the matrix conv is applied simultaneously to rm, rxm and rym. 
  !       The diagonal elements of conv are applied to rzm.
  ! 
  !       externals
  !       ---------
  !       subroutines:    tstamp
  !
  !       reference
  !       ---------       
  !       see manual
  !-----------------------------------------------------------------------
  !       Following the method by Walter Guelle (1997) convective removal of 
  !        tracers is added with some little changes:        
  !       a vector cvsfac is added containing scavenging efficiencies (0-1)
  !       per tracer.
  !                                   aj, knmi may 1998
  !      Maarten Krol, july 2000, implemented zoom version
  !      lmax_conv is now set as the maximum layer to which convection 
  !      is taken into account.
  !
  !       Second moments have been added: Bram Bregman, August 2004
  !--------------------------------------------------------------

  subroutine convec(region, status)

    use dims,          only : dx, dy, xref, yref, tref
    use dims,          only : isr, ier, jsr, jer
    use dims,          only : im, jm, lm
    use dims,          only : itau, okdebug, kdebug, nconv, lmax_conv
    use dims,          only : nregions
    use dims,          only : revert
    use global_data,   only : mass_dat,region_dat
    use global_data,   only : emis_data
    use global_data,   only : conv_dat
    use meteo         , only : entu_dat, entd_dat, detu_dat, detd_dat
    use meteo         , only : cp_dat
    use meteodata     , only : m_dat
    use chem_param,    only : ra, ntracet
#ifdef with_budgets    
    use budget_global, only : nzon_vg, apply_budget_global
#ifndef without_wet_deposition
    use wet_deposition,only : sum_wetdep, buddep_dat
#endif
#endif
    use toolbox,       only : lvlpress
#ifndef without_wet_deposition
    use wet_deposition,only : cvsfac, cp_scale
#endif
    use datetime,      only : tstamp
#ifdef with_tendencies  
    use tracer_data  , only : PLC_AddValue_t, plc_ipr_lwdep, plc_ipr_tcnvd, plc_kg_from_tm
#endif
#ifdef MPI
    use mpi_const, only: my_real, &
         mpi_min,mpi_max,com_trac,ierr
    use mpi_const, only: pe_first_tracer
#endif
    use ParTools, only : myid, root_t
    use ParTools, only : ntracetloc, ntracet_ar
    use omp_ParTools, only : my_omp_domain, my_omp_sum

    implicit none

    ! input/output
    integer,intent(in)     ::  region
    integer,intent(out)    ::  status

    ! local
    real,dimension(:,:,:,:), pointer  :: rm
#ifdef slopes
    real,dimension(:,:,:,:), pointer  :: rxm,rym,rzm
#ifdef secmom
    real,dimension(:,:,:,:), pointer  :: rxxm,rxym,rxzm,ryym,ryzm,rzzm
#endif
#endif
    real,dimension(:,:,:)  , pointer  :: m
    real,dimension(:,:,:)  , pointer  :: entu, entd, detu, detd
#ifndef without_diffusion
    real,dimension(:,:,:)  , pointer  :: dkg
#endif
    integer,dimension(:,:) , pointer  :: cloud_base, cloud_top, cloud_lfs, zoomed
    integer,dimension(3)              :: mxloc
    ! arrays for convection matrix
    ! conv1: convection matrix, lbdcv: wet removal matrix
    real,dimension(lmax_conv,lmax_conv)     :: conv1
    real,dimension(lmax_conv,lmax_conv)     :: lbdcv
    real,dimension(0:lmax_conv,lmax_conv)   :: f
    real,dimension(0:lmax_conv,lmax_conv)   :: fu
    real,dimension(0:lmax_conv)             :: amu,amd
    real,dimension(lmax_conv)               :: new_mass
    real,dimension(lmax_conv)               :: zwetrm, zcnvd
    real,dimension(lmax_conv)               :: zrm
    real,dimension(lmax_conv)               :: srm
#ifdef slopes
    real,dimension(lmax_conv)               :: zrxm,zrym
    real,dimension(lmax_conv)               :: srxm,srym,srzm
#ifdef secmom
    real,dimension(lmax_conv)               :: zrxxm,zrxym,zrxzm,zryym,zryzm
    real,dimension(lmax_conv)               :: srxxm,srxym,srxzm,sryym,sryzm,srzzm
#endif
#endif
    integer :: n,l,k,j,i
    integer :: imr, jmr, lmr
    integer :: info,ld,li,kk, lshallowtop, nzv
    real    :: scavf,dt,zxi,wetpct
    real    :: minvalue, mxval,mnval,mxvald
    real    :: minvalue_all, mxval_all,mnval_all,mxvald_all
    real    :: sceffdeep, cp_value
    integer :: nglob, offsetn
    ! Ompenmp parameters
    integer :: jss,jee
    real    :: l_sum_wet, g_sum_wet

    ! start
    imr=im(region) ; jmr=jm(region) ; lmr=lm(region)

    !WP! only PE's with nonzero ntracet
    if ( ntracetloc == 0 ) return  

    offsetn = sum(ntracet_ar(0:myid-1))

    sceffdeep = 1.0
    ! rain scale on this resolution: e.g. on 6x4 this read 1x1/(6x4)
    !    thus 1 mm/hr becomes 1mm/day
    !  cp_scaler = cp_scale*xref(region)*yref(region)/(dx*dy)  
    !CMK removed after test study: hardly resolution dependent;;

    m => m_dat(region)%data
    rm => mass_dat(region)%rm_t
#ifdef slopes
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t
#ifdef secmom
    rxxm => mass_dat(region)%rxxm_t
    rxym => mass_dat(region)%rxym_t
    rxzm => mass_dat(region)%rxzm_t
    ryym => mass_dat(region)%ryym_t
    ryzm => mass_dat(region)%ryzm_t
    rzzm => mass_dat(region)%rzzm_t
#endif
#endif
    zoomed => region_dat(region)%zoomed
    entu => entu_dat(region)%data
    detu => detu_dat(region)%data
    entd => entd_dat(region)%data
    detd => detd_dat(region)%data
    cloud_top => conv_dat(region)%cloud_top
    cloud_base => conv_dat(region)%cloud_base
    cloud_lfs => conv_dat(region)%cloud_lfs
#ifndef without_diffusion
    dkg => conv_dat(region)%dkg
#endif

    if ( okdebug ) call tstamp(kdebug,itau,'convec ')
    if ( okdebug ) print *,'convec: Convection called (region)', region

    !timestep and position in packed arrays
    dt = nconv/(2.*tref(region))     

    if ( okdebug ) then
       minvalue = 2.0
       mxvald = -1.0
       mxval = 0.0
       mnval = 2.0
    end if

    ! start main loop over surface cells
!
! PLS : re-incorporate the test on wet dep, from older cy3base/src.
!       It hides cvsfac and cp_scale when wet dep is off. 
!$OMP PARALLEL &
!$OMP  default (none) &
#ifdef with_budgets
!$OMP  shared  (nzon_vg, apply_budget_global) &
#endif
!$OMP  shared  (region, ntracetloc, sceffdeep, dt, offsetn,  &
!$OMP           g_sum_wet, cloud_base, cloud_top, cloud_lfs, cp_dat, &
!$OMP           zoomed, rm, rxm, rym, rzm, m, entu, entd, detu, detd, &
!$OMP           dkg, revert, isr, ier, jsr, jer) &
#ifndef without_wet_deposition
!$OMP shared   ( cvsfac,cp_scale ) &
#ifdef with_budgets
!$OMP shared   ( buddep_dat,ra ) &
#endif
#endif
!$OMP  private (i, j, k,kk,l,n, jss, jee, li, ld, nglob, nzv, zxi, scavf, &
!$OMP           l_sum_wet, conv1, lbdcv, f, fu,amu, amd, zrm, zrxm, &
!$OMP           zrym, zwetrm, zcnvd, srm, srxm, srym, srzm, cp_value)
    call my_omp_domain (jsr(region), jer(region), jss, jee)
    l_sum_wet = 0.0
    do  j=jss,jee
       do  i=isr(region),ier(region)
#ifdef with_zoom
          if (zoomed(i,j) /= region) cycle
#endif
          ! calculate convection matrix and directly apply to rm(i,j,*,1:ntracet)
          ! Set all amu's, amd's entu' etc to 0 for gridcells without clouds
          !
          ! updraft   - amu is positive upward
          !   
          f(:,:)=0.
          fu(:,:)=0.
          amu(:)=0.
          li = cloud_top(i,j)

          do k=1,li
             amu(k)=amu(k-1)+entu(i,j,k)-detu(i,j,k)
             if ( amu(k) > 0. ) then
                !
                ! we limit zxi from below at 0. in case 
                ! there are inconsistent en/detrainment
                ! rates. The case zxi>1 should not occur
                ! since eu and du are >=0.
                !
                zxi=max(0.,1.-detu(i,j,k)/(amu(k-1)+entu(i,j,k)))

             else
                ! set massflux at upper boundary to strictly 0.

                amu(k)=0.
                zxi=0.
             end if
             do kk=1,k-1
                fu(k,kk)=fu(k-1,kk)*zxi        !all previous levels
             end do
             !cmk: note: the fu(k-1,k) values are zero (see manual Heimann)
             fu(k,k)=entu(i,j,k)/m(i,j,k)*zxi
          end do !k
          !
          ! downdraft  -  amd is negative downward
          !
          amd(:)=0.
          ld = cloud_lfs(i,j)
          do  k=ld,2,-1   
             amd(k-1)=amd(k)-entd(i,j,k)+detd(i,j,k)
             if ( amd(k-1) < 0. ) then
                zxi=max(0.,1.+detd(i,j,k)/(amd(k)-entd(i,j,k)))
             else
                amd(k-1)=0.
                zxi=0.
             end if
             do  kk=k+1,ld   !
                f(k-1,kk)=f(k,kk)*zxi !f is here only fd downdraftmatrix
             end do !kk
             f(k-1,k)=-entd(i,j,k)/m(i,j,k)*zxi    !note the negatives for j  
          end do !k  

          !
          ! add coefficients from up and downdraft together
          !
          do k=1,lmax_conv-1
             do kk=1,lmax_conv
                f(k,kk)=fu(k,kk)+f(k,kk)
             end do !kk
             ! add coefficient from subsidence
             !CMK SH       f(k,k+1)=f(k,k+1)-(amu(k)+amd(k))/m(i,j,k+1)
             ! reported by Sander Houweling.........
             f(k,k+1)=f(k,k+1)-amu(k)/m(i,j,k+1)
             f(k,k  )=f(k,k  )-amd(k)/m(i,j,k  )
             !
#ifndef without_diffusion
             ! diffusion terms
             ! no diffusion in the stratosphere  k>lmax_conv?
             f(k,k  ) = f(k,k  ) + dkg(i,j,k)/m(i,j,k  )
             f(k,k+1) = f(k,k+1) - dkg(i,j,k)/m(i,j,k+1)
#endif
          end do
          
          ! >>> wet removal >>>

          ! initialize wet-removal
          lbdcv(:,:) = 0.0

#ifndef without_wet_deposition
          ! fill scaveging stuff if necessary

          !
          ! WARNING: For matrices f(k,kk), fu(k,kk), and fd(k,kk), index k 
          !          corresponds to the upper boundary of level k 
          !          and index j to level j.
          !          For matrices g(i,j,k,kk) and lbdcv(i,j,k,kk), 
          !          both indices k and j
          !          correspond to respective levels and not to boundaries.
          !
          cp_value = cp_dat(region)%data(i,j,1)   ! rain in m/s in region
          do k=2,lmax_conv
             do kk=1,k-1
                !cmk lbdcv(k,kk)=sceffdeep*cp_eff(region)%surf(i,j)* &
                !cmk             dt*(fu(k-1,kk) - fu(k,kk)) 
                ! scale factor between 0 - 1, with 1 for cp > cp_scaler
                !CMK OLD lbdcv(k,kk)=sceffdeep*min(cp_value, cp_scaler)/cp_scaler * &
                lbdcv(k,kk)=sceffdeep*(1.0 - exp(-cp_value/cp_scale))* &
                            dt*(fu(k-1,kk) - fu(k,kk)) 
             enddo !kk
          enddo !k
          ! CMK end
          !
          ! which ensures that the rate of scavenging is equal to dn/dt=-s Mu n
          !
          ! TvN (8 July 2004)

#endif
          
          ! <<<

          !---- WG end 21/04/98 -----

          ! generate forward matrix
          do k=1,lmax_conv
             do kk=1,lmax_conv
                conv1(k,kk)=-dt*(f(k-1,kk)-f(k,kk))
             end do !kk
             conv1(k,k) = conv1(k,k) + 1.0
          end do !k
          call fastminv(conv1,lmax_conv)

          ! some checks...can be removed later!

          !if(okdebug) then
          !    minvalue = min(minvalue,minval(conv1))
          !    do l=1,lmax_conv
          !        mxval = max(mxval,sum(conv1(:,l)))
          !        mnval = min(mnval,sum(conv1(:,l)))
          !    end do
          !
          !    ! test deviations for a uniform mixing ratio of 1.....ppb
          !
          !    do l=1,lmax_conv
          !       new_mass(l) = dot_product(conv1(l,:),m(i,j,1:lmax_conv)*1e-9)
          !    end do
          !    mxvald = max(mxvald, maxval(abs((m(i,j,1:lmax_conv)*1e-9-new_mass)
          !                 /(m(i,j,1:lmax_conv)*1e-9))))
          !end if   !okdebug
          !
          ! directly apply the matrix:
          ! copy i,j

          if(revert==-1) conv1 = transpose(conv1)

          ! loop over local tracers
          do n = 1, ntracetloc

                ! global tracer index:
                nglob=n+offsetn  

                ! copy tracer concentrations for this column in srm/srxm/etc:
                do l=1,lmax_conv
                   srm(l) = rm(i,j,l,n)
#ifdef slopes
                   srxm(l) = rxm(i,j,l,n)
                   srym(l) = rym(i,j,l,n)
                   srzm(l) = rzm(i,j,l,n)
#ifdef secmom
                   srxxm(l) = rxxm(i,j,l,n)
                   srxym(l) = rxym(i,j,l,n)
                   srxzm(l) = rxzm(i,j,l,n)
                   sryym(l) = ryym(i,j,l,n)
                   sryzm(l) = ryzm(i,j,l,n)
                   srzzm(l) = rzzm(i,j,l,n)
#endif
#endif
                end do

                ! init budgets to zero:
                do l=1,lmax_conv
                   zwetrm(l) = 0.0
                   zcnvd (l) = 0.0
                end do

                ! init column of concentrations to zero:
                do l=1,lmax_conv
                   zrm(l)=0.  
#ifdef slopes
                   zrxm(l) = 0.0
                   zrym(l) = 0.0
#ifdef secmom
                   zrxxm(l) = 0.0
                   zrxym(l) = 0.0
                   zrxzm(l) = 0.0
                   zryym(l) = 0.0
                   zryzm(l) = 0.0
#endif
#endif
                end do

                ! apply convection matrix:
                do l=1,lmax_conv
#ifdef slopes
                   srzm(l)=conv1(l,l)*srzm(l)   !rzm is directly calculated
#ifdef secmom
                   srzzm(l)=conv1(l,l)*srzzm(l)   !rzzm is directly calculated
#endif
#endif
                   do k=1,lmax_conv

                      ! species dependent removal....
#ifndef without_wet_deposition
                      scavf = min( lbdcv(l,k)*cvsfac(nglob), conv1(l,k) ) 
#else
                      scavf = 0.0
#endif

                      ! for budgets:
                      zcnvd (l) = zcnvd (l) +  conv1(l,k)        * srm(k)
                      zwetrm(l) = zwetrm(l) +             scavf  * srm(k)  ! loss is positive

                      ! store new concentrations in zrm/zrxm/etc
                      !  using old concentrations in srm/srxm/etc
                      zrm(l)    = zrm(l)    + (conv1(l,k)-scavf) * srm(k)
#ifdef slopes
                      zrxm (l) = zrxm (l) + (conv1(l,k)-scavf) * srxm(k)
                      zrym (l) = zrym (l) + (conv1(l,k)-scavf) * srym(k)
#ifdef secmom
                      zrxxm(l) = zrxxm(l) + (conv1(l,k)-scavf) * srxxm(k)
                      zrxym(l) = zrxym(l) + (conv1(l,k)-scavf) * srxym(k)
                      zrxzm(l) = zrxzm(l) + (conv1(l,k)-scavf) * srxzm(k)
                      zryym(l) = zryym(l) + (conv1(l,k)-scavf) * sryym(k)
                      zryzm(l) = zryzm(l) + (conv1(l,k)-scavf) * sryzm(k)
#endif
#endif
                   end do    !k

                   ! substract original tracer mass
                   ! VH, TvN (2007-11-16)
                   zcnvd(l)=zcnvd(l)-srm(l)

                end do      !l

                ! restore concentrations:
                do l=1,lmax_conv
                   rm(i,j,l,n)=zrm(l)
#ifdef slopes
                   rxm(i,j,l,n)=zrxm(l)
                   rym(i,j,l,n)=zrym(l)
                   rzm(i,j,l,n)=srzm(l)
#ifdef secmom
                    rxxm(i,j,l,n)=zrxxm(l)
                    rxym(i,j,l,n)=zrxym(l)
                    rxzm(i,j,l,n)=srxzm(l)
                    ryym(i,j,l,n)=zryym(l)
                    ryzm(i,j,l,n)=zryzm(l)
                    rzzm(i,j,l,n)=srzzm(l)
#endif
#endif
                end do

#ifdef with_budgets    
#ifndef without_wet_deposition
                  if (cvsfac(nglob).gt.0.) then
                     do l=1,lmax_conv
                        nzv=nzon_vg(l)
                        buddep_dat(region)%cp(i,j,nzv,nglob) = &
                             buddep_dat(region)%cp(i,j,nzv,nglob) + &
                             zwetrm(l)/ra(nglob)*1.e3       !kg=> moles
                        if ( nglob == 1 ) l_sum_wet = l_sum_wet + &
                             zwetrm(l)  !in kg!
                     end do !l
                  end if
#endif
#endif

#ifdef with_tendencies  
                ! Add conv/diff and wet deposition budgets in kg to tendencies:
                do l = 1, lmax_conv
                  call PLC_AddValue_t( region, plc_ipr_tcnvd, i, j, l, nglob, &
                                         zcnvd(l)*plc_kg_from_tm(nglob), status )
                  IF_NOTOK_RETURN(status=1)
                  call PLC_AddValue_t( region, plc_ipr_lwdep, i, j, l, nglob, &
                                         zwetrm(l)*plc_kg_from_tm(nglob), status )
                  IF_NOTOK_RETURN(status=1)
                end do
#endif


             end do          !ntracet
          end do   !i
       end do     !j

#if defined (with_budgets)
       call my_omp_sum (l_sum_wet, g_sum_wet)
#endif

!$OMP END PARALLEL


#ifndef without_wet_deposition
#if defined (with_budgets)
       sum_wetdep(region) = sum_wetdep(region) + g_sum_wet
#endif
#endif

       !if (okdebug)  then
       !
       !   call mpi_allreduce(minvalue,minvalue_all,1,my_real, &
       !        mpi_min,com_trac,ierr)
       !   call mpi_allreduce(mxval,mxval_all,1,my_real,mpi_max,com_trac,ierr)
       !   call mpi_allreduce(mnval,mnval_all,1,my_real,mpi_min,com_trac,ierr)
       !   call mpi_allreduce(mxvald,mxvald_all,1,my_real,mpi_max,com_trac,ierr)
       !   call barrier_t
       !   if (myid==root_t) then
       !    print* ,'-----convection matrix information----------------------'
       !    print* ,'Minimum value of convection matrix :',minvalue_all
       !    print *,'Max/Min along the columns  (?/=1)  :',mxval_all,mnval_all  
       !     !these should be one/zero to ensure mass-conservation...
       !    print *,'Maximum fractional deviation from uniform mixing ', &
       !       'ratio column:',mxvald_all  
       !    print* ,'------convection matrix information---------------------'
       !   end if
       !   call barrier_t
       !   if (mxvald_all > 0.1) &
       !      call exitus(1, ' Convection matrix appears wrong') 
       !end if


       nullify(m)
       nullify(rm)
#ifdef slopes
       nullify(rxm)
       nullify(rym)
       nullify(rzm)
#ifdef secmom
       nullify(rxxm)
       nullify(rxym)
       nullify(rxzm)
       nullify(ryym)
       nullify(ryzm)
       nullify(rzzm)
#endif
#endif
       nullify(zoomed)
       nullify(entu)
       nullify(detu)
       nullify(entd)
       nullify(detd)
       nullify(cloud_top)
       nullify(cloud_base)
       nullify(cloud_lfs)
#ifndef without_diffusion
       nullify(dkg)
#endif

     ! ok
     status = 0

     end subroutine convec


#ifdef with_lapack

     !-----------------------------------------------------------------------
     ! compute inverse matrix using LU decompositon and forward-backward
     ! substitution without pivoting (for easy vectorization)
     !
     ! LAPACK version, WP, Jan-2006
     !-----------------------------------------------------------------------

     subroutine fastminv(a,lm)

       use toolbox,       only: escape_tm

       ! input/output
       integer,intent(in)                  :: lm
       real,dimension(lm,lm),intent(inout) :: a

       ! local
       integer :: istat
       real,dimension(lm) :: work
       real,dimension(lm) :: ipiv

       call DGETRF( lm, lm, a, lm, ipiv, istat ) 
       if(istat.ne.0) call escape_tm('fastminv: Failed to make LU decomposition ')
       call DGETRI( lm, a, lm, ipiv, work, lm, istat )
       if(istat.ne.0) call escape_tm('fastminv: Failed to make inverse matrix')

     end subroutine fastminv

#else

     !-----------------------------------------------------------------------
     ! compute inverse matrix using LU decompositon and forward-backward
     ! substitution without pivoting (for easy vectorization)
     !
     ! Source: EMR and NR, sec. edition (p. 35ff).
     !
     !                                        mh, {Sat, Jun 4, 1994} 16:27:15 
     !-----------------------------------------------------------------------

     subroutine fastminv(a,lm)

       ! input/output
       integer,intent(in)                  :: lm
       real,dimension(lm,lm),intent(inout) :: a

       ! local
       real,dimension(lm)    :: y
       real,dimension(lm)    :: b
       real,dimension(lm,lm) :: bb
       integer :: i,j,k

       do j=1,lm

          do i=1,j
             bb(i,j)=a(i,j)
             do k=1,i-1
                bb(i,j)=bb(i,j)-bb(i,k)*bb(k,j)
             end do
          end do

          do i=j+1,lm
             bb(i,j)=a(i,j)
             do k=1,j-1
                bb(i,j)=bb(i,j)-bb(i,k)*bb(k,j)
             end do
             bb(i,j)=bb(i,j)/bb(j,j)
          end do
       end do

       do k=1,lm

          do i=1,lm
             b(i)=0.
          end do
          b(k)=1.
          do i=1,lm
             y(i)=b(i)
             do j=1,i-1
                y(i)=y(i)-bb(i,j)*y(j)
             end do
          end do

          do i=lm,1,-1
             a(i,k)=y(i)
             do j=i+1,lm
                a(i,k)=a(i,k)-bb(i,j)*a(j,k)
             end do
             a(i,k)=a(i,k)/bb(i,i)
          end do

       end do


     end subroutine fastminv

#endif


   end module convection
