C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE SFCDIF(LMHK,SM,THS,QS,UZ0,VZ0,THZ0,QZ0
     &,                 USTAR,WSTAR,Z0,ZEFF,AKMS,AKHS,HPBL,CT
     &,                 U10,V10,TH02,TH10,Q02,Q10
CGSM v100m
     &,                 TH100,Q100,U100,V100
CGSM v100m
     &,                 ULM,VLM,T,Q,APE,Z,PD,PT,TLM)
C     ******************************************************************
C     *                                                                *
C     *                        SURFACE LAYER                           *
C     *                                                                *
C     ******************************************************************
C     Ammended to use the "Effective roughness" of Mason (1986, see
C     Georgelin et al., MWR July 1994), by FM, RW, June 1995
C-----------------------------------------------------------------------
      INCLUDE "parmeta"
      INCLUDE "mpp.h"
#include "sp.h"
C-----------------------------------------------------------------------
                             P A R A M E T E R
     &(LP1=LM+1)
C-----------------------------------------------------------------------
                             P A R A M E T E R
     &(WWST=1.2,WWST2=WWST*WWST,G=9.8,USTFC=0.018/G
     &,VKRM=0.40,RIC=0.183,RFC=0.191,FHNEU=0.8
     &,RRIC=1.0/RIC,RFAC=RIC/(FHNEU*RFC*RFC),EXCM=0.001
     &,BETA=1./270.,BTG=BETA*G
     &,ELFC=VKRM*BTG,CNV=0.608*G/BTG
     &,WOLD=.15,WNEW=1.-WOLD,ITRMX=05
     &,PIHF=3.14159265/2.,PIFR=3.14159265/4.
C-----------------------------------------------------------------------
     &,EPSU2=1.E-4,EPSUST=0.07,EPSIT=1.E-4,EPSA=1.E-8
     &,ZTMIN=-5.
C-----------------------------------------------------------------------
CMM     &,SMALL=0.35, GLKBS=30.0,GLKBR=10.0,GRRS=GLKBR/GLKBS
     &,        GLKBS=30.0,GLKBR=10.0,GRRS=GLKBR/GLKBS                
CMM     &,CZIV=SMALL*GLKBS
     &,VISC=1.5E-5, TVISC=2.1E-5, QVISC=2.1E-5
     &,RVISC=1./VISC,RTVISC=1./TVISC,RQVISC=1./QVISC
     &,SQPR=0.84,SQSC=0.84,ZQRZT=SQSC/SQPR
     &,USTR=0.225,USTC=0.7
CMM     &,FZU1=CZIV*VISC,FZT1=RVISC *TVISC*SQPR,   FZQ1=RTVISC*QVISC*ZQRZT
     &,FZT1=RVISC *TVISC*SQPR, FZQ1=RTVISC*QVISC*ZQRZT
CMM     &,               FZT2=CZIV*GRRS*TVISC*SQPR,FZQ2=RTVISC*QVISC*ZQRZT
     &,      FZQ2=RTVISC*QVISC*ZQRZT, M=30.0
C-----------------------------------------------------------------------
     &,ZTFC=1.0
c     &,CZIL=.1000,SQVISC=258.2,ZILFC=-CZIL*VKRM*SQVISC
     &,CZIL=.2000,SQVISC=258.2,ZILFC=-CZIL*VKRM*SQVISC
     &,PQ0=379.90516,A2=17.2693882,A3=273.16,A4=35.86
     &,CAPA=0.28589641E0,H1M5=1.E-5)
C-----------------------------------------------------------------------
      parameter  (ZTMAX=2.)                     !! FRAN-URIX:2011
C-----------------------------------------------------------------------     
                             D I M E N S I O N
     & T     (LM),Q     (LM)
                             D I M E N S I O N
     & APE   (LM)
     &,Z     (LP1)
C-ZEFF-ZEFF-ZEFF-ZEFF
     &,ZEFF  (4)
