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

Last change on this file since 5159 was 5159, checked in by abarral, 3 months ago

Put dimensions.h and paramet.h into modules
  • 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 lmdz_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
34
35USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
36  USE lmdz_paramet
37  IMPLICIT NONE
38
39  !
40
41
42
43
44  !   Arguments:
45  !   ----------
46  REAL :: masse(ijb_u:ije_u, llm), pente_max
47  REAL :: pbaru(ijb_u:ije_u, llm), pbarv(ijb_v:ije_v, llm)
48  REAL :: q(ijb_u:ije_u, llm, nqtot)
49  REAL :: w(ijb_u:ije_u, llm), pdt
50  REAL :: p(ijb_u:ije_u, llmp1), teta(ijb_u:ije_u, llm)
51  REAL :: pk(ijb_u:ije_u, llm)
52
53  !  Local
54  !   ---------
55
56  INTEGER :: ij, l
57
58  REAL :: zzpbar, zzw
59
60  REAL :: qmin, qmax
61  DATA qmin, qmax/0., 1.e33/
62
63  !--pour rapport de melange saturant--
64
65  REAL :: rtt, retv, r2es, r3les, r3ies, r4les, r4ies, play
66  REAL :: ptarg, pdelarg, foeew, zdelta
67  REAL :: tempe(ijb_u:ije_u)
68  INTEGER :: ijb, ije, iq, iq2, ifils
69  LOGICAL, SAVE :: firstcall = .TRUE.
70  !$OMP THREADPRIVATE(firstcall)
71  type(request), SAVE :: MyRequest1
72  !$OMP THREADPRIVATE(MyRequest1)
73  type(request), SAVE :: MyRequest2
74  !$OMP THREADPRIVATE(MyRequest2)
75  !    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  ! Allocate variables depending on dynamic variable nqtot
90
91  IF (firstcall) THEN
92    firstcall = .FALSE.
93  END IF
94  !-- Calcul de Qsat en chaque point
95  !-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
96  !   pour eviter une exponentielle.
97
98  CALL SetTag(MyRequest1, 100)
99  CALL SetTag(MyRequest2, 101)
100
101  ijb = ij_begin - iip1
102  ije = ij_end + iip1
103  IF (pole_nord) ijb = ij_begin
104  IF (pole_sud) ije = ij_end
105
106  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
107  DO l = 1, llm
108    DO ij = ijb, ije
109      tempe(ij) = teta(ij, l) * pk(ij, l) / cpp
110    ENDDO
111    DO ij = ijb, ije
112      IF (adv_qsat_liq) THEN
113        zdelta = 0.
114      ELSE
115        zdelta = MAX(0., SIGN(1., rtt - tempe(ij)))
116      ENDIF
117      play = 0.5 * (p(ij, l) + p(ij, l + 1))
118      qsat(ij, l) = MIN(0.5, r2es * FOEEW(tempe(ij), zdelta) / play)
119      qsat(ij, l) = qsat(ij, l) / (1. - retv * qsat(ij, l))
120    ENDDO
121  ENDDO
122  !$OMP END DO NOWAIT
123  ! PRINT*,'Debut vlsplt version debug sans vlyqs'
124
125  zzpbar = 0.5 * pdt
126  zzw = pdt
127
128  ijb = ij_begin
129  ije = ij_end
130  IF (pole_nord) ijb = ijb + iip1
131  IF (pole_sud)  ije = ije - iip1
132
133  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
134  DO l = 1, llm
135    DO ij = ijb, ije
136      mu(ij, l) = pbaru(ij, l) * zzpbar
137    ENDDO
138  ENDDO
139  !$OMP END DO NOWAIT
140
141  ijb = ij_begin - iip1
142  ije = ij_end
143  IF (pole_nord) ijb = ij_begin
144  IF (pole_sud)  ije = ij_end - iip1
145
146  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
147  DO l = 1, llm
148    DO ij = ijb, ije
149      mv(ij, l) = pbarv(ij, l) * zzpbar
150    ENDDO
151  ENDDO
152  !$OMP END DO NOWAIT
153
154  ijb = ij_begin
155  ije = ij_end
156
157  DO iq = 1, nqtot
158    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
159    DO l = 1, llm
160      DO ij = ijb, ije
161        mw(ij, l, iq) = w(ij, l) * zzw
162      ENDDO
163    ENDDO
164    !$OMP END DO NOWAIT
165  ENDDO
166
167  DO iq = 1, nqtot
168    !$OMP MASTER
169    DO ij = ijb, ije
170      mw(ij, llm + 1, iq) = 0.
171    ENDDO
172    !$OMP END MASTER
173  ENDDO
174
175  ijb = ij_begin
176  ije = ij_end
177
178  DO iq = 1, nqtot
179    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
180    DO l = 1, llm
181      zq(ijb:ije, l, iq) = q(ijb:ije, l, iq)
182      zm(ijb:ije, l, iq) = masse(ijb:ije, l)
183    ENDDO
184    !$OMP END DO NOWAIT
185  ENDDO
186
187  ! verif temporaire
188  ijb = ij_begin
189  ije = ij_end
190  CALL check_isotopes(zq, ijb, ije, 'vlspltgen_loc 191')
191
192  !$OMP BARRIER
193  DO iq = 1, nqtot
194    ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
195    IF(tracers(iq)%parent /= 'air') CYCLE
196    !WRITE(*,*) 'vlspltgen 192: iq,iadv=',iq,tracers(iq)%iadv
197    SELECT CASE(tracers(iq)%iadv)
198    CASE(0); CYCLE
199    CASE(10)
200      CALL vlx_loc(zq, pente_max, zm, mu, &
201              ij_begin, ij_end, iq)
202
203      !$OMP MASTER
204      CALL VTb(VTHallo)
205      !$OMP END MASTER
206      CALL Register_Hallo_u(zq(:, :, iq), llm, 2, 2, 2, 2, MyRequest1)
207      CALL Register_Hallo_u(zm(:, :, iq), llm, 1, 1, 1, 1, MyRequest1)
208      ! CRisi
209      DO ifils = 1, tracers(iq)%nqDescen
210        iq2 = tracers(iq)%iqDescen(ifils)
211        CALL Register_Hallo_u(zq(:, :, iq2), llm, 2, 2, 2, 2, MyRequest1)
212        CALL Register_Hallo_u(zm(:, :, iq2), llm, 1, 1, 1, 1, MyRequest1)
213      enddo
214
215      !$OMP MASTER
216      CALL VTe(VTHallo)
217      !$OMP END MASTER
218    CASE(14)
219      CALL vlxqs_loc(zq, pente_max, zm, mu, &
220              qsat, ij_begin, ij_end, iq)
221
222      !$OMP MASTER
223      CALL VTb(VTHallo)
224      !$OMP END MASTER
225
226      CALL Register_Hallo_u(zq(:, :, iq), llm, 2, 2, 2, 2, MyRequest1)
227      CALL Register_Hallo_u(zm(:, :, iq), llm, 1, 1, 1, 1, MyRequest1)
228      DO ifils = 1, tracers(iq)%nqDescen
229        iq2 = tracers(iq)%iqDescen(ifils)
230        CALL Register_Hallo_u(zq(:, :, iq2), llm, 2, 2, 2, 2, MyRequest1)
231        CALL Register_Hallo_u(zm(:, :, iq2), llm, 1, 1, 1, 1, MyRequest1)
232      enddo
233
234      !$OMP MASTER
235      CALL VTe(VTHallo)
236      !$OMP END MASTER
237    CASE DEFAULT
238      CALL abort_gcm("vlspltgen_loc", "schema non parallelise", 1)
239    END SELECT
240
241  enddo !DO iq=1,nqtot
242
243
244  !$OMP BARRIER
245  !$OMP MASTER
246  CALL VTb(VTHallo)
247  !$OMP END MASTER
248
249  CALL SendRequest(MyRequest1)
250
251  !$OMP MASTER
252  CALL VTe(VTHallo)
253  !$OMP END MASTER
254  !$OMP BARRIER
255
256  ! verif temporaire
257  ijb = ij_begin - 2 * iip1
258  ije = ij_end + 2 * iip1
259  IF (pole_nord) ijb = ij_begin
260  IF (pole_sud)  ije = ij_end
261  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 280')
262
263  DO iq = 1, nqtot
264    IF(tracers(iq)%parent /= 'air') CYCLE
265    !WRITE(*,*) 'vlspltgen 279: iq=',iq
266
267    SELECT CASE(tracers(iq)%iadv)
268    CASE(0); CYCLE
269    CASE(10)
270    CASE(14)
271    CASE DEFAULT
272      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
273    END SELECT
274
275  enddo
276  !$OMP BARRIER
277  !$OMP MASTER
278  CALL VTb(VTHallo)
279  !$OMP END MASTER
280
281  ! CALL WaitRecvRequest(MyRequest1)
282  ! CALL WaitSendRequest(MyRequest1)
283  !$OMP BARRIER
284  CALL WaitRequest(MyRequest1)
285
286
287  !$OMP MASTER
288  CALL VTe(VTHallo)
289  !$OMP END MASTER
290  !$OMP BARRIER
291
292  IF(isoCheck) THEN
293    CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 326')
294    ijb = ij_begin - 2 * iip1
295    ije = ij_end + 2 * iip1
296    IF (pole_nord) ijb = ij_begin
297    IF (pole_sud)  ije = ij_end
298    CALL check_isotopes(zq, ijb, ije, 'vlspltgen_loc 336')
299  END IF
300
301  DO iq = 1, nqtot
302    IF(tracers(iq)%parent /= 'air') CYCLE
303
304    SELECT CASE(tracers(iq)%iadv)
305    CASE(0); CYCLE
306    CASE(10); CALL   vly_loc(zq, pente_max, zm, mv, iq)
307    CASE(14); CALL vlyqs_loc(zq, pente_max, zm, mv, qsat, iq)
308    CASE DEFAULT
309      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
310    END SELECT
311
312  enddo
313
314  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 357')
315
316  DO iq = 1, nqtot
317    IF(tracers(iq)%parent /= 'air') CYCLE
318    SELECT CASE(tracers(iq)%iadv)
319    CASE(0); CYCLE
320    CASE(10, 14)
321      !$OMP BARRIER
322      CALL vlz_loc(zq, pente_max, zm, mw, &
323              ij_begin, ij_end, iq)
324      !$OMP BARRIER
325
326      !$OMP MASTER
327      CALL VTb(VTHallo)
328      !$OMP END MASTER
329
330      CALL Register_Hallo_u(zq(:, :, iq), llm, 2, 2, 2, 2, MyRequest2)
331      CALL Register_Hallo_u(zm(:, :, iq), llm, 1, 1, 1, 1, MyRequest2)
332      ! CRisi
333      DO ifils = 1, tracers(iq)%nqDescen
334        iq2 = tracers(iq)%iqDescen(ifils)
335        CALL Register_Hallo_u(zq(:, :, iq2), llm, 2, 2, 2, 2, MyRequest2)
336        CALL Register_Hallo_u(zm(:, :, iq2), llm, 1, 1, 1, 1, MyRequest2)
337      enddo
338      !$OMP MASTER
339      CALL VTe(VTHallo)
340      !$OMP END MASTER
341      !$OMP BARRIER
342    CASE DEFAULT
343
344      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
345    END SELECT
346
347  enddo
348  !$OMP BARRIER
349
350  !$OMP MASTER
351  CALL VTb(VTHallo)
352  !$OMP END MASTER
353
354  CALL SendRequest(MyRequest2)
355
356  !$OMP MASTER
357  CALL VTe(VTHallo)
358  !$OMP END MASTER
359
360  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 429')
361
362  !$OMP BARRIER
363  DO iq = 1, nqtot
364    IF(tracers(iq)%parent /= 'air') CYCLE
365    !WRITE(*,*) 'vlspltgen 409: iq=',iq
366
367    SELECT CASE(tracers(iq)%iadv)
368    CASE(0); CYCLE
369    CASE(10, 14)
370      !$OMP BARRIER
371    CASE DEFAULT
372      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
373    END SELECT
374
375  enddo
376  !WRITE(*,*) 'vlspltgen_loc 476'
377
378  !$OMP BARRIER
379  !WRITE(*,*) 'vlspltgen_loc 477'
380  !$OMP MASTER
381  CALL VTb(VTHallo)
382  !$OMP END MASTER
383
384  ! CALL WaitRecvRequest(MyRequest2)
385  ! CALL WaitSendRequest(MyRequest2)
386  !$OMP BARRIER
387  CALL WaitRequest(MyRequest2)
388
389  !$OMP MASTER
390  CALL VTe(VTHallo)
391  !$OMP END MASTER
392  !$OMP BARRIER
393
394
395  !WRITE(*,*) 'vlspltgen_loc 494'
396  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 461')
397
398  DO iq = 1, nqtot
399    IF(tracers(iq)%parent /= 'air') CYCLE
400    SELECT CASE(tracers(iq)%iadv)
401    CASE(0); CYCLE
402    CASE(10); CALL   vly_loc(zq, pente_max, zm, mv, iq)
403    CASE(14); CALL vlyqs_loc(zq, pente_max, zm, mv, qsat, iq)
404    CASE DEFAULT
405      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
406    END SELECT
407
408  enddo !do iq=1,nqtot
409
410  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 493')
411
412  DO iq = 1, nqtot
413    IF(tracers(iq)%parent /= 'air') CYCLE
414    SELECT CASE(tracers(iq)%iadv)
415    CASE(0); CYCLE
416    CASE(10); CALL   vlx_loc(zq, pente_max, zm, mu, &
417            ij_begin, ij_end, iq)
418    CASE(14); CALL vlxqs_loc(zq, pente_max, zm, mu, &
419            qsat, ij_begin, ij_end, iq)
420    CASE DEFAULT
421      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
422    END SELECT
423
424  enddo !do iq=1,nqtot
425
426  !WRITE(*,*) 'vlspltgen 550: apres derniere serie de CALL vlx'
427  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 521')
428
429  ijb = ij_begin
430  ije = ij_end
431  !$OMP BARRIER
432
433  DO iq = 1, nqtot
434    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
435    DO l = 1, llm
436      DO ij = ijb, ije
437        ! print *,'zq-->',ij,l,iq,zq(ij,l,iq)
438        ! print *,'q-->',ij,l,iq,q(ij,l,iq)
439        q(ij, l, iq) = zq(ij, l, iq)
440      ENDDO
441    ENDDO
442    !$OMP END DO NOWAIT
443
444    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
445    DO l = 1, llm
446      DO ij = ijb, ije - iip1 + 1, iip1
447        q(ij + iim, l, iq) = q(ij, l, iq)
448      ENDDO
449    ENDDO
450    !$OMP END DO NOWAIT
451  ENDDO !DO iq=1,nqtot
452
453  CALL check_isotopes(q, ij_begin, ij_end, 'vlspltgen_loc 557')
454
455  !$OMP BARRIER
456
457  !c$OMP MASTER
458  ! CALL WaitSendRequest(MyRequest1)
459  ! CALL WaitSendRequest(MyRequest2)
460  !c$OMP END MASTER
461  !c$OMP BARRIER
462
463  !WRITE(*,*) 'vlspltgen 597: sortie'
464
465END SUBROUTINE vlspltgen_loc
Note: See TracBrowser for help on using the repository browser.