subroutine atm_diffus(lbase_min,lbase_max,dtime) 2,23
!@sum  atm_diffus updates u,v,t,q due to turbulent transport throughout
!@+    all GCM layers using a non-local turbulence model
!@auth Ye Cheng/G. Hartke (modifications by G. Schmidt)
!@ver  1.0 (from diffB347D6M20)
!@cont atm_diffus,getdz,dout,de_solver_main,de_solver_edge,l_gcm,k_gcm,
!@+    e_gcm,find_pbl_top,zze,ave_uv_to_agrid,
!@+    ave_uv_to_bgrid,ave_s_to_bgrid,apply_fluxes_to_atm
!@var lbase_min/max levels through which to apply turbulence (dummy)
!@var dtime time step

!@var qmin minimum value of specific humidity
!@var itest/jtest longitude/latitude at which dout may be called
!@var call_diag logical variable whether dout is called

      USE CONSTANT, only : grav,deltx,lhe,sha,by3,teeny,mb2kg
      USE MODEL_COM, only :
     *     im,jm,lm,sig,sige,u_3d=>u,v_3d=>v,t_3d=>t,q_3d=>q,itime,psf
     *     ,pmtop
cc      USE QUSDEF, only : nmom,zmoms,xymoms
cc      USE SOMTQ_COM, only : tmom,qmom
      USE GEOM, only : imaxj,kmaxj,ravj,idij,idjj,bydxyp,dxyp
      USE DYNAMICS, only : pk,pdsig,plij,pek,byam,am,dke
      USE DAGCOM, only : ajl,jl_trbhr,jl_damdc,jl_trbke,jl_trbdlht

      USE SOCPBL, only : b1,b123,prt,kappa,zgs
      USE PBLCOM, only : tsavg,qsavg,dclev,uflux,vflux,tflux,qflux
     *     ,e_3d=>egcm,w2_3d=>w2gcm !,t2_3d=>t2gcm
      USE FLUXES, only : uflux1,vflux1,tflux1,qflux1


      IMPLICIT NONE

      integer, intent(in) :: lbase_min,lbase_max
      real*8, intent(in) :: dtime

      real*8, parameter :: qmin=1.d-20
      integer, parameter :: itest= 1,jtest=11
      logical, parameter :: call_diag=.false.

      real*8, dimension(lm) :: u,v,t,q,e,u0,v0,t0,q0,e0
     &    ,dudz,dvdz,dtdz,dqdz,g_alpha,as2,an2
     &    ,rhoebydz,bydzerho,rho,rhoe,dz,dze
     &    ,km,kh,ke,wt_nl,wq_nl
     &    ,lscale,qturb,p3,p4,rhobydze,bydzrhoe,uw,vw,wt,wq
      real*8, dimension(lm+1) :: ze

      real*8, dimension(lm,im,jm) :: u_3d_old,rho_3d,rhoe_3d,dz_3d
     &    ,dze_3d,u_3d_agrid,v_3d_agrid,t_3d_virtual,km_3d,km_3d_bgrid
     &    ,dz_3d_bgrid,dze_3d_bgrid,rho_3d_bgrid,rhoe_3d_bgrid
     &    ,wt_3d,v_3d_old
      real*8, dimension(im,jm) :: tvsurf,uflux_bgrid,vflux_bgrid
cc      real*8, dimension(nmom,lm) :: tmomij,qmomij

      real*8 :: uflx,vflx,tvflx,qflx,tvs
     &   ,ustar2,t0ijl,tijl,rak,alpha1,ustar
     &   ,flux_bot,flux_top,x_surf
     &   ,wstar,dbl,lmonin,tpe0,tpe1,ediff
      integer :: idik,idjk,ldbl,kmax,
     &    i,j,l,k,n,iter !@i,j,l,k,n,iter loop variable


      ! Note that lbase_min/max are here for backwards compatibility
      ! with original drycnv. They are only used to determine where the
      ! routine has been called from.

      if (lbase_min.eq.2) return       ! quit if called from main

      !  convert input T to virtual T

!$OMP  PARALLEL DO PRIVATE (L,I,J)
      do j=1,jm
        do i=1,imaxj(j)
          !@var tvsurf(i,j) surface virtual temperature
          !@var tsavg(i,j) composite surface air temperature (k)
          tvsurf(i,j)=tsavg(i,j)*(1.d0+deltx*qsavg(i,j))
          do l=1,lm
            ! t_3d_virtual is virtual potential temp. referenced at 1 mb
            t_3d_virtual(l,i,j)=t_3d(i,j,l)*(1.d0+deltx*q_3d(i,j,l))
          end do
        end do
      end do
!$OMP  END PARALLEL DO



      ! integrate equations other than u,v at agrids

      ! get u_3d_agrid and v_3d_agrid
      call ave_uv_to_agrid(u_3d,v_3d,u_3d_agrid,v_3d_agrid,im,jm,lm)

      call getdz(t_3d_virtual,dz_3d,dze_3d,rho_3d,rhoe_3d,tvsurf
     &    ,im,jm,lm)

!$OMP  PARALLEL DO PRIVATE (L,I,J,u,v,t,q,e,rho,rhoe,t0,q0,e0,qturb,
!$OMP*   dze,dz,bydzerho,rhobydze,bydzrhoe,rhoebydz,tvs,uflx,vflx,
!$OMP*   qflx,tvflx,ustar,ustar2,alpha1,dudz,dvdz,dtdz,dqdz,g_alpha,
!$OMP*   an2,as2,ze,lscale,dbl,ldbl,wstar,kh,km,ke,wt,wq,uw,vw,
!$OMP*   wt_nl,wq_nl,lmonin,p3,p4,x_surf,flux_bot,flux_top,t0ijl,tijl,
!$OMP*   tpe0,tpe1,ediff

!$OMP*    ) SHARED(dtime

!$OMP*    ) SCHEDULE(DYNAMIC,2)

      loop_j_tq: do j=1,jm
        loop_i_tq: do i=1,imaxj(j)

          do l=1,lm
            u(l)=u_3d_agrid(l,i,j)
            v(l)=v_3d_agrid(l,i,j)
            ! virtual potential temp. referenced at 1 mb
            t(l)=t_3d_virtual(l,i,j)
            q(l)=q_3d(i,j,l)