C-ZEFF-ZEFF-ZEFF-ZEFF
C-----------------------------------------------------------------------
C morelli COMMON MOMENTO
      COMMON /MOMENTO/ UMFLX,VMFLX,AKMS10
C morelli end
      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.))
C
      PSPMU(XX)=-2.*ALOG((XX+1.)*0.5)-ALOG((XX*XX+1.)*0.5)+2.*ATAN(XX)
     &          -PIHF
      PSPMS(YY)=4.7*YY                       !! FRAN-URIX:2011
      PSPHU(XX)=-2.*ALOG((XX*XX+1.)*0.5)
      PSPHS(YY)=4.7*YY                             !! FRAN-URIX:2011
C***********************************************************************
      LMHP=LMHK+1
C
      THLM=T(LMHK)*APE(LMHK)
      QLM=Q(LMHK)
C-----------------------------------------------------------------------
      Z0=(1.-SM)*Z0+SM*AMAX1(USTFC*USTAR*USTAR,1.59E-5)
C--------------VISCOUS SUBLAYER-----------------------------------------
      IF(SM.GT.0.5.AND.USTAR.LT.USTC)THEN                !! OCEAN CASE
C-----------------------------------------------------------------------

C--------- Malagutti class experiments-------------
         Rr= USTAR*Z0*1/VISC
         SMALL=11/(M*(Rr**0.25))
         CZIV=SMALL*GLKBS
         FZU1=CZIV*VISC
         FZT2=CZIV*GRRS*TVISC*SQPR
C        WRITE(6,*)"SMALL1=" , SMALL
C-----------------------------------------------------------------------      
         IF(USTAR.LT.USTR)THEN
C
           ZU=FZU1*SQRT(SQRT(Z0*USTAR*RVISC))/USTAR
           WGHT=AKMS*ZU*RVISC
           RWGH=WGHT/(WGHT+1.)
           UZ0=(ULM*RWGH+UZ0)*0.5
           VZ0=(VLM*RWGH+VZ0)*0.5
C
           ZT=FZT1*ZU
           WGHT=AKHS*ZT*RTVISC
           THZ0=((WGHT*THLM+THS)/(WGHT+1.)+THZ0)*0.5
C
           ZQ=FZQ1*ZT
           WGHT=AKHS*ZQ*RQVISC
           QZ0 =((WGHT*QLM+QS)/(WGHT+1.)+QZ0)*0.5
C
         ENDIF
         IF(USTAR.GE.USTR.AND.USTAR.LT.USTC)THEN
C
           ZU=Z0
           UZ0=0.
           VZ0=0.
C
           ZT=FZT2*SQRT(SQRT(Z0*USTAR*RVISC))/USTAR
           WGHT=AKHS*ZT*RTVISC
           THZ0=((WGHT*THLM+THS)/(WGHT+1.)+THZ0)*0.5
C
           ZQ=FZQ2*ZT
           WGHT=AKHS*ZQ*RQVISC
           QZ0 =((WGHT*QLM+QS)/(WGHT+1.)+QZ0)*0.5
         ENDIF
C-----------------------------------------------------------------------
      ELSE                                             
C-----------------------------------------------------------------------
         ZU=Z0
C-ZEFF-ZEFF-ZEFF-ZEFF
         IF(SM.LE.0.5)THEN
           IF(ULM.EQ.0.) ULM=EPSU2
           ALPHA=ABS(ATAN(VLM/ULM)+PIHF-EPSA)
           X=ALPHA/PIFR
           ML=1+X
           ML=MIN(4,ML)
           MH=1+MOD(ML,4)
           WLOW=X-ML+1
           ZU=WLOW*ZEFF(ML)+(1.-WLOW)*ZEFF(MH)
         ENDIF
C-ZEFF-ZEFF-ZEFF-ZEFF
         UZ0=0.
         VZ0=0.
C
         ZT=Z0
         THZ0=THS
C
         ZQ=Z0
         QZ0=QS
C-----------------------------------------------------------------------
      ENDIF
