PROGRAM MAIN
C  This MAIN program is to demonstrate the workings of SUBROUTINE FLUX;
C  only the latter would be included in a user created subroutine to
C  supply surface forcing data.
      DO N=1,3
      write(6,'(/,''        DU        DV        DT        DQ      TAUX
     &      TAUY       SHF       SVF      CD        CH     '')')
      DT=-20+10.*FLOAT(N)
      do i=1,26,5
      DU=FLOAT(i-1)
      DU=MAX(DU,1.0)
      DV=0.0 
      DV=MAX(DV,1.0)
      DQ=0.
      CALL FLUX(CU,CT,DU,DV,DT,DQ,TAUX,TAUY,SHF,SVF)
      write(6,'(1x,8f10.4)') DU,DV,DT,DQ,TAUX,TAUY,SHF,SVF
      CD=CU*CU
      CH=CU*CT
      write(6,'(1x,8(10x),2f10.4)') CD,CH
      enddo
      ENDDO
      stop
      end

      SUBROUTINE FLUX(CU,CT,DU,DV,DT,DQ,TAUX,TAUY,SHF,SVF)
C
C      Subroutine to calculate ABL surface fluxes.
C
C      CU, CT = inverse law of the wall variables.
C        (Routine should benefit from a reasonable first guess)
C      DU, DV = velocity components relative to surface
C      TAUX,TAUY = momentum fluxes/air density
C      DT, DQ = temperature and specific humidity relative to surface
C      SHF, SVF = Sensible Heat/(air density*specific heat)
C           and Vapor flux/air density.
C      ITER = Number of iterations.
C
      REAL kappa,LINV
      DATA gee/9.807/
      DATA kappa/0.41/,BETAT/.00357/,BETAQ/0.608/
      DATA VISC/15.E-6/
      DATA PRT/0.9/
      ALPHA=0.0144
      ITER=4
C  Begin interation
      CU=.03
      CT=.03
      DO I=1,ITER
      UST=ABS(CU)*SQRT(DU**2+DV**2+.0001)
      SHF=CT*UST*DT
      SVF=CT*UST*DQ
       LINV=kappa*gee*(BETAT*SHF+BETAQ*SVF)/UST**3
      ZETA=10.*LINV
      CALL ZETAPSI(ZETA,PSI)
      Z0=MAX(ALPHA *UST**2/gee,0.14*VISC/UST)
      YKF=3.14*sqrt(ust*z0/VISC)*.60+2.11
      CU=kappa/(log(10/Z0)+PSI)
      CU=MAX(CU,0.)
      CT=kappa/(log(10/Z0)+kappa*YKF/PRT+PSI)/PRT
      CT=MAX(CT,0.)
C     write(6,'(1x,8F10.4)')
C    1    DU,DV,DT,DQ,ZETA,PSI,CU,CT
      ENDDO
C  END iteration
      UST=CU*SQRT(DU**2+DV**2)
      TAUX=CU*UST*DU
      TAUY=CU*UST*DV
      SHF=CT*UST*DT
      SVF=CT*UST*DQ
      Return
      end
c
      SUBROUTINE ZETAPSI( ZETA,PSI)
      DIMENSION PS(41)
      DATA PS/
     1 -2.0095,  -1.9915,  -1.9732,  -1.9544,  -1.9352,  -1.9155,
     2 -1.8953,  -1.8747,  -1.8534,  -1.8316,  -1.8092,  -1.7862,
     3 -1.7624,  -1.7380,  -1.7127,  -1.6867,  -1.6597,  -1.6318,
     4 -1.6029,  -1.5728,  -1.5416,  -1.5091,  -1.4751,  -1.4395,
     5 -1.4023,  -1.3631,  -1.3218,  -1.2781,  -1.2316,  -1.1821,
     6 -1.1291,  -1.0718,  -1.0097,  -0.9416,  -0.8663,  -0.7818,
     7 -0.6853,  -0.5723,  -0.4352,  -0.2583,   0.0000/
      I=1+(ZETA+4.)/.1
      IF(I.LT.1) PSI=-2.0095
      IF(I.GE.1.and.I.LT.41)
     1       PSI=PS(I)+(PS(I+1)-PS(I))*(ZETA+4.-.1*(I-1))/.1
      IF(I.GE.41) PSI=5.*ZETA
      RETURN
      END