
      PROGRAM MYMAIN
C
C  This main program is used to drive a Mellor ocean model 
C  through a diurnal cycle using an idealized, prescribed 
C  daily cycle of heating, or, an arbitrary (observed) time
C  series of surface flux data. 
C
C  The program assumes the following logical units:
C   5 is the terminal for input,
C   6 is the terminal for output,
C   8 may be opened to read air/sea flux data (optional),
C   9 may be opened to read the initial T, S profile (optional),
C  10 may be opened to store matrix data (disabled),
C  11 is opened to store lists from each time step, 
C  13 is opened to store profile data if requested.
C
C
C  The original code was written by Patrice Klein. The main program was
C  heavily revised by Jim Price in March, 1991 who is responsible for 
C  any and all errors.  This code may not be a standard 'Mellor' model,
C  so be very wary of your results.  Please direct any
C  comments or questions to Jim Price, WHOI, Woods Hole, MA
C  02543, USA, tel. 508-548-1400, x2526.
C
C
C
      CHARACTER*1 ANY
      CHARACTER*16 IFILE
C
      PARAMETER (NSOBS = 399)
C
C  The arrays below are used to store the initial Z,T,S.
C
      DIMENSION ZO(NSOBS), TO(NSOBS), SO(NSOBS), DO(NSOBS) 
C
C  The arrays below store times series of surface flux data.
C
      DIMENSION DAYA(NSOBS), QIA(NSOBS), QLA(NSOBS), TXA(NSOBS),
     X TYA(NSOBS), EMPA(NSOBS) 
C
C
C  Where arrays have had to be dimensioned here or in subroutines, 
C  they are set to 399 so that they can be easily found.
C
      PARAMETER (KB=150)
      DIMENSION ZP(KB),DZP(KB)
      COMMON/BLKT/TF(KB),T(KB),TB(KB),SF(KB),S(KB),SB(KB)
      COMMON/BLKU/UF(KB),U(KB),UB(KB),VF(KB),V(KB),VB(KB)
      COMMON/BLKQ/Q2F(KB),Q2(KB),Q2B(KB)
      COMMON/BLKR/RHO(KB),WT(KB)
      COMMON/BLKWU/WUSURF,WVSURF,WUBOT,WVBOT
      COMMON/BLKWT/WTSURF,WSSURF,PROFIL(KB)
      COMMON/BLKZ/Z(KB),ZZ(KB),DZ(KB),DZZ(KB),UMOL
      COMMON/BLKH/H,DT
      COMMON/BLKA/A(KB),C(KB),VH(KB),VHP(KB)
      COMMON/BLKK/KM(KB),KH(KB),KQ(KB),L(KB)
      COMMON/INDS/KBM1,KBM2
      COMMON/MODEL/GM(KB),GH(KB),SM(KB),SH(KB),KN(KB),
     1     SPROD(KB),BPROD(KB),PROD(KB),DTEF(KB)
      COMMON/CSTUR/A1,B1,A2,B2,C1,D1,D2,A1A16,A1A112,A1A29,A1A26,
     1A2B23,A1A212
      REAL KAPPA,KM,KH,KQ,L
      DIMENSION UU(KB),VV(KB)
C
c      PARAMETER (NT1=24)
c
c      DIMENSION UX(NT1,KB),VX(NT1,KB),TTX(NT1,KB)
c      DIMENSION DX(NT1,KB),EZ1(KB)
C
      DIMENSION TOX(10000),TOY(10000),RA(10000),HSURF(10000)
      DIMENSION TSEX(10000)
c
c
      DATA KAPPA/0.40/,EPS/1.E-10/,RPI/3.14159/
C
C  NZMAX should match the size of the initial profile arrays 
C  dimensioned above, and NMET should match the flux arrays.
C
      NZMAX = 399
      NMET = 399
C
      ISTART = 1
      IEND = 5000
C
C
        OPEN(10,FILE='MYMAT.DAT',STATUS='UNKNOWN',
     -      FORM='UNFORMATTED')
        OPEN(11,FILE='MYLIS.DAT',STATUS='UNKNOWN',
     -      FORM='FORMATTED')
        OPEN(13,FILE='MYPROF.DAT',STATUS='UNKNOWN',
     -      FORM='FORMATTED')

      BILEC1=0.
      BILEC2=0.
      BILECT=0.
      ECB=0.
C
c
      GRAV=9.806
      SMOTH=0.10
      UMOL=1.34E-5
C
      A1 = 0.78
      B1 = 15.0
      A2 = 0.79
      B2 = 8.0
      C1 = 0.056
C
      D1=A1*(1.-3.*C1)
      D2=A2
      A1A16=A1*A1*6.
      A1A112=A1*A1*12.
      A1A29=A1*A2*9.
      A1A26=A1*A2*6.
      A2B23=A2*B2*3.
      A1A212=A1*A2*12.

      WRITE (6,800)
  800 FORMAT (
     X 1X, 'This program provides a means to run the Mellor',/,
     X 1X, 'turbulence closure model. To run the model you must',/,
     X 1X, 'prescribe air/sea fluxes, the initial',/,
     X 1X, 'temperature and salinity profile, and several',/,
     X 1X, 'paprameters defining the site (i.e., latitude).',/,
     X 1X, 'Some of these (Ug,Vg) are relevent to the bottom',/,
     X 1X, 'boundary layer which this model includes, and are',/,
     X 1X, 'different from the PWP model. So watch out!',/,
     X 1X, 'The site and model parameters are as follows',/,
     X 1X, '(default values and units given in parens.):',/)
      WRITE (6,802)
  802 FORMAT (
     X 1X, 'RLAT,  latitude (31.)',/,
     X 1X, 'BETA1, longwave extinction coefficient (0.6 m)',/,
     X 1X, 'BETA2, shortwave extinction coefficient (20. m)',/,
     X 1X, 'Ug, Vg, the imposed geostrophic current (0., 0. m/sec)',/,

     X 1X,/,
     X 1X, 'DT,  time step (400. sec)',/,
     X 1X, 'DAYS, the number of days to run (1. day)',/,
     X 1X, 'H, the water column depth (80. m)',/,
     X 1X,/)
C
C  Initialize RLAT, the latitude.
C
      RLAT = 31.
