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

Last change on this file since 5122 was 5118, checked in by abarral, 4 months ago

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