source: LMDZ4/branches/V3_test/libf/dyn3d/pentes_ini.F @ 1839

Last change on this file since 1839 was 524, checked in by lmdzadmin, 20 years ago

Initial revision

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