source: LMDZ5/branches/testing/libf/dyn3d_common/pentes_ini.F @ 3946

Last change on this file since 3946 was 2641, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2593:2640 into testing branch

  • 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: 13.5 KB
RevLine 
[524]1!
2! $Header$
3!
4      SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode)
[2641]5     
6      USE comconst_mod, ONLY: pi, dtvr
7     
[524]8      IMPLICIT NONE
9
10c=======================================================================
11c   Adaptation LMDZ:  A.Armengaud (LGGE)
12c   ----------------
13c
14c   ********************************************************************
15c   Transport des traceurs par la methode des pentes
16c   ********************************************************************
17c   Reference possible : Russel. G.L., Lerner J.A.:
18c         A new Finite-Differencing Scheme for Traceur Transport
19c         Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
20c   ********************************************************************
21c   q,w,masse,pbaru et pbarv
22c                      sont des arguments d'entree  pour le s-pg ....
23c
24c=======================================================================
25
26
[2641]27      include "dimensions.h"
28      include "paramet.h"
29      include "comgeom2.h"
[524]30
31c   Arguments:
32c   ----------
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)
38c   Local:
39c   ------
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
47c  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
56c
57      REAL      SSUM
58      integer ismax,ismin,lati,latf
[1952]59      EXTERNAL  SSUM, ismin,ismax
[524]60      logical first
61      save first
62c   fin modif
63
64c      EXTERNAL masskg
65      EXTERNAL advx
66      EXTERNAL advy
67      EXTERNAL advz
68
69c  modif Fred 24 03 96
70      data first/.true./
71
72      limit = .TRUE.
73      pente_max=2
74c     if (mode.eq.1.or.mode.eq.3) then
75c     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
113c   Fin modif Fred
114
115c *** q contient les qqtes de traceur avant l'advection
116
117c *** Affectation des tableaux S a partir de Q
118c *** 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
131c      PRINT*,'----- S0 just before conversion -------'
132c      PRINT*,'S0(16,12,1)=',s0(16,12,1)
133c      PRINT*,'Q(16,12,1,4)=',q(16,12,1,4)
134
135c *** 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
145c *** On converti les champs S en atome (resp. kg)
146c *** Les routines d'advection traitent les champs
147c *** a advecter si ces derniers sont en atome (resp. kg)
148c *** 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
161c       ss0 = 0.
162c       DO l = 1,llm
163c        DO j = 1,jjp1
164c         DO i = 1,iim
165c            ss0 = ss0 + s0 ( i,j,l )
166c         ENDDO
167c        ENDDO
168c       ENDDO
169c       PRINT*, 'valeur tot s0 avant advection=',ss0
170
171c *** Appel des subroutines d'advection en X, en Y et en Z
172c *** Advection avec "time-splitting"
173     
174c-----------------------------------------------------------
175c      PRINT*,'----- S0 just before ADVX -------'
176c      PRINT*,'S0(16,12,1)=',s0(16,12,1)
177
178c-----------------------------------------------------------
179c      do l=1,llm
180c         do j=1,jjp1
181c          do i=1,iip1
182c             zq=s0(i,j,l)/sm(i,j,l)
183c            if(zq.lt.qmin)
184c    ,       print*,'avant advx1, s0(',i,',',j,',',l,')=',zq
185c          enddo
186c         enddo
187c      enddo
188CCC
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)
226c   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)
260c     call minmaxq(zq,1.e33,-1.e33,'avant advx     ')
261       call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
262c     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 )
275c     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
307c     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)
320c     call minmaxq(zq,1.e33,-1.e33,'apres advx     ')
321c      do l=1,llm
322c         do j=1,jjp1
323c          do i=1,iip1
324c             zq=s0(i,j,l)/sm(i,j,l)
325c            if(zq.lt.qmin)
326c    ,       print*,'apres advx2, s0(',i,',',j,',',l,')=',zq
327c          enddo
328c         enddo
329c      enddo
330c ***   On repasse les S dans la variable q directement 14/10/94
331c   On revient a des rapports de melange en divisant par la masse
332
333c 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
346c Traitements specifiques au pole
347
348      if(mode.ge.1) then
349      DO l=1,llm
350c   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
385c   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
428c 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
437c       PRINT*, ' SORTIE DE PENTES ---  ca peut glisser ....'
438
439        DO l = 1,llm
[2641]440             DO j = 1,jjp1
441              DO i = 1,iip1
[524]442                IF (q(i,j,l,0).lt.0.)  THEN
443c                    PRINT*,'------------ BIP-----------'
444c                    PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0)
445c                    PRINT*,'QX(',i,j,l,')=',q(i,j,l,1)
446c                    PRINT*,'QY(',i,j,l,')=',q(i,j,l,2)
447c                    PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3)
[2641]448c                            PRINT*,' PBL EN SORTIE DE PENTES'
[524]449                     q(i,j,l,0)=0.
450c                    STOP
451                 ENDIF
452          ENDDO
453         ENDDO
454        ENDDO
455
456c       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
Note: See TracBrowser for help on using the repository browser.