My Project
 All Classes Files Functions Variables Macros
exner_hyb_loc.F
Go to the documentation of this file.
1 c
2 c $Id$
3 c
4  SUBROUTINE exner_hyb_loc(ngrid, ps, p,alpha,beta, pks,pk,pkf)
5 c
6 c Auteurs : P.Le Van , Fr. Hourdin .
7 c ..........
8 c
9 c .... ngrid, ps,p sont des argum.d'entree au sous-prog ...
10 c .... alpha,beta, pks,pk,pkf sont des argum.de sortie au sous-prog ...
11 c
12 c ************************************************************************
13 c Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des
14 c couches . Pk(l) sera calcule aux milieux des couches l ,entre les
15 c pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
16 c ************************************************************************
17 c .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont
18 c la pression et la fonction d'Exner au sol .
19 c
20 c -------- z
21 c A partir des relations ( 1 ) p*dz(pk) = kappa *pk*dz(p) et
22 c ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
23 c ( voir note de Fr.Hourdin ) ,
24 c
25 c on determine successivement , du haut vers le bas des couches, les
26 c coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
27 c puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,
28 c pk(ij,l) donne par la relation (2), pour l = 2 a l = llm .
29 c
30 c
31  USE parallel
32  USE mod_filtreg_p
33  USE write_field_loc
34  IMPLICIT NONE
35 c
36 #include "dimensions.h"
37 #include "paramet.h"
38 #include "comconst.h"
39 #include "comgeom.h"
40 #include "comvert.h"
41 #include "serre.h"
42 
43  INTEGER ngrid
44  REAL p(ijb_u:ije_u,llmp1),pk(ijb_u:ije_u,llm)
45  REAL pkf(ijb_u:ije_u,llm)
46  REAL ps(ijb_u:ije_u),pks(ijb_u:ije_u)
47  REAL alpha(ijb_u:ije_u,llm),beta(ijb_u:ije_u,llm)
48 
49 c .... variables locales ...
50 
51  INTEGER l, ij
52  REAL unpl2k,dellta
53 
54  REAL ppn(iim),pps(iim)
55  REAL xpn, xps
56  REAL ssum
57  EXTERNAL ssum
58  INTEGER ije,ijb,jje,jjb
59  logical,save :: firstcall=.true.
60 !$OMP THREADPRIVATE(firstcall)
61  character(len=*),parameter :: modname="exner_hyb_loc"
62 c
63 c$OMP BARRIER
64 
65  ! Sanity check
66  if (firstcall) then
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 
82 c$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 !$OMP BARRIER
126  jjb=jj_begin
127  jje=jj_end
128  CALL filtreg_p( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm,
129  & 2, 1, .true., 1 )
130 
131  ! our work is done, exit routine
132  return
133  endif ! of if (llm.eq.1)
134 
135 
136  unpl2k = 1.+ 2.* kappa
137 c
138  ijb=ij_begin
139  ije=ij_end
140 
141 c$OMP DO SCHEDULE(STATIC)
142  DO ij = ijb, ije
143  pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
144  ENDDO
145 c$OMP ENDDO
146 c Synchro OPENMP ici
147 
148 c$OMP MASTER
149  if (pole_nord) then
150  DO ij = 1, iim
151  ppn(ij) = aire( ij ) * pks( ij )
152  ENDDO
153  xpn = ssum(iim,ppn,1) /apoln
154 
155  DO ij = 1, iip1
156  pks( ij ) = xpn
157  ENDDO
158  endif
159 
160  if (pole_sud) then
161  DO ij = 1, iim
162  pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
163  ENDDO
164  xps = ssum(iim,pps,1) /apols
165 
166  DO ij = 1, iip1
167  pks( ij+ip1jm ) = xps
168  ENDDO
169  endif
170 c$OMP END MASTER
171 c$OMP BARRIER
172 c
173 c
174 c .... Calcul des coeff. alpha et beta pour la couche l = llm ..
175 c
176 c$OMP DO SCHEDULE(STATIC)
177  DO ij = ijb,ije
178  alpha(ij,llm) = 0.
179  beta(ij,llm) = 1./ unpl2k
180  ENDDO
181 c$OMP ENDDO NOWAIT
182 c
183 c ... Calcul des coeff. alpha et beta pour l = llm-1 a l = 2 ...
184 c
185  DO l = llm -1 , 2 , -1
186 c
187 c$OMP DO SCHEDULE(STATIC)
188  DO ij = ijb, ije
189  dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
190  alpha(ij,l) = - p(ij,l+1) / dellta * alpha(ij,l+1)
191  beta(ij,l) = p(ij,l ) / dellta
192  ENDDO
193 c$OMP ENDDO NOWAIT
194 c
195  ENDDO
196 
197 c
198 c ***********************************************************************
199 c ..... Calcul de pk pour la couche 1 , pres du sol ....
200 c
201 c$OMP DO SCHEDULE(STATIC)
202  DO ij = ijb, ije
203  pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) ) /
204  * ( p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
205  ENDDO
206 c$OMP ENDDO NOWAIT
207 c
208 c ..... Calcul de pk(ij,l) , pour l = 2 a l = llm ........
209 c
210  DO l = 2, llm
211 c$OMP DO SCHEDULE(STATIC)
212  DO ij = ijb, ije
213  pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
214  ENDDO
215 c$OMP ENDDO NOWAIT
216  ENDDO
217 c
218 c
219 c CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 )
220  DO l = 1, llm
221 c$OMP DO SCHEDULE(STATIC)
222  DO ij = ijb, ije
223  pkf(ij,l)=pk(ij,l)
224  ENDDO
225 c$OMP ENDDO NOWAIT
226  ENDDO
227 
228 c$OMP BARRIER
229 
230  jjb=jj_begin
231  jje=jj_end
232 #ifdef DEBUG_IO
233  call writefield_u('pkf',pkf)
234 #endif
235  CALL filtreg_p( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm,
236  & 2, 1, .true., 1 )
237 #ifdef DEBUG_IO
238  call writefield_u('pkf',pkf)
239 #endif
240 
241  RETURN
242  END