/[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.12 - (hide annotations) (download)
Wed Jul 28 20:38:53 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
Changes since 1.11: +2 -2 lines
debugging

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

  ViewVC Help
Powered by ViewVC 1.1.22