/[lmdze]/trunk/dyn3d/vlspltqs.f
ViewVC logotype

Annotation of /trunk/dyn3d/vlspltqs.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 10 - (hide annotations)
Fri Apr 18 14:45:53 2008 UTC (16 years, 1 month ago) by guez
Original Path: trunk/libf/dyn3d/vlspltqs.f
File size: 20979 byte(s) 
Added NetCDF directory "/home/guez/include" in "g95.mk" and
"nag_tools.mk".

Added some "intent" attributes in "PVtheta", "advtrac", "caladvtrac",
"calfis", "diagedyn", "dissip", "vlspltqs", "aeropt", "ajsec",
"calltherm", "clmain", "cltrac", "cltracrn", "concvl", "conema3",
"conflx", "fisrtilp", "newmicro", "nuage", "diagcld1", "diagcld2",
"drag_noro", "lift_noro", "SUGWD", "physiq", "phytrac", "radlwsw", "thermcell".

Removed the case "ierr == 0" in "abort_gcm"; moved call to "histclo"
and messages for end of run from "abort_gcm" to "gcm"; replaced call
to "abort_gcm" in "leapfrog" by exit from outer loop.

In "calfis": removed argument "pp" and variable "unskap"; changed
"pksurcp" from scalar to rank 2; use "pressure_var"; rewrote
computation of "zplev", "zplay", "ztfi", "pcvgt" using "dyn_phy";
added computation of "pls".

Removed unused variable in "dynredem0".

In "exner_hyb": changed "dellta" from scalar to rank 1; replaced call
to "ssum" by call to "sum"; removed variables "xpn" and "xps";
replaced some loops by array expressions.

In "leapfrog": use "pressure_var"; deleted variables "p", "longcles".

Converted common blocks "YOECUMF", "YOEGWD" to modules.

Removed argument "pplay" in "cvltr", "diagetpq", "nflxtr".

Created module "raddimlw" from include file "raddimlw.h".

