source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/vlspltgen_loc.F90 @ 5116

Last change on this file since 5116 was 5116, checked in by abarral, 2 months ago

rename modules properly lmdz_*
move ismin, ismax, minmax into new lmdz_libmath.f90
(lint) uppercase fortran keywords
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 11.4 KB
Line 
1! $Header$
2
3SUBROUTINE vlspltgen_loc(q, pente_max, masse, w, pbaru, pbarv, &
4        pdt, p, pk, teta)
5
6  !
7  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron
8  !
9  !    ********************************************************************
10  !      Schema  d'advection " pseudo amont " .
11  !  + test sur humidite specifique: Q advecte< Qsat aval
12  !               (F. Codron, 10/99)
13  !    ********************************************************************
14  ! q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
15  !
16  ! pente_max facteur de limitation des pentes: 2 en general
17  !                                            0 pour un schema amont
18  ! pbaru,pbarv,w flux de masse en u ,v ,w
19  ! pdt pas de temps
20  !
21  ! teta temperature potentielle, p pression aux interfaces,
22  ! pk exner au milieu des couches necessaire pour calculer Qsat
23  !   --------------------------------------------------------------------
24  USE parallel_lmdz
25  USE mod_hallo
26  USE write_field_loc, ONLY: WriteField_u, WriteField_v
27  USE VAMPIR
28  ! CRisi: on rajoute variables utiles d'infotrac
29  USE infotrac, ONLY: nqtot, tracers, isoCheck
30  USE vlspltgen_mod
31  USE comconst_mod, ONLY: cpp
32  USE logic_mod, ONLY: adv_qsat_liq
33  IMPLICIT NONE
34
35  !
36  include "dimensions.h"
37  include "paramet.h"
38
39  !
40  !   Arguments:
41  !   ----------
42  REAL :: masse(ijb_u:ije_u, llm), pente_max
43  REAL :: pbaru(ijb_u:ije_u, llm), pbarv(ijb_v:ije_v, llm)
44  REAL :: q(ijb_u:ije_u, llm, nqtot)
45  REAL :: w(ijb_u:ije_u, llm), pdt
46  REAL :: p(ijb_u:ije_u, llmp1), teta(ijb_u:ije_u, llm)
47  REAL :: pk(ijb_u:ije_u, llm)
48  !
49  !  Local
50  !   ---------
51  !
52  INTEGER :: ij, l
53  !
54  REAL :: zzpbar, zzw
55
56  REAL :: qmin, qmax
57  DATA qmin, qmax/0., 1.e33/
58
59  !--pour rapport de melange saturant--
60
61  REAL :: rtt, retv, r2es, r3les, r3ies, r4les, r4ies, play
62  REAL :: ptarg, pdelarg, foeew, zdelta
63  REAL :: tempe(ijb_u:ije_u)
64  INTEGER :: ijb, ije, iq, iq2, ifils
65  LOGICAL, SAVE :: firstcall = .TRUE.
66  !$OMP THREADPRIVATE(firstcall)
67  type(request), SAVE :: MyRequest1
68  !$OMP THREADPRIVATE(MyRequest1)
69  type(request), SAVE :: MyRequest2
70  !$OMP THREADPRIVATE(MyRequest2)
71  !    fonction psat(T)
72
73  FOEEW (PTARG, PDELARG) = EXP (&
74          (R3LES * (1. - PDELARG) + R3IES * PDELARG) * (PTARG - RTT) &
75                  / (PTARG - (R4LES * (1. - PDELARG) + R4IES * PDELARG)))
76
77  r2es = 380.11733
78  r3les = 17.269
79  r3ies = 21.875
80  r4les = 35.86
81  r4ies = 7.66
82  retv = 0.6077667
83  rtt = 273.16
84
85  ! Allocate variables depending on dynamic variable nqtot
86
87  IF (firstcall) THEN
88    firstcall = .FALSE.
89  END IF
90  !-- Calcul de Qsat en chaque point
91  !-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
92  !   pour eviter une exponentielle.
93
94  CALL SetTag(MyRequest1, 100)
95  CALL SetTag(MyRequest2, 101)
96
97  ijb = ij_begin - iip1
98  ije = ij_end + iip1
99  if (pole_nord) ijb = ij_begin
100  if (pole_sud) ije = ij_end
101
102  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
103  DO l = 1, llm
104    DO ij = ijb, ije
105      tempe(ij) = teta(ij, l) * pk(ij, l) / cpp
106    ENDDO
107    DO ij = ijb, ije
108      IF (adv_qsat_liq) THEN
109        zdelta = 0.
110      ELSE
111        zdelta = MAX(0., SIGN(1., rtt - tempe(ij)))
112      ENDIF
113      play = 0.5 * (p(ij, l) + p(ij, l + 1))
114      qsat(ij, l) = MIN(0.5, r2es * FOEEW(tempe(ij), zdelta) / play)
115      qsat(ij, l) = qsat(ij, l) / (1. - retv * qsat(ij, l))
116    ENDDO
117  ENDDO
118  !$OMP END DO NOWAIT
119  ! PRINT*,'Debut vlsplt version debug sans vlyqs'
120
121  zzpbar = 0.5 * pdt
122  zzw = pdt
123
124  ijb = ij_begin
125  ije = ij_end
126  if (pole_nord) ijb = ijb + iip1
127  if (pole_sud)  ije = ije - iip1
128
129  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
130  DO l = 1, llm
131    DO ij = ijb, ije
132      mu(ij, l) = pbaru(ij, l) * zzpbar
133    ENDDO
134  ENDDO
135  !$OMP END DO NOWAIT
136
137  ijb = ij_begin - iip1
138  ije = ij_end
139  if (pole_nord) ijb = ij_begin
140  if (pole_sud)  ije = ij_end - iip1
141
142  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
143  DO l = 1, llm
144    DO ij = ijb, ije
145      mv(ij, l) = pbarv(ij, l) * zzpbar
146    ENDDO
147  ENDDO
148  !$OMP END DO NOWAIT
149
150  ijb = ij_begin
151  ije = ij_end
152
153  DO iq = 1, nqtot
154    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
155    DO l = 1, llm
156      DO ij = ijb, ije
157        mw(ij, l, iq) = w(ij, l) * zzw
158      ENDDO
159    ENDDO
160    !$OMP END DO NOWAIT
161  ENDDO
162
163  DO iq = 1, nqtot
164    !$OMP MASTER
165    DO ij = ijb, ije
166      mw(ij, llm + 1, iq) = 0.
167    ENDDO
168    !$OMP END MASTER
169  ENDDO
170
171  ijb = ij_begin
172  ije = ij_end
173
174  DO iq = 1, nqtot
175    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
176    DO l = 1, llm
177      zq(ijb:ije, l, iq) = q(ijb:ije, l, iq)
178      zm(ijb:ije, l, iq) = masse(ijb:ije, l)
179    ENDDO
180    !$OMP END DO NOWAIT
181  ENDDO
182
183  ! verif temporaire
184  ijb = ij_begin
185  ije = ij_end
186  CALL check_isotopes(zq, ijb, ije, 'vlspltgen_loc 191')
187
188  !$OMP BARRIER
189  DO iq = 1, nqtot
190    ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
191    IF(tracers(iq)%parent /= 'air') CYCLE
192    !WRITE(*,*) 'vlspltgen 192: iq,iadv=',iq,tracers(iq)%iadv
193    SELECT CASE(tracers(iq)%iadv)
194    CASE(0); CYCLE
195    CASE(10)
196      CALL vlx_loc(zq, pente_max, zm, mu, &
197              ij_begin, ij_end, iq)
198
199      !$OMP MASTER
200      CALL VTb(VTHallo)
201      !$OMP END MASTER
202      CALL Register_Hallo_u(zq(:, :, iq), llm, 2, 2, 2, 2, MyRequest1)
203      CALL Register_Hallo_u(zm(:, :, iq), llm, 1, 1, 1, 1, MyRequest1)
204      ! CRisi
205      do ifils = 1, tracers(iq)%nqDescen
206        iq2 = tracers(iq)%iqDescen(ifils)
207        CALL Register_Hallo_u(zq(:, :, iq2), llm, 2, 2, 2, 2, MyRequest1)
208        CALL Register_Hallo_u(zm(:, :, iq2), llm, 1, 1, 1, 1, MyRequest1)
209      enddo
210
211      !$OMP MASTER
212      CALL VTe(VTHallo)
213      !$OMP END MASTER
214    CASE(14)
215      CALL vlxqs_loc(zq, pente_max, zm, mu, &
216              qsat, ij_begin, ij_end, iq)
217
218      !$OMP MASTER
219      CALL VTb(VTHallo)
220      !$OMP END MASTER
221
222      CALL Register_Hallo_u(zq(:, :, iq), llm, 2, 2, 2, 2, MyRequest1)
223      CALL Register_Hallo_u(zm(:, :, iq), llm, 1, 1, 1, 1, MyRequest1)
224      do ifils = 1, tracers(iq)%nqDescen
225        iq2 = tracers(iq)%iqDescen(ifils)
226        CALL Register_Hallo_u(zq(:, :, iq2), llm, 2, 2, 2, 2, MyRequest1)
227        CALL Register_Hallo_u(zm(:, :, iq2), llm, 1, 1, 1, 1, MyRequest1)
228      enddo
229
230      !$OMP MASTER
231      CALL VTe(VTHallo)
232      !$OMP END MASTER
233    CASE DEFAULT
234      CALL abort_gcm("vlspltgen_loc", "schema non parallelise", 1)
235    END SELECT
236
237  enddo !DO iq=1,nqtot
238
239
240  !$OMP BARRIER
241  !$OMP MASTER
242  CALL VTb(VTHallo)
243  !$OMP END MASTER
244
245  CALL SendRequest(MyRequest1)
246
247  !$OMP MASTER
248  CALL VTe(VTHallo)
249  !$OMP END MASTER
250  !$OMP BARRIER
251
252  ! verif temporaire
253  ijb = ij_begin - 2 * iip1
254  ije = ij_end + 2 * iip1
255  if (pole_nord) ijb = ij_begin
256  if (pole_sud)  ije = ij_end
257  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 280')
258
259  do iq = 1, nqtot
260    IF(tracers(iq)%parent /= 'air') CYCLE
261    !WRITE(*,*) 'vlspltgen 279: iq=',iq
262
263    SELECT CASE(tracers(iq)%iadv)
264    CASE(0); CYCLE
265    CASE(10)
266    CASE(14)
267    CASE DEFAULT
268      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
269    END SELECT
270
271  enddo
272  !$OMP BARRIER
273  !$OMP MASTER
274  CALL VTb(VTHallo)
275  !$OMP END MASTER
276
277  ! CALL WaitRecvRequest(MyRequest1)
278  ! CALL WaitSendRequest(MyRequest1)
279  !$OMP BARRIER
280  CALL WaitRequest(MyRequest1)
281
282
283  !$OMP MASTER
284  CALL VTe(VTHallo)
285  !$OMP END MASTER
286  !$OMP BARRIER
287
288  IF(isoCheck) THEN
289    CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 326')
290    ijb = ij_begin - 2 * iip1
291    ije = ij_end + 2 * iip1
292    if (pole_nord) ijb = ij_begin
293    if (pole_sud)  ije = ij_end
294    CALL check_isotopes(zq, ijb, ije, 'vlspltgen_loc 336')
295  END IF
296
297  do iq = 1, nqtot
298    IF(tracers(iq)%parent /= 'air') CYCLE
299
300    SELECT CASE(tracers(iq)%iadv)
301    CASE(0); CYCLE
302    CASE(10); CALL   vly_loc(zq, pente_max, zm, mv, iq)
303    CASE(14); CALL vlyqs_loc(zq, pente_max, zm, mv, qsat, iq)
304    CASE DEFAULT
305      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
306    END SELECT
307
308  enddo
309
310  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 357')
311
312  do iq = 1, nqtot
313    IF(tracers(iq)%parent /= 'air') CYCLE
314    SELECT CASE(tracers(iq)%iadv)
315    CASE(0); CYCLE
316    CASE(10, 14)
317      !$OMP BARRIER
318      CALL vlz_loc(zq, pente_max, zm, mw, &
319              ij_begin, ij_end, iq)
320      !$OMP BARRIER
321
322      !$OMP MASTER
323      CALL VTb(VTHallo)
324      !$OMP END MASTER
325
326      CALL Register_Hallo_u(zq(:, :, iq), llm, 2, 2, 2, 2, MyRequest2)
327      CALL Register_Hallo_u(zm(:, :, iq), llm, 1, 1, 1, 1, MyRequest2)
328      ! CRisi
329      do ifils = 1, tracers(iq)%nqDescen
330        iq2 = tracers(iq)%iqDescen(ifils)
331        CALL Register_Hallo_u(zq(:, :, iq2), llm, 2, 2, 2, 2, MyRequest2)
332        CALL Register_Hallo_u(zm(:, :, iq2), llm, 1, 1, 1, 1, MyRequest2)
333      enddo
334      !$OMP MASTER
335      CALL VTe(VTHallo)
336      !$OMP END MASTER
337      !$OMP BARRIER
338    CASE DEFAULT
339
340      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
341    END SELECT
342
343  enddo
344  !$OMP BARRIER
345
346  !$OMP MASTER
347  CALL VTb(VTHallo)
348  !$OMP END MASTER
349
350  CALL SendRequest(MyRequest2)
351
352  !$OMP MASTER
353  CALL VTe(VTHallo)
354  !$OMP END MASTER
355
356  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 429')
357
358  !$OMP BARRIER
359  do iq = 1, nqtot
360    IF(tracers(iq)%parent /= 'air') CYCLE
361    !WRITE(*,*) 'vlspltgen 409: iq=',iq
362
363    SELECT CASE(tracers(iq)%iadv)
364    CASE(0); CYCLE
365    CASE(10, 14)
366      !$OMP BARRIER
367    CASE DEFAULT
368      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
369    END SELECT
370
371  enddo
372  !WRITE(*,*) 'vlspltgen_loc 476'
373
374  !$OMP BARRIER
375  !WRITE(*,*) 'vlspltgen_loc 477'
376  !$OMP MASTER
377  CALL VTb(VTHallo)
378  !$OMP END MASTER
379
380  ! CALL WaitRecvRequest(MyRequest2)
381  ! CALL WaitSendRequest(MyRequest2)
382  !$OMP BARRIER
383  CALL WaitRequest(MyRequest2)
384
385  !$OMP MASTER
386  CALL VTe(VTHallo)
387  !$OMP END MASTER
388  !$OMP BARRIER
389
390
391  !WRITE(*,*) 'vlspltgen_loc 494'
392  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 461')
393
394  do iq = 1, nqtot
395    IF(tracers(iq)%parent /= 'air') CYCLE
396    SELECT CASE(tracers(iq)%iadv)
397    CASE(0); CYCLE
398    CASE(10); CALL   vly_loc(zq, pente_max, zm, mv, iq)
399    CASE(14); CALL vlyqs_loc(zq, pente_max, zm, mv, qsat, iq)
400    CASE DEFAULT
401      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
402    END SELECT
403
404  enddo !do iq=1,nqtot
405
406  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 493')
407
408  do iq = 1, nqtot
409    IF(tracers(iq)%parent /= 'air') CYCLE
410    SELECT CASE(tracers(iq)%iadv)
411    CASE(0); CYCLE
412    CASE(10); CALL   vlx_loc(zq, pente_max, zm, mu, &
413            ij_begin, ij_end, iq)
414    CASE(14); CALL vlxqs_loc(zq, pente_max, zm, mu, &
415            qsat, ij_begin, ij_end, iq)
416    CASE DEFAULT
417      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
418    END SELECT
419
420  enddo !do iq=1,nqtot
421
422  !WRITE(*,*) 'vlspltgen 550: apres derniere serie de CALL vlx'
423  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 521')
424
425  ijb = ij_begin
426  ije = ij_end
427  !$OMP BARRIER
428
429  DO iq = 1, nqtot
430    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
431    DO l = 1, llm
432      DO ij = ijb, ije
433        ! print *,'zq-->',ij,l,iq,zq(ij,l,iq)
434        ! print *,'q-->',ij,l,iq,q(ij,l,iq)
435        q(ij, l, iq) = zq(ij, l, iq)
436      ENDDO
437    ENDDO
438    !$OMP END DO NOWAIT
439
440    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
441    DO l = 1, llm
442      DO ij = ijb, ije - iip1 + 1, iip1
443        q(ij + iim, l, iq) = q(ij, l, iq)
444      ENDDO
445    ENDDO
446    !$OMP END DO NOWAIT
447  ENDDO !DO iq=1,nqtot
448
449  CALL check_isotopes(q, ij_begin, ij_end, 'vlspltgen_loc 557')
450
451  !$OMP BARRIER
452
453  !c$OMP MASTER
454  ! CALL WaitSendRequest(MyRequest1)
455  ! CALL WaitSendRequest(MyRequest2)
456  !c$OMP END MASTER
457  !c$OMP BARRIER
458
459  !WRITE(*,*) 'vlspltgen 597: sortie'
460
461END SUBROUTINE vlspltgen_loc
Note: See TracBrowser for help on using the repository browser.