cc            qmomij(:,l)=qmom(:,i,j,l)
cc            tmomij(:,l)=tmom(:,i,j,l) ! vert. grad. should virtual ?
            if(            e(l)=e_3d(l,i,j) !e_3d was called egcM
            ! t2(l)=t2_3d(l,i,j)  ! not in use
            rho(l)=rho_3d(l,i,j)
            rhoe(l)=rhoe_3d(l,i,j)
            t0(l)=t(l)
            q0(l)=q(l)
            e0(l)=e(l)
            qturb(l)=sqrt(2.d0*e(l))
            dze(l)=dze_3d(l,i,j)
            dz(l)=dz_3d(l,i,j)
            bydzerho(l)=1.d0/(dze(l)*rho(l))
            rhobydze(l)=rho(l)/dze(l)

          end do
          bydzrhoe(1)=0.
          rhoebydz(1)=0.
          do l=1,lm-1
            bydzrhoe(l+1)=1.d0/(dz(l)*rhoe(l+1))
            rhoebydz(l+1)=rhoe(l+1)/dz(l)
          end do

          ! tvs is surface virtual potential temp. referenced at 1 mb
          tvs=tvsurf(i,j)/pek(1,i,j)
          uflx=uflux1(i,j)/rhoe(1)
          vflx=vflux1(i,j)/rhoe(1)
          qflx=qflux1(i,j)/rhoe(1)
          ! tvflx is virtual, potential temp. flux referenced at 1 mb
          tvflx=tflux1(i,j)*(1.d0+deltx*qsavg(i,j))/(rhoe(1)*pek(1,i,j))
     &         +deltx*tsavg(i,j)/pek(1,i,j)*qflx
          ! redefine uflux1,vflux1 for later use
          uflux1(i,j)=uflx
          vflux1(i,j)=vflx


          if(abs(uflx).lt.teeny) uflx=sign(teeny,uflx)
          if(abs(vflx).lt.teeny) vflx=sign(teeny,vflx)
          ustar=(uflx*uflx+vflx*vflx)**(0.25d0)
          ustar2=ustar*ustar
          alpha1=atan2(vflx,uflx)

          ! calculate z-derivatives at the surface

          ! @var zgs height of surface layer (m), imported from SOCPBL
          dudz(1)=ustar/(kappa*zgs)*cos(alpha1)
          dvdz(1)=ustar/(kappa*zgs)*sin(alpha1)
          dtdz(1)=(tvflx*prt/ustar)/(kappa*zgs)
          dqdz(1)=(qflx*prt/ustar)/(kappa*zgs)

          g_alpha(1)=grav/tvs

          ! calculate z-derivatives on the edges of the layers

          do l=2,lm
            dudz(l)=(            dvdz(l)=(            dtdz(l)=(            dqdz(l)=(            g_alpha(l)=grav*2.d0/(          end do

          !@var an2 brunt-vassala frequency
          !@var as2 shear number squared

          do l=1,lm
            an2(l)=g_alpha(l)*dtdz(l)
            as2(l)=dudz(l)*dudz(l)+dvdz(l)*dvdz(l)
          end do

          ! calculate turbulence length scale lscale
          call zze(dze,ze,lm)
          call l_gcm(ze,lscale,lm)

          call find_pbl_top(e,ze,dbl,ldbl,lm)
          if(tvflx.lt.0.) then ! convective case
            wstar=(-g_alpha(1)*tvflx*dbl)**by3
          else
            wstar=0.
          endif
          ! calculate turbulent diffusivities km,kh and ke
          call k_gcm(tvflx,qflx,ustar,wstar,dbl
     &        ,ze,lscale,e,qturb,an2,as2,dtdz,dqdz,dudz,dvdz
     &        ,kh,km,ke,wt,wq,uw,vw,wt_nl,wq_nl

     &        ,lm)

          call e_gcm(tvflx,wstar,ustar,dbl,lmonin,ze,g_alpha
     &              ,an2,as2,lscale,e,lm)

          ! integrate differential eqn for e
          p3(1)=0. ; p3(lm)=0.
          p4(1)=0. ; p4(lm)=0.
c         do l=2,lm-1
c             p3(l)=2.d0*(qturb(l)/(b1*lscale(l)))
c             p4(l)=-uw(l)*dudz(l)-vw(l)*dvdz(l)+g_alpha(l)*wt(l)
c             ! p4(l)=km(l)*as2(l)-kh(l)*an2(l)
c         end do
c         x_surf=0.5d0*b123*ustar2
c         call de_solver_edge(e,e0,ke,p3,p4,
c    &        rhobydze,bydzrhoe,x_surf,dtime,lm,.true.)

          do l=1,lm
              qturb(l)=sqrt(2.d0*e(l))
          end do

          ! integrate differential eqn for T
          do l=2,lm-1
              p4(l)=-(rhoe(l+1)*wt_nl(l+1)-rhoe(l)*wt_nl(l))
     &              *bydzerho(l)
          end do
          flux_bot=rhoe(1)*tvflx+rhoe(2)*wt_nl(2)
          flux_top=0.
          call de_solver_main(t,t0,kh,p4,
     &        rhoebydz,bydzerho,flux_bot,flux_top,dtime,lm,.false.)

C**** also diffuse moments
cc        call diff_mom(tmomij)

          ! integrate differential eqn for Q
          do l=2,lm-1
              p4(l)=-(rhoe(l+1)*wq_nl(l+1)-rhoe(l)*wq_nl(l))
     &              *bydzerho(l)
C**** check on physicality of non-local fluxes....
              if ( p4(l)*dtime+q0(l).lt.0 ) then
                p4(l)=-q(l)/dtime
                wq_nl(l+1)=(q0(l)/(dtime*bydzerho(l))+rhoe(l)
     *               *wq_nl(l))/rhoe(l+1)
              end if
          end do
          flux_bot=rhoe(1)*qflx+rhoe(2)*wq_nl(2)
C**** fix first layer for rare tracer problems
C**** Does this ever happen for q? (put this in just in case)
            if ( q0(1)-dtime*bydzerho(1)*flux_bot.lt.0 ) then
              flux_bot=-q0(1)/(dtime*bydzerho(1))
              wq_nl(2)=(flux_bot-rhoe(1)*qflx)/rhoe(2)
            end if
          flux_top=0.

          call de_solver_main(q,q0,kh,p4,
     &        rhoebydz,bydzerho,flux_bot,flux_top,dtime,lm,.true.)
          do l=1,lm
              if(          end do

C**** also diffuse moments
cc        call diff_mom(qmomij)


          call find_pbl_top(e,ze,dbl,ldbl,lm)
          dclev(i,j)=real(ldbl)

C**** calculate possible energy loss
          tpe0=-tflux1(i,j)*dtime*sha
          tpe1=0.
          do l=1,lm
            tpe0=tpe0+t_3d(i,j,l)*pk(l,i,j)*pdsig(l,i,j)*sha*mb2kg    
            tpe1=tpe1+t(l)*pk(l,i,j)*pdsig(l,i,j)*sha*mb2kg/(1.d0+deltx
     *           *q(l))
          end do
          ediff=(tpe1-tpe0)/((psf-pmtop)*sha*mb2kg)        ! C

          do l=1,lm
C**** correct virt pot t for energy error
            t(l) = t(l) - ediff*(1.d0+deltx*q(l))/pk(l,i,j)
            ! update 3-d t,q,e and km
            t0ijl=t_3d(i,j,l)
            tijl=t(l)/(1.d0+deltx*q(l))
            t_3d(i,j,l)=tijl
C**** moment variation to be added
cc            qmom(:,i,j,l)=qmomij(:,l)
cc            tmom(:,i,j,l)=tmomij(:,l)

            q_3d(i,j,l)=q(l)
            e_3d(l,i,j)=e(l)
            w2_3d(l,i,j)=2.*by3*e(l)
            ! t2_3d(l,i,j)=t2(l)  ! not in use
            km_3d(l,i,j)=km(l)
            ! ACCUMULATE DIAGNOSTICS for t and q
            AJL(J,L,JL_TRBHR)=AJL(J,L,JL_TRBHR)
     2                 +(tijl-t0ijl)*PK(L,I,J)*PLIJ(L,I,J)
            AJL(J,L,JL_TRBDLHT)=AJL(J,L,JL_TRBDLHT)
     2                 +(            AJL(J,L,JL_TRBKE)=AJL(J,L,JL_TRBKE)+e(l)

          end do

          ! Write out diagnostics if at selected grid point:

          if (call_diag.and.(i.eq.itest).and.(j.eq.jtest)) then
            call dout(ze,dz,u,v,t,q,ke,dtdz,dqdz,an2,as2
     &         ,wt_nl,wq_nl,kh,km,e,lscale
     &         ,uflx,vflx,tvflx,qflx,dbl,ldbl,i,j,lm)
          endif

        end do loop_i_tq
      end do loop_j_tq
!$OMP  END PARALLEL DO

      ! integrate U,V equations at bgrids

      call ave_uv_to_bgrid(uflux1,vflux1,uflux_bgrid,vflux_bgrid,
     &                     im,jm,1)
      call ave_s_to_bgrid(km_3d,dz_3d,dze_3d,rho_3d,rhoe_3d,
     &  km_3d_bgrid,dz_3d_bgrid,dze_3d_bgrid,rho_3d_bgrid,
     &  rhoe_3d_bgrid,im,jm,lm)

!$OMP  PARALLEL DO PRIVATE (L,I,J,u,v,rho,rhoe,u0,v0,km,dze,dz,
!$OMP*   bydzerho,rhoebydz,uflx,vflx,p4,flux_bot,flux_top)
!$OMP*    SHARED(dtime)
!$OMP*    SCHEDULE(DYNAMIC,2)
      loop_j_uv: do j=2,jm
        loop_i_uv: do i=1,im

          do l=1,lm
            u(l)=u_3d(i,j,l)
            v(l)=v_3d(i,j,l)
            rho(l)=rho_3d_bgrid(l,i,j)
            rhoe(l)=rhoe_3d_bgrid(l,i,j)
            u0(l)=u(l)
            v0(l)=v(l)
            u_3d_old(l,i,j)=u(l)
            v_3d_old(l,i,j)=v(l)
            km(l)=km_3d_bgrid(l,i,j)
            dze(l)=dze_3d_bgrid(l,i,j)
            dz(l)=dz_3d_bgrid(l,i,j)
            bydzerho(l)=1.d0/(dze(l)*rho(l))
          end do
          rhoebydz(1)=0.
          do l=1,lm-1
            rhoebydz(l+1)=rhoe(l+1)/dz(l)
          end do

          uflx=uflux_bgrid(i,j)
          vflx=vflux_bgrid(i,j)

          ! integrate differential eqns for U and V
          do l=2,lm-1
              p4(l)=0.d0
          end do
          flux_bot=rhoe(1)*uflx
          flux_top=0.d0
          call de_solver_main(u,u0,km,p4,
     &        rhoebydz,bydzerho,flux_bot,flux_top,dtime,lm,.false.)
          flux_bot=rhoe(1)*vflx
          call de_solver_main(v,v0,km,p4,
     &        rhoebydz,bydzerho,flux_bot,flux_top,dtime,lm,.false.)

          ! update u and v
          do l=1,lm
            u_3d(i,j,l)= u(l)
            v_3d(i,j,l)= v(l)
          end do

        end do loop_i_uv
      end do loop_j_uv
!$OMP  END PARALLEL DO

      ! ACCUMULATE DIAGNOSTICS for u and v

      DO J=1,JM
        KMAX=KMAXJ(J)
        DO I=1,IMAXJ(J)
          DO L=1,LM
            DO K=1,KMAX
              RAK=RAVJ(K,J)
              IDIK=IDIJ(K,I,J)
              IDJK=IDJJ(K,J)
              AJL(IDJK,L,JL_DAMDC)=AJL(IDJK,L,JL_DAMDC)
     &        +(U_3d(IDIK,IDJK,L)-u_3d_old(L,IDIK,IDJK))*PLIJ(L,I,J)*RAK
            ENDDO
          ENDDO
        ENDDO
      ENDDO

C**** Save additional changes in KE for addition as heat later
!$OMP  PARALLEL DO PRIVATE (L,I,J)
      DO L=1,LM
      DO J=2,JM
      DO I=1,IM
        DKE(I,J,L)=DKE(I,J,L)+0.5*(U_3d(I,J,L)*U_3d(I,J,L)
     *       +V_3d(I,J,L)*V_3d(I,J,L)-U_3d_old(L,I,J)*U_3d_old(L,I,J)
     *       -V_3d_old(L,I,J)*V_3d_old(L,I,J))
      END DO
      END DO
      END DO
!$OMP  END PARALLEL DO

      return
      end subroutine atm_diffus


      subroutine getdz(tv,dz,dze,rho,rhoe,tvsurf,im,jm,lm) 1,5
!@sum  getdz computes the 3-d finite difference dz and dze
!@+    as well as the 3-d density rho and rhoe
!@+    called at the primary grid (A-grid)
!@auth Ye Cheng/G. Hartke
!@ver  1.0
!@var  tv virtual potential temp. referenced at 1 mb
!@var  dz main grid spacing
!@var  dze edge grid spacing
!@var  rho,rhoe air density at the main/edge grids
!@var  tvsurf surface virtual temperature
!@var  im,jm,lm 3-d grids

      !
      !     Grids:
      !
      !                   -------------------------  lm+1
      !                lm  - - - - - - - - - - - - -
      !                   -------------------------  lm
      !                l+1 - - - - - - - - - - - - -
      !                    -------------------------  l+1
      !     (main)     l   - - - - - - - - - - - - -
      !                    -------------------------  l     (edge)
      !                l-1 - - - - - - - - - - - - -
      !                    -------------------------  l-1
      !                2   - - - - - - - - - - - - -
      !                    -------------------------    2
      !                1   - - - - - - - - - - - - -
      !                    -------------------------    1
      !           dz(l,i,j) z(l+1,i,j) - z(l,i,j)
      !           dze(l,i,j) ze(l+1,i,j) - ze(l,i,j)
      !           rhoe(l+1,i,j)=100.d0*(pl-pl1)/(grav*dz(l,i,j))
      !           rho(l,i,j)=100.d0*(ple-pl1e)/(grav*dze(l,i,j))
      !
      !     at main: u,v,tv,q,ke
      !     at edge: e,lscale,km,kh,gm,gh
      !
      USE CONSTANT, only : grav,rgas
      USE GEOM, only : imaxj
      USE DYNAMICS, only : pmid,pk,pedn

      implicit none

      integer, intent(in) :: im,jm,lm
      real*8, dimension(im,jm), intent(in) :: tvsurf
      real*8, dimension(lm,im,jm), intent(in) :: tv
      real*8, dimension(lm,im,jm), intent(out) :: rho,rhoe,dz,dze

      real*8 :: temp0,temp1,temp1e,pl1,pl,pl1e,ple,plm1e
      integer :: i,j,l  !@var i,j,l loop variable

      !@ temp0 virtual temperature (K) at (i,j) and SIG(l)
      !@ temp1 virtual temperature (K) at (i,j) and SIG(l+1)
      !@ temp1e average of temp0 and temp1
!$OMP  PARALLEL DO PRIVATE (J,I,L,pl1,pl,pl1e,ple,temp0,temp1,temp1e,
!$OMP*  plm1e)
!$OMP*    SCHEDULE(DYNAMIC,2)
      do j=1,jm
        do i=1,imaxj(j)
          do l=1,lm-1

            pl1 =pmid(l+1,i,j)
            pl  =pmid(l,i,j)
            pl1e=pedn(l+1,i,j)
            ple =pedn(l,i,j)
            temp0 =tv(l,i,j)*pk(l,i,j)
            temp1 =tv(l+1,i,j)*pk(l+1,i,j)
            temp1e=0.5d0*(temp0+temp1)
            dz(l,i,j)    =-(rgas/grav)*temp1e*log(pl1/pl)
            dze(l,i,j)=-(rgas/grav)*temp0 *log(pl1e/ple)
            rhoe(l+1,i,j)=100.d0*(pl-pl1)/(grav*dz(l,i,j))
            rho(l,i,j)=100.d0*(ple-pl1e)/(grav*dze(l,i,j))
            if(l.eq.1) then
              rhoe(1,i,j)=100.d0*ple/(tvsurf(i,j)*rgas)
            endif
            if(l.eq.lm-1) then
              plm1e=pedn(lm+1,i,j)
              dze(lm,i,j)=-(rgas/grav)*temp1 *log(plm1e/pl1e)
              dz(lm,i,j)=0.
              rho(lm,i,j)=100.d0*(pl1e-plm1e)/(grav*dze(lm,i,j))
            endif

          end do
        end do
      end do
!$OMP  END PARALLEL DO

      return
      end subroutine getdz


      subroutine dout(ze,dz,u,v,t,q,ke,dtdz,dqdz,an2,as2 1,1
     &      ,wt_nl,wq_nl,kh,km,e,lscale
     &      ,uflx,vflx,tvflx,qflx,dbl,ldbl,i,j,n)
!@sum dout writes out diagnostics at (i,j)
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
!@var p  pressure at main grid z
!@var pe  pressure at secondary grid ze
!@var u  west-east   velocity component
!@var v  south-north velocity component
!@var t  virt. pot. temperature (referenced to 1mb)
!@var q  specific humidity
!@var e  turbulent kinetic energy
!@var ze hight at edge lyer
!@var dz(j) z(j+1)-z(j)
!@var dze(j) ze(j+1)-ze(j)
!@var as2 dudz^2+dvdz^2
!@var dtdz z-derivative of t at edge grid ze
!@var dqdz z-derivative of q at edge grid ze
!@var an2 g_alpha*dtdz, brunt-vasala frequency
!@var km turbulent viscosity for u and v equations
!@var kh turbulent diffusivity for t,q
!@var ke turbulent diffusivity for e
!@var lscale  turbulent length scale
!@var uflx momentun flux -uw at surface, ze(1)
!@var vflx momentun flux -vw at surface, ze(1)
!@var tvflx heat flux -wt at surface, ze(1)
!@var qflx moisture flux -wq at surface, ze(1)
!@var i/j horizontal location at which the output is written
!@var n number of vertical main layers

      USE DYNAMICS, only : pmid,pedn,pk,pek

      implicit none

      integer, intent(in) :: ldbl,i,j,n 
      real*8, dimension(n), intent(in) ::
     &   dz,u,v,t,q,ke,dtdz,dqdz,an2,as2
     &   ,wt_nl,wq_nl,kh,km,e,lscale
      real*8, dimension(n+1), intent(in) :: ze
      real*8, intent(in) :: uflx,vflx,tvflx,qflx,dbl

      real*8, dimension(n) :: p,pe
      real*8 :: z,ri,wt_lcl,wq_lcl
      integer :: l  !@var l loop variable

      Write (67,1100) "i=",i,"j=",j
      Write (67,1200) "uflx=",uflx,"vflx=",vflx
      Write (67,1200) "tvflx=",tvflx*pek(1,i,j),"qflx=",qflx
      Write (67,1250) "ldbl=",ldbl,"dbl=",dbl

      ! pressure at main and edge layers

      do l=1,n
          p(l)=100.*pmid(l,i,j)
          pe(l)=100.*pedn(l,i,j)
      end do

      ! fields at layer center

      write (67,1300)
      do l=1,n
        z=.5d0*(ze(l)+ze(l+1))
        write (67,2000) l,p(l),z,dz(l),u(l),v(l),t(l)*pk(l,i,j),q(l)
     &                  ,ke(l)
      end do
      write (67,*)

      ! fields at layer edge

      write (67,1400)
      do l=1,n
        wt_lcl=-kh(l)*dtdz(l)*pek(l,i,j)
        wq_lcl=-kh(l)*dqdz(l)
        ri=an2(l)/as2(l)
        write (67,2000) l,pe(l),ze(l),wt_lcl,wt_nl(l)*pek(l,i,j)
     &    ,wq_lcl,wq_nl(l),kh(l),km(l),lscale(l),ri,e(l)
      end do
      write (67,1500)

      return
1100  format(1200  format(1250  format(2x,a,i6,2x,a,1pe14.4)
1300  format (' ',' l',1x,
     1            '     p     ',1x,
     2            '     z     ',1x,'     dz    ',1x,'     u     ',1x,
     3            '     v     ',1x,'     t     ',1x,'     q     ',1x,
     4            '     ke    ',1x,'           ',1x,'           ',1x,
     5            '           ',/)
        write (67,2000) l,pe(l),ze(l),wt_lcl,wt_nl(l),wq_lcl,wq_nl(l)
     2    ,kh(l),km(l),lscale(l),ri,e(l)
1400  format (' ',' l',1x,
     2            '  p (edge) ',1x,
     2            '  z (edge) ',1x,'  wt_lcl   ',1x,'   wt_nl   ',1x,
     3            '   wq_lcl  ',1x,'  wq_nl    ',1x,'     kh    ',1x,
     4            '     km    ',1x,'  lscale   ',1x,'     ri    ',1x,
     5            '     e     ',/)
1500  format(132('-'))
2000  format (1x,i2,1x,1pe11.4,9(1x,1pe11.4),1x,1pe10.3)
      end subroutine dout


      subroutine de_solver_main(x,x0,p1,p4, 4,1
     &    rhoebydz,bydzerho,flux_bot,flux_top,dtime,n,qlimit)
!@sum differential eqn solver for x using tridiagonal method
!@+   d/dt x = d/dz (P1 d/dz x) + P4
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
!@var x the unknown to be solved (at main drid)
!@var x0 x at previous time step
!@var p1,p4 coeff. of the d.e.
!@var rhoebydz(j) rhoe(j+1)/dz(j)
!@var bydzerho(j) 1/(dze(j)*rho(j))
!@var flux_bot flux from the bottom
!@var flux_top flux at the top
!@var dtime time step
!@var n number of main grid
!@var qlimit true if tracer must be positive definite
      implicit none

      integer, intent(in) :: n
      real*8, dimension(n), intent(in) ::
     &        x0,p1,p4,rhoebydz,bydzerho
      real*8, intent(in) :: flux_bot,flux_top,dtime
      real*8, dimension(n), intent(out) :: x

      real*8, dimension(n) :: sub,dia,sup,rhs
      real*8 :: alpha
      integer :: j  !@var j loop variable
      logical, intent(in) :: qlimit 

      !sub(j)*x_jm1_kp1+dia(j)*x_j_kp1+sup(j)*x_jp1_kp1 = rhs(j)
      !k refers to time step, j refers to main grid
      !except p1(j) which is defined on the edge grid

      do j=2,n-1
          sub(j)=-dtime*p1(j)*rhoebydz(j)*bydzerho(j)
          sup(j)=-dtime*p1(j+1)*rhoebydz(j+1)*bydzerho(j)
          dia(j)=1.d0-(sub(j)+sup(j))
          rhs(j)=x0(j)+dtime*p4(j)
          if (qlimit.and. rhs(j).lt.0) rhs(j)=0.  ! prevent roundoff error
      end do

      ! Lower boundary conditions(x=T) :
      ! d/dt T = -(1/rho)*d/dz(rho*wt)
      ! d/dt T = (      ! -d/dz(rho*wt)=-(rhoe(2)*wt(2)-rhoe(1)*wt(1))/dze(1)
      ! wt(2)=-p1(2)*(      ! wt(1)=-tvflx, therefore,flux_bot(the quantity used below)
      !       flux_bot=rhoe(1)*tvflx+rhoe(2)*wt_nl(2)
      
      alpha=dtime*p1(2)*rhoebydz(2)*bydzerho(1)
      dia(1)=1.d0+alpha
      sup(1)=-alpha
      rhs(1)=x0(1)-dtime*bydzerho(1)*flux_bot

      ! Upper boundary conditions:

      ! Upper boundary conditions(x=T) :
      ! d/dt T = -(1/rho)*d/dz(rho*wt)
      ! d/dt T = (      ! -d/dz(rho*wt)=-(rhoe(n+1)*wt(n+1)-rhoe(n)*wt(n))/dze(n)
      ! wt(n)=-p1(n)*(      ! wt(n+1)=0, therefore,flux_top
      ! flux_top=0.

      alpha=dtime*p1(n)*rhoebydz(n)*bydzerho(n)
      sub(n)=-alpha
      dia(n)=1.d0+alpha
      rhs(n)=x0(n)+dtime*bydzerho(n)*flux_top

      call tridiag(sub,dia,sup,rhs,x,n)

      return
      end subroutine de_solver_main


      subroutine de_solver_edge(x,x0,p1,p3,p4,,2
     &    rhobydze,bydzrhoe,x_surf,dtime,n)
!@sum differential eqn solver for x using tridiagonal method
!@+   d/dt x = d/dz (P1 d/dz x) - P3 x + P4
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
!@var x the unknown to be solved (at edge drid)
!@var x0 x at previous time step
!@var p1,p3,p4 coeff. of the d.e.
!@var rhobydze(j) rho(j)/dze(j)
!@var bydzrhoe(j) 1/(dz(j-1)*rhoe(j))
!@var x_surf surface value of x
!@var dtime time step
!@var n number of vertical edge grid

      use CONSTANT, only : teeny
      implicit none

      integer, intent(in) :: n
      real*8, dimension(n), intent(in) ::
     &        x0,p1,p3,p4,rhobydze,bydzrhoe
      real*8, intent(in) :: x_surf,dtime
      real*8, dimension(n), intent(out) :: x

      real*8, dimension(n) :: sub,dia,sup,rhs
      integer :: j  !@var j loop variable

      ! sub(j)*x_jm1_kp1+dia(j)*x_j_kp1+sup(j)*x_jp1_kp1 = rhs(j)
      ! k refers to time step, j refers to edge grid
      ! except p1(j) which is defined on the main grid

      do j=2,n-1
          sub(j)=-dtime*p1(j-1)*rhobydze(j-1)*bydzrhoe(j)
          sup(j)=-dtime*p1(j)*rhobydze(j)*bydzrhoe(j)
          dia(j)=1.d0-(sub(j)+sup(j))+dtime*p3(j)
          rhs(j)=x0(j)+dtime*p4(j)
      end do

      ! Boundary conditions:

      dia(1)=1.d0
      sup(1)=0.d0
      rhs(1)=x_surf

      sub(n)=-1.d0
      dia(n)=1.d0
      rhs(n)=0.d0

      call tridiag(sub,dia,sup,rhs,x,n)
      ! assuming x is a non-negative scalar:
      do j=1,n
         if(      end do

      return
      end subroutine de_solver_edge


      subroutine ave_uv_to_agrid(u,v,u_a,v_a,im,jm,lm) 1,1
!@sum Computes u_a,v_a from u and v
!@var u x-component at secondary grids (B_grid)
!@var v y-component at secondary grids (B_grid)
!@var u_a x-component at primary grids (A_grid)
!@var v_a y-component at primary grids (A_grid)
!@auth Ye Cheng
!@ver  1.0
      USE GEOM, only : imaxj,idij,idjj,kmaxj,rapj,cosiv,siniv
      implicit none

      integer, intent(in) :: im,jm,lm
      real*8, dimension(im,jm,lm), intent(in) ::  u,v
      real*8, dimension(lm,im,jm), intent(out) :: u_a,v_a

      real*8, dimension(im) :: ra
      integer, dimension(im) :: idj
      real*8 :: HEMI,u_t,v_t,rak,ck,sk,uk,vk
      integer :: i,j,l,k,idik,idjk,kmax

!     polar boxes
      DO J=1,JM,JM-1
        KMAX=KMAXJ(J)
        HEMI=1.
        IF(J.LE.JM/2) HEMI=-1.
!$OMP  PARALLEL DO PRIVATE (I,L,u_t,v_t,K,IDIK,IDJK,RAK,ck,sk,uk,vk)
        DO I=1,IMAXJ(J)
          DO L=1,LM
            u_t=0.d0; v_t=0.d0
            DO K=1,KMAX
              IDIK=IDIJ(K,I,J)
              IDJK=IDJJ(K,J)
              RAK=RAPJ(K,J)
              ck=cosiv(k)
              sk=siniv(k)
              uk=u(idik,idjk,L)
              vk=v(idik,idjk,L)
              u_t=u_t+rak*(uk*ck-hemi*vk*sk)
              v_t=v_t+rak*(vk*ck+hemi*uk*sk)
            END DO
            u_a(l,i,j)=u_t
            v_a(l,i,j)=v_t
          END DO
        END DO
!$OMP  END PARALLEL DO
      END DO

!     non polar boxes
!$OMP  PARALLEL DO PRIVATE (J,I,L,u_t,v_t,K,KMAX,IDJ,RA,IDIK,IDJK,RAK)
!$OMP*    SCHEDULE(DYNAMIC,2)
      DO J=2,JM-1
        KMAX=KMAXJ(J)
        DO K=1,KMAX
          IDJ(K)=IDJJ(K,J)
          RA(K)=RAPJ(K,J)
        END DO
        DO I=1,IMAXJ(J)
          DO L=1,LM
            u_t=0.d0; v_t=0.d0
            DO K=1,KMAX
              IDIK=IDIJ(K,I,J)
              IDJK=IDJ(K)
              RAK=RA(K)
              u_t=u_t+u(IDIK,IDJK,L)*RAK
              v_t=v_t+v(IDIK,IDJK,L)*RAK
            END DO
            u_a(l,i,j)=u_t
            v_a(l,i,j)=v_t
          END DO
        END DO
      END DO
!$OMP  END PARALLEL DO
C****
      return
      end subroutine ave_uv_to_agrid


      subroutine ave_uv_to_bgrid(u_a,v_a,u,v,im,jm,lm) 1,1
!@sum Computes u and v from u_a and v_a
!@var u_a x-component of wind at primary grids (A_grid)
!@var v_a y-component of wind at primary grids (A_grid)
!@var u x-component of wind at secondary grids (B_grid)
!@var v y-component of wind at secondary grids (B_grid)
!@auth Ye Cheng
!@ver  1.0

      USE GEOM, only : imaxj,idij,idjj,kmaxj,ravj,cosiv,siniv

      implicit none

      integer, intent(in) :: im,jm,lm
      real*8, dimension(lm,im,jm), intent(in) ::  u_a,v_a
      real*8, dimension(im,jm,lm), intent(out) :: u,v

      real*8, dimension(im) :: ra
      integer, dimension(im) :: idj
      real*8 :: HEMI,rak,ck,sk,uk,vk
      integer :: i,j,l,k,idik,idjk,kmax

      u=0.d0; v=0.d0
!     polar boxes
      DO J=1,JM,JM-1
        KMAX=KMAXJ(J)
        HEMI=1.
        IF(J.LE.JM/2) HEMI=-1.
        DO I=1,IMAXJ(J)
        DO L=1,LM
        DO K=1,KMAX
              IDIK=IDIJ(K,I,J)
              IDJK=IDJJ(K,J)
              RAK=RAVJ(K,J)
              ck=cosiv(k)
              sk=siniv(k)
              uk=u_a(L,I,J)
              vk=v_a(L,I,J)
          U(IDIK,IDJK,L)=U(IDIK,IDJK,L)+RAK*(UK*CK+VK*SK*HEMI)
          V(IDIK,IDJK,L)=V(IDIK,IDJK,L)+RAK*(VK*CK-UK*SK*HEMI)
        END DO
        END DO
        END DO
      END DO
!     non polar boxes
      DO J=2,JM-1
        KMAX=KMAXJ(J)
        DO K=1,KMAX
          IDJ(K)=IDJJ(K,J)
          RA(K)=RAVJ(K,J)
        END DO
        DO I=1,IMAXJ(J)
        DO L=1,LM
        DO K=1,KMAX
          IDIK=IDIJ(K,I,J)
          IDJK=IDJ(K)
          RAK=RA(K)
          U(IDIK,IDJK,L)=U(IDIK,IDJK,L)+RAK*U_A(L,I,J)
          V(IDIK,IDJK,L)=V(IDIK,IDJK,L)+RAK*V_A(L,I,J)
        END DO
        END DO
        END DO
      END DO

      return
      end subroutine ave_uv_to_bgrid


      subroutine ave_s_to_bgrid(s1,s2,s3,s4,s5, 1,1
     &  sb1,sb2,sb3,sb4,sb5,im,jm,lm)
!@sum Computes s_b from s
!@var s  scalar at primary grid (A_grid)
!@var sb scalar at secondary grid (B_grid)
!@auth Ye Cheng
!@ver  1.0

      USE GEOM, only : imaxj,idij,idjj,kmaxj,ravj

      implicit none

      integer, intent(in) :: im,jm,lm
      real*8, dimension(lm,im,jm), intent(in) ::  s1,s2,s3,s4,s5
      real*8, dimension(lm,im,jm), intent(out) :: sb1,sb2,sb3,sb4,sb5

      real*8, dimension(im) :: ra
      integer, dimension(im) :: idj
      integer :: idik,idjk,kmax
      integer :: i,j,l,k
      real*8 :: rak

      SB1=0.;SB2=0.;SB3=0.;SB4=0.;SB5=0.
      DO J=1,JM
        KMAX=KMAXJ(J)
        DO K=1,KMAX
          IDJ(K)=IDJJ(K,J)
          RA(K)=RAVJ(K,J)
        END DO
        DO I=1,IMAXJ(J)
        DO L=1,LM
        DO K=1,KMAX
          IDIK=IDIJ(K,I,J)
          IDJK=IDJ(K)
          RAK=RA(K)
          SB1(L,IDIK,IDJK)=SB1(L,IDIK,IDJK)+RAK*S1(L,I,J)
          SB2(L,IDIK,IDJK)=SB2(L,IDIK,IDJK)+RAK*S2(L,I,J)
          SB3(L,IDIK,IDJK)=SB3(L,IDIK,IDJK)+RAK*S3(L,I,J)
          SB4(L,IDIK,IDJK)=SB4(L,IDIK,IDJK)+RAK*S4(L,I,J)
          SB5(L,IDIK,IDJK)=SB5(L,IDIK,IDJK)+RAK*S5(L,I,J)
        END DO
        END DO
        END DO
      END DO
      return
      end subroutine ave_s_to_bgrid


      subroutine apply_fluxes_to_atm 1
!@sum dummy subroutine - replaces the real one needed by DRYCNV
!@auth I. Aleinov
!@ver  1.0
      return
      end subroutine apply_fluxes_to_atm


      subroutine zze(dze,ze,n) 1,1
!@sum finds the layer edge height ze
!@var  ze vertical coordinate associated with SIGE(l)
!@var  dze(l) ze(l+1) - ze(l)
!@var  n number of layers

      USE SOCPBL, only : zgs

      implicit none

      integer, intent(in) :: n
      real*8, dimension(n), intent(in) :: dze
      real*8, dimension(n+1), intent(out) :: ze
      integer :: j  !@var j loop variable

      ze(1)=zgs
      do j=2,n+1
          ze(j)=ze(j-1)+dze(j-1)
      end do
      return
      end subroutine zze



      subroutine l_gcm(ze,lscale,n) 1,2
!@sum l_gcm calculates the turbulent length scale
!@auth Ye Cheng/G. Hartke
!@ver  1.0
!@var ze height (meters) of layer edge
!@var lscale turbulent length scale
!@var n number of layers

      USE CONSTANT, only : teeny
      USE SOCPBL, only : kappa

      implicit none

      integer, intent(in) :: n
      real*8, dimension(n+1), intent(in) :: ze
      real*8, dimension(n), intent(out) :: lscale

      real*8 :: l0,kz
      integer :: j  !@var j loop variable

      !@ Ref: Holtslag and Boville 1993, J. Climate, 6, 1825-1842.
      do j=1,n
        l0=30.+270.*exp(1.-1.d-3*ze(j))
        kz=kappa*ze(j)
        lscale(j)=l0*kz/(l0+kz)
      end do

      return
      end subroutine l_gcm


      subroutine k_gcm(tvflx,qflx,ustar,wstar,dbl 1,2
     &  ,ze,lscale,e,qturb,an2,as2,dtdz,dqdz,dudz,dvdz
     &  ,kh,km,ke,wt,wq,uw,vw,wt_nl,wq_nl

     &  ,n)

!@sum k_gcm computes the turbulent stability functions Km, Kc
!@+   and the non-local part of the fluxes
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
!@var tvflx virtual potential temperature flux at surface
!@var qflx moisture flux at surface
!@var ustar friction velocity
!@var wstar Deardorff convective velocity scale
!@var dbl the height of the PBL (real*8, in meters)
!@var kh turbulent diffusivity for scalars (heat, moisture,...)
!@var ze height (in meters) of layer edges
!@var lscale turbulent length scale
!@var e turbulent kinetic energy
!@var qturb sqrt(2*e)
!@var an2 g_alpha*dTdz
!@var as2 (dUdz)**2+(dVdz)**2
!@var dtdz dt/dz
!@var dqdz dq/dz
!@var dudz du/dz
!@var dvdz dv/dz
!@var km turbulent diffusivity for momentun
!@var ke turbulent diffusivity for e
!@var wt,wq,uw,vw turbulent fluxes
!@var wt_nl non-local part of heat flux wt
!@var wq_nl non-local part of moisture flux wq
!@var wc_nl non-local part of tracer flux
!@var n number of layers

      USE CONSTANT, only : teeny,by3
      USE SOCPBL, only : kappa,prt,ghmin,ghmax,d1,d2,d3,d4,d5
     &                  ,s0,s1,s2,s4,s5,s6,b1,g5

      implicit none

      integer, intent(in) :: n


      real*8, intent(in) :: tvflx,qflx,ustar,wstar,dbl
      real*8, dimension(n), intent(in) :: lscale,e,qturb,an2,as2
     &        ,dtdz,dqdz,dudz,dvdz
      real*8, dimension(n+1), intent(in) :: ze
      real*8, dimension(n), intent(out) ::
     &  kh,km,ke,wt,wq,uw,vw,wt_nl,wq_nl

      real*8, parameter :: kmmin=1.5d-5,khmin=2.5d-5,kmax=600.d0
      real*8 :: tmp0,tmp,tau,gm,gh,gmmax,byden,sm,sh
     &    ,ustar2,wstar3,zzi,tau_pt,w2j
      integer :: j  !@var j loop variable
      
      !@ Ref: Holtslag and Moeng 1991, JAS, 48, 1690-1698.
      ustar2=ustar*ustar
      wstar3=wstar*wstar*wstar

      tmp0=-wstar/(g5*dbl)
      do j=1,n
          tau=b1*lscale(j)/(qturb(j)+teeny)
          gh=tau*tau*an2(j)
          gm=tau*tau*as2(j)
          if(gh.lt.ghmin) gh=ghmin
          if(gh.gt.ghmax) gh=ghmax
          gmmax=(1.+d1*gh+d3*gh*gh)/(d2+d4*gh)
          if(gm.gt.gmmax) gm=gmmax
          byden=1./(1.+d1*gh+d2*gm+d3*gh*gh+d4*gh*gm+d5*gm*gm)
          sm=(s0+s1*gh+s2*gm)*byden
          km(j)=min(max(tau*e(j)*sm,kmmin),kmax)
          if(ze(j).le.dbl) then
              tau_pt=tau/g5
              zzi=ze(j)/dbl
              tmp=(1.6d0*ustar2*(1.-zzi)+teeny)**1.5d0
     &           +1.2d0*wstar3*zzi*(1.-.9d0*zzi)**1.5d0
              w2j=tmp**(2.*by3)
              kh(j)=min(max(.5d0*w2j*tau_pt,khmin),kmax)
              tmp=tmp0*tau
              wt_nl(j)=tmp*tvflx
              wq_nl(j)=tmp*qflx

          else
              sh=(s4+s5*gh+s6*gm)*byden
              kh(j)=min(max(tau*e(j)*sh,khmin),kmax)
              wt_nl(j)=0.
              wq_nl(j)=0.

          endif
          ke(j)=5.*km(j)
          wt(j) = -kh(j)*dtdz(j)+wt_nl(j)
          wq(j) = -kh(j)*dqdz(j)+wq_nl(j)
          uw(j) = -km(j)*dudz(j)
          vw(j) = -km(j)*dvdz(j)
      end do

      return
      end subroutine k_gcm


      subroutine e_gcm(tvflx,wstar,ustar,dbl,lmonin,ze,g_alpha 1,2
     &    ,an2,as2,lscale,e,n)
!@sum e_gcm finds e according to the parameterization of les data
!@+   (Moeng and Sullivan 1994) for ze<=dbl and using the giss
!@+   soc model (level 2) for ze>dbl
!@auth  Ye Cheng
!@ver   1.0
!@var (see subroutine k_gcm)
!@var lmonin Monin-Obukov length
!@var g_alpha grav*alpha
      USE CONSTANT, only : teeny,by3
      USE SOCPBL, only : kappa,emax,rimax,b1,c1,c2,c3,c4,c5

      implicit none

      integer, intent(in) :: n   !@var n  array dimension
      real*8, intent(in) :: tvflx,wstar,ustar,dbl
      real*8, dimension(n), intent(in)   :: g_alpha,an2,as2,lscale
      real*8, dimension(n+1), intent(in) :: ze
      real*8, dimension(n), intent(out) :: e
      real*8, intent(out) :: lmonin
      integer :: j !@var j loop variable
      real*8 :: ri,gm,aa,bb,cc,phi_m,tmp,wstar3,ustar3,zj,zeta

      ustar3=ustar*ustar*ustar
      wstar3=wstar*wstar*wstar
      lmonin=ustar3/(kappa*g_alpha(1)*tvflx)  ! tvflx=-(wtv)_0
      do j=1,n   ! Dyer 1974
        zj=ze(j)
        if(zj.le.dbl) then
          zeta=zj/lmonin
          if(zeta.ge.0.) then ! stable or neutral
            if(zeta.le.1.) then
              phi_m=1.+5.*zeta
            else
             phi_m=5.+zeta
            endif
          else                ! unstable
            phi_m=(1.-15.*zeta)**(-.25d0)
          endif
          tmp=.4d0*wstar3+ustar3*(dbl-zj)*phi_m/(kappa*zj)
          e(j)=min(max(tmp**(2.*by3),teeny),emax)
        else
          ri=an2(j)/max(as2(j),teeny)
          if(ri.lt.rimax) then
            aa=c1*ri*ri-c2*ri+c3
            bb=c4*ri+c5
            cc=2.d0
            if(abs(aa).lt.1d-8) then
              gm= -cc/bb
            else
              tmp=bb*bb-4.*aa*cc
              gm=(-bb-sqrt(tmp))/(2.*aa)
            endif
            tmp=0.5d0*(b1*lscale(j))**2*as2(j)/max(gm,teeny)
            e(j)=min(max(tmp,teeny),emax)
          else
            e(j)=teeny
          endif
        endif
      end do
      return
      end subroutine e_gcm


      subroutine find_pbl_top(e,ze,dbl,ldbl,n) 2
!@sum find_pbl_top finds the pbl top (at main level)
!@auth  Ye Cheng
!@ver   1.0
!@var e turbulent kinetic energy
!@var ze height at the edge level (meters)
!@var ldbl the (main) layer corresponding to top of pbl
!@var dbl the height (in meters) of the pbl, at main layer
!@+   this dbl is different from the dbl in module socpbl,
!@+   the latter is itype dependent
!@var n number of layers

      implicit none

      integer, intent(in) :: n
      real*8, dimension(n), intent(in) :: e
      real*8, dimension(n+1), intent(in) :: ze
      real*8, intent(out) :: dbl
      integer, intent(out) :: ldbl

      real*8, parameter :: dbl_max=3000.,fraction = 0.1d0
      real*8 :: e1p    ! a fraction of e(1)
      integer :: l

      e1p=fraction*e(1)
      do l=2,n
        if (      end do
      ldbl=l-1
      dbl=.5d0*(ze(l-1)+ze(l))  ! dbl is at main layer (l-1)
      dbl=min(dbl,dbl_max)
    
      return
      end subroutine find_pbl_top