C
C  Initialize the solar radiation absorption;
C  BETA1 is the  non-penetrating part of the insolation, 
C  BETA2 is the penetrating (shortwave) part.
C  Values below are appropriate for fairly clear waters.
C
      BETA1 = 0.60
      BETA2 = 20.
C
      UG = 0.
      VG = 0.
C
      WRITE (6,888)
  888 FORMAT (1X,'ENTER RLAT, BETA1, BETA2, Ug, Vg')
      READ (5,*) RLAT, BETA1, BETA2, UG, VG
C
C
      XM=.62
      ETA1=1./BETA1
      ETA2=1./BETA2
C
      COR = 2.*7.292E-5*SIN(RLAT*3.141/180.)
C
      PGX = COR*VG
      PGY = -COR*UG
C
C  Define the temporal and vertical resolution. In this model, the vertical
c  resolution (dz) is given by the water column depth, H, divided by the number
c  of grid cells, KB.
c
      DT1 = 400.
      DAYS = 1.
      H = 80.
C
      WRITE (6,706) 
  706 FORMAT (1X,/,1X,'ENTER DT, DAYS, H')
      READ (5,*) DT1, DAYS, H
C
      DTD = DT1/8.64E4
C
C
      KBM1=KB-1
      KBM2=KB-2
      CALL DEPTH(Z,ZZ,DZ,DZZ,KB,KBM1)
      TIME=0.
      IDT1=INT(DT1)
      DT2=2.*DT1
C
C  Initialize the Z, T, S, U, V profiles.
C
C  Read in the data used to define the initial T,S profiles
C  from logical unit NUNIT (read in below). The data
C  is presumed to be free format, with each record being a 
C  triple (Z, T, S).  The first record should be the surface
C  value, i.e., 0., 20., 36. is a reasonable surface value, 
C  and the last true value should be deep enough to avoid
C  entering into the caluclation (unless intentional).
C  The end of file is defined either by the end
C  of the data (if read from a disk file), or by a z value less
C  than -10. These data are then linearly interpolated to give
C  profiles at DZ resolution needed by the model.
C
      WRITE (6,458)
  458 FORMAT (1X,/,
     X 1X,'The temperature and salinity  profile can',/,
     X 1X,'be defaulted to a profile that is hardwired into the',/,
     X 1X,'program (ENTER 1, the default), or enter Z, T, S data'/,
     X 1X,'from the terminal (ENTER 2), or read the data from a',/,
     X 1X,'disk file (ENTER 3).')
C
      MZTS = 1
      READ (5,*) MZTS
C
C
      IF(MZTS.EQ.1) THEN
      NTS = 3
      ZO(1) = 0.
      ZO(2) = 20.
      ZO(3) = 200.
      TO(1) = 20.
      TO(2) = 20.
      TO(3) = 5.
      SO(1) = 36.
      SO(2) = 36.
      SO(3) = 36.
      END IF
      IF(MZTS.EQ.1) GO TO 53
C
      IF(MZTS.EQ.3) THEN
      WRITE (6,3290)
 3290 FORMAT (1X,/,1X,'ENTER THE NAME OF THE Z,T,S DATA FILE')
      READ (5,116) IFILE
      OPEN (UNIT=9,FILE=IFILE)
      END IF
C
C
      DO 50 J=1,NZMAX
      IF(MZTS.EQ.2) THEN
      WRITE (6,876)
  876 FORMAT (1X,'ENTER A Z,T,S (END = Z < -10)')
      READ (5,*) ZO(J), TO(J), SO(J)
      END IF
C
      IF(MZTS.EQ.3) THEN
      READ(9,*,END=53) ZO(J), TO(J), SO(J)
      END IF
C
      IF(ZO(J).LT.-9.) GO TO 53
      NTS = NTS + 1
C     IF(MZTS.EQ.3) WRITE (6,503) ZO(J),TO(J),SO(J)
  503 FORMAT (1X,'INITIAL Z, T, S ARE', 3F10.2)
   50 CONTINUE
   53 CONTINUE
C
C  Interpolate the initial data to get a profile at DZ resolution.
C  The arrays ZP and DZP are first computed to find the grid levels.
C
      DO 52 J=1,KB
      ZP(J) = -1.*H*ZZ(J) 
      DZP(J) = -DZZ(J)*H
      CALL INTX(ZO,TO,NSOBS,NTS,ZP(J),TB(J))
      CALL INTX(ZO,SO,NSOBS,NTS,ZP(J),SB(J))
c      WRITE (6,54) J,ZP(J),ZZ(J),TB(J),SB(J),RHO(J)
   54 FORMAT (1X,'DENSITY PROF.',I4,2X,7F10.3)
   52 CONTINUE
C

      DO 5 K=1,KB
      UB(K)=UG
      U(K)=UB(K)
      UF(K)=U(K)
      VB(K)=VG
      V(K)=VB(K)
      VF(K)=V(K)
      WT(K)=0.
      T(K)=TB(K)
      TF(K)=T(K)
      S(K)=SB(K)
      SF(K)=S(K)
C
      Q2B(K)=1.E-5
      Q2(K)=Q2B(K)
      L(K)=1.0
      KM(K)=0.
      KH(K)=0.
      KQ(K)=0
    5 CONTINUE
C
      CALL DENS(T,S,RHO)
C
C
C  Initialization of Z, T, S, U, V is finished.
C
      DVT = 0.
      VVT = 10.
C
C  Set up to input air/sea flux data.
C
      WRITE (6,457)
  457 FORMAT (1X,//,
     X 1X,'You can prescribe the air/sea flux data by',/,
     X 1X,'specifying the parameters of an idealized diurnal',/,
     X 1X,'cycle (ENTER 1 (default)), or, read the air/sea',/,
     X 1X,'flux data from a disk file (ENTER 2).')
C
      METU = 1
      READ (5,*) METU
C
      IF(METU.EQ.1) THEN
