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

Last change on this file since 5246 was 5246, checked in by abarral, 23 hours ago

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

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