GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: dyn3d_common/pentes_ini.F Lines: 0 207 0.0 %
Date: 2023-06-30 12:51:15 Branches: 0 128 0.0 %

Line Branch Exec Source
1
!
2
! $Header$
3
!
4
      SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode)
5
6
      USE comconst_mod, ONLY: pi, dtvr
7
8
      IMPLICIT NONE
9
10
c=======================================================================
11
c   Adaptation LMDZ:  A.Armengaud (LGGE)
12
c   ----------------
13
c
14
c   ********************************************************************
15
c   Transport des traceurs par la methode des pentes
16
c   ********************************************************************
17
c   Reference possible : Russel. G.L., Lerner J.A.:
18
c         A new Finite-Differencing Scheme for Traceur Transport
19
c         Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
20
c   ********************************************************************
21
c   q,w,masse,pbaru et pbarv
22
c                      sont des arguments d'entree  pour le s-pg ....
23
c
24
c=======================================================================
25
26
27
      include "dimensions.h"
28
      include "paramet.h"
29
      include "comgeom2.h"
30
31
c   Arguments:
32
c   ----------
33
      integer mode
34
      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
35
      REAL q( iip1,jjp1,llm,0:3)
36
      REAL w( ip1jmp1,llm )
37
      REAL masse( iip1,jjp1,llm)
38
c   Local:
39
c   ------
40
      LOGICAL limit
41
      REAL sm ( iip1,jjp1, llm )
42
      REAL s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )
43
      REAL sy( iip1,jjp1,llm ),  sz( iip1,jjp1,llm )
44
      real masn,mass,zz
45
      INTEGER i,j,l,iq
46
47
c  modif Fred 24 03 96
48
49
      real sinlon(iip1),sinlondlon(iip1)
50
      real coslon(iip1),coslondlon(iip1)
51
      save sinlon,coslon,sinlondlon,coslondlon
52
      real dyn1,dyn2,dys1,dys2
53
      real qpn,qps,dqzpn,dqzps
54
      real smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)
55
      real qmin,zq,pente_max
56
c
57
      REAL      SSUM
58
      integer ismax,ismin,lati,latf
59
      EXTERNAL  SSUM, ismin,ismax
60
      logical first
61
      save first
62
c   fin modif
63
64
c      EXTERNAL masskg
65
      EXTERNAL advx
66
      EXTERNAL advy
67
      EXTERNAL advz
68
69
c  modif Fred 24 03 96
70
      data first/.true./
71
72
      limit = .TRUE.
73
      pente_max=2
74
c     if (mode.eq.1.or.mode.eq.3) then
75
c     if (mode.eq.1) then
76
      if (mode.ge.1) then
77
        lati=2
78
        latf=jjm
79
      else
80
        lati=1
81
        latf=jjp1
82
      endif
83
84
      qmin=0.4995
85
      qmin=0.
86
      if(first) then
87
         print*,'SCHEMA AMONT NOUVEAU'
88
         first=.false.
89
         do i=2,iip1
90
            coslon(i)=cos(rlonv(i))
91
            sinlon(i)=sin(rlonv(i))
92
            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
93
            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
94
            print*,coslondlon(i),sinlondlon(i)
95
         enddo
96
         coslon(1)=coslon(iip1)
97
         coslondlon(1)=coslondlon(iip1)
98
         sinlon(1)=sinlon(iip1)
99
         sinlondlon(1)=sinlondlon(iip1)
100
         print*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1)
101
         print*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1)
102
        DO l = 1,llm
103
        DO j = 1,jjp1
104
         DO i = 1,iip1
105
         q ( i,j,l,1 )=0.
106
         q ( i,j,l,2 )=0.
107
         q ( i,j,l,3 )=0.
108
         ENDDO
109
         ENDDO
110
        ENDDO
111
112
      endif
113
c   Fin modif Fred
114
115
c *** q contient les qqtes de traceur avant l'advection
116
117
c *** Affectation des tableaux S a partir de Q
118
c *** Rem : utilisation de SCOPY ulterieurement
119
120
       DO l = 1,llm
121
        DO j = 1,jjp1
122
         DO i = 1,iip1
123
             s0( i,j,llm+1-l ) = q ( i,j,l,0 )
