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

Last change on this file since 5276 was 5274, checked in by abarral, 23 hours ago

Replace yomcst.h by existing module

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