C-----------------------------------------------------------------------
      ZSL=(Z(LMHK)-Z(LMHP))*0.5
C-ZEFF-ZEFF-ZEFF-ZEFF
      ZU=AMIN1(ZU,0.5*ZSL)
C-ZEFF-ZEFF-ZEFF-ZEFF
      RDZ=1./ZSL
      CXCH=EXCM*RDZ
C-----------------------------------------------------------------------
      IF(SM.GT.0.5)THEN
        DTHV=(0.608*QLM+1.)*THLM-(0.608*QZ0+1.)*THZ0
      ELSE
        DTHV=(QLM-QZ0)*CNV+THLM-THZ0
        ZT=Z0*ZTFC
      ENDIF
C
      DU2=AMAX1((ULM-UZ0)**2+(VLM-VZ0)**2,EPSU2)
C-----------------------------------------------------------------------
      RIB=BTG*DTHV*ZSL/DU2
C--------------BELJARS CORRECTION OF USTAR------------------------------
      BTGH=BTG*HPBL
      WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.)
      USTAR=AMAX1(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST)
C--------- Malagutti class experiments-------------
         Rr= USTAR*Z0*1/VISC
         SMALL=11/(M*(Rr**0.25))
         CZIV=SMALL*GLKBS
         FZU1=CZIV*VISC
         FZT2=CZIV*GRRS*TVISC*SQPR
C        WRITE(6,*)"SMALL2=",SMALL
C        WRITE(6,*)"ZU2=",ZU,"ZT2=",ZT
C--------------ZILITINKEVITCH FIX FOR ZT--------------------------------
      IF(SM.LT.0.5)ZT=EXP(ZILFC*SQRT(USTAR*Z0))*Z0  
C-----------------------------------------------------------------------
      
      IF (SM.GT.0.5.AND.RIB.GE.RIC)THEN      

          AKMS=AMAX1( VISC*RDZ,CXCH)        !! NO TURBULENCE - MAR
          AKHS=AMAX1(TVISC*RDZ,CXCH)
C-----------------------------------------------------------------------
      ELSE
C-----------------------------------------------------------------------
        ZSLU=ZSL+ZU                         !! TURBULENCE
        ZSLT=ZSL+ZT
C
        RLOGU=ALOG(ZSLU/ZU)
        RLOGT=ALOG(ZSLT/ZT)
C
        RLMO=ELFC*AKHS*DTHV/USTAR**3

        IF(SM.GT.0.5)THEN                      !! SEA POINTS FIRST

           DO ITR=1,ITRMX
C--------------1./MONIN-OBUKKHOV LENGTH-SCALE---------------------------
               ZETALT=AMAX1(ZSLT*RLMO,ZTMIN)
               RLMO=ZETALT/ZSLT
               ZETALU=ZSLU*RLMO
C
               ZETAU=ZU*RLMO
               ZETAT=ZT*RLMO
C--------------LL FUNCTIONS OVER SEA------------------------------------
               IF(RLMO.LT.0.)THEN
                 PSMZ=PSLMU(ZETAU)
                 SIMM=       PSLMU(ZETALU)-PSMZ+RLOGU
                 PSHZ=PSLHU(ZETAT)
                 SIMH=FHNEU*(PSLHU(ZETALT)-PSHZ+RLOGT)
               ELSE
                 PSMZ=PSLMS(ZETAU)
                 SIMM=       PSLMS(ZETALU)-PSMZ+RLOGU
                 PSHZ=PSLHS(ZETAT)
                 SIMH=FHNEU*(PSLHS(ZETALT)-PSHZ+RLOGT)
               ENDIF
C--------------BELJAARS CORRECTION FOR USTAR----------------------------
               USTAR=AMAX1(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST)
C--------- Malagutti class experiments-------------
               Rr= USTAR*Z0*1/VISC
               SMALL=11/(M*(Rr**0.25))
               CZIV=SMALL*GLKBS
               FZU1=CZIV*VISC
               FZT2=CZIV*GRRS*TVISC*SQPR
