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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 12 - (show annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 9 months ago) by guez
File size: 26057 byte(s) 
-- Minor modification of input/output:

Created procedure "read_logic". Variables of module "logic" are read
by "read_logic" instead of "conf_gcm". Variable "offline" of module
"conf_gcm" is read from namelist instead of "*.def".

Deleted arguments "dtime", "co2_ppm_etat0", "solaire_etat0",
"tabcntr0" and local variables "radpas", "tab_cntrl" of
"phyetat0". "phyetat0" does not read "controle" in "startphy.nc" any
longer. "phyetat0" now reads global attribute "itau_phy" from
"startphy.nc". "phyredem" does not create variable "controle" in
"startphy.nc" any longer. "phyredem" now writes global attribute
"itau_phy" of "startphy.nc". Deleted argument "tabcntr0" of
"printflag". Removed diagnostic messages written by "printflag" for
comparison of the variable "controle" of "startphy.nc" and the
variables read from "*.def" or namelist input.

-- Removing unwanted functionality:

Removed variable "lunout" from module "iniprint", replaced everywhere
by standard output.

Removed case "ocean == 'couple'" in "clmain", "interfsurf_hq" and
"physiq". Removed procedure "interfoce_cpl".

-- Should not change anything at run time:

Automated creation of graphs in documentation. More documentation on
input files.

Converted Fortran files to free format: "phyredem.f90", "printflag.f90".

Split module "clesphy" into "clesphys" and "clesphys2".

Removed variables "conser", "leapf", "forward", "apphys", "apdiss" and
"statcl" from module "logic". Added arguments "conser" to "advect",
"leapf" to "integrd". Added local variables "forward", "leapf",
"apphys", "conser", "apdiss" in "leapfrog".

Added intent attributes.

Deleted arguments "dtime" of "phyredem", "pdtime" of "flxdtdq", "sh"
of "phytrac", "dt" of "yamada".