Corrected call to "new_unit" in "test_disvert".
1 guez 3 !
2     ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/vlspltqs.F,v 1.2 2005/02/24 12:16:57 fairhead Exp $
3     !
4     SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
5     , p,pk,teta )
6     c
7     c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron
8     c
9     c ********************************************************************
10     c Shema d'advection " pseudo amont " .
11     c + test sur humidite specifique: Q advecte< Qsat aval
12     c (F. Codron, 10/99)
13     c ********************************************************************
14     c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
15     c
16     c pente_max facteur de limitation des pentes: 2 en general
17     c 0 pour un schema amont
18     c pbaru,pbarv,w flux de masse en u ,v ,w
19     c pdt pas de temps
20     c
21     c teta temperature potentielle, p pression aux interfaces,
22     c pk exner au milieu des couches necessaire pour calculer Qsat
23     c --------------------------------------------------------------------
24     use dimens_m
25     use paramet_m
26     use comconst
27     use comvert
28     use logic
29     IMPLICIT NONE
30     c
31    
32     c
33     c Arguments:
34     c ----------
35     REAL masse(ip1jmp1,llm),pente_max
36     REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
37     REAL q(ip1jmp1,llm)
38     REAL w(ip1jmp1,llm),pdt
39 guez 10 REAL, intent(in):: p(ip1jmp1,llmp1)
40     real teta(ip1jmp1,llm),pk(ip1jmp1,llm)
41 guez 3 c
42     c Local
43     c ---------
44     c
45     INTEGER i,ij,l,j,ii
46     c
47     REAL qsat(ip1jmp1,llm)
48     REAL zm(ip1jmp1,llm)
49     REAL mu(ip1jmp1,llm)
50     REAL mv(ip1jm,llm)
51     REAL mw(ip1jmp1,llm+1)
52     REAL zq(ip1jmp1,llm)
53     REAL temps1,temps2,temps3
54     REAL zzpbar, zzw
55     LOGICAL testcpu
56     SAVE testcpu
57     SAVE temps1,temps2,temps3
58    
59     REAL qmin,qmax
60     DATA qmin,qmax/0.,1.e33/
61     DATA testcpu/.false./
62     DATA temps1,temps2,temps3/0.,0.,0./
63    
64     c--pour rapport de melange saturant--
65    
66     REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
67     REAL ptarg,pdelarg,foeew,zdelta
68     REAL tempe(ip1jmp1)
69    
70     c fonction psat(T)
71    
72     FOEEW ( PTARG,PDELARG ) = EXP (
73     * (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)
74     * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
75    
76     r2es = 380.11733
77     r3les = 17.269
78     r3ies = 21.875
79     r4les = 35.86
80     r4ies = 7.66
81     retv = 0.6077667
82     rtt = 273.16
83    
84     c-- Calcul de Qsat en chaque point
85     c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
86     c pour eviter une exponentielle.
87     DO l = 1, llm
88     DO ij = 1, ip1jmp1
89     tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
90     ENDDO
91     DO ij = 1, ip1jmp1
92     zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
93     play = 0.5*(p(ij,l)+p(ij,l+1))
94     qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
95     qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
96     ENDDO
97     ENDDO
98    
99     c PRINT*,'Debut vlsplt version debug sans vlyqs'
100    
101     zzpbar = 0.5 * pdt
102     zzw = pdt
103     DO l=1,llm
104     DO ij = iip2,ip1jm
105     mu(ij,l)=pbaru(ij,l) * zzpbar
106     ENDDO
107     DO ij=1,ip1jm
108     mv(ij,l)=pbarv(ij,l) * zzpbar
109     ENDDO
110     DO ij=1,ip1jmp1
111     mw(ij,l)=w(ij,l) * zzw
112     ENDDO
113     ENDDO
114    
115     DO ij=1,ip1jmp1
116     mw(ij,llm+1)=0.
117     ENDDO
118    
119     CALL SCOPY(ijp1llm,q,1,zq,1)
120     CALL SCOPY(ijp1llm,masse,1,zm,1)
121    
122     c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
123     call vlxqs(zq,pente_max,zm,mu,qsat)
124    
125    
126     c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
127    
128     call vlyqs(zq,pente_max,zm,mv,qsat)
129    
130    
131     c call minmaxq(zq,qmin,qmax,'avant vlz ')
132    
133     call vlz(zq,pente_max,zm,mw)
134    
135    
136     c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
137     c call minmaxq(zm,qmin,qmax,'M avant vlyqs ')
138    
139     call vlyqs(zq,pente_max,zm,mv,qsat)
140    
141    
142     c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
143     c call minmaxq(zm,qmin,qmax,'M avant vlxqs ')
144    
145     call vlxqs(zq,pente_max,zm,mu,qsat)
146    
147     c call minmaxq(zq,qmin,qmax,'apres vlxqs ')
148     c call minmaxq(zm,qmin,qmax,'M apres vlxqs ')
149    
150    
151     DO l=1,llm
152     DO ij=1,ip1jmp1
153     q(ij,l)=zq(ij,l)
154     ENDDO
155     DO ij=1,ip1jm+1,iip1
156     q(ij+iim,l)=q(ij,l)
157     ENDDO
158     ENDDO
159    
160     RETURN
161     END
162     SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
163     c
164     c Auteurs: P.Le Van, F.Hourdin, F.Forget
165     c
166     c ********************************************************************
167     c Shema d'advection " pseudo amont " .
168     c ********************************************************************
169     c
170     c --------------------------------------------------------------------
171     use dimens_m
172     use paramet_m
173     use comconst
174     use comvert
175     use logic
176     IMPLICIT NONE
177     c
178     c
179     c
180     c Arguments:
181     c ----------
182     REAL masse(ip1jmp1,llm),pente_max
183     REAL u_m( ip1jmp1,llm )
184     REAL q(ip1jmp1,llm)
185     REAL qsat(ip1jmp1,llm)
186     c
187     c Local
188     c ---------
189     c
190     INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
191     INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
192     c
193     REAL new_m,zu_m,zdum(ip1jmp1,llm)
194     REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
195     REAL zz(ip1jmp1)
196     REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
197     REAL u_mq(ip1jmp1,llm)
198    
199     Logical first,testcpu
200     SAVE first,testcpu
201    
202     REAL SSUM
203     REAL temps0,temps1,temps2,temps3,temps4,temps5
204     SAVE temps0,temps1,temps2,temps3,temps4,temps5
205    
206    
207     DATA first,testcpu/.true.,.false./
208    
209     IF(first) THEN
210     temps1=0.
211     temps2=0.
212     temps3=0.
213     temps4=0.
214     temps5=0.
215     first=.false.
216     ENDIF
217    
218     c calcul de la pente a droite et a gauche de la maille
219    
220    
221     IF (pente_max.gt.-1.e-5) THEN
222     c IF (pente_max.gt.10) THEN
223    
224     c calcul des pentes avec limitation, Van Leer scheme I:
225     c -----------------------------------------------------
226    
227     c calcul de la pente aux points u
228     DO l = 1, llm
229     DO ij=iip2,ip1jm-1
230     dxqu(ij)=q(ij+1,l)-q(ij,l)
231     c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
232     c sigu(ij)=u_m(ij,l)/masse(ij,l)
233     ENDDO
234     DO ij=iip1+iip1,ip1jm,iip1
235     dxqu(ij)=dxqu(ij-iim)
236     c sigu(ij)=sigu(ij-iim)
237     ENDDO
238    
239     DO ij=iip2,ip1jm
240     adxqu(ij)=abs(dxqu(ij))
241     ENDDO
242    
243     c calcul de la pente maximum dans la maille en valeur absolue
244    
245     DO ij=iip2+1,ip1jm
246     dxqmax(ij,l)=pente_max*
247     , min(adxqu(ij-1),adxqu(ij))
248     c limitation subtile
249     c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
250    
251    
252     ENDDO
253    
254     DO ij=iip1+iip1,ip1jm,iip1
255     dxqmax(ij-iim,l)=dxqmax(ij,l)
256     ENDDO
257    
258     DO ij=iip2+1,ip1jm
259     IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
260     dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
261     ELSE
262     c extremum local
263     dxq(ij,l)=0.
264     ENDIF
265     dxq(ij,l)=0.5*dxq(ij,l)
266     dxq(ij,l)=
267     , sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
268     ENDDO
269    
270     ENDDO ! l=1,llm
271    
272     ELSE ! (pente_max.lt.-1.e-5)
273    
274     c Pentes produits:
275     c ----------------
276    
277     DO l = 1, llm
278     DO ij=iip2,ip1jm-1
279     dxqu(ij)=q(ij+1,l)-q(ij,l)
280     ENDDO
281     DO ij=iip1+iip1,ip1jm,iip1
282     dxqu(ij)=dxqu(ij-iim)
283     ENDDO
284    
285     DO ij=iip2+1,ip1jm
286     zz(ij)=dxqu(ij-1)*dxqu(ij)
287     zz(ij)=zz(ij)+zz(ij)
288     IF(zz(ij).gt.0) THEN
289     dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
290     ELSE
291     c extremum local
292     dxq(ij,l)=0.
293     ENDIF
294     ENDDO
295    
296     ENDDO
297    
298     ENDIF ! (pente_max.lt.-1.e-5)
299    
300     c bouclage de la pente en iip1:
301     c -----------------------------
302    
303     DO l=1,llm
304     DO ij=iip1+iip1,ip1jm,iip1
305     dxq(ij-iim,l)=dxq(ij,l)
306     ENDDO
307    
308     DO ij=1,ip1jmp1
309     iadvplus(ij,l)=0
310     ENDDO
311    
312     ENDDO
313    
314    
315     c calcul des flux a gauche et a droite
316    
317     c on cumule le flux correspondant a toutes les mailles dont la masse
318     c au travers de la paroi pENDant le pas de temps.
319     c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
320     DO l=1,llm
321     DO ij=iip2,ip1jm-1
322     IF (u_m(ij,l).gt.0.) THEN
323     zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
324     u_mq(ij,l)=u_m(ij,l)*
325     $ min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
326     ELSE
327     zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
328     u_mq(ij,l)=u_m(ij,l)*
329     $ min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
330     ENDIF
331     ENDDO
332     ENDDO
333    
334    
335     c detection des points ou on advecte plus que la masse de la
336     c maille
337     DO l=1,llm
338     DO ij=iip2,ip1jm-1
339     IF(zdum(ij,l).lt.0) THEN
340     iadvplus(ij,l)=1
341     u_mq(ij,l)=0.
342     ENDIF
343     ENDDO
344     ENDDO
345     DO l=1,llm
346     DO ij=iip1+iip1,ip1jm,iip1
347     iadvplus(ij,l)=iadvplus(ij-iim,l)
348     ENDDO
349     ENDDO
350    
351    
352    
353     c traitement special pour le cas ou on advecte en longitude plus que le
354     c contenu de la maille.
355     c cette partie est mal vectorisee.
356    
357     c pas d'influence de la pression saturante (pour l'instant)
358    
359     c calcul du nombre de maille sur lequel on advecte plus que la maille.
360    
361     n0=0
362     DO l=1,llm
363     nl(l)=0
364     DO ij=iip2,ip1jm
365     nl(l)=nl(l)+iadvplus(ij,l)
366     ENDDO
367     n0=n0+nl(l)
368     ENDDO
369    
370     IF(n0.gt.0) THEN
371     ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
372     ccc & ,'contenu de la maille : ',n0
373    
374     DO l=1,llm
375     IF(nl(l).gt.0) THEN
376     iju=0
377     c indicage des mailles concernees par le traitement special
378     DO ij=iip2,ip1jm
379     IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
380     iju=iju+1
381     indu(iju)=ij
382     ENDIF
383     ENDDO
384     niju=iju
385     c PRINT*,'niju,nl',niju,nl(l)
386    
387     c traitement des mailles
388     DO iju=1,niju
389     ij=indu(iju)
390     j=(ij-1)/iip1+1
391     zu_m=u_m(ij,l)
392     u_mq(ij,l)=0.
393     IF(zu_m.gt.0.) THEN
394     ijq=ij
395     i=ijq-(j-1)*iip1
396     c accumulation pour les mailles completements advectees
397     do while(zu_m.gt.masse(ijq,l))
398     u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
399     zu_m=zu_m-masse(ijq,l)
400     i=mod(i-2+iim,iim)+1
401     ijq=(j-1)*iip1+i
402     ENDDO
403     c ajout de la maille non completement advectee
404     u_mq(ij,l)=u_mq(ij,l)+zu_m*
405     & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
406     ELSE
407     ijq=ij+1
408     i=ijq-(j-1)*iip1
409     c accumulation pour les mailles completements advectees
410     do while(-zu_m.gt.masse(ijq,l))
411     u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
412     zu_m=zu_m+masse(ijq,l)
413     i=mod(i,iim)+1
414     ijq=(j-1)*iip1+i
415     ENDDO
416     c ajout de la maille non completement advectee
417     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
418     & 0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
419     ENDIF
420     ENDDO
421     ENDIF
422     ENDDO
423     ENDIF ! n0.gt.0
424    
425    
426    
427     c bouclage en latitude
428    
429     DO l=1,llm
430     DO ij=iip1+iip1,ip1jm,iip1
431     u_mq(ij,l)=u_mq(ij-iim,l)
432     ENDDO
433     ENDDO
434    
435    
436     c calcul des tendances
437    
438     DO l=1,llm
439     DO ij=iip2+1,ip1jm
440     new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
441     q(ij,l)=(q(ij,l)*masse(ij,l)+
442     & u_mq(ij-1,l)-u_mq(ij,l))
443     & /new_m
444     masse(ij,l)=new_m
445     ENDDO
446     c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
447     DO ij=iip1+iip1,ip1jm,iip1
448     q(ij-iim,l)=q(ij,l)
449     masse(ij-iim,l)=masse(ij,l)
450     ENDDO
451     ENDDO
452    
453     c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
454     c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
455    
456    
457     RETURN
458     END
459     SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
460     c
461     c Auteurs: P.Le Van, F.Hourdin, F.Forget
462     c
463     c ********************************************************************
464     c Shema d'advection " pseudo amont " .
465     c ********************************************************************
466     c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
467     c qsat est un argument de sortie pour le s-pg ....
468     c
469     c
470     c --------------------------------------------------------------------
471     c
472     use dimens_m
473     use paramet_m
474     use comconst
475     use comvert
476     use logic
477     use comgeom
478     IMPLICIT NONE
479     c
480     c
481     c Arguments:
482     c ----------
483     REAL masse(ip1jmp1,llm),pente_max
484     REAL masse_adv_v( ip1jm,llm)
485     REAL q(ip1jmp1,llm)
486     REAL qsat(ip1jmp1,llm)
487     c
488     c Local
489     c ---------
490     c
491     INTEGER i,ij,l
492     c
493     REAL airej2,airejjm,airescb(iim),airesch(iim)
494     REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
495     REAL adyqv(ip1jm),dyqmax(ip1jmp1)
496     REAL qbyv(ip1jm,llm)
497    
498     REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
499     c REAL newq,oldmasse
500     Logical first,testcpu
501     REAL temps0,temps1,temps2,temps3,temps4,temps5
502     SAVE temps0,temps1,temps2,temps3,temps4,temps5
503     SAVE first,testcpu
504    
505     REAL convpn,convps,convmpn,convmps
506     REAL sinlon(iip1),sinlondlon(iip1)
507     REAL coslon(iip1),coslondlon(iip1)
508     SAVE sinlon,coslon,sinlondlon,coslondlon
509     SAVE airej2,airejjm
510     c
511     c
512     REAL SSUM
513    
514     DATA first,testcpu/.true.,.false./
515     DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
516    
517     IF(first) THEN
518     PRINT*,'Shema Amont nouveau appele dans Vanleer '
519     first=.false.
520     do i=2,iip1
521     coslon(i)=cos(rlonv(i))
522     sinlon(i)=sin(rlonv(i))
523     coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
524     sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
525     ENDDO
526     coslon(1)=coslon(iip1)
527     coslondlon(1)=coslondlon(iip1)
528     sinlon(1)=sinlon(iip1)
529     sinlondlon(1)=sinlondlon(iip1)
530     airej2 = SSUM( iim, aire(iip2), 1 )
531     airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
532     ENDIF
533    
534     c
535    
536    
537     DO l = 1, llm
538     c
539     c --------------------------------
540     c CALCUL EN LATITUDE
541     c --------------------------------
542    
543     c On commence par calculer la valeur du traceur moyenne sur le premier cercle
544     c de latitude autour du pole (qpns pour le pole nord et qpsn pour
545     c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
546    
547     DO i = 1, iim
548     airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
549     airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
550     ENDDO
551     qpns = SSUM( iim, airescb ,1 ) / airej2
552     qpsn = SSUM( iim, airesch ,1 ) / airejjm
553    
554     c calcul des pentes aux points v
555    
556     DO ij=1,ip1jm
557     dyqv(ij)=q(ij,l)-q(ij+iip1,l)
558     adyqv(ij)=abs(dyqv(ij))
559     ENDDO
560    
561     c calcul des pentes aux points scalaires
562    
563     DO ij=iip2,ip1jm
564     dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
565     dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
566     dyqmax(ij)=pente_max*dyqmax(ij)
567     ENDDO
568    
569     c calcul des pentes aux poles
570    
571     DO ij=1,iip1
572     dyq(ij,l)=qpns-q(ij+iip1,l)
573     dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
574     ENDDO
575    
576     c filtrage de la derivee
577     dyn1=0.
578     dys1=0.
579     dyn2=0.
580     dys2=0.
581     DO ij=1,iim
582     dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
583     dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
584     dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
585     dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
586     ENDDO
587     DO ij=1,iip1
588     dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
589     dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
590     ENDDO
591    
592     c calcul des pentes limites aux poles
593    
594     fn=1.
595     fs=1.
596     DO ij=1,iim
597     IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
598     fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
599     ENDIF
600     IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
601     fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
602     ENDIF
603     ENDDO
604     DO ij=1,iip1
605     dyq(ij,l)=fn*dyq(ij,l)
606     dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
607     ENDDO
608    
609     CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
610     C En memoire de dIFferents tests sur la
611     C limitation des pentes aux poles.
612     CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
613     C PRINT*,dyq(1)
614     C PRINT*,dyqv(iip1+1)
615     C apn=abs(dyq(1)/dyqv(iip1+1))
616     C PRINT*,dyq(ip1jm+1)
617     C PRINT*,dyqv(ip1jm-iip1+1)
618     C aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
619     C DO ij=2,iim
620     C apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
621     C aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
622     C ENDDO
623     C apn=min(pente_max/apn,1.)
624     C aps=min(pente_max/aps,1.)
625     C
626     C
627     C cas ou on a un extremum au pole
628     C
629     C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
630     C & apn=0.
631     C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
632     C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
633     C & aps=0.
634     C
635     C limitation des pentes aux poles
636     C DO ij=1,iip1
637     C dyq(ij)=apn*dyq(ij)
638     C dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
639     C ENDDO
640     C
641     C test
642     C DO ij=1,iip1
643     C dyq(iip1+ij)=0.
644     C dyq(ip1jm+ij-iip1)=0.
645     C ENDDO
646     C DO ij=1,ip1jmp1
647     C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
648     C ENDDO
649     C
650     C changement 10 07 96
651     C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
652     C & THEN
653     C DO ij=1,iip1
654     C dyqmax(ij)=0.
655     C ENDDO
656     C ELSE
657     C DO ij=1,iip1
658     C dyqmax(ij)=pente_max*abs(dyqv(ij))
659     C ENDDO
660     C ENDIF
661     C
662     C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
663     C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
664     C &THEN
665     C DO ij=ip1jm+1,ip1jmp1
666     C dyqmax(ij)=0.
667     C ENDDO
668     C ELSE
669     C DO ij=ip1jm+1,ip1jmp1
670     C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
671     C ENDDO
672     C ENDIF
673     C fin changement 10 07 96
674     CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
675    
676     c calcul des pentes limitees
677    
678     DO ij=iip2,ip1jm
679     IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
680     dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
681     ELSE
682     dyq(ij,l)=0.
683     ENDIF
684     ENDDO
685    
686     ENDDO
687    
688     DO l=1,llm
689     DO ij=1,ip1jm
690     IF( masse_adv_v(ij,l).GT.0. ) THEN
691     qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) +
692     , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
693     ELSE
694     qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
695     , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
696     ENDIF
697     qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
698     ENDDO
699     ENDDO
700    
701    
702     DO l=1,llm
703     DO ij=iip2,ip1jm
704     newmasse=masse(ij,l)
705     & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
706     q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
707     & /newmasse
708     masse(ij,l)=newmasse
709     ENDDO
710     c.-. ancienne version
711     convpn=SSUM(iim,qbyv(1,l),1)/apoln
712     convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
713     DO ij = 1,iip1
714     newmasse=masse(ij,l)+convmpn*aire(ij)
715     q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
716     & newmasse
717     masse(ij,l)=newmasse
718     ENDDO
719     convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
720     convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
721     DO ij = ip1jm+1,ip1jmp1
722     newmasse=masse(ij,l)+convmps*aire(ij)
723     q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
724     & newmasse
725     masse(ij,l)=newmasse
726     ENDDO
727     c.-. fin ancienne version
728    
729     c._. nouvelle version
730     c convpn=SSUM(iim,qbyv(1,l),1)
731     c convmpn=ssum(iim,masse_adv_v(1,l),1)
732     c oldmasse=ssum(iim,masse(1,l),1)
733     c newmasse=oldmasse+convmpn
734     c newq=(q(1,l)*oldmasse+convpn)/newmasse
735     c newmasse=newmasse/apoln
736     c DO ij = 1,iip1
737     c q(ij,l)=newq
738     c masse(ij,l)=newmasse*aire(ij)
739     c ENDDO
740     c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
741     c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
742     c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
743     c newmasse=oldmasse+convmps
744     c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
745     c newmasse=newmasse/apols
746     c DO ij = ip1jm+1,ip1jmp1
747     c q(ij,l)=newq
748     c masse(ij,l)=newmasse*aire(ij)
749     c ENDDO
750     c._. fin nouvelle version
751     ENDDO
752    
753     RETURN
754     END
  ViewVC Help
Powered by ViewVC 1.1.21