C              WRITE(6,*)"SMALL3=",SMALL
C-----------------------------------------------------------------------
               USTARK=USTAR*VKRM
               AKMS=AMAX1(USTARK/SIMM,CXCH)
               AKHS=AMAX1(USTARK/SIMH,CXCH)
C-----------------------------------------------------------------------
               WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.)
               RLMN=ELFC*AKHS*DTHV/USTAR**3
C-----------------------------------------------------------------------
               RLMP=RLMO
               RLMA=RLMO*WOLD+RLMN*WNEW
               RLMO=RLMA
           END DO
        
        ELSE                         !! NOW LAND POINTS ...

          DO ITR=1,ITRMX
C--------------1./MONIN-OBUKKHOV LENGTH-SCALE---------------------------
             ZETALT=AMAX1(ZSLT*RLMO,ZTMIN)
             RLMO=ZETALT/ZSLT
             ZETALU=ZSLU*RLMO
C
             ZETAU=ZU*RLMO
             ZETAT=ZT*RLMO

C--------------Businger 1981 FUNCTIONS OVER LAND W RAD. SKIN T-----------
             IF(RLMO.LT.0.)THEN
               XLU4=1.-15.*ZETALU                        !! FRAN-URIX:2011
               XLT4=1.-9.*ZETALT                         !! FRAN-URIX:2011
               XU4 =1.-15.*ZETAU                         !! FRAN-URIX:2011
               XT4 =1.-9.*ZETAT                          !! FRAN-URIX:2011
C
               XLU=SQRT(SQRT(XLU4))
               XLT=SQRT(XLT4)                            !! FRAN-URIX:2011
               XU =SQRT(SQRT(XU4))
               XT =SQRT(XT4)                             !! FRAN-URIX:2011
C
               PSMZ=PSPMU(XU)
               SIMM=PSPMU(XLU)-PSMZ+RLOGU
               PSHZ=PSPHU(XT)
               SIMH=PSPHU(XLT)-PSHZ+RLOGT
             ELSE
               ZETAU=AMIN1(ZETAU,ZTMAX)
               ZETAT=AMIN1(ZETAT,ZTMAX)
               ZETALU=AMIN1(ZETALU,ZTMAX)
               ZETALT=AMIN1(ZETALT,ZTMAX)
               PSMZ=PSPMS(ZETAU)
               SIMM=PSPMS(ZETALU)-PSMZ+RLOGU
               PSHZ=PSPHS(ZETAT)
               SIMH=PSPHS(ZETALT)-PSHZ+RLOGT
            ENDIF

C--------------BELJAARS CORRECTION FOR USTAR----------------------------
            USTAR=AMAX1(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST)
C--------------ZILITINKEVITCH FIX FOR ZT--------------------------------
            ZT=EXP(ZILFC*SQRT(USTAR*Z0))*Z0
            ZSLT=ZSL+ZT
            RLOGT=ALOG(ZSLT/ZT)
C-----------------------------------------------------------------------
            USTARK=USTAR*VKRM
            AKMS=AMAX1(USTARK/SIMM,CXCH)
            AKHS=AMAX1(USTARK/SIMH,CXCH)
C-----------------------------------------------------------------------
            WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.)
            RLMN=ELFC*AKHS*DTHV/USTAR**3
C-----------------------------------------------------------------------
            RLMP=RLMO
            RLMA=RLMO*WOLD+RLMN*WNEW
            RLMO=RLMA
         END DO      

       ENDIF              !! END OF LAND POINT PROCESSING AND SEA-LAND BRANCHIN

      ENDIF      !! END OF TURBULENCE-NO TURBULENCE BRANCHING

C--------------COUNTERGRADIENT FIX--------------------------------------
C      HV=-AKHS*DTHV
C      IF(HV.GT.0.)THEN
C        FCT=-10.*(BTG)**(-1./3.)
C        CT=FCT*(HV/(HPBL*HPBL))**(2./3.)
C      ELSE
        CT=0.
