MODULE SURFACE_PROCESS
	
	USE ICOSA
	USE dimphys_mod
	USE RADIATION
	DATA lmixmin,emin_turb,karman/100.,1.e-8,.4/

	contains

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	      subroutiNE vdif(ngrid,nlay,ptime,
     $                ptimestep,pcapcal,pz0,
     $                pplay,pplev,pzlay,pzlev,
     $                pu,pv,ph,ptsrf,pemis,
     $                pdufi,pdvfi,pdhfi,pfluxsrf,
     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l,
     $                lwrite)
      IMPLICIT NONE

c=======================================================================
c
c   Diffusion verticale
c   Shema implicite
c   On commence par rajouter au variables x la tendance physique
c   et on resoult en fait:
c      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
c
c   !!! attention :
c   pour utilisation sur une machine sans allocation dynamique de 
c   memoires (sur SUN par exemple) il faut que ngrid soit egal
c   a ngrid.
c
c   arguments:
c   ----------
c
c   entree:
c   -------
c
c
c=======================================================================

c-----------------------------------------------------------------------
c   declarations:
c   -------------
c
c   arguments:
c   ----------

      INTEGER ngrid,nlay
      REAL ptime,ptimestep
      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
      REAL ptsrf(ngrid),pemis(ngrid)
      REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
      REAL pfluxsrf(ngrid)
      REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
      REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid)
      REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)
      LOGICAL lwrite
c
c   local:
c   ------

      INTEGER ilev,ig,ilay,nlev
      INTEGER unit,ierr,it1,it2,icount
      SAVE icount
      INTEGER cluvdb,putdat,putvdim,setname,setvdim
      REAL z4st,zdplanck(ngrid),zu2
      REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)
      REAL zcdv(ngrid),zcdh(ngrid)
      REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
      REAL zh(ngrid,nlayermx)
      REAL ztsrf2(ngrid)
      REAL z1(ngrid),z2(ngrid)
      REAL za(ngrid,nlayermx),zb(ngrid,nlayermx)
      REAL zb0(ngrid,nlayermx)
      REAL zc(ngrid,nlayermx),zd(ngrid,nlayermx)
      REAL zout_dyn(iim+1,jjm+1,nlayermx+1),zout_fi(ngrid,nlayermx+1)
      REAL zcst1
      REAL karman

      EXTERNAL coefdifv
      EXTERNAL SSUM
      REAL SSUM
      SAVE karman

      DATA karman/0.4/
      DATA icount/0/
c
c-----------------------------------------------------------------------
c   initialisations:
c   ----------------

      nlev=nlay+1

      IF(ngrid.NE.ngrid) THEN
         PRINT*,'STOP dans coefdifv'
         PRINT*,'probleme de dimensions :'
         PRINT*,'ngrid  =',ngrid
         PRINT*,'ngrid  =',ngrid
         STOP
      ENDIF

c   computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp:
c   with rho=p/RT=p/ (R Theta) (p/ps)**kappa
c   ---------------------------------

      DO ilay=1,nlay
         DO ig=1,ngrid
            za(ig,ilay)=
     s       (pplev(ig,ilay)-pplev(ig,ilay+1))/g
         ENDDO
      ENDDO

      zcst1=4.*g*ptimestep/(kappa*cpp)**2 
      DO ilev=2,nlev-1
         DO ig=1,ngrid
            zb0(ig,ilev)=pplev(ig,ilev)*
     s      (pplev(ig,1)/pplev(ig,ilev))**kappa /
     s      (ph(ig,ilev-1)+ph(ig,ilev))
            zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/
     s      (pplay(ig,ilev-1)-pplay(ig,ilev))
         ENDDO
      ENDDO
      DO ig=1,ngrid
         zb0(ig,1)=ptimestep*pplev(ig,1)/(kappa*cpp*ptsrf(ig))
      ENDDO
      IF(lwrite) THEN
         ig=ngrid/2+1
         PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
         DO ilay=1,nlay
            WRITE(*,*) .01*pplay(ig,ilay),.001*pzlay(ig,ilay),
     s      pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
         ENDDO
         PRINT*,'Pression (mbar) ,altitude (km),zb'
         DO ilev=1,nlay
            WRITE(*,*) .01*pplev(ig,ilev),.001*pzlev(ig,ilev),
     s      zb0(ig,ilev)
         ENDDO
      ENDIF

