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

Last change on this file since 5195 was 5159, checked in by abarral, 5 months ago

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