124
             sx( i,j,llm+1-l ) = q ( i,j,l,1 )
125
             sy( i,j,llm+1-l ) = q ( i,j,l,2 )
126
             sz( i,j,llm+1-l ) = q ( i,j,l,3 )
127
         ENDDO
128
        ENDDO
129
       ENDDO
130
131
c      PRINT*,'----- S0 just before conversion -------'
132
c      PRINT*,'S0(16,12,1)=',s0(16,12,1)
133
c      PRINT*,'Q(16,12,1,4)=',q(16,12,1,4)
134
135
c *** On calcule la masse d'air en kg
136
137
       DO  l = 1,llm
138
         DO  j = 1,jjp1
139
           DO  i = 1,iip1
140
            sm ( i,j,llm+1-l)=masse( i,j,l )
141
          ENDDO
142
         ENDDO
143
       ENDDO
144
145
c *** On converti les champs S en atome (resp. kg)
146
c *** Les routines d'advection traitent les champs
147
c *** a advecter si ces derniers sont en atome (resp. kg)
148
c *** A optimiser !!!
149
150
       DO  l = 1,llm
151
         DO  j = 1,jjp1
152
           DO  i = 1,iip1
153
               s0(i,j,l) = s0(i,j,l) * sm ( i,j,l )
154
               sx(i,j,l) = sx(i,j,l) * sm ( i,j,l )
155
               sy(i,j,l) = sy(i,j,l) * sm ( i,j,l )
156
               sz(i,j,l) = sz(i,j,l) * sm ( i,j,l )
157
           ENDDO
158
         ENDDO
159
       ENDDO
160
161
c       ss0 = 0.
162
c       DO l = 1,llm
163
c        DO j = 1,jjp1
164
c         DO i = 1,iim
165
c            ss0 = ss0 + s0 ( i,j,l )
166
c         ENDDO
167
c        ENDDO
168
c       ENDDO
169
c       PRINT*, 'valeur tot s0 avant advection=',ss0
170
171
c *** Appel des subroutines d'advection en X, en Y et en Z
172
c *** Advection avec "time-splitting"
173
174
c-----------------------------------------------------------
175
c      PRINT*,'----- S0 just before ADVX -------'
176
c      PRINT*,'S0(16,12,1)=',s0(16,12,1)
177
178
c-----------------------------------------------------------
179
c      do l=1,llm
180
c         do j=1,jjp1
181
c          do i=1,iip1
182
c             zq=s0(i,j,l)/sm(i,j,l)
183
c            if(zq.lt.qmin)
184
c    ,       print*,'avant advx1, s0(',i,',',j,',',l,')=',zq
185
c          enddo
186
c         enddo
187
c      enddo
188
CCC
189
       if(mode.eq.2) then
190
          do l=1,llm
191
            s0s=0.
192
            s0n=0.
193
            dyn1=0.
194
            dys1=0.
195
            dyn2=0.
196
            dys2=0.
197
            smn=0.
198
            sms=0.
199
            do i=1,iim
200
               smn=smn+sm(i,1,l)
201
               sms=sms+sm(i,jjp1,l)
202
               s0n=s0n+s0(i,1,l)
203
               s0s=s0s+s0(i,jjp1,l)
204
               zz=sy(i,1,l)/sm(i,1,l)
205
               dyn1=dyn1+sinlondlon(i)*zz
206
               dyn2=dyn2+coslondlon(i)*zz
207
               zz=sy(i,jjp1,l)/sm(i,jjp1,l)
208
               dys1=dys1+sinlondlon(i)*zz
209
               dys2=dys2+coslondlon(i)*zz
210
            enddo
211
            do i=1,iim
212
               sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i)
213
               sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i)
214
            enddo
215
            do i=1,iim
216
               s0(i,1,l)=s0n/smn+sy(i,1,l)
217
               s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l)
218
            enddo
219
220
            s0(iip1,1,l)=s0(1,1,l)
221
            s0(iip1,jjp1,l)=s0(1,jjp1,l)
222
223
            do i=1,iim
224
               sxn(i)=s0(i+1,1,l)-s0(i,1,l)
225
               sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l)
226
c   on rerentre les masses
227
            enddo
228
            do i=1,iim
229
               sy(i,1,l)=sy(i,1,l)*sm(i,1,l)
