/[MITgcm]/MITgcm_contrib/jscott/igsm/src/eddypa.F
ViewVC logotype

Annotation of /MITgcm_contrib/jscott/igsm/src/eddypa.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Fri Aug 11 19:35:30 2006 UTC (18 years, 11 months ago) by jscott
Branch: MAIN
CVS Tags: HEAD
atm2d package

1 jscott 1.1 #include "ctrparam.h"
2    
3     ! ==========================================================
4     !
5     ! EDDYPA_NEW.F: THIS ROUTINE DEALS WITH THE PARAMETERIZATION OF
6     ! EDDY TRANSPORTS AND VARIANCES
7     ! Modified version, includes real averaging
8     !
9     !
10     SUBROUTINE EDDYPA (UT,VT,TT,PT,QT,U,V,T,P,Q,DTE,EDDFLX) 2181.
11     C 2181.5
12     C THIS ROUTINE DEALS WITH THE PARAMETERIZATION OF 2182.
13     C EDDY TRANSPORTS AND VARIANCES 2182.5
14     C 2183.
15     #if ( defined CPL_CHEM )
16     !
17     #include "chem_para"
18     #include "chem_com"
19     !
20     #endif
21    
22     #include "BD2G04.COM"
23    
24     COMMON/WORK3/PHI(IM0,JM0,LM0) 2184.
25     COMMON/WORK4/FD(IM0,JM0) 2184.5
26     DIMENSION UT(IM0,JM0,LM0),VT(IM0,JM0,LM0),TT(IM0,JM0,LM0),
27     & PT(IM0,JM0),
28     * QT(IM0,JM0,LM0) 2185.5
29     DIMENSION SINL(JM0),COSL(JM0), 2186.
30     * DXU(JM0),DYU(JM0),DXYU(JM0),PHIS(IO0,JM0),UX(IO0,JM0,2*LM0) 2186.5
31     EQUIVALENCE (SINL(1),SINP(1)),(COSL(1),COSP(1)), 2187.
32     * (DXU(1),DXV(1)),(DYU(1),DYV(1)),(DXYU(1),DXYV(1)), 2187.5
33     * (PHIS,FDATA) 2188.
34     c REAL*4 KAPAP1 2188.5
35     REAL KAPAP1,MOMCOEF
36     COMMON/EPARA/VTH(JM0,LM0),WTH(JM0,LM0),VU(JM0,LM0),VV(JM0,LM0)
37     & ,DQSDT(JM0,LM0) 2189.
38     * ,DWV(JM0),PHIT(JM0,LM0),TPRIM2(JM0,LM0),WU(JM0,LM0),CKS,CKN 2189.5
39     * ,WQ(JM0,LM0),VQ(JM0,LM0) 2190.
40     DIMENSION PHITAV(LM0),TPRIM0(JM0)
41     COMMON/SPEC2/KM,KINC,COEK,C3LAND(IO0,JM0),C3OICE(IO0,JM0) 2190.5
42     * ,C3LICE(IO0,JM0),WMGE(IO0,JM0),TSSFC(1,JM0,4) 2191.
43     DIMENSION TSZM(JM0),THSZM(JM0),FO(JM0),WUDATA(JM0,LM0 )
44     & ,VUDATA(JM0,LM0), 2191.5
45     & VSH(JM0,LM0),WSH(JM0,LM0),TEMZ(JM0,LM0),COEKY(JM0)
46     & ,DTHDY(JM0,LM0),TEMJL(JM0,LM0+1) 2192.
47     * ,DSHDY(JM0,LM0),DSHDP(JM0,LM0+1),DTHDZ(JM0,LM0+1)
48     & ,DSHDZ(JM0,LM0+1) 2192.5
49     * ,SHSAT(JM0,LM0),DSHDYS(JM0,LM0) 2193.
50     DIMENSION DQDYJ(JM0),DQDY(JM0,LM0),UZM(JM0),PSROK(JM0,LM0)
51     & ,DQSDT1(JM0,LM0), 2193.5
52     * DUDP(JM0,LM0),TZA(JM0,LM0),DTHDP(JM0,LM0+1),PRESR(JM0,LM0+1)
53     & ,PSRK(JM0,LM0+1) 2194.
54     * ,PSI(JM0,LM0),PVOR(JM0,LM0),VQZ(JM0),UNEW(JM0,LM0),UOLD(JM0,LM0)
55     & ,FMIX(JM0,LM0) 2194.5
56     * ,SCLH(JM0,LM0),GAMA(JM0),SHJ(JM0),HTJ(JM0),ZSTAR(JM0,LM0)
57     & ,BVF(JM0) 2195.
58     * ,PRSRU(JM0,LM0),DTDYJ(JM0),DTDZJ(JM0),DWVJL(JM0,LM0)
59     & ,BVFJL(JM0,LM0), 2195.5
60     * DUDZ(JM0,LM0),DUDZJ(JM0),COEJL(JM0,LM0),DQDY1(JM0,LM0)
61     & ,COEJL1(JM0,LM0) 2196.
62     * ,UVAR(JM0),TEMJ(JM0),WTDATA(JM0,LM0),WQDATA(JM0,LM0)
63     & ,VTDATA(JM0,LM0), 2196.5
64     * VQDATA(JM0,LM0),VVDATA(JM0,LM0),TTDATA(JM0,LM0),GAMAJ(JM0)
65     & ,DTHDZJ(JM0) 2197.
66     * ,BVFM(JM0,LM0),BVFMJ(JM0),DQDZ(JM0,LM0),VLAT(JM0),VVN(9)
67     & ,VVS(9)
68     & ,VVNEWA(LM0),FOV(JM0),BETAV(JM0),BETAP(JM0),TEMKGP(JM0)
69     & ,TEMKGV(JM0),RDWVERR(JM0),ADWVERR(JM0)
70     & ,TEMKV(JM0),TEMKP(JM0),GAMAP(JM0),DWVP(JM0),HTJP(JM0)
71     & ,DWVA(JM0),DWVPA(JM0)
72     dimension TZN(JM0,LM0),QZN(JM0,LM0),UZN(JM0,LM0),VZN(JM0,LM0)
73     common/tprmtg/tprmg(JM0),ntprmg(JM0)
74     common/Dscale/DWAV0(JM0)
75     logical first,uvtr,wvtr,uttr,wttr,uqtr,wqtr,EDDFLX
76     &,OLDUVVERT,HPRNT
77     SAVE
78     data first/.true./
79     DATA VVN/37.6,39.5,57.8,94.8,150.6,236.3,329.5,114.8,43.8/ 2261.5
80     DATA VVS/32.7,28.3,33.7,44.9,79.4,151.5,193.6,90.9,15.7/ 2262.
81     C 2262.5
82     QSAT(TM,PR,QL)=3.797915*EXP(QL*(7.93252E-6-2.166847E-3/TM))/PR 2263.
83     C 2263.5
84     c HPRNT=TAU.gt.50.0.and.TAU.lt.75.0
85     HPRNT=.FALSE.
86     c HPRNT=TAU.ge.17520.0.and.TAU.lt.17521.0
87     c HPRNT=.TRUE.
88     JPR=24
89     uvtr=EDDFLX
90     wvtr=EDDFLX
91     uttr=EDDFLX
92     wttr=EDDFLX
93     uqtr=EDDFLX
94     wqtr=EDDFLX
95     if(first)then
96     LZ=LM 2264.
97     LZM1=LZ-1 2264.5
98     LZM2=LZ-2 2265.
99     JHALF=JM/2 2265.5
100     JHAFM1=JHALF-1 2266.
101     JHAFP2=JHALF+2 2266.5
102     JMM2=JM-2 2267.
103     IX=IM 2267.5
104     IS=IM 2268.
105     FIS=IS 2268.5
106     CKS=1. 2268.6
107     CKN=1. 2268.7
108     CP=RGAS/KAPA 2269.5
109     SHA=RGAS/KAPA 2270.
110     CONSTA=2.0
111     CONSTA=1.5
112     CONSTA=1.0
113     RCZZ=1.35
114     c RCZZ=1.83
115     MOMCOEF=1.
116     c MOMCOEF=4.*0.5*(0.57+0.71)/CONSTA
117     KAPAP1=KAPA+1. 2271.
118     HLAT=LHE 2271.5
119     CLH=HLAT/CP 2272.
120     DQDTX=.622*HLAT/RGAS 2272.5
121     P1000=1000.**KAPA 2273.
122     DTDYMN=1.E-6/P1000 2273.5
123     DQDYMN=1.E-09 2274.
124     DSCALE=626.*1.e3
125     c DSCALE=323.*1.e3
126     c DSCALE=1000.*1.e3
127     FDR=exp(-DLAT*RADIUS/DSCALE)
128     OLDUVVERT=.false.
129     DWVCRI=1000.
130     DWVCRI=100.
131     JINT=6
132     IF(JM.EQ.46)JINT=12
133     c JINT=JM/2-2
134     JEQ=JM/2+1
135     JEQP1=JEQ+1
136     JEQM1=JEQ-1
137     JEQM2=JEQ-2
138     LBGN=1
139     ALPHA=1.
140     NVAIT=2
141     NVAIT=3
142     WMGMIN=3.0
143     WMGMIN=0.0
144     print *,' Mom. flux param. for dry atmosph.'
145     print *,' New vertical avereging'
146     print *,' Change in YCUT calculations'
147     print *,' with tg difference'
148     print *,' CONSTA=',CONSTA
149     print *,' RCZZ=',RCZZ
150     print *,' DWVCRI=',DWVCRI
151     print *,' MRCH=',MRCH
152     print *,' FDR=',FDR,' L=',DSCALE
153     print *,' JINT=',JINT
154     print *,' NVAIT=',NVAIT
155     print *,' PSF=',PSF
156     print *,' WMGMIN=',WMGMIN
157     print *,' With A for momentum=',CONSTA*MOMCOEF
158     print *,' UVTR=',uvtr,' WVTR=',wvtr
159     print *,' UTTR=',uttr,' WTTR=',wttr
160     print *,' UQTR=',uqtr,' WQTR=',wqtr
161     if(OLDUVVERT)then
162     print *,' Vertical structure of UV flux from GISS
163     & 2D model'
164     endif
165     DO 3 J=1,JM 2274.5
166     3 FO(J)=F(J)/DXYP(J) 2275.
167     DO 4 J=2,JM 2275.5
168     VLAT(J)=(LAT(J)+LAT(J-1))*.5 2276.
169     FOV(J)=0.5*(FO(J)+FO(J-1))
170     BETAV(J)=(FO(J)-FO(J-1))/DYU(J)
171     if(j.NE.JM)BETAP(J)=.5*(FO(J+1)-FO(J-1))/DYU(J)
172     TEMY=RADIUS*(.25*TWOPI-ABS(VLAT(J)))
173     Xk=1./(RADIUS*COSV(j))
174     Yk=TWOPI/(4.*TEMY)
175     TEMKGV(j)=SQRT(Xk**2+Yk**2)
176     4 CONTINUE
177     TPR00=1.0
178     TPR90=2.5
179     TPR11=0.5*(TPR90+TPR00)
180     TPR12=0.5*(TPR90-TPR00)
181     print *,'TPR00=',TPR00,' TPR90=',TPR90
182     do j=1,jm
183     TEMY=RADIUS*(.25*TWOPI-ABS(LAT(J)))
184     Xk=1./(RADIUS*COSP(j))
185     Yk=TWOPI/(4.*TEMY)
186     TEMKGP(j)=SQRT(Xk**2+Yk**2)
187     tprimrad = 3.14159*(float(j-1)/22.5)
188     TPRIM0(j)=5.5+(4.5*cos(tprimrad)) ! 2323
189     TPRIM0(j)=0.0 ! 2324
190     ! TPRIM0(j)=1.5+(1.0*cos(tprimrad)) ! 2327 2332
191     TPRIM0(j)=TPR11+(TPR12*cos(tprimrad)) ! 2327 2332
192     ! TPRIM0(j)=2.0 ! 2331
193     enddo
194     c print *,'LAT'
195     c print *,LAT
196     c print *,'VLAT'
197     c print *,VLAT
198     c print *,'COSV'
199     c print *,COSV
200     print *,'TPRIM0'
201     print *,TPRIM0
202     first=.false.
203     endif
204     C 2276.5
205     #if ( defined CPL_CHEM )
206     !
207     do 1 j=1,nlat
208     beta1(j)=0.0
209     beta2(j)=0.0
210     do 2 k=1,nlev
211     beta3(j,k)=0.0
212     beta4(j,k)=0.0
213     fkt (j,k)=0.0
214     2 continue
215     1 continue
216     !
217     #endif
218    
219     C reassign T,Q,U and V to internal variables
220    
221     do l=1,lm
222     do j=1,jm
223     TZN(j,l)=T(1,j,l)
224     QZN(j,l)=Q(1,j,l)
225     UZN(j,l)=U(1,j,l)
226     VZN(j,l)=V(1,j,l)
227     enddo
228     enddo
229     if(HPRNT)then
230     print *,' Eddypa UZN before filter'
231     print *,'LZM2=',LZM2
232     print *,(UZN(13,L),L=1,LM)
233     print *,(UZN(14,L),L=1,LM)
234     endif
235     if(JM.eq.46)then
236     CALL FLTR(TZN,JM,LM,1)
237     CALL FLTR(QZN,JM,LM,1)
238     CALL FLTR(UZN,JM,LM,2)
239     endif
240    
241     DO 6 J=1,JM 2277.
242     DWV(J)=0. 2277.5
243     DO 6 I=1,IM 2278.
244     6 WMGE(I,J)=WMGMIN 2278.5
245     DO 10 L=1,LM 2279.
246     DO 10 J=1,JM 2279.5
247     VU(J,L)=0. 2280.
248     VQ(J,L)=0. 2280.5
249     VTH(J,L)=0. 2281.
250     WTH(J,L)=0. 2281.5
251     WU(J,L)=0. 2282.
252     WQ(J,L)=0. 2282.5
253     VV(J,L)=0. 2283.
254     !
255     !CAS Applying latitudinal depedence of convection params
256     !
257     ! tprimrad = 3.14159*(float(j-1)/22.5)
258     ! TPRIMCO=6.0+(4.0*cos(tprimrad))
259     TPRIM2(J,L)=TPRIM0(j)/(EXPBYK(SIG(L)*P(1,J)+PTOP))**2.
260     !CAS TPRIM2(J,L)=2./(EXPBYK(SIG(L)*P(1,J)+PTOP))**2. 2283.5
261     ! AJL(J,L,59)=AJL(J,L,59)+2.*PSF
262     AJL(J,L,59)=AJL(J,L,59)+TPRIM0(j)*PSF
263     c TPRIM2(J,L)=0
264     PHIT(J,L)=PHI(1,J,L) 2284.
265     10 CONTINUE 2284.5
266     DO L=1,LM
267     PHITAV(L)=0.
268     CCC=0.
269     DO J=1,JM
270     PHITAV(L)=PHITAV(L)+PHIT(J,L)*COSP(J)
271     CCC=CCC+COSP(J)
272     ENDDO
273     PHITAV(L)=PHITAV(L)/CCC
274     ENDDO
275     C 2285.
276     C COMPUTE PRESSURE, P**KAPA, STATIC STABILITY AND SCALE HEIGHT 2285.5
277     C 2286.
278     DO 24 L=1,LZ 2286.5
279     DO 24 J=1,JM 2287.
280     PRESR(J,L+1)= P(1,J)*SIGE(L+1)+PTOP 2287.5
281     PSRK(J,L+1)=PRESR(J,L+1)**KAPA 2288.
282     IF(J.GE.2) PRSRU(J,L)=.5*(P(1,J)+P(1,J-1))*SIG(L)+PTOP 2288.5
283     IF(J.GE.2) PSROK(J,L)=PRSRU(J,L)**KAPA 2289.
284     TEM=( P(1,J)*SIG(L)+PTOP)**KAPA 2289.5
285     c24 TZA(J,L)= T(1,J,L)*TEM 2290.
286     24 TZA(J,L)= TZN(J,L)*TEM
287     DO 23 J=1,JM 2290.5
288     DO 22 L=1,LZM1 2291.
289     c TEM= T(1,J,L+1)- T(1,J,L) 2291.5
290     TEM= TZN(J,L+1)- TZN(J,L)
291     IF(TEM.LT..1) TEM=.1 2292.
292     DTHDP(J,L)=-TEM/(P(1,J)*DSIGO(L)) 2292.5
293     c DSHDP(J,L)=(Q(1,J,L)-Q(1,J,L+1))/(P(1,J)*DSIGO(L)) 2293.
294     DSHDP(J,L)=(QZN(J,L)-QZN(J,L+1))/(P(1,J)*DSIGO(L))
295     TEMJL(J,L+1)=-TEM/(SIG(L+1)-SIG(L))*RGAS* 2293.5
296     * PSRK(J,L+1)/PRESR(J,L+1)/ P(1,J) 2294.
297     TEMX=GRAV*PRESR(J,L+1)/(.5*RGAS*(TZA(J,L)+TZA(J,L+1))) 2294.5
298    
299     #if ( defined CPL_CHEM )
300     !
301     deltap(j,l)=1./(p(1,j)*dsigo(l))
302     dp2dz (j,l)=-temx
303     !
304     #endif
305    
306    
307     DTHDZ(J,L+1)=-DTHDP(J,L)*TEMX 2295.
308     DSHDZ(J,L+1)=-DSHDP(J,L)*TEMX 2295.5
309     22 CONTINUE 2296.
310     TEMJL(J,1)=(TEMJL(J,2)*(DSIG(1)+DSIG(2))-TEMJL(J,3)*DSIG(1))/ 2296.5
311     * DSIG(2) 2297.
312     IF(TEMJL(J,1).LE..1*TEMJL(J,2)) TEMJL(J,1)=.1*TEMJL(J,2) 2297.5
313     TEMJL(J,LZ+1)=TEMJL(J,LZ) 2298.
314     PRESR(J,1)= P(1,J)+PTOP 2298.5
315     PSRK(J,1)=PRESR(J,1)**KAPA 2299.
316     DTHDZ(J,1)=DTHDZ(J,2) 2299.5
317     DSHDZ(J,1)=DSHDZ(J,2) 2300.
318     DTHDZ(J,LZ+1)=DTHDZ(J,LZ) 2300.5
319     DSHDZ(J,LZ+1)=DSHDZ(J,LZ) 2301.
320     23 CONTINUE 2301.5
321     DO 27 L=1,LZ 2302.
322     DO 26 J=2,JMM1 2302.5
323     26 TEMZ(J,L)=0.25*(TEMJL(J,L)*2.+TEMJL(J+1,L)+TEMJL(J-1,L)) 2303.
324     TEMJL(1,L)=(2.*TEMJL(1,L)+TEMJL(2,L))/3. 2303.5
325     TEMJL(JM,L)=(2.*TEMJL(JM,L)+TEMJL(JMM1,L))/3. 2304.
326     DO 27 J=2,JMM1 2304.5
327     TEMJL(J,L)=TEMZ(J,L) 2305.
328     27 CONTINUE 2305.5
329     DO 21 L=1,LZM1 2306.
330     DO 21 J=1,JM 2306.5
331     TEM=2.*PRESR(J,L+1)*GRAV/RGAS/(TZA(J,L)+TZA(J,L+1)) 2307.
332     BVFJL(J,L)=TEM*SQRT(TEMJL(J,L+1)) 2307.5
333     C DRY
334     c BVFM(J,L)=BVFJL(J,L)*BVFJL(J,L)+2.*DSHDZ(J,L+1)*CLH*GRAV/ 2308.
335     c * (TZA(J,L)+TZA(J,L+1)) 2308.5
336     BVFM(J,L)=BVFJL(J,L)*BVFJL(J,L)
337     c RVFM is now (Nd)**2 not (Nm)**2
338     C DRY
339     SCLH(J,L)=RGAS*.5*(TZA(J,L)+TZA(J,L+1))/GRAV 2309.
340     21 CONTINUE 2309.5
341     C 2310.
342     2310.5
343     C COMPUTE DEL QS/DEL T 2311.
344     C 2311.5
345     DO 30 L=1,LZM1 2312.
346     DO 30 J=1,JM 2312.5
347     TEMZ(J,L)=0. 2313.
348     DQSDT1(J,L)=0. 2313.5
349     DO 25 I=1,IS 2314.
350     PR=(P(I,J)*SIG(L)+PTOP) 2314.5
351     TM=TZA(J,L) 2315.
352     TEM=QSAT(TM,PR,HLAT) 2315.5
353     SHSAT(J,L)=TEM 2316.
354     c TEM=DQDTX/TM/TM*Q(1,J,L) 2316.5
355     c TEM1=DQDTX*Q(1,J,L)/TM/TM 2317.
356     TEM=DQDTX/TM/TM*QZN(J,L)
357     TEM1=DQDTX*QZN(J,L)/TM/TM
358     TEMZ(J,L)=TEMZ(J,L)+TEM 2317.5
359     DQSDT1(J,L)=DQSDT1(J,L)+TEM1 2318.
360     25 CONTINUE 2318.5
361     DQSDT1(J,L)=DQSDT1(J,L)*CLH/FIS 2319.
362     30 DQSDT(J,L)=TEMZ(J,L)/FIS 2319.5
363     C 2320.
364     C COMPUTE MERIDIONAL TEMPERATURE AND MOISTURE GRADIENT 2320.5
365     C 2321.
366     DO 3140 L=1,LZ 2321.5
367     DO 3140 J=2,JM 2322.
368     c DTHDY(J,L)=( T(1,J,L)- T(1,J-1,L))/DYU(J) 2322.5
369     DTHDY(J,L)=( TZN(J,L)- TZN(J-1,L))/DYU(J)
370     3140 CONTINUE 2323.
371     DO 3141 L=1,LZ 2323.5
372     DO 3142 J=3,JMM1 2324.
373     3142 TEMZ(J,L)=.25*(2.*DTHDY(J,L)+DTHDY(J+1,L)+DTHDY(J-1,L)) 2324.5
374     DTHDY(2,L)=(2.*DTHDY(2,L)+DTHDY(3,L))/3. 2325.
375     DTHDY(JM,L)=(2.*DTHDY(JM,L)+DTHDY(JMM1,L))/3. 2325.5
376     DO 3141 J=3,JMM1 2326.
377     3141 DTHDY(J,L)=TEMZ(J,L) 2326.5
378     DO 3144 L=1,LZ 2327.
379     DO 3144 J=2,JM 2327.5
380     c DSHDY(J,L)=( Q(1,J,L)- Q(1,J-1,L))/DYU(J) 2328.
381     DSHDY(J,L)=( QZN(J,L)- QZN(J-1,L))/DYU(J)
382     3144 CONTINUE 2328.5
383     DO 3148 L=1,LZ 2329.
384     DO 3148 J=2,JM 2329.5
385     DSHDYS(J,L)=(SHSAT(J,L)-SHSAT(J-1,L))/DYU(J) 2330.
386     3148 CONTINUE 2330.5
387     C DO 3147 L=1,LZ 2331.
388     C DO 3146 J=3,JMM1 2331.5
389     C3146 TEMZ(J,L)=.25*(2.*DSHDY(J,L)+DSHDY(J+1,L)+DSHDY(J-1,L)) 2332.
390     C DSHDY(2,L)=(2.*DSHDY(2,L)+DSHDY(3,L))/3. 2332.5
391     C DSHDY(JM,L)=(2.*DSHDY(JM,L)+DSHDY(JMM1,L))/3. 2333.
392     C DO 3147 J=3,JMM1 2333.5
393     C3147 DSHDY(J,L)=TEMZ(J,L) 2334.
394     C 2334.5
395     C COMPUTE WAVE DEPTH 2335.
396     C 2335.5
397     DO 3145 J=2,JM 2336.5
398     B45=BETAV(J)
399     FCONT=FOV(J)
400     DO 3145 L=1,LZM2 2338.
401     TEM1=.5*(BVFJL(J,L)+BVFJL(J-1,L)) 2338.5
402     TEM=.5*(SCLH(J,L)+SCLH(J-1,L))*TEM1*TEM1*B45*.25* 2339.
403     c * (T(1,J,L)+T(1,J,L+1)+T(1,J-1,L+1)+ 2339.5
404     * (TZN(J,L)+TZN(J,L+1)+TZN(J-1,L+1)+
405     c * T(1,J-1,L))/(FCONT*GRAV*.5*(DTHDY(J,L)+DTHDY(J,L+1))+1.E-20) 2340.
406     * TZN(J-1,L))/(FCONT*GRAV*.5*(DTHDY(J,L)+DTHDY(J,L+1))+1.E-20)
407     c DWVJL(J,L)= .5*(SCLH(J,L)+SCLH(J-1,L))/ 2340.5
408     c * (1.48*ABS(TEM)+.48) 2341.
409     DWVJL(J,L)= .5*(SCLH(J,L)+SCLH(J-1,L))/
410     & (2./RCZZ*(1.+ABS(TEM))-1.)
411     IF(DWVJL(J,L).LT.DWVCRI) DWVJL(J,L)=DWVCRI 2341.5
412     3145 CONTINUE 2342.
413     C 2342.5
414     C DETERMINE THE LAT.FOR CUTTING OFF THE WAVE PROPAGATION 2343.
415     C 2343.5
416     DO 3150 J=2,JM 2345.
417     USUM=0. 2345.5
418     DO 3151 L=1,LM 2346.
419     c3151 USUM=USUM+U(1,J,L)*DSIG(L) 2346.5
420     3151 USUM=USUM+UZN(J,L)*DSIG(L)
421     3150 UZM(J)=USUM 2347.
422     JB=JEQ+1
423     JE=JB+JINT 2348.
424     YCUTN=0. 2348.5
425     DO 3155 J=JB,JE 2349.
426     JR=JB+JE-J 2349.5
427     3155 IF(UZM(JR).LE.0.) GO TO 3156 2350.
428     GO TO 3157 2350.5
429     3156 YCUTN=(UZM(JR)*VLAT(JR+1)-UZM(JR+1)*VLAT(JR))/
430     * (UZM(JR)-UZM(JR+1)-1.E-20)
431     3157 CONTINUE
432     JE=JEQ-1
433     JB=JE-JINT 2352.5
434     YCUTS=0. 2353.
435     DO 3158 JR=JB,JE 2353.5
436     3158 IF(UZM(JR).LE.0.) GO TO 3159 2354.
437     GO TO 3160 2354.5
438     3159 YCUTS=(UZM(JR)*VLAT(JR-1)-UZM(JR-1)*VLAT(JR))/
439     * (UZM(JR)-UZM(JR-1)-1.E-20)
440     3160 CONTINUE 2356.
441     c print *,'YCUTN=',180./TWOPI*YCUTN,' YCUTS=',180./TWOPI*YCUTS
442     C 2356.5
443     C 2356.5
444     C NEW WEIGHTING VARIABLES BY WAVE DEPTH
445     C
446     c print *,' DWAV0'
447     c print *,DWAV0
448     DO J=1,JM
449     DWV(J)=DWAV0(J)
450     ENDDO
451     DO II=1,NVAIT
452     DO J=1,JM
453     DO L=1,LZM1
454     TEMZ(J,L)=(TZA(J,L+1)-TZA(J,L))*GRAV/(PHI(1,J,L+1)-PHI(1,J,L))
455     ENDDO ! L
456     ENDDO ! J
457     CALL VWEIGHAV(GRAV,TEMZ,DTDZJ,PHI,PHIS,DWV,1,LZM2,1,JM ,PTOP,
458     * SCLH,P,SIGE,DSIGO,CKS,CKN,IO,IM,JM,LM)
459     DO J=2,JM
460     DO L=1,LZM2
461     TEMZ(J,L)=.5*(DTHDY(J,L)+DTHDY(J,L+1))
462     ENDDO ! L
463     ENDDO ! J
464     CALL VWEI1AV(GRAV,TEMZ ,DTDYJ,PHI,PHIS,DWV,1,LZM2,2,JM,PTOP,
465     * SCLH,P,SIGE,DSIGO,CKS,CKN,IO,IM,JM,LM)
466     DO L=1,LZM1
467     DO J=1,JM
468     TEMZ(J,L)=DTHDZ(J,L+1)
469     ENDDO ! L
470     ENDDO ! J
471     CALL VWEIGHAV(GRAV,TEMZ,DTHDZJ,PHI,PHIS,DWV,1,LZM2,1,JM ,PTOP,
472     * SCLH,P,SIGE,DSIGO,CKS,CKN,IO,IM,JM,LM)
473     c Vertical averaging of N**2
474     CALL VWEIGHAV(GRAV,BVFM,BVFMJ,PHI,PHIS,DWV,1,LZM2,1,JM ,PTOP,
475     * SCLH,P,SIGE,DSIGO,CKS,CKN,IO,IM,JM,LM)
476     c Nw=sqrt((N**2)w)
477     DO j=1,JM
478     BVF(j)=sqrt(BVFMJ(j))
479     ENDDO
480     c Nw=sqrt((N**2)w)
481     CALL VWEIGHAV(GRAV,SCLH,SHJ, PHI,PHIS,DWV,1,LZM2,1,JM ,PTOP,
482     * SCLH,P,SIGE,DSIGO,CKS,CKN,IO,IM,JM,LM)
483     DO L=1,LZM1
484     DO J=1,JM
485     c TEMZ(J,L)=.5*(T(1,J,L)+T(1,J,L+1))
486     TEMZ(J,L)=.5*(TZN(J,L)+TZN(J,L+1))
487     ENDDO ! L
488     ENDDO ! J
489     CALL VWEIGHAV(GRAV,TEMZ,THSZM, PHI,PHIS,DWV,1,LZM2,1,JM ,PTOP,
490     * SCLH,P,SIGE,DSIGO,CKS,CKN,IO,IM,JM,LM)
491     DO L=1,LZM2
492     DO J=2,JM
493     c TEMZ(J,L)=.5*(U(1,J,L)+U(1,J,L+1))
494     TEMZ(J,L)=.5*(UZN(J,L)+UZN(J,L+1))
495     ENDDO ! L
496     ENDDO ! J
497     if(HPRNT)then
498     J=JM-1
499     print *,'LZ=',LZ
500     print *,' Eddypa UZN before VWEI1AV'
501     print *,(UZN(13,L),L=1,LM)
502     print *,(UZN(14,L),L=1,LM)
503     print *,DWV(J+1),DWV(J-1),DWV(J)
504     endif
505     CALL VWEI1AV(GRAV,TEMZ,UZM,PHI,PHIS,DWV,1,LZM2,2,JM,
506     * PTOP,SCLH,P,SIGE,DSIGO,CKS,CKN,IO,IM,JM,LM)
507     if(HPRNT)then
508     print *,' Eddypa UZM after VWEI1AV'
509     J=14
510     print *,UZM(J+1),UZM(J-1),UZM(J)
511     endif
512     GAMAD=GRAV/CP
513     DO J=1,JM
514     DO L=1,LZM1
515     TEDGE=.5*(TZA(J,L+1)+TZA(J,L))
516     QSE=QSAT(TEDGE,PRESR(J,L+1),LHE)
517     TEMZ(J,L)=GAMAD*(1.+LHE*QSE/(RGAS*TEDGE))/
518     * (1.+.622*LHE*LHE*QSE/(CP*RGAS*TEDGE*TEDGE))
519     ENDDO ! L
520     ENDDO ! J
521     CALL VWEIGHAV(GRAV,TEMZ,GAMAJ,PHI,PHIS,DWV,1,LZM2,1,JM ,PTOP,
522     * SCLH,P,SIGE,DSIGO,CKS,CKN,IO,IM,JM,LM)
523     DO J=2,JM
524     B45=BETAV(J)
525     FCONT=FOV(J)
526     c
527     c BVFMJ=(N**2)w
528     TEM1=.5*(BVFMJ(J)+BVFMJ(J-1))
529     TEM=.5*(SHJ(J)+SHJ(J-1))*TEM1*B45* .5*(THSZM(J)+THSZM(J-1))/
530     * (FCONT*GRAV*DTDYJ(J)+1.E-20)
531     c
532     GAMA(J)=ABS(TEM)
533     HTJ(J)=.5*(SHJ(J)+SHJ(J-1))/(1.+GAMA(J))
534     IF(HTJ(J).LT.DWVCRI) HTJ(J)=DWVCRI
535     c DWV(J)=.5*(SHJ(J)+SHJ(J-1))/(1.48*GAMA(J)+.48)
536     DWV(J)=.5*(SHJ(J)+SHJ(J-1))/(2./RCZZ*(1.+GAMA(J))-1.)
537     IF(DWV(J).LT.DWVCRI) DWV(J)=DWVCRI
538     IF(BVFMJ(J).LT.1.E-5) BVFMJ(J)=1.E-5
539     C
540     DWV(J)=0.5*(DWV(j)+DWAV0(J))
541     C
542     ENDDO ! J
543     IF(BVFMJ(1).LT.1.E-5) BVFMJ(1)=1.E-5
544     c print *,' NEW II=',ii
545     c print *,' DWV'
546     c print *,DWV
547     AERR=0.0
548     RERR=0.0
549     AERRMAX=0.0
550     RERRMAX=0.0
551     DO J=2,JM
552     ADWVERR(J)=abs(DWV(J)-DWAV0(J))
553     RDWVERR(J)=ADWVERR(J)/DWAV0(J)
554     AERR=AERR+ADWVERR(J)
555     RERR=RERR+RDWVERR(J)
556     DWAV0(J)=DWV(J)
557     if(ADWVERR(J).gt.AERRMAX)then
558     AERRMAX=ADWVERR(J)
559     JAMAX=j
560     endif
561     if(RDWVERR(J).gt.RERRMAX)then
562     RERRMAX=RDWVERR(J)
563     JRMAX=j
564     endif
565     ENDDO
566     c print *,'Eddy iteration ii=',ii
567     c print *,'AERR=',AERR/(jm-1.0),' RERR=',RERR/(jm-1.0)
568     c print * ,'MAXAERR=',AERRMAX,' JAMAX=',JAMAX
569     c print * ,'MAXRERR=',RERRMAX,' JRMAX=',JRMAX
570     c print *,' DWV D'
571     c print *,DWV
572     c print *,' HTJ d'
573     c print *,HTJ
574     c print *,' SHJ H'
575     c print *,SHJ
576     if(RERRMAX.lt.1.e-3) go to 7788
577     ENDDO ! iterations
578     7788 continue
579     c if(ii.gt.20)print *,' NIT=',II
580     DO J=1,JM
581     DWAV0(J)=DWV(J)
582     ENDDO
583     c print *,' after DAAVO calculation'
584     c print *,' GAMA '
585     c print *,GAMA
586     c print *,' DWV'
587     c print *,DWV
588     c print *,' HTJ'
589     c print *,HTJ
590     C
591     C DWV (D) and CAMA on P-grid
592     do j=2,jmm1
593     TEM=SHJ(J)*BVFMJ(J)*BETAP(J)*THSZM(J)/
594     & (FO(J)*GRAV*0.5*(DTDYJ(J)+DTDYJ(J+1)))
595     GAMAP(j)=abs(TEM)
596     DWVP(J)=SHJ(J)/(2./RCZZ*(1.+GAMAP(J))-1.)
597     IF(DWVP(J).LT.DWVCRI) DWVP(J)=DWVCRI
598     HTJP(J)=SHJ(J)/(1.+GAMAP(J))
599     IF(HTJP(J).LT.DWVCRI) HTJP(J)=DWVCRI
600     enddo
601     c print *,' GAMAP '
602     c print *,(GAMAP(j),j=2,jm-1)
603     c print *,' '
604     c print *,(0.5*(GAMA(j)+GAMA(J+1)),j=2,jm-1)
605     c print *,' DWVP'
606     c print *,(DWVP(j),j=2,jm-1)
607     c print *,' '
608     c print *,(0.5*(DWV(j)+DWV(J+1)),j=2,jm-1)
609     C NEW
610     C 2387.
611     c Adjust wave length near pole
612     DO j=2,JM
613     if(j.ne.JEQ)then
614     TEM1=.5*(BVF(J)+BVF(J-1))
615     TEM=.5*(SHJ(J)+SHJ(J-1))*TEM1/ABS(FOV(J))
616     TEMH=TEMKGV(j)*TEM
617     TEMKV(J)=sqrt(((1.+GAMA(J))/RCZZ)**2-0.25)
618     C
619     if(TEMKV(J).lt.TEMH)then
620     DWVA(J)=.5*(SHJ(J)+SHJ(J-1))/(sqrt(4.*TEMH**2+1.)-1.)
621     else
622     DWVA(J)=DWV(J)
623     endif
624     endif
625     ENDDO ! J
626     do j=2,jmm1
627     TEM=BVF(J)*SHJ(J)/abs(FO(J))
628     TEMH=TEMKGP(j)*TEM
629     TEMKP(J)=sqrt(((1.+GAMAP(J))/RCZZ)**2-0.25)
630     C
631     if(TEMKP(J).lt.TEMH)then
632     DWVPA(J)=SHJ(J)/(sqrt(4.*TEMH**2+1.)-1.)
633     else
634     DWVPA(J)=DWVP(J)
635     endif
636     enddo
637     c print *,' After adjustment'
638     c print *,' DWV'
639     c print *,DWV
640     c print *,' DWVA'
641     c print *,DWVA
642     c print *,' DWVP'
643     c print *,DWVP
644     c print *,' DWVPA'
645     c print *,DWVPA
646     c print *,' HTJ'
647     c print *,HTJ
648     c stop
649     C COMPUTE TEMPERATURE VARIANCES 2387.5
650     C 2388.
651     DO 4013 J=2,JMM1 2388.5
652     COEKC=CKS 2389.
653     IF(J.GT.JHALF) COEKC=CKN 2389.5
654     c TEM=BVF(J)*.25*(HTJ(J)+HTJ(J+1))*(DTDYJ(J)+DTDYJ(J+1))/ 2390.
655     TEM=BVF(J)*.5*HTJP(J)*(DTDYJ(J)+DTDYJ(J+1))/
656     * (FO(J)+1.E-20) 2390.5
657     IF ((DTDYJ(J)+DTDYJ(J+1))*FO(J).GT.0.) TEM=0. 2390.6
658     c IF(FO(J).EQ.0.) TEM=0. 2391.
659     if(abs(FO(J)).lt.1.e-15) TEM=0.
660     PLAND=FDATA(1,J,2) 2395.
661     PWATER=1.-PLAND 2395.5
662     PLICE=FDATA(1,J,3)*PLAND 2396.
663     PEARTH=PLAND-PLICE 2396.5
664     POICE=ODATA(1,J,2)*PWATER 2397.
665     POCEAN=PWATER-POICE 2397.5
666     if(POCEAN.LE.1.E-5)then
667     POCEAN=0.
668     POICE=PWATER
669     endif
670     tland=GDATA(1,J,4)
671     toice=GDATA(1,J,3)
672     tlice=GDATA(1,J,13)
673     tocean=ODATA(1,J,1)
674     TSS=POCEAN*tocean+POICE*toice+PLICE*tlice+PEARTH*tland
675     TPR1=POCEAN*(tocean-TSS)**2
676     TPR2=POICE*(toice-TSS)**2
677     TPR3=PLICE*(tlice-TSS)**2
678     TPR4=PEARTH*(tland-TSS)**2
679     tprmg(J)=tprmg(J)+(TPR1+TPR2+TPR3+TPR4)
680     ntprmg(J)=ntprmg(J)+1
681     c TEMH=.5*(DWV(J)+DWV(J+1)) 2400.5
682     TEMH=DWVPA(J)
683     DO 4013 L=1,LZ 2400.
684     ZPHI=(PHI(1,J,L)-PHIS(1,J))/(GRAV*TEMH+1.E-20) 2401.5
685     C TEM1=TEM*(1.-.5*ZPHI) 2402.
686     TEM1=TEM*(1.-.5*ZPHI)*(1.+DQSDT1(J,L)) 2402.5
687     C for analisys
688     ! TEMAN=TEM*(1.-.5*ZPHI)
689     C for analisys
690     TPRIM2(J,L)=TPRIM2(J,L)+1.*TEM1*TEM1*EXP(-ZPHI/COEKC)/COEKC 2403.
691     ! 2319
692     ! TPRIM2(J,L)=max(TPRIM2(J,L),1.*TEM1*TEM1*EXP(-ZPHI/COEKC)/COEKC)
693     PK2=(EXPBYK(SIG(L)*P(1,J)+PTOP))**2. 2403.5
694     IF(L.EQ.1) THEN
695     TPRIM2(J,L)=TPRIM2(J,L)+.6065* 2404.
696     & (TPR1+TPR2+TPR3+TPR4)/PK2 2404.5
697     AJL(J,L,59)=AJL(J,L,59)+.6065*
698     & (TPR1+TPR2+TPR3+TPR4)*P(1,J)
699     ELSEIF(L.EQ.2) THEN
700     TPRIM2(J,L)=TPRIM2(J,L)+.1653* 2405.
701     & (TPR1+TPR2+TPR3+TPR4)/PK2 2405.5
702     AJL(J,L,59)=AJL(J,L,59)+.1653*
703     & (TPR1+TPR2+TPR3+TPR4)*P(1,J)
704     ENDIF
705     AJL(J,L,59)=AJL(J,L,59)+TPRIM2(J,L)
706     4013 CONTINUE 2406.
707     tprmg(JM)=0.
708     tprmg(1)=0.
709     ntprmg(JM)=1.
710     ntprmg(1)=1.
711     CCC ************
712     C 2406.5
713     C COMPUTE WIND VARIANCE AND EDDY TRANSFER COEFFICIENTS 2407.
714     C 2407.5
715     DO 4015 I=1,IM 2408.
716     DO 4015 J=2,JM 2408.5
717     COEKD=CKS 2409.
718     IF(J.GT.JHALF) COEKD=CKN 2409.5
719     TEM2=FOV(J)
720     TEM3=.5*(THSZM(J)+THSZM(J-1)) 2410.5
721     TEM=HTJ(J)*DTDYJ(J)*GRAV/(TEM2*TEM3+1.E-20) 2411.
722     IF (DTDYJ(J)*TEM2.GT.0.) TEM=0. 2411.1
723     c IF(TEM2.EQ.0.) TEM=0. 2411.5
724     if(abs(TEM2).lt.1.e-15) TEM=0.
725     C TEML=3.*1.4142*.5*(BVF(J)+BVF(J-1))*DWV(J)/(TEM2+1.E-20) 2412.
726     C TEML=ABS(TEML) 2412.5
727     C TEMY=RADIUS*(.25*TWOPI-.5*ABS(LAT(J)+LAT(J-1))) 2413.
728     C RATIO=TEMY/(TEML+1.E-20) 2413.5
729     C IF(TEML.GT.TEMY) TEM=TEM*RATIO 2414.
730     UVAR(J)=.5*TEM*TEM/COEKD 2414.5
731     WMGETEMP=.33*TEM*TEM/COEKD
732     c WMGE(I,J)=max(WMGETEMP,WMGMIN) 2415.
733     WMGE(I,J)=WMGETEMP+WMGMIN 2415.1
734     4015 continue
735     IF(DMOD(TAU,25.).EQ.0..AND.IABS(MRCH).EQ.1) 2415.5
736     * WRITE(6,905)(WMGE(1,J),J=2,JM) 2416.
737     905 FORMAT(' WMGE=',24F5.1) 2416.5
738     C 2417.
739     DO 2230 J=2,JM 2417.5
740     COEKC=CKS 2418.
741     IF(J.GT.JHALF) COEKC=CKN 2418.5
742     TEM=CONSTA*(BVF(J)+BVF(J-1))*ABS(DTDYJ(J))*HTJ(J)*HTJ(J) 2419.
743     * *GRAV/(THSZM(J)+THSZM(J-1)) 2419.5
744     COEKY(J)=0. 2420.
745     FVTH=FOV(J)
746     c IF(FVTH.NE.0.) COEKY(J)=TEM/(FVTH*FVTH)/COEKC 2421.
747     if(abs(FVTH).gt.1.e-15)COEKY(J)=TEM/(FVTH*FVTH)/COEKC
748     IF(DTDYJ(J)*FVTH.GT.0.) COEKY(J)=0. 2421.1
749     c TEMH=DWV(J) 2425.5
750     TEMH=DWVA(J) 2425.5
751     c
752     C Old restriction
753     c
754     c TEML=3.*1.4142*(BVF(J)+BVF(J-1))*DWV(J)/((FO(J)+FO(J-1))+1.E-20) 2421.5
755     c TEML=ABS(TEML) 2422.
756     c TEMY=RADIUS*(.25*TWOPI-.5*ABS(LAT(J)+LAT(J-1))) 2422.5
757     c TEMK=.7071*.5*TWOPI*.25*(BVF(J)+BVF(J-1))* 2423.
758     c * (SHJ(J)+SHJ(J-1))/(TEMY*ABS(FVTH)+1.E-20) 2423.5
759     c if(TEML.GT.TEMY) then
760     c if(TEMK.lt.1.e19)then
761     c TEMH=.5*(SHJ(J)+SHJ(J-1))/(SQRT(4.*TEMK*TEMK+1.)-1.)
762     c else
763     c TEMH=.5*(SHJ(J)+SHJ(J-1))/(2.*TEMK-1.)
764     c end if
765     c end if
766     c
767     C Old restriction
768     c
769    
770     DO 2230 L=1,LZM1 2425.
771     ZSTAR(J,L)=.5*(PHI(1,J,L)+PHI(1,J-1,L)-PHIS(1,J)-PHIS(1,J-1)) 2427.
772     * /(GRAV*TEMH+1.E-20) 2427.5
773     EXPFAC=EXP(-ZSTAR(J,L)/COEKC) 2428.
774     C TEMV=5.+45.*(5.*SIG(L)*EXP(1.-5.*SIG(L)))**4. 2428.5
775     C * +TEMV 2429.5
776     COEJL(J,L)=COEKY(J)*2. 2430.
777     COEJL1(J,L)=COEKY(J)*EXPFAC 2430.5
778     2230 CONTINUE 2431.
779     C 2431.5
780     C PARA. MERIDIONAL EDDY MOMENTUM TRANSPORT 2432.
781     C 2432.5
782     DO 35 J=2,JM 2436.
783     BETA=BETAV(J)
784     FAVE=FOV(J)
785     C DRY
786     TEM1=.5*(BVFMJ(J)+BVFMJ(J-1)) 2437.5
787     c Nd instead of Nm Eq 23 from Part II.
788     c See change in calculation of BVFM above
789     c BVFMJ=(Nd)**2, while BVF=Nd
790     C DRY
791     TEM=FAVE*GRAV*DTDYJ(J)*2./(1. *TEM1*(THSZM(J)+THSZM(J-1)) 2438.
792     * *DWV(J)+1.E-20) 2438.5
793     if(HPRNT)then
794     if(J.eq.16)then
795     print *,' Eddypa 2350'
796     print *,FAVE,DTDYJ(J),TEM1
797     print *,THSZM(J),THSZM(J-1),DWV(J)
798     print *,j,TEM
799     endif
800     endif
801     IF (J.GT.2.AND.J.LT.JM) 2439.
802     * DUDY2=(UZM(J+1)+UZM(J-1)-2.*UZM(J))/DYP(J)/DYP(J) 2439.5
803     * -(SINL(J)+SINL(J-1))/(COSL(J)+COSL(J-1))*.5* 2440.
804     * (UZM(J+1)-UZM(J-1))/(RADIUS*DYU(J)) 2440.5
805     * -UZM(J)/(RADIUS*.5*(COSL(J)+COSL(J-1)))**2. 2441.
806     IF (J.EQ.2) DUDY2=(UZM(3)-2.*UZM(2))/DYP(2)**2. 2441.5
807     * -(SINL(2)+SINL(1))/(COSL(2)+COSL(1))* 2442.
808     * (UZM(3)-UZM(2))/(RADIUS*DYP(2)) 2442.5
809     * -UZM(2)/(RADIUS*.5*(COSL(2)+COSL(1)))**2. 2443.
810     IF (J.EQ.JM) DUDY2=(UZM(JMM1)-2.*UZM(JM))/DYP(JM)**2. 2443.5
811     * -(SINL(JMM1)+SINL(JM))/(COSL(JMM1)+COSL(JM))* 2444.
812     * (UZM(JM)-UZM(JMM1))/(RADIUS*DYP(JMM1)) 2444.5
813     * -UZM(JM)/(RADIUS*.5*(COSL(JMM1)+COSL(JM)))**2. 2445.
814     if(HPRNT)then
815     if(J.eq.14)then
816     print *,' Eddypa 35 J=',J
817     print *,BETAV(J),DUDY2
818     print *,UZM(J+1),UZM(J-1),UZM(J)
819     endif
820     endif
821     DO 35 L=1,LZM1 2445.5
822     DQDY(J,L)= BETAV(J)-DUDY2
823     C DRY
824     c
825     c35 DQDY1(J,L)=TEM*(1.+.5*(DQSDT1(J,L)+DQSDT1(J-1,L)))+DQDY(J,L) 2446.5
826     c
827     c M=0. Eq 23 from Part II.
828     c
829     DQDY1(J,L)=TEM+DQDY(J,L)
830     CCCCC check + or *
831     c
832     C DRY
833     35 CONTINUE
834     DO 2350 J=2,JM 2447.
835     TEM=0.25*(COSL(J)+COSL(J-1))**2. 2447.5
836     DO 2380 L=1,LZM1 2448.
837     TEMZ(J,L)=DQDY(J,L)*TEM 2448.5
838     DQDY1(J,L)=COEJL1(J,L)*DQDY1(J,L)*TEM 2449.
839     if(HPRNT)then
840     if(J.eq.14)then
841     print *,' Eddypa 2350'
842     print *,L,TEM
843     print *,DQDY(J,L),COEJL1(J,L),DQDY1(J,L)
844     endif
845     endif
846     2380 CONTINUE 2449.5
847     2350 CONTINUE 2450.
848     C**** TAKING VERTICAL AVERAGES 2450.5
849     DO 2365 J=2,JM 2451.
850     SUMSIG=0. 2451.5
851     TEMJ(J)=0. 2452.
852     TSZM(J)=0. 2452.5
853     DO 2366 L=LBGN,LZM1 2453.
854     SUMSIG=SUMSIG+DSIG(L) 2453.5
855     TSZM(J)=TSZM(J)+DQDY1(J,L)*DSIG(L) 2454.
856     2366 TEMJ(J)=TEMJ(J)+TEMZ(J,L)*DSIG(L) 2454.5
857     TSZM(J)=TSZM(J)/SUMSIG 2455.
858     TEMJ(J) = TEMJ(J)/SUMSIG 2455.5
859     TSZM(J)=TSZM(J)*MOMCOEF
860     2365 CONTINUE 2456.
861     if(HPRNT)then
862     print *,' Eddypa 2365'
863     print *,'TEMJ'
864     print *,TEMJ
865     print *,'TSZM'
866     print *,TSZM
867     endif
868     C**** CALCULATE THE ADJUSTMENT COEFFICIENT FOR MOMENTUM BALANCING 2456.5
869     PFRAC=.0 2457.
870     YSCALE=400000./RADIUS 2457.5
871     DO 2355 JH=1,2 2458.5
872     JB=2 2459.
873     JE=JM 2459.5
874     IF(JH.EQ.1) JE=JEQ-1 2460.
875     IF(JH.EQ.2) JB=JEQ+1 2460.5
876     YZERO=YCUTN 2461.
877     IF(JH.EQ.1) YZERO=YCUTS 2461.5
878     TEM1=0. 2462.
879     TEM2=0. 2462.5
880     DO 2353 J=JB,JE 2463.
881     FKB=1. 2463.5
882     DELTY=VLAT(J)-YZERO 2464.
883     c NLAT=ABS(DELTY/DLAT)+1. 2464.5
884     IF(JH.EQ.1.AND.DELTY.GT.0.) 2465.
885     * FKB=EXP(-RADIUS*DELTY/DSCALE)
886     IF(JH.EQ.2.AND.DELTY.LT.0.) 2466.
887     * FKB=EXP(RADIUS*DELTY/DSCALE)
888     SP=.5*(P(1,J)+P(1,J-1)) 2467.
889     TEM1=TEM1+TEMJ(J)*SP*FKB 2467.5
890     TEM2=TEM2+TSZM(J)*SP 2468.
891     2353 CONTINUE 2468.5
892     DO 2354 J=JB,JE 2469.
893     SP=.5*(P(1,J)+P(1,J-1)) 2469.5
894     2354 TSZM(J)=TSZM(J)-PFRAC*TEM2/(SP*(JE-JB+1)) 2470.
895     COEKB=-TEM2*(1.-PFRAC)/(TEM1+1.E-20) 2470.5
896     C IF(JH.EQ.1) CKS=1./COEKB 2471.
897     C IF(JH.EQ.2) CKN=1./COEKB 2471.5
898     IF(DMOD(TAU,25.).EQ.0..AND.IABS(MRCH).EQ.1) 2472.
899     * WRITE(6,909) JH,COEKB,TAU 2472.5
900     DO 2358 J=JB,JE 2473.
901     FKB=1. 2473.5
902     DELTY=VLAT(J)-YZERO 2474.
903     c NLAT=ABS(DELTY/DLAT)+1. 2474.5
904     IF(JH.EQ.1.AND.DELTY.GT.0.) 2475.
905     * FKB=EXP(-RADIUS*DELTY/DSCALE)
906     IF(JH.EQ.2.AND.DELTY.LT.0.) 2476.
907     * FKB=EXP(RADIUS*DELTY/DSCALE)
908     CONKB=COEKB*FKB 2477.
909     2358 TEMJ(J)=DYU(J)*(CONKB*TEMJ(J)+1.000*TSZM(J)) 2477.5
910     2355 CONTINUE 2478.
911     909 FORMAT(10X,'JH =',I5,8X,'KB =',E16.5,5X,'TAU =',F10.2) 2478.5
912     C TOTLJ=JM-3 2479.
913     DO 2451 J=1,JM 2479.5
914     2451 FD(1,J)=PT(1,J)*DXYP(J) 2480.
915     FD(1,1)=2.*FD(1,1) 2480.5
916     FD(1,JM)=2.*FD(1,JM) 2481.
917     DO 2452 J=2,JM 2481.5
918     RFDU=2./(FD(1,J)+FD(1,J-1)) 2482.
919     DO 2452 L=1,LZM1 2482.5
920     2452 UOLD(J,L)=UT(1,J,L)*RFDU 2483.
921     c
922     if(uvtr)then
923     c
924     C**** CALCULATE U'V' CONTRIBUTION IN U-MOMENTUM EQUATION 2483.5
925     C TDISG=1.-SIGE(LZ) 2484.
926     C TEM2=.36-SIGE(LZ)*SIGE(LZ) 2484.5
927     C TEM3=.216-SIGE(LZ)**3. 2485.
928     C XCONT=TDISG/(TEM2*.3-TEM3/3.) 2485.5
929     if(OLDUVVERT)then
930     C Old
931     VVSUM=0. 2486.
932     DO 2384 L=1,LM 2486.5
933     2384 VVSUM=VVSUM+VVN(L)+VVS(L) 2487.
934     FSUM=0. 2487.5
935     DO 2383 L=1,LM 2488.
936     2383 FSUM=FSUM+(VVN(L)+VVS(L))*DSIG(L)/VVSUM 2488.5
937     else
938     C New
939     VVNEWAV=0.
940     ZZD=250.
941     CALL VFUNCT(VVD,ZZD)
942     DO L=1,LM
943     ZZZ=PHITAV(L)/GRAV+250.
944     if(L.eq.LM)then
945     ZZU=30000.
946     else
947     ZZU=0.5*(PHITAV(L)+PHITAV(L+1))/GRAV+250.
948     endif
949     CALL VFUNCT(VVNEW,ZZZ)
950     CALL VFUNCT(VVU,ZZU)
951     VVNEWA(L)=0.25*(VVD+2.*VVNEW+VVU)
952     VVNEWAV=VVNEWAV+VVNEWA(L)*DSIG(L)
953     VVD=VVU
954     ZZD=ZZU
955     ENDDO
956     if(abs(1.-VVNEWAV).gt.1.e-1)then
957     print *,' VVNEWAV=',VVNEWAV
958     endif
959     endif ! OLDUVVERT
960     DO 2385 L=1,LZ 2489.
961     if(OLDUVVERT)then
962     TEM=(VVN(L)+VVS(L))/(VVSUM*FSUM) 2489.5
963     C IF(SIG(L).LT.0.6) TEM=XCONT*SIG(L)*(0.6-SIG(L)) 2490.
964     else
965     TEM=VVNEWA(L)/VVNEWAV
966     endif ! OLDUVVERT
967     C New Aexp 049.98
968     JB=2
969     JE=JMM1
970     DO 2360 J=JB,JE
971     VU(J,L)=VU(J-1,L)*P(1,J-1)*(COSL(J-1)/COSL(J))**2.+
972     * TEMJ(J)*TEM*.5*(P(1,J)+P(1,J-1))/COSL(J)**2.
973     if(HPRNT)then
974     c if(J.eq.JPR-1.and.L.eq.4)then
975     if(L.eq.4)then
976     print *,' Eddypa 2360'
977     print *,VU(J-1,L),TEMJ(J),P(1,J-1)
978     c print *,TEMJ(J),TEM,P(1,J)
979     endif
980     endif
981     2360 VU(J,L)=VU(J,L)/P(1,J)
982     C New
983     DO 2370 J=JB,JE 2495.
984     DO 2370 I=1,IX 2495.5
985     if(HPRNT)then
986     if(J.eq.JPR-1.and.L.eq.4)then
987     print *,' Eddypa 2370'
988     print *,DTE,DXP(J),VU(J,L),P(1,J)
989     print *,FLUX,UT(I,J+1,L)
990     endif
991     endif
992     FLUX=DTE*VU(J,L)*P (1,J)*DXP(J) 2496.
993     UT(I,J+1,L)=UT(I,J+1,L)+FLUX 2496.5
994     2370 UT(I,J,L)=UT(I,J,L)-FLUX 2497.
995     2385 CONTINUE 2497.5
996     c
997     if(HPRNT)then
998     print *,' Eddypa 2385'
999     print *,'UT'
1000     print *,(UT(1,JPR,L),L=1,LM)
1001     endif ! PRNT
1002     endif
1003     c end uvtr
1004     C 2498.
1005     if(wvtr)then
1006     c
1007     C PRESCRIBE W'U' BASED ON THE SEMI-SPECTRAL MODEL 2498.5
1008     C 2499.
1009     C**** ADJUSTMENT OF U TO MAINTAIN TOTAL KINETIC ENERGY 2499.5
1010     C 2500.
1011     DO 2454 J=2,JM 2500.5
1012     RFDU=2./(FD(1,J)+FD(1,J-1)) 2501.
1013     DO 2454 L=1,LZM1 2501.5
1014     2454 UNEW(J,L)=UT(1,J,L)*RFDU 2502.
1015     DO 2456 L=1,LZM1 2502.5
1016     DO 2456 JH=1,2 2503.
1017     JB=2 2503.5
1018     JE=JM 2504.
1019     IF(JH.EQ.1) JE=JEQM1 2504.5
1020     IF(JH.EQ.2) JB=JEQP1 2505.
1021     SUM1=0. 2505.5
1022     SUM2=0. 2506.
1023     DO 2458 J=JB,JE 2506.5
1024     TEM=0.5*(FD(1,J)+FD(1,J-1)) 2507.
1025     SUM1=SUM1+TEM*UOLD(J,L)*UOLD(J,L) 2507.5
1026     2458 SUM2=SUM2+TEM*UNEW(J,L)*UNEW(J,L) 2508.
1027     DIFF=(SUM2-SUM1)/SUM2 2508.5
1028     DO 2460 J=JB,JE 2509.
1029     TEM=0.5*(FD(1,J)+FD(1,J-1)) 2509.5
1030     TEM1=TEM*UNEW(J,L)*UNEW(J,L)*(1.-DIFF) 2510.
1031     IF(TEM1.LE.0.) GO TO 2465 2510.5
1032     TEM2=SQRT(TEM1/TEM) 2511.
1033     IF(UNEW(J,L).GE.0.) UNEW(J,L)=TEM2 2511.5
1034     IF(UNEW(J,L).LT.0.) UNEW(J,L)=-TEM2 2512.
1035     GO TO 2460 2512.5
1036     2465 UNEW(J,L)=0. 2513.
1037     2460 CONTINUE 2513.5
1038     2456 CONTINUE 2514.
1039     DO 2462 J=2,JM 2514.5
1040     FDU=0.5*(FD(1,J)+FD(1,J-1)) 2515.
1041     DO 2462 L=1,LZM1 2515.5
1042     2462 UT(1,J,L)=UNEW(J,L)*FDU 2516.
1043     if(HPRNT)then
1044     print *,' Eddypa 2462'
1045     print *,'UT'
1046     print *,(UT(1,JPR,L),L=1,LM)
1047     endif ! PRNT
1048     endif
1049     c end wvtr
1050     C 2516.5
1051     if(uvtr)then
1052     C**** CALCULATE OTHER U'V' AND V'V' CONTRIBUTIONS IN MOMENTUM EQUATIONS 2517.
1053     C 2517.5
1054     DO 2490 L=1,LZM1 2518.
1055     DO 2490 J=2,JMM1 2518.5
1056     FLUX=DTE*P(1,J)*DXP(J)*.5*(VV(J,L)+VV(J+1,L)) 2519.
1057     VT(1,J+1,L)=VT(1,J+1,L)+FLUX 2519.5
1058     2490 VT(1,J,L)=VT(1,J,L)-FLUX 2520.
1059     DO 2492 L=1,LZM1 2520.5
1060     DO 2494 J=2,JM 2521.
1061     TEM=DTE*(DXP(J)-DXP(J-1))*(P(1,J)+P(1,J-1))*.5 2521.5
1062     FLUXU=.5*TEM*(VU(J,L)+VU(J-1,L)) 2522.
1063     FLUXV=TEM*VV(J,L) 2522.5
1064     UT(1,J,L)=UT(1,J,L)-FLUXU 2523.
1065     2494 VT(1,J,L)=VT(1,J,L)+FLUXV 2523.5
1066     2492 CONTINUE 2524.
1067     if(HPRNT)then
1068     print *,' Eddypa 2492'
1069     print *,'UT'
1070     print *,(UT(1,JPR,L),L=1,LM)
1071     endif ! PRNT
1072     endif
1073     c end uvtr
1074     C 2524.5
1075     if(uttr)then
1076     C PARAMETERIZATION OF MERIDIONAL EDDY HEAT TRANSPORTS 2525.
1077     C 2525.5
1078     COEKC=1. 2526.
1079     DO 2245 L=1,LZ 2526.5
1080     DO 2240 J=2,JM 2527.
1081     COEKD=CKS 2527.5
1082     IF(J.GT.JHALF) COEKD=CKN 2528.
1083     TEM=.5*(PHI(1,J,L)+PHI(1,J-1,L)-PHIS(1,J)-PHIS(1,J-1)) 2528.5
1084     * /GRAV/(DWV(J)+1.E-20) 2529.
1085     EXPFAC=EXP(-TEM/COEKD) 2529.5
1086     COEJL(J,L)=COEJL(J,L)*EXPFAC 2530.
1087     IF(L.EQ.1) COEJL(J,L)=COEJL(J,L)*.3935 2530.5
1088     IF(L.EQ.2) COEJL(J,L)=COEJL(J,L)*.8347 2531.
1089     VTH(J,L)=-.5*COEJL(J,L)*DTDYJ(J) 2531.5
1090     #if ( defined CPL_CHEM )
1091     !
1092     ! --- Get the Kt
1093     !
1094     ! xxx =(T(1,j,l)-T(1,j-1,l))/dyv(j)
1095     xxx =dthdy(j,l)
1096     fkt(j,l)=-vth(j,l)/(xxx+1.e-20)
1097     c fkt(j,l)=dmax1(0.0,fkt(j,l))
1098     c
1099     c 012896: increase midlatitude eddy
1100     c here 0.0 can be changed to minimum
1101     c diffussion constant too
1102     c
1103     c fkt(j,l)=max(2.e6,fkt(j,l)*1.5)
1104     c fkt(j,l)=max(0.0,fkt(j,l)*1.5)
1105     c
1106     c 111596: climate model change vth to 1.5 time large already:
1107     c
1108     fkt(j,l)=max(0.0,fkt(j,l))
1109     !
1110     #endif
1111    
1112     FLUXT=DTE*VTH(J,L)*( P(1,J)+ P(1,J-1))*DXU(J)*.5 2532.
1113     TT(1,J,L)=TT(1,J,L)+FLUXT 2532.5
1114     TT(1,J-1,L)=TT(1,J-1,L)-FLUXT 2533.
1115     2240 CONTINUE 2533.5
1116     2245 CONTINUE 2534.
1117     endif
1118     c end uttr
1119     C
1120     C 2534.5
1121     if(uqtr)then
1122     C PARAMETERIZATION OF MERIDIONAL EDDY MOISTURE TRANSPORTS 2535.
1123     C 2535.5
1124     COEKC=1. 2536.
1125     DO 2345 L=1,LZ 2536.5
1126     DO 2340 J=2,JM 2537.
1127     IF(MRCH.EQ.0 ) GO TO 2340 2537.5
1128     VQ(J,L)=VTH(J,L)*P1000*.5*(DQSDT(J,L)+DQSDT(J-1,L)) 2538.
1129     IF(DSHDY(J,L)*DTHDY(J,L).LT.0.) VQ(J,L)=0. 2538.1
1130     FLUXSH=VQ(J,L)*DTE*.5*(P(1,J)+P(1,J-1))*DXV(J) 2538.5
1131     IF(FLUXSH.GT..5* QT(1,J-1,L)) FLUXSH=.5* QT(1,J-1,L) 2539.
1132     IF(FLUXSH.LT.-.5* QT(1,J,L)) FLUXSH=-.5* QT(1,J,L) 2539.5
1133     QT(1,J,L)= QT(1,J,L)+FLUXSH 2540.
1134     QT(1,J-1,L)= QT(1,J-1,L)-FLUXSH 2540.5
1135     2340 CONTINUE 2541.
1136     2345 CONTINUE 2541.5
1137     endif
1138     c end uqtr
1139     C 2542.
1140     if(wttr) then
1141     C PARAMETERIZE W'T' 2542.5
1142     C 2543.
1143     COEF=.25 2543.5
1144     DO 4150 J=2,JMM1 2544.
1145     COEKD=1. 2544.5
1146     C IF(J.GT.JHALF) COEKD=CKN 2545.
1147     FCONT=ABS(FO(J)) 2545.5
1148     c TEM2=.5*(HTJ(J)+HTJ(J+1)) 2546.
1149     TEM2=HTJP(J)
1150     TEM3=ABS(.5*(DTDYJ(J)+DTDYJ(J+1))) 2546.5
1151     TEM=TEM3**3.*TEM2*TEM2*GRAV*GRAV/ 2547.
1152     * (BVF(J)*THSZM(J)*THSZM(J)*FCONT*FCONT+1.E-20) 2547.5
1153     IF ((DTDYJ(J)+DTDYJ(J+1))*FO(J).GT.0.) TEM=0. 2547.6
1154     RFAC=(GAMAJ(J)/GAMAD)*(DTDZJ(J)+GAMAJ(J))/(DTDZJ(J)+GAMAD) 2551.
1155     IF(RFAC.LT.0.) RFAC=0. 2551.01
1156     RFAC=RFAC+(TEM3*GRAV/(THSZM(J)*FCONT*BVF(J)+1.E-20))**2 2551.1
1157     IF(RFAC.LT.0.01) RFAC=0.01 2551.5
1158     ALAM=.573*SQRT(RFAC)/(1.-.427*RFAC**.302) 2552.
1159     c TEMH=DWV(J)+DWV(J+1) 2553.
1160     TEMH=DWVPA(J)
1161     DO 4150 L=1,LZM2 2552.5
1162     TEMP=(PHI(1,J,L)+PHI(1,J,L+1)-2.*PHIS(1,J))/ 2554.
1163     * (GRAV*TEMH+1.E-20) 2554.5
1164     WTH(J,L)=-CONSTA*TEM*EXP(-TEMP/COEKD)*(P(1,J)*SIGE(L+1)+PTOP) 2555.
1165     * /SCLH(J,L)/P(1,J)*TEMP*(1.-.25*TEMP)/COEKD 2555.5
1166     WTH(J,L)=WTH(J,L)*(1.+ALAM/RFAC)/(1.+ALAM) 2556.
1167     C IF(WTH(J,L).GT.0.) WTH(J,L)=0. 2556.5
1168     4150 CONTINUE 2557.
1169     DO 4160 J=2,JMM1 2557.5
1170     SP= P(1,J) 2558.
1171     DO 4160 I=1,IS 2558.5
1172     FLUXDT=0. 2559.
1173     DO 4155 L=1,LZM2 2559.5
1174     FLUXT=DTE*SP*DXYP(J)*WTH(J,L) 2560.
1175     TT(I,J,L)=TT(I,J,L)+(FLUXT-FLUXDT)/DSIG(L) 2560.5
1176     FLUXDT=FLUXT 2561.
1177     4155 CONTINUE 2561.5
1178     TT(I,J,LZM1)=TT(I,J,LZM1)-FLUXDT/DSIG(LZM1) 2562.
1179     4160 CONTINUE 2562.5
1180     endif
1181     c end wttr
1182     C
1183     C 2563.
1184     if(wqtr) then
1185     C PARAMETERIZE W'Q' 2563.5
1186     C 2564.
1187     IF(MRCH.ne.0) then
1188     DO 3170 J=2,JMM1 2564.5
1189     SP=P(1,J) 2565.
1190     BETA=BETAP(J)
1191     FCONT=ABS(FO(J)) 2565.6
1192     TEM3=ABS(.5*(DTDYJ(J)+DTDYJ(J+1))) 2565.7
1193     ccc GAMAP=.5*(GAMA(J)+GAMA(J+1)) 2566.
1194     C FUNK=.55*(1.+GAMAP)*(1.+GAMAP)-.25 2566.5
1195     c FUND=SHJ(J)/(1.48*GAMAP+.48) 2567.
1196     FUND=SHJ(J)/(2./RCZZ*(1.+GAMAP(j))-1.)
1197     RFAC=(GAMAJ(J)/GAMAD)*(DTDZJ(J)+GAMAJ(J))/(DTDZJ(J)+GAMAD) 2567.5
1198     IF(RFAC.LT.0.) RFAC=0. 2567.51
1199     RFAC=RFAC+(TEM3*GRAV/(THSZM(J)*FCONT*BVF(J)+1.E-20))**2 2567.6
1200     c IF(RFAC.LT.0.01) RFAC=0.01 2568.
1201     ALAM=.573*SQRT(RFAC)/(1.-.427*RFAC**.302) 2568.5
1202     TEMPB=BETA*SHJ(J)/(FO(J)*GAMAP(j)+1.E-20) 2569.
1203     IF ((DTDYJ(J)+DTDYJ(J+1))*FO(J).GT.0.) TEMPB=0. 2569.1
1204    
1205     #if ( defined CPL_CHEM )
1206     !
1207     ! --- Get parameters:
1208     !
1209     beta1(j)=tempb
1210     beta2(j)=rfac
1211     !
1212     #endif
1213    
1214     DO 3170 I=1,IS 2569.5
1215     FLUXDT=0. 2570.
1216     FLUXDS=0. 2570.5
1217     DO 3175 L=1,LZM2 2571.
1218     ZHEI=.5*(PHI(1,J,L)+PHI(1,J,L+1)-2.*PHIS(1,J))/GRAV 2572.
1219     ZOVD=ZHEI/FUND 2572.5
1220     TEMDN1=.25*(DSHDY(J,L)+DSHDY(J,L+1)+DSHDY(J+1,L)+DSHDY(J+1,L+1)) 2573.
1221     IF(J.GT.JHALF.AND.TEMDN1.GT.-DQDYMN) TEMDN1=-DQDYMN 2573.5
1222     IF(J.LE.JHALF.AND.TEMDN1.LT.DQDYMN) TEMDN1=DQDYMN 2574.
1223     TERM1=ZOVD+.25*ZOVD*ZOVD*TEMPB*DSHDZ(J,L+1)/TEMDN1 2574.5
1224     DSHDZE=CP*(DTDZJ(J)+GAMAJ(J))*(1.-GAMAJ(J)/GAMAD)/HLAT 2575.
1225     IF(DSHDZE.LT.0.) DSHDZE=0. 2575.1
1226     TEMDN2=.25*(DSHDYS(J,L)+DSHDYS(J,L+1)+DSHDYS(J+1,L)+ 2575.5
1227     * DSHDYS(J+1,L+1)) 2576.
1228     IF(J.GT.JHALF.AND.TEMDN2.GT.-DQDYMN) TEMDN2=-DQDYMN 2576.5
1229     IF(J.LE.JHALF.AND.TEMDN2.LT.DQDYMN) TEMDN2=DQDYMN 2577.
1230     TERM2=ZOVD+.25*ZOVD*ZOVD*TEMPB*DSHDZE/RFAC/TEMDN2 2577.5
1231     TEM=(TEMPB/(1.+ALAM))*(TERM1+ALAM*TERM2/RFAC) 2578.
1232     WQ(J,L)=-.25*(VQ(J,L)+VQ(J,L+1)+VQ(J+1,L)+VQ(J+1,L+1))*TEM* 2578.5
1233     * (P(1,J)*SIGE(L+1)+PTOP)/SCLH(J,L)/P(1,J) 2579.
1234    
1235     #if ( defined CPL_CHEM )
1236     !
1237     ! --- Get parameters:
1238     !
1239     beta3(j,l+1)=zovd
1240     c 112696 added /p(1,j)
1241     beta4(j,l+1)=(p(1,j)*sige(l+1)+ptop)/sclh(j,l)/p(1,j)
1242     !
1243     #endif
1244    
1245     FLUXS=DTE*SP*DXYP(J)*WQ(J,L) 2579.5
1246     IF(FLUXS.GT..5* QT(I,J,L+1)*DSIG(L+1)) FLUXS=.5* QT(I,J,L+1)* 2580.
1247     * DSIG(L+1) 2580.5
1248     IF(FLUXS.LT.-.5* QT(I,J,L)*DSIG(L)) FLUXS=-.5* QT(I,J,L)*DSIG(L) 2581.
1249     QT(I,J,L)= QT(I,J,L)+(FLUXS-FLUXDS)/DSIG(L) 2581.5
1250     FLUXDS=FLUXS 2582.
1251     3175 CONTINUE 2582.5
1252     QT(I,J,LZM1)= QT(I,J,LZM1)-FLUXDS/DSIG(LZM1)
1253     3170 CONTINUE 2583.5
1254     endif ! MRCH
1255     endif
1256     c end wqtr
1257     C 2584.
1258     c print *,' TPRIM from eddypa'
1259     c print *,TPRIM2
1260     c do j=1,jm
1261     c do l=1,lm
1262     c print *,TPRIM2(J,L)*(EXPBYK(SIG(L)*P(1,J)+PTOP))**2.
1263     c enddo
1264     c enddo
1265    
1266     RETURN 2584.5
1267     END 2585.
1268     subroutine FLTR(X,JM,LM,JIN)
1269     dimension X(JM,LM),Y(JM,LM)
1270     do l=1,LM
1271     do j=1,JM
1272     Y(j,l)=X(j,l)
1273     enddo
1274     enddo
1275     do l=1,LM
1276     do j=JIN+1,JM-1
1277     X(j,l)=0.25*(X(j-1,l)+2.0*X(j,l)+X(j+1,l))
1278     enddo
1279     X(JIN,l)=(2.0*X(JIN,l)+X(JIN+1,l))/3.0
1280     X(JM,l)=(2.0*X(JM,l)+X(JM-1,l))/3.0
1281     enddo
1282     return
1283     end

  ViewVC Help
Powered by ViewVC 1.1.22