C
      WRITE (6,2790)
 2790 FORMAT (1X,/,1X,'In order to prescribe an idealized',/,
     X 1X, 'diurnal cycle, the program will now ',/,
     X 1X, 'query for the values of the following variables, all ',/,
     X 1X, 'of which may be defaulted to values in parens.',/,
     X 1X, 'The default values are chosen to match day 131',/,
     X 1X, 'of the PWP JGR paper (91) C7, 8411-8427, 1986.',/)
      WRITE (6,801)
  801 FORMAT (
     X 1X, 'QIMAX, noon amplitude of insolation (978. W/m**2)',/,
     X 1X, 'QL, steady heat loss (-126. W/m**2)',/,
     X 1X, 'TX, steady east wind stress (0.07 Pa)',/,
     X 1X, 'EMP, evaporation minus precipitation (0. m/sec)',/,
     X 1X, 'PQFAC,  duration of insolation (10. hours)',/)
C
C  Initialize the amplitude of solar insolation (QI), the
C  heat loss (QL), and the wind stress (TX, uses east only).
C
      QIMAX = 978.
      QL = -126.
      TX = 0.07
      EMP = 0.
C
C  Initialize PQFAC, the duration of daylight in hours.
C
      PQFAC = 10.
C
C  
      WRITE (6,777)
  777 FORMAT (1X,'ENTER QIMAX, QL, TX, EMP, PQFAC')
      READ (5,*) QIMAX, QL, TX, EMP, PQFAC
C
C  Convert PQFAC to a non-d form.
C
      IF(PQFAC.LT.0.001) PQFAC = 0.001
      PQFAC = PQFAC/12.
      END IF 
C
C  Come here if you intend to read a time series of flux data from
C  a disk file.
C
      IF(METU.EQ.2) THEN
C
      WRITE (6,117) 
  117 FORMAT (1X,/,1X,'ENTER THE NAME OF THE THE AIR/SEA FLUX FILE')
      READ (5,116) IFILE
  116 FORMAT (A16)
      OPEN (UNIT=8,FILE=IFILE)
C
C  In the file read below, the array EMPA is the integrated evaporation
C  minus precipitation (m), not the rate as above. This is done to insure
C  that the net E-P is correctly applied to the model regardless of the
C  time step.
C
      DO 233 J=1,NMET
      READ (8,*,END=459) DAYA(J), QIA(J), QLA(J), TXA(J), TYA(J),
     X EMPA(J)
C
      IF(DAYA(J).LT.-9.) GO TO 459
C
      NMETO = NMETO + 1
C
      IF(NMETO.LT.25) THEN
      WRITE (6,8838) NMETO, DAYA(J), QIA(J), QLA(J), TXA(J), TYA(J),
     X EMPA(J)
 8838 FORMAT (1X,'FIRST 25 FLUX DATA', I4, F7.2,2F7.0,2F7.2,E10.3)
      END IF
C
C
  233 CONTINUE
  459 CONTINUE
C
C  Set NMET equal to the actual number of met observations input
C
      NMET = NMETO
C
      WRITE (6,8833)
 8833 FORMAT (1X,/,1X,'ENTER THE STARTING TIME (DAY)')
      READ (5,*) DAY0
C
      END IF
C
C  Finished with specifying air/sea flux data.
C
C
      WSSURF=0.
      WUBOT=0.
      WVBOT=0.
C
      ROCP=4.2E+6
      DO 9 I=1,KBM1
      PROFIL(I)=(XM*(EXP(ETA1*Z(I)*H)-EXP(ETA1*Z(I+1)*H))
     1+(1.-XM)*(EXP(ETA2*Z(I)*H)-EXP(ETA2*Z(I+1)*H)))/DZ(I)/H/ROCP
  9   CONTINUE
C
C
      WRITE (6,8899)
 8899 FORMAT (1X,/,1X,'DO YOU WANT A SCREEN LISTING OF DATA ?',/,
     X 1X,'(0 = NO, N = EVERY Nth STEP (DEFAULT = 6)')
      IWRT = 6
      READ (5,*) IWRT
C
C
C  Section below stores some data onto a disk file for later
C  analysis or plotting.
C
      ISTOR = 1 
      WRITE (6,3186)
 3186 FORMAT (1X,/,1X,'DO YOU WANT TO STORE THE DATA ON DISK ?',/,1X,
     X '(0 = NO, 1 = YES (DEFAULT)')
      READ (5,*) ISTOR
C
C
      ITFREQ = 20
      IZFREQ = 3
C
      IF(ISTOR.EQ.1) THEN
      WRITE (6,3187) 
 3187 FORMAT (1X,/,1X,'ENTER THE DECIMATION RATE FOR TIME AND DEPTH',
     X /,1X,'(DEFAULT IS 20 AND 3)')
      READ (5,*) ITFREQ, IZFREQ
      END IF
C
C
C*********************************************************************
C      BEGIN TIME MARCH
C*********************************************************************
C
       DO 9000 ITER=ISTART,IEND
C
C
C
      TIME=TIME+DT1/3600.
C
C  Compute the solar insolation as a function of the time,
C  where TIMED is time since start in days.
C
C
      TIMED = TIME/24. 
      DTD = DT1/8.64E4
      IF(TIMED.GT.DAYS) GO TO 499
C
      IF(METU.EQ.1) THEN
      P2 = 3.141*2.
      TIMS = AMOD(TIMED,1.)
      TIMS = TIMS - 0.50
      QI = 0.
      IF(COS(P2*TIMS).GT.0.) THEN
      QI = QIMAX*COS(TIMS*P2/PQFAC)
      IF(QI.LT.0.) QI = 0.
      END IF
      EMP = 0.
      END IF
C
      IF(METU.EQ.2) THEN
C
      DAY = TIMED + DAY0
      CALL INTX (DAYA,QIA,NSOBS,NMET,DAY,QI)
      CALL INTX (DAYA,QLA,NSOBS,NMET,DAY,QL)
      CALL INTX (DAYA,TXA,NSOBS,NMET,DAY,TX)
      CALL INTX (DAYA,TYA,NSOBS,NMET,DAY,TY)
C
C  Compute the E-P rate from the integrated E-P in EMPA.
C 
      DAYP = DAY + DTD
      CALL INTX (DAYA,EMPA,NSOBS,NMET,DAYP,EMP2)
      CALL INTX (DAYA,EMPA,NSOBS,NMET,DAY,EMP1)
      EMP = (EMP2 - EMP1)/(DTD*8.64E4)
