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

Contents of /trunk/libf/dyn3d/vlspltqs.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 10 - (show annotations)
Fri Apr 18 14:45:53 2008 UTC (16 years, 1 month ago) by guez
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 !
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 REAL, intent(in):: p(ip1jmp1,llmp1)
40 real teta(ip1jmp1,llm),pk(ip1jmp1,llm)
41 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