source: trunk/LMDZ.GENERIC/libf/dyn3d/ps_amontF @ 1422

Last change on this file since 1422 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

File size: 13.3 KB
Line 
1      SUBROUTINE ps_amont(nq,iq,q,w,pbaru,pbarv,dq)
2c
3c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
4c
5c    ********************************************************************
6c     Shema  d'advection " pseudo amont " .
7c    ********************************************************************
8c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
9c     dq               sont des arguments de sortie pour le s-pg ....
10c
11c
12c   --------------------------------------------------------------------
13      IMPLICIT NONE
14c
15#include "dimensions.h"
16#include "paramet.h"
17#include "comgeom.h"
18c
19c
20c   Arguments:
21c   ----------
22      INTEGER nq,iq
23      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
24      REAL q(ip1jmp1,llm,nq), dq( ip1jmp1,llm,nq )
25      REAL w(ip1jmp1,llm)
26c
27c      Local
28c   ---------
29c
30      INTEGER i,ij,l
31c
32      REAL airej2,airejjm,airescb(iim),airesch(iim)
33      REAL pente(ip1jmp1),xg(ip1jmp1),xd(ip1jmp1),xs(ip1jmp1) ,
34     *     xn(ip1jmp1),xb(ip1jmp1),xh(ip1jmp1)
35      REAL qg(ip1jmp1),qd(ip1jmp1),qs(ip1jmp1) ,
36     *     qn(ip1jmp1),qb(ip1jmp1,llm),qh(ip1jmp1,llm)
37      REAL qbyv(ip1jm,llm), qbxu(ip1jmp1,llm), ww,dqh(ip1jmp1,llm)
38      REAL qpns,qpsn
39      logical first,extrpn,extrps
40      save first
41c
42c
43      REAL      SSUM,CVMGP,CVMGT
44      EXTERNAL  SSUM, convflu
45
46      data first/.true./
47
48      if(first) then
49         print*,'SCHEMA AMONT NOUVEAU'
50         first=.false.
51      endif
52
53
54c
55c
56      IF( forward.OR.leapf )   THEN
57c
58c
59      DO  100  l = 1, llm
60c
61c   ...  Boucle sur les  llm niveaux verticaux ...
62c
63c
64c  --------------------------------------------------------------
65c  --------------------------------------------------------------
66c  .............      Traitement en longitude     ...............
67c  --------------------------------------------------------------
68c  --------------------------------------------------------------
69c
70c     
71c        |          |            |           |
72c        |   q(i-1) |     q(i)   |   q(i+1)  |
73c        |          |qg(i)  qd(i)|qg(i+1)    |
74c
75c
76c      En  longitude ,
77c      Pour chaque maille ( i ) avec q(i,j,l,iq), on cherche a determiner
78c      avec une methode de ' pente' les valeurs qg(i) et qd(i) qui se trouvent
79c      au bord gauche et droite de cette maille .
80c
81c      Si ( q(i+1)-q(i) ) * ( q(i)-q(i-1)) < 0. ,on a qg(i)=qd(i)=q(i)
82c              Sinon
83c      qg(i)= q(i) - 1/4 * ( q(i+1) - q(i-1))
84c      qd(i)= q(i) + 1/4 * ( q(i+1) - q(i-1) )
85c
86c      On utilisera la meme methode pour determiner les valeurs qs(i) et qn(i)
87c      en latitude , ainsi que les valeurs qb(i) et qh(i) en  altitude .
88c
89c
90c
91       DO      ij    = 1,iip1
92       qg(ij)        = 0.
93       qd(ij)        = 0.
94       qg(ij+ ip1jm) = 0.
95       qd(ij+ ip1jm) = 0.
96       ENDDO
97c
98c     ....  calculs pour les lignes j= 2 a j = jjm  ....
99c
100       DO   ij    = iip2, ip1jm -1
101       pente(ij)  =( q(ij+1,l,iq)-q(ij,l,iq)) *(q(ij,l,iq)-q(ij-1,l,iq))
102       xg(ij)     = q(ij,l,iq) - 0.25 * ( q(ij+1,l,iq) - q(ij-1,l,iq) )
103       xd(ij)     = q(ij,l,iq) + 0.25 * ( q(ij+1,l,iq) - q(ij-1,l,iq) )
104       qg(ij)     = CVMGP( xg(ij), q(ij,l,iq) ,pente(ij) )
105       qd(ij)     = CVMGP( xd(ij), q(ij,l,iq), pente(ij) )
106       ENDDO
107
108c    ... Correction aux points  ( i= 1,  j )  .....
109c
110       DO   ij   = iip2, ip1jm, iip1
111       pente(ij) = ( q(ij+1,l,iq) - q(ij,l,iq) )     *
112     *                                ( q(ij,l,iq)  - q(ij+iim-1,l,iq) )
113       xg(ij)    = q(ij,l,iq) - 0.25* ( q(ij+1,l,iq)- q(ij+iim-1,l,iq) )
114       xd(ij)    = q(ij,l,iq) + 0.25* ( q(ij+1,l,iq)- q(ij+iim-1,l,iq) )
115       qg(ij)    = CVMGP( xg(ij), q(ij,l,iq) ,pente(ij) )
116       qd(ij)    = CVMGP( xd(ij), q(ij,l,iq), pente(ij) )
117       ENDDO
118c
119c    ... Correction aux points ( i= iip1,  j ) .....
120c
121       DO     ij      = iip2, ip1jm, iip1
122       qg( ij+ iim  ) = qg( ij )
123       qd( ij+ iim )  = qd( ij )
124       ENDDO
125c
126c   .............................................................
127c   .........    Limitation des pentes a gauche des boites  .....
128c
129c    Si (q(i)-qg(i))*(qg(i)-q(i-1)) < 0.  , on a  qg(i)=q(i-1) 
130c                                         et qd(i)=2*q(i)-qg(i)
131c   .............................................................
132c
133       DO   ij  = iip2,ip1jm -1
134       pente(ij)= ( qg(ij) -q(ij-1,l,iq))*(q(ij,l,iq)-qg(ij) )
135       qg(ij)   = CVMGP( qg(ij), q(ij-1,l,iq) ,pente(ij) )
136       qd(ij)   = CVMGP( qd(ij),
137     *                q(ij,l,iq)+ q(ij,l,iq) -qg(ij) , pente(ij) )
138       ENDDO
139c
140c      .....   Correction aux points ( i= 1,  j )  ......
141c
142       DO    ij  = iip2 ,ip1jm, iip1
143       qg(ij)    =  qg(ij+ iim)
144       qd(ij)    =  qd(ij+ iim)
145       ENDDO
146c   
147c      ...............................................................
148c      ...... Limitation des pentes a droite des boites      .........
149c      Si (q(i+1)-qd(i))*(qd(i)-q(i)) < 0. , on a qd(i)=q(i+1)
150c                                              et qg(i)=2*q(i)-qd(i) .
151c      ...............................................................
152c
153       DO ij = iip2, ip1jm -1
154       pente(ij) = ( qd(ij)-q(ij,l,iq) )*(q(ij+1,l,iq)-qd(ij) )
155       qd(ij)    =  CVMGP( qd(ij), q(ij+1,l,iq), pente(ij) )
156       qg(ij)    =  CVMGP( qg(ij),
157     *                q(ij,l,iq)+ q(ij,l,iq) -qd(ij) , pente(ij) )
158       ENDDO
159c
160c      ....  Correction aux points ( i = iip1, j )  .....
161c
162       DO     ij     = iip2, ip1jm, iip1
163       qg( ij+ iim ) = qg( ij )
164       qd( ij+ iim ) = qd( ij )
165       ENDDO
166c
167
168c     -------------------------------------------------------------
169c     -------------------------------------------------------------
170c     .............    Traitement en  latitude    .................
171c     -------------------------------------------------------------
172c     -------------------------------------------------------------
173c
174c
175c      q(j=1)   PN
176c   --------------
177c     --------- 
178c   --------------
179c
180c      q(j-1)
181c
182c   --------------
183c            qn(j)
184c      q(j)
185c            qs(j)
186c   --------------
187c       
188c      q(j+1)
189c
190c   --------------
191c
192c      q(jjp1)  PS
193c
194c   --------------
195c
196c
197c     ......   operations pour les lignes j= 2 a j= jjm  .......
198c
199       DO    ij    = iip2, ip1jm
200       pente(ij)   = ( q(ij-iip1,l,iq)- q(ij,l,iq)  )  *
201     *               ( q(ij,l,iq) - q(ij+iip1,l,iq) )
202       xs(ij) = q(ij,l,iq) - 0.25 * ( q(ij-iip1,l,iq) -q(ij+iip1,l,iq) )
203       xn(ij) = q(ij,l,iq) + 0.25 * ( q(ij-iip1,l,iq) -q(ij+iip1,l,iq) )
204       qs(ij) = CVMGP( xs(ij), q(ij,l,iq), pente(ij) )
205       qn(ij) = CVMGP( xn(ij), q(ij,l,iq), pente(ij) )
206       ENDDO
207c
208c
209c     ......    Calculs aux  poles   .............................
210c     ............................................................
211c
212c     On n'a pas besoin des valeurs de  qn au pole Nord ( j=    1) ,
213c      ainsi que de celles de           qs au pole Sud  ( j= jjp1)
214c
215c
216      airej2 = SSUM( iim, aire(iip2), 1 )
217      airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
218      DO i = 1, iim
219      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
220      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
221      ENDDO
222      qpns   = SSUM( iim,  airescb ,1 ) / airej2
223      qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
224c
225c     qpns , val.moyenne de q sur la ligne j= 2
226c     qpsn , val.moyenne de q sur la ligne j= jjm
227c
228c
229c
230c  on cherche si on a un extremum au pole
231c
232      extrpn=.true.
233      extrps=.true.
234      DO ij=2,iim
235         if((q(iip1+i,l,iq)-q(1,l,iq))*(q(iip2,l,iq)-q(1,l,iq)).lt.0.)
236     .   extrpn=.false.
237         if((q(ip1jm-iip1+i,l,iq)-q(1,l,iq))*
238     .     (q(ip1jm-iim,l,iq)-q(1,l,iq)).lt.0.)
239     .      extrps=.false.
240      ENDDO
241
242c   calcul des pentes au pole
243
244      if(extrpn) then
245         DO    ij     = 1, iip1
246            qs(ij)= q(ij,l,iq)
247         ENDDO
248      else
249         DO    ij     = 1, iip1
250            qs(ij)= q(ij,l,iq) + 0.5 * ( q(ij+ iip1,l,iq) - qpns )
251         ENDDO
252      endif
253
254      if(extrps) then
255         DO    ij     = 1, iip1
256            qn(ij+ip1jm) = q(ij+ip1jm,l,iq)
257         ENDDO
258      else
259         DO    ij     = 1, iip1
260            qn(ij+ip1jm) = q(ij+ip1jm,l,iq) +   0.5 *
261     *                    ( q(ij+ip1jm-iip1,l,iq) - qpsn )
262         ENDDO
263      endif
264
265c
266c
267c    .........................................................
268c    ......    Limitation des pentes au sud des boites   .....
269c    .........................................................   
270c
271      DO   ij     = 1, ip1jm
272      pente(ij)   = ( qs(ij) - q (ij+iip1,l,iq) )   *
273     *              ( q( ij,l,iq) - qs( ij ) )
274      qs(ij) = CVMGP( qs(ij) , q(ij+iip1,l,iq), pente(ij) )
275      qn(ij) = CVMGP( qn(ij) ,
276     *                q(ij,l,iq)+ q(ij,l,iq) -qs(ij), pente(ij) )
277      ENDDO
278c
279c
280c     .......................................................
281c     .... Limitation des pentes au nord des boites .........
282c     .......................................................
283c
284      DO   ij     = iip2, ip1jmp1
285      pente(ij)   = ( qn( ij  ) -  q(ij,l,iq) )  *
286     *              (  q(ij-iip1,l,iq) - qn(ij) )
287      qn(ij)      = CVMGP( qn(ij), q(ij-iip1,l,iq), pente(ij) )
288      qs(ij)      = CVMGP( qs(ij),
289     *                q(ij,l,iq)+ q(ij,l,iq) -qn(ij) , pente(ij) )
290      ENDDO
291c
292c
293c    .............................................................
294c    .....   Calculs des flux de  q  sur le plan horizontal ......
295c    .............................................................
296c
297c
298c
299c     
300c      ....  Selon   X  ....
301c
302         DO  ij     = iip2, ip1jm - 1
303c
304         qbxu( ij,l ) =  pbaru( ij,l )   *
305     *            CVMGT( qd(ij), qg(ij +1), pbaru(ij,l).GT.0. )
306         ENDDO
307c
308c     ..... correction  pour  qbxu(iip1,j,l)   .....
309c     ...   qbxu(iip1,j,l)= qbxu(1,j,l)  ...
310c
311c           &&&CDIR$ IVDEP
312         DO    ij     = iip1 +iip1, ip1jm, iip1
313         qbxu( ij,l ) = qbxu( ij - iim, l )
314         ENDDO
315c
316c      .... Selon Y  .....
317c
318         DO    ij     = 1, ip1jm
319         qbyv( ij,l ) =  pbarv( ij,l )    *
320     *             CVMGT( qn(ij+iip1), qs(ij), pbarv(ij,l).GT.0. )
321         ENDDO
322c
323c
324c
325  100  CONTINUE
326c
327c     ..........................................................
328c     ( ... fin des traitements en longitude et latitude  ... )
329c
330c
331c   stockage dans  dqh  de la convergence horiz.du flux d'humidite  .
332c                   ....
333c
334c
335                 CALL convflu( qbxu, qbyv, llm, dqh )
336c
337c
338c
339c    ----------------------------------------------------------------
340c    ----------------------------------------------------------------
341c    .............      Traitement  en  altitude        .............
342c    ----------------------------------------------------------------
343c    ----------------------------------------------------------------
344c
345c
346c     -------------
347c        q (llm)
348c
349c     -------------
350c     -------------
351c
352c        q(l+1)
353c
354c     -------------
355c             qh(l)
356c        q(l)
357c             qb(l)
358c     -------------
359c
360c        q(l-1)
361c
362c     -------------
363c     -------------
364c        q(1)
365c     -------------
366c
367c
368c    ... Calculs pour les niveaux 2 a llm -1  ...
369c
370c
371      DO  200   l = 2, llm -1
372
373       DO   ij     = 1, ip1jmp1
374       pente(ij)   = ( q(ij, l+1 ,iq)  - q(ij , l  ,   iq) )  *
375     *               ( q(ij,   l ,iq)  - q(ij ,l-1 ,   iq) )
376       xb(ij)      = q(ij,l,iq) - 0.25* ( q(ij,l+1,iq) - q(ij,l-1,iq) )
377       xh(ij)      = q(ij,l,iq) + 0.25* ( q(ij,l+1,iq) - q(ij,l-1,iq) )
378       qb(ij,l)    = CVMGP( xb(ij), q(ij,l,iq), pente(ij) )
379       qh(ij,l)    = CVMGP( xh(ij), q(ij,l,iq), pente(ij) )
380       ENDDO
381c
382c     ........................................................
383c     ......  Limitation des pentes en bas des boites   ......
384c     ........................................................
385c
386       DO   ij     = 1, ip1jmp1
387       pente(ij)   = ( qb(ij,l) - q ( ij,l-1,iq) )   *
388     *                ( q (ij,l,iq) - qb( ij,l ) )
389       qb(ij,l)    = CVMGP( qb(ij,l), q(ij,l-1,iq), pente(ij) )
390       qh(ij,l)    = CVMGP( qh(ij,l),
391     *               q(ij,l,iq) + q(ij,l,iq) -qb(ij,l), pente(ij) )
392       ENDDO
393c
394c
395c     ........................................................
396c     ......  Limitation des pentes en haut des boites   ......
397c     ........................................................
398c
399       DO   ij     = 1, ip1jmp1
400       pente(ij)   = ( qh(ij,l) - q ( ij,l+1,iq) )   *
401     *                ( q (ij,l,iq) - qh( ij,l ) )
402       qh(ij,l)    = CVMGP( qh(ij,l), q(ij,l+1,iq), pente(ij) )
403       qb(ij,l)    = CVMGP( qb(ij,l),
404     *               q(ij,l,iq) + q(ij,l,iq) -qh(ij,l), pente(ij) )
405       ENDDO
406c
407c
408  200  CONTINUE
409c
410c
411c    ............................................................
412c    .....     Calculs pour les niveaux l= 1 et l= llm  .........
413c    ............................................................
414c
415      DO   ij    = 1, ip1jmp1
416      qb(ij,1)   = q(ij, 1 , iq)
417      qb(ij,llm) = q(ij,llm, iq)
418      qh(ij,1)   = q(ij, 1 , iq)
419      qh(ij,llm) = q(ij,llm, iq)
420      ENDDO
421c
422
423c ---------------------------------------------------------------
424c   .... calcul des termes d'advection verticale  .......
425c ---------------------------------------------------------------
426
427c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dqh pour calculer dq
428c
429
430       DO 300 l = 1,llmm1
431c
432         DO  ij = 1,ip1jmp1
433          ww= - w( ij,l+1 ) *
434     *        CVMGT ( qh(ij,l), qb(ij,l+1), w(ij,l+1).LT.0.)
435
436          dq (ij, l ,iq ) = dqh(ij, l )   - dsig1( l ) * ww
437          dqh(ij,l+1    ) = dqh(ij,l+1)   + dsig1(l+1) * ww
438         ENDDO
439c
440  300   CONTINUE
441c
442c
443c
444        DO  ij = 1,ip1jmp1
445          dq( ij,llm,iq ) = dqh( ij,llm )
446        END DO
447c
448c
449      END IF
450c
451      RETURN
452      END
Note: See TracBrowser for help on using the repository browser.