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

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

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