C      ENDIF
C--------------DIAGNOSTIC BLOCK-----------------------------------------
      WSTAR=SQRT(WSTAR2)/WWST
C
      UMFLX=AKMS*(ULM -UZ0 )
      VMFLX=AKMS*(VLM -VZ0 )
      HSFLX=AKHS*(THLM-THZ0)
      HLFLX=AKHS*(QLM -QZ0 )
C-----------------------------------------------------------------------
      IF(SM.GT.0.5.AND.RIB.GE.RIC)THEN
C-----------------------------------------------------------------------
Cmp        AKMS10=AMAX1( VISC/10.,CXCH)
Cmp        AKHS02=AMAX1(TVISC/02.,CXCH)
Cmp        AKHS10=AMAX1(TVISC/10.,CXCH)
        AKMS10=AMAX1( VISC/10.,4.*CXCH)
CGSM v100m
        AKMS100=AMAX1( VISC/100.,4.*CXCH)
CGSM v100m
        AKHS02=AMAX1(TVISC/02.,4.*CXCH)
        
        AKHS10=AMAX1(TVISC/10.,4.*CXCH)
CGSM v100m
        AKHS100=AMAX1(TVISC/100.,4.*CXCH)
CGSM v100m        
C-----------------------------------------------------------------------
      ELSE
C-----------------------------------------------------------------------
        ZU10=ZU+10.
CGSM v100m
        ZU100=ZU+100.
CGSM v100m
        ZT02=ZT+02.
        ZT10=ZT+10.
CGSM v100m
        ZT100=ZT+100.
CGSM v100m
C
        RLNU10=ALOG(ZU10/ZU)
CGSM v100
        RLNU100=ALOG(ZU100/ZU)
CGSM v100
        RLNT02=ALOG(ZT02/ZT)
        RLNT10=ALOG(ZT10/ZT)
CGSM v100
        RLNT100=ALOG(ZT100/ZT)
CGSM v100m
C
        ZTAU10=ZU10*RLMP
CGSM v100m
        ZTAU100=ZU100*RLMP
CGSM v100m
        ZTAT02=ZT02*RLMP
        ZTAT10=ZT10*RLMP
CGSM v100m
        ZTAT100=ZT100*RLMP
CGSM v100m

C--------------LL FUNCTIONS OVER SEA------------------------------------
        IF (SM.GT.0.5)THEN
           IF(RLMP.LT.0.)THEN                          !! CGSM v100m
             SIMM10=       PSLMU(ZTAU10)-PSMZ+RLNU10           
             SIMM100=       PSLMU(ZTAU100)-PSMZ+RLNU100
             SIMH02=FHNEU*(PSLHU(ZTAT02)-PSHZ+RLNT02)
             SIMH10=FHNEU*(PSLHU(ZTAT10)-PSHZ+RLNT10)
             SIMH100=FHNEU*(PSLHU(ZTAT100)-PSHZ+RLNT100)
           ELSE
             SIMM10=       PSLMS(ZTAU10)-PSMZ+RLNU10
             SIMM100=       PSLMS(ZTAU100)-PSMZ+RLNU100
             SIMH02=FHNEU*(PSLHS(ZTAT02)-PSHZ+RLNT02)
             SIMH10=FHNEU*(PSLHS(ZTAT10)-PSHZ+RLNT10)
             SIMH100=FHNEU*(PSLHS(ZTAT100)-PSHZ+RLNT100)
           ENDIF
C--------------Businger 1981 FUNCTIONS OVER LAND W RAD. SKIN T-----------
        ELSE
           IF(RLMP.LT.0.)THEN             !! CGSM v100                  
             XLU104=1.-15.*ZTAU10	                     !! FRAN-URIX:2011                       
             XLU1004=1.-15.*ZTAU100                      !! FRAN-URIX:2011
             XLT024=1.-9.*ZTAT02                         !! FRAN-URIX:2011
             XLT104=1.-9.*ZTAT10                         !! FRAN-URIX:2011
             XLT1004=1.-9.*ZTAT100                       !! FRAN-URIX:2011
