source: LMDZ6/trunk/libf/dyn3d_common/pentes_ini.f90 @ 5450

Last change on this file since 5450 was 5285, checked in by abarral, 2 months ago

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