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

Last change on this file since 773 was 135, checked in by aslmd, 14 years ago

CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET

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