C
      SUBROUTINE E4(Q,PQ,ALPHA,TAU,ANG,F,H,DTDZ,Z,U,V)
C
C  COMPUTES THE STRATIFIED EKMAN LAYER CURRENT AT DEPTH Z.
C  THIS MODEL IS MOST APPROPRIATE IF THERE IS SIGNIFICANT DIURNAL CYCLING
C  AT THE SITE. IF THERE IS NOT, THEN THE RESULT IS ESSENTIALLY THAT FOR
C  A HOMOGENEOUS (SLAB) LAYER OF DEPTH H (PRETTY DULL).
C
C  THE MODEL REQUIRES THE FOLLOWING LONG LIST OF PARAMETERS THAT 
C  CHARACTERIZE THE SITE AND THE SURFACE FLUXES: 
C            
C  INPUT:
C  Q IS THE NOON (MAXIMUM) NET SURFACE HEAT FLUX (W/m*m),
C  PQ IS THE TIME INTERVAL DURING WHICH THE NET HEAT FLUX IS POSITIVE (s,)
C  ALPHA IS THE THERMAL EXPANSION COEFFICIENT, kg/(m*m*m*C), 
C  TAU IS THE WIND STRESS MAGNITUDE (Pa),
C  ANG IS THE WIND STRESS DIRECTION TOWARD (deg),
C  F IS THE CORIOLIS PARAMETER (1/s),
C  H IS THE DEPTH OF THE SEMI-PERMANENT STRATIFICATION (m),
C  DTDZ IS THE TEMPERATURE VERTICAL GRADIENT JUST BELOW H (C/m),
C  Z IS THE SINGLE DEPTH WHERE THE SOLUTION IS TO BE EVALUATED (m),
C        
C  OUTPUT:
C  U AND V ARE THE CURRENT COMPONENTS IN THE SAME FRAME THE WIND STRESS (cm/s).
C
C  Jim Price, January 1996.  THIS IS WORK IN PROGRESS.
C
C  SOME COMMENTS:
C  DENSITY IS ASSUMED TO DEPEND UPON TEMPERATURE ONLY (S IS IGNORED). 
C  THE LIMIT OF VANISHING Q IS FINE BUT Q = 0 IS AVOIDED.  THE LIMIT
C  OF VANISHING TAU OR F ARE SIMILAR.  
C
C  COMPARED TO THE TWO LAYER MODEL OF 1993, THIS VERSION ADDS A THIRD
C  'TRANSITION' LAYER BETWEEN THE DIURNAL LAYER AND THE NIGHTTIME MIXED LAYER.
C  THEN IT ADDS A FOURTH TRANSITION LAYER AT THE TOP OF THE MAIN THERMOCLINE.
C
C  SOME FURTHER EMBELLISHMENT COULD INCLUDE USING THE FULL PWP FORMULA
C  FOR TRAPPING DEPTH (I.E., INCLUDE OPTICAL ABSORPTION) AND ADD SALINITY.
C
C  THE BACKGROUND FOR THIS CALCULATION IS STARTED IN PRICE ET AL., (1986, JGR).
C
      G = 9.8
      DW = 1023.
      CPW = 4000.
      U = 0.
      V = 0.
C
C  CHECK SOME VARIABLES AND RESET IF NEED BE TO AVOID MATH ERRORS.
C  
      IF(Q.LT.0.01) Q = 0.01
      IF(ABS(F).LT.1.E-9) F = SIGN(1.E-9,F)
      IF(TAU.LT.0.00001) TAU = 0.00001
      DTDZ = ABS(DTDZ)
      IF(DTDZ.LT.0.001) DTDZ = 0.001
      IF(ABS(ALPHA).LT.0.001) ALPHA = 0.001
C
C
C  THE CALCULATION IS DONE AS IF F > 0, AND THE RESULT IS ROTATED AT THE END IF
C  THIS IS A SOUTHERN HEMISPHERE CASE.
C
C  NOTE THAT PQ IN THE ORIGINAL PWP FORMULAE WAS THE HALF THE INTERVAL DURING
C  WHICH THE NET HEAT FLUX IS POSITIVE. HERE IT IS THE FULL INTERVAL, AND 
C  THUS THE FACTOR OF 2 THAT APPEARS BELOW IN THE VSTAR, ETC.
C
      USTAR = SQRT(TAU/DW)
      QSTAR = G*ABS(ALPHA)*Q/(CPW*(DW**2))
      VSTAR = SQRT(QSTAR*PQ/2.)
