source: LMDZ6/trunk/libf/phylmdiso/calwake.F90 @ 4657

Last change on this file since 4657 was 4594, checked in by fhourdin, 18 months ago

Correction pour phylmdiso

  • Property svn:keywords set to Id
File size: 14.0 KB
Line 
1
2! $Id: calwake.F90 4594 2023-06-29 16:09:08Z fhourdin $
3
4SUBROUTINE calwake(iflag_wake_tend, paprs, pplay, dtime, &
5    t, q, omgb, &
6    dt_dwn, dq_dwn, m_dwn, m_up, dt_a, dq_a, wgen, &
7    sigd, Cin, &
8    wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens, &
9    wake_dth, wake_h, &
10    wake_pe, wake_fip, wake_gfl, &
11    dt_wake, dq_wake, wake_k, t_x, q_x, wake_omgbdth, &
12    wake_dp_omgb, &
13    wake_dtke, wake_dqke, &
14    wake_omg, wake_dp_deltomg, &
15    wake_spread, wake_cstar, wake_d_deltat_gw, &
16    wake_ddeltat, wake_ddeltaq, wake_ds, awake_ddens, wake_ddens &
17#ifdef ISO
18     &               ,xt,dxt_dwn,dxt_a &
19     &               ,wake_deltaxt,dxt_wake,xt_x,wake_ddeltaxt &
20#endif     
21     &   )
22  ! **************************************************************
23  ! *
24  ! CALWAKE                                                     *
25  ! interface avec le schema de calcul de la poche    *
26  ! froide                                            *
27  ! *
28  ! written by   : CHERUY Frederique, 13/03/2000, 10.31.05      *
29  ! modified by :  ROEHRIG Romain,    01/30/2007                *
30  ! **************************************************************
31
32  USE dimphy
33  USE phys_state_var_mod, ONLY: pctsrf
34  USE indice_sol_mod, ONLY: is_oce
35  USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
36#ifdef ISO
37  USE infotrac_phy, ONLY : ntraciso=>ntiso
38#ifdef ISOVERIF
39  USE isotopes_mod, ONLY: iso_eau
40  USE isotopes_verif_mod
41  USE lmdz_wake, ONLY : wake
42
43#endif
44#endif
45  IMPLICIT NONE
46  ! ======================================================================
47  include "YOMCST.h"
48
49  ! Arguments
50  ! ----------
51  ! Input
52  ! ----
53  INTEGER,                       INTENT (IN)         :: iflag_wake_tend
54  REAL,                          INTENT (IN)         :: dtime
55  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: pplay
56  REAL, DIMENSION(klon, klev+1), INTENT (IN)         :: paprs
57  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: t, q, omgb
58  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: dt_dwn, dq_dwn
59  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: m_up, m_dwn
60  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: dt_a, dq_a
61  REAL, DIMENSION(klon),         INTENT (IN)         :: wgen
62  REAL, DIMENSION(klon),         INTENT (IN)         :: sigd
63  REAL, DIMENSION(klon),         INTENT (IN)         :: Cin
64#ifdef ISO
65  REAL, DIMENSION(ntraciso,klon, klev),   INTENT (IN) :: xt,dxt_dwn,dxt_a
66#endif
67  ! Input/Output
68  ! ------------
69  REAL, DIMENSION(klon, klev),   INTENT (INOUT)      :: wake_deltat, wake_deltaq
70  REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_s
71  REAL, DIMENSION(klon),         INTENT (INOUT)      :: awake_dens, wake_dens
72#ifdef ISO
73  REAL, DIMENSION(ntraciso,klon, klev),   INTENT (INOUT)      :: wake_deltaxt
74#endif
75  ! Output
76  ! ------
77  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: dt_wake, dq_wake
78!!jyg  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_k
79  INTEGER, DIMENSION(klon),      INTENT (OUT)        :: wake_k
80  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_d_deltat_gw
81  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_h
82  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_dth
83  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_pe, wake_fip, wake_gfl
84  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: t_x, q_x
85  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_omgbdth, wake_dp_omgb
86  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_dtke, wake_dqke
87  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_omg, wake_dp_deltomg
88  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_spread
89  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_cstar
90  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_ddeltat, wake_ddeltaq
91  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_ds, awake_ddens, wake_ddens
92#ifdef ISO
93  REAL, DIMENSION(ntraciso,klon, klev),   INTENT (OUT)        :: dxt_wake
94  REAL, DIMENSION(ntraciso,klon, klev),   INTENT (OUT)        :: xt_x
95  REAL, DIMENSION(ntraciso,klon, klev),   INTENT (OUT)        :: wake_ddeltaxt
96  REAL, DIMENSION(ntraciso,klon, klev)                        :: wake_dxtke ! pas besoin de la sortir
97#endif
98
99
100  ! Variable internes
101  ! -----------------
102  LOGICAL, SAVE                                      :: first = .TRUE.
103  !$OMP THREADPRIVATE(first)
104  INTEGER                                            :: i, l
105  INTEGER, DIMENSION(klon)                           :: znatsurf    ! 0 if pctsrf(is_oce)>0.1; 1 else.
106  REAL                                               :: aire
107  REAL, DIMENSION(klon, klev)                        :: p,  pi
108  REAL, DIMENSION(klon, klev+1)                      ::  ph
109  REAL, DIMENSION(klon, klev)                        ::  omgbe
110  REAL, DIMENSION(klon, klev)                        :: te, qe
111  REAL, DIMENSION(klon, klev)                        :: dtdwn, dqdwn
112  REAL, DIMENSION(klon, klev)                        :: dta, dqa
113  REAL, DIMENSION(klon, klev)                        :: amdwn, amup
114  REAL, DIMENSION(klon, klev)                        :: dtw, dqw, dth
115  REAL, DIMENSION(klon, klev)                        :: dtls, dqls
116  REAL, DIMENSION(klon, klev)                        :: tx, qx
117  REAL, DIMENSION(klon)                              :: hw, wape, fip, gfl
118  REAL, DIMENSION(klon)                              :: sigmaw, awdens, wdens
119  REAL, DIMENSION(klon, klev)                        :: omgbdth
120  REAL, DIMENSION(klon, klev)                        :: dp_omgb
121  REAL, DIMENSION(klon, klev)                        :: dtke, dqke
122  REAL, DIMENSION(klon, klev)                        :: omg
123  REAL, DIMENSION(klon, klev)                        :: dp_deltomg, spread
124  REAL, DIMENSION(klon)                              :: cstar
125  REAL, DIMENSION(klon)                              :: sigd0
126  INTEGER, DIMENSION(klon)                           :: ktopw
127  REAL, DIMENSION(klon, klev)                        :: d_deltat_gw
128  REAL, DIMENSION(klon, klev)                        :: d_deltatw, d_deltaqw
129  REAL, DIMENSION(klon)                              :: d_sigmaw, d_awdens, d_wdens
130#ifdef ISO
131  REAL, DIMENSION(ntraciso,klon, klev)                        :: xte
132  REAL, DIMENSION(ntraciso,klon, klev)                        :: dxtdwn
133  REAL, DIMENSION(ntraciso,klon, klev)                        :: dxta
134  REAL, DIMENSION(ntraciso,klon, klev)                        :: dxtw
135  REAL, DIMENSION(ntraciso,klon, klev)                        :: dxtls
136  REAL, DIMENSION(ntraciso,klon, klev)                        :: xtx
137  REAL, DIMENSION(ntraciso,klon, klev)                        :: dxtke
138  REAL, DIMENSION(ntraciso,klon, klev)                        :: d_deltaxtw
139  INTEGER                                            :: ixt
140#endif
141
142  REAL                                               :: rdcp
143
144  IF (prt_level >= 10) THEN
145    print *, '-> calwake, wake_s, wgen input ', wake_s(1), wgen(1)
146  ENDIF
147
148  rdcp = 1./3.5
149
150  znatsurf(:) = 0
151  DO i = 1,klon
152    IF (pctsrf(i,is_oce) < 0.1) znatsurf(i) = 1
153  ENDDO
154
155
156  ! -----------------------------------------------------------
157  ! IM 290108     DO 999 i=1,klon   ! a vectoriser
158  ! ----------------------------------------------------------
159
160
161  DO l = 1, klev
162    DO i = 1, klon
163      p(i, l) = pplay(i, l)
164      ph(i, l) = paprs(i, l)
165      pi(i, l) = (pplay(i,l)/100000.)**rdcp
166
167      te(i, l) = t(i, l)
168      qe(i, l) = q(i, l)
169      omgbe(i, l) = omgb(i, l)
170
171      dtdwn(i, l) = dt_dwn(i, l)
172      dqdwn(i, l) = dq_dwn(i, l)
173      dta(i, l) = dt_a(i, l)
174      dqa(i, l) = dq_a(i, l)
175#ifdef ISO
176      do ixt=1,ntraciso
177        xte(ixt,i, l) = xt(ixt,i, l)
178        dxtdwn(ixt,i, l) = dxt_dwn(ixt,i, l)
179        dxta(ixt,i, l) = dxt_a(ixt,i, l)
180      enddo
181#endif
182    END DO
183  END DO
184
185!----------------------------------------------------------------
186!         Initialize tendencies to zero
187!----------------------------------------------------------------
188dtls(:,:) = 0.
189dqls(:,:) = 0.
190d_deltat_gw(:,:) = 0.
191d_deltatw(:,:) = 0.
192d_deltaqw(:,:) = 0.
193d_sigmaw(:) = 0.
194d_awdens(:) = 0.
195d_wdens(:) = 0.
196#ifdef ISO
197dxtls(:,:,:) = 0.
198d_deltaxtw(:,:,:) = 0.
199#endif
200!
201
202  DO i = 1, klon
203    sigd0(i) = sigd(i)
204  END DO
205  ! print*, 'sigd0,sigd', sigd0, sigd(i)
206  DO i = 1, klon
207    ph(i, klev+1) = 0.
208  END DO
209
210!!jyg!  DO i = 1, klon                 
211!!jyg!    ktopw(i) = NINT(wake_k(i))   
212!!jyg!  END DO                         
213
214  DO i = 1, klon
215    hw(i) = wake_h(i)
216  END DO
217!
218!    Make a copy of state variables
219  DO l = 1, klev
220    DO i = 1, klon
221      dtw(i, l) = wake_deltat(i, l)
222      dqw(i, l) = wake_deltaq(i, l)
223#ifdef ISO
224      dxtw(:,i, l) = wake_deltaxt(:,i, l)
225#endif
226    END DO
227  END DO
228
229  DO i = 1, klon
230    sigmaw(i) = wake_s(i)
231  END DO
232
233  DO i = 1, klon
234    awdens(i) = max(0., awake_dens(i))
235    wdens(i) = max(0., wake_dens(i))
236  END DO
237
238  ! fkc les flux de masses sont evalues aux niveaux et valent 0 a la surface
239  ! fkc  on veut le flux de masse au milieu des couches
240
241  DO l = 1, klev - 1
242    DO i = 1, klon
243      amdwn(i, l) = 0.5*(m_dwn(i,l)+m_dwn(i,l+1))
244      amdwn(i, l) = (m_dwn(i,l+1))
245    END DO
246  END DO
247
248  ! au sommet le flux de masse est nul
249
250  DO i = 1, klon
251    amdwn(i, klev) = 0.5*m_dwn(i, klev)
252  END DO
253
254  DO l = 1, klev
255    DO i = 1, klon
256      amup(i, l) = m_up(i, l)
257    END DO
258  END DO
259
260#ifdef ISOVERIF
261! verif des inputs des wakes
262        write(*,*) 'calwake 257: verif des inputs des wakes'
263  DO l = 1, klev
264    DO i = 1, klon
265      if (iso_eau.gt.0) then
266        call iso_verif_egalite(qe(i,l),xte(iso_eau,i,l),'calwake 261a')
267        call iso_verif_egalite(dqa(i,l),dxta(iso_eau,i,l),'calwake 261b')
268        call iso_verif_egalite(dqdwn(i,l),dxtdwn(iso_eau,i,l),'calwake 261c')
269        call iso_verif_egalite(dqw(i,l),dxtw(iso_eau,i,l),'calwake 261d')
270      endif
271    END DO
272  END DO
273#endif
274
275  CALL wake(znatsurf, p, ph, pi, dtime, &
276    te, qe, omgbe, &
277    dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
278    sigd0, Cin, &
279    dtw, dqw, sigmaw, awdens, wdens, &                                   ! state variables
280    dth, hw, wape, fip, gfl, &
281    dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
282    dtke, dqke, omg, dp_deltomg, spread, cstar, &
283    d_deltat_gw, &
284    d_deltatw, d_deltaqw, d_sigmaw, d_awdens, d_wdens &                     ! tendencies
285#ifdef ISO
286     , xte,dxtdwn,dxta,dxtw &
287     , dxtls,xtx,dxtke,d_deltaxtw &
288#endif
289    )
290#ifdef ISOVERIF
291! verif des ouputs des wakes
292        write(*,*) 'calwake 286: verif des outputs des wakes'
293  DO l = 1, klev
294    DO i = 1, klon
295      if (iso_eau.gt.0) then
296        call iso_verif_egalite(dqls(i,l),dxtls(iso_eau,i,l),'calwake 290a')
297        call iso_verif_egalite(qx(i,l),xtx(iso_eau,i,l),'calwake 290b')
298        call iso_verif_egalite(dqke(i,l),dxtke(iso_eau,i,l),'calwake 290c')
299        call iso_verif_egalite(d_deltaqw(i,l),d_deltaxtw(iso_eau,i,l),'calwake 290d')
300        call iso_verif_egalite(dqw(i,l),dxtw(iso_eau,i,l),'calwake 290e')
301      endif
302    END DO
303  END DO
304#endif
305
306!
307  DO l = 1, klev
308    DO i = 1, klon
309      IF (ktopw(i)>0) THEN
310        wake_d_deltat_gw(i, l) = d_deltat_gw(i, l)
311        wake_omgbdth(i, l) = omgbdth(i, l)
312        wake_dp_omgb(i, l) = dp_omgb(i, l)
313        wake_dtke(i, l) = dtke(i, l)
314        wake_dqke(i, l) = dqke(i, l)
315        wake_omg(i, l) = omg(i, l)
316        wake_dp_deltomg(i, l) = dp_deltomg(i, l)
317        wake_spread(i, l) = spread(i, l)
318        wake_dth(i, l) = dth(i, l)
319        dt_wake(i, l) = dtls(i, l)*dtime         ! derivative -> tendency
320        dq_wake(i, l) = dqls(i, l)*dtime         ! derivative -> tendency
321        t_x(i, l) = tx(i, l)
322        q_x(i, l) = qx(i, l)
323#ifdef ISO
324        do ixt=1,ntraciso
325          wake_dxtke(ixt,i, l) = dxtke(ixt,i, l)
326          dxt_wake(ixt,i, l) = dxtls(ixt,i, l)*dtime
327          xt_x(ixt,i, l) = xtx(ixt,i, l)
328        enddo
329#endif
330      ELSE
331        wake_d_deltat_gw(i, l) = 0.
332        wake_omgbdth(i, l) = 0.
333        wake_dp_omgb(i, l) = 0.
334        wake_dtke(i, l) = 0.
335        wake_dqke(i, l) = 0.
336        wake_omg(i, l) = 0.
337        wake_dp_deltomg(i, l) = 0.
338        wake_spread(i, l) = 0.
339        wake_dth(i, l) = 0.
340        dt_wake(i, l) = 0.
341        dq_wake(i, l) = 0.
342        t_x(i, l) = te(i, l)
343        q_x(i, l) = qe(i, l)
344#ifdef ISO
345        do ixt=1,ntraciso
346          wake_dxtke(ixt,i, l) = 0.
347          dxt_wake(ixt,i, l) = 0.
348          xt_x(ixt,i, l) = xte(ixt,i, l)
349        enddo
350#endif
351      END IF
352    END DO
353  END DO
354
355  DO i = 1, klon
356    wake_h(i) = hw(i)
357    wake_pe(i) = wape(i)
358    wake_fip(i) = fip(i)
359    wake_gfl(i) = gfl(i)
360    wake_k(i) = ktopw(i)
361    wake_cstar(i) = cstar(i)
362  END DO
363
364!  Tendencies of state variables
365  DO l = 1, klev
366    DO i = 1, klon
367      IF (ktopw(i)>0) THEN
368        wake_ddeltat(i, l) = d_deltatw(i, l)*dtime
369        wake_ddeltaq(i, l) = d_deltaqw(i, l)*dtime
370#ifdef ISO
371        do ixt=1,ntraciso
372          wake_ddeltaxt(ixt,i, l) = d_deltaxtw(ixt,i, l)*dtime
373        enddo
374#endif
375      ELSE
376        wake_ddeltat(i, l) = -wake_deltat(i, l)
377        wake_ddeltaq(i, l) = -wake_deltaq(i, l)
378#ifdef ISO
379        do ixt=1,ntraciso
380          wake_ddeltaxt(ixt,i, l) = -wake_deltaxt(ixt,i, l)
381        enddo
382#endif
383      END IF
384    END DO
385  END DO
386  DO i = 1, klon
387    IF (ktopw(i)>0) THEN
388      wake_ds(i) = d_sigmaw(i)*dtime
389      awake_ddens(i) = d_awdens(i)*dtime
390      wake_ddens(i) = d_wdens(i)*dtime
391    ELSE
392      wake_ds(i)   = -wake_s(i)
393      wake_ddens(i)= -wake_dens(i)
394    END IF
395  END DO
396!
397
398!jyg< 
399  IF (iflag_wake_tend .EQ. 0) THEN
400!  Update State variables
401    DO l = 1, klev
402      DO i = 1, klon
403        IF (ktopw(i)>0) THEN
404          wake_deltat(i, l) = dtw(i, l)
405          wake_deltaq(i, l) = dqw(i, l)
406#ifdef ISO
407        do ixt=1,ntraciso
408          wake_deltaxt(ixt,i, l) = dxtw(ixt,i, l)
409        enddo
410#endif
411        ELSE
412          wake_deltat(i, l) = 0.
413          wake_deltaq(i, l) = 0.
414#ifdef ISO
415        do ixt=1,ntraciso
416          wake_deltaxt(ixt,i, l) = 0.
417        enddo
418#endif
419        END IF
420      END DO
421    END DO
422    DO i = 1, klon
423      wake_s(i) = sigmaw(i)
424      awake_dens(i) = awdens(i)
425      wake_dens(i) = wdens(i)
426    END DO
427  ENDIF  ! (iflag_wake_tend .EQ. 0)
428!
429  IF (first) THEN
430    DO i = 1,klon
431      IF (wake_dens(i) < -1.) THEN
432        wake_dens(i) = wdens(i)
433      ENDIF
434    ENDDO
435    first=.false.
436  ENDIF  ! (first)
437!>jyg
438
439  RETURN
440END SUBROUTINE calwake
441
Note: See TracBrowser for help on using the repository browser.