230
               sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l)
231
               s0(i,1,l)=s0(i,1,l)*sm(i,1,l)
232
               s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l)
233
            enddo
234
            sxn(iip1)=sxn(1)
235
            sxs(iip1)=sxs(1)
236
            do i=1,iim
237
               sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l)
238
               sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l)
239
            enddo
240
            s0(iip1,1,l)=s0(1,1,l)
241
            s0(iip1,jjp1,l)=s0(1,jjp1,l)
242
            sy(iip1,1,l)=sy(1,1,l)
243
            sy(iip1,jjp1,l)=sy(1,jjp1,l)
244
            sx(1,1,l)=sx(iip1,1,l)
245
            sx(1,jjp1,l)=sx(iip1,jjp1,l)
246
          enddo
247
      endif
248
249
      if (mode.eq.4) then
250
         do l=1,llm
251
            do i=1,iip1
252
               sx(i,1,l)=0.
253
               sx(i,jjp1,l)=0.
254
               sy(i,1,l)=0.
255
               sy(i,jjp1,l)=0.
256
            enddo
257
         enddo
258
      endif
259
      call limx(s0,sx,sm,pente_max)
260
c     call minmaxq(zq,1.e33,-1.e33,'avant advx     ')
261
       call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
262
c     call minmaxq(zq,1.e33,-1.e33,'avant advy     ')
263
      if (mode.eq.4) then
264
         do l=1,llm
265
            do i=1,iip1
266
               sx(i,1,l)=0.
267
               sx(i,jjp1,l)=0.
268
               sy(i,1,l)=0.
269
               sy(i,jjp1,l)=0.
270
            enddo
271
         enddo
272
      endif
273
       call   limy(s0,sy,sm,pente_max)
274
       call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
275
c     call minmaxq(zq,1.e33,-1.e33,'avant advz     ')
276
       do j=1,jjp1
277
          do i=1,iip1
278
             sz(i,j,1)=0.
279
             sz(i,j,llm)=0.
280
          enddo
281
       enddo
282
       call limz(s0,sz,sm,pente_max)
283
       call advz( limit,dtvr,w,sm,s0,sx,sy,sz )
284
      if (mode.eq.4) then
285
         do l=1,llm
286
            do i=1,iip1
287
               sx(i,1,l)=0.
288
               sx(i,jjp1,l)=0.
289
               sy(i,1,l)=0.
290
               sy(i,jjp1,l)=0.
291
            enddo
292
         enddo
293
      endif
294
        call limy(s0,sy,sm,pente_max)
295
       call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
296
       do l=1,llm
297
          do j=1,jjp1
298
             sm(iip1,j,l)=sm(1,j,l)
299
             s0(iip1,j,l)=s0(1,j,l)
300
             sx(iip1,j,l)=sx(1,j,l)
301
             sy(iip1,j,l)=sy(1,j,l)
302
             sz(iip1,j,l)=sz(1,j,l)
303
          enddo
304
       enddo
305
306
307
c     call minmaxq(zq,1.e33,-1.e33,'avant advx     ')
308
      if (mode.eq.4) then
309
         do l=1,llm
310
            do i=1,iip1
311
               sx(i,1,l)=0.
312
               sx(i,jjp1,l)=0.
313
               sy(i,1,l)=0.
314
               sy(i,jjp1,l)=0.
315
            enddo
316
         enddo
317
      endif
318
       call limx(s0,sx,sm,pente_max)
319
       call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
320
c     call minmaxq(zq,1.e33,-1.e33,'apres advx     ')
321
c      do l=1,llm
322
c         do j=1,jjp1
323
c          do i=1,iip1
324
c             zq=s0(i,j,l)/sm(i,j,l)
325
c            if(zq.lt.qmin)
326
c    ,       print*,'apres advx2, s0(',i,',',j,',',l,')=',zq
327
c          enddo
328
c         enddo
329
c      enddo
330
c ***   On repasse les S dans la variable q directement 14/10/94
331
c   On revient a des rapports de melange en divisant par la masse
332
333
c En dehors des poles:
334
335
       DO  l = 1,llm
336
        DO  j = 1,jjp1
337
         DO  i = 1,iim
338
             q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l)
339
             q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l)
340
             q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l)
341
             q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l)
