/[MITgcm]/MITgcm/pkg/fizhi/fizhi_utils.F
ViewVC logotype

Contents of /MITgcm/pkg/fizhi/fizhi_utils.F

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


Revision 1.5 - (show annotations) (download)
Wed Jul 14 17:31:58 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
Changes since 1.4: +1 -0 lines
Debugging

1 function minval (q,im)
2 implicit none
3 #include "CPP_EEOPTIONS.h"
4 integer im, i
5 _RL q(im), minval
6 minval = 1.e15
7 do i=1,im
8 if( q(i).lt.minval ) minval = q(i)
9 enddo
10 return
11 end
12 FUNCTION ERRF (ARG)
13 C***********************************************************************
14 C FUNCTION ERRF
15 C PURPOSE
16 C COMPUTES ERROR FUNCTION OF ARGUMENT
17 C USAGE
18 C CALLED BY TRBFLX
19 C DESCRIPTION OF PARAMETERS
20 C ARG - INPUTED ARGUMENT
21 C REMARKS:
22 C USED TO COMPUTE FRACTIONAL CLOUD COVER AND LIQUID WATER CONTENT
23 C FROM TURBULENCE STATISTICS
24 C **********************************************************************
25
26 PARAMETER ( AA1 = 0.254829592 )
27 PARAMETER ( AA2 = -0.284496736 )
28 PARAMETER ( AA3 = 1.421413741 )
29 PARAMETER ( AA4 = -1.453152027 )
30 PARAMETER ( AA5 = 1.061405429 )
31 PARAMETER ( PP = 0.3275911 )
32 PARAMETER ( X2 = AA2 / AA1 )
33 PARAMETER ( X3 = AA3 / AA2 )
34 PARAMETER ( X4 = AA5 / AA3 )
35 PARAMETER ( X5 = AA5 / AA4 )
36
37 ERRF = 1.
38 AARG=ABS(ARG)
39
40 IF ( AARG .LT. 4.0 ) THEN
41 TT = 1./(1.+PP*AARG)
42 ERRF = 1. -
43 1 (AA1*TT*(1.+X2*TT*(1.+X3*TT*(1.+X4*TT*(1.+X5*TT)))))
44 2 * EXP(-AARG*AARG)
45 ENDIF
46
47 IF ( ARG .LT. 0.0 ) ERRF = -ERRF
48
49 RETURN
50 END
51
52 SUBROUTINE STRIP(A,B,IA,IB,L,K)
53 #include "CPP_EEOPTIONS.h"
54 _RL A(IA,L), B(IB,L)
55 INTEGER OFFSET
56
57 OFFSET = IB*(K-1)
58 LEN = MIN(IB,IA-OFFSET)
59 OFFSET = OFFSET+1
60
61 IF(LEN.EQ.IB) THEN
62 DO 100 J=1,L
63 DO 100 I=1,LEN
64 B(I,J) = A(I+OFFSET-1,J)
65 100 CONTINUE
66 ELSE
67 DO 200 J=1,L
68 DO 300 I=1,LEN
69 B(I,J) = A(I+OFFSET-1,J)
70 300 CONTINUE
71 DO 400 I=1,IB-LEN
72 B(LEN+I,J) = A(LEN+OFFSET-1,J)
73 400 CONTINUE
74 200 CONTINUE
75 ENDIF
76
77 RETURN
78 END
79 SUBROUTINE PASTE(B,A,IB,IA,L,K)
80 #include "CPP_EEOPTIONS.h"
81 _RL A(IA,L), B(IB,L)
82 INTEGER OFFSET
83
84 OFFSET = IB*(K-1)
85 LEN = MIN(IB,IA-OFFSET)
86 OFFSET = OFFSET+1
87
88 DO 100 J=1,L
89 DO 100 I=1,LEN
90 A(I+OFFSET-1,J) = B(I,J)
91 100 CONTINUE
92
93 RETURN
94 END
95 SUBROUTINE PSTBMP(B,A,IB,IA,L,K)
96 #include "CPP_EEOPTIONS.h"
97 _RL A(IA,L), B(IB,L)
98 INTEGER OFFSET
99
100 OFFSET = IB*(K-1)
101 LEN = MIN(IB,IA-OFFSET)
102 OFFSET = OFFSET+1
103
104 DO 100 J=1,L
105 DO 100 I=1,LEN
106 A(I+OFFSET-1,J) = A(I+OFFSET-1,J) + B(I,J)
107 100 CONTINUE
108 C
109 RETURN
110 END
111 SUBROUTINE STRINT(A,B,IA,IB,L,K)
112 INTEGER A(IA,L), B(IB,L)
113 INTEGER OFFSET
114
115 OFFSET = IB*(K-1)
116 LEN = MIN(IB,IA-OFFSET)
117 OFFSET = OFFSET+1
118
119 IF(LEN.EQ.IB) THEN
120 DO 100 J=1,L
121 DO 100 I=1,LEN
122 B(I,J) = A(I+OFFSET-1,J)
123 100 CONTINUE
124 ELSE
125 DO 200 J=1,L
126 DO 300 I=1,LEN
127 B(I,J) = A(I+OFFSET-1,J)
128 300 CONTINUE
129 DO 400 I=1,IB-LEN
130 B(LEN+I,J) = A(LEN+OFFSET-1,J)
131 400 CONTINUE
132 200 CONTINUE
133 ENDIF
134
135 RETURN
136 END
137 SUBROUTINE QSAT (TT,P,Q,DQDT,LDQDT)
138 C***********************************************************************
139 C
140 C PURPOSE:
141 C ========
142 C Compute Saturation Specific Humidity
143 C
144 C INPUT:
145 C ======
146 C TT ......... Temperature (Kelvin)
147 C P .......... Pressure (mb)
148 C LDQDT ...... Logical Flag to compute QSAT Derivative
149 C
150 C OUTPUT:
151 C =======
152 C Q .......... Saturation Specific Humidity
153 C DQDT ....... Saturation Specific Humidity Derivative wrt Temperature
154 C
155 C
156 C***********************************************************************
157 C* GODDARD LABORATORY FOR ATMOSPHERES *
158 C***********************************************************************
159
160 IMPLICIT NONE
161 #include "CPP_EEOPTIONS.h"
162 _RL TT, P, Q, DQDT
163 LOGICAL LDQDT
164 _RL AIRMW, H2OMW
165
166 PARAMETER ( AIRMW = 28.97 )
167 PARAMETER ( H2OMW = 18.01 )
168
169 _RL ESFAC, ERFAC
170 PARAMETER ( ESFAC = H2OMW/AIRMW )
171 PARAMETER ( ERFAC = (1.0-ESFAC)/ESFAC )
172
173 _RL aw0, aw1, aw2, aw3, aw4, aw5, aw6
174 _RL bw0, bw1, bw2, bw3, bw4, bw5, bw6
175 _RL ai0, ai1, ai2, ai3, ai4, ai5, ai6
176 _RL bi0, bi1, bi2, bi3, bi4, bi5, bi6
177
178 _RL d0, d1, d2, d3, d4, d5, d6
179 _RL e0, e1, e2, e3, e4, e5, e6
180 _RL f0, f1, f2, f3, f4, f5, f6
181 _RL g0, g1, g2, g3, g4, g5, g6
182
183 c ********************************************************
184 c *** Polynomial Coefficients WRT Water (Lowe, 1977) ****
185 c *** (Valid +50 C to -50 C) ****
186 c ********************************************************
187
188 parameter ( aw0 = 6.107799961e+00 * esfac )
189 parameter ( aw1 = 4.436518521e-01 * esfac )
190 parameter ( aw2 = 1.428945805e-02 * esfac )
191 parameter ( aw3 = 2.650648471e-04 * esfac )
192 parameter ( aw4 = 3.031240396e-06 * esfac )
193 parameter ( aw5 = 2.034080948e-08 * esfac )
194 parameter ( aw6 = 6.136820929e-11 * esfac )
195
196 parameter ( bw0 = +4.438099984e-01 * esfac )
197 parameter ( bw1 = +2.857002636e-02 * esfac )
198 parameter ( bw2 = +7.938054040e-04 * esfac )
199 parameter ( bw3 = +1.215215065e-05 * esfac )
200 parameter ( bw4 = +1.036561403e-07 * esfac )
201 parameter ( bw5 = +3.532421810e-10 * esfac )
202 parameter ( bw6 = -7.090244804e-13 * esfac )
203
204
205 c ********************************************************
206 c *** Polynomial Coefficients WRT Ice (Lowe, 1977) ****
207 c *** (Valid +0 C to -50 C) ****
208 c ********************************************************
209
210 parameter ( ai0 = +6.109177956e+00 * esfac )
211 parameter ( ai1 = +5.034698970e-01 * esfac )
212 parameter ( ai2 = +1.886013408e-02 * esfac )
213 parameter ( ai3 = +4.176223716e-04 * esfac )
214 parameter ( ai4 = +5.824720280e-06 * esfac )
215 parameter ( ai5 = +4.838803174e-08 * esfac )
216 parameter ( ai6 = +1.838826904e-10 * esfac )
217
218 parameter ( bi0 = +5.030305237e-01 * esfac )
219 parameter ( bi1 = +3.773255020e-02 * esfac )
220 parameter ( bi2 = +1.267995369e-03 * esfac )
221 parameter ( bi3 = +2.477563108e-05 * esfac )
222 parameter ( bi4 = +3.005693132e-07 * esfac )
223 parameter ( bi5 = +2.158542548e-09 * esfac )
224 parameter ( bi6 = +7.131097725e-12 * esfac )
225
226
227 c ********************************************************
228 c *** Polynomial Coefficients WRT Ice ****
229 c *** Starr and Cox (1985) (Valid -40 C to -70 C) ****
230 c ********************************************************
231
232
233 parameter ( d0 = 0.535098336e+01 * esfac )
234 parameter ( d1 = 0.401390832e+00 * esfac )
235 parameter ( d2 = 0.129690326e-01 * esfac )
236 parameter ( d3 = 0.230325039e-03 * esfac )
237 parameter ( d4 = 0.236279781e-05 * esfac )
238 parameter ( d5 = 0.132243858e-07 * esfac )
239 parameter ( d6 = 0.314296723e-10 * esfac )
240
241 parameter ( e0 = 0.469290530e+00 * esfac )
242 parameter ( e1 = 0.333092511e-01 * esfac )
243 parameter ( e2 = 0.102164528e-02 * esfac )
244 parameter ( e3 = 0.172979242e-04 * esfac )
245 parameter ( e4 = 0.170017544e-06 * esfac )
246 parameter ( e5 = 0.916466531e-09 * esfac )
247 parameter ( e6 = 0.210844486e-11 * esfac )
248
249
250 c ********************************************************
251 c *** Polynomial Coefficients WRT Ice ****
252 c *** Starr and Cox (1985) (Valid -65 C to -95 C) ****
253 c ********************************************************
254
255 parameter ( f0 = 0.298152339e+01 * esfac )
256 parameter ( f1 = 0.191372282e+00 * esfac )
257 parameter ( f2 = 0.517609116e-02 * esfac )
258 parameter ( f3 = 0.754129933e-04 * esfac )
259 parameter ( f4 = 0.623439266e-06 * esfac )
260 parameter ( f5 = 0.276961083e-08 * esfac )
261 parameter ( f6 = 0.516000335e-11 * esfac )
262
263 parameter ( g0 = 0.312654072e+00 * esfac )
264 parameter ( g1 = 0.195789002e-01 * esfac )
265 parameter ( g2 = 0.517837908e-03 * esfac )
266 parameter ( g3 = 0.739410547e-05 * esfac )
267 parameter ( g4 = 0.600331350e-07 * esfac )
268 parameter ( g5 = 0.262430726e-09 * esfac )
269 parameter ( g6 = 0.481960676e-12 * esfac )
270
271 _RL TMAX, TICE
272 PARAMETER ( TMAX=323.15, TICE=273.16)
273
274 _RL T, D, W, QX, DQX
275 T = MIN(TT,TMAX) - TICE
276 DQX = 0.
277 QX = 0.
278
279 c Fitting for temperatures above 0 degrees centigrade
280 c ---------------------------------------------------
281 if(t.gt.0.) then
282 qx = aw0+T*(aw1+T*(aw2+T*(aw3+T*(aw4+T*(aw5+T*aw6)))))
283 if (ldqdt) then
284 dqx = bw0+T*(bw1+T*(bw2+T*(bw3+T*(bw4+T*(bw5+T*bw6)))))
285 endif
286 endif
287
288 c Fitting for temperatures between 0 and -40
289 c ------------------------------------------
290 if( t.le.0. .and. t.gt.-40.0 ) then
291 w = (40.0 + t)/40.0
292 qx = w *(aw0+T*(aw1+T*(aw2+T*(aw3+T*(aw4+T*(aw5+T*aw6))))))
293 . + (1.-w)*(ai0+T*(ai1+T*(ai2+T*(ai3+T*(ai4+T*(ai5+T*ai6))))))
294 if (ldqdt) then
295 dqx = w *(bw0+T*(bw1+T*(bw2+T*(bw3+T*(bw4+T*(bw5+T*bw6))))))
296 . + (1.-w)*(bi0+T*(bi1+T*(bi2+T*(bi3+T*(bi4+T*(bi5+T*bi6))))))
297 endif
298 endif
299
300 c Fitting for temperatures between -40 and -70
301 c --------------------------------------------
302 if( t.le.-40.0 .and. t.ge.-70.0 ) then
303 qx = d0+T*(d1+T*(d2+T*(d3+T*(d4+T*(d5+T*d6)))))
304 if (ldqdt) then
305 dqx = e0+T*(e1+T*(e2+T*(e3+T*(e4+T*(e5+T*e6)))))
306 endif
307 endif
308
309 c Fitting for temperatures less than -70
310 c --------------------------------------
311 if(t.lt.-70.0) then
312 qx = f0+t*(f1+t*(f2+t*(f3+t*(f4+t*(f5+t*f6)))))
313 if (ldqdt) then
314 dqx = g0+t*(g1+t*(g2+t*(g3+t*(g4+t*(g5+t*g6)))))
315 endif
316 endif
317
318 c Compute Saturation Specific Humidity
319 c ------------------------------------
320 D = (P-ERFAC*QX)
321 IF(D.LT.0.) THEN
322 Q = 1.0
323 IF (LDQDT) DQDT = 0.
324 ELSE
325 D = 1.0 / D
326 Q = MIN(QX * D,1.0)
327 IF (LDQDT) DQDT = (1.0 + ERFAC*Q) * D * DQX
328 ENDIF
329 RETURN
330 END
331 subroutine vqsat (tt,p,q,dqdt,ldqdt,n)
332 implicit none
333 #include "CPP_EEOPTIONS.h"
334 integer i,n
335 logical ldqdt
336 _RL tt(n), p(n), q(n), dqdt(n)
337 #if CRAY
338 #if f77
339 cfpp$ expand (qsat)
340 #endif
341 #endif
342 do i=1,n
343 call qsat ( tt(i),p(i),q(i),dqdt(i),ldqdt )
344 enddo
345 return
346 end
347
348 subroutine stripit(a,b,irun,ia,ib,l,k)
349 implicit none
350 #include "CPP_EEOPTIONS.h"
351 integer ia,ib,irun,l,k
352 _RL a(ia,l), b(ib,l)
353 integer i,j,len,offset
354
355 offset = ib*(k-1)
356 len = min(ib,irun-offset)
357 offset = offset+1
358
359 if(len.eq.ib) then
360 do 100 j=1,l
361 do 100 i=1,len
362 b(i,j) = a(i+offset-1,j)
363 100 continue
364 else
365 do 200 j=1,l
366 do 300 i=1,len
367 b(i,j) = a(i+offset-1,j)
368 300 continue
369 do 400 i=1,ib-len
370 b(len+i,j) = a(len+offset-1,j)
371 400 continue
372 200 continue
373 endif
374 return
375 end
376
377 subroutine stripitint(a,b,irun,ia,ib,l,k)
378 implicit none
379 #include "CPP_EEOPTIONS.h"
380 integer ia,ib,irun,l,k,a(ia,l),b(ib,l)
381 integer i,j,len,offset
382
383 offset = ib*(k-1)
384 len = min(ib,irun-offset)
385 offset = offset+1
386
387 if(len.eq.ib) then
388 do 100 j=1,l
389 do 100 i=1,len
390 b(i,j) = a(i+offset-1,j)
391 100 continue
392 else
393 do 200 j=1,l
394 do 300 i=1,len
395 b(i,j) = a(i+offset-1,j)
396 300 continue
397 do 400 i=1,ib-len
398 b(len+i,j) = a(len+offset-1,j)
399 400 continue
400 200 continue
401 endif
402 return
403 end
404
405 subroutine pastit(b,a,ib,ia,irun,L,k)
406 implicit none
407 #include "CPP_EEOPTIONS.h"
408 _RL a(ia,l), b(ib,l)
409 integer ib,ia,L,k,irun,len,offset
410 integer i,j
411
412 offset = ib*(k-1)
413 len = min(ib,irun-offset)
414 offset = offset+1
415
416 do 100 j=1,L
417 do 100 i=1,len
418 a(i+offset-1,j) = b(i,j)
419 100 continue
420 return
421 end
422
423 subroutine pstbitint(b,a,ib,ia,irun,l,k)
424 implicit none
425 #include "CPP_EEOPTIONS.h"
426 _RL a(ia,l)
427 integer b(ib,l)
428 integer ib,ia,L,k,irun,len,offset
429 integer i,j
430
431 offset = ib*(k-1)
432 len = min(ib,irun-offset)
433 offset = offset+1
434
435 do 100 j=1,L
436 do 100 i=1,len
437 a(i+offset-1,j) = a(i+offset-1,j) + float(b(i,j))
438 100 continue
439 return
440 end
441
442
443 subroutine pstbmpit(b,a,ib,ia,irun,l,k)
444 implicit none
445 #include "CPP_EEOPTIONS.h"
446 _RL a(ia,l), b(ib,l)
447 integer ib,ia,L,k,irun,len,offset
448 integer i,j
449
450 offset = ib*(k-1)
451 len = min(ib,irun-offset)
452 offset = offset+1
453
454 do 100 j=1,L
455 do 100 i=1,len
456 a(i+offset-1,j) = a(i+offset-1,j) + b(i,j)
457 100 continue
458 return
459 end
460
461 subroutine strip2tile(a,index,b,irun,ia,ib,levs,npeice)
462 c-----------------------------------------------------------------------
463 c subroutine strip2tile - extract one processors worth of grid points
464 c from a grid space array to a stripped tile
465 c space array
466 c
467 c input:
468 c a - array to be stripped FROM [ia,levs]
469 c index - array of horizontal indeces of grid points to convert to
470 c tile space
471 c irun - number of points in array a that need to be stripped
472 c ia - inner of dimension of source array
473 c ib - inner dimension of target array AND the number of points
474 c in the target array to be filled
475 c levs - number of vertical levels AND outer dimension of arrays
476 c npeice - the current strip number to be filled
477 c output:
478 c b - array to be filled, ie, one processors field [ib,levs]
479 c-----------------------------------------------------------------------
480 implicit none
481 #include "CPP_EEOPTIONS.h"
482 integer ia,ib,irun,levs,npeice
483 _RL a(ia,levs), b(ib,levs)
484 integer index(irun)
485 integer i,k,len,offset
486
487 offset = ib*(npeice-1)
488 len = min(ib,irun-offset)
489 offset = offset+1
490
491 if(len.eq.ib) then
492 do 100 k=1,levs
493 do 100 i=1,len
494 b(i,k) = a(index(i+offset-1),k)
495 100 continue
496 else
497 do 200 k=1,levs
498 do 300 i=1,len
499 b(i,k) = a(index(i+offset-1),k)
500 300 continue
501 do 400 i=1,ib-len
502 b(len+i,k) = a(index(len+offset-1),k)
503 400 continue
504 200 continue
505 endif
506 return
507 end
508
509 subroutine paste2grd_old(b,index,chfr,ib,numpts,a,ia,levs,npeice)
510 c-----------------------------------------------------------------------
511 c subroutine paste2grd - paste one processors worth of grid points
512 c from a stripped tile array to a grid
513 c space array
514 c
515 c input:
516 c b - array to be pasted back into grid space [ib,levs]
517 c index - array of horizontal indeces of grid points to convert to
518 c tile space[numpts]
519 c chfr - fractional area covered by the tile [ib]
520 c ib - inner dimension of source array AND number of points in
521 c array a that need to be pasted
522 c numpts - total number of points which were stripped
523 c ia - inner of dimension of target array
524 c levs - number of vertical levels AND outer dimension of arrays
525 c npeice - the current strip number to be filled
526 c output:
527 c a - grid space array to be filled [ia,levs]
528 c
529 c IMPORTANT NOTE:
530 c
531 c This routine will result in roundoff differences if called from
532 c within a parallel region.
533 c-----------------------------------------------------------------------
534
535 implicit none
536 #include "CPP_EEOPTIONS.h"
537 integer ia,ib,levs,numpts,npeice
538 integer index(numpts)
539 _RL a(ia,levs), b(ib,levs), chfr(ib)
540 real b(ib,levs)
541
542 integer i,L,offset,len
543
544 offset = ib*(npeice-1)
545 len = min(ib,numpts-offset)
546 offset = offset+1
547
548 do L = 1,levs
549 do i=1,len
550 a(index(i+offset-1),L) = a(index(i+offset-1),L) + b(i,L)*chfr(i)
551 enddo
552 enddo
553 return
554 end
555 subroutine paste2grd (b,index,chfr,ib,numpts,a,ia,levs,npeice,check)
556 c-----------------------------------------------------------------------
557 c subroutine paste2grd - paste one processors worth of grid points
558 c from a stripped tile array to a grid
559 c space array
560 c
561 c input:
562 c b - array to be pasted back into grid space [ib,levs]
563 c index - array of horizontal indeces of grid points to convert to
564 c tile space[numpts]
565 c chfr - fractional area covered by the tile [ib]
566 c ib - inner dimension of source array AND number of points in
567 c array a that need to be pasted
568 c numpts - total number of points which were stripped
569 c ia - inner of dimension of target array
570 c levs - number of vertical levels AND outer dimension of arrays
571 c npeice - the current strip number to be filled
572 c check - logical to check for undefined values
573 c output:
574 c a - grid space array to be filled [ia,levs]
575 c
576 c IMPORTANT NOTE:
577 c
578 c This routine will result in roundoff differences if called from
579 c within a parallel region.
580 c-----------------------------------------------------------------------
581
582 implicit none
583 #include "CPP_EEOPTIONS.h"
584 integer ia,ib,levs,numpts,npeice
585 integer index(numpts)
586 _RL a(ia,levs), b(ib,levs), chfr(ib)
587 logical check
588
589 integer i,L,offset,len
590 _RL undef,getcon
591
592 offset = ib*(npeice-1)
593 len = min(ib,numpts-offset)
594 offset = offset+1
595
596 if( check ) then
597 undef = getcon('UNDEF')
598 do L= 1,levs
599 do i= 1,len
600 if( a(index(i+offset-1),L).eq.undef .or. b(i,L).eq.undef ) then
601 a(index(i+offset-1),L) = undef
602 else
603 a(index(i+offset-1),L) = a(index(i+offset-1),L) + b(i,L)*chfr(i)
604 endif
605 enddo
606 enddo
607 else
608 do L= 1,levs
609 do i= 1,len
610 a(index(i+offset-1),L) = a(index(i+offset-1),L) + b(i,L)*chfr(i)
611 enddo
612 enddo
613 endif
614
615 return
616 end
617 SUBROUTINE GRD2MSC(A,IM,JM,IGRD,B,MXCHPS,NCHP)
618
619 implicit none
620 #include "CPP_EEOPTIONS.h"
621 integer im,jm,mxchps,nchp
622 integer igrd(mxchps)
623 _RL A(IM,JM), B(MXCHPS)
624
625 integer i
626
627 IF(NCHP.GE.0) THEN
628 DO I=1,NCHP
629 B(I) = A(IGRD(I),1)
630 ENDDO
631 ELSE
632 PRINT *, 'ERROR IN GRD2MSC'
633 ENDIF
634
635 RETURN
636 END
637
638 SUBROUTINE MSC2GRD(IGRD,CHFR,B,MXCHPS,NCHP,FRACG,A,IM,JM)
639
640 implicit none
641 #include "CPP_EEOPTIONS.h"
642 _RL zero,one
643 parameter ( one = 1.)
644 parameter (zero = 0.)
645 integer im,jm,mxchps,nchp
646 integer igrd(mxchps)
647 _RL A(IM,JM), B(MXCHPS), CHFR(MXCHPS), FRACG(IM,JM)
648
649 _RL VT1(IM,JM)
650 integer i
651
652 IF(NCHP.GE.0) THEN
653 DO I=1,IM*JM
654 VT1(I,1) = ZERO
655 ENDDO
656
657 DO I=1,NCHP
658 VT1(IGRD(I),1) = VT1(IGRD(I),1) + B(I)*CHFR(I)
659 ENDDO
660
661 DO I=1,IM*JM
662 A(I,1) = A(I,1)*(ONE-FRACG(I,1)) + VT1(I,1)
663 ENDDO
664 ELSE
665 PRINT *, 'ERROR IN MSC2GRD'
666 ENDIF
667
668 RETURN
669 END
670
671 subroutine chpprm(nymd,nhms,mxchps,nchp,chlt,ityp,alai,
672 1 agrn,zoch,z2ch,cdrc,cdsc,sqsc,ufac,rsl1,rsl2,rdcs)
673
674 implicit none
675 #include "CPP_EEOPTIONS.h"
676
677 integer nymd,nhms,nchp,mxchps,ityp(mxchps)
678 _RL chlt(mxchps)
679 _RL alai(mxchps),agrn(mxchps)
680 _RL zoch(mxchps), z2ch(mxchps), cdrc(mxchps), cdsc(mxchps)
681 _RL sqsc(mxchps), ufac(mxchps), rsl1(mxchps), rsl2(mxchps)
682 _RL rdcs(mxchps)
683
684 C*********************************************************************
685 C********************* SUBROUTINE CHPPRM ****************************
686 C********************** 14 JUNE 1991 ******************************
687 C*********************************************************************
688
689 integer ntyps
690 parameter (ntyps=10)
691
692 _RL pblzet
693 parameter (pblzet = 50.)
694 integer k1,k2,nymd1,nhms1,nymd2,nhms2,i
695 _RL getcon,vkrm,rootl,vroot,dum1,dum2,alphaf
696 _RL facm,facp
697 _RL scat,d
698
699 _RL
700 & vgdd(12, ntyps), vgz0(12, ntyps),
701 & vgrd(12, ntyps), vgrt(12, ntyps),
702
703 & vgrf11(ntyps), vgrf12(ntyps),
704 & vgtr11(ntyps), vgtr12(ntyps),
705 & vgroca(ntyps), vgrotd(ntyps),
706 & vgrdrs(ntyps), vgz2 (ntyps)
707
708
709 data vgz0 /
710 1 2.6530, 2.6530, 2.6530, 2.6530, 2.6530, 2.6530, 2.6530,
711 1 2.6530, 2.6530, 2.6530, 2.6530, 2.6530,
712 2 0.5200, 0.5200, 0.6660, 0.9100, 1.0310, 1.0440, 1.0420,
713 2 1.0370, 1.0360, 0.9170, 0.6660, 0.5200,
714 3 1.1120, 1.1030, 1.0880, 1.0820, 1.0760, 1.0680, 1.0730,
715 3 1.0790, 1.0820, 1.0880, 1.1030, 1.1120,
716 4 0.0777, 0.0778, 0.0778, 0.0779, 0.0778, 0.0771, 0.0759,
717 4 0.0766, 0.0778, 0.0779, 0.0778, 0.0778,
718 5 0.2450, 0.2450, 0.2270, 0.2000, 0.2000, 0.2000, 0.2000,
719 5 0.267, 0.292, 0.280, 0.258, 0.2450,
720 6 0.0752, 0.0752, 0.0752, 0.0752, 0.0752, 0.0757, 0.0777,
721 6 0.0778, 0.0774, 0.0752, 0.0752, 0.0752,
722 7 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
723 7 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
724 8 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
725 8 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
726 9 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
727 9 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
728 1 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
729 1 0.0112, 0.0112, 0.0112, 0.0112, 0.0112
730 & /
731
732
733 data vgrd /
734 1 285.87, 285.87, 285.87, 285.87, 285.87, 285.87, 285.87,
735 1 285.87, 285.87, 285.87, 285.87, 285.87,
736 2 211.32, 211.32, 218.78, 243.40, 294.87, 345.90, 355.18,
737 2 341.84, 307.22, 244.84, 218.78, 211.32,
738 3 565.41, 587.05, 623.46, 638.13, 652.86, 675.04, 660.24,
739 3 645.49, 638.13, 623.46, 587.05, 565.41,
740 4 24.43, 24.63, 24.80, 24.96, 25.72, 27.74, 30.06,
741 4 28.86, 25.90, 25.11, 24.80, 24.63,
742 5 103.60, 103.60, 102.35, 100.72, 100.72, 100.72, 100.72,
743 5 105.30, 107.94, 106.59, 104.49, 103.60,
744 6 22.86, 22.86, 22.86, 22.86, 22.86, 23.01, 24.36,
745 6 24.69, 24.04, 22.86, 22.86, 22.86,
746 7 23.76, 23.76, 23.76, 23.76, 23.76, 23.76, 23.76,
747 7 23.76, 23.76, 23.76, 23.76, 23.76,
748 8 23.76, 23.76, 23.76, 23.76, 23.76, 23.76, 23.76,
749 8 23.76, 23.76, 23.76, 23.76, 23.76,
750 9 23.76, 23.76, 23.76, 23.76, 23.76, 23.76, 23.76,
751 9 23.76, 23.76, 23.76, 23.76, 23.76,
752 1 23.76, 23.76, 23.76, 23.76, 23.76, 23.76, 23.76,
753 1 23.76, 23.76, 23.76, 23.76, 23.76
754 & /
755
756 data vgrt /
757 1 19737.8, 19737.8, 19737.8, 19737.8, 19737.8, 19737.8, 19737.8,
758 1 19737.8, 19737.8, 19737.8, 19737.8, 19737.8,
759 2 5010.0, 5010.0, 5270.0, 6200.0, 8000.0, 9700.0, 9500.0,
760 2 8400.0, 6250.0, 5270.0, 5010.0, 5010.0,
761 3 9000.0, 9200.0, 9533.3, 9666.7, 9800.0, 9866.7, 9733.3,
762 3 9666.7, 9533.3, 9200.0, 9000.0, 9000.0,
763 4 5500.0, 5625.0, 5750.0, 5875.0, 6625.0, 8750.0, 9375.0,
764 4 6875.0, 6000.0, 5750.0, 5625.0, 5500.0,
765 5 6500.0, 6000.0, 5500.0, 5500.0, 5500.0, 5500.0, 5500.0,
766 5 7500.0, 8500.0, 7000.0, 6500.0, 6500.0,
767 6 10625.0, 10625.0, 10625.0, 10625.0, 10625.0, 11250.0, 18750.0,
768 6 17500.0, 10625.0, 10625.0, 10625.0, 10625.0,
769 7 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
770 7 1.0, 1.0, 1.0, 1.0, 1.0,
771 8 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
772 8 1.0, 1.0, 1.0, 1.0, 1.0,
773 9 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
774 9 1.0, 1.0, 1.0, 1.0, 1.0,
775 1 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
776 1 1.0, 1.0, 1.0, 1.0, 1.0
777 & /
778
779
780 data vgdd /
781 1 27.37, 27.37, 27.37, 27.37, 27.37, 27.37, 27.37,
782 1 27.37, 27.37, 27.37, 27.37, 27.37,
783 2 13.66, 13.66, 14.62, 15.70, 16.33, 16.62, 16.66,
784 2 16.60, 16.41, 15.73, 14.62, 13.66,
785 3 13.76, 13.80, 13.86, 13.88, 13.90, 13.93, 13.91,
786 3 13.89, 13.88, 13.86, 13.80, 13.76,
787 4 0.218, 0.227, 0.233, 0.239, 0.260, 0.299, 0.325,
788 4 0.313, 0.265, 0.244, 0.233, 0.227,
789 5 2.813, 2.813, 2.662, 2.391, 2.391, 2.391, 2.391,
790 5 2.975, 3.138, 3.062, 2.907, 2.813,
791 6 0.10629, 0.10629, 0.10629, 0.10629, 0.10629, 0.12299, 0.21521,
792 6 0.22897, 0.19961, 0.10629, 0.10629, 0.10629,
793 7 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
794 7 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
795 8 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
796 8 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
797 9 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
798 9 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
799 1 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
800 1 0.0001, 0.0001, 0.0001, 0.0001, 0.0001
801 & /
802
803
804 data vgrf11 /0.10,0.10,0.07,0.105,0.10,0.10,.001,.001,.001,.001/
805
806 data vgrf12 /0.16,0.16,0.16,0.360,0.16,0.16,.001,.001,.001,.001/
807
808 data vgtr11 /0.05,0.05,0.05,0.070,0.05,0.05,.001,.001,.001,.001/
809
810 data vgtr12 /.001,.001,.001, .220,.001,.001,.001,.001,.001,.001/
811
812 data vgroca /
813 & 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6,
814 & .1E-6, .1E-6, .1E-6, .1E-6 /
815
816 data vgrotd /1.00,1.00,0.50,0.50,0.50,0.20,0.10,0.10,0.10,0.10/
817
818 data vgrdrs /
819 & 0.75E13, 0.75E13, 0.75E13, 0.40E13, 0.75E13, 0.75E13,
820 & 0.10E13, 0.10E13, 0.10E13, 0.10E13 /
821
822 data vgz2 /35.0, 20.0, 17.0, 0.6, 5.0, 0.6, 0.1, 0.1, 0.1, 0.1/
823
824 vkrm = GETCON('VON KARMAN')
825
826 call time_bound ( nymd,nhms, nymd1,nhms1, nymd2,nhms2, k1,k2 )
827 call interp_time ( nymd,nhms, nymd1,nhms1, nymd2,nhms2, facm,facp )
828
829 do i=1,nchp
830
831 zoch(i) = vgz0(k2,ityp(i))*facp + vgz0(k1,ityp(i))*facm
832 rdcs(i) = vgrd(k2,ityp(i))*facp + vgrd(k1,ityp(i))*facm
833
834 rootl = vgrt(k2,ityp(i))*facp + vgrt(k1,ityp(i))*facm
835
836 vroot = rootl * vgroca(ityp (i))
837 dum1 = log (vroot / (1. - vroot))
838 dum2 = 1. / (8. * 3.14159 * rootl)
839 alphaf = dum2 * (vroot - 3. -2. * dum1)
840
841 rsl1(i) = vgrdrs (ityp (i)) / (rootl * vgrotd (ityp (i)))
842 rsl2(i) = alphaf / vgrotd (ityp (i))
843
844 scat = agrn(i) *(vgtr11(ityp(i)) + vgrf11(ityp(i)))
845 & + (1. - agrn(i))*(vgtr12(ityp(i)) + vgrf12(ityp(i)))
846 sqsc(i) = sqrt (1. - scat)
847
848 d = vgdd(k2,ityp(i))*facp + vgdd(k1,ityp(i))*facm
849 ufac(i) = log( (vgz2(ityp(i)) - d) / zoch(i) )
850 * / log( pblzet / zoch(i) )
851
852 z2ch(i) = vgz2(ityp (i))
853
854 cdsc(i) = pblzet/zoch(i)+1.
855 cdrc(i) = vkrm/log(cdsc(i))
856 cdrc(i) = cdrc(i)*cdrc(i)
857 cdsc(i) = sqrt(cdsc(i))
858 cdsc(i) = cdrc(i)*cdsc(i)
859
860 enddo
861
862 return
863 end
864
865 subroutine pkappa (ple,pkle,im,jm,lm,pkz)
866 C***********************************************************************
867 C Purpose
868 C Calculate Phillips P**Kappa
869 C
870 C Arguments
871 C PLE .... edge-level pressure
872 C PKLE ... edge-level pressure**kappa
873 C IM ..... longitude dimension
874 C JM ..... latitude dimension
875 C LM ..... vertical dimension
876 C PKZ .... mid-level pressure**kappa
877 C***********************************************************************
878 implicit none
879
880 integer im,jm,lm
881 real ple (im,jm,lm+1)
882 real pkle(im,jm,lm+1)
883 real pkz (im,jm,lm)
884
885 real akap1,getcon
886 integer i,j,L
887
888 akap1 = 1.0 + getcon('KAPPA')
889
890 do L = 1,lm
891 do j = 1,jm
892 do i = 1,im
893 pkz(i,j,L) = ( ple (i,j,l+1)*pkle(i,j,l+1)
894 . - ple (i,j,l)*pkle(i,j,l) )
895 . / ( akap1* (ple (i,j,l+1)-ple (i,j,l)) )
896 enddo
897 enddo
898 enddo
899
900 return
901 end

  ViewVC Help
Powered by ViewVC 1.1.22