SUBROUTINE 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)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !
  !  surface layer 
  !
  !  ammended to use the "effective roughness" of mason (1986, see
  !  georgelin et al., mwr july 1994), by fm, rw, june 1995
  !
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INTEGER, PARAMETER :: lp1=lm+1
  !
  INTEGER, DIMENSION(im,jm)            , INTENT(IN   ) :: lmh
  REAL   , DIMENSION(im,jm)            , INTENT(IN   ) :: sm
  REAL   , DIMENSION(im,jm)            , INTENT(IN   ) :: ths
  REAL   , DIMENSION(im,jm)            , INTENT(IN   ) :: qs
  REAL   , DIMENSION(0:im+1,0:jm+1)    , INTENT(INOUT) :: uz0h
  REAL   , DIMENSION(0:im+1,0:jm+1)    , INTENT(INOUT) :: vz0h
  REAL   , DIMENSION(im,jm)            , INTENT(INOUT) :: thz0
  REAL   , DIMENSION(im,jm)            , INTENT(INOUT) :: qz0
  REAL   , DIMENSION(im,jm)            , INTENT(INOUT) :: ustar
  REAL   , DIMENSION(im,jm)            , INTENT(INOUT) :: z0
  REAL   , DIMENSION(im,jm,4)          , INTENT(IN   ) :: zeffij
  REAL   , DIMENSION(0:im+1,0:jm+1)    , INTENT(INOUT) :: akms
  REAL   , DIMENSION(im,jm)            , INTENT(INOUT) :: akhs
  REAL   , DIMENSION(im,jm)            , INTENT(IN   ) :: hpblwork
  REAL   , DIMENSION(im,jm)            , INTENT(OUT  ) :: ct
  REAL   , DIMENSION(im,jm)            , INTENT(OUT  ) :: u10
  REAL   , DIMENSION(im,jm)            , INTENT(OUT  ) :: v10
  REAL   , DIMENSION(im,jm)            , INTENT(OUT  ) :: tshltr
  REAL   , DIMENSION(im,jm)            , INTENT(OUT  ) :: th10
  REAL   , DIMENSION(im,jm)            , INTENT(OUT  ) :: qshltr
  REAL   , DIMENSION(im,jm)            , INTENT(OUT  ) :: q10
  REAL   , DIMENSION(im,jm,lp1)        , INTENT(INOUT) :: uwork
  REAL   , DIMENSION(im,jm,lp1)        , INTENT(IN   ) :: vwork
  REAL   , DIMENSION(0:im+1,0:jm+1,lm) , INTENT(IN   ) :: t
  REAL   , DIMENSION(0:im+1,0:jm+1,lm) , INTENT(IN   ) :: q
  REAL   , DIMENSION(im,jm,lm)         , INTENT(IN   ) :: ape
  REAL   , DIMENSION(0:im+1,0:jm+1,lp1), INTENT(IN   ) :: zint
  !
  REAL, PARAMETER :: wwst=1.2
  REAL, PARAMETER :: wwst2=wwst*wwst
  REAL, PARAMETER :: g=9.8
  REAL, PARAMETER :: ustfc=0.018/g
  REAL, PARAMETER :: vkrm=0.40
  REAL, PARAMETER :: rfc=0.191
  REAL, PARAMETER :: fhneu=0.8
  REAL, PARAMETER :: ric=0.183
  REAL, PARAMETER :: rric=1.0/ric
  REAL, PARAMETER :: rfac=ric/(fhneu*rfc*rfc)
  REAL, PARAMETER :: excm=0.001
  REAL, PARAMETER :: beta=1./270.
  REAL, PARAMETER :: btg=beta*g
  REAL, PARAMETER :: elfc=vkrm*btg
  REAL, PARAMETER :: cnv=0.608*g/btg
  REAL, PARAMETER :: wold=.15
  REAL, PARAMETER :: wnew=1.-wold
  REAL, PARAMETER :: itrmx=05
  REAL, PARAMETER :: pihf=3.14159265/2.
  REAL, PARAMETER :: pifr=3.14159265/4.
  !
  REAL, PARAMETER :: epsu2=1.e-4
  REAL, PARAMETER :: epsust=0.07
  REAL, PARAMETER :: epsit=1.e-4
  REAL, PARAMETER :: epsa=1.e-8
  REAL, PARAMETER :: ztmin=-5.
  REAL, PARAMETER :: ztmax=1.
  !
  REAL, PARAMETER :: small=0.35
  REAL, PARAMETER :: glkbs=30.0
  REAL, PARAMETER :: glkbr=10.0
  REAL, PARAMETER :: grrs=glkbr/glkbs
  REAL, PARAMETER :: cziv=small*glkbs
  REAL, PARAMETER :: visc=1.5e-5
  REAL, PARAMETER :: tvisc=2.1e-5
  REAL, PARAMETER :: qvisc=2.1e-5
  REAL, PARAMETER :: rvisc=1./visc
  REAL, PARAMETER :: rtvisc=1./tvisc
  REAL, PARAMETER :: rqvisc=1./qvisc
  REAL, PARAMETER :: sqpr=0.84
  REAL, PARAMETER :: sqsc=0.84
  REAL, PARAMETER :: zqrzt=sqsc/sqpr
  REAL, PARAMETER :: ustr=0.225
  REAL, PARAMETER :: ustc=0.7
  REAL, PARAMETER :: fzu1=cziv*visc
  REAL, PARAMETER :: fzt1=rvisc *tvisc*sqpr
  REAL, PARAMETER :: fzq1=rtvisc*qvisc*zqrzt
  REAL, PARAMETER :: fzt2=cziv*grrs*tvisc*sqpr
  REAL, PARAMETER :: fzq2=rtvisc*qvisc*zqrzt
  !
  REAL, PARAMETER :: ztfc=1.0
  REAL, PARAMETER :: czil=.1000
  REAL, PARAMETER :: sqvisc=258.2
  REAL, PARAMETER :: zilfc=-czil*vkrm*sqvisc
  !
  INTEGER :: itr, i, j
  LOGICAL, DIMENSION(im,jm) :: ExtIfTrue
  LOGICAL, DIMENSION(im,jm) :: IntIfTrue
  INTEGER, DIMENSION(im,jm) :: ml
  INTEGER, DIMENSION(im,jm) :: mh
  REAL   , DIMENSION(im,jm) :: thlm
  REAL   , DIMENSION(im,jm) :: qlm
  REAL   , DIMENSION(im,jm) :: zu
  REAL   , DIMENSION(im,jm) :: wght
  REAL   , DIMENSION(im,jm) :: rwgh
  REAL   , DIMENSION(im,jm) :: dthv
  REAL   , DIMENSION(im,jm) :: du2
  REAL   , DIMENSION(im,jm) :: rib
  REAL   , DIMENSION(im,jm) :: zt
  REAL   , DIMENSION(im,jm) :: zq
  REAL   , DIMENSION(im,jm) :: alpha
  REAL   , DIMENSION(im,jm) :: x
  REAL   , DIMENSION(im,jm) :: wlow
  REAL   , DIMENSION(im,jm) :: zsl
  REAL   , DIMENSION(im,jm) :: zslu
  REAL   , DIMENSION(im,jm) :: zslt
  REAL   , DIMENSION(im,jm) :: btgh
  REAL   , DIMENSION(im,jm) :: wstar2
  REAL   , DIMENSION(im,jm) :: rdz
  REAL   , DIMENSION(im,jm) :: cxch
  REAL   , DIMENSION(im,jm) :: rlmo
  REAL   , DIMENSION(im,jm) :: rlmp
  REAL   , DIMENSION(im,jm) :: rlmn
  REAL   , DIMENSION(im,jm) :: rlma
  REAL   , DIMENSION(im,jm) :: rlogu
  REAL   , DIMENSION(im,jm) :: rlogt
  REAL   , DIMENSION(im,jm) :: zetalt
  REAL   , DIMENSION(im,jm) :: zetalu
  REAL   , DIMENSION(im,jm) :: zetau
  REAL   , DIMENSION(im,jm) :: zetat
  REAL   , DIMENSION(im,jm) :: psmz
  REAL   , DIMENSION(im,jm) :: simm
  REAL   , DIMENSION(im,jm) :: pshz
  REAL   , DIMENSION(im,jm) :: simh
  REAL   , DIMENSION(im,jm) :: ustark
  REAL   , DIMENSION(im,jm) :: xlu
  REAL   , DIMENSION(im,jm) :: xlu4
  REAL   , DIMENSION(im,jm) :: xlt
  REAL   , DIMENSION(im,jm) :: xlt4
  REAL   , DIMENSION(im,jm) :: xu
  REAL   , DIMENSION(im,jm) :: xu4
  REAL   , DIMENSION(im,jm) :: xt
  REAL   , DIMENSION(im,jm) :: xt4
  REAL   , DIMENSION(im,jm) :: umflx
  REAL   , DIMENSION(im,jm) :: vmflx
  REAL   , DIMENSION(im,jm) :: hsflx
  REAL   , DIMENSION(im,jm) :: hlflx
  REAL   , DIMENSION(im,jm) :: akms10
  REAL   , DIMENSION(im,jm) :: akhs02
  REAL   , DIMENSION(im,jm) :: akhs10
  REAL   , DIMENSION(im,jm) :: zu10
  REAL   , DIMENSION(im,jm) :: zt02
  REAL   , DIMENSION(im,jm) :: zt10
  REAL   , DIMENSION(im,jm) :: rlnu10
  REAL   , DIMENSION(im,jm) :: rlnt02
  REAL   , DIMENSION(im,jm) :: rlnt10
  REAL   , DIMENSION(im,jm) :: ztau10
  REAL   , DIMENSION(im,jm) :: ztat02
  REAL   , DIMENSION(im,jm) :: ztat10
  REAL   , DIMENSION(im,jm) :: simm10
  REAL   , DIMENSION(im,jm) :: simh02
  REAL   , DIMENSION(im,jm) :: simh10
  REAL   , DIMENSION(im,jm) :: xlu104
  REAL   , DIMENSION(im,jm) :: xlt024
  REAL   , DIMENSION(im,jm) :: xlt104
  REAL   , DIMENSION(im,jm) :: xlu10
  REAL   , DIMENSION(im,jm) :: xlt02
  REAL   , DIMENSION(im,jm) :: xlt10
  !
  REAL    :: xx, yy, zz
  REAL    :: pslmu, pslms, pslhu, pslhs, pspmu, pspms, psphu, psphs
  !
  pslmu(zz)=-0.96*alog(1.0-4.5*zz)
  pslms(zz)=zz*rric-2.076*(1.-1./(zz+1.))
  pslhu(zz)=-0.96*alog(1.0-4.5*zz)
  pslhs(zz)=zz*rfac-2.076*(1.-1./(zz+1.))
  pspmu(xx)=-2.*alog((xx+1.)*0.5)-alog((xx*xx+1.)*0.5)+&
       2.*ATAN(xx)-pihf
  pspms(yy)=5.*yy
  psphu(xx)=-2.*alog((xx*xx+1.)*0.5)
  psphs(yy)=5.*yy


  DO j = 3, jm-2
     DO i = 1, im

        thlm(i,j)=t(i,j,lmh(i,j))*ape(i,j,lmh(i,j))
        qlm(i,j)=q(i,j,lmh(i,j))

        z0(i,j)=(1.-sm(i,j))*z0(i,j)+&
             sm(i,j)*amax1(ustfc*ustar(i,j)*ustar(i,j),1.59e-5)

        !*** viscous sublayer

        IF(sm(i,j).GT.0.5.AND.ustar(i,j).LT.ustc)THEN

           IF(ustar(i,j).LT.ustr)THEN

              zu(i,j)=fzu1*SQRT(SQRT(z0(i,j)*ustar(i,j)*rvisc))/ustar(i,j)

              wght(i,j)=akms(i,j)*zu(i,j)*rvisc
              rwgh(i,j)=wght(i,j)/(wght(i,j)+1.)
              uz0h(i,j)=(uwork(i,j,lmh(i,j))*rwgh(i,j)+uz0h(i,j))*0.5
              vz0h(i,j)=(vwork(i,j,lmh(i,j))*rwgh(i,j)+vz0h(i,j))*0.5

              zt(i,j)=fzt1*zu(i,j)
              wght(i,j)=akhs(i,j)*zt(i,j)*rtvisc
              thz0(i,j)=((wght(i,j)*thlm(i,j)+ths(i,j))/&
                   (wght(i,j)+1.)+thz0(i,j))*0.5

              zq(i,j)=fzq1*zt(i,j)
              wght(i,j)=akhs(i,j)*zq(i,j)*rqvisc
              qz0(i,j) =((wght(i,j)*qlm(i,j)+qs(i,j))&
                   /(wght(i,j)+1.)+qz0(i,j))*0.5

           ENDIF
           IF(ustar(i,j).GE.ustr.AND.ustar(i,j).LT.ustc)THEN

              zu(i,j)=z0(i,j)
              uz0h(i,j)=0.
              vz0h(i,j)=0.

              zt(i,j)=fzt2*SQRT(SQRT(z0(i,j)*ustar(i,j)*rvisc))/ustar(i,j)
              wght(i,j)=akhs(i,j)*zt(i,j)*rtvisc
              thz0(i,j)=((wght(i,j)*thlm(i,j)+ths(i,j))/&
                   (wght(i,j)+1.)+thz0(i,j))*0.5

              zq(i,j)=fzq2*zt(i,j)
              wght(i,j)=akhs(i,j)*zq(i,j)*rqvisc
              qz0(i,j) =((wght(i,j)*qlm(i,j)+qs(i,j))/&
                   (wght(i,j)+1.)+qz0(i,j))*0.5
           ENDIF

        ELSE

           zu(i,j)=z0(i,j)

           IF(sm(i,j).LE.0.5)THEN
              IF(uwork(i,j,lmh(i,j)).EQ.0.) uwork(i,j,lmh(i,j))=epsu2
              alpha(i,j)=ABS(ATAN(vwork(i,j,lmh(i,j))/&
                   uwork(i,j,lmh(i,j)))+pihf-epsa)
              x(i,j)=alpha(i,j)/pifr
              ml(i,j)=1+x(i,j)
              mh(i,j)=1+MOD(ml(i,j),4)
              wlow(i,j)=x(i,j)-ml(i,j)+1
              zu(i,j)=wlow(i,j)*zeffij(i,j,ml(i,j))+&
                   (1.-wlow(i,j))*zeffij(i,j,mh(i,j))
           ENDIF

           uz0h(i,j)=0.
           vz0h(i,j)=0.

           zt(i,j)=z0(i,j)
           thz0(i,j)=ths(i,j)

           zq(i,j)=z0(i,j)
           qz0(i,j)=qs(i,j)

        ENDIF
     END DO
  END DO

  DO j = 3, jm-2
     DO i = 1, im

        zsl(i,j)=(zint(i,j,lmh(i,j))-zint(i,j,lmh(i,j)+1))*0.5

        zu(i,j)=amin1(zu(i,j),0.5*zsl(i,j))

        rdz(i,j)=1./zsl(i,j)
        cxch(i,j)=excm*rdz(i,j)

        IF(sm(i,j).GT.0.5)THEN
           dthv(i,j)=(0.608*qlm(i,j)+1.)*&
                thlm(i,j)-(0.608*qz0(i,j)+1.)*thz0(i,j)
        ELSE
           dthv(i,j)=(qlm(i,j)-qz0(i,j))*cnv+thlm(i,j)-thz0(i,j)
           zt(i,j)=z0(i,j)*ztfc
        ENDIF

        du2(i,j)=amax1((uwork(i,j,lmh(i,j))-uz0h(i,j))**2+&
             (vwork(i,j,lmh(i,j))-vz0h(i,j))**2,epsu2)

        rib(i,j)=btg*dthv(i,j)*zsl(i,j)/du2(i,j)

        !*** beljars correction of ustar

        btgh(i,j)=btg*hpblwork(i,j)
        wstar2(i,j)=wwst2*ABS(btgh(i,j)*akhs(i,j)*dthv(i,j))**(2./3.)
        ustar(i,j)=amax1(SQRT(akms(i,j)*SQRT(du2(i,j)+wstar2(i,j))),epsust)

        !*** zilitinkevitch fix for zt

        IF(sm(i,j).LT.0.5)zt(i,j)=&
             EXP(zilfc*SQRT(ustar(i,j)*z0(i,j)))*z0(i,j)

     END DO
  END DO

  DO j = 3, jm-2
     DO i = 1, im
        IntIfTrue(i,j) = sm(i,j)>0.5
        ExtIfTrue(i,j) = (sm(i,j)>0.5) .AND. (rib(i,j).GE.ric)
     END DO
  END DO

  DO j = 3, jm-2
     DO i = 1, im
        IF(ExtIfTrue(i,j))THEN
           akms(i,j)=amax1( visc*rdz(i,j),cxch(i,j))
           akhs(i,j)=amax1(tvisc*rdz(i,j),cxch(i,j))
        ELSE
           zslu(i,j)=zsl(i,j)+zu(i,j)
           zslt(i,j)=zsl(i,j)+zt(i,j)
           rlogu(i,j)=alog(zslu(i,j)/zu(i,j))
           rlogt(i,j)=alog(zslt(i,j)/zt(i,j))
           rlmo(i,j)=elfc*akhs(i,j)*dthv(i,j)/ustar(i,j)**3
        END IF
     END DO
  END DO

  !*** sea points first ... 

  DO itr=1,itrmx
     DO j = 3, jm-2
        DO i = 1, im
           IF ((.NOT.ExtIfTrue(i,j)) .AND. IntIfTrue(i,j)) THEN

              !*** 1./monin-obukkhov length-scale

              zetalt(i,j)=amax1(zslt(i,j)*rlmo(i,j),ztmin)
              rlmo(i,j)=zetalt(i,j)/zslt(i,j)
              zetalu(i,j)=zslu(i,j)*rlmo(i,j)

              zetau(i,j)=zu(i,j)*rlmo(i,j)
              zetat(i,j)=zt(i,j)*rlmo(i,j)

              !*** ll functions over sea

              IF(rlmo(i,j).LT.0.)THEN
                 psmz(i,j)=pslmu(zetau(i,j))
                 simm(i,j)=pslmu(zetalu(i,j))-psmz(i,j)+rlogu(i,j)
                 pshz(i,j)=pslhu(zetat(i,j))
                 simh(i,j)=fhneu*&
                      (pslhu(zetalt(i,j))-pshz(i,j)+rlogt(i,j))
              ELSE
                 psmz(i,j)=pslms(zetau(i,j))
                 simm(i,j)=pslms(zetalu(i,j))-psmz(i,j)+rlogu(i,j)
                 pshz(i,j)=pslhs(zetat(i,j))
                 simh(i,j)=fhneu*&
                      (pslhs(zetalt(i,j))-pshz(i,j)+rlogt(i,j))
              ENDIF

              !*** beljaars correction for ustar

              ustar(i,j)=amax1(&
                   SQRT(akms(i,j)*SQRT(du2(i,j)+wstar2(i,j))),epsust)

              ustark(i,j)=ustar(i,j)*vkrm
              akms(i,j)=amax1(ustark(i,j)/simm(i,j),cxch(i,j))
              akhs(i,j)=amax1(ustark(i,j)/simh(i,j),cxch(i,j))

              wstar2(i,j)=wwst2*&
                   ABS(btgh(i,j)*akhs(i,j)*dthv(i,j))**(2./3.)
              rlmn(i,j)=elfc*akhs(i,j)*dthv(i,j)/ustar(i,j)**3

              rlmp(i,j)=rlmo(i,j)
              rlma(i,j)=rlmo(i,j)*wold+rlmn(i,j)*wnew

              rlmo(i,j)=rlma(i,j)
           END IF
        END DO
     END DO
  END DO

  !*** now land points ...

  DO itr=1,itrmx
     DO j = 3, jm-2
        DO i = 1, im
           IF ((.NOT.ExtIfTrue(i,j)) .AND. (.NOT.IntIfTrue(i,j))) THEN

              !*** 1./monin-obukkhov length-scale

              zetalt(i,j)=amax1(zslt(i,j)*rlmo(i,j),ztmin)
              rlmo(i,j)=zetalt(i,j)/zslt(i,j)
              zetalu(i,j)=zslu(i,j)*rlmo(i,j)

              zetau(i,j)=zu(i,j)*rlmo(i,j)
              zetat(i,j)=zt(i,j)*rlmo(i,j)

              !*** paulson 1970 functions over land w rad. skin t

              IF(rlmo(i,j).LT.0.)THEN
                 xlu4(i,j)=1.-16.*zetalu(i,j)
                 xlt4(i,j)=1.-16.*zetalt(i,j)
                 xu4(i,j) =1.-16.*zetau(i,j)
                 xt4(i,j) =1.-16.*zetat(i,j)

                 xlu(i,j)=SQRT(SQRT(xlu4(i,j)))
                 xlt(i,j)=SQRT(SQRT(xlt4(i,j)))
                 xu(i,j) =SQRT(SQRT(xu4(i,j)))
                 xt(i,j) =SQRT(SQRT(xt4(i,j)))

                 psmz(i,j)=pspmu(xu(i,j))
                 simm(i,j)=pspmu(xlu(i,j))-psmz(i,j)+rlogu(i,j)
                 pshz(i,j)=psphu(xt(i,j))
                 simh(i,j)=psphu(xlt(i,j))-pshz(i,j)+rlogt(i,j)
              ELSE
                 zetau(i,j)=amin1(zetau(i,j),ztmax)
                 zetat(i,j)=amin1(zetat(i,j),ztmax)
                 zetalu(i,j)=amin1(zetalu(i,j),ztmax)
                 zetalt(i,j)=amin1(zetalt(i,j),ztmax)
                 psmz(i,j)=pspms(zetau(i,j))
                 simm(i,j)=pspms(zetalu(i,j))-psmz(i,j)+rlogu(i,j)
                 pshz(i,j)=psphs(zetat(i,j))
                 simh(i,j)=psphs(zetalt(i,j))-pshz(i,j)+rlogt(i,j)
              ENDIF

              !*** beljaars correction for ustar

              ustar(i,j)=amax1(&
                   SQRT(akms(i,j)*SQRT(du2(i,j)+wstar2(i,j))),epsust)

              !*** zilitinkevitch fix for zt

              zt(i,j)=EXP(zilfc*SQRT(ustar(i,j)*z0(i,j)))*z0(i,j)
              zslt(i,j)=zsl(i,j)+zt(i,j)
              rlogt(i,j)=alog(zslt(i,j)/zt(i,j))

              ustark(i,j)=ustar(i,j)*vkrm
              akms(i,j)=amax1(ustark(i,j)/simm(i,j),cxch(i,j))
              akhs(i,j)=amax1(ustark(i,j)/simh(i,j),cxch(i,j))

              wstar2(i,j)=wwst2*&
                   ABS(btgh(i,j)*akhs(i,j)*dthv(i,j))**(2./3.)
              rlmn(i,j)=elfc*akhs(i,j)*dthv(i,j)/ustar(i,j)**3

              rlmp(i,j)=rlmo(i,j)
              rlma(i,j)=rlmo(i,j)*wold+rlmn(i,j)*wnew

              rlmo(i,j)=rlma(i,j)

           END IF
        END DO
     END DO
  END DO


  DO j = 3, jm-2
     DO i = 1, im
        ct(i,j)=0.

        !*** diagnostic block

        !
        !*** Jairo: comentado pois wstar, antigo argumento de saida
        !    desta subrotina,  nao eh utilizado em turbl
        !    wstar=SQRT(wstar2)/wwst

        umflx(i,j)=akms(i,j)*(uwork(i,j,lmh(i,j)) -uz0h(i,j) )
        vmflx(i,j)=akms(i,j)*(vwork(i,j,lmh(i,j)) -vz0h(i,j) )
        hsflx(i,j)=akhs(i,j)*(thlm(i,j)-thz0(i,j))
        hlflx(i,j)=akhs(i,j)*(qlm(i,j) -qz0(i,j) )

        IF(sm(i,j).GT.0.5.AND.rib(i,j).GE.ric)THEN

           akms10(i,j)=amax1( visc/10.,cxch(i,j))
           akhs02(i,j)=amax1(tvisc/02.,cxch(i,j))
           akhs10(i,j)=amax1(tvisc/10.,cxch(i,j))

        ELSE

           zu10(i,j)=zu(i,j)+10.
           zt02(i,j)=zt(i,j)+02.
           zt10(i,j)=zt(i,j)+10.

           rlnu10(i,j)=alog(zu10(i,j)/zu(i,j))
           rlnt02(i,j)=alog(zt02(i,j)/zt(i,j))
           rlnt10(i,j)=alog(zt10(i,j)/zt(i,j))

           ztau10(i,j)=zu10(i,j)*rlmp(i,j)
           ztat02(i,j)=zt02(i,j)*rlmp(i,j)
           ztat10(i,j)=zt10(i,j)*rlmp(i,j)

           !*** ll functions over sea

           IF(sm(i,j).GT.0.5)THEN

              IF(rlmp(i,j).LT.0.)THEN
                 simm10(i,j)=       &
                      pslmu(ztau10(i,j))-psmz(i,j)+rlnu10(i,j)
                 simh02(i,j)=fhneu*&
                      (pslhu(ztat02(i,j))-pshz(i,j)+rlnt02(i,j))
                 simh10(i,j)=fhneu*&
                      (pslhu(ztat10(i,j))-pshz(i,j)+rlnt10(i,j))
              ELSE
                 simm10(i,j)=       &
                      pslms(ztau10(i,j))-psmz(i,j)+rlnu10(i,j)
                 simh02(i,j)=fhneu*&
                      (pslhs(ztat02(i,j))-pshz(i,j)+rlnt02(i,j))
                 simh10(i,j)=fhneu*&
                      (pslhs(ztat10(i,j))-pshz(i,j)+rlnt10(i,j))
              ENDIF

              !*** paulson 1970 functions over land w rad. skin t

           ELSE

              IF(rlmp(i,j).LT.0.)THEN
                 xlu104(i,j)=1.-16.*ztau10(i,j)
                 xlt024(i,j)=1.-16.*ztat02(i,j)
                 xlt104(i,j)=1.-16.*ztat10(i,j)

                 xlu10(i,j)=SQRT(SQRT(xlu104(i,j)))
                 xlt02(i,j)=SQRT(SQRT(xlt024(i,j)))
                 xlt10(i,j)=SQRT(SQRT(xlt104(i,j)))

                 simm10(i,j)=pspmu(xlu10(i,j))-psmz(i,j)+rlnu10(i,j)
                 simh02(i,j)=psphu(xlt02(i,j))-pshz(i,j)+rlnt02(i,j)
                 simh10(i,j)=psphu(xlt10(i,j))-pshz(i,j)+rlnt10(i,j)
              ELSE
                 ztau10(i,j)=amin1(ztau10(i,j),ztmax)
                 ztat02(i,j)=amin1(ztat02(i,j),ztmax)
                 ztat10(i,j)=amin1(ztat10(i,j),ztmax)

                 simm10(i,j)=pspms(ztau10(i,j))-psmz(i,j)+rlnu10(i,j)
                 simh02(i,j)=psphs(ztat02(i,j))-pshz(i,j)+rlnt02(i,j)
                 simh10(i,j)=psphs(ztat10(i,j))-pshz(i,j)+rlnt10(i,j)
              ENDIF

           ENDIF

           akms10(i,j)=amax1(ustark(i,j)/simm10(i,j),cxch(i,j))
           akhs02(i,j)=amax1(ustark(i,j)/simh02(i,j),cxch(i,j))
           akhs10(i,j)=amax1(ustark(i,j)/simh10(i,j),cxch(i,j))

        ENDIF

     END DO
  END DO


  DO j = 3, jm-2
     DO i = 1, im
        u10(i,j) =umflx(i,j)/akms10(i,j)+uz0h(i,j)
        v10(i,j) =vmflx(i,j)/akms10(i,j)+vz0h(i,j)
        tshltr(i,j)=hsflx(i,j)/akhs02(i,j)+thz0(i,j)
        th10(i,j)=hsflx(i,j)/akhs10(i,j)+thz0(i,j)
        qshltr(i,j) =hlflx(i,j)/akhs02(i,j)+qz0(i,j)
        q10(i,j) =hlflx(i,j)/akhs10(i,j)+qz0(i,j)

        IF(qshltr(i,j).LT.0.)qshltr(i,j)=qlm(i,j)
        IF(q10(i,j).LT.0.)q10(i,j)=qlm(i,j)

     END DO

  END DO
END SUBROUTINE sfcdif
