source: LMDZ6/branches/Amaury_dev/libf/phylmdiso/calwake.F90 @ 5105

Last change on this file since 5105 was 5105, checked in by abarral, 8 weeks ago

Replace 1DUTILS.h by module lmdz_1dutils.f90
Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files)
Convert *.F to *.f90
Fix gradsdef.h formatting
Remove unnecessary "RETURN" at the end of functions/subroutines

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