My Project
 All Classes Files Functions Variables Macros
exner_hyb_p.F
Go to the documentation of this file.
1 !
2 ! $Id $
3 !
4  SUBROUTINE exner_hyb_p ( 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  IMPLICIT NONE
33 c
34 #include "dimensions.h"
35 #include "paramet.h"
36 #include "comconst.h"
37 #include "comgeom.h"
38 #include "comvert.h"
39 #include "serre.h"
40 
41  INTEGER ngrid
42  REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
43  REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
44 
45 c .... variables locales ...
46 
47  INTEGER l, ij
48  REAL unpl2k,dellta
49 
50  REAL ppn(iim),pps(iim)
51  REAL xpn, xps
52  REAL ssum
53  EXTERNAL ssum
54  INTEGER ije,ijb,jje,jjb
55  logical,save :: firstcall=.true.
56 !$OMP THREADPRIVATE(firstcall)
57  character(len=*),parameter :: modname="exner_hyb_p"
58 c
59 
60  ! Sanity check
61  if (firstcall) then
62  ! sanity checks for Shallow Water case (1 vertical layer)
63  if (llm.eq.1) then
64  if (kappa.ne.1) then
65  call abort_gcm(modname,
66  & "kappa!=1 , but running in Shallow Water mode!!",42)
67  endif
68  if (cpp.ne.r) then
69  call abort_gcm(modname,
70  & "cpp!=r , but running in Shallow Water mode!!",42)
71  endif
72  endif ! of if (llm.eq.1)
73 
74  firstcall=.false.
75  endif ! of if (firstcall)
76 
77 c$OMP BARRIER
78 
79 ! Specific behaviour for Shallow Water (1 vertical layer) case
80  if (llm.eq.1) then
81 
82  ! Compute pks(:),pk(:),pkf(:)
83  ijb=ij_begin
84  ije=ij_end
85 !$OMP DO SCHEDULE(STATIC)
86  DO ij=ijb, ije
87  pks(ij)=(cpp/preff)*ps(ij)
88  pk(ij,1) = .5*pks(ij)
89  pkf(ij,1)=pk(ij,1)
90  ENDDO
91 !$OMP ENDDO
92 
93 !$OMP MASTER
94  if (pole_nord) then
95  DO ij = 1, iim
96  ppn(ij) = aire( ij ) * pks( ij )
97  ENDDO
98  xpn = ssum(iim,ppn,1) /apoln
99 
100  DO ij = 1, iip1
101  pks( ij ) = xpn
102  pk(ij,1) = .5*pks(ij)
103  pkf(ij,1)=pk(ij,1)
104  ENDDO
105  endif
106 
107  if (pole_sud) then
108  DO ij = 1, iim
109  pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
110  ENDDO
111  xps = ssum(iim,pps,1) /apols
112 
113  DO ij = 1, iip1
114  pks( ij+ip1jm ) = xps
115  pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
116  pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
117  ENDDO
118  endif
119 !$OMP END MASTER
120 !$OMP BARRIER
121  jjb=jj_begin
122  jje=jj_end
123  CALL filtreg_p( pkf,jjb,jje, jmp1, llm, 2, 1, .true., 1 )
124 
125  ! our work is done, exit routine
126  return
127  endif ! of if (llm.eq.1)
128 
129 !!!! General case:
130 
131  unpl2k = 1.+ 2.* kappa
132 c
133  ijb=ij_begin
134  ije=ij_end
135 
136 c$OMP DO SCHEDULE(STATIC)
137  DO ij = ijb, ije
138  pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
139  ENDDO
140 c$OMP ENDDO
141 c Synchro OPENMP ici
142 
143 c$OMP MASTER
144  if (pole_nord) then
145  DO ij = 1, iim
146  ppn(ij) = aire( ij ) * pks( ij )
147  ENDDO
148  xpn = ssum(iim,ppn,1) /apoln
149 
150  DO ij = 1, iip1
151  pks( ij ) = xpn
152  ENDDO
153  endif
154 
155  if (pole_sud) then
156  DO ij = 1, iim
157  pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
158  ENDDO
159  xps = ssum(iim,pps,1) /apols
160 
161  DO ij = 1, iip1
162  pks( ij+ip1jm ) = xps
163  ENDDO
164  endif
165 c$OMP END MASTER
166 c$OMP BARRIER
167 c
168 c
169 c .... Calcul des coeff. alpha et beta pour la couche l = llm ..
170 c
171 c$OMP DO SCHEDULE(STATIC)
172  DO ij = ijb,ije
173  alpha(ij,llm) = 0.
174  beta(ij,llm) = 1./ unpl2k
175  ENDDO
176 c$OMP ENDDO NOWAIT
177 c
178 c ... Calcul des coeff. alpha et beta pour l = llm-1 a l = 2 ...
179 c
180  DO l = llm -1 , 2 , -1
181 c
182 c$OMP DO SCHEDULE(STATIC)
183  DO ij = ijb, ije
184  dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
185  alpha(ij,l) = - p(ij,l+1) / dellta * alpha(ij,l+1)
186  beta(ij,l) = p(ij,l ) / dellta
187  ENDDO
188 c$OMP ENDDO NOWAIT
189 c
190  ENDDO
191 
192 c
193 c ***********************************************************************
194 c ..... Calcul de pk pour la couche 1 , pres du sol ....
195 c
196 c$OMP DO SCHEDULE(STATIC)
197  DO ij = ijb, ije
198  pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) ) /
199  * ( p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
200  ENDDO
201 c$OMP ENDDO NOWAIT
202 c
203 c ..... Calcul de pk(ij,l) , pour l = 2 a l = llm ........
204 c
205  DO l = 2, llm
206 c$OMP DO SCHEDULE(STATIC)
207  DO ij = ijb, ije
208  pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
209  ENDDO
210 c$OMP ENDDO NOWAIT
211  ENDDO
212 c
213 c
214 c CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 )
215  DO l = 1, llm
216 c$OMP DO SCHEDULE(STATIC)
217  DO ij = ijb, ije
218  pkf(ij,l)=pk(ij,l)
219  ENDDO
220 c$OMP ENDDO NOWAIT
221  ENDDO
222 
223 c$OMP BARRIER
224 
225  jjb=jj_begin
226  jje=jj_end
227  CALL filtreg_p( pkf,jjb,jje, jmp1, llm, 2, 1, .true., 1 )
228 
229 
230  RETURN
231  END