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

Last change on this file since 4026 was 4004, checked in by Laurent Fairhead, 3 years ago

Added properties on phylmdiso directory
LF

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