source: trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F @ 1032

Last change on this file since 1032 was 787, checked in by aslmd, 13 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

  • nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this
File size: 19.6 KB
Line 
1      SUBROUTINE vdif_kc(ngrid,dt,g,zlev,zlay,u,v,teta,cd,q2,km,kn)
2      IMPLICIT NONE
3c.......................................................................
4#include "dimensions.h"
5#include "dimphys.h"
6c.......................................................................
7c
8c dt : pas de temps
9c g  : g
10c zlev : altitude a chaque niveau (interface inferieure de la couche
11c        de meme indice)
12c zlay : altitude au centre de chaque couche
13c u,v : vitesse au centre de chaque couche
14c       (en entree : la valeur au debut du pas de temps)
15c teta : temperature potentielle au centre de chaque couche
16c        (en entree : la valeur au debut du pas de temps)
17c cd : cdrag
18c      (en entree : la valeur au debut du pas de temps)
19c q2 : $q^2$ au bas de chaque couche
20c      (en entree : la valeur au debut du pas de temps)
21c      (en sortie : la valeur a la fin du pas de temps)
22c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
23c      couche)
24c      (en sortie : la valeur a la fin du pas de temps)
25c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
26c      (en sortie : la valeur a la fin du pas de temps)
27c
28c.......................................................................
29      INTEGER ngrid
30      REAL dt,g
31      REAL zlev(ngrid,nlayermx+1)
32      REAL zlay(ngrid,nlayermx)
33      REAL u(ngrid,nlayermx)
34      REAL v(ngrid,nlayermx)
35      REAL teta(ngrid,nlayermx)
36      REAL cd(ngrid)
37      REAL q2(ngrid,nlayermx+1)
38      REAL km(ngrid,nlayermx+1)
39      REAL kn(ngrid,nlayermx+1)
40c.......................................................................
41c
42c nlay : nombre de couches       
43c nlev : nombre de niveaux
44c ngrid : nombre de points de grille       
45c unsdz : 1 sur l'epaisseur de couche
46c unsdzdec : 1 sur la distance entre le centre de la couche et le
47c            centre de la couche inferieure
48c q : echelle de vitesse au bas de chaque couche
49c     (valeur a la fin du pas de temps)
50c
51c.......................................................................
52      INTEGER nlay,nlev
53      REAL unsdz(ngrid,nlayermx)
54      REAL unsdzdec(ngrid,nlayermx+1)
55      REAL q(ngrid,nlayermx+1)
56c.......................................................................
57c
58c kmpre : km au debut du pas de temps
59c qcstat : q : solution stationnaire du probleme couple
60c          (valeur a la fin du pas de temps)
61c q2cstat : q2 : solution stationnaire du probleme couple
62c           (valeur a la fin du pas de temps)
63c
64c.......................................................................
65      REAL kmpre(ngrid,nlayermx+1)
66      REAL qcstat
67      REAL q2cstat
68c.......................................................................
69c
70c long : longueur de melange calculee selon Blackadar
71c
72c.......................................................................
73      REAL long(ngrid,nlayermx+1)
74c.......................................................................
75c
76c kmq3 : terme en q^3 dans le developpement de km
77c        (valeur au debut du pas de temps)
78c kmcstat : valeur de km solution stationnaire du systeme {q2 ; du/dz}
79c           (valeur a la fin du pas de temps)
80c knq3 : terme en q^3 dans le developpement de kn
81c mcstat : valeur de m solution stationnaire du systeme {q2 ; du/dz}
82c          (valeur a la fin du pas de temps)
83c m2cstat : valeur de m2 solution stationnaire du systeme {q2 ; du/dz}
84c           (valeur a la fin du pas de temps)
85c m : valeur a la fin du pas de temps
86c mpre : valeur au debut du pas de temps
87c m2 : valeur a la fin du pas de temps
88c n2 : valeur a la fin du pas de temps
89c
90c.......................................................................
91      REAL kmq3
92      REAL kmcstat
93      REAL knq3
94      REAL mcstat
95      REAL m2cstat
96      REAL m(ngrid,nlayermx+1)
97      REAL mpre(ngrid,nlayermx+1)
98      REAL m2(ngrid,nlayermx+1)
99      REAL n2(ngrid,nlayermx+1)
100c.......................................................................
101c
102c gn : intermediaire pour les coefficients de stabilite
103c gnmin : borne inferieure de gn (-0.23 ou -0.28)
104c gnmax : borne superieure de gn (0.0233)
105c gninf : vrai si gn est en dessous de sa borne inferieure
106c gnsup : vrai si gn est en dessus de sa borne superieure
107c gm : drole d'objet bien utile
108c ri : nombre de Richardson
109c sn : coefficient de stabilite pour n
110c snq2 : premier terme du developement limite de sn en q2
111c sm : coefficient de stabilite pour m
112c smq2 : premier terme du developement limite de sm en q2
113c
114c.......................................................................
115      REAL gn
116      REAL gnmin
117      REAL gnmax
118      LOGICAL gninf
119      LOGICAL gnsup
120      REAL gm
121c      REAL ri(ngrid,nlayermx+1)
122      REAL sn(ngrid,nlayermx+1)
123      REAL snq2(ngrid,nlayermx+1)
124      REAL sm(ngrid,nlayermx+1)
125      REAL smq2(ngrid,nlayermx+1)
126c.......................................................................
127c
128c kappa : consatnte de Von Karman (0.4)
129c long0 : longueur de reference pour le calcul de long (160)
130c a1,a2,b1,b2,c1 : constantes d'origine pour les  coefficients
131c                  de stabilite (0.92/0.74/16.6/10.1/0.08)
132c cn1,cn2 : constantes pour sn
133c cm1,cm2,cm3,cm4 : constantes pour sm
134c
135c.......................................................................
136      REAL kappa
137      REAL long0
138      REAL a1,a2,b1,b2,c1
139      REAL cn1,cn2
140      REAL cm1,cm2,cm3,cm4
141c.......................................................................
142c
143c termq : termes en $q$ dans l'equation de q2
144c termq3 : termes en $q^3$ dans l'equation de q2
145c termqm2 : termes en $q*m^2$ dans l'equation de q2
146c termq3m2 : termes en $q^3*m^2$ dans l'equation de q2
147c
148c.......................................................................
149      REAL termq
150      REAL termq3
151      REAL termqm2
152      REAL termq3m2
153c.......................................................................
154c
155c q2min : borne inferieure de q2
156c q2max : borne superieure de q2
157c
158c.......................................................................
159      REAL q2min
160      REAL q2max
161c.......................................................................
162c knmin : borne inferieure de kn
163c kmmin : borne inferieure de km
164c.......................................................................
165      REAL knmin
166      REAL kmmin
167c.......................................................................
168      INTEGER ilay,ilev,igrid
169      REAL tmp1,tmp2
170c.......................................................................
171      PARAMETER (kappa=0.4E+0)
172      PARAMETER (long0=160.E+0)
173      PARAMETER (gnmin=-10.E+0)
174      PARAMETER (gnmax=0.0233E+0)
175      PARAMETER (a1=0.92E+0)
176      PARAMETER (a2=0.74E+0)
177      PARAMETER (b1=16.6E+0)
178      PARAMETER (b2=10.1E+0)
179      PARAMETER (c1=0.08E+0)
180      PARAMETER (knmin=1.E-5)
181      PARAMETER (kmmin=1.E-5)
182      PARAMETER (q2min=1.E-3)
183      PARAMETER (q2max=1.E+2)
184      PARAMETER (nlay=nlayermx)
185      PARAMETER (nlev=nlayermx+1)
186c
187      PARAMETER (
188     &  cn1=a2*(1.E+0 -6.E+0 *a1/b1)
189     &          )
190      PARAMETER (
191     &  cn2=-3.E+0 *a2*(6.E+0 *a1+b2)
192     &          )
193      PARAMETER (
194     &  cm1=a1*(1.E+0 -3.E+0 *c1-6.E+0 *a1/b1)
195     &          )
196      PARAMETER (
197     &  cm2=a1*(-3.E+0 *a2*((b2-3.E+0 *a2)*(1.E+0 -6.E+0 *a1/b1)
198     &          -3.E+0 *c1*(b2+6.E+0 *a1)))
199     &          )
200      PARAMETER (
201     &  cm3=-3.E+0 *a2*(6.E+0 *a1+b2)
202     &          )
203      PARAMETER (
204     &  cm4=-9.E+0 *a1*a2
205     &          )
206c.......................................................................
207c  traitment des valeur de q2 en entree
208c.......................................................................
209c
210      DO ilev=1,nlev
211                                                 DO igrid=1,ngrid
212        q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
213        q(igrid,ilev)=sqrt(q2(igrid,ilev))
214                                                 ENDDO
215      ENDDO
216c
217                                                 DO igrid=1,ngrid
218      tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
219      q2(igrid,1)=b1**(2.E+0/3.E+0)*tmp1
220      q2(igrid,1)=amax1(q2(igrid,1),q2min)
221      q(igrid,1)=sqrt(q2(igrid,1))
222                                                 ENDDO
223c
224c.......................................................................
225c  les increments verticaux
226c.......................................................................
227c
228c!!!!! allerte !!!!!c
229c!!!!! zlev n'est pas declare a nlev !!!!!c
230c!!!!! ---->
231c                                                     DO igrid=1,ngrid
232c           zlev(igrid,nlev)=zlay(igrid,nlay)
233c    &             +( zlay(igrid,nlay) - zlev(igrid,nlev-1) )
234c                                                     ENDDO           
235c!!!!! <----
236c!!!!! allerte !!!!!c
237c
238      DO ilay=1,nlay
239                                                 DO igrid=1,ngrid
240        unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
241                                                 ENDDO
242      ENDDO
243                                                 DO igrid=1,ngrid
244      unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
245                                                 ENDDO
246      DO ilay=2,nlay
247                                                 DO igrid=1,ngrid
248        unsdzdec(igrid,ilay)=1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
249                                                 ENDDO
250      ENDDO
251                                                 DO igrid=1,ngrid
252      unsdzdec(igrid,nlay+1)=1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
253                                                 ENDDO
254c
255c.......................................................................
256c  le cisaillement et le gradient de temperature
257c.......................................................................
258c
259                                                 DO igrid=1,ngrid
260      m2(igrid,1)=(unsdzdec(igrid,1)
261     &                   *u(igrid,1))**2
262     &                 +(unsdzdec(igrid,1)
263     &                   *v(igrid,1))**2
264      m(igrid,1)=sqrt(m2(igrid,1))
265      mpre(igrid,1)=m(igrid,1)
266                                                 ENDDO
267c
268c-----------------------------------------------------------------------
269      DO ilev=2,nlev-1
270                                                 DO igrid=1,ngrid
271c-----------------------------------------------------------------------
272c
273        n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
274     &                   *(teta(igrid,ilev)-teta(igrid,ilev-1))
275     &                   /(teta(igrid,ilev)+teta(igrid,ilev-1)) *2.E+0
276c
277c --->
278c       on ne sais traiter que les cas stratifies. et l'ajustement
279c       convectif est cense faire en sorte que seul des configurations
280c       stratifiees soient rencontrees en entree de cette routine.
281c       mais, bon ... on sait jamais (meme on sait que n2 prends
282c       quelques valeurs negatives ... parfois) alors :
283c<---
284c
285        IF (n2(igrid,ilev).lt.0.E+0) THEN
286          n2(igrid,ilev)=0.E+0
287        ENDIF
288c
289        m2(igrid,ilev)=(unsdzdec(igrid,ilev)
290     &                     *(u(igrid,ilev)-u(igrid,ilev-1)))**2
291     &                   +(unsdzdec(igrid,ilev)
292     &                     *(v(igrid,ilev)-v(igrid,ilev-1)))**2
293        m(igrid,ilev)=sqrt(m2(igrid,ilev))
294        mpre(igrid,ilev)=m(igrid,ilev)
295c
296c-----------------------------------------------------------------------
297                                                 ENDDO
298      ENDDO
299c-----------------------------------------------------------------------
300c
301                                                 DO igrid=1,ngrid
302      m2(igrid,nlev)=m2(igrid,nlev-1)
303      m(igrid,nlev)=m(igrid,nlev-1)
304      mpre(igrid,nlev)=m(igrid,nlev)
305                                                 ENDDO
306c
307c.......................................................................
308c  calcul des fonctions de stabilite
309c.......................................................................
310c
311c-----------------------------------------------------------------------
312      DO ilev=2,nlev-1
313                                                 DO igrid=1,ngrid
314c-----------------------------------------------------------------------
315c
316        tmp1=kappa*(zlev(igrid,ilev)-zlev(igrid,1))
317        long(igrid,ilev)=tmp1/(1.E+0 + tmp1/long0)
318        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
319     &                                           * n2(igrid,ilev)
320        gm=long(igrid,ilev)**2 / q2(igrid,ilev)
321     &                                           * m2(igrid,ilev)
322c
323        gninf=.false.
324        gnsup=.false.
325        long(igrid,ilev)=long(igrid,ilev)
326        long(igrid,ilev)=long(igrid,ilev)
327c
328        IF (gn.lt.gnmin) THEN
329          gninf=.true.
330          gn=gnmin
331        ENDIF
332c
333        IF (gn.gt.gnmax) THEN
334          gnsup=.true.
335          gn=gnmax
336        ENDIF
337c
338        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
339        sm(igrid,ilev)=
340     &    (cm1+cm2*gn)
341     &   /( (1.E+0 +cm3*gn)
342     &     *(1.E+0 +cm4*gn) )
343c
344        IF ((gninf).or.(gnsup)) THEN
345          snq2(igrid,ilev)=0.E+0
346          smq2(igrid,ilev)=0.E+0
347        ELSE
348          snq2(igrid,ilev)=
349     &     -gn
350     &     *(-cn1*cn2/(1.E+0 +cn2*gn)**2 )
351          smq2(igrid,ilev)=
352     &     -gn
353     &     *( cm2*(1.E+0 +cm3*gn)
354     &           *(1.E+0 +cm4*gn)
355     &       -( cm3*(1.E+0 +cm4*gn)
356     &         +cm4*(1.E+0 +cm3*gn) )
357     &       *(cm1+cm2*gn)            )
358     &     /( (1.E+0 +cm3*gn)
359     &       *(1.E+0 +cm4*gn) )**2
360        ENDIF
361c
362c --->
363c       la decomposition de Taylor en q2 n'a de sens que
364c       dans les cas stratifies ou sn et sm sont quasi
365c       proportionnels a q2. ailleurs on laisse le meme
366c       algorithme car l'ajustement convectif fait le travail.
367c       mais c'est delirant quand sn et snq2 n'ont pas le meme
368c       signe : dans ces cas, on ne fait pas la decomposition.
369c<---
370c
371        IF (snq2(igrid,ilev)*sn(igrid,ilev).le.0.E+0)
372     &      snq2(igrid,ilev)=0.E+0
373        IF (smq2(igrid,ilev)*sm(igrid,ilev).le.0.E+0)
374     &      smq2(igrid,ilev)=0.E+0
375c
376c-----------------------------------------------------------------------
377                                                  ENDDO
378      ENDDO
379c-----------------------------------------------------------------------
380c
381c.......................................................................
382c  calcul de km et kn au debut du pas de temps
383c.......................................................................
384c
385                                                 DO igrid=1,ngrid
386      kn(igrid,1)=knmin
387      km(igrid,1)=kmmin
388      kmpre(igrid,1)=km(igrid,1)
389                                                 ENDDO
390c
391c-----------------------------------------------------------------------
392      DO ilev=2,nlev-1
393                                                 DO igrid=1,ngrid
394c-----------------------------------------------------------------------
395c
396        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
397     &                                         *sn(igrid,ilev)
398        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
399     &                                         *sm(igrid,ilev)
400        kmpre(igrid,ilev)=km(igrid,ilev)
401c
402c-----------------------------------------------------------------------
403                                                 ENDDO
404      ENDDO
405c-----------------------------------------------------------------------
406c
407                                                 DO igrid=1,ngrid
408      kn(igrid,nlev)=kn(igrid,nlev-1)
409      km(igrid,nlev)=km(igrid,nlev-1)
410      kmpre(igrid,nlev)=km(igrid,nlev)
411                                                 ENDDO
412c
413c.......................................................................
414c  boucle sur les niveaux 2 a nlev-1
415c.......................................................................
416c
417c---->
418      DO 10001 ilev=2,nlev-1
419c---->
420      DO 10002 igrid=1,ngrid
421c
422c.......................................................................
423c
424c  calcul des termes sources et puits de l'equation de q2
425c  ------------------------------------------------------
426c
427        knq3=kn(igrid,ilev)*snq2(igrid,ilev)
428     &                                    /sn(igrid,ilev)
429        kmq3=km(igrid,ilev)*smq2(igrid,ilev)
430     &                                    /sm(igrid,ilev)
431c
432        termq=0.E+0
433        termq3=0.E+0
434        termqm2=0.E+0
435        termq3m2=0.E+0
436c
437        tmp1=dt*2.E+0 *km(igrid,ilev)*m2(igrid,ilev)
438        tmp2=dt*2.E+0 *kmq3*m2(igrid,ilev)
439        termqm2=termqm2
440     &    +dt*2.E+0 *km(igrid,ilev)*m2(igrid,ilev)
441     &    -dt*2.E+0 *kmq3*m2(igrid,ilev)
442        termq3m2=termq3m2
443     &    +dt*2.E+0 *kmq3*m2(igrid,ilev)
444c
445        termq=termq
446     &    -dt*2.E+0 *kn(igrid,ilev)*n2(igrid,ilev)
447     &    +dt*2.E+0 *knq3*n2(igrid,ilev)
448        termq3=termq3
449     &    -dt*2.E+0 *knq3*n2(igrid,ilev)
450c
451        termq3=termq3
452     &    -dt*2.E+0 *q(igrid,ilev)**3 / (b1*long(igrid,ilev))
453c
454c.......................................................................
455c
456c  resolution stationnaire couplee avec le gradient de vitesse local
457c  -----------------------------------------------------------------
458c
459c  -----{on cherche le cisaillement qui annule l'equation de q^2
460c        supposee en q3}
461c
462        tmp1=termq+termq3
463        tmp2=termqm2+termq3m2
464        m2cstat=m2(igrid,ilev)
465     &      -(tmp1+tmp2)/(dt*2.E+0*km(igrid,ilev))
466        mcstat=sqrt(m2cstat)
467c
468c  -----{puis on ecrit la valeur de q qui annule l'equation de m
469c        supposee en q3}
470c
471        IF (ilev.eq.2) THEN
472          kmcstat=1.E+0 / mcstat
473     &    *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
474     &                        *mpre(igrid,ilev+1)
475     &      +unsdz(igrid,ilev-1)
476     &              *cd(igrid)
477     &              *( sqrt(u(igrid,3)**2+v(igrid,3)**2)
478     &                -mcstat/unsdzdec(igrid,ilev)
479     &                -mpre(igrid,ilev+1)/unsdzdec(igrid,ilev+1) )**2)
480     &      /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
481        ELSE
482          kmcstat=1.E+0 / mcstat
483     &    *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
484     &                        *mpre(igrid,ilev+1)
485     &      +unsdz(igrid,ilev-1)*kmpre(igrid,ilev-1)
486     &                          *mpre(igrid,ilev-1) )
487     &      /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
488        ENDIF
489        tmp2=kmcstat
490     &      /( sm(igrid,ilev)/q2(igrid,ilev) )
491     &      /long(igrid,ilev)
492        qcstat=tmp2**(1.E+0/3.E+0)
493        q2cstat=qcstat**2
494c
495c.......................................................................
496c
497c  choix de la solution finale
498c  ---------------------------
499c
500          q(igrid,ilev)=qcstat
501          q2(igrid,ilev)=q2cstat
502          m(igrid,ilev)=mcstat
503          m2(igrid,ilev)=m2cstat
504c
505c --->
506c       pour des raisons simples q2 est minore
507c<---
508c
509        IF (q2(igrid,ilev).lt.q2min) THEN
510          q2(igrid,ilev)=q2min
511          q(igrid,ilev)=sqrt(q2min)
512        ENDIF
513c
514c.......................................................................
515c
516c  calcul final de kn et km
517c  ------------------------
518c
519        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
520     &                                           * n2(igrid,ilev)
521        IF (gn.lt.gnmin) gn=gnmin
522        IF (gn.gt.gnmax) gn=gnmax
523        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
524        sm(igrid,ilev)=
525     &    (cm1+cm2*gn)
526     &   /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )
527        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
528     &                 *sn(igrid,ilev)
529        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
530     &                 *sm(igrid,ilev)
531c
532c.......................................................................
533c
53410002 CONTINUE
535c
53610001 CONTINUE
537c
538c.......................................................................
539c
540c
541                                                 DO igrid=1,ngrid
542      kn(igrid,1)=knmin
543      km(igrid,1)=kmmin
544      q2(igrid,nlev)=q2(igrid,nlev-1)
545      q(igrid,nlev)=q(igrid,nlev-1)
546      kn(igrid,nlev)=kn(igrid,nlev-1)
547      km(igrid,nlev)=km(igrid,nlev-1)
548                                                 ENDDO
549
550c
551c.......................................................................
552c
553      RETURN
554      END
Note: See TracBrowser for help on using the repository browser.