source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/vlsplt_loc.f90 @ 5133

Last change on this file since 5133 was 5128, checked in by abarral, 5 months ago

Correct bug in vlspltqs_loc.f90 from r2270 where we call SSUM with incorrect arguments.
Merge the three different versions of abort_gcm into one
Fix seq, para 3D compilation broken from r5107 onwards
(lint) usual + Remove uneeded fixed-form continuations

  • 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:keywords set to Id
File size: 31.7 KB
Line 
1! $Id: vlsplt_loc.f90 5128 2024-07-25 15:47:25Z abarral $
2
3SUBROUTINE vlx_loc(q, pente_max, masse, u_m, ijb_x, ije_x, iq)
4
5  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
6  !
7  !    ********************************************************************
8  ! Shema  d'advection " pseudo amont " .
9  !    ********************************************************************
10  ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
11  !
12  !
13  !   --------------------------------------------------------------------
14  USE parallel_lmdz
15  USE infotrac, ONLY: nqtot, tracers, & ! CRisi                 &
16          min_qParent, min_qMass, min_ratio ! MVals et CRisi
17  USE lmdz_iniprint, ONLY: lunout, prt_level
18  IMPLICIT NONE
19  !
20  include "dimensions.h"
21  include "paramet.h"
22  !
23  !
24  !   Arguments:
25  !   ----------
26  REAL :: masse(ijb_u:ije_u, llm, nqtot), pente_max
27  REAL :: u_m(ijb_u:ije_u, llm), pbarv(iip1, jjb_v:jje_v, llm)
28  REAL :: q(ijb_u:ije_u, llm, nqtot) ! CRisi: ajout dimension nqtot
29  REAL :: w(ijb_u:ije_u, llm)
30  INTEGER :: iq ! CRisi
31  !
32  !  Local
33  !   ---------
34  !
35  INTEGER :: ij, l, j, i, iju, ijq, indu(ijnb_u), niju
36  INTEGER :: n0, iadvplus(ijb_u:ije_u, llm), nl(llm)
37  !
38  REAL :: new_m, zu_m, zdum(ijb_u:ije_u, llm)
39  REAL :: sigu(ijb_u:ije_u), dxq(ijb_u:ije_u, llm), dxqu(ijb_u:ije_u)
40  REAL :: zz(ijb_u:ije_u)
41  REAL :: adxqu(ijb_u:ije_u), dxqmax(ijb_u:ije_u, llm)
42  REAL :: u_mq(ijb_u:ije_u, llm)
43
44  REAL :: Ratio(ijb_u:ije_u, llm, nqtot) ! CRisi
45  INTEGER :: ifils, iq2 ! CRisi
46
47  Logical :: extremum
48
49  REAL :: z1, z2, z3
50
51  INTEGER :: ijb, ije, ijb_x, ije_x
52
53  !WRITE(*,*) 'vlsplt 58: entree dans vlx_loc, iq,ijb_x=',
54  ! &   iq,ijb_x
55  !   calcul de la pente a droite et a gauche de la maille
56
57  ijb = ijb_x
58  ije = ije_x
59
60  IF (pole_nord.AND.ijb==1) ijb = ijb + iip1
61  IF (pole_sud.AND.ije==ip1jmp1)  ije = ije - iip1
62
63  IF (pente_max>-1.e-5) THEN
64    ! IF (pente_max.gt.10) THEN
65
66    !   calcul des pentes avec limitation, Van Leer scheme I:
67    !   -----------------------------------------------------
68    ! on a besoin de q entre ijb et ije
69    !   calcul de la pente aux points u
70    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
71    DO l = 1, llm
72
73      DO ij = ijb, ije - 1
74        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
75        ! IF(u_m(ij,l).lt.0.) stop 'limx n admet pas les U<0'
76        ! sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
77      ENDDO
78      DO ij = ijb + iip1 - 1, ije, iip1
79        dxqu(ij) = dxqu(ij - iim)
80        ! sigu(ij)=sigu(ij-iim)
81      ENDDO
82
83      DO ij = ijb, ije
84        adxqu(ij) = abs(dxqu(ij))
85      ENDDO
86
87      !   calcul de la pente maximum dans la maille en valeur absolue
88
89      DO ij = ijb + 1, ije
90        dxqmax(ij, l) = pente_max * &
91                min(adxqu(ij - 1), adxqu(ij))
92        ! limitation subtile
93        !    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
94
95      ENDDO
96
97      DO ij = ijb + iip1 - 1, ije, iip1
98        dxqmax(ij - iim, l) = dxqmax(ij, l)
99      ENDDO
100
101      DO ij = ijb + 1, ije
102        IF(dxqu(ij - 1) * dxqu(ij)>0) THEN
103          dxq(ij, l) = dxqu(ij - 1) + dxqu(ij)
104        ELSE
105          !   extremum local
106          dxq(ij, l) = 0.
107        ENDIF
108        dxq(ij, l) = 0.5 * dxq(ij, l)
109        dxq(ij, l) = &
110                sign(min(abs(dxq(ij, l)), dxqmax(ij, l)), dxq(ij, l))
111      ENDDO
112
113    ENDDO ! l=1,llm
114    !$OMP END DO NOWAIT
115    ! PRINT*,'Ok calcul des pentes'
116
117  ELSE ! (pente_max.lt.-1.e-5)
118
119    !   Pentes produits:
120    !   ----------------
121    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
122    DO l = 1, llm
123      DO ij = ijb, ije - 1
124        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
125      ENDDO
126      DO ij = ijb + iip1 - 1, ije, iip1
127        dxqu(ij) = dxqu(ij - iim)
128      ENDDO
129
130      DO ij = ijb + 1, ije
131        zz(ij) = dxqu(ij - 1) * dxqu(ij)
132        zz(ij) = zz(ij) + zz(ij)
133        IF(zz(ij)>0) THEN
134          dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij))
135        ELSE
136          !   extremum local
137          dxq(ij, l) = 0.
138        ENDIF
139      ENDDO
140
141    ENDDO
142    !$OMP END DO NOWAIT
143  ENDIF ! (pente_max.lt.-1.e-5)
144
145  !WRITE(*,*) 'vlx 156: iq,ijb_x=',iq,ijb_x
146
147  !   bouclage de la pente en iip1:
148  !   -----------------------------
149  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
150  DO l = 1, llm
151    DO ij = ijb + iip1 - 1, ije, iip1
152      dxq(ij - iim, l) = dxq(ij, l)
153    ENDDO
154    DO ij = ijb, ije
155      iadvplus(ij, l) = 0
156    ENDDO
157
158  ENDDO
159  !$OMP END DO NOWAIT
160  ! PRINT*,'Bouclage en iip1'
161
162  !   calcul des flux a gauche et a droite
163
164
165  !   on cumule le flux correspondant a toutes les mailles dont la masse
166  !   au travers de la paroi pENDant le pas de temps.
167  ! PRINT*,'Cumule ....'
168  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
169  ! on a besoin de masse entre ijb et ije
170  DO l = 1, llm
171    DO ij = ijb, ije - 1
172      ! PRINT*,'masse(',ij,')=',masse(ij,l,iq)
173      IF (u_m(ij, l)>0.) THEN
174        zdum(ij, l) = 1. - u_m(ij, l) / masse(ij, l, iq)
175        u_mq(ij, l) = u_m(ij, l) * (q(ij, l, iq) &
176                + 0.5 * zdum(ij, l) * dxq(ij, l))
177      ELSE
178        zdum(ij, l) = 1. + u_m(ij, l) / masse(ij + 1, l, iq)
179        u_mq(ij, l) = u_m(ij, l) * (q(ij + 1, l, iq) &
180                - 0.5 * zdum(ij, l) * dxq(ij + 1, l))
181      ENDIF
182    ENDDO
183  ENDDO
184  !$OMP END DO NOWAIT
185
186  !   detection des points ou on advecte plus que la masse de la
187  !   maille
188  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
189  DO l = 1, llm
190    DO ij = ijb, ije - 1
191      IF(zdum(ij, l)<0) THEN
192        iadvplus(ij, l) = 1
193        u_mq(ij, l) = 0.
194      ENDIF
195    ENDDO
196  ENDDO
197  !$OMP END DO NOWAIT
198  ! PRINT*,'Ok test 1'
199
200  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
201  DO l = 1, llm
202    DO ij = ijb + iip1 - 1, ije, iip1
203      iadvplus(ij, l) = iadvplus(ij - iim, l)
204    ENDDO
205  ENDDO
206  !$OMP END DO NOWAIT
207  ! PRINT*,'Ok test 2'
208
209
210  !   traitement special pour le cas ou on advecte en longitude plus que le
211  !   contenu de la maille.
212  !   cette partie est mal vectorisee.
213
214  !  calcul du nombre de maille sur lequel on advecte plus que la maille.
215
216  n0 = 0
217  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
218  DO l = 1, llm
219    nl(l) = 0
220    DO ij = ijb, ije
221      nl(l) = nl(l) + iadvplus(ij, l)
222    ENDDO
223    n0 = n0 + nl(l)
224  ENDDO
225  !$OMP END DO NOWAIT
226  !ym      IF(n0.gt.1) THEN
227  !ym      IF(n0.gt.0) THEN
228
229  ! PRINT*,'Nombre de points pour lesquels on advect plus que le'
230  ! &       ,'contenu de la maille : ',n0
231  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
232
233  DO l = 1, llm
234    IF(nl(l)>0) THEN
235      iju = 0
236      !   indicage des mailles concernees par le traitement special
237      DO ij = ijb, ije
238        IF(iadvplus(ij, l)==1.AND.mod(ij, iip1)/=0) THEN
239          iju = iju + 1
240          indu(iju) = ij
241        ENDIF
242      ENDDO
243      niju = iju
244      !PRINT*,'vlx 278, niju,nl',niju,nl(l)
245
246      !  traitement des mailles
247      DO iju = 1, niju
248        ij = indu(iju)
249        j = (ij - 1) / iip1 + 1
250        zu_m = u_m(ij, l)
251        u_mq(ij, l) = 0.
252        IF(zu_m>0.) THEN
253          ijq = ij
254          i = ijq - (j - 1) * iip1
255          !   accumulation pour les mailles completements advectees
256          do while(zu_m>masse(ijq, l, iq))
257            u_mq(ij, l) = u_mq(ij, l) &
258                    + q(ijq, l, iq) * masse(ijq, l, iq)
259            zu_m = zu_m - masse(ijq, l, iq)
260            i = mod(i - 2 + iim, iim) + 1
261            ijq = (j - 1) * iip1 + i
262          ENDDO
263          !   ajout de la maille non completement advectee
264          u_mq(ij, l) = u_mq(ij, l) + zu_m * &
265                  (q(ijq, l, iq) + 0.5 * &
266                          (1. - zu_m / masse(ijq, l, iq)) * dxq(ijq, l))
267        ELSE
268          ijq = ij + 1
269          i = ijq - (j - 1) * iip1
270          !   accumulation pour les mailles completements advectees
271          do while(-zu_m>masse(ijq, l, iq))
272            u_mq(ij, l) = u_mq(ij, l) - q(ijq, l, iq) &
273                    * masse(ijq, l, iq)
274            zu_m = zu_m + masse(ijq, l, iq)
275            i = mod(i, iim) + 1
276            ijq = (j - 1) * iip1 + i
277          ENDDO
278          !   ajout de la maille non completement advectee
279          u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l, iq) - &
280                  0.5 * (1. + zu_m / masse(ijq, l, iq)) * dxq(ijq, l))
281        ENDIF
282      ENDDO
283    ENDIF
284  ENDDO
285  !$OMP END DO NOWAIT
286  !ym      ENDIF  ! n0.gt.0
287
288  !   bouclage en latitude
289  ! PRINT*,'Avant bouclage en latitude'
290  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
291  DO l = 1, llm
292    DO ij = ijb + iip1 - 1, ije, iip1
293      u_mq(ij, l) = u_mq(ij - iim, l)
294    ENDDO
295  ENDDO
296  !$OMP END DO NOWAIT
297
298  ! CRisi: appel récursif de l'advection sur les fils.
299  ! Il faut faire ça avant d'avoir mis à jour q et masse
300
301  do ifils = 1, tracers(iq)%nqDescen
302    ! attention: comme Ratio est utilisé comme q dans l'appel
303    ! recursif, il doit contenir à lui seul tous les indices de tous
304    ! les descendants!
305    iq2 = tracers(iq)%iqDescen(ifils)
306    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
307    DO l = 1, llm
308      DO ij = ijb, ije
309        ! On a besoin de q et masse seulement entre ijb et ije. On ne
310        ! les calcule donc que de ijb à ije
311        !MVals: veiller a ce qu'on n'ait pas de denominateur nul
312        masse(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
313        IF (q(ij, l, iq)>min_qParent) then ! modif 13 nov 2020
314          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
315        else
316          Ratio(ij, l, iq2) = min_ratio
317        endif
318      enddo
319    enddo
320    !$OMP END DO NOWAIT
321  enddo !do ifils=1,tracers(iq)%nqDescen
322  do ifils = 1, tracers(iq)%nqChildren
323    iq2 = tracers(iq)%iqDescen(ifils)
324    CALL vlx_loc(Ratio, pente_max, masse, u_mq, ijb_x, ije_x, iq2)
325  enddo
326  ! end CRisi
327
328
329  !   calcul des tENDances
330  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
331  DO l = 1, llm
332    DO ij = ijb + 1, ije
333      !MVals: veiller a ce qu'on n'ait pas de denominateur nul
334      new_m = max(masse(ij, l, iq) + u_m(ij - 1, l) - u_m(ij, l), min_qMass)
335      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + &
336              u_mq(ij - 1, l) - u_mq(ij, l)) &
337              / new_m
338      masse(ij, l, iq) = new_m
339    ENDDO
340    !   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
341    DO ij = ijb + iip1 - 1, ije, iip1
342      q(ij - iim, l, iq) = q(ij, l, iq)
343      masse(ij - iim, l, iq) = masse(ij, l, iq)
344    ENDDO
345  ENDDO
346  !$OMP END DO NOWAIT
347
348  ! retablir les fils en rapport de melange par rapport a l'air:
349  ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio
350  ! puis on boucle en longitude
351  do ifils = 1, tracers(iq)%nqDescen
352    iq2 = tracers(iq)%iqDescen(ifils)
353    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
354    DO l = 1, llm
355      DO ij = ijb + 1, ije
356        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
357      enddo
358      DO ij = ijb + iip1 - 1, ije, iip1
359        q(ij - iim, l, iq2) = q(ij, l, iq2)
360      enddo
361    enddo
362    !$OMP END DO NOWAIT
363  enddo
364
365END SUBROUTINE vlx_loc
366
367
368SUBROUTINE vly_loc(q, pente_max, masse, masse_adv_v, iq)
369  !
370  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
371  !
372  !    ********************************************************************
373  ! Shema  d'advection " pseudo amont " .
374  !    ********************************************************************
375  ! q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
376  ! dq         sont des arguments de sortie pour le s-pg ....
377  !
378  !
379  !   --------------------------------------------------------------------
380  USE parallel_lmdz
381  USE infotrac, ONLY: nqtot, tracers, & ! CRisi                 &
382          min_qParent, min_qMass, min_ratio ! MVals et CRisi
383  USE comconst_mod, ONLY: pi
384  USE lmdz_ssum_scopy, ONLY: ssum
385
386  IMPLICIT NONE
387  !
388  include "dimensions.h"
389  include "paramet.h"
390  include "comgeom.h"
391  !
392  !
393  !   Arguments:
394  !   ----------
395  REAL :: masse(ijb_u:ije_u, llm, nqtot), pente_max
396  REAL :: masse_adv_v(ijb_v:ije_v, llm)
397  REAL :: q(ijb_u:ije_u, llm, nqtot), dq(ijb_u:ije_u, llm)
398  INTEGER :: iq ! CRisi
399  !
400  !  Local
401  !   ---------
402  !
403  INTEGER :: i, ij, l
404  !
405  REAL :: airej2, airejjm, airescb(iim), airesch(iim)
406  REAL :: dyq(ijb_u:ije_u, llm), dyqv(ijb_v:ije_v), zdvm(ijb_u:ije_u, llm)
407  REAL :: adyqv(ijb_v:ije_v), dyqmax(ijb_u:ije_u)
408  REAL :: qbyv(ijb_v:ije_v, llm)
409
410  REAL :: qpns, qpsn, appn, apps, dyn1, dys1, dyn2, dys2, newmasse, fn, fs
411  ! REAL newq,oldmasse
412  Logical :: extremum, first
413  REAL :: temps0, temps1, temps2, temps3, temps4, temps5, second
414  SAVE temps0, temps1, temps2, temps3, temps4, temps5
415  !$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
416  SAVE first
417  !$OMP THREADPRIVATE(first)
418
419  REAL :: convpn, convps, convmpn, convmps
420  REAL :: massepn, masseps, qpn, qps
421  REAL :: sinlon(iip1), sinlondlon(iip1)
422  REAL :: coslon(iip1), coslondlon(iip1)
423  SAVE sinlon, coslon, sinlondlon, coslondlon
424  !$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
425  SAVE airej2, airejjm
426  !$OMP THREADPRIVATE(airej2,airejjm)
427
428  REAL :: Ratio(ijb_u:ije_u, llm, nqtot) ! CRisi
429  INTEGER :: ifils, iq2 ! CRisi
430
431  DATA first/.TRUE./
432  DATA temps0, temps1, temps2, temps3, temps4, temps5/0., 0., 0., 0., 0., 0./
433  INTEGER :: ijb, ije
434  INTEGER :: ijbm, ijem
435
436  ijb = ij_begin - 2 * iip1
437  ije = ij_end + 2 * iip1
438  IF (pole_nord) ijb = ij_begin
439  IF (pole_sud)  ije = ij_end
440
441  IF(first) THEN
442    PRINT*, 'Shema  Amont nouveau  appele dans  Vanleer   '
443    first = .FALSE.
444    do i = 2, iip1
445      coslon(i) = cos(rlonv(i))
446      sinlon(i) = sin(rlonv(i))
447      coslondlon(i) = coslon(i) * (rlonu(i) - rlonu(i - 1)) / pi
448      sinlondlon(i) = sinlon(i) * (rlonu(i) - rlonu(i - 1)) / pi
449    ENDDO
450    coslon(1) = coslon(iip1)
451    coslondlon(1) = coslondlon(iip1)
452    sinlon(1) = sinlon(iip1)
453    sinlondlon(1) = sinlondlon(iip1)
454    airej2 = SSUM(iim, aire(iip2), 1)
455    airejjm = SSUM(iim, aire(ip1jm - iim), 1)
456  ENDIF
457
458  !
459  ! PRINT*,'CALCUL EN LATITUDE'
460
461  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
462  DO l = 1, llm
463    !
464    !   --------------------------------
465    !  CALCUL EN LATITUDE
466    !   --------------------------------
467
468    !   On commence par calculer la valeur du traceur moyenne sur le premier cercle
469    !   de latitude autour du pole (qpns pour le pole nord et qpsn pour
470    !    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
471
472    IF (pole_nord) THEN
473      DO i = 1, iim
474        airescb(i) = aire(i + iip1) * q(i + iip1, l, iq)
475      ENDDO
476      qpns = SSUM(iim, airescb, 1) / airej2
477    endif
478
479    IF (pole_sud) THEN
480      DO i = 1, iim
481        airesch(i) = aire(i + ip1jm - iip1) * q(i + ip1jm - iip1, l, iq)
482      ENDDO
483      qpsn = SSUM(iim, airesch, 1) / airejjm
484    endif
485
486    !   calcul des pentes aux points v
487
488    ijb = ij_begin - 2 * iip1
489    ije = ij_end + iip1
490    IF (pole_nord) ijb = ij_begin
491    IF (pole_sud)  ije = ij_end - iip1
492
493    ! on a besoin de q entre ij_begin-2*iip1 et ij_end+2*iip1
494    ! Si pole sud, entre ij_begin-2*iip1 et ij_end
495    ! Si pole Nord, entre ij_begin et ij_end+2*iip1
496    DO ij = ijb, ije
497      dyqv(ij) = q(ij, l, iq) - q(ij + iip1, l, iq)
498      adyqv(ij) = abs(dyqv(ij))
499    ENDDO
500
501
502    !   calcul des pentes aux points scalaires
503    ijb = ij_begin - iip1
504    ije = ij_end + iip1
505    IF (pole_nord) ijb = ij_begin + iip1
506    IF (pole_sud)  ije = ij_end - iip1
507
508    DO ij = ijb, ije
509      dyq(ij, l) = .5 * (dyqv(ij - iip1) + dyqv(ij))
510      dyqmax(ij) = min(adyqv(ij - iip1), adyqv(ij))
511      dyqmax(ij) = pente_max * dyqmax(ij)
512    ENDDO
513
514    !   calcul des pentes aux poles
515    IF (pole_nord) THEN
516      DO ij = 1, iip1
517        dyq(ij, l) = qpns - q(ij + iip1, l, iq)
518      ENDDO
519
520      dyn1 = 0.
521      dyn2 = 0.
522      DO ij = 1, iim
523        dyn1 = dyn1 + sinlondlon(ij) * dyq(ij, l)
524        dyn2 = dyn2 + coslondlon(ij) * dyq(ij, l)
525      ENDDO
526      DO ij = 1, iip1
527        dyq(ij, l) = dyn1 * sinlon(ij) + dyn2 * coslon(ij)
528      ENDDO
529
530      DO ij = 1, iip1
531        dyq(ij, l) = 0.
532      ENDDO
533      ! ym tout cela ne sert pas a grand chose
534    ENDIF
535
536    IF (pole_sud) THEN
537
538      DO ij = 1, iip1
539        dyq(ip1jm + ij, l) = q(ip1jm + ij - iip1, l, iq) - qpsn
540      ENDDO
541
542      dys1 = 0.
543      dys2 = 0.
544
545      DO ij = 1, iim
546        dys1 = dys1 + sinlondlon(ij) * dyq(ip1jm + ij, l)
547        dys2 = dys2 + coslondlon(ij) * dyq(ip1jm + ij, l)
548      ENDDO
549
550      DO ij = 1, iip1
551        dyq(ip1jm + ij, l) = dys1 * sinlon(ij) + dys2 * coslon(ij)
552      ENDDO
553
554      DO ij = 1, iip1
555        dyq(ip1jm + ij, l) = 0.
556      ENDDO
557      ! ym tout cela ne sert pas a grand chose
558    ENDIF
559
560    !   filtrage de la derivee
561
562    !   calcul des pentes limites aux poles
563    ! ym partie inutile
564    ! goto 8888
565    ! fn=1.
566    ! fs=1.
567    ! DO ij=1,iim
568    !    IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
569    !       fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
570    !    ENDIF
571    ! IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
572    !    fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
573    !    ENDIF
574    ! ENDDO
575    ! DO ij=1,iip1
576    !    dyq(ij,l)=fn*dyq(ij,l)
577    !    dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
578    ! ENDDO
579    ! 8888    continue
580
581
582    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
583    !  En memoire de dIFferents tests sur la
584    !  limitation des pentes aux poles.
585    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
586    ! PRINT*,dyq(1)
587    ! PRINT*,dyqv(iip1+1)
588    ! appn=abs(dyq(1)/dyqv(iip1+1))
589    ! PRINT*,dyq(ip1jm+1)
590    ! PRINT*,dyqv(ip1jm-iip1+1)
591    ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
592    ! DO ij=2,iim
593    !    appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
594    !    apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
595    ! ENDDO
596    ! appn=min(pente_max/appn,1.)
597    ! apps=min(pente_max/apps,1.)
598    !
599    !
600    !   cas ou on a un extremum au pole
601    !
602    ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
603    !    &   appn=0.
604    ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
605    !    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
606    !    &   apps=0.
607    !
608    !   limitation des pentes aux poles
609    ! DO ij=1,iip1
610    !    dyq(ij)=appn*dyq(ij)
611    !    dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
612    ! ENDDO
613    !
614    !   test
615    !  DO ij=1,iip1
616    !     dyq(iip1+ij)=0.
617    !     dyq(ip1jm+ij-iip1)=0.
618    !  ENDDO
619    !  DO ij=1,ip1jmp1
620    !     dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
621    !  ENDDO
622    !
623    ! changement 10 07 96
624    ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
625    !    &   THEN
626    !    DO ij=1,iip1
627    !       dyqmax(ij)=0.
628    !    ENDDO
629    ! ELSE
630    !    DO ij=1,iip1
631    !       dyqmax(ij)=pente_max*abs(dyqv(ij))
632    !    ENDDO
633    ! ENDIF
634    !
635    ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
636    !    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
637    !    &THEN
638    !    DO ij=ip1jm+1,ip1jmp1
639    !       dyqmax(ij)=0.
640    !    ENDDO
641    ! ELSE
642    !    DO ij=ip1jm+1,ip1jmp1
643    !       dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
644    !    ENDDO
645    ! ENDIF
646    !   fin changement 10 07 96
647    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
648
649    !   calcul des pentes limitees
650    ijb = ij_begin - iip1
651    ije = ij_end + iip1
652    IF (pole_nord) ijb = ij_begin + iip1
653    IF (pole_sud)  ije = ij_end - iip1
654
655    DO ij = ijb, ije
656      IF(dyqv(ij) * dyqv(ij - iip1)>0.) THEN
657        dyq(ij, l) = sign(min(abs(dyq(ij, l)), dyqmax(ij)), dyq(ij, l))
658      ELSE
659        dyq(ij, l) = 0.
660      ENDIF
661    ENDDO
662
663  ENDDO
664  !$OMP END DO NOWAIT
665
666  ijb = ij_begin - iip1
667  ije = ij_end
668  IF (pole_nord) ijb = ij_begin
669  IF (pole_sud)  ije = ij_end - iip1
670
671  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
672  DO l = 1, llm
673    DO ij = ijb, ije
674      IF(masse_adv_v(ij, l)>0) THEN
675        qbyv(ij, l) = q(ij + iip1, l, iq) + dyq(ij + iip1, l) * &
676                0.5 * (1. - masse_adv_v(ij, l) &
677                / masse(ij + iip1, l, iq))
678      ELSE
679        qbyv(ij, l) = q(ij, l, iq) - dyq(ij, l) * &
680                0.5 * (1. + masse_adv_v(ij, l) / masse(ij, l, iq))
681      ENDIF
682      qbyv(ij, l) = masse_adv_v(ij, l) * qbyv(ij, l)
683    ENDDO
684  ENDDO
685  !$OMP END DO NOWAIT
686
687  ! CRisi: appel récursif de l'advection sur les fils.
688  ! Il faut faire ça avant d'avoir mis à jour q et masse
689  ! WRITE(*,*)'vly 689: iq,nqChildren(iq)=',iq,tracers(iq)%nqChildren
690
691  ijb = ij_begin - 2 * iip1
692  ije = ij_end + 2 * iip1
693  ijbm = ij_begin - iip1
694  ijem = ij_end + iip1
695  IF (pole_nord) ijb = ij_begin
696  IF (pole_sud)  ije = ij_end
697  IF (pole_nord) ijbm = ij_begin
698  IF (pole_sud)  ijem = ij_end
699
700  do ifils = 1, tracers(iq)%nqDescen
701    iq2 = tracers(iq)%iqDescen(ifils)
702    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
703    DO l = 1, llm
704      ! modif des bornes: CRisi 16 nov 2020
705      ! d'abord masse avec bornes corrigées
706      DO ij = ijbm, ijem
707        !MVals: veiller a ce qu'on n'ait pas de denominateur nul
708        masse(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
709      enddo
710
711      ! ensuite Ratio avec anciennes bornes
712      DO ij = ijb, ije
713        !MVals: veiller a ce qu'on n'ait pas de denominateur nul
714        IF (q(ij, l, iq)>min_qParent) then ! modif 13 nov 2020
715          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
716        else
717          Ratio(ij, l, iq2) = min_ratio
718        endif
719      enddo !DO ij=ijbm,ijem
720    enddo !DO l=1,llm
721    !$OMP END DO NOWAIT
722  enddo
723
724  do ifils = 1, tracers(iq)%nqChildren
725    iq2 = tracers(iq)%iqDescen(ifils)
726    CALL vly_loc(Ratio, pente_max, masse, qbyv, iq2)
727  enddo
728  ! end CRisi
729
730  ijb = ij_begin
731  ije = ij_end
732  IF (pole_nord) ijb = ij_begin + iip1
733  IF (pole_sud)  ije = ij_end - iip1
734
735  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
736  DO l = 1, llm
737    DO ij = ijb, ije
738      newmasse = masse(ij, l, iq) &
739              + masse_adv_v(ij, l) - masse_adv_v(ij - iip1, l)
740
741      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + qbyv(ij, l) &
742              - qbyv(ij - iip1, l)) / newmasse
743
744      masse(ij, l, iq) = newmasse
745
746    ENDDO
747
748    IF (pole_nord) THEN
749      convpn = SSUM(iim, qbyv(1, l), 1)
750      convmpn = ssum(iim, masse_adv_v(1, l), 1)
751      massepn = ssum(iim, masse(1, l, iq), 1)
752      qpn = 0.
753      do ij = 1, iim
754        qpn = qpn + masse(ij, l, iq) * q(ij, l, iq)
755      enddo
756      qpn = (qpn + convpn) / (massepn + convmpn)
757      do ij = 1, iip1
758        q(ij, l, iq) = qpn
759      enddo
760    endif
761
762    IF (pole_sud) THEN
763      convps = -SSUM(iim, qbyv(ip1jm - iim, l), 1)
764      convmps = -ssum(iim, masse_adv_v(ip1jm - iim, l), 1)
765      masseps = ssum(iim, masse(ip1jm + 1, l, iq), 1)
766      qps = 0.
767      do ij = ip1jm + 1, ip1jmp1 - 1
768        qps = qps + masse(ij, l, iq) * q(ij, l, iq)
769      enddo
770      qps = (qps + convps) / (masseps + convmps)
771      do ij = ip1jm + 1, ip1jmp1
772        q(ij, l, iq) = qps
773      enddo
774    endif
775  ENDDO
776  !$OMP END DO NOWAIT
777
778  ! retablir les fils en rapport de melange par rapport a l'air:
779  ijb = ij_begin
780  ije = ij_end
781  ! if (pole_nord) ijb=ij_begin
782  ! if (pole_sud)  ije=ij_end
783
784  do ifils = 1, tracers(iq)%nqDescen
785    iq2 = tracers(iq)%iqDescen(ifils)
786    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
787    DO l = 1, llm
788      DO ij = ijb, ije
789        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
790      enddo
791    enddo
792    !$OMP END DO NOWAIT
793  enddo
794
795END SUBROUTINE vly_loc
796
797
798SUBROUTINE vlz_loc(q, pente_max, masse, w, ijb_x, ije_x, iq)
799  !
800  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
801  !
802  !    ********************************************************************
803  ! Shema  d'advection " pseudo amont " .
804  !    ********************************************************************
805  !    q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
806  ! dq         sont des arguments de sortie pour le s-pg ....
807  !
808  !
809  !   --------------------------------------------------------------------
810  USE parallel_lmdz
811  USE vlz_mod
812  USE infotrac, ONLY: nqtot, tracers, & ! CRisi                 &
813          min_qParent, min_qMass, min_ratio ! MVals et CRisi
814  USE lmdz_iniprint, ONLY: lunout, prt_level
815
816
817  IMPLICIT NONE
818  !
819  include "dimensions.h"
820  include "paramet.h"
821  !
822  !
823  !   Arguments:
824  !   ----------
825  REAL :: masse(ijb_u:ije_u, llm, nqtot), pente_max
826  REAL :: q(ijb_u:ije_u, llm, nqtot)
827  REAL :: w(ijb_u:ije_u, llm + 1, nqtot)
828  INTEGER :: iq
829  !
830  !  Local
831  !   ---------
832  !
833  INTEGER :: i, ij, l, j, ii
834
835  REAL, DIMENSION(ijb_u:ije_u, llm + 1) :: wresi, morig, qorig, dzqorig
836  INTEGER, DIMENSION(ijb_u:ije_u, llm + 1) :: lorig
837  INTEGER, SAVE :: countcfl
838  !$OMP THREADPRIVATE(countcfl)
839  !
840  REAL :: newmasse
841
842  REAL :: dzqmax
843  REAL :: sigw
844
845  REAL :: temps0, temps1, temps2, temps3, temps4, temps5, second
846  SAVE temps0, temps1, temps2, temps3, temps4, temps5
847  !$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
848
849  DATA temps0, temps1, temps2, temps3, temps4, temps5/0., 0., 0., 0., 0., 0./
850  INTEGER :: ijb, ije, ijb_x, ije_x
851  LOGICAL, SAVE :: first = .TRUE.
852  !$OMP THREADPRIVATE(first)
853
854  !REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
855  ! Ces varibles doivent être déclarées en pointer et en save dans
856  ! vlz_loc si on veut qu'elles soient vues par tous les threads.
857  INTEGER :: ifils, iq2 ! CRisi
858
859  IF (first) THEN
860    first = .FALSE.
861  ENDIF
862  !    On oriente tout dans le sens de la pression c'est a dire dans le
863  !    sens de W
864
865  !WRITE(*,*) 'vlsplt 926: entree dans vlz_loc, iq=',iq
866
867  ijb = ijb_x
868  ije = ije_x
869
870  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
871  DO l = 2, llm
872    DO ij = ijb, ije
873      dzqw(ij, l) = q(ij, l - 1, iq) - q(ij, l, iq)
874      adzqw(ij, l) = abs(dzqw(ij, l))
875    ENDDO
876  ENDDO
877  !$OMP END DO
878
879  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
880  DO l = 2, llm - 1
881    DO ij = ijb, ije
882      IF(dzqw(ij, l) * dzqw(ij, l + 1)>0.) THEN
883        dzq(ij, l) = 0.5 * (dzqw(ij, l) + dzqw(ij, l + 1))
884      ELSE
885        dzq(ij, l) = 0.
886      ENDIF
887      dzqmax = pente_max * min(adzqw(ij, l), adzqw(ij, l + 1))
888      dzq(ij, l) = sign(min(abs(dzq(ij, l)), dzqmax), dzq(ij, l))
889    ENDDO
890  ENDDO
891  !$OMP END DO NOWAIT
892
893  !$OMP MASTER
894  DO ij = ijb, ije
895    dzq(ij, 1) = 0.
896    dzq(ij, llm) = 0.
897  ENDDO
898  !$OMP END MASTER
899  !$OMP BARRIER
900
901  !--------------------------------------------------------
902  ! On repere les points qui violent le CFL (|w| > masse)
903  !--------------------------------------------------------
904
905  countcfl = 0
906  ! PRINT*,'vlz nouveau'
907  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
908  DO l = 2, llm
909    DO ij = ijb, ije
910      IF((w(ij, l, iq)>0.AND.w(ij, l, iq)>masse(ij, l, iq)) &
911              .OR. (w(ij, l, iq)<=0.AND.ABS(w(ij, l, iq))>masse(ij, l - 1, iq))) &
912              countcfl = countcfl + 1
913    ENDDO
914  ENDDO
915  !$OMP END DO NOWAIT
916
917  ! ---------------------------------------------------------------
918  !  Identification des mailles ou on viole le CFL : w > masse
919  ! ---------------------------------------------------------------
920
921  IF (countcfl==0) THEN
922
923    ! ---------------------------------------------------------------
924    !   .... calcul des termes d'advection verticale  .......
925    ! Dans le cas où le  |w| < masse partout.
926    ! Version d'origine
927    ! Pourrait etre enleve si on voit que le code plus general
928    ! est aussi rapide
929    ! ---------------------------------------------------------------
930
931    ! calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
932
933    !  !WRITE(*,*) 'vlz 982,ijb,ije=',ijb,ije
934    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
935    DO l = 1, llm - 1
936      do  ij = ijb, ije
937        IF(w(ij, l + 1, iq)>0.) THEN
938          sigw = w(ij, l + 1, iq) / masse(ij, l + 1, iq)
939          wq(ij, l + 1, iq) = w(ij, l + 1, iq) * (q(ij, l + 1, iq) &
940                  + 0.5 * (1. - sigw) * dzq(ij, l + 1))
941        ELSE
942          sigw = w(ij, l + 1, iq) / masse(ij, l, iq)
943          wq(ij, l + 1, iq) = w(ij, l + 1, iq) * (q(ij, l, iq) &
944                  - 0.5 * (1. + sigw) * dzq(ij, l))
945        ENDIF
946      ENDDO
947    ENDDO
948    !$OMP END DO NOWAIT
949    !WRITE(*,*) 'vlz 1001'
950
951  ELSE ! countcfl>=1
952
953    IF (prt_level>9) THEN
954      WRITE(lunout, *)'vlz passage dans le non local'
955    ENDIF
956    ! ---------------------------------------------------------------
957    !  Debut du traitement du cas ou on viole le CFL : w > masse
958    ! ---------------------------------------------------------------
959
960    ! Initialisation
961
962    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
963    DO l = 2, llm
964      DO ij = ijb, ije
965        wresi(ij, l) = w(ij, l, iq)
966        wq(ij, l, iq) = 0.
967        IF(w(ij, l, iq)>0.) THEN
968          lorig(ij, l) = l
969          morig(ij, l) = masse(ij, l, iq)
970          qorig(ij, l) = q(ij, l, iq)
971          dzqorig(ij, l) = dzq(ij, l)
972        ELSE
973          lorig(ij, l) = l - 1
974          morig(ij, l) = masse(ij, l - 1, iq)
975          qorig(ij, l) = q(ij, l - 1, iq)
976          dzqorig(ij, l) = dzq(ij, l - 1)
977        ENDIF
978      ENDDO
979    ENDDO
980    !$OMP END DO NOWAIT
981
982    ! Reindicage vertical en accumulant les flux sur
983    !  les mailles qui viollent le CFL
984    !  on itère jusqu'à ce que tous les poins satisfassent
985    !  le critère
986    DO WHILE (countcfl>=1)
987      IF (prt_level>9) THEN
988        WRITE(lunout, *)'On viole le CFL Vertical sur ', countcfl, ' pts'
989      ENDIF
990      countcfl = 0
991
992      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
993      DO l = 2, llm
994        DO ij = ijb, ije
995          IF (ABS(wresi(ij, l))>morig(ij, l)) THEN
996            countcfl = countcfl + 1
997            ! rm : les 8 lignes ci dessous pourraient sans doute s'ecrire
998            ! avec la fonction sign
999            IF(w(ij, l, iq)>0.) THEN
1000              wresi(ij, l) = wresi(ij, l) - morig(ij, l)
1001              wq(ij, l, iq) = wq(ij, l, iq) + morig(ij, l) * qorig(ij, l)
1002              lorig(ij, l) = lorig(ij, l) + 1
1003            ELSE
1004              wresi(ij, l) = wresi(ij, l) + morig(ij, l)
1005              wq(ij, l, iq) = wq(ij, l, iq) - morig(ij, l) * qorig(ij, l)
1006              lorig(ij, l) = lorig(ij, l) - 1
1007            ENDIF
1008            ! CRisi 24nov2020: ajout d'un message d'erreur clair au lieu d'un plantage
1009            ! pour seg fault
1010            IF (lorig(ij, l)==0) THEN
1011              CALL abort_gcm("vlz in vlsplt_loc", &
1012                      "unfixable violation of CFL", 1)
1013            endif
1014            morig(ij, l) = masse(ij, lorig(ij, l), iq)
1015            qorig(ij, l) = q(ij, lorig(ij, l), iq)
1016            dzqorig(ij, l) = dzq(ij, lorig(ij, l))
1017          ENDIF
1018        ENDDO
1019      ENDDO
1020      !$OMP END DO NOWAIT
1021
1022    ENDDO ! WHILE (countcfl>=1)
1023
1024    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1025    DO l = 2, llm
1026      do  ij = ijb, ije
1027        sigw = wresi(ij, l) / morig(ij, l)
1028        IF(w(ij, l, iq)>0.) THEN
1029          wq(ij, l, iq) = wq(ij, l, iq) + wresi(ij, l) * (qorig(ij, l) &
1030                  + 0.5 * (1. - sigw) * dzqorig(ij, l))
1031        ELSE
1032          wq(ij, l, iq) = wq(ij, l, iq) + wresi(ij, l) * (qorig(ij, l) &
1033                  - 0.5 * (1. + sigw) * dzqorig(ij, l))
1034        ENDIF
1035      ENDDO
1036    ENDDO
1037    !$OMP END DO NOWAIT
1038
1039  ENDIF ! councfl=0
1040
1041
1042
1043  !$OMP MASTER
1044  DO ij = ijb, ije
1045    wq(ij, llm + 1, iq) = 0.
1046    wq(ij, 1, iq) = 0.
1047  ENDDO
1048  !$OMP END MASTER
1049  !$OMP BARRIER
1050
1051  ! CRisi: appel récursif de l'advection sur les fils.
1052  ! Il faut faire ça avant d'avoir mis à jour q et masse
1053  ! WRITE(*,*)'vlsplt 942: iq,nqChildren(iq)=',iq,tracers(iq)%nqChildren
1054  do ifils = 1, tracers(iq)%nqDescen
1055    iq2 = tracers(iq)%iqDescen(ifils)
1056    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1057    DO l = 1, llm
1058      DO ij = ijb, ije
1059        !MVals: veiller a ce qu'on n'ait pas de denominateur nul
1060        masse(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
1061        IF (q(ij, l, iq)>min_qParent) THEN
1062          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
1063        else
1064          Ratio(ij, l, iq2) = min_ratio
1065        endif
1066        !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015
1067        w(ij, l, iq2) = wq(ij, l, iq)
1068      enddo
1069    enddo
1070    !$OMP END DO NOWAIT
1071  enddo
1072  !$OMP BARRIER
1073
1074  do ifils = 1, tracers(iq)%nqChildren
1075    iq2 = tracers(iq)%iqDescen(ifils)
1076    CALL vlz_loc(Ratio, pente_max, masse, w, ijb_x, ije_x, iq2)
1077  enddo
1078  ! end CRisi
1079
1080  ! CRisi: On rajoute ici une barrière car on veut être sur que tous les
1081  ! wq soient synchronisés
1082
1083  !$OMP BARRIER
1084  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1085  DO l = 1, llm
1086    DO ij = ijb, ije
1087      newmasse = masse(ij, l, iq) + w(ij, l + 1, iq) - w(ij, l, iq)
1088      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) &
1089              + wq(ij, l + 1, iq) - wq(ij, l, iq)) &
1090              / newmasse
1091      masse(ij, l, iq) = newmasse
1092    ENDDO
1093  ENDDO
1094  !$OMP END DO NOWAIT
1095
1096
1097  ! retablir les fils en rapport de melange par rapport a l'air:
1098  do ifils = 1, tracers(iq)%nqDescen
1099    iq2 = tracers(iq)%iqDescen(ifils)
1100    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1101    DO l = 1, llm
1102      DO ij = ijb, ije
1103        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
1104      enddo
1105    enddo
1106    !$OMP END DO NOWAIT
1107  enddo
1108
1109END SUBROUTINE vlz_loc
1110
1111
1112
1113
Note: See TracBrowser for help on using the repository browser.