source: LMDZ4/trunk/libf/dyn3d/exner_hyb.F @ 701

Last change on this file since 701 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.3 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE  exner_hyb ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
5c
6c     Auteurs :  P.Le Van  , Fr. Hourdin  .
7c    ..........
8c
9c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
10c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
11c
12c   ************************************************************************
13c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des
14c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
15c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
16c   ************************************************************************
17c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
18c    la pression et la fonction d'Exner  au  sol  .
19c
20c                                 -------- z                                   
21c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
22c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
23c    ( voir note de Fr.Hourdin )  ,
24c
25c    on determine successivement , du haut vers le bas des couches, les
26c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
27c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
28c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
29c
30c
31      IMPLICIT NONE
32c
33#include "dimensions.h"
34#include "paramet.h"
35#include "comconst.h"
36#include "comgeom.h"
37#include "comvert.h"
38#include "serre.h"
39
40      INTEGER  ngrid
41      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
42      REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
43
44c    .... variables locales   ...
45
46      INTEGER l, ij
47      REAL unpl2k,dellta
48
49      REAL ppn(iim),pps(iim)
50      REAL xpn, xps
51      REAL SSUM
52c
53     
54      unpl2k    = 1.+ 2.* kappa
55c
56      DO   ij  = 1, ngrid
57        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
58      ENDDO
59
60      DO  ij   = 1, iim
61        ppn(ij) = aire(   ij   ) * pks(  ij     )
62        pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
63      ENDDO
64      xpn      = SSUM(iim,ppn,1) /apoln
65      xps      = SSUM(iim,pps,1) /apols
66
67      DO ij   = 1, iip1
68        pks(   ij     )  =  xpn
69        pks( ij+ip1jm )  =  xps
70      ENDDO
71c
72c
73c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
74c
75      DO     ij      = 1, ngrid
76       alpha(ij,llm) = 0.
77       beta (ij,llm) = 1./ unpl2k
78      ENDDO
79c
80c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
81c
82      DO l = llm -1 , 2 , -1
83c
84        DO ij = 1, ngrid
85        dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
86        alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
87        beta (ij,l)  =   p(ij,l  ) / dellta   
88        ENDDO
89c
90      ENDDO
91c
92c  ***********************************************************************
93c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
94c
95      DO   ij   = 1, ngrid
96       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
97     *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
98      ENDDO
99c
100c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
101c
102      DO l = 2, llm
103        DO   ij   = 1, ngrid
104         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
105        ENDDO
106      ENDDO
107c
108c
109      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
110      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
111     
112
113      RETURN
114      END
Note: See TracBrowser for help on using the repository browser.