source: LMDZ5/trunk/libf/dyn3dmem/exner_milieu_loc.F @ 1632

Last change on this file since 1632 was 1632, checked in by Laurent Fairhead, 12 years ago

Import initial du répertoire dyn3dmem

Attention! ceci n'est qu'une version préliminaire du code "basse mémoire":
le code contenu dans ce répertoire est basé sur la r1320 et a donc besoin
d'être mis à jour par rapport à la dynamique parallèle d'aujourd'hui.
Ce code est toutefois mis à disposition pour circonvenir à des problèmes
de mémoire que certaines configurations du modèle pourraient rencontrer.
Dans l'état, il compile et tourne sur vargas et au CCRT


Initial import of dyn3dmem

Warning! this is just a preliminary version of the memory light code:
it is based on r1320 of the code and thus needs to be updated before
it can replace the present dyn3dpar code. It is nevertheless put at your
disposal to circumvent some memory problems some LMDZ configurations may
encounter. In its present state, it will compile and run on vargas and CCRT

File size: 5.9 KB
Line 
1!
2! $Id $
3!
4      SUBROUTINE  exner_milieu_loc ( ngrid, ps, p,beta, pks, pk, pkf )
5c
6c     Auteurs :  F. Forget , Y. Wanherdrick
7c P.Le Van  , Fr. Hourdin  .
8c    ..........
9c
10c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
11c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
12c
13c   ************************************************************************
14c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
15c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
16c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
17c   ************************************************************************
18c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
19c    la pression et la fonction d'Exner  au  sol  .
20c
21c     WARNING : CECI est une version speciale de exner_hyb originale
22c               Utilise dans la version martienne pour pouvoir
23c               tourner avec des coordonnees verticales complexe
24c              => Il ne verifie PAS la condition la proportionalite en
25c              energie totale/ interne / potentielle (F.Forget 2001)
26c    ( voir note de Fr.Hourdin )  ,
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(ijb_u:ije_u,llmp1),pk(ijb_u:ije_u,llm)
40      REAL pkf(ijb_u:ije_u,llm)
41      REAL ps(ijb_u:ije_u),pks(ijb_u:ije_u)
42      REAL alpha(ijb_u:ije_u,llm),beta(ijb_u:ije_u,llm)
43
44c    .... variables locales   ...
45
46      INTEGER l, ij
47      REAL dum1
48
49      REAL ppn(iim),pps(iim)
50      REAL xpn, xps
51      REAL SSUM
52      EXTERNAL SSUM
53      INTEGER ije,ijb,jje,jjb
54      logical,save :: firstcall=.true.
55!$OMP THREADPRIVATE(firstcall)
56      character(len=*),parameter :: modname="exner_milieu_p"
57
58      ! Sanity check
59      if (firstcall) then
60        ! check that vertical discretization is compatible
61        ! with this routine
62        if (disvert_type.ne.2) then
63          call abort_gcm(modname,
64     &     "this routine should only be called if disvert_type==2",42)
65        endif
66       
67        ! sanity checks for Shallow Water case (1 vertical layer)
68        if (llm.eq.1) then
69          if (kappa.ne.1) then
70            call abort_gcm(modname,
71     &      "kappa!=1 , but running in Shallow Water mode!!",42)
72          endif
73          if (cpp.ne.r) then
74            call abort_gcm(modname,
75     &      "cpp!=r , but running in Shallow Water mode!!",42)
76          endif
77        endif ! of if (llm.eq.1)
78
79        firstcall=.false.
80      endif ! of if (firstcall)
81     
82c$OMP BARRIER
83
84! Specific behaviour for Shallow Water (1 vertical layer) case
85      if (llm.eq.1) then
86             
87        ! Compute pks(:),pk(:),pkf(:)
88        ijb=ij_begin
89        ije=ij_end
90!$OMP DO SCHEDULE(STATIC)
91        DO ij=ijb, ije
92          pks(ij)=(cpp/preff)*ps(ij)
93          pk(ij,1) = .5*pks(ij)
94          pkf(ij,1)=pk(ij,1)
95        ENDDO
96!$OMP ENDDO
97
98!$OMP MASTER
99      if (pole_nord) then
100        DO  ij   = 1, iim
101          ppn(ij) = aire(   ij   ) * pks(  ij     )
102        ENDDO
103        xpn      = SSUM(iim,ppn,1) /apoln
104 
105        DO ij   = 1, iip1
106          pks(   ij     )  =  xpn
107          pk(ij,1) = .5*pks(ij)
108          pkf(ij,1)=pk(ij,1)
109        ENDDO
110      endif
111     
112      if (pole_sud) then
113        DO  ij   = 1, iim
114          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
115        ENDDO
116        xps      = SSUM(iim,pps,1) /apols
117 
118        DO ij   = 1, iip1
119          pks( ij+ip1jm )  =  xps
120          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
121          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
122        ENDDO
123      endif
124!$OMP END MASTER
125
126        jjb=jj_begin
127        jje=jj_end
128        CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
129
130        ! our work is done, exit routine
131        return
132      endif ! of if (llm.eq.1)
133
134!!!! General case:
135
136c     -------------
137c     Calcul de pks
138c     -------------
139   
140      ijb=ij_begin
141      ije=ij_end
142
143c$OMP DO SCHEDULE(STATIC)
144      DO   ij  = ijb, ije
145        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
146      ENDDO
147c$OMP ENDDO
148c Synchro OPENMP ici
149
150c$OMP MASTER
151      if (pole_nord) then
152        DO  ij   = 1, iim
153          ppn(ij) = aire(   ij   ) * pks(  ij     )
154        ENDDO
155        xpn      = SSUM(iim,ppn,1) /apoln
156 
157        DO ij   = 1, iip1
158          pks(   ij     )  =  xpn
159        ENDDO
160      endif
161     
162      if (pole_sud) then
163        DO  ij   = 1, iim
164          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
165        ENDDO
166        xps      = SSUM(iim,pps,1) /apols
167 
168        DO ij   = 1, iip1
169          pks( ij+ip1jm )  =  xps
170        ENDDO
171      endif
172c$OMP END MASTER
173c
174c
175c    .... Calcul de pk  pour la couche l
176c    --------------------------------------------
177c
178      dum1 = cpp * (2*preff)**(-kappa)
179      DO l = 1, llm-1
180c$OMP DO SCHEDULE(STATIC)
181        DO   ij   = ijb, ije
182         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
183        ENDDO
184c$OMP ENDDO NOWAIT       
185      ENDDO
186
187c    .... Calcul de pk  pour la couche l = llm ..
188c    (on met la meme distance (en log pression)  entre Pk(llm)
189c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
190
191c$OMP DO SCHEDULE(STATIC)
192      DO   ij   = ijb, ije
193         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
194      ENDDO
195c$OMP ENDDO NOWAIT       
196
197
198c    calcul de pkf
199c    -------------
200c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
201      DO l = 1, llm
202c$OMP DO SCHEDULE(STATIC)
203         DO   ij   = ijb, ije
204           pkf(ij,l)=pk(ij,l)
205         ENDDO
206c$OMP ENDDO NOWAIT             
207      ENDDO
208
209c$OMP BARRIER
210     
211      jjb=jj_begin
212      jje=jj_end
213      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
214     
215c    EST-CE UTILE ?? : calcul de beta
216c    --------------------------------
217      DO l = 2, llm
218c$OMP DO SCHEDULE(STATIC)
219        DO   ij   = ijb, ije
220          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
221        ENDDO
222c$OMP ENDDO NOWAIT             
223      ENDDO
224
225      RETURN
226      END
Note: See TracBrowser for help on using the repository browser.