source: LMDZ4/branches/LMDZ4-dev-20091210/libf/dyn3dpar/exner_hyb_p.F @ 1959

Last change on this file since 1959 was 985, checked in by Laurent Fairhead, 16 years ago

Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.1 KB
Line 
1      SUBROUTINE  exner_hyb_p ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
2c
3c     Auteurs :  P.Le Van  , Fr. Hourdin  .
4c    ..........
5c
6c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
7c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
8c
9c   ************************************************************************
10c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des
11c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
12c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
13c   ************************************************************************
14c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
15c    la pression et la fonction d'Exner  au  sol  .
16c
17c                                 -------- z                                   
18c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
19c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
20c    ( voir note de Fr.Hourdin )  ,
21c
22c    on determine successivement , du haut vers le bas des couches, les
23c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
24c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
25c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
26c
27c
28      USE parallel
29      IMPLICIT NONE
30c
31#include "dimensions.h"
32#include "paramet.h"
33#include "comconst.h"
34#include "comgeom.h"
35#include "comvert.h"
36#include "serre.h"
37
38      INTEGER  ngrid
39      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
40      REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
41
42c    .... variables locales   ...
43
44      INTEGER l, ij
45      REAL unpl2k,dellta
46
47      REAL ppn(iim),pps(iim)
48      REAL xpn, xps
49      REAL SSUM
50      EXTERNAL SSUM
51      INTEGER ije,ijb,jje,jjb
52c
53c$OMP BARRIER           
54      unpl2k    = 1.+ 2.* kappa
55c
56      ijb=ij_begin
57      ije=ij_end
58
59c$OMP DO SCHEDULE(STATIC)
60      DO   ij  = ijb, ije
61        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
62      ENDDO
63c$OMP ENDDO
64c Synchro OPENMP ici
65
66c$OMP MASTER
67      if (pole_nord) then
68        DO  ij   = 1, iim
69          ppn(ij) = aire(   ij   ) * pks(  ij     )
70        ENDDO
71        xpn      = SSUM(iim,ppn,1) /apoln
72 
73        DO ij   = 1, iip1
74          pks(   ij     )  =  xpn
75        ENDDO
76      endif
77     
78      if (pole_sud) then
79        DO  ij   = 1, iim
80          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
81        ENDDO
82        xps      = SSUM(iim,pps,1) /apols
83 
84        DO ij   = 1, iip1
85          pks( ij+ip1jm )  =  xps
86        ENDDO
87      endif
88c$OMP END MASTER
89c
90c
91c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
92c
93c$OMP DO SCHEDULE(STATIC)
94      DO     ij      = ijb,ije
95       alpha(ij,llm) = 0.
96       beta (ij,llm) = 1./ unpl2k
97      ENDDO
98c$OMP ENDDO NOWAIT
99c
100c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
101c
102      DO l = llm -1 , 2 , -1
103c
104c$OMP DO SCHEDULE(STATIC)
105        DO ij = ijb, ije
106        dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
107        alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
108        beta (ij,l)  =   p(ij,l  ) / dellta   
109        ENDDO
110c$OMP ENDDO NOWAIT
111c
112      ENDDO
113
114c
115c  ***********************************************************************
116c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
117c
118c$OMP DO SCHEDULE(STATIC)
119      DO   ij   = ijb, ije
120       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
121     *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
122      ENDDO
123c$OMP ENDDO NOWAIT
124c
125c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
126c
127      DO l = 2, llm
128c$OMP DO SCHEDULE(STATIC)
129        DO   ij   = ijb, ije
130         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
131        ENDDO
132c$OMP ENDDO NOWAIT       
133      ENDDO
134c
135c
136c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
137      DO l = 1, llm
138c$OMP DO SCHEDULE(STATIC)
139         DO   ij   = ijb, ije
140           pkf(ij,l)=pk(ij,l)
141         ENDDO
142c$OMP ENDDO NOWAIT             
143      ENDDO
144
145c$OMP BARRIER
146     
147      jjb=jj_begin
148      jje=jj_end
149      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
150     
151
152      RETURN
153      END
Note: See TracBrowser for help on using the repository browser.