c-----------------------------------------------------------------------
c   2. ajout des tendances physiques:
c   ------------------------------

      DO ilev=1,nlay
         DO ig=1,ngrid
            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
            zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
         ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   3. calcul de  cd :
c   ----------------
c
      CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)

c     CALL my_25(ptimestep,g,pzlev,pzlay,pu,pv,ph,zcdv,
c    a          pq2,pq2l,zkv,zkh)

        CALL vdif_k(ngrid,nlay,
     s   ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh,pq2,pq2l)

      DO ig=1,ngrid
         zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
         zcdv(ig)=zcdv(ig)*sqrt(zu2)
         zcdh(ig)=zcdh(ig)*sqrt(zu2)
      ENDDO

      IF(lwrite) THEN
         PRINT*
         PRINT*,'Diagnostique diffusion verticale'
         PRINT*,'coefficients Cd pour v et h'
         PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
         PRINT*,'coefficients K pour v et h'
         DO ilev=1,nlay
            PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
         ENDDO
      ENDIF

c-----------------------------------------------------------------------
c   integration verticale pour u:
c   -----------------------------
c
      CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
      CALL multipl(ngrid,zcdv,zb0,zb)
      DO ig=1,ngrid
         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
         zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
      ENDDO

      DO ilay=nlay-1,1,-1
         DO ig=1,ngrid
            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
            zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
         ENDDO
      ENDDO

      DO ig=1,ngrid
         zu(ig,1)=zc(ig,1)
      ENDDO
      DO ilay=2,nlay
         DO ig=1,ngrid
            zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
         ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   integration verticale pour v:
c   -----------------------------
c
      DO ig=1,ngrid
         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
         zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
      ENDDO

      DO ilay=nlay-1,1,-1
         DO ig=1,ngrid
            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
            zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
         ENDDO
      ENDDO

      DO ig=1,ngrid
         zv(ig,1)=zc(ig,1)
      ENDDO
      DO ilay=2,nlay
         DO ig=1,ngrid
            zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
         ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   integration verticale pour h:
c   -----------------------------
c
      CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
      CALL multipl(ngrid,zcdh,zb0,zb)
      DO ig=1,ngrid
         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
         zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
      ENDDO

      DO ilay=nlay-1,1,-1
         DO ig=1,ngrid
            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
            zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
         ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   rajout eventuel de planck dans le shema implicite:
c   --------------------------------------------------

      z4st=4.*5.67e-8*ptimestep
c     z4st=0.
      DO ig=1,ngrid
         zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
      ENDDO

c-----------------------------------------------------------------------
c   calcul le l'evolution de la temperature du sol':
c   -----------------------------------------------

      DO ig=1,ngrid
         z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1)
     s     +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
         z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
         ztsrf2(ig)=z1(ig)/z2(ig)
         zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
      ENDDO

c-----------------------------------------------------------------------
c   integration verticale finale:
c   -----------------------------

      DO ilay=2,nlay
         DO ig=1,ngrid
            zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
         ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   calcul final des tendances de la diffusion verticale:
c   -----------------------------------------------------

      DO ilev = 1, nlay
         DO ig=1,ngrid
            pdudif(ig,ilev)=(    zu(ig,ilev)-
     $      (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep)    )/ptimestep
            pdvdif(ig,ilev)=(    zv(ig,ilev)-
     $      (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep)    )/ptimestep
            pdhdif(ig,ilev)=(    zh(ig,ilev)-
     $      (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep)    )/ptimestep
         ENDDO
      ENDDO

      IF(lwrite) THEN
         PRINT*
         PRINT*,'Diagnostique de la diffusion verticale'
         PRINT*,'h avant et apres diffusion verticale'
         PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
         DO 3110 ilev=1,nlay
            PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev)
3110     CONTINUE
      ENDIF
c---------------------------------------------------------------------
      RETURN
      END SUBROUTINE vdif
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	      SUBROUTINE convadj(ngrid,nlay,ptimestep,
     S                   pplay,pplev,ppopsk,
     $                   pu,pv,ph,
     $                   pdufi,pdvfi,pdhfi,
     $                   pduadj,pdvadj,pdhadj)
      IMPLICIT NONE

c=======================================================================
c
c   ajustement convectif sec
c   on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement
c'
c=======================================================================

c-----------------------------------------------------------------------
c   declarations:
c   -------------
c   arguments:
c   ----------

      INTEGER ngrid,nlay
      REAL ptimestep
      REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
      REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
      REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)