C
      END IF
C
C
      TOX(ITER) = TX/1023.
      TOY(ITER) = TY/1023.
      RA(ITER) = QI
      HSURF(ITER) = QL
c      WRITE (6,3344) ITER, TIMS, RA(ITER), HSURF(ITER), TOX(ITER),
c     c   TOY(ITER)
 3344 FORMAT (1X,'FLUX DATA' I6, 5F10.5)
c
c 
      WUSURF=TOX(ITER)
      WVSURF=TOY(ITER)
      WWSURF=(TOX(ITER)**2.+TOY(ITER)**2.)**.5
      WTSURF=HSURF(ITER)/ROCP
      WSSURF=EMP*S(1)
C
      CALL PROFQ(DT2)
C
      DO 325 K=1,KB
      Q2(K)=Q2(K)+.5*SMOTH*(Q2F(K)+Q2B(K)-2.*Q2(K))
      Q2B(K)=Q2(K)
 325  Q2(K)=Q2F(K)
C
C
C
      CALL PROFT(TF,TB,WTSURF,DT2,RA(ITER))
      CALL PROFT(SF,SB,WSSURF,DT2,0.)
C
      DO 345 k=1,KB
      T(K)=T(K)+.5*SMOTH*(TF(K)+TB(K)-2.*T(K))
      TB(K)=T(K)
      T(K)=TF(K)
      S(K)=S(K)+.5*SMOTH*(SF(K)+SB(K)-2.*S(K))
      SB(K)=S(K)
      S(K)=SF(K)
  345 CONTINUE
C
C
      CALL DENS(T,S,RHO)
C
      DO 380 K=1,KBM1
      UU(K) = UB(K) + DT2*(COR*V(K) - PGX)
  380 VV(K) = VB(K) + DT2*(-1.*COR*U(K) - PGY)
C
C
      CBC=AMAX1(.0025,.16/ALOG((ZZ(KBM1)-Z(KB))*H/.01)**2)
      CDB = CBC
      CBC=CBC*SQRT(UB(KBM1)**2+VB(KBM1)**2)
      CDB2 = CBC
C
C  TO DELETE BOTTOM B.L. SET CBC = 0. AS BELOW                                
C      CBC = 0.
C
      CALL PROFC(UU,UF,WUSURF,CBC,DT2)
      WUBOT1 = WUBOT
      WVBOT1 = WVBOT
      CALL PROFC(VV,VF,WVSURF,CBC,DT2)
      WUBOT2 = WUBOT
      WVBOT2 = WVBOT
C
C STORE BOTTOM STRESS
C 
      WUBOT = WUBOT1
      WVBOT = WUBOT2
C
C
C      WRITE (11,7733) CDB,CDB2,WUBOT1,WVBOT1,WUBOT2,WVBOT2
C 7733 FORMAT (1X,'CDS, ETC', 6E10.3)
C
      DO 382 K=1,KB
      U(K)=U(K)+.5*SMOTH*(UF(K)+UB(K)-2.*U(K))
      V(K)=V(K)+.5*SMOTH*(VF(K)+VB(K)-2.*V(K))
      UB(K)=U(K)
      U(K)=UF(K)
      VB(K)=V(K)
  382 V(K)=VF(K)
C
C  AS A TEST, EXRAPOLATE ALL VARIABLES TO THE BOTTOM
C  (SEEMS FINE)
C
      T(KB) = T(KBM1)
      S(KB) = S(KBM1)
      Q2(KB) = Q2(KBM1)
      U(KB) = U(KBM1)
      V(KB) = V(KBM1)
C
C  CHECK MOMENTUM BUDGET FOR BOTTOM LAYER IN CASE WHERE F = 0.
C  (WORKS FINE)
C 
C      KBB = KB - 15
C      UBSUM = 0.
C      DO 4445 J=KBB,KB
C      UBSUM = UBSUM + (UG - U(J))*DZP(J)
C 4445 CONTINUE
C
C      TXBSUM = TXBSUM + DT1*WUBOT
C
C  Evaluate the energy budget. Seems to make sense only with surface
C  stress (for now).
C
C  ecf is the kinetic energy
C  work is the work by the wind stress on the surface current
C  shrp is shear production
C  ecb is the kinetic energy at the previous time step
C
      IF(ITER.LE.2) GO TO 5555
      ECF=0.
      PGWORK = 0.
      DO 601 K=1,KB
      UE = 0.5*(UF(K) + UB(K))
      VE = 0.5*(VF(K) + VB(K))
C
      IF(K.EQ.(KB-2)) THEN
      UEB = UE
      VEB = VE
      END IF
C
      ECF = ECF + (UE**2 + VE**2)*DZ(K)*H/2.
      PGWORK = PGWORK - DT1*(UE*PGX + VE*PGY)*DZ(K)*H
601   CONTINUE
      WORK = DT1*(U(1)*TOX(ITER)+V(1)*TOY(ITER))
      WORK = WORK + DT1*(WUBOT*UEB + WVBOT*VEB)
C
      SHRP=0.
      BUOYP = 0.
C
      DO 602 K=2,KB
      AUD = (U(K-1) - U(K))/(H*DZ(K-1))
      AVD = (V(K-1) - V(K))/(H*DZ(K-1))
      RHODZ = (RHO(K-1) - RHO(K))/(H*DZ(K-1))
      BUOYP = BUOYP - 9.8*KH(K)*RHODZ*(H*DZ(K-1))
      SHRP = SHRP - KM(K)*(AUD**2 + AVD**2)*(H*DZ(K-1))
 602  CONTINUE
      SHRP = SHRP*DT1
      BUOYP = BUOYP*DT1
C
      BILEC1 = ECF - ECB
      BILEC2 = BILEC1 + (WORK + SHRP + PGWORK)
      BILECT = BILECT + DT1*BILEC2 
      ECB = ECF
 5555 CONTINUE
C
 603  CONTINUE