342
         ENDDO
343
        ENDDO
344
      ENDDO
345
346
c Traitements specifiques au pole
347
348
      if(mode.ge.1) then
349
      DO l=1,llm
350
c   filtrages aux poles
351
         masn=ssum(iim,sm(1,1,l),1)
352
         mass=ssum(iim,sm(1,jjp1,l),1)
353
         qpn=ssum(iim,s0(1,1,l),1)/masn
354
         qps=ssum(iim,s0(1,jjp1,l),1)/mass
355
         dqzpn=ssum(iim,sz(1,1,l),1)/masn
356
         dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
357
         do i=1,iip1
358
            q( i,1,llm+1-l,3)=dqzpn
359
            q( i,jjp1,llm+1-l,3)=dqzps
360
            q( i,1,llm+1-l,0)=qpn
361
            q( i,jjp1,llm+1-l,0)=qps
362
         enddo
363
         if(mode.eq.3) then
364
            dyn1=0.
365
            dys1=0.
366
            dyn2=0.
367
            dys2=0.
368
            do i=1,iim
369
               dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l)
370
               dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l)
371
               dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
372
               dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
373
            enddo
374
            do i=1,iim
375
               q(i,1,llm+1-l,2)=
376
     s          (sinlon(i)*dyn1+coslon(i)*dyn2)
377
               q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
378
               q(i,jjp1,llm+1-l,2)=
379
     s          (sinlon(i)*dys1+coslon(i)*dys2)
380
               q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
381
     s         -q(i,jjp1,llm+1-l,2)
382
            enddo
383
         endif
384
         if(mode.eq.1) then
385
c   on filtre les valeurs au bord de la "grande maille pole"
386
            dyn1=0.
387
            dys1=0.
388
            dyn2=0.
389
            dys2=0.
390
            do i=1,iim
391
               zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
392
               dyn1=dyn1+sinlondlon(i)*zz
393
               dyn2=dyn2+coslondlon(i)*zz
394
               zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
395
               dys1=dys1+sinlondlon(i)*zz
396
               dys2=dys2+coslondlon(i)*zz
397
            enddo
398
            do i=1,iim
399
               q(i,1,llm+1-l,2)=
400
     s          (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
401
               q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
402
               q(i,jjp1,llm+1-l,2)=
403
     s          (sinlon(i)*dys1+coslon(i)*dys2)/2.
404
               q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
405
     s         -q(i,jjp1,llm+1-l,2)
406
            enddo
407
            q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
408
            q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
409
410
            do i=1,iim
411
               sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
412
               sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
413
            enddo
414
            sxn(iip1)=sxn(1)
415
            sxs(iip1)=sxs(1)
416
            do i=1,iim
417
               q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
418
               q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
419
            enddo
420
            q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
421
            q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1)
422
423
         endif
424
425
       ENDDO
426
       endif
427
428
c bouclage en longitude
429
      do iq=0,3
430
         do l=1,llm
431
            do j=1,jjp1
432
               q(iip1,j,l,iq)=q(1,j,l,iq)
433
            enddo
434
         enddo
435
      enddo
436
437
c       PRINT*, ' SORTIE DE PENTES ---  ca peut glisser ....'
438
439
        DO l = 1,llm
440
             DO j = 1,jjp1
441
              DO i = 1,iip1
442
                IF (q(i,j,l,0).lt.0.)  THEN
443
c                    PRINT*,'------------ BIP-----------'
444
c                    PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0)
445
c                    PRINT*,'QX(',i,j,l,')=',q(i,j,l,1)
446
c                    PRINT*,'QY(',i,j,l,')=',q(i,j,l,2)
447
c                    PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3)
448
c                            PRINT*,' PBL EN SORTIE DE PENTES'
449
                     q(i,j,l,0)=0.
450
c                    STOP
451
                 ENDIF
452
          ENDDO
453
         ENDDO
454
        ENDDO
455
456
c       PRINT*, '-------------------------------------------'
457
458
       do l=1,llm
459
          do j=1,jjp1
460
           do i=1,iip1
461
             if(q(i,j,l,0).lt.qmin)
462
     ,       print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)
463
           enddo
464
          enddo
465
       enddo
466
      RETURN
467
      END
468
469
470
471
472
473
474
475
476
477
478
479