My Project
 All Classes Files Functions Variables Macros
vlspltgen_p.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4  SUBROUTINE vlspltgen_p( q,iadv,pente_max,masse,w,pbaru,pbarv,pdt,
5  , p,pk,teta )
6 c
7 c auteurs: p.le van, f.hourdin, f.forget, f.codron
8 c
9 c ********************************************************************
10 c shema d
11 
12 
13 
14 'advection " pseudo amont " .c + test sur humidite specifique: Q advecte< Qsat avalc (F. Codron, 10/99)c ********************************************************************c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
15 c
16 c pente_max facteur de limitation des pentes: 2 en general
17 c 0 pour un schema amont
18 c pbaru,pbarv,w flux de masse en u ,v ,w
19 c pdt pas de temps
20 c
21 c teta temperature potentielle, p pression aux interfaces,
22 c pk exner au milieu des couches necessaire pour calculer qsat
23 c --------------------------------------------------------------------
24  USE parallel
25  USE mod_hallo
26  USE write_field_p
27  USE vampir
28  USE infotrac, ONLY : nqtot
29  IMPLICIT NONE
30 
31 c
32 #include "dimensions.h"
33 #include "paramet.h"
34 #include "logic.h"
35 #include "comvert.h"
36 #include "comconst.h"
37 
38 c
39 c arguments:
40 c ----------
41  INTEGER iadv(nqtot)
42  REAL masse(ip1jmp1,llm),pente_max
43  REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
44  REAL q(ip1jmp1,llm,nqtot)
45  REAL w(ip1jmp1,llm),pdt
46  REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
47 c
48 c local
49 c ---------
50 c
51  INTEGER ij,l
52 c
53  REAL,SAVE :: qsat(ip1jmp1,llm)
54  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zm
55  REAL,SAVE :: mu(ip1jmp1,llm)
56  REAL,SAVE :: mv(ip1jm,llm)
57  REAL,SAVE :: mw(ip1jmp1,llm+1)
58  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zq
59  REAL zzpbar, zzw
60 
61  REAL qmin,qmax
62  DATA qmin,qmax/0.,1.e33/
63 
64 c--pour rapport de melange saturant--
65 
66  REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
67  REAL ptarg,pdelarg,foeew,zdelta
68  REAL tempe(ip1jmp1)
69  INTEGER ijb,ije,iq
70  LOGICAL, SAVE :: firstcall=.true.
71 !$OMP THREADPRIVATE(firstcall)
72  type(request) :: myrequest1
73  type(request) :: myrequest2
74 
75 c fonction psat(t)
76 
77  foeew( ptarg,pdelarg ) = exp(
78  * (r3les*(1.-pdelarg)+r3ies*pdelarg) * (ptarg-rtt)
79  * / (ptarg-(r4les*(1.-pdelarg)+r4ies*pdelarg)) )
80 
81  r2es = 380.11733
82  r3les = 17.269
83  r3ies = 21.875
84  r4les = 35.86
85  r4ies = 7.66
86  retv = 0.6077667
87  rtt = 273.16
88 
89 c Allocate variables depending on dynamic variable nqtot
90 
91  IF (firstcall) THEN
92  firstcall=.false.
93 !$OMP MASTER
94  ALLOCATE(zm(ip1jmp1,llm,nqtot))
95  ALLOCATE(zq(ip1jmp1,llm,nqtot))
96 !$OMP END MASTER
97 !$OMP BARRIER
98  END IF
99 c-- calcul de qsat en chaque point
100 c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
101 c pour eviter une exponentielle.
102 
103  call settag(myrequest1,100)
104  call settag(myrequest2,101)
105 
106 
107  ijb=ij_begin-iip1
108  ije=ij_end+iip1
109  if (pole_nord) ijb=ij_begin
110  if (pole_sud) ije=ij_end
111 
112 c$omp DO schedule(static,omp_chunk)
113  DO l = 1, llm
114  DO ij = ijb, ije
115  tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
116  ENDDO
117  DO ij = ijb, ije
118  zdelta = max( 0., sign(1., rtt - tempe(ij)) )
119  play = 0.5*(p(ij,l)+p(ij,l+1))
120  qsat(ij,l) = min(0.5, r2es* foeew(tempe(ij),zdelta) / play )
121  qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
122  ENDDO
123  ENDDO
124 c$omp END DO nowait
125 c print*,'Debut vlsplt version debug sans vlyqs'
126 
127  zzpbar = 0.5 * pdt
128  zzw = pdt
129 
130  ijb=ij_begin
131  ije=ij_end
132  if (pole_nord) ijb=ijb+iip1
133  if (pole_sud) ije=ije-iip1
134 
135 c$omp DO schedule(static,omp_chunk)
136  DO l=1,llm
137  DO ij = ijb,ije
138  mu(ij,l)=pbaru(ij,l) * zzpbar
139  ENDDO
140  ENDDO
141 c$omp END DO nowait
142 
143  ijb=ij_begin-iip1
144  ije=ij_end
145  if (pole_nord) ijb=ij_begin
146  if (pole_sud) ije=ij_end-iip1
147 
148 c$omp DO schedule(static,omp_chunk)
149  DO l=1,llm
150  DO ij=ijb,ije
151  mv(ij,l)=pbarv(ij,l) * zzpbar
152  ENDDO
153  ENDDO
154 c$omp END DO nowait
155 
156  ijb=ij_begin
157  ije=ij_end
158 
159 c$omp DO schedule(static,omp_chunk)
160  DO l=1,llm
161  DO ij=ijb,ije
162  mw(ij,l)=w(ij,l) * zzw
163  ENDDO
164  ENDDO
165 c$omp END DO nowait
166 
167 c$omp master
168  DO ij=ijb,ije
169  mw(ij,llm+1)=0.
170  ENDDO
171 c$omp END master
172 
173 c CALL scopy(ijp1llm,q,1,zq,1)
174 c CALL scopy(ijp1llm,masse,1,zm,1)
175 
176  ijb=ij_begin
177  ije=ij_end
178 
179  DO iq=1,nqtot
180 c$omp DO schedule(static,omp_chunk)
181  DO l=1,llm
182  zq(ijb:ije,l,iq)=q(ijb:ije,l,iq)
183  zm(ijb:ije,l,iq)=masse(ijb:ije,l)
184  ENDDO
185 c$omp END DO nowait
186  ENDDO
187 
188 
189 c$omp barrier
190  DO iq=1,nqtot
191 
192  if(iadv(iq) == 0) then
193 
194  cycle
195 
196  else if (iadv(iq)==10) then
197 
198 #ifdef _ADV_HALO
199  call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
200  & ij_begin,ij_begin+2*iip1-1)
201  call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
202  & ij_end-2*iip1+1,ij_end)
203 #else
204  call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
205  & ij_begin,ij_end)
206 #endif
207 
208 c$omp master
209  call vtb(vthallo)
210 c$omp END master
211  call register_hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,myrequest1)
212  call register_hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,myrequest1)
213 
214 c$omp master
215  call vte(vthallo)
216 c$omp END master
217  else if (iadv(iq)==14) then
218 
219 #ifdef _ADV_HALO
220  call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
221  & ij_begin,ij_begin+2*iip1-1)
222  call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
223  & ij_end-2*iip1+1,ij_end)
224 #else
225 
226  call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
227  & ij_begin,ij_end)
228 #endif
229 
230 c$omp master
231  call vtb(vthallo)
232 c$omp END master
233 
234  call register_hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,myrequest1)
235  call register_hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,myrequest1)
236 
237 c$omp master
238  call vte(vthallo)
239 c$omp END master
240  else
241 
242  stop 'vlspltgen_p : schema non parallelise'
243 
244  endif
245 
246  enddo
247 
248 
249 c$omp barrier
250 c$omp master
251  call vtb(vthallo)
252 c$omp END master
253 
254  call sendrequest(myrequest1)
255 
256 c$omp master
257  call vte(vthallo)
258 c$omp END master
259 c$omp barrier
260  do iq=1,nqtot
261 
262  if(iadv(iq) == 0) then
263 
264  cycle
265 
266  else if (iadv(iq)==10) then
267 
268 #ifdef _ADV_HALLO
269  call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
270  & ij_begin+2*iip1,ij_end-2*iip1)
271 #endif
272  else if (iadv(iq)==14) then
273 #ifdef _ADV_HALLO
274  call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
275  & ij_begin+2*iip1,ij_end-2*iip1)
276 #endif
277  else
278 
279  stop 'vlspltgen_p : schema non parallelise'
280 
281  endif
282 
283  enddo
284 c$omp barrier
285 c$omp master
286  call vtb(vthallo)
287 c$omp END master
288 
289 ! call WaitRecvRequest(MyRequest1)
290 ! call WaitSendRequest(MyRequest1)
291 c$omp barrier
292  call waitrequest(myrequest1)
293 
294 
295 c$omp master
296  call vte(vthallo)
297 c$omp END master
298 c$omp barrier
299 
300  do iq=1,nqtot
301 
302  if(iadv(iq) == 0) then
303 
304  cycle
305 
306  else if (iadv(iq)==10) then
307 
308  call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
309 
310  else if (iadv(iq)==14) then
311 
312  call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
313 
314  else
315 
316  stop 'vlspltgen_p : schema non parallelise'
317 
318  endif
319 
320  enddo
321 
322 
323  do iq=1,nqtot
324 
325  if(iadv(iq) == 0) then
326 
327  cycle
328 
329  else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
330 
331 c$omp barrier
332 #ifdef _ADV_HALLO
333  call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
334  & ij_begin,ij_begin+2*iip1-1)
335  call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
336  & ij_end-2*iip1+1,ij_end)
337 #else
338  call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
339  & ij_begin,ij_end)
340 #endif
341 c$omp barrier
342 
343 c$omp master
344  call vtb(vthallo)
345 c$omp END master
346 
347  call register_hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,myrequest2)
348  call register_hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,myrequest2)
349 
350 c$omp master
351  call vte(vthallo)
352 c$omp END master
353 c$omp barrier
354  else
355 
356  stop 'vlspltgen_p : schema non parallelise'
357 
358  endif
359 
360  enddo
361 c$omp barrier
362 
363 c$omp master
364  call vtb(vthallo)
365 c$omp END master
366 
367  call sendrequest(myrequest2)
368 
369 c$omp master
370  call vte(vthallo)
371 c$omp END master
372 
373 c$omp barrier
374  do iq=1,nqtot
375 
376  if(iadv(iq) == 0) then
377 
378  cycle
379 
380  else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
381 c$omp barrier
382 
383 #ifdef _ADV_HALLO
384  call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
385  & ij_begin+2*iip1,ij_end-2*iip1)
386 #endif
387 
388 c$omp barrier
389  else
390 
391  stop 'vlspltgen_p : schema non parallelise'
392 
393  endif
394 
395  enddo
396 
397 c$omp barrier
398 c$omp master
399  call vtb(vthallo)
400 c$omp END master
401 
402 ! call WaitRecvRequest(MyRequest2)
403 ! call WaitSendRequest(MyRequest2)
404 c$omp barrier
405  CALL waitrequest(myrequest2)
406 
407 c$omp master
408  call vte(vthallo)
409 c$omp END master
410 c$omp barrier
411 
412 
413  do iq=1,nqtot
414 
415  if(iadv(iq) == 0) then
416 
417  cycle
418 
419  else if (iadv(iq)==10) then
420 
421  call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
422 
423  else if (iadv(iq)==14) then
424 
425  call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
426 
427  else
428 
429  stop 'vlspltgen_p : schema non parallelise'
430 
431  endif
432 
433  enddo
434 
435  do iq=1,nqtot
436 
437  if(iadv(iq) == 0) then
438 
439  cycle
440 
441  else if (iadv(iq)==10) then
442 
443  call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
444  & ij_begin,ij_end)
445 
446  else if (iadv(iq)==14) then
447 
448  call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
449  & ij_begin,ij_end)
450 
451  else
452 
453  stop 'vlspltgen_p : schema non parallelise'
454 
455  endif
456 
457  enddo
458 
459 
460  ijb=ij_begin
461  ije=ij_end
462 c$omp barrier
463 
464 
465  DO iq=1,nqtot
466 
467 c$omp DO schedule(static,omp_chunk)
468  DO l=1,llm
469  DO ij=ijb,ije
470 c print *,'zq-->',ij,l,iq,zq(ij,l,iq)
471 c print *,'q-->',ij,l,iq,q(ij,l,iq)
472  q(ij,l,iq)=zq(ij,l,iq)
473  ENDDO
474  ENDDO
475 c$omp END DO nowait
476 
477 c$omp DO schedule(static,omp_chunk)
478  DO l=1,llm
479  DO ij=ijb,ije-iip1+1,iip1
480  q(ij+iim,l,iq)=q(ij,l,iq)
481  ENDDO
482  ENDDO
483 c$omp END DO nowait
484 
485  ENDDO
486 
487 
488 c$omp barrier
489 
490 cc$omp master
491 c call waitsendrequest(myrequest1)
492 c call waitsendrequest(myrequest2)
493 cc$omp END master
494 cc$omp barrier
495 
496  RETURN
497  END