SUBROUTINE turbl
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INTEGER, PARAMETER :: ktmq2=1
  REAL   , PARAMETER :: capa=0.28589641
  REAL   , PARAMETER :: g=9.8
  REAL   , PARAMETER :: rg=1./g
  REAL   , PARAMETER :: rog=287.04/g
  REAL   , PARAMETER ::  epsz=1.e-4
  REAL   , PARAMETER :: epsq2=0.2
  REAL   , PARAMETER :: d608=0.608
  INTEGER, PARAMETER :: imjm=im*jm-jm/2
  INTEGER, PARAMETER :: lm1=lm-1
  INTEGER, PARAMETER :: lp1=lm+1
  INTEGER, PARAMETER :: jam=6+2*(jm-10)
  INTEGER, PARAMETER :: nsoil=2
  INTEGER, PARAMETER :: itb=76
  INTEGER, PARAMETER :: jtb=134
  INTEGER, PARAMETER :: itbq=152
  INTEGER, PARAMETER :: jtbq=440
  INTEGER, PARAMETER :: la=13
  INTEGER, PARAMETER :: lscrch=4*lm+1+la+1
  INTEGER, PARAMETER :: l1=la+1
  INTEGER, PARAMETER :: l2=la+lm+1
  INTEGER, PARAMETER :: l3=la+2*lm+1
  INTEGER, PARAMETER :: l4=la+3*lm+1
  INCLUDE "ctlblk.f90"
  INCLUDE "loops.f90"
  INCLUDE "masks.f90"
  INCLUDE "dynam.f90"
  INCLUDE "phys2.f90"
  INCLUDE "vrbls.f90"
  INCLUDE "pvrbls.f90"
  INCLUDE "srfdsp.f90"
  INCLUDE "indx.f90"
  INCLUDE "soil.f90"
  INCLUDE "z0efft.f90"
  !
  INTEGER i, j, l
  INTEGER :: lmhk, lmhp, lmhm, lpbl, lmvk
  REAL :: apests, wmsk, rwmsk, hpbl, ulm, vlm, wstar
  REAL pdsl(im,jm),cklq(im,jm),ct(im,jm)
  REAL ape (im,jm,lm),akh (im,jm,lm1)
  REAL zint(0:im+1,0:jm+1,lp1),akm (0:im+1,0:jm+1,lm1)
  REAL uz0h(0:im+1,0:jm+1)    ,vz0h(0:im+1,0:jm+1)
  REAL tcol (lm),qcol (lm),q2col(lm)
  REAL ucol (lm),vcol (lm)
  REAL gm  (lm1),gh  (lm1),el  (lm1)
  REAL zcol(lp1)
  REAL akhcol(lm1),akmcol(lm1),apecol(lm)
  REAL zeff  (  4)
  INTEGER lpblwork(im,jm)
  REAL    hpblwork(im,jm)
  REAL, DIMENSION(im,jm,lm1) :: gmwork, ghwork, elwork

  !***  areas alocadas para vetorizacao

  REAL, DIMENSION(im, jm, lm1) :: akhwork, akmwork
  REAL, DIMENSION(im, jm, lp1) :: zcolwork
  REAL, DIMENSION(im, jm     ) :: akmsvwork
  REAL, DIMENSION(im, jm, lp1) :: uwork, vwork

  CALL zero3(zint,lp1)
  CALL zero3(akm,lm1)
  CALL zero2(uz0h)
  CALL zero2(vz0h)

  IF(ntsd.EQ.1)THEN
     DO j=1,jm
        DO i=1,im
           ustar(i,j)=0.1
           thz0(i,j)=ths(i,j)
           qz0(i,j)=qs(i,j)
           uz0(i,j)=0.
           vz0(i,j)=0.
        ENDDO
     ENDDO
  ENDIF

  DO l=1,lm1
     DO j=1,jm
        DO i=1,im
           akh(i,j,l)=0.
           akm(i,j,l)=0.
        ENDDO
     ENDDO
  ENDDO

  !***  compute the heights of the layer interfaces and the exner function

  DO  j=1,jm
     DO  i=1,im
        pdsl(i,j)=pd(i,j)*res(i,j)
        zint(i,j,lp1)=epsz 
        IF(sigma)zint(i,j,lp1)=rg*fis(i,j)

        DO  l=lm,1,-1
           apests=pdsl(i,j)*aeta(l)+pt
           zint(i,j,l)=t(i,j,l)/apests*pdsl(i,j)*deta(l)*rog &
                *(q(i,j,l)*d608+1.)+zint(i,j,l+1)
           zint(i,j,l)=(zint(i,j,l)-dfrlg(l))*htm(i,j,l)+dfrlg(l)
        END DO

        DO  l=lm,1,-1
           apests=pdsl(i,j)*aeta(l)+pt
           ape(i,j,l)=(1.e5/apests)**capa
        END DO
     END DO
  END DO

  !***  remove negative q2

  DO l=1,lm
     DO j=1,jm
        DO i=1,im
           q2(i,j,l)=amax1(q2(i,j,l)*hbm2(i,j),epsq2)
        END DO
     END DO
  END DO

  DO j=3,jm-2
     DO i=1,im
        uz0h(i,j)=(uz0(i+ihe(j),j)+uz0(i+ihw(j),j) &
             +uz0(i,j+1)+uz0(i,j-1))*hbm2(i,j)*0.25
        vz0h(i,j)=(vz0(i+ihe(j),j)+vz0(i+ihw(j),j) &
             +vz0(i,j+1)+vz0(i,j-1))*hbm2(i,j)*0.25
     END DO
  END DO

  !***                      the 1st grand loop
  !***               prepare the exchange coefficients

  DO l=1,lm
     DO j=3,jm-2
        DO i=1,im
           IF (l <= lmh(i,j)) THEN
              wmsk=vtm(i+ihe(j),j,l)+vtm(i+ihw(j),j,l) &
                   +vtm(i,j+1,l)+vtm(i,j-1,l)
              IF(wmsk.GT.0.)THEN
                 rwmsk=1./wmsk
                 uwork(i,j,l)=(u(i+ihe(j),j,l)*vtm(i+ihe(j),j,l) &
                      +u(i+ihw(j),j,l)*vtm(i+ihw(j),j,l) &
                      +u(i,j+1,l)*vtm(i,j+1,l)+u(i,j-1,l)*vtm(i,j-1,l))*rwmsk
                 vwork(i,j,l)=(v(i+ihe(j),j,l)*vtm(i+ihe(j),j,l) &
                      +v(i+ihw(j),j,l)*vtm(i+ihw(j),j,l) &
                      +v(i,j+1,l)*vtm(i,j+1,l)+v(i,j-1,l)*vtm(i,j-1,l))*rwmsk
              ELSE
                 uwork(i,j,l)=0.
                 vwork(i,j,l)=0.
              END IF
           ENDIF
        END DO
     END DO
  END DO

  !***  find the mixing length

  CALL mixlen(lmh, lpblwork, hpblwork, uwork, vwork,&
       t, q, q2, ape, zint, gmwork, ghwork, elwork)

  !***  solve for the production/dissipation of 
  !***  the turbulent kinetic energy

  CALL prodq2(lmh, dtq2, ustar, gmwork, ghwork, elwork, q2)

  !***  find the exchange coefficients in the free atmosphere

  CALL difcof(lmh, gmwork, ghwork, elwork, q2, zint, akm, akh)

  !***  carry out the vertical diffusion of
  !***  turbulent kinetic energy

  CALL vdifq  (lmh, ktmq2, dtq2, q2, elwork, zint)

  !***  find the surface exchange coefficients

  CALL sfcdif (lmh, sm, ths, qs, uz0h, vz0h, thz0, qz0, &
       ustar, z0, zeffij, akms, akhs, hpblwork, ct, &
       u10, v10, tshltr, th10, qshltr, q10, &
       uwork, vwork, t, q, ape, zint)

  !***                      end of 1st grand loop
  !***  average uz0 and vz0 back to v points

  DO  j=3,jm-2
     DO  i=1,im
        uz0(i,j)=(uz0h(i+ive(j),j)*hbm2(i+ive(j),j) &
             +uz0h(i+ivw(j),j)*hbm2(i+ivw(j),j) &
             +uz0h(i,j+1)*hbm2(i,j+1)+uz0h(i,j-1)*hbm2(i,j-1))*0.25
        vz0(i,j)=(vz0h(i+ive(j),j)*hbm2(i+ive(j),j) &
             +vz0h(i+ivw(j),j)*hbm2(i+ivw(j),j) &
             +vz0h(i,j+1)*hbm2(i,j+1)+vz0h(i,j-1)*hbm2(i,j-1))*0.25
     END DO
  END DO

  !***  execute the ground processes

  CALL surfce(ape,zint,cklq)

  !***                      the 2nd grand loop
  !***                 execute the vertical exchange

  DO j=3,jm-2
     DO i=1,im
        thz0(i,j)=(1.-sm(i,j))*ths(i,j)+sm(i,j)*thz0(i,j)
        qz0 (i,j)=(1.-sm(i,j))*qs (i,j)+sm(i,j)*qz0 (i,j)
     END DO
  END DO

  DO l=1,lm1
     DO j=3,jm-2
        DO i=1,im
           akhwork(i,j,l)=akh(i,j,l)*htm(i,j,l)*hbm2(i,j)
        END DO
     END DO
  END DO

  !***  carry out the vertical diffusion of
  !***  temperature and water vapor

  CALL vdifh(lmh, ktmq2, dtq2, thz0, qz0, akhs, ct, cklq, &
       t, q, akhwork, ape, zint)

  !***  carry out the vertical diffusion of
  !***  velocity components

  DO l=1,lm1
     DO j=3,jm-2
        DO i=1,im
           akmwork(i,j,l)=(akm(i+ive(j),j,l)*hbm2(i+ive(j),j) &
                +akm(i+ivw(j),j,l)*hbm2(i+ivw(j),j) &
                +akm(i,j+1,l)*hbm2(i,j+1)+akm(i,j-1,l)*hbm2(i,j-1)) &
                *vtm(i,j,l)*vbm2(i,j)*0.25
        END DO
     END DO
  END DO

  DO j=3,jm-2
     DO i=1,im
        akmsvwork(i,j)=(akms(i+ive(j),j)*hbm2(i+ive(j),j) &
             +akms(i+ivw(j),j)*hbm2(i+ivw(j),j) &
             +akms(i,j+1)*hbm2(i,j+1)+akms(i,j-1)*hbm2(i,j-1)) &
             *vbm2(i,j)*0.25
     END DO
  END DO

  DO l=1,lp1
     DO j=3,jm-2
        DO i=1,im
           zcolwork(i,j,l)=0.25*(zint(i+ive(j),j,l)+zint(i+ivw(j),j,l) &
                +zint(i,j+1,l)+zint(i,j-1,l))
        END DO
     END DO
  END DO

  CALL vdifv (lmv, ktmq2, dtq2,uz0, vz0, &
       akmsvwork, u, v, akmwork, zcolwork)

  !*** end of 2nd grand loop

  !***  remove negative q2

  DO l=1,lm
     DO j=1,jm
        DO i=1,im
           q2(i,j,l)=amax1(q2(i,j,l)*hbm2(i,j),epsq2)
        END DO
     END DO
  END DO
END SUBROUTINE turbl