Deleted local variables "dtime", "co2_ppm_etat0", "solaire_etat0",
"length", "tabcntr0" in "physiq". Replaced all references to "dtime"
by references to "pdtphys".
1 !
2 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advn.F,v 1.1.1.1 2004/05/19 12:53:06 lmdzadmin Exp $
3 !
4 SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)
5 c
6 c Auteur : F. Hourdin
7 c
8 c ********************************************************************
9 c Shema d'advection " pseudo amont " .
10 c ********************************************************************
11 c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
12 c
13 c pbaru,pbarv,w flux de masse en u ,v ,w
14 c pdt pas de temps
15 c
16 c --------------------------------------------------------------------
17 use dimens_m
18 use paramet_m
19 use comconst
20 use comvert
21 use logic
22 use comgeom
23 use iniprint
24 IMPLICIT NONE
25 c
26
27 c
28 c Arguments:
29 c ----------
30 integer mode
31 real masse(ip1jmp1,llm)
32 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
33 REAL q(ip1jmp1,llm)
34 REAL w(ip1jmp1,llm),pdt
35 c
36 c Local
37 c ---------
38 c
39 INTEGER i,ij,l,j,ii
40 integer ijlqmin,iqmin,jqmin,lqmin
41 integer ismin
42 c
43 real zm(ip1jmp1,llm),newmasse
44 real mu(ip1jmp1,llm)
45 real mv(ip1jm,llm)
46 real mw(ip1jmp1,llm+1)
47 real zq(ip1jmp1,llm),zz,qpn,qps
48 real zqg(ip1jmp1,llm),zqd(ip1jmp1,llm)
49 real zqs(ip1jmp1,llm),zqn(ip1jmp1,llm)
50 real zqh(ip1jmp1,llm),zqb(ip1jmp1,llm)
51 real temps0,temps1,temps2,temps3
52 real ztemps1,ztemps2,ztemps3,ssum
53 logical testcpu
54 save testcpu
55 save temps1,temps2,temps3
56 real zzpbar,zzw
57
58 real qmin,qmax
59 data qmin,qmax/0.,1./
60 data testcpu/.false./
61 data temps1,temps2,temps3/0.,0.,0./
62
63 zzpbar = 0.5 * pdt
64 zzw = pdt
65
66 DO l=1,llm
67 DO ij = iip2,ip1jm
68 mu(ij,l)=pbaru(ij,l) * zzpbar
69 ENDDO
70 DO ij=1,ip1jm
71 mv(ij,l)=pbarv(ij,l) * zzpbar
72 ENDDO
73 DO ij=1,ip1jmp1
74 mw(ij,l)=w(ij,l) * zzw
75 ENDDO
76 ENDDO
77
78 DO ij=1,ip1jmp1
79 mw(ij,llm+1)=0.
80 ENDDO
81
82 do l=1,llm
83 qpn=0.
84 qps=0.
85 do ij=1,iim
86 qpn=qpn+q(ij,l)*masse(ij,l)
87 qps=qps+q(ip1jm+ij,l)*masse(ip1jm+ij,l)
88 enddo
89 qpn=qpn/ssum(iim,masse(1,l),1)
90 qps=qps/ssum(iim,masse(ip1jm+1,l),1)
91 do ij=1,iip1
92 q(ij,l)=qpn
93 q(ip1jm+ij,l)=qps
94 enddo
95 enddo
96
97 do ij=1,ip1jmp1
98 mw(ij,llm+1)=0.
99 enddo
100 do l=1,llm
101 do ij=1,ip1jmp1
102 zq(ij,l)=q(ij,l)
103 zm(ij,l)=masse(ij,l)
104 enddo
105 enddo
106
107 c call minmaxq(zq,qmin,qmax,'avant vlx ')
108 call advnqx(zq,zqg,zqd)
109 call advnx(zq,zqg,zqd,zm,mu,mode)
110 call advnqy(zq,zqs,zqn)
111 call advny(zq,zqs,zqn,zm,mv)
112 call advnqz(zq,zqh,zqb)
113 call advnz(zq,zqh,zqb,zm,mw)
114 c call vlz(zq,0.,zm,mw)
115 call advnqy(zq,zqs,zqn)
116 call advny(zq,zqs,zqn,zm,mv)
117 call advnqx(zq,zqg,zqd)
118 call advnx(zq,zqg,zqd,zm,mu,mode)
119 c call minmaxq(zq,qmin,qmax,'apres vlx ')
120
121 do l=1,llm
122 do ij=1,ip1jmp1
123 q(ij,l)=zq(ij,l)
124 enddo
125 do ij=1,ip1jm+1,iip1
126 q(ij+iim,l)=q(ij,l)
127 enddo
128 enddo
129
130 RETURN
131 END
132
133 SUBROUTINE advnqx(q,qg,qd)
134 c
135 c Auteurs: Calcul des valeurs de q aux point u.
136 c
137 c --------------------------------------------------------------------
138 use dimens_m
139 use paramet_m
140 use iniprint
141 IMPLICIT NONE
142 c
143 c
144 c
145 c Arguments:
146 c ----------
147 real q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
148 c
149 c Local
150 c ---------
151 c
152 INTEGER ij,l
153 c
154 real dxqu(ip1jmp1),zqu(ip1jmp1)
155 real zqmax(ip1jmp1),zqmin(ip1jmp1)
156 logical extremum(ip1jmp1)
157
158 integer mode
159 save mode
160 data mode/1/
161
162 c calcul des pentes en u:
163 c -----------------------
164 if (mode.eq.0) then
165 do l=1,llm
166 do ij=1,ip1jm
167 qd(ij,l)=q(ij,l)
168 qg(ij,l)=q(ij,l)
169 enddo
170 enddo
171 else
172 do l = 1, llm
173 do ij=iip2,ip1jm-1
174 dxqu(ij)=q(ij+1,l)-q(ij,l)
175 zqu(ij)=0.5*(q(ij+1,l)+q(ij,l))
176 enddo
177 do ij=iip1+iip1,ip1jm,iip1
178 dxqu(ij)=dxqu(ij-iim)
179 zqu(ij)=zqu(ij-iim)
180 enddo
181 do ij=iip2,ip1jm-1
182 zqu(ij)=zqu(ij)-dxqu(ij+1)/12.
183 enddo
184 do ij=iip1+iip1,ip1jm,iip1
185 zqu(ij)=zqu(ij-iim)
186 enddo
187 do ij=iip2+1,ip1jm
188 zqu(ij)=zqu(ij)+dxqu(ij-1)/12.
189 enddo
190 do ij=iip1+iip1,ip1jm,iip1
191 zqu(ij-iim)=zqu(ij)
192 enddo
193
194 c calcul des valeurs max et min acceptees aux interfaces
195
196 do ij=iip2,ip1jm-1
197 zqmax(ij)=max(q(ij+1,l),q(ij,l))
198 zqmin(ij)=min(q(ij+1,l),q(ij,l))
199 enddo
200 do ij=iip1+iip1,ip1jm,iip1
201 zqmax(ij)=zqmax(ij-iim)
202 zqmin(ij)=zqmin(ij-iim)
203 enddo
204 do ij=iip2+1,ip1jm
205 extremum(ij)=dxqu(ij)*dxqu(ij-1).le.0.
206 enddo
207 do ij=iip1+iip1,ip1jm,iip1
208 extremum(ij-iim)=extremum(ij)
209 enddo
210 do ij=iip2,ip1jm
211 zqu(ij)=min(max(zqmin(ij),zqu(ij)),zqmax(ij))
212 enddo
213 do ij=iip2+1,ip1jm
214 if(extremum(ij)) then
215 qg(ij,l)=q(ij,l)
216 qd(ij,l)=q(ij,l)
217 else
218 qd(ij,l)=zqu(ij)
219 qg(ij,l)=zqu(ij-1)
220 endif
221 enddo
222 do ij=iip1+iip1,ip1jm,iip1
223 qd(ij-iim,l)=qd(ij,l)
224 qg(ij-iim,l)=qg(ij,l)
225 enddo
226
227 goto 8888
228
229 do ij=iip2+1,ip1jm
230 if(extremum(ij).and..not.extremum(ij-1))
231 s qd(ij-1,l)=q(ij,l)
232 enddo
233
234 do ij=iip1+iip1,ip1jm,iip1
235 qd(ij-iim,l)=qd(ij,l)
236 enddo
237 do ij=iip2,ip1jm-1
238 if (extremum(ij).and..not.extremum(ij+1))
239 s qg(ij+1,l)=q(ij,l)
240 enddo
241
242 do ij=iip1+iip1,ip1jm,iip1
243 qg(ij,l)=qg(ij-iim,l)
244 enddo
245 8888 continue
246 enddo
247 endif
248 RETURN
249 END
250 SUBROUTINE advnqy(q,qs,qn)
251 c
252 c Auteurs: Calcul des valeurs de q aux point v.
253 c
254 c --------------------------------------------------------------------
255 use dimens_m
256 use paramet_m
257 use iniprint
258 IMPLICIT NONE
259 c
260 c
261 c
262 c Arguments:
263 c ----------
264 real q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
265 c
266 c Local
267 c ---------
268 c
269 INTEGER ij,l
270 c
271 real dyqv(ip1jm),zqv(ip1jm,llm)
272 real zqmax(ip1jm),zqmin(ip1jm)
273 logical extremum(ip1jmp1)
274
275 integer mode
276 save mode
277 data mode/1/
278
279 if (mode.eq.0) then
280 do l=1,llm
281 do ij=1,ip1jmp1
282 qn(ij,l)=q(ij,l)
283 qs(ij,l)=q(ij,l)
284 enddo
285 enddo
286 else
287
288 c calcul des pentes en u:
289 c -----------------------
290 do l = 1, llm
291 do ij=1,ip1jm
292 dyqv(ij)=q(ij,l)-q(ij+iip1,l)
293 enddo
294
295 do ij=iip2,ip1jm-iip1
296 zqv(ij,l)=0.5*(q(ij+iip1,l)+q(ij,l))
297 zqv(ij,l)=zqv(ij,l)+(dyqv(ij+iip1)-dyqv(ij-iip1))/12.
298 enddo
299
300 do ij=iip2,ip1jm
301 extremum(ij)=dyqv(ij)*dyqv(ij-iip1).le.0.
302 enddo
303
304 c Pas de pentes aux poles
305 do ij=1,iip1
306 zqv(ij,l)=q(ij,l)
307 zqv(ip1jm-iip1+ij,l)=q(ip1jm+ij,l)
308 extremum(ij)=.true.
309 extremum(ip1jmp1-iip1+ij)=.true.
310 enddo
311
312 c calcul des valeurs max et min acceptees aux interfaces
313 do ij=1,ip1jm
314 zqmax(ij)=max(q(ij+iip1,l),q(ij,l))
315 zqmin(ij)=min(q(ij+iip1,l),q(ij,l))
316 enddo
317
318 do ij=1,ip1jm
319 zqv(ij,l)=min(max(zqmin(ij),zqv(ij,l)),zqmax(ij))
320 enddo
321
322 do ij=iip2,ip1jm
323 if(extremum(ij)) then
324 qs(ij,l)=q(ij,l)
325 qn(ij,l)=q(ij,l)
326 c if (.not.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
327 c if (.not.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
328 else
329 qs(ij,l)=zqv(ij,l)
330 qn(ij,l)=zqv(ij-iip1,l)
331 endif
332 enddo
333
334 do ij=1,iip1
335 qs(ij,l)=q(ij,l)
336 qn(ij,l)=q(ij,l)
337 qs(ip1jm+ij,l)=q(ip1jm+ij,l)
338 qn(ip1jm+ij,l)=q(ip1jm+ij,l)
339 enddo
340
341 enddo
342 endif
343 RETURN
344 END
345
346 SUBROUTINE advnqz(q,qh,qb)
347 c
348 c Auteurs: Calcul des valeurs de q aux point v.
349 c
350 c --------------------------------------------------------------------
351 use dimens_m
352 use paramet_m
353 use iniprint
354 IMPLICIT NONE
355 c
356 c
357 c
358 c Arguments:
359 c ----------
360 real q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
361 c
362 c Local
363 c ---------
364 c
365 INTEGER ij,l
366 c
367 real dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)
368 real zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)
369 logical extremum(ip1jmp1,llm)
370
371 integer mode
372 save mode
373
374 data mode/1/
375
376 c calcul des pentes en u:
377 c -----------------------
378
379 if (mode.eq.0) then
380 do l=1,llm
381 do ij=1,ip1jmp1
382 qb(ij,l)=q(ij,l)
383 qh(ij,l)=q(ij,l)
384 enddo
385 enddo
386 else
387 do l = 2, llm
388 do ij=1,ip1jmp1
389 dzqw(ij,l)=q(ij,l-1)-q(ij,l)
390 zqw(ij,l)=0.5*(q(ij,l-1)+q(ij,l))
391 enddo
392 enddo
393 do ij=1,ip1jmp1
394 dzqw(ij,1)=0.
395 dzqw(ij,llm+1)=0.
396 enddo
397 do l=2,llm
398 do ij=1,ip1jmp1
399 zqw(ij,l)=zqw(ij,l)+(dzqw(ij,l+1)-dzqw(ij,l-1))/12.
400 enddo
401 enddo
402 do l=2,llm-1
403 do ij=1,ip1jmp1
404 extremum(ij,l)=dzqw(ij,l)*dzqw(ij,l+1).le.0.
405 enddo
406 enddo
407
408 c Pas de pentes en bas et en haut
409 do ij=1,ip1jmp1
410 zqw(ij,2)=q(ij,1)
411 zqw(ij,llm)=q(ij,llm)
412 extremum(ij,1)=.true.
413 extremum(ij,llm)=.true.
414 enddo
415
416 c calcul des valeurs max et min acceptees aux interfaces
417 do l=2,llm
418 do ij=1,ip1jmp1
419 zqmax(ij,l)=max(q(ij,l-1),q(ij,l))
420 zqmin(ij,l)=min(q(ij,l-1),q(ij,l))
421 enddo
422 enddo
423
424 do l=2,llm
425 do ij=1,ip1jmp1
426 zqw(ij,l)=min(max(zqmin(ij,l),zqw(ij,l)),zqmax(ij,l))
427 enddo
428 enddo
429
430 do l=2,llm-1
431 do ij=1,ip1jmp1
432 if(extremum(ij,l)) then
433 qh(ij,l)=q(ij,l)
434 qb(ij,l)=q(ij,l)
435 else
436 qh(ij,l)=zqw(ij,l+1)
437 qb(ij,l)=zqw(ij,l)
438 endif
439 enddo
440 enddo
441 c do l=2,llm-1
442 c do ij=1,ip1jmp1
443 c if(extremum(ij,l)) then
444 c if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
445 c if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
446 c endif
447 c enddo
448 c enddo
449
450 do ij=1,ip1jmp1
451 qb(ij,1)=q(ij,1)
452 qh(ij,1)=q(ij,1)
453 qb(ij,llm)=q(ij,llm)
454 qh(ij,llm)=q(ij,llm)
455 enddo
456
457 endif
458
459 RETURN
460 END
461
462 SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)
463 c
464 c Auteur : F. Hourdin
465 c
466 c ********************************************************************
467 c Shema d'advection " pseudo amont " .
468 c ********************************************************************
469 c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
470 c
471 c
472 c --------------------------------------------------------------------
473 use dimens_m
474 use paramet_m
475 use comconst
476 use comvert
477 use logic
478 use iniprint
479 IMPLICIT NONE
480 c
481 c
482 c
483 c Arguments:
484 c ----------
485 integer mode
486 real masse(ip1jmp1,llm)
487 real u_m( ip1jmp1,llm )
488 real q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)
489 c
490 c Local
491 c ---------
492 c
493 INTEGER i,j,ij,l,indu(ip1jmp1),niju,iju,ijq
494 integer n0,nl(llm)
495 c
496 real new_m,zu_m,zdq,zz
497 real zsigg(ip1jmp1,llm),zsigd(ip1jmp1,llm),zsig
498 real u_mq(ip1jmp1,llm)
499
500 real zm,zq,zsigm,zsigp,zqm,zqp,zu
501
502 logical ladvplus(ip1jmp1,llm)
503
504 real prec
505 save prec
506
507 data prec/1.e-15/
508
509 do l=1,llm
510 do ij=iip2,ip1jm
511 zdq=qd(ij,l)-qg(ij,l)
512 c if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) then
513 c print*,'probleme au point ij=',ij,' l=',l
514 c print*,qd(ij,l),q(ij,l),qg(ij,l)
515 c qd(ij,l)=q(ij,l)
516 c qg(ij,l)=q(ij,l)
517 c endif
518 if(abs(zdq).gt.prec) then
519 zsigd(ij,l)=(q(ij,l)-qg(ij,l))/zdq
520 zsigg(ij,l)=1.-zsigd(ij,l)
521 c if(.not.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .and.
522 c s zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) then
523 c print*,'probleme au point ij=',ij,' l=',l
524 c print*,'sigg=',zsigg(ij,l),' sigd=',zsigd(ij,l)
525 c print*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
526 c stop
527 c endif
528 else
529 zsigd(ij,l)=0.5
530 zsigg(ij,l)=0.5
531 qd(ij,l)=q(ij,l)
532 qg(ij,l)=q(ij,l)
533 endif
534 enddo
535 enddo
536
537 c calcul de la pente maximum dans la maille en valeur absolue
538
539 do l=1,llm
540 do ij=iip2,ip1jm-1
541 if (u_m(ij,l).ge.0.) then
542 zsigp=zsigd(ij,l)
543 zsigm=zsigg(ij,l)
544 zqp=qd(ij,l)
545 zqm=qg(ij,l)
546 zm=masse(ij,l)
547 zq=q(ij,l)
548 else
549 zsigm=zsigd(ij+1,l)
550 zsigp=zsigg(ij+1,l)
551 zqm=qd(ij+1,l)
552 zqp=qg(ij+1,l)
553 zm=masse(ij+1,l)
554 zq=q(ij+1,l)
555 endif
556 zu=abs(u_m(ij,l))
557 ladvplus(ij,l)=zu.gt.zm
558 zsig=zu/zm
559 if(zsig.eq.0.) zsigp=0.1
560 if (mode.eq.1) then
561 if (zsig.le.zsigp) then
562 u_mq(ij,l)=u_m(ij,l)*zqp
563 else if (mode.eq.1) then
564 u_mq(ij,l)=
565 s sign(zm,u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm)
566 endif
567 else
568 if (zsig.le.zsigp) then
569 u_mq(ij,l)=u_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
570 else
571 zz=0.5*(zsig-zsigp)/zsigm
572 u_mq(ij,l)=sign(zm,u_m(ij,l))*( 0.5*(zq+zqp)*zsigp
573 s +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
574 endif
575 endif
576 c if(zsig.lt.0.) then
577 c print*,'au point ij=',ij,' l=',l,' sig=',zsig
578 c stop
579 c endif
580 enddo
581 enddo
582
583 do l=1,llm
584 do ij=iip1+iip1,ip1jm,iip1
585 u_mq(ij,l)=u_mq(ij-iim,l)
586 ladvplus(ij,l)=ladvplus(ij-iim,l)
587 enddo
588 enddo
589
590 c=================================================================
591 C SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
592 c=================================================================
593 c tris des regions a traiter
594 n0=0
595 do l=1,llm
596 nl(l)=0
597 do ij=iip2,ip1jm
598 if(ladvplus(ij,l)) then
599 nl(l)=nl(l)+1
600 u_mq(ij,l)=0.
601 endif
602 enddo
603 n0=n0+nl(l)
604 enddo
605
606 if(n0.gt.1) then
607 IF (prt_level > 9) print *,
608 & 'Nombre de points pour lesquels on advect plus que le'
609 & ,'contenu de la maille : ',n0
610
611 do l=1,llm
612 if(nl(l).gt.0) then
613 iju=0
614 c indicage des mailles concernees par le traitement special
615 do ij=iip2,ip1jm
616 if(ladvplus(ij,l).and.mod(ij,iip1).ne.0) then
617 iju=iju+1
618 indu(iju)=ij
619 endif
620 enddo
621 niju=iju
622 c print*,'niju,nl',niju,nl(l)
623
624 c traitement des mailles
625 do iju=1,niju
626 ij=indu(iju)
627 j=(ij-1)/iip1+1
628 zu_m=u_m(ij,l)
629 u_mq(ij,l)=0.
630 if(zu_m.gt.0.) then
631 ijq=ij
632 i=ijq-(j-1)*iip1
633 c accumulation pour les mailles completements advectees
634 do while(zu_m.gt.masse(ijq,l))
635 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
636 zu_m=zu_m-masse(ijq,l)
637 i=mod(i-2+iim,iim)+1
638 ijq=(j-1)*iip1+i
639 enddo
640 c MODIFS SPECIFIQUES DU SCHEMA
641 c ajout de la maille non completement advectee
642 zsig=zu_m/masse(ijq,l)
643 if(zsig.le.zsigd(ijq,l)) then
644 u_mq(ij,l)=u_mq(ij,l)+zu_m*(qd(ijq,l)
645 s -0.5*zsig/zsigd(ijq,l)*(qd(ijq,l)-q(ijq,l)))
646 else
647 c u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
648 c goto 8888
649 zz=0.5*(zsig-zsigd(ijq,l))/zsigg(ijq,l)
650 if(.not.(zz.gt.0..and.zz.le.0.5)) then
651 print *,'probleme2 au point ij=',ij,
652 s ' l=',l
653 print *,'zz=',zz
654 stop
655 endif
656 u_mq(ij,l)=u_mq(ij,l)+masse(ijq,l)*(
657 s 0.5*(q(ijq,l)+qd(ijq,l))*zsigd(ijq,l)
658 s +(zsig-zsigd(ijq,l))*(q(ijq,l)+zz*(qg(ijq,l)-q(ijq,l))) )
659 endif
660 else
661 ijq=ij+1
662 i=ijq-(j-1)*iip1
663 c accumulation pour les mailles completements advectees
664 do while(-zu_m.gt.masse(ijq,l))
665 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
666 zu_m=zu_m+masse(ijq,l)
667 i=mod(i,iim)+1
668 ijq=(j-1)*iip1+i
669 enddo
670 c ajout de la maille non completement advectee
671 c 2eme MODIF SPECIFIQUE
672 zsig=-zu_m/masse(ij+1,l)
673 if(zsig.le.zsigg(ijq,l)) then
674 u_mq(ij,l)=u_mq(ij,l)+zu_m*(qg(ijq,l)
675 s -0.5*zsig/zsigg(ijq,l)*(qg(ijq,l)-q(ijq,l)))
676 else
677 c u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
678 c goto 9999
679 zz=0.5*(zsig-zsigg(ijq,l))/zsigd(ijq,l)
680 if(.not.(zz.gt.0..and.zz.le.0.5)) then
681 print *,'probleme22 au point ij=',ij
682 s ,' l=',l
683 print *,'zz=',zz
684 stop
685 endif
686 u_mq(ij,l)=u_mq(ij,l)-masse(ijq,l)*(
687 s 0.5*(q(ijq,l)+qg(ijq,l))*zsigg(ijq,l)
688 s +(zsig-zsigg(ijq,l))*
689 s (q(ijq,l)+zz*(qd(ijq,l)-q(ijq,l))) )
690 endif
691 c fin de la modif
692 endif
693 enddo
694 endif
695 enddo
696 endif ! n0.gt.0
697
698 c bouclage en latitude
699 do l=1,llm
700 do ij=iip1+iip1,ip1jm,iip1
701 u_mq(ij,l)=u_mq(ij-iim,l)
702 enddo
703 enddo
704
705 c=================================================================
706 c CALCUL DE LA CONVERGENCE DES FLUX
707 c=================================================================
708
709 do l=1,llm
710 do ij=iip2+1,ip1jm
711 new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
712 q(ij,l)=(q(ij,l)*masse(ij,l)+
713 & u_mq(ij-1,l)-u_mq(ij,l))
714 & /new_m
715 masse(ij,l)=new_m
716 enddo
717 c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
718 do ij=iip1+iip1,ip1jm,iip1
719 q(ij-iim,l)=q(ij,l)
720 masse(ij-iim,l)=masse(ij,l)
721 enddo
722 enddo
723
724 RETURN
725 END
726 SUBROUTINE advny(q,qs,qn,masse,v_m)
727 c
728 c Auteur : F. Hourdin
729 c
730 c ********************************************************************
731 c Shema d'advection " pseudo amont " .
732 c ********************************************************************
733 c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
734 c
735 c
736 c --------------------------------------------------------------------
737 use dimens_m
738 use paramet_m
739 use comgeom
740 use iniprint
741 IMPLICIT NONE
742 c
743 c
744 c
745 c Arguments:
746 c ----------
747 real masse(ip1jmp1,llm)
748 real v_m( ip1jm,llm )
749 real q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)
750 c
751 c Local
752 c ---------
753 c
754 INTEGER ij,l
755 c
756 real new_m,zdq,zz
757 real zsigs(ip1jmp1),zsign(ip1jmp1),zsig
758 real v_mq(ip1jm,llm)
759 real convpn,convps,convmpn,convmps,massen,masses
760 real zm,zq,zsigm,zsigp,zqm,zqp
761 real ssum
762 real prec
763 save prec
764
765 data prec/1.e-15/
766 do l=1,llm
767 do ij=1,ip1jmp1
768 zdq=qn(ij,l)-qs(ij,l)
769 c if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) then
770 c print*,'probleme au point ij=',ij,' l=',l,' advnqx'
771 c print*,qn(ij,l),q(ij,l),qs(ij,l)
772 c qn(ij,l)=q(ij,l)
773 c qs(ij,l)=q(ij,l)
774 c endif
775 if(abs(zdq).gt.prec) then
776 zsign(ij)=(q(ij,l)-qs(ij,l))/zdq
777 zsigs(ij)=1.-zsign(ij)
778 c if(.not.(zsign(ij).ge.0..and.zsign(ij).le.1. .and.
779 c s zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) then
780 c print*,'probleme au point ij=',ij,' l=',l
781 c print*,'sigs=',zsigs(ij),' sign=',zsign(ij)
782 c stop
783 c endif
784 else
785 zsign(ij)=0.5
786 zsigs(ij)=0.5
787 endif
788 enddo
789
790 c calcul de la pente maximum dans la maille en valeur absolue
791
792 do ij=1,ip1jm
793 if (v_m(ij,l).ge.0.) then
794 zsigp=zsign(ij+iip1)
795 zsigm=zsigs(ij+iip1)
796 zqp=qn(ij+iip1,l)
797 zqm=qs(ij+iip1,l)
798 zm=masse(ij+iip1,l)
799 zq=q(ij+iip1,l)
800 else
801 zsigm=zsign(ij)
802 zsigp=zsigs(ij)
803 zqm=qn(ij,l)
804 zqp=qs(ij,l)
805 zm=masse(ij,l)
806 zq=q(ij,l)
807 endif
808 zsig=abs(v_m(ij,l))/zm
809 if(zsig.eq.0.) zsigp=0.1
810 if (zsig.le.zsigp) then
811 v_mq(ij,l)=v_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
812 else
813 zz=0.5*(zsig-zsigp)/zsigm
814 v_mq(ij,l)=sign(zm,v_m(ij,l))*( 0.5*(zq+zqp)*zsigp
815 s +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
816 endif
817 enddo
818 enddo
819
820 do l=1,llm
821 do ij=iip2,ip1jm
822 new_m=masse(ij,l)
823 & +v_m(ij,l)-v_m(ij-iip1,l)
824 q(ij,l)=(q(ij,l)*masse(ij,l)+v_mq(ij,l)-v_mq(ij-iip1,l))
825 & /new_m
826 masse(ij,l)=new_m
827 enddo
828 c.-. ancienne version
829 convpn=SSUM(iim,v_mq(1,l),1)
830 convmpn=ssum(iim,v_m(1,l),1)
831 massen=ssum(iim,masse(1,l),1)
832 new_m=massen+convmpn
833 q(1,l)=(q(1,l)*massen+convpn)/new_m
834 do ij = 1,iip1
835 q(ij,l)=q(1,l)
836 masse(ij,l)=new_m*aire(ij)/apoln
837 enddo
838
839 convps=-SSUM(iim,v_mq(ip1jm-iim,l),1)
840 convmps=-ssum(iim,v_m(ip1jm-iim,l),1)
841 masses=ssum(iim,masse(ip1jm+1,l),1)
842 new_m=masses+convmps
843 q(ip1jm+1,l)=(q(ip1jm+1,l)*masses+convps)/new_m
844 do ij = ip1jm+1,ip1jmp1
845 q(ij,l)=q(ip1jm+1,l)
846 masse(ij,l)=new_m*aire(ij)/apols
847 enddo
848 enddo
849
850 RETURN
851 END
852 SUBROUTINE advnz(q,qh,qb,masse,w_m)
853 c
854 c Auteurs: F.Hourdin
855 c
856 c ********************************************************************
857 c Shema d'advection " pseudo amont " .
858 c b designe le bas et h le haut
859 c il y a une correspondance entre le b en z et le d en x
860 c ********************************************************************
861 c
862 c
863 c --------------------------------------------------------------------
864 use dimens_m
865 use paramet_m
866 use comgeom
867 use iniprint
868 IMPLICIT NONE
869 c
870 c
871 c
872 c Arguments:
873 c ----------
874 real masse(ip1jmp1,llm)
875 real w_m( ip1jmp1,llm+1)
876 real q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)
877
878 c
879 c Local
880 c ---------
881 c
882 INTEGER ij,l
883 c
884 real new_m,zdq,zz
885 real zsigh(ip1jmp1,llm),zsigb(ip1jmp1,llm),zsig
886 real w_mq(ip1jmp1,llm+1)
887 real zm,zq,zsigm,zsigp,zqm,zqp
888 real prec
889 save prec
890
891 data prec/1.e-13/
892
893 do l=1,llm
894 do ij=1,ip1jmp1
895 zdq=qb(ij,l)-qh(ij,l)
896 c if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) then
897 c print*,'probleme au point ij=',ij,' l=',l
898 c print*,qh(ij,l),q(ij,l),qb(ij,l)
899 c qh(ij,l)=q(ij,l)
900 c qb(ij,l)=q(ij,l)
901 c endif
902
903 if(abs(zdq).gt.prec) then
904 zsigb(ij,l)=(q(ij,l)-qh(ij,l))/zdq
905 zsigh(ij,l)=1.-zsigb(ij,l)
906 zsigb(ij,l)=min(max(zsigb(ij,l),0.),1.)
907 else
908 zsigb(ij,l)=0.5
909 zsigh(ij,l)=0.5
910 endif
911 enddo
912 enddo
913
914 c print*,'ok1'
915 c calcul de la pente maximum dans la maille en valeur absolue
916 do l=2,llm
917 do ij=1,ip1jmp1
918 if (w_m(ij,l).ge.0.) then
919 zsigp=zsigb(ij,l)
920 zsigm=zsigh(ij,l)
921 zqp=qb(ij,l)
922 zqm=qh(ij,l)
923 zm=masse(ij,l)
924 zq=q(ij,l)
925 else
926 zsigm=zsigb(ij,l-1)
927 zsigp=zsigh(ij,l-1)
928 zqm=qb(ij,l-1)
929 zqp=qh(ij,l-1)
930 zm=masse(ij,l-1)
931 zq=q(ij,l-1)
932 endif
933 zsig=abs(w_m(ij,l))/zm
934 if(zsig.eq.0.) zsigp=0.1
935 if (zsig.le.zsigp) then
936 w_mq(ij,l)=w_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
937 else
938 zz=0.5*(zsig-zsigp)/zsigm
939 w_mq(ij,l)=sign(zm,w_m(ij,l))*( 0.5*(zq+zqp)*zsigp
940 s +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
941 endif
942 enddo
943 enddo
944
945 do ij=1,ip1jmp1
946 w_mq(ij,llm+1)=0.
947 w_mq(ij,1)=0.
948 enddo
949
950 do l=1,llm
951 do ij=1,ip1jmp1
952 new_m=masse(ij,l)+w_m(ij,l+1)-w_m(ij,l)
953 q(ij,l)=(q(ij,l)*masse(ij,l)+w_mq(ij,l+1)-w_mq(ij,l))
954 & /new_m
955 masse(ij,l)=new_m
956 enddo
957 enddo
958 c print*,'ok3'
959 RETURN
960 END
  ViewVC Help
Powered by ViewVC 1.1.21