C
C      WRITE (11,9344) BILEC1, WORK, PGWORK, SHRP, BUOYP, WUBOT, UE,
C     X  UBSUM, TXBSUM
C 9344 FORMAT (1X,/,1X,'E', 9E8.2)
C
C
C      NPLT = 9
C      ICON=1
C      MNPLT = ITER/NPLT
C      IF (ITER .EQ. MNPLT*NPLT) THEN
C      DO 401 K=1,KB
C        IF (ICON .EQ. 1) THEN
C         UX(MNPLT,K)=U(KB-K+1)
C         VX(MNPLT,K)=V(KB-K+1)
C         DX(MNPLT,K)=RHO(KB-K+1)
C         TTX(MNPLT,K)=T(KB-K+1)
C         EZ1(K)=FLOAT(K-1)*0.1
C        ELSE
C         UX(MNPLT,K)=U(K)*50.+FLOAT(MNPLT)*2.
C         VX(MNPLT,K)=V(K)*50.+FLOAT(MNPLT)*2.
C         DX(MNPLT,K)=RHO(k)*50.+FLOAT(MNPLT)*2.
C         TTX(MNPLT,K)=T(K)*50.+FLOAT(MNPLT)*2.
C         EZ1(K)=FLOAT(K-1)*.5
C        END IF
C 401  CONTINUE
C      END IF
C
C  Write out a few things on every IWRTth time step.
C
      IF(IWRT.LT.1) GO TO 1339
C
      IF(ITER.EQ.1) WRITE (6,444) 
  444 FORMAT (1X,/,1X,'MODEL SOLUTION VARIABLES ARE:',/,1X,
     X '   M TIMED    QI     QL      TX   T(1)   S(1)   U(1) ',
     X '  V(1)  DML   TLD',/,1X)
C
      CALL MLDPTH(ZZ,RHO,Q2,KBM1,ZZMLD,TLD)
      ZZDMLD=-ZZMLD*H
      TLD = -TLD*H 
C
      IF(MOD(ITER,IWRT).NE.1) GO TO 1339
      WRITE (6,133) ITER,TIMED,QI,QL,TX,T(1),S(1),
     C U(1),V(1),ZZDMLD,TLD
  133 FORMAT (1X,I4,F6.2,2F7.0,F7.2,4F7.2,3F6.0)
C
 1339 CONTINUE
C
C  WRITE OUT DATA TO UNIT 11 EVERY TIME STEP
C
      TXB = 1023.*WUBOT
      TYB = 1023.*WVBOT
      TXS = 1023.*TOX(ITER)
      TYS = 1023.*TOY(ITER)
C
C
      IF(ITER.EQ.1) WRITE (11,305)
  305 FORMAT (1X,'ITER',2X,'TIMED',5X,'QI',5X,'QL',5X,'TX',5X,
     C 'TY',6X,'T',6X,'S',6X,'U',6X,'V',4X,'MLD',
     C 4X,'TXB',4X,'TYB',6X,'DEKE',4X,'BLWORK',4X,'PGWORK',4X,
     C 'SHRPRO')
      WRITE (11,500)   ITER,TIMED,RA(ITER), HSURF(ITER),
     C TXS,TYS,T(1),S(1),
     C U(1),V(1), ZZDMLD, TXB, TYB, BILEC1, WORK, PGWORK, SHRP
  500 FORMAT (1X,I4,F7.2, 2F7.0, 6F7.2, F7.0, 2F7.2, 4E10.3)
9001  CONTINUE
C
C  Store profile data on logical unit 13, if requested.
C
      IF(ISTOR.EQ.1) THEN
      IF(MOD(ITER,ITFREQ).EQ.1) THEN
      WRITE (13,3183)
 3183 FORMAT (1X)
      WRITE (13,3182) 
 3182 FORMAT (1X,5X,'TIMED',9X,'Z',9X,'T',9X,'S',9X,'U',9X,
     X 'V',8X,'Q2',8X,'KH',8X,'KM',8X,'KQ',8X,'SH',
     X 8X,'SM',5X,'BPROD',5X,'SPROD',9X,'Z')
      DO 3189 K=1,KB,IZFREQ
      WRITE (13,3188) TIMED,ZP(K),T(K),S(K),U(K),V(K),
     X Q2(K), KH(K), KM(K), KQ(K), SH(K), SM(K), 
     X BPROD(K), SPROD(K), ZP(K)
 3188 FORMAT (1X,6F10.2,8E10.3,F10.2)
 3189 CONTINUE
      END IF
      END IF
C
C
 9009 CONTINUE
C
C  Store the density profile if this is the first step.
C
      IF(ITER.EQ.1) THEN
      DO 2278 J=1,KB
      DO(J) = RHO(J)
 2278 CONTINUE
      END IF
C
C  This is the place to compute diagnostics.
C
C  Find max/min values on the first day.
C
      IF(ITER.EQ.1) DMLMIN = 100.
      IF(ITER.EQ.1) TMIN = T(1)
      IF(TIMED.LT.1.) THEN
      IF(T(1).LT.TMIN) TMIN = T(1)
      IF(T(1).GT.TMAX) TMAX = T(1)
      SPD = SQRT(U(1)**2 + V(1)**2)
      IF(SPD.GT.SPDMAX) SPDMAX = SPD
      END IF
C
C  Compute the net heat flux as a check on heat conservation.
C
      QNET = QNET + DT1*(RA(ITER) + HSURF(ITER))
C     
C
 9000 CONTINUE
  499 CONTINUE
C
C  Compute heat and potential energy diagnostics.
C  DPE will be the change in potential energy over the
C  integration period, and HEAT is the change in heat
C  content. HEAT should nearly equal the time-integrated
C  surface heat flux, QNET, if heat has been conserved.
C
      DPE = 0.
      HEAT = 0.
C
      DO 8376 J=1,KBM1
      CALL INTX(ZO,TO,NSOBS,NTS,ZP(J),TI)
      TANOM = T(J) - TI
      HEAT = HEAT + TANOM*DZP(J)
      IF(J.EQ.1) TANOMS = TANOM
C
      DANOM = RHO(J) - DO(J)
      DPE = DPE + 9.8*1000.*DANOM*ZP(J)*DZP(J)
C
C      WRITE (13,4465) J,DO(J),RHO(J),DANOM,DPE
 4465 FORMAT (1X,'DO,D,DANOM,DPE',I4,4E12.4)
 8376 CONTINUE
