source: LMDZ4/branches/LMDZ4_AR5/libf/dyn3dpar/exner_hyb_p.F @ 4056

Last change on this file since 4056 was 1558, checked in by yann meurdesoif, 13 years ago

Ajout d'une synchronisation OpenMP. Devrait corriger les problemes de reproductibilites en OpenMP.

YM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.1 KB
RevLine 
[630]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
[764]50      EXTERNAL SSUM
[630]51      INTEGER ije,ijb,jje,jjb
52c
[985]53c$OMP BARRIER           
[630]54      unpl2k    = 1.+ 2.* kappa
55c
56      ijb=ij_begin
57      ije=ij_end
[774]58
[985]59c$OMP DO SCHEDULE(STATIC)
[630]60      DO   ij  = ijb, ije
61        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
62      ENDDO
[985]63c$OMP ENDDO
64c Synchro OPENMP ici
[630]65
[985]66c$OMP MASTER
[630]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
[985]88c$OMP END MASTER
[1558]89c$OMP BARRIER
[630]90c
91c
92c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
93c
[985]94c$OMP DO SCHEDULE(STATIC)
[630]95      DO     ij      = ijb,ije
96       alpha(ij,llm) = 0.
97       beta (ij,llm) = 1./ unpl2k
98      ENDDO
[985]99c$OMP ENDDO NOWAIT
[630]100c
101c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
102c
103      DO l = llm -1 , 2 , -1
104c
[985]105c$OMP DO SCHEDULE(STATIC)
[630]106        DO ij = ijb, ije
107        dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
108        alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
109        beta (ij,l)  =   p(ij,l  ) / dellta   
110        ENDDO
[985]111c$OMP ENDDO NOWAIT
[630]112c
113      ENDDO
[774]114
[630]115c
116c  ***********************************************************************
117c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
118c
[985]119c$OMP DO SCHEDULE(STATIC)
[630]120      DO   ij   = ijb, ije
121       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
122     *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
123      ENDDO
[985]124c$OMP ENDDO NOWAIT
[630]125c
126c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
127c
128      DO l = 2, llm
[985]129c$OMP DO SCHEDULE(STATIC)
[630]130        DO   ij   = ijb, ije
131         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
132        ENDDO
[985]133c$OMP ENDDO NOWAIT       
[630]134      ENDDO
135c
136c
137c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
[985]138      DO l = 1, llm
139c$OMP DO SCHEDULE(STATIC)
140         DO   ij   = ijb, ije
141           pkf(ij,l)=pk(ij,l)
142         ENDDO
143c$OMP ENDDO NOWAIT             
144      ENDDO
145
[774]146c$OMP BARRIER
[630]147     
148      jjb=jj_begin
149      jje=jj_end
150      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
151     
152
153      RETURN
154      END
Note: See TracBrowser for help on using the repository browser.