!
!  $Author: pkubota $
!  $Date: 2007/02/15 17:03:10 $
!  $Revision: 1.5 $
!
MODULE Shall_Tied

USE Constants, ONLY :  &
       cp                 , &
       hl                 , &
       gasr               , &
       grav               , &
       rmwmd              , &
       rmwmdi             , &
       e0c                , &
       delq               , &
       p00                , &
       tov                , &
       r8
USE Options, ONLY :       &
       rccmbl            , &
       sthick            , &
       sacum             , &
       acum0             , &
       tbase             , &
       ki                , &
       mlrg              , &
       is                , &
       doprec            , &
       cflric            , &
       ifilt             , &
       dt                , &
       kt                , &
       ktp               , &
       jdt               , &
       nfprt

IMPLICIT NONE
  PRIVATE
PUBLIC :: InitShall_Tied
PUBLIC :: shalv2

  REAL(KIND=r8) :: aa(15)
  REAL(KIND=r8) :: ad(15)
  REAL(KIND=r8) :: ac(15)
  REAL(KIND=r8) :: actop
  REAL(KIND=r8) :: thetae(151,181)
  REAL(KIND=r8) :: tfmthe(431,241)
  REAL(KIND=r8) :: qfmthe(431,241)
  REAL(KIND=r8) :: ess
  INTEGER :: kbase
  INTEGER :: kcr
  REAL(KIND=r8), ALLOCATABLE :: dels  (:)
  REAL(KIND=r8), ALLOCATABLE :: gams  (:)
  REAL(KIND=r8), ALLOCATABLE :: gammod(:)
  REAL(KIND=r8), ALLOCATABLE :: delmod(:)
  REAL(KIND=r8) :: rlocp
  REAL(KIND=r8) :: rgrav
  REAL(KIND=r8) :: rlrv
  REAL(KIND=r8) :: const1
  REAL(KIND=r8) :: const2
  REAL(KIND=r8) :: xx1

  REAL(KIND=r8), PARAMETER :: xkapa=0.2857143_r8


CONTAINS

SUBROUTINE InitShall_Tied(si, del, sl, cl, kmax)
    INTEGER, INTENT(IN) :: kmax
    REAL(KIND=r8),    INTENT(IN) :: si (kmax+1)
    REAL(KIND=r8),    INTENT(IN) :: del(kmax  )
    REAL(KIND=r8),    INTENT(IN) :: sl (kmax  )
    REAL(KIND=r8),    INTENT(IN) :: cl (kmax  )
    CALL InitShalv2(si, del, sl, cl, kmax)

