ecearth3/sources/tm5mp/base/convection.F90

!### 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
  !   17 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !---------------------------------------------------------------------------

CONTAINS
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CONVECTION_INIT
  !
  ! !DESCRIPTION: set %used flag of releveant meteo to True
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CONVECTION_INIT( status )
    !
    ! !USES:
    !
    use dims       , only : nregions
    use Meteo      , only : Set
    use Meteodata      , only : entu_dat, entd_dat, detu_dat, detd_dat
    use Meteodata      , only : mfw_dat, omega_dat
#ifndef without_wet_deposition
    use Meteodata      , only : cp_dat
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  ::  status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/Convection_Init'
    
    integer           ::  region
    
    ! 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
  !EOC

  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CONVECTION_DONE
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CONVECTION_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)           ::  status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/Convection_Done'

    ! nothing to do

    status = 0

  END SUBROUTINE CONVECTION_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CONVEC
  !
  ! !DESCRIPTION: 
  !
  !****   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
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CONVEC( region, status )
    !
    ! !USES:
    !
    use dims,          only : tref !, dx, dy, xref, yref 
    use dims,          only : itau, okdebug, kdebug, nconv, lmax_conv
    use dims,          only : nregions
    use dims,          only : revert
    use global_data,   only : region_dat, conv_dat, mass_dat
    use meteodata,     only : entu_dat, entd_dat, detu_dat, detd_dat
    use meteodata,     only : cp_dat
    use meteodata,     only : m_dat
    use chem_param,    only : ra, ntracet
#ifdef with_budgets    
#ifndef without_wet_deposition
    use budget_global, only : nzon_vg
    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, plc_ipr_lwdep, plc_ipr_tcnvd, plc_kg_from_tm
#endif
    use tm5_distgrid,  only : dgrid, Get_DistGrid, reduce
    use partools,      only : isRoot
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in)     ::  region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,intent(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !   17 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter  ::  rname = mname//'/CONVEC'

    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, i1, i2, j1, j2
    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

    ! for global wetdep budget of tracer #1
    real              :: g_sum_wet
    real, allocatable :: l_sum_wet(:,:)

    !====================== start ==================================
    
    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;;

    !-------- local indices
    
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

#ifdef with_budgets    
#ifndef without_wet_deposition
    allocate(l_sum_wet(i1:i2,j1:j2))
    l_sum_wet = 0.0
#endif
#endif    
    
    !-------- shortcuts
    
    m    => m_dat(region)%data
    rm   => mass_dat(region)%rm

#ifdef slopes
    rxm  => mass_dat(region)%rxm
    rym  => mass_dat(region)%rym
    rzm  => mass_dat(region)%rzm
#ifdef secmom
    rxxm => mass_dat(region)%rxxm
    rxym => mass_dat(region)%rxym
    rxzm => mass_dat(region)%rxzm
    ryym => mass_dat(region)%ryym
    ryzm => mass_dat(region)%ryzm
    rzzm => mass_dat(region)%rzzm
#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) &
#endif
    !$OMP  shared  (region, ntracetloc, sceffdeep, dt,   &
    !$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, li, ld, n, 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)
    
    do j = j1, j2
    do i = i1, i2

#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
          TRACER : do n = 1, ntracet

             ! 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
                srxm(l)=conv1(l,l)*srxm(l) !no convection -> no multip.  !WP-Emma van der Veen
                srym(l)=conv1(l,l)*srym(l) !no convection -> no multip.  !WP-Emma van der Veen
#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(n), 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
!WP-vdV!           zrxm (l) = zrxm (l) + (conv1(l,k)-scavf) * srxm(k)
!WP-vdV!           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
!WP-vdV!                rxm(i,j,l,n)=zrxm(l)
!WP-vdV!                rym(i,j,l,n)=zrym(l)
                rxm(i,j,l,n)=srxm(l) !VDV
                rym(i,j,l,n)=srym(l) !VDv
                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(n).gt.0.) then
                do l=1,lmax_conv
                   nzv=nzon_vg(l)
                   buddep_dat(region)%cp(i,j,nzv,n) = buddep_dat(region)%cp(i,j,nzv,n) + &
                        zwetrm(l)/ra(n)*1.e3       !kg=> moles
                   if ( n == 1 ) l_sum_wet(i,j) = l_sum_wet(i,j) + 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( region, plc_ipr_tcnvd, i, j, l, n, &
                     zcnvd(l)*plc_kg_from_tm(n), status )
                IF_NOTOK_RETURN(status=1)
                call PLC_AddValue( region, plc_ipr_lwdep, i, j, l, n, &
                     zwetrm(l)*plc_kg_from_tm(n), status )
                IF_NOTOK_RETURN(status=1)
             end do
#endif


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

    !$OMP END PARALLEL


#ifndef without_wet_deposition
#ifdef with_budgets

    call REDUCE( dgrid(region), l_sum_wet, 0, 'SUM', g_sum_wet, status)
    IF_NOTOK_RETURN(status=1)
    
    if ( isRoot ) sum_wetdep(region) = sum_wetdep(region) + g_sum_wet

    deallocate( l_sum_wet )
#endif
#endif

    !--- Cleanup
    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
  !EOC



#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