C
             XLU10=SQRT(SQRT(XLU104))
             XLU100=SQRT(SQRT(XLU1004))
             XLT02=SQRT(XLT024)                          !! FRAN-URIX:2011
             XLT10=SQRT(XLT104)                          !! FRAN-URIX:2011
             XLT100=SQRT(XLT1004)                        !! FRAN-URIX:2011
CGSM v100
C
             SIMM10=PSPMU(XLU10)-PSMZ+RLNU10
CGSM v100	    
             SIMM100=PSPMU(XLU100)-PSMZ+RLNU100
CGSM v100
             SIMH02=PSPHU(XLT02)-PSHZ+RLNT02
             SIMH10=PSPHU(XLT10)-PSHZ+RLNT10
CGSM v100	    
             SIMH100=PSPHU(XLT100)-PSHZ+RLNT100
CGSM v100
           ELSE
             ZTAU10=AMIN1(ZTAU10,ZTMAX)
CGSM v100m	    
             ZTAU100=AMIN1(ZTAU100,ZTMAX) 
CGSM v100m
             ZTAT02=AMIN1(ZTAT02,ZTMAX)
             ZTAT10=AMIN1(ZTAT10,ZTMAX)
CGSM v100m	    
             ZTAT100=AMIN1(ZTAT100,ZTMAX)
CGSM v100m
C
             SIMM10=PSPMS(ZTAU10)-PSMZ+RLNU10
CGSM v100m	    
             SIMM100=PSPMS(ZTAU100)-PSMZ+RLNU100
CGSM v100m
             SIMH02=PSPHS(ZTAT02)-PSHZ+RLNT02
             SIMH10=PSPHS(ZTAT10)-PSHZ+RLNT10
CGSM v100m	    
             SIMH100=PSPHS(ZTAT100)-PSHZ+RLNT100
CGSM v100m
          ENDIF
C-----------------------------------------------------------------------
        ENDIF
C-----------------------------------------------------------------------
        AKMS10=AMAX1(USTARK/SIMM10,CXCH)
CGSM v100m
        AKMS100=AMAX1(USTARK/SIMM100,CXCH)
CGSM v100m
        AKHS02=AMAX1(USTARK/SIMH02,CXCH)
        AKHS10=AMAX1(USTARK/SIMH10,CXCH)
CGSM v100m	
        AKHS100=AMAX1(USTARK/SIMH100,CXCH)
CGSM v100m
C-----------------------------------------------------------------------
      ENDIF
C-----------------------------------------------------------------------
      U10 =UMFLX/AKMS10+UZ0
CGSM v100m
      U100 =UMFLX/AKMS100+UZ0
      V100 =VMFLX/AKMS100+VZ0
CGSM v100m
      V10 =VMFLX/AKMS10+VZ0
      TH02=HSFLX/AKHS02+THZ0
      TH10=HSFLX/AKHS10+THZ0
CGSM v100m      
      TH100=HSFLX/AKHS100+THZ0
CGSM v100m

C GSM  changed this section in response to problem with 2-m
C     dew point occasionally being greater than 2-m temperature
C     and similar problem at 10-m.   Now, a saturation Q is
C     calculated at each level, and the Q is constrained to
C     be no higher than the saturation value. 

      PDS=PD+PT
      TERM1=-0.068283/TLM
      PSHLTR=PDS*EXP(TERM1)
      T02=TH02*(PSHLTR*H1M5)**CAPA
      QSAT2 = PQ0/PSHLTR*EXP(A2*(T02-A3)/(T02-A4))
      Q02 =HLFLX/AKHS02+QZ0
      IF (Q02.LT.0.) THEN
        IF (QLM .GT. 0.) THEN
          Q02=QLM
        ELSE
          Q02=0.0001
        ENDIF
      ENDIF
      IF (Q02.GT.QSAT2)THEN
        Q02 = QSAT2
      ENDIF

      T10=TH10*(PSHLTR*H1M5)**CAPA
      QSAT10 = PQ0/PSHLTR*EXP(A2*(T10-A3)/(T10-A4))
      Q10 =HLFLX/AKHS10+QZ0
      IF (Q10.LT.0.) THEN
        IF (QLM .GT. 0.) THEN
          Q10=QLM
        ELSE
          Q10=0.0001
        ENDIF
      ENDIF
      IF (Q10.GT.QSAT10)THEN
        Q10 = QSAT10
      ENDIF

