source: LMDZ4/branches/V3_test/libf/dyn3dpar/pentes_ini.F @ 5415

Last change on this file since 5415 was 709, checked in by Laurent Fairhead, 18 years ago

Nouvelles versions de la dynamique YM
LF

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