c   local:
c   ------

      INTEGER ig,i,l,l1,l2,jj
      INTEGER jcnt, jadrs(ngrid)

      REAL*8 sig(nlayermx+1),sdsig(nlayermx),dsig(nlayermx)
      REAL*8 zu(ngrid,nlayermx),zv(ngrid,nlayermx)
      REAL*8 zh(ngrid,nlayermx)
      REAL*8 zu2(ngrid,nlayermx),zv2(ngrid,nlayermx)
      REAL*8 zh2(ngrid,nlayermx)
      REAL*8 zhm,zsm,zum,zvm,zalpha

      LOGICAL vtest(ngrid),down

c
c-----------------------------------------------------------------------
c   initialisation:
c   ---------------
c
      IF(ngrid.NE.ngrid) THEN
         PRINT*
         PRINT*,'STOP dans convadj'
         PRINT*,'ngrid    =',ngrid
         PRINT*,'ngrid  =',ngrid
      ENDIF
c
c-----------------------------------------------------------------------
c   detection des profils a modifier:
c   ---------------------------------
c   si le profil est a modifier
c   (i.e. ph(niv_sup) < ph(niv_inf) )
c   alors le tableau "vtest" est mis a .TRUE. ;
c   sinon, il reste a sa valeur initiale (.FALSE.)
c   cette operation est vectorisable
c   On en profite pour copier la valeur initiale de "ph"
c   dans le champ de travail "zh"


      DO 1010 l=1,nlay
         DO 1015 ig=1,ngrid
            zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
            zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
            zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
1015     CONTINUE
1010  CONTINUE

      zu2(:,:)=zu(:,:)
      zv2(:,:)=zv(:,:)
      zh2(:,:)=zh(:,:)

      DO 1020 ig=1,ngrid
        vtest(ig)=.FALSE.
 1020 CONTINUE
c
      DO 1040 l=2,nlay
        DO 1060 ig=1,ngrid
CRAY      vtest(ig)=CVMGM(.TRUE. , vtest(ig),
CRAY .                      zh2(ig,l)-zh2(ig,l-1))
          IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE.
 1060   CONTINUE
 1040 CONTINUE
c
CRAY  CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt)
      jcnt=0
      DO 1070 ig=1,ngrid
         IF(vtest(ig)) THEN
            jcnt=jcnt+1
            jadrs(jcnt)=ig
         ENDIF
1070  CONTINUE


c-----------------------------------------------------------------------
c  Ajustement des "jcnt" profils instables indices par "jadrs":
c  ------------------------------------------------------------
c
      DO 1080 jj = 1, jcnt
c
          i = jadrs(jj)
c
c   Calcul des niveaux sigma sur cette colonne
          DO l=1,nlay+1
            sig(l)=pplev(i,l)/pplev(i,1)
         ENDDO
         DO l=1,nlay
            dsig(l)=sig(l)-sig(l+1)
            sdsig(l)=ppopsk(i,l)*dsig(l)
         ENDDO
          l2 = 1
c
c      -- boucle de sondage vers le haut
c
cins$     Loop
 8000     CONTINUE
c
            l2 = l2 + 1
c
cins$       Exit
            IF (l2 .GT. nlay) Goto 8001
c
            IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN
c
c          -- l2 est le niveau le plus haut de la colonne instable
c
              l1 = l2 - 1
              l  = l1
              zsm = sdsig(l2)
              zhm = zh2(i, l2)
c
c          -- boucle de sondage vers le bas
c
cins$         Loop
 8020         CONTINUE
c
                zsm = zsm + sdsig(l)
                zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
c
c            -- doit on etendre la colonne vers le bas ?
c
c_EC (M1875) 20/6/87 : AND -> AND THEN
c
                down = .FALSE.
                IF (l1 .NE. 1) THEN    !--  and then
                  IF (zhm .LT. zh2(i, l1-1)) THEN
                    down = .TRUE.
                  END IF
                END IF
c
                IF (down) THEN
c
                  l1 = l1 - 1
                  l  = l1
c
                ELSE
c
c              -- peut on etendre la colonne vers le haut ?
c
cins$             Exit
                  IF (l2 .EQ. nlay) Goto 8021
c
cins$             Exit
                  IF (zh2(i, l2+1) .GE. zhm) Goto 8021
c
                  l2 = l2 + 1
                  l  = l2
c
                END IF
