/[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.2 - (hide annotations) (download)
Thu Jun 10 20:17:17 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
Changes since 1.1: +58 -130 lines
Developing

1 molod 1.1 function minval (q,im)
2     implicit none
3 molod 1.2 #include "CPP_OPTIONS.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.2 #include "CPP_OPTIONS.h"
54     _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.2 #include "CPP_OPTIONS.h"
81     _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.2 #include "CPP_OPTIONS.h"
97     _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.2 #include "CPP_OPTIONS.h"
162     _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.2 #include "CPP_OPTIONS.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.2 #include "CPP_OPTIONS.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.2 #include "CPP_OPTIONS.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.2 #include "CPP_OPTIONS.h"
408     _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.2 #include "CPP_OPTIONS.h"
426     _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.2 #include "CPP_OPTIONS.h"
446     _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.2 #include "CPP_OPTIONS.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.2 #include "CPP_OPTIONS.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.1
541     integer i,L,offset,len
542    
543     offset = ib*(npeice-1)
544     len = min(ib,numpts-offset)
545     offset = offset+1
546    
547     do L = 1,levs
548     do i=1,len
549     a(index(i+offset-1),L) = a(index(i+offset-1),L) + b(i,L)*chfr(i)
550     enddo
551     enddo
552     return
553     end
554     subroutine paste2grd (b,index,chfr,ib,numpts,a,ia,levs,npeice,check)
555     c-----------------------------------------------------------------------
556     c subroutine paste2grd - paste one processors worth of grid points
557     c from a stripped tile array to a grid
558     c space array
559     c
560     c input:
561     c b - array to be pasted back into grid space [ib,levs]
562     c index - array of horizontal indeces of grid points to convert to
563     c tile space[numpts]
564     c chfr - fractional area covered by the tile [ib]
565     c ib - inner dimension of source array AND number of points in
566     c array a that need to be pasted
567     c numpts - total number of points which were stripped
568     c ia - inner of dimension of target array
569     c levs - number of vertical levels AND outer dimension of arrays
570     c npeice - the current strip number to be filled
571     c check - logical to check for undefined values
572     c output:
573     c a - grid space array to be filled [ia,levs]
574     c
575     c IMPORTANT NOTE:
576     c
577     c This routine will result in roundoff differences if called from
578     c within a parallel region.
579     c-----------------------------------------------------------------------
580    
581     implicit none
582 molod 1.2 #include "CPP_OPTIONS.h"
583 molod 1.1 integer ia,ib,levs,numpts,npeice
584     integer index(numpts)
585 molod 1.2 _RL a(ia,levs), b(ib,levs), chfr(ib)
586 molod 1.1 logical check
587    
588     integer i,L,offset,len
589 molod 1.2 _RL undef,getcon
590 molod 1.1
591     offset = ib*(npeice-1)
592     len = min(ib,numpts-offset)
593     offset = offset+1
594    
595     if( check ) then
596     undef = getcon('UNDEF')
597     do L= 1,levs
598     do i= 1,len
599     if( a(index(i+offset-1),L).eq.undef .or. b(i,L).eq.undef ) then
600     a(index(i+offset-1),L) = undef
601     else
602     a(index(i+offset-1),L) = a(index(i+offset-1),L) + b(i,L)*chfr(i)
603     endif
604     enddo
605     enddo
606     else
607     do L= 1,levs
608     do i= 1,len
609     a(index(i+offset-1),L) = a(index(i+offset-1),L) + b(i,L)*chfr(i)
610     enddo
611     enddo
612     endif
613    
614     return
615     end
616     SUBROUTINE GRD2MSC(A,IM,JM,IGRD,B,MXCHPS,NCHP)
617    
618     implicit none
619 molod 1.2 #include "CPP_OPTIONS.h"
620 molod 1.1 integer im,jm,mxchps,nchp
621     integer igrd(mxchps)
622 molod 1.2 _RL A(IM,JM), B(MXCHPS)
623 molod 1.1
624     integer i
625    
626     IF(NCHP.GE.0) THEN
627     DO I=1,NCHP
628     B(I) = A(IGRD(I),1)
629     ENDDO
630     ELSE
631     PRINT *, 'ERROR IN GRD2MSC'
632     ENDIF
633    
634     RETURN
635     END
636    
637     SUBROUTINE MSC2GRD(IGRD,CHFR,B,MXCHPS,NCHP,FRACG,A,IM,JM)
638    
639     implicit none
640 molod 1.2 #include "CPP_OPTIONS.h"
641     _RL zero,one
642 molod 1.1 parameter ( one = 1.)
643     parameter (zero = 0.)
644     integer im,jm,mxchps,nchp
645     integer igrd(mxchps)
646 molod 1.2 _RL A(IM,JM), B(MXCHPS), CHFR(MXCHPS), FRACG(IM,JM)
647 molod 1.1
648 molod 1.2 _RL VT1(IM,JM)
649 molod 1.1 integer i
650    
651     IF(NCHP.GE.0) THEN
652     DO I=1,IM*JM
653     VT1(I,1) = ZERO
654     ENDDO
655    
656     DO I=1,NCHP
657     VT1(IGRD(I),1) = VT1(IGRD(I),1) + B(I)*CHFR(I)
658     ENDDO
659    
660     DO I=1,IM*JM
661     A(I,1) = A(I,1)*(ONE-FRACG(I,1)) + VT1(I,1)
662     ENDDO
663     ELSE
664     PRINT *, 'ERROR IN MSC2GRD'
665     ENDIF
666    
667     RETURN
668     END
669    
670     subroutine chpprm(nymd,nhms,mxchps,nchp,chlt,ityp,alai,
671     1 agrn,zoch,z2ch,cdrc,cdsc,sqsc,ufac,rsl1,rsl2,rdcs)
672    
673     implicit none
674 molod 1.2 #include "CPP_OPTIONS.h"
675 molod 1.1
676     integer nymd,nhms,nchp,mxchps,ityp(mxchps)
677 molod 1.2 _RL chlt(mxchps)
678     _RL alai(mxchps),agrn(mxchps)
679     _RL zoch(mxchps), z2ch(mxchps), cdrc(mxchps), cdsc(mxchps)
680     _RL sqsc(mxchps), ufac(mxchps), rsl1(mxchps), rsl2(mxchps)
681     _RL rdcs(mxchps)
682 molod 1.1
683     C*********************************************************************
684     C********************* SUBROUTINE CHPPRM ****************************
685     C********************** 14 JUNE 1991 ******************************
686     C*********************************************************************
687    
688     integer ntyps
689     parameter (ntyps=10)
690    
691 molod 1.2 _RL pblzet
692 molod 1.1 parameter (pblzet = 50.)
693     integer k1,k2,nymd1,nhms1,nymd2,nhms2,i
694 molod 1.2 _RL getcon,vkrm,rootl,vroot,dum1,dum2,alphaf
695     _RL facm,facp
696     _RL scat,d
697 molod 1.1
698 molod 1.2 _RL
699 molod 1.1 & vgdd(12, ntyps), vgz0(12, ntyps),
700     & vgrd(12, ntyps), vgrt(12, ntyps),
701    
702     & vgrf11(ntyps), vgrf12(ntyps),
703     & vgtr11(ntyps), vgtr12(ntyps),
704     & vgroca(ntyps), vgrotd(ntyps),
705     & vgrdrs(ntyps), vgz2 (ntyps)
706    
707    
708     data vgz0 /
709     1 2.6530, 2.6530, 2.6530, 2.6530, 2.6530, 2.6530, 2.6530,
710     1 2.6530, 2.6530, 2.6530, 2.6530, 2.6530,
711     2 0.5200, 0.5200, 0.6660, 0.9100, 1.0310, 1.0440, 1.0420,
712     2 1.0370, 1.0360, 0.9170, 0.6660, 0.5200,
713     3 1.1120, 1.1030, 1.0880, 1.0820, 1.0760, 1.0680, 1.0730,
714     3 1.0790, 1.0820, 1.0880, 1.1030, 1.1120,
715     4 0.0777, 0.0778, 0.0778, 0.0779, 0.0778, 0.0771, 0.0759,
716     4 0.0766, 0.0778, 0.0779, 0.0778, 0.0778,
717     5 0.2450, 0.2450, 0.2270, 0.2000, 0.2000, 0.2000, 0.2000,
718     5 0.267, 0.292, 0.280, 0.258, 0.2450,
719     6 0.0752, 0.0752, 0.0752, 0.0752, 0.0752, 0.0757, 0.0777,
720     6 0.0778, 0.0774, 0.0752, 0.0752, 0.0752,
721     7 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
722     7 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
723     8 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
724     8 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
725     9 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
726     9 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
727     1 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
728     1 0.0112, 0.0112, 0.0112, 0.0112, 0.0112
729     & /
730    
731    
732     data vgrd /
733     1 285.87, 285.87, 285.87, 285.87, 285.87, 285.87, 285.87,
734     1 285.87, 285.87, 285.87, 285.87, 285.87,
735     2 211.32, 211.32, 218.78, 243.40, 294.87, 345.90, 355.18,
736     2 341.84, 307.22, 244.84, 218.78, 211.32,
737     3 565.41, 587.05, 623.46, 638.13, 652.86, 675.04, 660.24,
738     3 645.49, 638.13, 623.46, 587.05, 565.41,
739     4 24.43, 24.63, 24.80, 24.96, 25.72, 27.74, 30.06,
740     4 28.86, 25.90, 25.11, 24.80, 24.63,
741     5 103.60, 103.60, 102.35, 100.72, 100.72, 100.72, 100.72,
742     5 105.30, 107.94, 106.59, 104.49, 103.60,
743     6 22.86, 22.86, 22.86, 22.86, 22.86, 23.01, 24.36,
744     6 24.69, 24.04, 22.86, 22.86, 22.86,
745     7 23.76, 23.76, 23.76, 23.76, 23.76, 23.76, 23.76,
746     7 23.76, 23.76, 23.76, 23.76, 23.76,
747     8 23.76, 23.76, 23.76, 23.76, 23.76, 23.76, 23.76,
748     8 23.76, 23.76, 23.76, 23.76, 23.76,
749     9 23.76, 23.76, 23.76, 23.76, 23.76, 23.76, 23.76,
750     9 23.76, 23.76, 23.76, 23.76, 23.76,
751     1 23.76, 23.76, 23.76, 23.76, 23.76, 23.76, 23.76,
752     1 23.76, 23.76, 23.76, 23.76, 23.76
753     & /
754    
755     data vgrt /
756     1 19737.8, 19737.8, 19737.8, 19737.8, 19737.8, 19737.8, 19737.8,
757     1 19737.8, 19737.8, 19737.8, 19737.8, 19737.8,
758     2 5010.0, 5010.0, 5270.0, 6200.0, 8000.0, 9700.0, 9500.0,
759     2 8400.0, 6250.0, 5270.0, 5010.0, 5010.0,
760     3 9000.0, 9200.0, 9533.3, 9666.7, 9800.0, 9866.7, 9733.3,
761     3 9666.7, 9533.3, 9200.0, 9000.0, 9000.0,
762     4 5500.0, 5625.0, 5750.0, 5875.0, 6625.0, 8750.0, 9375.0,
763     4 6875.0, 6000.0, 5750.0, 5625.0, 5500.0,
764     5 6500.0, 6000.0, 5500.0, 5500.0, 5500.0, 5500.0, 5500.0,
765     5 7500.0, 8500.0, 7000.0, 6500.0, 6500.0,
766     6 10625.0, 10625.0, 10625.0, 10625.0, 10625.0, 11250.0, 18750.0,
767     6 17500.0, 10625.0, 10625.0, 10625.0, 10625.0,
768     7 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
769     7 1.0, 1.0, 1.0, 1.0, 1.0,
770     8 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
771     8 1.0, 1.0, 1.0, 1.0, 1.0,
772     9 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
773     9 1.0, 1.0, 1.0, 1.0, 1.0,
774     1 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
775     1 1.0, 1.0, 1.0, 1.0, 1.0
776     & /
777    
778    
779     data vgdd /
780     1 27.37, 27.37, 27.37, 27.37, 27.37, 27.37, 27.37,
781     1 27.37, 27.37, 27.37, 27.37, 27.37,
782     2 13.66, 13.66, 14.62, 15.70, 16.33, 16.62, 16.66,
783     2 16.60, 16.41, 15.73, 14.62, 13.66,
784     3 13.76, 13.80, 13.86, 13.88, 13.90, 13.93, 13.91,
785     3 13.89, 13.88, 13.86, 13.80, 13.76,
786     4 0.218, 0.227, 0.233, 0.239, 0.260, 0.299, 0.325,
787     4 0.313, 0.265, 0.244, 0.233, 0.227,
788     5 2.813, 2.813, 2.662, 2.391, 2.391, 2.391, 2.391,
789     5 2.975, 3.138, 3.062, 2.907, 2.813,
790     6 0.10629, 0.10629, 0.10629, 0.10629, 0.10629, 0.12299, 0.21521,
791     6 0.22897, 0.19961, 0.10629, 0.10629, 0.10629,
792     7 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
793     7 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
794     8 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
795     8 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
796     9 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
797     9 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
798     1 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
799     1 0.0001, 0.0001, 0.0001, 0.0001, 0.0001
800     & /
801    
802    
803     data vgrf11 /0.10,0.10,0.07,0.105,0.10,0.10,.001,.001,.001,.001/
804    
805     data vgrf12 /0.16,0.16,0.16,0.360,0.16,0.16,.001,.001,.001,.001/
806    
807     data vgtr11 /0.05,0.05,0.05,0.070,0.05,0.05,.001,.001,.001,.001/
808    
809     data vgtr12 /.001,.001,.001, .220,.001,.001,.001,.001,.001,.001/
810    
811     data vgroca /
812     & 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6,
813     & .1E-6, .1E-6, .1E-6, .1E-6 /
814    
815     data vgrotd /1.00,1.00,0.50,0.50,0.50,0.20,0.10,0.10,0.10,0.10/
816    
817     data vgrdrs /
818     & 0.75E13, 0.75E13, 0.75E13, 0.40E13, 0.75E13, 0.75E13,
819     & 0.10E13, 0.10E13, 0.10E13, 0.10E13 /
820    
821     data vgz2 /35.0, 20.0, 17.0, 0.6, 5.0, 0.6, 0.1, 0.1, 0.1, 0.1/
822    
823     vkrm = GETCON('VON KARMAN')
824    
825     call time_bound ( nymd,nhms, nymd1,nhms1, nymd2,nhms2, k1,k2 )
826     call interp_time ( nymd,nhms, nymd1,nhms1, nymd2,nhms2, facm,facp )
827    
828     do i=1,nchp
829    
830     zoch(i) = vgz0(k2,ityp(i))*facp + vgz0(k1,ityp(i))*facm
831     rdcs(i) = vgrd(k2,ityp(i))*facp + vgrd(k1,ityp(i))*facm
832    
833     rootl = vgrt(k2,ityp(i))*facp + vgrt(k1,ityp(i))*facm
834    
835     vroot = rootl * vgroca(ityp (i))
836     dum1 = log (vroot / (1. - vroot))
837     dum2 = 1. / (8. * 3.14159 * rootl)
838     alphaf = dum2 * (vroot - 3. -2. * dum1)
839    
840     rsl1(i) = vgrdrs (ityp (i)) / (rootl * vgrotd (ityp (i)))
841     rsl2(i) = alphaf / vgrotd (ityp (i))
842    
843     scat = agrn(i) *(vgtr11(ityp(i)) + vgrf11(ityp(i)))
844     & + (1. - agrn(i))*(vgtr12(ityp(i)) + vgrf12(ityp(i)))
845     sqsc(i) = sqrt (1. - scat)
846    
847     d = vgdd(k2,ityp(i))*facp + vgdd(k1,ityp(i))*facm
848     ufac(i) = log( (vgz2(ityp(i)) - d) / zoch(i) )
849     * / log( pblzet / zoch(i) )
850    
851     z2ch(i) = vgz2(ityp (i))
852    
853     cdsc(i) = pblzet/zoch(i)+1.
854     cdrc(i) = vkrm/log(cdsc(i))
855     cdrc(i) = cdrc(i)*cdrc(i)
856     cdsc(i) = sqrt(cdsc(i))
857     cdsc(i) = cdrc(i)*cdsc(i)
858    
859     enddo
860    
861     return
862     end

  ViewVC Help
Powered by ViewVC 1.1.22