CGSM v100m
      P100=PDS*EXP(-100.0*G/(287.04*TLM))
      T100=TH100*(P100*H1M5)**CAPA
      QSAT100 = PQ0/P100*EXP(A2*(T100-A3)/(T100-A4))
      Q100 =HLFLX/AKHS100+QZ0
      IF (Q100.LT.0.) THEN
        IF (QLM .GT. 0.) THEN
          Q100=QLM
        ELSE
          Q100=0.0001
        ENDIF
      ENDIF
      IF (Q100.GT.QSAT100)THEN
        Q100 = QSAT100
      ENDIF
CGSM v100m

c new calculation of 10-m winds
C-----------------------------------------------------
      U10E=U10
      V10E=V10
CGSM v100m
      U100E=U100
      V100E=V100
CGSM v100m
C-----------------------------------------------------
      IF(SM.LT.0.5)  THEN
c choose the equivalent z0 here:
czj        ZU=0.01
ckm     zu=zu*0.1
C
        zuuz=amin1(zu*0.50,0.10)
        zu=amax1(zu*0.10,zuuz)
        ZU10=ZU+10.
CGSM v100m
        ZU100=ZU+100.
CGSM v100m
        RLNU10=ALOG(ZU10/ZU)
CGSM v100m 
        RLNU100=ALOG(ZU100/ZU)
CGSM v100m
        ZTAU=ZU*RLMP
        ZTAU10=ZU10*RLMP
CGSM v100m
        ZTAU100=ZU100*RLMP
CGSM v100m
c--------------------------------------------------------
        IF(RLMP.LT.0)THEN
          XLU104=1.-16.*ZTAU10
CGSM v100m
          XLU1004=1.-16.*ZTAU100
CGSM v100m
          XU104 =1.-16.*ZTAU
          XLU10=SQRT(SQRT(XLU104))
CGSM v100m	  
          XLU100=SQRT(SQRT(XLU1004))
CGSM v100m
          XU10 =SQRT(SQRT(XU104))
          SIMM10=PSPMU(XLU10)-PSPMU(XU10)+RLNU10
CGSM v100m	  
          SIMM100=PSPMU(XLU100)-PSPMU(XU10)+RLNU100
CGSM v100m
        ELSE
          ZTAU10=AMIN1(ZTAU10,ZTMAX)
CGSM v100m	  
          ZTAU100=AMIN1(ZTAU100,ZTMAX)
CGSM v100m
          SIMM10=PSPMS(ZTAU10)-PSPMS(ZTAU)+RLNU10
CGSM v100m
          SIMM100=PSPMS(ZTAU100)-PSPMS(ZTAU)+RLNU100
CGSM v100m
        ENDIF
c-----------------------------------------------------
        EKMS10=AMAX1(USTARK/SIMM10,CXCH)
CGSM v100m
        EKMS100=AMAX1(USTARK/SIMM100,CXCH)
CGSM v100m
        U10E=UMFLX/EKMS10+UZ0
        V10E=VMFLX/EKMS10+VZ0
CGSM v100m
        U100E=UMFLX/EKMS100+UZ0
        V100E=VMFLX/EKMS100+VZ0
CGSM v100m
      ENDIF
c----------------------------------------------------------------
cs.morelll     U10=U10E
cs.morelli     V10=V10E
CGSM v100m
C      U100=U100E
C      V100=V100E
CGSM v100m
C
C-----------------------------------------------------------------------
                           RETURN
                           END