c
cins$         End Loop
              GO TO 8020
 8021         CONTINUE
c
c          -- nouveau profil : constant (valeur moyenne)
c
              zalpha=0.
              zum=0.
              zvm=0.
              DO 1100 l = l1, l2
                zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
                zh2(i, l) = zhm
                zum=zum+dsig(l)*zu(i,l)
                zvm=zvm+dsig(l)*zv(i,l)
 1100         CONTINUE
              zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
              zum=zum/(sig(l1)-sig(l2+1))
              zvm=zvm/(sig(l1)-sig(l2+1))
              IF(zalpha.GT.1.) THEN
                 PRINT*,'WARNING dans convadj zalpha=',zalpha
                 if(ig.eq.1) then
                    print*,'Au pole nord'
                 elseif (ig.eq.ngrid) then
                    print*,'Au pole sud'
                 else
                    print*,'Point i=',
     .              ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1
                 endif
!                 STOP !problem with icosa pole 
                 zalpha=1.
              ELSE
c                IF(zalpha.LT.0.) STOP'zalpha=0'
                 IF(zalpha.LT.1.e-5) zalpha=1.e-5
              ENDIF
              DO l=l1,l2
                 zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
                 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
              ENDDO
 
              l2 = l2 + 1
c
            END IF
c
cins$     End Loop
          GO TO 8000
 8001     CONTINUE
c
 1080 CONTINUE
c
      DO 4000 l=1,nlay
        DO 4020 ig=1,ngrid
          pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
          pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
          pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
 4020   CONTINUE
 4000 CONTINUE
c
      RETURN
      END SUBROUTINE convadj 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	      SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i,
     s          ptimestep,ptsrf,ptsoil,
     s          pcapcal,pfluxgrd)
      IMPLICIT NONE

c=======================================================================
c
c   Auteur:  Frederic Hourdin     30/01/92
c   -------
c
c   objet:  computation of : the soil temperature evolution
c   ------                   the surfacic heat capacity "Capcal"
c                            the surface conduction flux pcapcal
c
c
c   Method: implicit time integration
c   -------
c   Consecutive ground temperatures are related by:
c           T(k+1) = C(k) + D(k)*T(k)  (1)
c   the coefficients C and D are computed at the t-dt time-step.
c   Routine structure:
c   1)new temperatures are computed  using (1)
c   2)C and D coefficients are computed from the new temperature
c     profile for the t+dt time-step
c   3)the coefficients A and B are computed where the diffusive
c     fluxes at the t+dt time-step is given by
c            Fdiff = A + B Ts(t+dt)
c     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
c            with F0 = A + B (Ts(t))
c                 Capcal = B*dt
c           
c   Interface:
c   ----------
c
c   Arguments:
c   ----------
c   ngird               number of grid-points
c   ptimestep              physical timestep (s)
c   pto(ngrid,nsoil)     temperature at time-step t (K)
c   ptn(ngrid,nsoil)     temperature at time step t+dt (K)
c   pcapcal(ngrid)      specific heat (W*m-2*s*K-1)
c   pfluxgrd(ngrid)      surface diffusive flux from ground (Wm-2)
c   
c=======================================================================
c   declarations:
c   -------------
c-----------------------------------------------------------------------
c  arguments
c  ---------

      INTEGER ngrid,nsoil
      REAL ptimestep
      REAL ptsrf(ngrid),ptsoil(ngrid,nsoilmx),ptherm_i(ngrid)
      REAL pcapcal(ngrid),pfluxgrd(ngrid)
      LOGICAL firstcall

c-----------------------------------------------------------------------
c  local arrays
c  ------------

      INTEGER ig,jk
      REAL za(ngrid),zb(ngrid)
      REAL zdz2(nsoilmx),z1(ngrid)
      REAL min_period,dalph_soil

c   local saved variables:
c   ----------------------
      REAL dz1(nsoilmx),dz2(nsoilmx)
      REAL zc(ngrid,nsoilmx),zd(ngrid,nsoilmx)
      REAL lambda

!!!!!!!! SARVESH !!!!!!! SAVE ATTRIBUTE 
!!      SAVE dz1,dz2,zc,zd,lambda

c-----------------------------------------------------------------------
c   Depthts:
c   --------

      REAL fz,rk,fz1,rk1,rk2
      fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)

      IF (firstcall) THEN