END SUBROUTINE InitShall_Tied


  ! shalv2 :subgrid scale shallow cumulus parameterization - enhanced
  !         vertical temperature and moisture diffusion returning adjusted
  !         temperature and specific humidity.


  
  SUBROUTINE shalv2(si, sl, tin, qin, ps, deltim, &
       ktop, plcl, kuo, kmaxp, kctop1, kcbot1, noshal1, &
       newr, ncols, kmax)
    !
    !**************************************************************************
  !
    !         use:
    !         sr to compute the parameterized effects of shallow
    !         convective clouds (vector version).  routine used in
    !         gwater  moist convection kuolcl.
    !         
    !         purpose:
    !         sub-grid-scale shallow convective cloud parameterization.
    !         this routine computes the effects of shallow convection
    !         based on tiedtke (1984), ecmwf workshop on convection in
    !         large-scale numerical models.
    !         tapered k profile in cloud  developed by caplan and long.
    !         
    !         srs called ... none
    !
    !         
    !         jpg modifications nov 6 85
    ! 
    ! input
    !     
    !     si    : p/ps at sigma level interfaces (array from surface up)
    !     del   : del p positive across sigma lyr  (array from surface up)
    !     sl    : p/ps at sigma layers               (array from surface up)
    !     cl    : 1 - sl                     (array from surface up)
    !     tin   : temperature       (longitude, height array, deg k)
    !     qin   : specific humidity (longitude, height array, gm/gm)
    !     ps    : surface pressure  (longitude array        , cb   )
    !     deltim: timestep (  sec  )
    !     ktop  : cloud tops (sigma layer index) passed from conkuo
    !     plcl  : pressure   (cb) at lcl   passed from conkuo
    !     kuo   : flag to indicate that deep convection was done
    !             kuo, ktop and plcl are longitude arrays
    !     newr  : flag for ccm3 based cloud radiation
    !     output
    !
    !     qin   : updated specific humidity    (gm/gm)
    !     tin   : updated temperature          (deg k)
    !
    !
    !     external variables
    !**************************************************************************
    ! ncols.......Number of grid points on a gaussian latitude circle
    ! kmax.......Number of sigma levels
    ! kmaxp......kmaxp=kmax+1
    ! deltim.....deltim=dt where dt time interval,usually =delt,but changes
    !            in nlnmi (dt=1.) and at dead start(delt/4,delt/2)
    ! si.........si(l)=1.0-ci(l).
    ! ci.........sigma value at each level.
    ! del........sigma spacing for each layer computed in routine "setsig".
    ! ps.........surface pressure
    ! sl.........sigma value at midpoint of
    !                                         each layer : (k=287/1005)
    !
    !                                                                     1
    !                                             +-                   + ---
    !                                             !     k+1         k+1!  k
    !                                             !si(l)   - si(l+1)   !
    !                                     sl(l) = !--------------------!
    !                                             !(k+1) (si(l)-si(l+1)!
    !                                             +-                  -+
    ! kcbot1....
    ! kctop1....
    ! noshal1...
    ! cp........Specific heat of air           (j/kg/k)
    ! hl........heat of evaporation of water     (j/kg)
    ! gasr......gas constant of dry air        (j/kg/k)
    ! grav......grav   gravity constant        (m/s**2)
    !**************************************************************************
    INTEGER, INTENT(in   ) :: ncols
    INTEGER, INTENT(in   ) :: kmax
    INTEGER, INTENT(in   ) :: kmaxp
    REAL(KIND=r8),    INTENT(in   ) :: deltim
    REAL(KIND=r8),    INTENT(in   ) :: si     ( kmaxp )
    REAL(KIND=r8),    INTENT(in   ) :: ps     ( ncols )
    REAL(KIND=r8),    INTENT(in   ) :: sl     ( kmax  )
    REAL(KIND=r8),    INTENT(inout) :: plcl   ( ncols )
    INTEGER, INTENT(in   ) :: ktop   ( ncols )
    INTEGER, INTENT(in   ) :: kuo    ( ncols )
    REAL(KIND=r8),    INTENT(inout) :: tin    ( ncols , kmax )
    REAL(KIND=r8),    INTENT(inout) :: qin    ( ncols , kmax )
    LOGICAL, INTENT(in   ) :: newr
    INTEGER, INTENT(inout) :: kcbot1 (ncols)
    INTEGER, INTENT(inout  ) :: kctop1 (ncols)
    INTEGER, INTENT(inout  ) :: noshal1(ncols)

    REAL(KIND=r8)    :: dk(    ncols,kmax-1)
    REAL(KIND=r8)    :: terp  (ncols,kmax-1)
    REAL(KIND=r8)    :: f     (ncols,kmax  )
    REAL(KIND=r8)    :: ff    (ncols,kmax  )
    REAL(KIND=r8)    :: g     (ncols,kmax  )
    REAL(KIND=r8)    :: gg    (ncols,kmax  )
    REAL(KIND=r8)    :: tnew  (ncols,kmax  )
    REAL(KIND=r8)    :: qnew  (ncols,kmax  )
    REAL(KIND=r8)    :: a     (ncols,kmax  )
    REAL(KIND=r8)    :: b     (ncols,kmax  )
    REAL(KIND=r8)    :: c     (ncols,kmax  )
    REAL(KIND=r8)    :: ud    (ncols,kmax-1)
    REAL(KIND=r8)    :: rec   (ncols)
    INTEGER :: icheck(ncols)
    INTEGER :: kctop (ncols)
    INTEGER :: kcbot (ncols)
    INTEGER :: noshal(ncols)
    LOGICAL :: searching(ncols)

    INTEGER ::  n
    INTEGER ::  i
    INTEGER ::  l
    INTEGER ::  k

    REAL(KIND=r8)          :: dt2
    INTEGER       :: kmaxm
    INTEGER       :: levhm1
    INTEGER       :: loleh1
    INTEGER       :: lolem1
    INTEGER       :: lonlev
    INTEGER       :: lonleh

    dt2 = deltim* 2.0e0_r8
    kmaxm = kmax  -1
    levhm1= kmax -1
    loleh1= ncols *levhm1
    lolem1 =  ncols * kmaxm
    lonlev =  ncols * kmax
    lonleh =  ncols * kmax

    !
    !     zero  arrays
    !
    DO i=1,ncols
       noshal(i)=0
       icheck(i)=0
    END DO

    dk=0.0_r8
    DO i=1,ncols
       IF(kuo(i) .EQ. 1) THEN
          noshal(i)=1
       END IF
    END DO
    !
    !     get cloud base .. overwrite plcl with normalized value
    !
    DO i=1,ncols
       plcl(i)=plcl(i)/ps(i)
    END DO

    !    DO  l=1, ncols
    !       DO n=2, kmax
    !          kcbot(l)=n-1
    !          IF ( plcl(l) .GE. sl(n) ) THEN
    !             EXIT
    !          END IF
    !       END DO
    !    END DO

    searching = .TRUE.
    kcbot = kmax-1
    DO n=2, kmax
       DO  l=1, ncols
          IF (searching(l) .AND. plcl(l) >= sl(n)) THEN
             searching(l) = .FALSE.
             kcbot(l)=n-1
          END IF
       END DO
    END DO
    IF (newr) THEN
       kcbot1(1:ncols) = kcbot(1:ncols)
    END IF
    !
    !     set cloud tops
    !
    DO i=1,ncols
       kctop(i) = MIN(kcr,ktop(i))
       IF(newr)kctop1(i)=kctop(i)
    END DO
    !
    !     check for too high cloud base
    !     remem. that plcl has been divided by ps
    !
    print*, 'Emarg: check too high cloud base',(plcl(1).LT.sl(kbase)),plcl(:),sl(kbase)
    DO i=1,ncols
       IF(plcl(i).LT.sl(kbase)) THEN
          noshal(i) = 1
       END IF
    END DO

    DO i=1,ncols
       IF(kcbot(i).GE.kctop(i)) THEN
          noshal(i) = 1
          kcbot (i) = kctop(i)
       END IF
    END DO
    !
    !     test for moist convective instability
    !
    DO k=1,kmax-1
       DO i=1,ncols
          a(i,k+1) =  tin(i,k+1) - tin(i,k)
          b(i,k+1) =  qin(i,k+1) - qin(i,k)
          c(i,k+1) =  tin(i,k+1) + tin(i,k)
          a(i,k  ) =  &
               (cp/(sl(k+1)-sl(k)))*a(i,k+1)-&
               (0.5e0_r8*gasr/si(k+1))*c(i,k+1)+&
               (hl/(sl(k+1)-sl(k)))*b(i,k+1)
       END DO
    END DO

    print*, 'Emarg: kcbot,kctop,a=', kcbot(:),kctop(:),a(:,:)
    DO k = 1, kmax
       DO i=1, ncols
          IF (k >= kcbot(i) .AND. k <= kctop(i) .AND. a(i,k) > 0.0e0_r8) THEN
             icheck(i)=1
          END IF
       END DO
    END DO

    DO i=1,ncols
       IF(icheck(i).EQ.0) THEN
          noshal(i)=1
       END IF
    END DO
    !
    !     set mixing coefficient dk (m**2/s) profile
    !     dk(n) is value at top of layer n
    !     n.b.  dk(kctop) .ne. 0
    !
    !DO l=1, ncols
    !   IF (noshal(l) .NE. 1) THEN
    !      dk(l,kctop(l)) = 1.0e0_r8
    !      dk(l,kcbot(l)) = 1.5e0_r8
    !      kbetw= kctop(l)-(kcbot(l)+1)
    !      IF(kbetw .GE. 1) dk(l,kctop(l)-1)=3.0e0_r8
    !      IF(kbetw .GT. 1) THEN
    !         DO k=kcbot(l)+1, kctop(l)-2
    !            dk(l,k)=5.0e0_r8
    !         END DO
    !      END IF
    !   END IF
    !END DO
    DO l=1, ncols
       IF (noshal(l) /= 1) THEN
          dk(l,kctop(l)) = 1.0e0_r8
          dk(l,kcbot(l)) = 1.5e0_r8
       END IF
    END DO
    DO l=1, ncols
       IF (noshal(l) /= 1 .AND. kctop(l) >= kcbot(l)+2) THEN
          dk(l,kctop(l)-1)=3.0e0_r8
       END IF
    END DO
    DO k = 1, kmax
       DO l = 1, ncols
          IF (&
               noshal(l) /= 1         .AND. &
               kctop(l) >  kcbot(l)+2 .AND. &
               k >= kcbot(l)+1        .AND. &
               k <= kctop(l)-2               ) THEN
             dk(l,k)=5.0e0_r8
          END IF
       END DO
    END DO
    !
    !     compute adiabatic lapse rate terms for temperature eq.
    !
    DO k=1,kmaxm
       DO i=1,ncols
          terp(i,k) = 1.0e0_r8/(tin(i,k)+tin(i,k+1))
       END DO
    END DO
    !
    !     n.b. terp(i,n) is valid at top of layer n
    !
    DO k=1, kmaxm
       DO i=1,ncols
          ff(i,k) = (dt2*delmod(k))*terp(i,k)*dk(i,k)
          gg(i,k) = (dt2*gammod(k))/(dt2*delmod(k))*ff(i,k)
       END DO

    END DO
    !
    !     compute elements of tridiagonal matrix a,b,c
    !
    a=0.0_r8 ! call reset(a,lonlev)
    b=0.0_r8 ! call reset(b,lonlev)
    c=0.0_r8 ! call reset(c,lonlev)

    DO k=1,kmaxm
       DO i=1,ncols
          terp(i,k) = terp(i,k)*terp(i,k)*dk(i,k)
       END DO
    END DO

    DO n=1, kmaxm
       DO i=1,ncols
          c(i,n) = (-dels(n) *dt2)*terp(i,n)
       END DO
    END DO

    DO n=2, kmax
       DO i=1,ncols
          a(i,n) = (-gams(n-1)*dt2) *terp(i,n-1)
       END DO
    END DO

    DO k=1,kmax
       DO i=1,ncols
          b(i,k)=  1.0e0_r8-a(i,k)-c(i,k)
       END DO
    END DO
    !
    !     compute forcing terms f for temperature and
    !     g for water vapor
    !
    DO i=1,ncols
       f(i,   1)=tin(i,   1)+ff(i,    1)
       f(i,kmax)=tin(i,kmax)-gg(i,kmaxm)
    END DO

    DO n=2, kmaxm
       DO i=1,ncols
          f(i,n)=tin(i,n)+ff(i,n)-gg(i,n-1)
       END DO
    END DO

    DO k=1,kmax
       DO i=1,ncols
          g(i,k)=qin(i,k)
       END DO
    END DO
   !
    !     solution of tridiagonal problems
    !     forward part
    !
    DO i=1,ncols
       rec (i)=1.0e0_r8/b(i,1)
       ud(i,1)=c(i,1)*rec(i)
       ff(i,1)=f(i,1)*rec(i)
       gg(i,1)=g(i,1)*rec(i)
    END DO

    DO k=2, kmaxm
       DO i=1,ncols
          rec(i)=1.0e0_r8/( b(i,k)-a(i,k)*ud(i,k-1))
          ud(i,k)=c(i,k)*rec(i)
          ff(i,k)=rec(i)*(f(i,k)-a(i,k)*ff(i,k-1))
       END DO

       IF(k .LE. levhm1)THEN
          DO i=1,ncols
             gg(i,k)=rec(i)*(g(i,k)-a(i,k)*gg(i,k-1))
          END DO
       END IF

    END DO
    !
    !     now determine solutions
    !
    DO i=1,ncols
       rec(i)=1.0e0_r8/(b(i,kmax)-a(i,kmax)*ud(i,kmaxm))
       tnew(i,kmax)=rec(i)*(f(i,kmax)-a(i,kmax)*ff(i,kmaxm))
    END DO

    DO n= kmaxm ,1,-1
       DO i=1,ncols
          tnew(i,n)=ff(i,n)-tnew(i,n+1)*ud(i,n)
       END DO
    END DO
    !     water vapor solution
    !
    DO i=1,ncols
       rec(i)=1.0e0_r8/(b(i,kmax)-a(i,kmax)*ud(i,levhm1))
       qnew(i,kmax)=rec(i)*(g(i,kmax)-a(i,kmax)*gg(i,levhm1))
       IF(newr)noshal1(i)=noshal(i)
    END DO

    DO k=kmax-1,1,-1
       DO i=1,ncols
          qnew(i,k)=gg(i,k)-ud(i,k)*qnew(i,k+1)
       END DO
    END DO
    !
    !     store solutions
    !
    DO k=1,kmax
       DO i=1,ncols
          IF(noshal(i).EQ.0)THEN
             tin(i,k)=tnew(i,k)
          END IF
       END DO
    END DO

    DO k=1,kmax
       DO i=1,ncols
          IF(noshal(i).EQ.0)THEN
             qin(i,k)=qnew(i,k)
          END IF
       END DO
    END DO
    !
    !     end of computation of shallow convection
    !
!999 FORMAT(' PARAMETER SETTING IS WRONG IN SUBR.SHALV2')
  END SUBROUTINE shalv2

SUBROUTINE InitShalv2(si, del, sl, cl, kmax)
    INTEGER, INTENT(IN) :: kmax
    REAL(KIND=r8),    INTENT(IN) :: si (kmax+1)
    REAL(KIND=r8),    INTENT(IN) :: del(kmax  )
    REAL(KIND=r8),    INTENT(IN) :: sl (kmax  )
    REAL(KIND=r8),    INTENT(IN) :: cl (kmax  )
    
    INTEGER :: n 
    REAL(KIND=r8)    :: gr
    REAL(KIND=r8)    :: ggrr
    REAL(KIND=r8)    :: dryl
  
    ALLOCATE (dels  (kmax-1))
    ALLOCATE (gams  (kmax-1))
    ALLOCATE (gammod(kmax-1))
    ALLOCATE (delmod(kmax-1))

    gr=2.0e0_r8*grav/gasr
    ggrr=gr*gr
    dryl=grav/cp
    DO n = 1, kmax-1
       dels(n)=ggrr*si(n+1)*si(n+1)/ ( &
            del(n)*( cl(n+1) - cl(n) )    )
       delmod(n)=gr*dryl*si(n+1)/del(n)
       gams(n)=del(n)*dels(n)/del(n+1)
       gammod(n)=del(n)*delmod(n)/del(n+1)
       IF(sl(n).GT.0.7e0_r8) THEN
          kbase=n
          kcr=n
       END IF
    END DO
  END SUBROUTINE InitShalv2

END MODULE Shall_Tied