C
      QNET = QNET/(1.023E3*4.183E3)
      WRITE (6,3877) HEAT, QNET, DPE
 3877 FORMAT (1X,/,1X,'HEAT, QNET, DPE ARE', 3F9.2)
C
C  End of diagnostics.
C
C
C          WRITE (10) UX
C          WRITE (10) VX
C          WRITE (10) DX
C          WRITE (10) TTX
C
C
   40 CONTINUE
C
C
      WRITE (6,870) TMIN,TMAX,DMLMIN,SPDMAX
  870 FORMAT (1X,/,1X,'TMIN, TMAX, DMLMIN, SPDMAX ARE',4F8.2)
C
C
      STOP
      END
C
C
      SUBROUTINE INTX(R,S,KB,N,RI,SI)
C
C  This subroutine interploates the array S to find the value 
C  at the point RI, where R is the independent variable.
C  The interpolation is linear, and assumes that R increases
C  montonically, and that RI is within the range of R. If
C  RI is not within the range of R, SI is set equal to the
C  appropriate endpoint of S, and a warning is printed. 

C
      DIMENSION R(KB), S(KB)
C 
C  Check for out of range RI.
C
      IF(RI.LT.R(1)) THEN
      SI = S(1)
      WRITE (6,20) RI, R(1)
  20  FORMAT (1X,/,1X,'WARNING, RI WAS BELOW R(1) IN CALL TO INTX',/,
     C 'RI AND R(1) WERE', 2E12.4)
      END IF
C
C
      IF(RI.GT.R(N)) THEN
      SI = S(N)
      WRITE (6,21) RI, R(N)
  21  FORMAT (1X,/,1X,'WARNING, RI WAS ABOVE R(N) IN CALL TO INTX',/,
     C 'RI AND R(N) WERE', 2E12.4)
      END IF    
C 
C  Data are within bounds for an interpolation.
C
      NM = N - 1
      DO 2 K=1,NM
      IF(RI.GE.R(K).AND.RI.LE.R(K+1)) GO TO 1
    2 CONTINUE
C
      GO TO 9
C
    1 CONTINUE
      SI = (S(K)*(R(K+1) - RI) + S(K+1)*(RI - R(K)))/
     X (R(K+1) - R(K))
C
C
    9 CONTINUE
      RETURN
      END
C
C
c----------------------------------------------------------------------
      SUBROUTINE MLDPTH(ZZ,R,Q2,KBM1,ZZMLD,TLD)
      PARAMETER(KB=150)
      DIMENSION ZZ(KB),R(KB),Q2(KB)
C
      REPS = 0.00002
      QEPS = 0.1
      KML = 1
      ZZMLD = 1.0
C
      DO 100 K=1,KBM1
      IF(ABS(R(1) - R(K)).GT.REPS) GO TO 100
      ZZMLD = ZZ(K) 
      KML = K
  100 CONTINUE
C
      QS = 0.
      DO 1 J=1,KML
      QS = QS + Q2(J)
    1 CONTINUE
      QA = QS/FLOAT(KML)
C
      DO 2 J=2,KBM1
      IF(Q2(J).LT.(QEPS*QA)) GO TO 3
      TLD = ZZ(J)
    2 CONTINUE
    3 CONTINUE
C
      RETURN
      END
c----------------------------------------------------------------------
      SUBROUTINE PROFQ(DT2)
      PARAMETER(KB=150)
      COMMON/BLKT/TF(KB),T(KB),TB(KB),SF(KB),S(KB),SB(KB)
      COMMON/BLKU/UF(KB),U(KB),UB(KB),VF(KB),V(KB),VB(KB)
      COMMON/BLKQ/Q2F(KB),Q2(KB),Q2B(KB)
      COMMON/BLKR/RHO(KB),WT(KB)
      COMMON/BLKWU/WUSURF,WVSURF,WUBOT,WVBOT
      COMMON/BLKWT/WTSURF,WSSURF,PROFIL(KB)
      COMMON/BLKZ/Z(KB),ZZ(KB),DZ(KB),DZZ(KB),UMOL
      COMMON/BLKH/H,DT
      COMMON/BLKA/A(KB),C(KB),VH(KB),VHP(KB)
      COMMON/BLKK/KM(KB),KH(KB),KQ(KB),L(KB)
      COMMON/INDS/KBM1,KBM2
      COMMON/MODEL/GM(KB),GH(KB),SM(KB),SH(KB),KN(KB),
     1     SPROD(KB),BPROD(KB),PROD(KB),DTEF(KB)
      COMMON/CSTUR/A1,B1,A2,B2,C1,D1,D2,A1A16,A1A112,A1A29,A1A26,
     1A2B23,A1A212

      REAL KM,KH,KQ,L
      REAL KAPPA,KN
      DATA KAPPA/0.40/,SQ/0.2/
      DATA GEE/9.806/
      DATA SMALL/1.E-10/,EPS/1.E-30/
C      DATA SM/KB*0.39/,SH/KB*0.49/,GM/KB*0.154/,GH/KB*0.154/
C
      DO 8877 J=1,KB
      SM(J) = 0.39
      SH(J) = 0.49
      GM(J) = 0.154
      GH(J) = 0.154
 8877 CONTINUE
C
      DO 100 K=2,KBM1
      A(K)=-DT2*(KQ(K+1)+KQ(K)+2.*UMOL)*.5/(DZZ(K-1)*DZ(K)*H*H)
      C(K)=-DT2*(KQ(K-1)+KQ(K)+2.*UMOL)*.5/(DZZ(K-1)*DZ(K-1)*H*H)
  100 CONTINUE
