PROGRAM E4MAIN 
C
C
C  COMPUTES EKMAN LAYER CURRENTS FROM INPUT PROVIDED THROUGH THE KEYBOARD.
C  Jim Price, WHOI, Sep 93. Last modified in January 1996. This can be used
C  to check the output from e4.f (SEL4)
C  Last modified in Jan 1996. 
C
      WRITE (6,72) 
   72 FORMAT (1X,//,
     C ' This program calculates Ekman layer currents using ',/,
     C ' a solution for a four layer stratifed Ekman layer. Some '
     C ,/,' data are hardwired in the FORTAN program, others are ',/,
     C ' input from the keyboard. See the listing for details. Data ',
     C /,' input from the keyboard and their default values are:',//,
     C ' QMAX is the noon maximum heat flux (630. W/m*m)',/,
     C ' TAU is the wind stress maginitude (0.07 Pa),',/,
     C ' ANG is the wind stress direction (0. = northward),',/,
     C ' RLAT is the latitude (35. deg), and ',/,
     C ' H is the depth of semi-permanent stratification (50. m)',/)
C
C
C
C  SET SOME VARIABLES TO THEIR DEFAULT VALUES (CLOSE TO THE MEAN
C   CONDITIONS FOUND AT THE LOTUS SITE IN SUMMER (SEE PRICE ET AL. (1987,
C   SCIENCE).
C
C  PQH IS THE DURATION OF HEATING (hours)
C  QMAX IS THE NOON (MAXIMUM) NET HEAT FLUX (W/m*m)
C  TAU IS THE WIND STRESS (Pa) 
C  ANG IS THE STRESS DIRECTION TOWARD, I.E., 90 = EASTWARD (deg)
C  RLAT IS THE LATITUDE (deg)
C  H IS THE DEPTH OF SEMI-PERMANENT STRATIFICATION (m)
C  ALPHA IS THE THERMAL EXPANSION COEFFICIENT (kg/m*m*m*C)
C  DTDZ IS THE STRATIFICATION JUST BELOW H (C/m).
c
      PQH = 12.
      QMAX = 630.
      TAU = 0.07
      ANG = 0.
      RLAT = 35.
      H = 50.
      ALPHA = -0.3
      DTDZ = 1./30.
C
      NPASS = 0
   60 CONTINUE
C 
      NPASS = NPASS + 1
      WRITE (6,88)  
   88 FORMAT (1X,
     C'Enter  QMAX, TAU, ANG, RLAT, H (to default',
     C ' to previous values, enter ,,,,)',/)
      READ (5,*) QMAX, TAU, ANG, RLAT, H
      IF(QMAX.LT.-1.) GO TO 999      
C
      QI = 1.2*QMAX
      QL = -0.2*QMAX
C
      IF(ABS(RLAT).LT.0.25) THEN
      RLAT = SIGN(0.25,RLAT)
      WRITE (6,3349) 
 3349 FORMAT (1X,/,1x,'WARNING, RLAT WAS RESET TO 0.25')
      END IF
C
      F = 2.*7.29E-5*SIN(RLAT/57.)
C      
      ektran = tau/(1025*f)
      write (6,888) ektran
  888 format (1x,/,'   The ekman transport is ', f9.2,/)
C  
      PQ = PQH*3.6E3
C
      IF(NPASS.EQ.1) WRITE (6,70)
   70 FORMAT (1X,//,
     C ' Now, get ready for a lot of exciting results:',/,
     C ' (U,V) is the Ekman layer current (cm/s).')
C
      DO 90 K=1,20
      ZK = FLOAT(K-1)*3.
C
      CALL E4(QMAX,PQ,ALPHA,TAU,ANG,F,H,DTDZ,ZK,U,V)
C
C  AC IS THE ANGLE TO THE RIGHT (F > 0) OF THE WIND (DEGREES)
C
      AC =  90. - (180./3.1415)*ATAN2(V,U)     
      IF(K.EQ.1) THEN
      WRITE (6,803) QMAX,TAU,ANG,RLAT,H
  803 FORMAT (1X,/,1X,' QMAX, TAU, ANG, RLAT and H are: ',F6.0, F6.2,
     C F6.0, 2X, F6.2,  F6.0)
      WRITE (6,801)
  801 FORMAT (1X,/,'     Depth(m)  U      V(cm/s)  ANG') 
      END IF
C
C  CONVERT TO CM/SEC
C 
      CMSEC = 100.
      U = CMSEC*U
      V = CMSEC*V
C
      WRITE (6,800) ZK,U,V,AC
  800 FORMAT (1X, F8.0,2X,3(2F7.1,4X))
   90 CONTINUE
      WRITE (6,92)
   92 FORMAT (1X,//)
C
      GO TO 60
C
  999 CONTINUE
c
      STOP
      END
C