/[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.6 - (show annotations) (download)
Fri Jul 16 16:11:36 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
Changes since 1.5: +82 -66 lines
Debugging

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

  ViewVC Help
Powered by ViewVC 1.1.22