C*********************************************************************
C
C        THE FOLLOWING SECTION SOLVES THE EQUATION
C        DT2*(KQ*Q2')' - Q2*(2.*DT2*DTEF+1.) =Q2B
C
C********************************************************************
      CONST1=16.6**.6666667
      VH (1)=0.0
      VHP(1)=SQRT(WUSURF**2+WVSURF**2)*CONST1
      Q2F(KB)=SQRT(WUBOT**2+WVBOT**2)*CONST1
      DO 101 K=2,KBM1
      DTEF(K)=SQRT(Q2B(K))/(B1*L(K)+EPS)
      SPROD(K)=KM(K)*((U(K)-U(K-1))**2+(V(K)-V(K-1))**2)/
     1(DZZ(K-1)*H)**2
      BPROD(K)=-KH(K)*GEE*(RHO(K)-RHO(K-1))/(DZZ(K-1)*H)
  101 PROD(K)=SPROD(K)+BPROD(K)
      DO 102 K=2,KBM1
      VHP(K)=1./(A(K)+C(K)*(1.-VH(K-1))-(2.*DT2*DTEF(K)+1.))
      VH(K)=A(K)*VHP(K)
      VHP(K)=(-2.*DT2*PROD(K)+C(K)*VHP(K-1)-Q2B(K))*VHP(K)
  102 CONTINUE
      DO 103 K=1,KBM1
      KI=KB-K
      Q2F(KI)=VH(KI)*Q2F(KI+1)+VHP(KI)
  103 CONTINUE
      DO 112 K=2,KBM1
      IF(Q2F(K).GT.SMALL) GO TO 112
      Q2F(K)=SMALL
  112 CONTINUE
C******************************************************************
C      THE FOLLOWING SECTION SOLVES FOR KM AND KH
C******************************************************************
      a1l=0.
      a2l=0.
      do 211 k=1,kb
      a1l=a1l+sqrt(q2(k))*dz(k)*h*z(k)*h
211   a2l=a2l+sqrt(q2(k))*dz(k)*h
      rl0=-.20*a1l/a2l
      do 210 k=1,kb
      rl1=1.0
c longueur variable en enlevant c ligne suivante
c      rl1=.40*z(k)*(-h)/(rl0+.4*z(k)*(-h))
      l(k)=rl0*rl1
210   continue
      DO 212 K=2,KBM1
  212 GH(K)=L(K)**2/Q2(K)*((U(K)-U(K-1))**2+(V(K)-V(K-1))**2)/
     1(-DZZ(K-1)*H)**2

C OCCASIONALLY SOLUTION SPLITTING W.R.T. Z OCCURS.THIS IS REMOVED BY TH
C WHICH DOES NOT SEEM TO EFFECT MIXED LAYER DEEPENING.
      DO 213 K=2,KBM1
  213 GM(K)=.25*GH(K-1)+.5*GH(K)+.25*GH(K+1)
      DO 214 K=2,KBM1
      GH(K)=L(K)**2/Q2(K)/(-DZZ(K-1)*H)*GEE*(RHO(K)-RHO(K-1))
      GH(K)=AMIN1(GH(K),.032)
      GM(K)=AMIN1(GM(K),.48-15.*GH(K))
      C11=1.+A1A16*GM(K)-A1A29*GH(K)
      C12=(-A1A112-A1A29)*GH(K)
      C21=A1A26*GM(K)
      C22=1.-(A2B23+A1A212)*GH(K)
      DENOM=1./(C11*C22-C12*C21)
      SM(K)=(C22*D1-C12*D2)*DENOM
      SH(K)=(C11*D2-C21*D1)*DENOM
  214 CONTINUE
      DO 216 K=2,KBM1
      KN(K)=L(K)*SQRT(Q2(K))
C***  KQ(K)=(KN(K)*SQ+KQ(K))*.5
      KQ(K)=(KN(K)*.41*SM(K)+KQ(K))*.5
      KM(K)=(KN(K)*SM(K)+KM(K))*.5
      KH(K)=(KN(K)*SH(K)+KH(K))*.5
  216 CONTINUE

      RETURN
      END
c----------------------------------------------------------------------
      SUBROUTINE PROFT(FF,FB,WFSURF,DT2,ra)
      PARAMETER(KB=150)
      COMMON/BLKT/TF(KB),T(KB),TB(KB),SF(KB),S(KB),SB(KB)
      COMMON/BLKWT/WTSURF,WSSURF,profil(kb)
      COMMON/BLKZ/Z(KB),ZZ(KB),DZ(KB),DZZ(KB),UMOL
      COMMON/BLKH/H,DT
      COMMON/BLKA/A(KB),C(KB),VH(KB),VHP(KB)
      COMMON/BLKK/KM(KB),KH(KB),KQ(KB),L(KB)
      COMMON/INDS/KBM1,KBM2
      REAL KM,KH,KQ,L
      REAL KAPPA
      DIMENSION FF(KB),FB(KB)
      DATA KAPPA/0.4/
      UMOLPR=UMOL
C************************************************************************************
C        THE FOLLOWING SECTION SOLVES THE EQUATION
C        DT2*(KH*FF')' -FF = FB
C*************************************************************************************
      do 96 k=1,kbm1
      fb(k)=fb(k)+profil(k)*ra*dt2
96    continue
      FF(KBM1)=FB(KBM1)
      DO 100 K=2,KBM1
      A(K-1)=-DT2*(KH(K)+UMOLPR)/(DZ(K-1)*DZZ(K-1)*H*H)
      C(K)=-DT2*(KH(K)+UMOLPR)/(DZ(K)*DZZ(K-1)*H*H)
  100 CONTINUE
C*************************************************************************************
C   INSTEAD OF HEAT FLUX ONE CAN IMPOSE A SURFACE TEMPERATURE
C   BY REPLACYNG THE NEXT TWO LINES BY VH(1)=0 and VHP(1)=F(1).
C   F(1) SHOULD BE SPECIFIED IN THE MAIN PROGAM
C***************************************************************************************
      VH(1)=A(1)/(A(1)-1.)
      VHP(1)=(-DT2*WFSURF/(DZ(1)*H)-FB(1))/(A(1)-1.)
   98 CONTINUE
      DO 101 K=2,KBM2
      VHP(K)=1./(A(K)+C(K)*(1.-VH(K-1))-1.)
      VH(K)=A(K)*VHP(K)
      VHP(K)=(C(K)*VHP(k-1)-FB(K))*VHP(K)
  101 CONTINUE
C************************************************************************************
C   INSTEAD OF MATCHING SOLUTION TO A LOWER LAYER VALUE OF F(KB),
C   ONE MAY IMPOSE AN ADIABATIC BOTTOM B.C. BY REMOVING C'S FROM
C   COL. 1 OF THE NEXT TWO LINES.
C*************************************************************************************
      FF(KBM1)=(C(KBM1)*VHP(KBM2)-FB(KBM1))
     1   /(C(KBM1)*(1.-VH(KBM2))-1.)
   99 CONTINUE
      DO 102 K=2,KBM1
      KI=KB-K
      FF(KI)=VH(KI)*FF(KI+1)+VHP(KI)
  102 CONTINUE
      RETURN
      END
c----------------------------------------------------------------------
      SUBROUTINE PROFC(CB,CF,WSURF,CBC,DT2)
      PARAMETER(KB=150)
      COMMON/BLKU/UF(KB),U(KB),UB(KB),VF(KB),V(KB),VB(KB)
      COMMON/BLKWU/WUSURF,WVSURF,WUBOT,WVBOT
      COMMON/BLKZ/Z(KB),ZZ(KB),DZ(KB),DZZ(KB),UMOL
      COMMON/BLKH/H,DT
      COMMON/BLKA/A(KB),C(KB),VH(KB),VHP(KB)
      COMMON/BLKK/KM(KB),KH(KB),KQ(KB),L(KB)
      COMMON/INDS/KBM1,KBM2
      REAL KM,KH,KQ,L
      DIMENSION  CF(KB),CB(KB)
      DO 100 K=2,KBM1
      A(K-1)=-DT2*(KM(K)+UMOL  )/(DZ(K-1)*DZZ(K-1)*H*H)
      C(K)=-DT2*(KM(K)+UMOL  )/(DZ(K)*DZZ(K-1)*H*H)
  100 CONTINUE
      VH(1)=A(1)/(A(1)-1.)
      VHP(1)=(-DT2*WSURF/(DZ(1)*H)-CB(1))/(A(1)-1.)
      DO 101 K=2,KBM2
      VHP(K)=1./(A(K)+C(K)*(1.-VH(K-1))-1.)
      VH(K)=A(K)*VHP(K)
      VHP(K)=(C(K)*VHP(k-1)-CB(K))*VHP(K)
  101 CONTINUE
      CF(KBM1)=(C(KBM1)*VHP(KBM2)-CB(KBM1))/(CBC
     1     *DT2/(-DZ(KBM1)*H)-1.-(VH(KBM2)-1.)*C(KBM1))
      DO 103 K=2,KBM1
      KI=KB-K
      CF(KI)=VH(KI)*CF(KI+1)+VHP(KI)
  103 CONTINUE
   92 WUBOT=-CBC*CF(KBM1)
      RETURN
      END

c----------------------------------------------------------------------

      SUBROUTINE DENS(TI,SI,RHOO)
      PARAMETER(KB=150)
      DIMENSION TI(KB),SI(KB),RHOO(KB)
C
C   ......... THIS SUBROUTINE COMPUTES (DENSITY-1)*0.001 .............
C
      DO 61 K=1,KB
      RHOO(K)=SGT(TI(K),SI(K),SG)*1.E-3
 61   CONTINUE
      RETURN
      END
c----------------------------------------------------------------------

      SUBROUTINE DEPTH(Z,ZZ,DZ,DZZ,KB,KBM1)
      DIMENSION Z(KB),ZZ(KB),DZ(KB),DZZ(KB)
      KL1=3
      KL2=KB-2
C**********************************************************************
C   THIS SUBROUTINE ESTABLISHES THE VERTICAL RESOLUTION WITH LOG
C   DISTRIBUTIONS AT THE TOP AND BOTTOM AND A LINEAR DISTRIBUTION
C   BETWEEN. CHOOSE KL1 AND KL2. THE DEFAULT KL1 = .3*KB AND KL2 = KB-2
C   YIELDS A LOG DISTRIBUTION AT THE TOP AND NONE AT THE BOTTOM.
C   KL1=3 AND KL2=KB-2 YIELDS NONE LOG AT THE TOP AND NONE AT BOT.
C**********************************************************************
      BB=FLOAT(KL2-KL1)+4.
      CC=FLOAT(KL1)-2.
      DEL1=2./BB/EXP(.693147*FLOAT(KL1-2))
      DEL2=2./BB/EXP(.693147*FLOAT(KB-KL2-1))
      Z(1)=0.
      ZZ(1)=-DEL1/2.
      DO 3 K=2,KL1
      Z(K)=-DEL1*EXP(.693147*FLOAT(K-2))
    3 ZZ(K)=-DEL1*EXP(.693147*(FLOAT(K)-1.5))
      DO 4 K=KL1,KL2
      Z(K)=-(FLOAT(K)-CC)/BB
    4 ZZ(K)=-(FLOAT(K)-CC+0.5)/BB
      DO 5 K=KL2,KBM1
      Z(K)=(1.0-DEL2*EXP(.693147*FLOAT(KB-K-1)))*(-1.)
    5 ZZ(K)=(1.0-DEL2*EXP(.693147*(FLOAT(KB-K)-1.5)))*(-1.)
      Z(KB)=-1.0
      ZZ(KBM1)=-1.*(1.-DEL2/2.)
      ZZ(KB)=-1.*(1.+DEL2/2.)
      DO 6 K=1,KBM1
      DZ(K)=Z(K)-Z(K+1)
    6 DZZ(K)=ZZ(K)-ZZ(K+1)
      RETURN
      END
c----------------------------------------------------------------------
C
C
      FUNCTION SG0(S) 
C 
C  A sigma-0 subroutine neede by the sigma-t subroutine
C  taken from SEAPROP.
C
C  SIGMA-0 KNUDSEN 
      SG0 = ((6.76786136E-6*S-4.8249614E-4)*S+0.814876577)*S  
     X -0.0934458632
      RETURN
      END
C
C
      FUNCTION SGT(T,S,SG)
C
C  A sigma-t subroutine taken from SEAPROP.
C
C SIGMA-T KNUDSEN 
      SG = SG0(S) 
   20 SGT = ((((-1.43803061E-7*T-1.98248399E-3)*T-0.545939111)*T
     X +4.53168426)*T)/(T+67.26)+((((1.667E-8*T-8.164E-7)*T 
     X +1.803E-5)*T)*SG+((-1.0843E-6*T+9.8185E-5)*T-4.7867E-3)*T
     X +1.0)*SG 
      RETURN
      END