c-----------------------------------------------------------------------
c   ground levels 
c   grnd=z/l where l is the skin depth of the diurnal cycle:
c   --------------------------------------------------------

         min_period=20000.
         dalph_soil=2.

         OPEN(99,file='soil.def',status='old',form='formatted',err=9999)
         READ(99,*) min_period
         READ(99,*) dalph_soil
         PRINT*,'Discretization for the soil model'
         PRINT*,'First level e-folding depth',min_period,
     s   '   dalph',dalph_soil
         CLOSE(99)
9999     CONTINUE

c   la premiere couche represente un dixieme de cycle diurne
         fz1=sqrt(min_period/3.14)

         DO jk=1,nsoil
            rk1=jk
            rk2=jk-1
            dz2(jk)=fz(rk1)-fz(rk2)
         ENDDO
         DO jk=1,nsoil-1
            rk1=jk+.5
            rk2=jk-.5
            dz1(jk)=1./(fz(rk1)-fz(rk2))
         ENDDO
         lambda=fz(.5)*dz1(1)
         PRINT*,'full layers, intermediate layers (secoonds)'
         DO jk=1,nsoil
            rk=jk
            rk1=jk+.5
            rk2=jk-.5
            PRINT*,fz(rk1)*fz(rk2)*3.14,
     s      fz(rk)*fz(rk)*3.14
         ENDDO

c   Initialisations:
c   ----------------

      ELSE
c-----------------------------------------------------------------------
c   Computation of the soil temperatures using the Cgrd and Dgrd
c  coefficient computed at the previous time-step:
c  -----------------------------------------------

c    surface temperature
         DO ig=1,ngrid
            ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/
     s      (lambda*(1.-zd(ig,1))+1.)
         ENDDO

c   other temperatures
         DO jk=1,nsoil-1
            DO ig=1,ngrid
               ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
            ENDDO
         ENDDO

      ENDIF
c-----------------------------------------------------------------------
c   Computation of the Cgrd and Dgrd coefficient for the next step:
c   ---------------------------------------------------------------

      DO jk=1,nsoil
         zdz2(jk)=dz2(jk)/ptimestep
      ENDDO

      DO ig=1,ngrid
         z1(ig)=zdz2(nsoil)+dz1(nsoil-1)
         zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig)
         zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig)
      ENDDO

      DO jk=nsoil-1,2,-1
         DO ig=1,ngrid
            z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
            zc(ig,jk-1)=
     s      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)
            zd(ig,jk-1)=dz1(jk-1)*z1(ig)
         ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   computation of the surface diffusive flux from ground and
c   calorific capacity of the ground:
c   ---------------------------------

      DO ig=1,ngrid
         pfluxgrd(ig)=ptherm_i(ig)*dz1(1)*
     s   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
         pcapcal(ig)=ptherm_i(ig)*
     s   (dz2(1)+ptimestep*(1.-zd(ig,1))*dz1(1))
         z1(ig)=lambda*(1.-zd(ig,1))+1.
         pcapcal(ig)=pcapcal(ig)/z1(ig)
         pfluxgrd(ig)=pfluxgrd(ig)
     s   +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig))
     s   /ptimestep
      ENDDO

      RETURN
      END SUBROUTINE SOIL
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	      SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
      IMPLICIT NONE
c=======================================================================
c
c   Subject: computation of the surface drag coefficient using the
c   -------  approch developed by Loui for ECMWF.
c
c   Author: Frederic Hourdin  15 /10 /93
c   -------
c
c   Arguments:
c   ----------
c
c   inputs:
c   ------
c     ngrid            size of the horizontal grid
c     pg               gravity (m s -2)
c     pz(ngrid)        height of the first atmospheric layer
c     pu(ngrid)        u component of the wind in that layer
c     pv(ngrid)        v component of the wind in that layer
c     pts(ngrid)       surfacte temperature
c     ph(ngrid)        potential temperature T*(p/ps)^kappa
c
c   outputs:
c   --------
c     pcdv(ngrid)      Cd for the wind
c     pcdh(ngrid)      Cd for potential temperature
c
c=======================================================================
c
c-----------------------------------------------------------------------
c   Declarations:
c   -------------

c   Arguments:
c   ----------

      INTEGER ngrid,nlay
      REAL pz0(ngrid)
      REAL pg,pz(ngrid)
      REAL pu(ngrid),pv(ngrid)
      REAL pts(ngrid),ph(ngrid)
      REAL pcdv(ngrid),pcdh(ngrid)

