source: LMDZ6/trunk/libf/dyn3d_common/advn.F90 @ 5272

Last change on this file since 5272 was 5272, checked in by abarral, 23 hours ago

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