C
      PHAS = ABS(F)*PQ/2.
      IF(PHAS.GT.3.1415)  PHAS = 3.1415
      PTAU = (1./ABS(F))*SQRT(2. - 2.*COS(PHAS)) 
      TRAPD = (PTAU*USTAR**2)/SQRT(QSTAR*PQ/2.)
C
C      WRITE (5,887) Q, PQ, ALPHA, TAU
C  887 FORMAT (1X,'Q, PQ, ALPHA, TAU', 6E10.3)
c      WRITE (6,888) USTAR, QSTAR, VSTAR, PTAU, TRAPD
c  888 FORMAT (1X,'USTAR, QSTAR, VSTAR,PTAU, TRAPD', 8E10.3)
C
C  MAKE SURE THAT THE TRAPPING DEPTH IS LESS THAN H 
C
      IF(TRAPD.GT.H) TRAPD = H 
C      
C  AND THE DIURNAL CYCLING EFFECT WILL BE OF NO CONSEQUENCE FOR THE EKMAN LAYER.      
C
C  NOW, EVALUATE THE TWO-LAYER STRATIFIED EKMAN LAYER MODEL. 
C
      PHAS = ABS(F)*PQ
      UEK = TAU/(ABS(F)*DW*H)
      UEK1 = TAU/(ABS(F)*DW*TRAPD)
      UDIFF = UEK1*(PQ/8.64E4)*(1. - SIN(PHAS)/PHAS)
      VDIFF = UEK1*(PQ/8.64E4)*(1. - COS(PHAS))/PHAS
C
      V1 = VDIFF*(H - TRAPD)/H
      V2 = V1 - VDIFF
      U1 = (UEK*H + UDIFF*(H - TRAPD))/H
      U2 = U1 - UDIFF
C
C EVALUATE THE CURRENT AT THE DEPTH Z. 
C
      IF(Z.GT.TRAPD) THEN
        U = U2
        V = V2
      END IF  
C      
      IF(Z.LT.TRAPD) THEN
        U = U1
        V = V1
      END IF    
C
C  THE THIRD LAYER EMBELLISHMENT OF THE TWO LAYER MODEL, IN WHICH THE
C  VELOCITY IS INTERPOLATED BETWEEN THE DIURNAL AND NIGHTTIME LAYERS.
C  THE TRANSITION LAYER THICKNESS, TTRAN, IS ASSUMED EQUAL TO TRAPD. 
C
C  FIRST, CHECK TO INSURE THAT THE DEEPEST EXTENT OF THE TRANSITION LAYER  
C  WILL NOT EXCEED H. IF IT DOES, THEN REDUCE THE THICKNESS OF THE
C  TRANSITION LAYER SO THAT IT DOES NOT.
C
      TTRAN = TRAPD
      IF((1.5*TRAPD).GT.H) THEN
        TTRAN = 2.*(H - TRAPD)
      END IF
C
      I3L = 1
      IF(I3L.EQ.1) THEN
        IF(Z.GT.(0.5*TTRAN).AND.Z.LT.(1.5*TTRAN)) THEN 
          ZL = 1.5*TTRAN - Z
          ZU = Z - 0.5*TTRAN
          U = (U1*ZL + U2*ZU)/TTRAN
          V = (V1*ZL + V2*ZU)/TTRAN
        END IF
      END IF     
C
C  THE FOURTH LAYER EMBELLISHMENT IS HERE. YIKES. WHEN WILL THIS EVER END!!!
C  THIS PART IS PROVISIONAL, FOR SURE.
C
      I4L = 1
      IF(I4L.EQ.1) THEN
C
C  COMPUTE THE THICKNESS OF THE FOURTH LAYER, D4.        
C
      SPD2 = 2.*SQRT(U2**2 + V2**2) 
      D4 = SPD2*SQRT(DW/(2.*G*ABS(ALPHA*DTDZ)))
C    
      IF(Z.GT.(H + D4/2.)) THEN
        U = 0.
        V = 0.
        RETURN
      END IF
      IF(Z.GT.(H - D4/2.)) THEN
        DEPL4 = Z - (H - D4/2.)      
        U = U2*(D4 - DEPL4)/D4
        V = V2*(D4 - DEPL4)/D4
      END IF 
C
      END IF
C      
C      
      IF(F.LT.0.) CALL ROTE2(U,V,180.)
C      
      RETURN
      END
C
      SUBROUTINE ROTE2(U,V,ANG)
C
C  ROTATE THROUGH THE ANGLE ANG (DEG)
C
      CA = COS(ANG*3.1415/180.)
      SA = SIN(ANG*3.1415/180.)
      U1 = U
      V1 = V
      U = U1*CA + V1*SA
      V = V1*CA - U1*SA
      RETURN
      END