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