source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advn.f90 @ 5117

Last change on this file since 5117 was 5117, checked in by abarral, 3 months ago

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