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

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

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


Revision 1.5 - (hide 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 molod 1.1 function minval (q,im)
2     implicit none
3 molod 1.3 #include "CPP_EEOPTIONS.h"
4 molod 1.1 integer im, i
5 molod 1.2 _RL q(im), minval
6 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
54 molod 1.2 _RL A(IA,L), B(IB,L)
55 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
81 molod 1.2 _RL A(IA,L), B(IB,L)
82 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
97 molod 1.2 _RL A(IA,L), B(IB,L)
98 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
162 molod 1.2 _RL TT, P, Q, DQDT
163 molod 1.1 LOGICAL LDQDT
164 molod 1.2 _RL AIRMW, H2OMW
165 molod 1.1
166     PARAMETER ( AIRMW = 28.97 )
167     PARAMETER ( H2OMW = 18.01 )
168    
169 molod 1.2 _RL ESFAC, ERFAC
170 molod 1.1 PARAMETER ( ESFAC = H2OMW/AIRMW )
171     PARAMETER ( ERFAC = (1.0-ESFAC)/ESFAC )
172    
173 molod 1.2 _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 molod 1.1
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 molod 1.2 _RL TMAX, TICE
272 molod 1.1 PARAMETER ( TMAX=323.15, TICE=273.16)
273    
274 molod 1.2 _RL T, D, W, QX, DQX
275 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
334 molod 1.1 integer i,n
335     logical ldqdt
336 molod 1.2 _RL tt(n), p(n), q(n), dqdt(n)
337 molod 1.1 #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 molod 1.3 #include "CPP_EEOPTIONS.h"
351 molod 1.1 integer ia,ib,irun,l,k
352 molod 1.2 _RL a(ia,l), b(ib,l)
353 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
380 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
408 molod 1.2 _RL a(ia,l), b(ib,l)
409 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
426 molod 1.2 _RL a(ia,l)
427 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
446 molod 1.2 _RL a(ia,l), b(ib,l)
447 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
482 molod 1.1 integer ia,ib,irun,levs,npeice
483 molod 1.2 _RL a(ia,levs), b(ib,levs)
484 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
537 molod 1.1 integer ia,ib,levs,numpts,npeice
538     integer index(numpts)
539 molod 1.2 _RL a(ia,levs), b(ib,levs), chfr(ib)
540 molod 1.5 real b(ib,levs)
541 molod 1.1
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 molod 1.3 #include "CPP_EEOPTIONS.h"
584 molod 1.1 integer ia,ib,levs,numpts,npeice
585     integer index(numpts)
586 molod 1.2 _RL a(ia,levs), b(ib,levs), chfr(ib)
587 molod 1.1 logical check
588    
589     integer i,L,offset,len
590 molod 1.2 _RL undef,getcon
591 molod 1.1
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 molod 1.3 #include "CPP_EEOPTIONS.h"
621 molod 1.1 integer im,jm,mxchps,nchp
622     integer igrd(mxchps)
623 molod 1.2 _RL A(IM,JM), B(MXCHPS)
624 molod 1.1
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 molod 1.3 #include "CPP_EEOPTIONS.h"
642 molod 1.2 _RL zero,one
643 molod 1.1 parameter ( one = 1.)
644     parameter (zero = 0.)
645     integer im,jm,mxchps,nchp
646     integer igrd(mxchps)
647 molod 1.2 _RL A(IM,JM), B(MXCHPS), CHFR(MXCHPS), FRACG(IM,JM)
648 molod 1.1
649 molod 1.2 _RL VT1(IM,JM)
650 molod 1.1 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 molod 1.3 #include "CPP_EEOPTIONS.h"
676 molod 1.1
677     integer nymd,nhms,nchp,mxchps,ityp(mxchps)
678 molod 1.2 _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 molod 1.1
684     C*********************************************************************
685     C********************* SUBROUTINE CHPPRM ****************************
686     C********************** 14 JUNE 1991 ******************************
687     C*********************************************************************
688    
689     integer ntyps
690     parameter (ntyps=10)
691    
692 molod 1.2 _RL pblzet
693 molod 1.1 parameter (pblzet = 50.)
694     integer k1,k2,nymd1,nhms1,nymd2,nhms2,i
695 molod 1.2 _RL getcon,vkrm,rootl,vroot,dum1,dum2,alphaf
696     _RL facm,facp
697     _RL scat,d
698 molod 1.1
699 molod 1.2 _RL
700 molod 1.1 & 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 molod 1.4
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