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

Last change on this file since 5423 was 5159, checked in by abarral, 5 months ago

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