PRESUB.inc
References to this file elsewhere.
1 SUBROUTINE CRH2SH( WV, T, PS, A, B, IJMAX, LMAX, LARHM,
2 \ QMIN, ITERMX )
3
4 DIMENSION
5 1 WV(IJMAX,LMAX), T(IJMAX,LMAX), PS(IJMAX), A(LMAX), B(LMAX)
6
7 C INPUT WV : RELATIVE HUMIDITY( 0-1 ) ( L=<LARHM ) : CHANGED
8 C : SPECIFIC HUMIDITY( KG/KG ) ( L >LARHM ) : UNCHANGED
9 C T : VIRTUAL TEMPERATURE(K) : CHANGED
10 C PS : SURFACE PRESSURE(HPA)
11 C OUTPUT WV : SPECIFIC HUMIDITY( KG/KG ) ( FOR ALL LEVELS )
12 C T : REAL TEMPERATURE(K)
13
14 EE = 1.D0/(0.608D0*t_kelvin)
15 FF = 1.D0/(0.608D0*35.9D0)
16 GG = 6.11D0*0.622D0
17 CC = t_kelvin/35.9D0*7.5D0*LOG(10.D0)
18
19 DO 100 IJ = 1, IJMAX
20
21 DO 110 L = 1, LARHM
22
23 VT = T (IJ,L)
24 RH = WV(IJ,L)
25 IF( L.LE.LMAX-1 ) THEN
26 PHF1 = A(L )+B(L )*PS(IJ)
27 PHF2 = A(L+1)+B(L+1)*PS(IJ)
28 PL = EXP( ( PHF1*LOG(PHF1)-PHF2*LOG(PHF2) )
29 \ /( PHF1-PHF2 )-1.0 )
30 ELSE
31 PL = 0.5*( A(L)+B(L)*PS(IJ) )
32 END IF
33
34 AA = GG*RH/PL
35 BB = EE*(VT-t_kelvin)
36 DD = FF*(VT-35.9)
37 HH = CC*(BB-DD)
38
39 C 繰り返し計算開始 : QQ = INITIAL VALUE OF SPECIFIC HUMIDITY
40
41 QQ = AA*EXP(CC*BB/DD)
42
43 DO 120 ITER = 1, ITERMX
44 QV = AA*EXP(CC*(QQ-BB)/(QQ-DD))
45 FQ = QV-QQ/(1.0+0.608*QQ)
46 DFQ = HH*QV/((QQ-DD)**2)-1.0/(1.0+0.608*QQ)**2
47 QQ = QQ-FQ/DFQ
48 IF( QQ.LT.QMIN ) QQ = QMIN
49 120 CONTINUE
50
51 C 繰り返し計算終了→実温度の計算
52
53 TT = VT/(1.0+0.608*QQ)
54 T (IJ,L) = TT
55 WV(IJ,L) = QQ
56
57 110 CONTINUE
58
59 IF( LARHM.LT.LMAX ) THEN
60 DO 130 L = LARHM+1, LMAX
61 T(IJ,L) = T(IJ,L)/(1.0+0.608*WV(IJ,L))
62 130 CONTINUE
63 ENDIF
64
65 100 CONTINUE
66 C
67 RETURN
68 END SUBROUTINE CRH2SH
69 C**********************************************************************
70 SUBROUTINE CRH2SHA
71 I(IJMAX,KMAX,PS,A,B,GRAV,GASR,TLAPS,QCONS,QMIN,KST,ITERMX,
72 *-->
73 I IDX, LARHM,
74 *<--
75 O WV,T)
76 C
77 C*** CALCULATE SPECIFIC HUMIDITY AND REAL TEMPERATURE
78 C FROM RELATIVE HUMIDITY ,REAL TEMPERATURE AND VIRTUAL TEMPERATURE
79 C CORRECTION
80 C SPECIFIC HUMIDITY = CONSTANT IN STRATOSPHERE
81 C
82 C*** (ARRAYS) (INPUT) ******************* (OUTPUT) ******************
83 C WV RELATIVE HUMIDITY(NOT %) SPECIFIC HUMIDITY(KG/KG)
84 C T VIRTUAL TEMPERATURE(K) REAL TEMPERATURE(K)
85 C PS SURFACE PRESSURE (MB) (UNCHANGED)
86 C A,B A+B*PS DEFINES INTER-LAYER LEVEL
87 C***********************************************************************
88 DIMENSION WV(IJMAX,KMAX),T(IJMAX,KMAX),PS(IJMAX),A(KMAX),B(KMAX)
89 DATA Q90/2.7E-6/ ! SAME AS IN WV300M
90 COMMON/CTETEN/TABLE(25000)
91 COMMON/DTETEN/DTABLE(25000)
92 REAL*8 TABLE,DTABLE,X
93 LOGICAL FIRST/.TRUE./
94 C
95 IF(FIRST) THEN
96 GBYR=GRAV/GASR
97 ENDIF
98 C
99 DO 100 K=1,KMAX
100 DO 100 I= 1,IJMAX
101 VT=T (I,K)
102 IF( K.LT.LARHM ) THEN
103 RH=WV (I,K)
104 ENDIF
105 C
106 IF(K.LE.KMAX-1) THEN
107 PHF1=A(K )+B(K )*PS(I)
108 PHF2=A(K+1)+B(K+1)*PS(I)
109 PL =EXP( (PHF1*LOG(PHF1)-PHF2*LOG(PHF2))/(PHF1-PHF2)-1.0 )
110 ELSE
111 PL =0.5*(A(K)+B(K)*PS(I))
112 END IF
113 C
114 C
115 T0 = VT
116 YI = (T0-123.2D0)*100.0D0
117 IY = YI
118 IY = MAX( IY, 1 ) ! 96/01/31
119 IY = MIN( IY, 24999 ) ! 96/01/31
120 X = YI - IY
121 QSAT = ((1.0D0-X)*TABLE(IY)+X*TABLE(IY+1))/PL
122 CMM IF((IDX.EQ.1).AND.(PL.LE.90.)) THEN
123 CMM QQ=MIN(Q90,QSAT*0.9)
124 CMM ELSE
125 IF( K.LT.LARHM ) THEN
126 DQSAT = ((1.0D0-X)*DTABLE(IY)+X*DTABLE(IY+1))/PL
127 FT = VT-T0-0.608*T0*QSAT*RH
128 DFT = -1.0-0.608*QSAT*RH-0.608*T0*DQSAT*RH
129 T0 = T0 - FT/DFT
130 C
131 YI = (T0-123.2D0)*100.0D0
132 IY = YI
133 IY = MAX( IY, 1 ) ! 96/01/31
134 IY = MIN( IY, 24999 ) ! 96/01/31
135 X = YI - IY
136 QSAT = ((1.0D0-X)*TABLE(IY)+X*TABLE(IY+1))/PL
137 DQSAT = ((1.0D0-X)*DTABLE(IY)+X*DTABLE(IY+1))/PL
138 FT = VT-T0-0.608*T0*QSAT*RH
139 DFT = -1.0-0.608*QSAT*RH-0.608*T0*DQSAT*RH
140 T0 = T0 - FT/DFT
141 C
142 YI = (T0-123.2D0)*100.0D0
143 IY = YI
144 IY = MAX( IY, 1 ) ! 96/01/31
145 IY = MIN( IY, 24999 ) ! 96/01/31
146 X = YI - IY
147 QQ = ((1.0D0-X)*TABLE(IY)+X*TABLE(IY+1))/PL*RH
148 ELSE
149 QQ = WV(I,K)
150 ENDIF
151 C
152 IF((IDX.EQ.1).AND.(PL.LE.90.)) THEN
153 QQ=MIN(QQ,QSAT*0.9)
154 ENDIF
155 C
156 C +++ CALCULATE REAL TEMPERATURE +++
157 C
158 TT=VT/(1.+0.608*QQ)
159 C
160 T (I,K)=TT
161 WV(I,K)=QQ
162 C
163 100 CONTINUE
164 C
165 RETURN
166 END SUBROUTINE CRH2SHA
167 C ---------------------------------------------------------------------
168 SUBROUTINE SPLDIF( ZOUT, POUT, LMXOUT, ZIN, PIN, LMXIN )
169
170 DIMENSION ZOUT(LMXOUT), POUT(LMXOUT), ZIN(LMXIN), PIN(LMXIN)
171 DIMENSION SM(40), H(40), AL(40), AM(40), AP(40), C(40)
172
173 C INPUT / ZIN (L), PIN (L), LMXIN : INPUT.DATA, PRES(LOG), NUMBER
174 C OUTPUT/ ZOUT(L), POUT(L), LMXOUT : OUTPUT-VAL, PRES(LOG), NUMBER
175
176 LM1 = LMXIN-1
177 GR = -gravity/gas_constant
178
179 DO 110 L = 2, LMXIN
180 H(L) = PIN(L)-PIN(L-1)
181 110 CONTINUE
182
183 DO 120 L = 2, LM1
184 AL(L) = 0.5*H(L+1)/(H(L)+H(L+1))
185 AM(L) = 0.5-AL(L)
186 120 CONTINUE
187
188 C 終端条件( LAPSE RATE IS CONSTANT )
189 C SM(1) = SM(2) ; SM(LMXIN-1) = SM(LMXIN) : SECOND DERIVATIVE
190
191 AL(1) = -1.0
192 AM(LMXIN) = -1.0
193 AL(LMXIN) = 0.0
194
195 DO 130 L = 2, LMXIN
196 AP(L) = 1.0/(1.0-AL(L-1)*AM(L))
197 AL(L) = AL(L)*AP(L)
198 130 CONTINUE
199
200 C(1) = 0.0
201 C(LMXIN) = 0.0
202 DO 160 L = 2, LM1
203 C(L) = 3.0*((ZIN(L+1)-ZIN(L))/H(L+1) - (ZIN(L)-ZIN(L-1))/H(L))
204 \ /(H(L)+H(L+1))
205 160 CONTINUE
206
207 C FORWARD SUBSTITUTION
208
209 DO 200 L = 2, LMXIN
210 C(L) = (C(L)-C(L-1)*AM(L))*AP(L)
211 200 CONTINUE
212 SM(LMXIN) = C(LMXIN)
213
214 C BACKWARD SUBSTUTUTION
215
216 DO 220 K = 1, LM1
217 L = LMXIN-K
218 SM(L) = C(L)-AL(L)*SM(L+1)
219 220 CONTINUE
220
221 C INTERPOLATION
222
223 LB = 2
224
225 DO 500 LOUT = 1, LMXOUT
226
227 X = POUT(LOUT)
228 DO 300 L = LB, LMXIN
229 IF( X.GE.PIN(L) ) GO TO 310
230 300 CONTINUE
231 L = LMXIN
232 310 LB = L
233
234 C ３次スプライン関数の微分
235
236 ZOUT(LOUT) = SM(L-1)*( -(PIN(L)-X)**2 /(2.0*H(L))+H(L)/6.0 )
237 \ + SM(L) *( (X-PIN(L-1))**2 /(2.0*H(L))-H(L)/6.0 )
238 \ + ( ZIN(L)-ZIN(L-1) )/H(L)
239 ZOUT(LOUT) = ZOUT(LOUT)*GR
240
241 500 CONTINUE
242
243 RETURN
244 END SUBROUTINE SPLDIF
245 C ---------------------------------------------------------------------
246 SUBROUTINE SPLDIF3( ZOUT, POUT, LMXOUT, ZIN, PIN, LMXIN, IJMAX,
247 W SM, H, AL, AM, AP, C )
248 DIMENSION ZOUT(IJMAX,LMXOUT), POUT(IJMAX,LMXOUT),
249 1 ZIN (IJMAX,LMXIN), PIN (IJMAX,LMXIN)
250 CMM DIMENSION ZOUT(LMXOUT), POUT(LMXOUT), ZIN(LMXIN), PIN(LMXIN)
251 DIMENSION SM(IJMAX,LMXIN), H(IJMAX,LMXIN), AL(IJMAX,LMXIN),
252 1 AM(IJMAX,LMXIN), AP(IJMAX,LMXIN), C(IJMAX,LMXIN)
253 CMM DIMENSION SM(40), H(40), AL(40), AM(40), AP(40), C(40)
254
255 C INPUT / ZIN (L), PIN (L), LMXIN : INPUT.DATA, PRES(LOG), NUMBER
256 C OUTPUT/ ZOUT(L), POUT(L), LMXOUT : OUTPUT-VAL, PRES(LOG), NUMBER
257
258 LM1 = LMXIN-1
259 GR = -gravity/gas_constant
260
261 DO 110 L = 2, LMXIN
262 DO 110 I = 1, IJMAX
263 H(I,L) = PIN(I,L)-PIN(I,L-1)
264 110 CONTINUE
265
266 DO 120 L = 2, LM1
267 DO 120 I = 1, IJMAX
268 AL(I,L) = 0.5*H(I,L+1)/(H(I,L)+H(I,L+1))
269 AM(I,L) = 0.5-AL(I,L)
270 120 CONTINUE
271
272 C 終端条件( LAPSE RATE IS CONSTANT )
273 C SM(1) = SM(2) ; SM(LMXIN-1) = SM(LMXIN) : SECOND DERIVATIVE
274 DO 125 I = 1, IJMAX
275 AL(I,1) = -1.0
276 AM(I,LMXIN) = -1.0
277 AL(I,LMXIN) = 0.0
278 125 CONTINUE
279
280 DO 130 L = 2, LMXIN
281 DO 130 I = 1, IJMAX
282 C
283 AP(I,L) = 1.0/(1.0-AL(I,L-1)*AM(I,L))
284 AL(I,L) = AL(I,L)*AP(I,L)
285 130 CONTINUE
286
287 DO 155 I = 1, IJMAX
288 C(I,1) = 0.0
289 C(I,LMXIN) = 0.0
290 155 CONTINUE
291 C
292 DO 160 L = 2, LM1
293 DO 160 I = 1, IJMAX
294 C(I,L) = 3.0*((ZIN(I,L+1)-ZIN(I,L))/H(I,L+1)
295 1 - (ZIN(I,L)-ZIN(I,L-1))/H(I,L))
296 2 /(H(I,L)+H(I,L+1))
297 160 CONTINUE
298
299 C FORWARD SUBSTITUTION
300
301 DO 200 L = 2, LMXIN
302 DO 200 I = 1, IJMAX
303 C(I,L) = (C(I,L)-C(I,L-1)*AM(I,L))*AP(I,L)
304 200 CONTINUE
305 DO 205 I = 1, IJMAX
306 SM(I,LMXIN) = C(I,LMXIN)
307 205 CONTINUE
308
309 C BACKWARD SUBSTUTUTION
310
311 DO 220 K = 1, LM1
312 DO 220 I = 1, IJMAX
313 L = LMXIN-K
314 SM(I,L) = C(I,L)-AL(I,L)*SM(I,L+1)
315 220 CONTINUE
316
317 C INTERPOLATION
318
319 C LB = 2
320
321 DO 500 LOUT = 1, LMXOUT
322 DO 500 I = 1, IJMAX
323
324 X = POUT(I,LOUT)
325 L = LOUT ! FOR ONLY PIN.EQ.POUT
326 IF( L.LT.2 ) L=2
327 CM DO 300 L = LB, LMXIN
328 CM IF( X.GE.PIN(L) ) GO TO 310
329 CM300 CONTINUE
330 CM L = LMXIN
331 CM310 LB = L
332
333 ZOUT(I,LOUT) = SM(I,L-1)*( -(PIN(I,L)-X)**2
334 1 /(2.0*H(I,L))+H(I,L)/6.0 )
335 2 + SM(I,L) *( (X-PIN(I,L-1))**2
336 3 /(2.0*H(I,L))-H(I,L)/6.0 )
337 \ + ( ZIN(I,L)-ZIN(I,L-1) )/H(I,L)
338 ZOUT(I,LOUT) = ZOUT(I,LOUT)*GR
339
340 500 CONTINUE
341
342 RETURN
343 END SUBROUTINE SPLDIF3
344 C======================================================================
345 SUBROUTINE SETWHT
346 1(IMAX,JMAX,DP,IP,DCOSCL,DGW,COCLT)
347 DIMENSION DP(4,IMAX,JMAX)
348 INTEGER*2 IP(2,IMAX,JMAX)
349 REAL*8 DCOSCL(JMAX),DGW(JMAX),COCLT(JMAX)
350 C
351 Crizvi Already defined in module_wave2grid_kma
352 C JMAXHF=JMAX/2
353 DELX=360.0/FLOAT(IMAX)
354 CALL GAUSS(DCOSCL,DGW,JMAX)
355 C
356 R2D=180.0/pi
357 DO 100 J=1,JMAXHF
358 COCLT( J)= R2D*ACOS(DCOSCL(J))
359 COCLT(JMAX+1-J)=180.0-R2D*ACOS(DCOSCL(J))
360 100 CONTINUE
361 C WRITE(6,*) 'CHECK OF COLATTITUDE',COCLT
362 C
363 DO 120 J=1,JMAX
364 PX=COCLT(J)+1.5
365 LAT=PX
366 DLAT=PX-FLOAT(LAT)
367 DO 120 I=1,IMAX
368 QX=DELX*FLOAT(I-1)+1.5
369 LON=QX
370 DLON=QX-FLOAT(LON)
371 IP(1,I,J)=LON
372 IP(2,I,J)=LAT
373 DP(1,I,J)=(1.0-DLAT)*(1.0-DLON)
374 DP(2,I,J)=(1.0-DLAT)*DLON
375 DP(3,I,J)=DLAT*DLON
376 DP(4,I,J)=DLAT*(1.0-DLON)
377 120 CONTINUE
378 C DO 1000 J=1,JMAX
379 C WRITE(6,1010) J,IP(1,1,J),IP(2,1,J),(DP(I,1,J),I=1,4)
380 C1000 CONTINUE
381 C DO 1100 J=1,JMAX
382 C WRITE(6,1010) J,IP(1,IMAX,J),IP(2,IMAX,J),(DP(I,IMAX,J),I=1,4)
383 C1100 CONTINUE
384 C DO 1200 J=1,IMAX
385 C WRITE(6,1020) J,IP(1,J,1),IP(2,J,1),(DP(I,J,1),I=1,4)
386 C1200 CONTINUE
387 C DO 1300 J=1,IMAX
388 C WRITE(6,1020) J,IP(1,J,JMAX),IP(2,J,JMAX),(DP(I,J,JMAX),I=1,4)
389 C1300 CONTINUE
390 C1010 FORMAT(1H ,'J=',I3,5X,2I5,5X,4F10.3)
391 C1020 FORMAT(1H ,'I=',I3,5X,2I5,5X,4F10.3)
392 C
393 RETURN
394 END SUBROUTINE SETWHT
395 SUBROUTINE INTERP
396 1(WORK,DATA,IMAX,JMAX,DP,IP)
397 C
398 DIMENSION WORK(362,182),DATA(IMAX,JMAX),DP(4,IMAX,JMAX)
399 INTEGER*2 IP(2,IMAX,JMAX)
400 C
401 DO 100 J=1,JMAX
402 DO 100 I=1,IMAX
403 II=IP(1,I,J)
404 JJ=IP(2,I,J)
405 DATA(I,J)=DP(1,I,J)*WORK(II,JJ )+DP(2,I,J)*WORK(II+1,JJ )
406 * +DP(4,I,J)*WORK(II,JJ+1)+DP(3,I,J)*WORK(II+1,JJ+1)
407 100 CONTINUE
408 C
409 RETURN
410 END SUBROUTINE INTERP