source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/pentes_ini.f90 @ 5136

Last change on this file since 5136 was 5136, checked in by abarral, 8 weeks ago

Put comgeom.h, comgeom2.h into modules

  • 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: 12.2 KB
Line 
1
2! $Header$
3
4SUBROUTINE pentes_ini(q,w,masse,pbaru,pbarv,mode)
5
6  USE comconst_mod, ONLY: pi, dtvr
7  USE lmdz_ssum_scopy, ONLY: ssum
8  USE lmdz_comgeom2
9
10  IMPLICIT NONE
11
12  !=======================================================================
13  !   Adaptation LMDZ:  A.Armengaud (LGGE)
14  !   ----------------
15  !
16  !   ********************************************************************
17  !   Transport des traceurs par la methode des pentes
18  !   ********************************************************************
19  !   Reference possible : Russel. G.L., Lerner J.A.:
20  !     A new Finite-Differencing Scheme for Traceur Transport
21  !     Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
22  !   ********************************************************************
23  !   q,w,masse,pbaru et pbarv
24  !                  sont des arguments d'entree  pour le s-pg ....
25  !
26  !=======================================================================
27
28
29  INCLUDE "dimensions.h"
30  INCLUDE "paramet.h"
31
32  !   Arguments:
33  !   ----------
34  INTEGER :: mode
35  REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
36  REAL :: q( iip1,jjp1,llm,0:3)
37  REAL :: w( ip1jmp1,llm )
38  REAL :: masse( iip1,jjp1,llm)
39  !   Local:
40  !   ------
41  LOGICAL :: limit
42  REAL :: sm ( iip1,jjp1, llm )
43  REAL :: s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )
44  REAL :: sy( iip1,jjp1,llm ),  sz( iip1,jjp1,llm )
45  REAL :: masn,mass,zz
46  INTEGER :: i,j,l,iq
47
48  !  modif Fred 24 03 96
49
50  REAL :: sinlon(iip1),sinlondlon(iip1)
51  REAL :: coslon(iip1),coslondlon(iip1)
52  save sinlon,coslon,sinlondlon,coslondlon
53  REAL :: dyn1,dyn2,dys1,dys2
54  REAL :: qpn,qps,dqzpn,dqzps
55  REAL :: smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)
56  REAL :: qmin,zq,pente_max
57  !
58  INTEGER :: lati,latf
59  LOGICAL :: first
60  save first
61  !   fin modif
62
63   ! EXTERNAL masskg
64  EXTERNAL advx
65  EXTERNAL advy
66  EXTERNAL advz
67
68  !  modif Fred 24 03 96
69  data first/.TRUE./
70
71  limit = .TRUE.
72  pente_max=2
73  ! if (mode.EQ.1.OR.mode.EQ.3) THEN
74  ! if (mode.EQ.1) THEN
75  IF (mode>=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
112  !   Fin modif Fred
113
114  ! *** q contient les qqtes de traceur avant l'advection
115
116  ! *** Affectation des tableaux S a partir de Q
117  ! *** 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
130   ! PRINT*,'----- S0 just before conversion -------'
131   ! PRINT*,'S0(16,12,1)=',s0(16,12,1)
132   ! PRINT*,'Q(16,12,1,4)=',q(16,12,1,4)
133
134  ! *** 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
144  ! *** On converti les champs S en atome (resp. kg)
145  ! *** Les routines d'advection traitent les champs
146  ! *** a advecter si ces derniers sont en atome (resp. kg)
147  ! *** 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
160    ! ss0 = 0.
161    ! DO l = 1,llm
162    !  DO j = 1,jjp1
163    !   DO i = 1,iim
164    !      ss0 = ss0 + s0 ( i,j,l )
165    !   ENDDO
166    !  ENDDO
167    ! ENDDO
168    ! PRINT*, 'valeur tot s0 avant advection=',ss0
169
170  ! *** Appel des subroutines d'advection en X, en Y et en Z
171  ! *** Advection avec "time-splitting"
172
173  !-----------------------------------------------------------
174   ! PRINT*,'----- S0 just before ADVX -------'
175   ! PRINT*,'S0(16,12,1)=',s0(16,12,1)
176
177  !-----------------------------------------------------------
178   ! do l=1,llm
179   !    do j=1,jjp1
180   !     do i=1,iip1
181   !        zq=s0(i,j,l)/sm(i,j,l)
182   !       IF(zq.lt.qmin)
183  !    ,       PRINT*,'avant advx1, s0(',i,',',j,',',l,')=',zq
184   !     enddo
185   !    enddo
186   ! enddo
187  !CC
188   IF(mode==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)
225  !   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==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)
259  ! CALL minmaxq(zq,1.e33,-1.e33,'avant advx     ')
260   CALL advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
261  ! CALL minmaxq(zq,1.e33,-1.e33,'avant advy     ')
262  IF (mode==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 )
274  ! 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==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
306  ! CALL minmaxq(zq,1.e33,-1.e33,'avant advx     ')
307  IF (mode==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)
319  ! CALL minmaxq(zq,1.e33,-1.e33,'apres advx     ')
320  !  do l=1,llm
321  !     do j=1,jjp1
322  !      do i=1,iip1
323  !         zq=s0(i,j,l)/sm(i,j,l)
324  !        IF(zq.lt.qmin)
325  !    ,       PRINT*,'apres advx2, s0(',i,',',j,',',l,')=',zq
326  !      enddo
327  !     enddo
328  !  enddo
329  ! ***   On repasse les S dans la variable q directement 14/10/94
330  !   On revient a des rapports de melange en divisant par la masse
331
332  ! 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
345  ! Traitements specifiques au pole
346
347  IF(mode>=1) THEN
348  DO l=1,llm
349  !   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==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                 (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                 (sinlon(i)*dys1+coslon(i)*dys2)
379           q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) &
380                 -q(i,jjp1,llm+1-l,2)
381        enddo
382     endif
383     IF(mode==1) THEN
384  !   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                 (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                 (sinlon(i)*dys1+coslon(i)*dys2)/2.
403           q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) &
404                 -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
427  ! 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
436    ! 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)<0.)  THEN
442                 ! PRINT*,'------------ BIP-----------'
443                 ! PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0)
444                 ! PRINT*,'QX(',i,j,l,')=',q(i,j,l,1)
445                 ! PRINT*,'QY(',i,j,l,')=',q(i,j,l,2)
446                 ! PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3)
447                 !         PRINT*,' PBL EN SORTIE DE PENTES'
448                 q(i,j,l,0)=0.
449                 ! STOP
450             ENDIF
451      ENDDO
452     ENDDO
453    ENDDO
454
455    ! PRINT*, '-------------------------------------------'
456
457   do l=1,llm
458      do j=1,jjp1
459       do i=1,iip1
460         IF(q(i,j,l,0)<qmin) &
461               PRINT*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)
462       enddo
463      enddo
464   enddo
465  RETURN
466END SUBROUTINE pentes_ini
467
468
469
470
471
472
473
474
475
476
477
478
Note: See TracBrowser for help on using the repository browser.