c   Local:
c   ------

      INTEGER ig

      REAL zu2,z1,zri,zcd0,zz

      REAL karman,b,c,d,c2b,c3bc,c3b,z0,umin2
      LOGICAL firstcal
      DATA karman,b,c,d,umin2/.4,5.,5.,5.,1.e-12/
      DATA firstcal/.true./
      SAVE b,c,d,karman,c2b,c3bc,c3b,firstcal,umin2

c-----------------------------------------------------------------------
c   couche de surface:
c   ------------------

c     DO ig=1,ngrid
c        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
c        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
c        pcdh(ig)=pcdv(ig)
c     ENDDO
c     RETURN

      IF (firstcal) THEN
         c2b=2.*b
         c3bc=3.*b*c
         c3b=3.*b
         firstcal=.false.
      ENDIF

c!!!! WARNING, verifier la formule originale de Louis!
      DO ig=1,ngrid
         zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
         zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)
         z1=1.+pz(ig)/pz0(ig)
         zcd0=karman/log(z1)
         zcd0=zcd0*zcd0*sqrt(zu2)
         IF(zri.LT.0.) THEN
            z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))
            pcdv(ig)=zcd0*(1.-2.*z1)
            pcdh(ig)=zcd0*(1.-3.*z1)
         ELSE
            zz=sqrt(1.+d*zri)
            pcdv(ig)=zcd0/(1.+c2b*zri/zz)
            pcdh(ig)=zcd0/(1.+c3b*zri*zz)
         ENDIF
      ENDDO

c-----------------------------------------------------------------------

      RETURN
      END SUBROUTINE vdif_cd
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        SUBROUTINE vdif_k(ngrid,nlay,
     s   ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pcdv,pkv,pkh,pq2,pq2l)

      IMPLICIT NONE

      INTEGER ngrid,nlay

      REAL ptimestep
      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
      REAL pz0(ngrid)
      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
      REAL pg,pcdv(ngrid)
      REAL pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)
	REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1) !!!! SARVESH ADDED to 

      INTEGER ig,il
      REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix

      REAL karman
      SAVE karman
!      DATA lmixmin,emin_turb,karman/100.,1.e-8,.4/
!!!!! SARVESH !!!!!!
!Error: Host associated variable 'lmixmin' may not be in the DATA statement 

!      print*,'LMIXMIN',lmixmin
      DO ig=1,ngrid
         pkv(ig,1)=0.
         pkh(ig,1)=0.
         pkv(ig,nlay+1)=0.
         pkh(ig,nlay+1)=0.
      ENDDO
c    s     ' zdu,zdv,zdz,zdovdz2,ph(ig,il)+ph(ig,il-1)'
      DO il=2,nlay
         DO ig=1,ngrid
            z1=pzlev(ig,il)+pz0(ig)
            lmix=karman*z1/(1.+karman*z1/lmixmin)
c           lmix=lmixmin
c WARNING test lmix=lmixmin
            zdu=pu(ig,il)-pu(ig,il-1)
            zdv=pv(ig,il)-pv(ig,il-1)
            zdz=pzlay(ig,il)-pzlay(ig,il-1)
            zdvodz2=(zdu*zdu+zdv*zdv)/(zdz*zdz)
            IF(zdvodz2.LT.1.e-5) THEN
                pkv(ig,il)=lmix*sqrt(emin_turb)
            ELSE
               zri=2.*pg*(ph(ig,il)-ph(ig,il-1))
     s         / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2  )
               pkv(ig,il)=
     s         lmix*sqrt(MAX(lmix*lmix*zdvodz2*(1-zri/.4),emin_turb))
            ENDIF
            pkh(ig,il)=pkv(ig,il)
c           IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il),
c    s      zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1),
c    s      lmix*lmix*zdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1),
c    s      ph(ig,il),ph(ig,il-1)
         ENDDO
      ENDDO

      RETURN
      END SUBROUTINE vdif_k
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	      SUBROUTINE multipl(n,x1,x2,y)
      IMPLICIT NONE
c====================================================================
c
c   multiplication de deux vecteurs
c
c=======================================================================
c
      INTEGER n,i
      REAL x1(n),x2(n),y(n)
c
      DO 10 i=1,n
	 y(i)=x1(i)*x2(i)
10    CONTINUE
c
      RETURN
      END SUBROUTINE multipl
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        END MODULE SURFACE_PROCESS