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

Last change on this file since 5300 was 5285, checked in by abarral, 4 days ago

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