GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/lmdz_thermcell_alp.F90 Lines: 144 152 94.7 %
Date: 2023-06-30 12:56:34 Branches: 123 142 86.6 %

Line Branch Exec Source
1
MODULE lmdz_thermcell_alp
2
! $Id: thermcell_main.F90 2351 2015-08-25 15:14:59Z emillour $
3
!
4
CONTAINS
5
6
11377955
      SUBROUTINE thermcell_alp(ngrid,nlay,ptimestep  &                         ! in
7
288
     &                  ,pplay,pplev  &                                        ! in
8
     &                  ,fm0,entr0,lmax  &                                     ! in
9
288
     &                  ,pbl_tke,pctsrf,omega,airephy &                        ! in
10
     &                  ,zw2,fraca &                                           ! in
11
     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &  ! in
12
!
13
288
     &                  ,ale_bl,alp_bl,lalim_conv,wght_th &                    ! out
14
288
     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &   ! out
15
     &                  ,n2,s2,ale_bl_stat &                                   ! out
16
288
     &                  ,therm_tke_max,env_tke_max &                           ! out
17
     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &          ! out
18
     &                  ,alp_bl_conv,alp_bl_stat &                             ! out
19
     &)
20
21
      USE indice_sol_mod
22
      USE lmdz_thermcell_main, ONLY : thermcell_tke_transport
23
      IMPLICIT NONE
24
25
!=======================================================================
26
!
27
!   Auteurs: Catherine Rio
28
!   Modifications :
29
!   Nicolas Rochetin et Jean-Yves Grandpeix
30
!         pour la fermeture stochastique. 2012
31
!   Fr�d�ric Hourdin :
32
!         netoyage informatique. 2022
33
!
34
!=======================================================================
35
!-----------------------------------------------------------------------
36
!   declarations:
37
!   -------------
38
39
      INCLUDE "YOMCST.h"
40
      INCLUDE "YOETHF.h"
41
      INCLUDE "FCTTRE.h"
42
      INCLUDE "alpale.h"
43
44
!   arguments:
45
!   ----------
46
47
!------Entrees
48
      integer, intent(in) :: ngrid,nlay
49
      real, intent(in) :: ptimestep
50
      real, intent(in) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
51
      integer, intent(in), dimension(ngrid) ::lmax,lalim
52
      real, intent(in), dimension(ngrid) :: zmax
53
      real, intent(in), dimension(ngrid,nlay+1) :: zw2
54
      real, intent(in), dimension(ngrid,nlay+1) :: fraca
55
      real, intent(in), dimension(ngrid,nlay) :: wth3
56
      real, intent(in), dimension(ngrid,nlay) :: rhobarz
57
      real, intent(in), dimension(ngrid) :: wmax_sec
58
      real, intent(in), dimension(ngrid,nlay) :: entr0
59
      real, intent(in), dimension(ngrid,nlay+1) :: fm0,fm
60
      real, intent(in), dimension(ngrid) :: pcon
61
      real, intent(in), dimension(ngrid,nlay) :: alim_star
62
      real, intent(in), dimension(ngrid,nlay+1,nbsrf) :: pbl_tke
63
      real, intent(in), dimension(ngrid,nbsrf) :: pctsrf
64
      real, intent(in), dimension(ngrid,nlay) :: omega
65
      real, intent(in), dimension(ngrid) :: airephy
66
!------Sorties
67
      real, intent(out), dimension(ngrid) :: ale_bl,alp_bl
68
      real, intent(out), dimension(ngrid,nlay) :: wght_th
69
      integer, intent(out), dimension(ngrid) :: lalim_conv
70
      real, intent(out), dimension(ngrid) :: zlcl,fraca0,w0,w_conv
71
      real, intent(out), dimension(ngrid) :: therm_tke_max0,env_tke_max0,n2,s2,ale_bl_stat
72
      real, intent(out), dimension(ngrid,nlay) :: therm_tke_max,env_tke_max
73
      real, intent(out), dimension(ngrid) :: alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke
74
      real, intent(out), dimension(ngrid) :: alp_bl_conv,alp_bl_stat
75
76
!=============================================================================================
77
!------Local
78
!=============================================================================================
79
80
      REAL susqr2pi, reuler
81
      INTEGER ig,k,l
82
      integer nsrf
83
576
      real rhobarz0(ngrid)                    ! Densit� au LCL
84
576
      logical ok_lcl(ngrid)                   ! Existence du LCL des thermiques
85
576
      integer klcl(ngrid)                     ! Niveau du LCL
86
576
      real interp(ngrid)                      ! Coef d'interpolation pour le LCL
87
!--Triggering
88
      real, parameter :: su_cst=4e4              ! Surface unite: celle d'un updraft �l�mentaire
89
      real, parameter :: hcoef=1             ! Coefficient directeur pour le calcul de s2
90
      real, parameter :: hmincoef=0.3        ! Coefficient directeur pour l'ordonn�e � l'origine pour le calcul de s2
91
      real, parameter :: eps1=0.3            ! Fraction de surface occup�e par la population 1 : eps1=n1*s1/(fraca0*Sd)
92
576
      real, dimension(ngrid) :: hmin         ! Ordonn�e � l'origine pour le calcul de s2
93
576
      real, dimension(ngrid) :: zmax_moy     ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
94
      real, parameter :: zmax_moy_coef=0.33
95
576
      real, dimension(ngrid) :: depth        ! Epaisseur moyenne du cumulus
96
576
      real, dimension(ngrid) ::  w_max                 ! Vitesse max statistique
97
576
      real, dimension(ngrid) ::  s_max(ngrid)
98
!--Closure
99
576
      real, dimension(ngrid,nlay) :: pbl_tke_max       ! Profil de TKE moyenne
100
576
      real, dimension(ngrid) :: pbl_tke_max0           ! TKE moyenne au LCL
101
576
      real, dimension(ngrid,nlay) :: w_ls              ! Vitesse verticale grande �chelle (m/s)
102
      real, parameter :: coef_m=1.            ! On consid�re un rendement pour alp_bl_fluct_m
103
      real, parameter :: coef_tke=1.          ! On consid�re un rendement pour alp_bl_fluct_tke
104
      real :: zdp
105
576
      real, dimension(ngrid) :: alp_int,dp_int
106
288
      real, dimension(ngrid) :: fm_tot
107
108
!------------------------------------------------------------
109
!  Initialize output arrays related to stochastic triggering
110
!------------------------------------------------------------
111
286560
  DO ig = 1,ngrid
112
286272
     zlcl(ig) = 0.
113
286272
     fraca0(ig) = 0.
114
286272
     w0(ig) = 0.
115
286272
     w_conv(ig) = 0.
116
286272
     therm_tke_max0(ig) = 0.
117
286272
     env_tke_max0(ig) = 0.
118
286272
     n2(ig) = 0.
119
286272
     s2(ig) = 0.
120
286272
     ale_bl_stat(ig) = 0.
121
286272
     alp_bl_det(ig) = 0.
122
286272
     alp_bl_fluct_m(ig) = 0.
123
286272
     alp_bl_fluct_tke(ig) = 0.
124
286272
     alp_bl_conv(ig) = 0.
125
286560
     alp_bl_stat(ig) = 0.
126
  ENDDO
127
11520
  DO l = 1,nlay
128
11176128
    DO ig = 1,ngrid
129
11164608
     therm_tke_max(ig,l) = 0.
130
11175840
     env_tke_max(ig,l) = 0.
131
    ENDDO
132
  ENDDO
133
134
!------------Test sur le LCL des thermiques
135
286560
    do ig=1,ngrid
136
286272
      ok_lcl(ig)=.false.
137

286560
      if ( (pcon(ig) .gt. pplay(ig,nlay-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
138
    enddo
139
140
!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
141
11232
    do l=1,nlay-1
142
10889568
      do ig=1,ngrid
143
10889280
        if (ok_lcl(ig)) then
144
!ATTENTION,zw2 calcule en pplev
145
!          if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
146
!          klcl(ig)=l
147
!          interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
148
!          endif
149

8041522
          if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then
150
211619
          klcl(ig)=l
151
211619
          interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(ig,klcl(ig)))
152
          endif
153
        endif
154
      enddo
155
    enddo
156
157
286560
    do ig =1,ngrid
158
!CR:REHABILITATION ZMAX CONTINU
159
286560
     if (ok_lcl(ig)) then
160
      rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
161
211619
 &               -rhobarz(ig,klcl(ig)))*interp(ig)
162
211619
      zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
163
211619
      zlcl(ig)=min(zlcl(ig),zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
164
     else
165
74653
      rhobarz0(ig)=0.
166
74653
      zlcl(ig)=zmax(ig)
167
     endif
168
    enddo
169
!!jyg fin
170
171
!------------Calcul des propri�t�s du thermique au LCL
172

288
  IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
173
174
  !-----Initialisation de la TKE moyenne
175
11520
   do l=1,nlay
176
11176128
    do ig=1,ngrid
177
11175840
     pbl_tke_max(ig,l)=0.
178
    enddo
179
   enddo
180
181
!-----Calcul de la TKE moyenne
182
1440
   do nsrf=1,nbsrf
183
46368
    do l=1,nlay
184
44704512
     do ig=1,ngrid
185
44703360
     pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
186
     enddo
187
    enddo
188
   enddo
189
190
!-----Initialisations des TKE dans et hors des thermiques
191
11520
   do l=1,nlay
192
11176128
    do ig=1,ngrid
193
11164608
    therm_tke_max(ig,l)=pbl_tke_max(ig,l)
194
11175840
    env_tke_max(ig,l)=pbl_tke_max(ig,l)
195
    enddo
196
   enddo
197
198
!-----Calcul de la TKE transport�e par les thermiques : therm_tke_max
199
   call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &  ! in
200
288
  &           rg,pplev,therm_tke_max)                               ! out
201
!   print *,' thermcell_tke_transport -> '   !!jyg
202
203
!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande �chelle : W_ls
204
11520
   do l=1,nlay
205
11176128
    do ig=1,ngrid
206
11164608
     pbl_tke_max(ig,l)=fraca(ig,l)*therm_tke_max(ig,l)+(1.-fraca(ig,l))*env_tke_max(ig,l)         !  Recalcul de TKE moyenne apr�s transport de TKE_TH
207
11164608
     env_tke_max(ig,l)=(pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l))       !  Recalcul de TKE dans  l'environnement apr�s transport de TKE_TH
208
11175840
     w_ls(ig,l)=-1.*omega(ig,l)/(RG*rhobarz(ig,l))                                                !  Vitesse verticale de grande �chelle
209
    enddo
210
   enddo
211
!    print *,' apres w_ls = '   !!jyg
212
213
286560
  do ig=1,ngrid
214
286560
   if (ok_lcl(ig)) then
215
     fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
216
211619
 &             -fraca(ig,klcl(ig)))*interp(ig)
217
     w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
218
211619
 &         -zw2(ig,klcl(ig)))*interp(ig)
219
     w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
220
211619
 &             -w_ls(ig,klcl(ig)))*interp(ig)
221
     therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
222
211619
 &                     +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
223
     env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
224
211619
 &                   -env_tke_max(ig,klcl(ig)))*interp(ig)
225
     pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
226
211619
 &                   -pbl_tke_max(ig,klcl(ig)))*interp(ig)
227
211619
     if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
228
211619
     if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
229
211619
     if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
230
   else
231
74653
     fraca0(ig)=0.
232
74653
     w0(ig)=0.
233
!!jyg le 27/04/2012
234
!!     zlcl(ig)=0.
235
!!
236
   endif
237
  enddo
238
239
  ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
240
!  print *,'ENDIF  ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) '    !!jyg
241
242
!------------Triggering------------------
243
288
  IF (iflag_trig_bl.ge.1) THEN
244
245
!-----Initialisations
246
286560
   depth(:)=0.
247
286560
   n2(:)=0.
248
286560
   s2(:)=100. ! some low value, arbitrary
249
286560
   s_max(:)=0.
250
251
!-----Epaisseur du nuage (depth) et d�termination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
252
286560
   do ig=1,ngrid
253
286272
     zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
254
286272
     depth(ig)=zmax_moy(ig)-zlcl(ig)
255
286272
     hmin(ig)=hmincoef*zlcl(ig)
256
286560
     if (depth(ig).ge.10.) then
257
66778
       s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
258
66778
       n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
259
!!
260
!!jyg le 27/04/2012
261
!!       s_max(ig)=s2(ig)*log(n2(ig))
262
!!       if (n2(ig) .lt. 1) s_max(ig)=0.
263
66778
       s_max(ig)=s2(ig)*log(max(n2(ig),1.))
264
!!fin jyg
265
     else
266
219494
       n2(ig)=0.
267
219494
       s_max(ig)=0.
268
     endif
269
   enddo
270
!   print *,'avant Calcul de Wmax '    !!jyg
271
272
288
   susqr2pi=su_cst*sqrt(2.*Rpi)
273
   reuler=exp(1.)
274
286560
   do ig=1,ngrid
275

286560
     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*reuler) ) then
276
63570
      w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
277
63570
      ale_bl_stat(ig)=0.5*w_max(ig)**2
278
     else
279
222702
      w_max(ig)=0.
280
222702
      ale_bl_stat(ig)=0.
281
     endif
282
   enddo
283
284
  ENDIF ! iflag_trig_bl
285
!  print *,'ENDIF  iflag_trig_bl'    !!jyg
286
287
!------------Closure------------------
288
289
288
  IF (iflag_clos_bl.ge.2) THEN
290
291
!-----Calcul de ALP_BL_STAT
292
  do ig=1,ngrid
293
  alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
294
  alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
295
 &                   (w0(ig)**2)
296
  alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
297
 &                    +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
298
    if (iflag_clos_bl.ge.2) then
299
    alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
300
 &                   (w0(ig)**2)
301
    else
302
    alp_bl_conv(ig)=0.
303
    endif
304
  alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
305
  enddo
306
307
!-----S�curit� ALP infinie
308
  do ig=1,ngrid
309
   if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
310
  enddo
311
312
  ENDIF ! (iflag_clos_bl.ge.2)
313
314
!!! fin nrlmd le 10/04/2012
315
316
!      print*,'avant calcul ale et alp'
317
!calcul de ALE et ALP pour la convection
318
286560
      alp_bl(:)=0.
319
286560
      ale_bl(:)=0.
320
!          print*,'ALE,ALP ,l,zw2(ig,l),ale_bl(ig),alp_bl(ig)'
321
11520
      do l=1,nlay
322
11176128
      do ig=1,ngrid
323
11164608
           alp_bl(ig)=max(alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
324
11175840
           ale_bl(ig)=max(ale_bl(ig),0.5*zw2(ig,l)**2)
325
!          print*,'ALE,ALP',l,zw2(ig,l),ale_bl(ig),alp_bl(ig)
326
      enddo
327
      enddo
328
329
! ale sec (max de wmax/2 sous la zone d'inhibition) dans
330
! le cas iflag_trig_bl=3
331

288
      IF (iflag_trig_bl==3) ale_bl(:)=0.5*wmax_sec(:)**2
332
333
!test:calcul de la ponderation des couches pour KE
334
!initialisations
335
336
286560
      fm_tot(:)=0.
337

11176128
      wght_th(:,:)=1.
338
286560
      lalim_conv(:)=lalim(:)
339
340
11520
      do k=1,nlay
341
11176128
         do ig=1,ngrid
342
11175840
            if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
343
         enddo
344
      enddo
345
346
! assez bizarre car, si on est dans la couche d'alim et que alim_star et
347
! plus petit que 1.e-10, on prend wght_th=1.
348
11520
      do k=1,nlay
349
11176128
         do ig=1,ngrid
350

11175840
            if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
351
341067
               wght_th(ig,k)=alim_star(ig,k)
352
            endif
353
         enddo
354
      enddo
355
356
!      print*,'apres wght_th'
357
!test pour prolonger la convection
358
286560
      do ig=1,ngrid
359
!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
360
286560
      if ((alim_star(ig,1).lt.1.e-10)) then
361
148935
      lalim_conv(ig)=1
362
148935
      wght_th(ig,1)=1.
363
!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
364
      endif
365
      enddo
366
367
!------------------------------------------------------------------------
368
! Modif CR/FH 20110310 : alp integree sur la verticale.
369
! Integrale verticale de ALP.
370
! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
371
! couches
372
!------------------------------------------------------------------------
373
374
286560
      alp_int(:)=0.
375
286560
      dp_int(:)=0.
376
11232
      do l=2,nlay
377
10889568
        do ig=1,ngrid
378
10889280
           if(l.LE.lmax(ig)) THEN
379
742187
           zdp=pplay(ig,l-1)-pplay(ig,l)
380
742187
           alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)*zdp
381
742187
           dp_int(ig)=dp_int(ig)+zdp
382
           endif
383
        enddo
384
      enddo
385
386
288
      if (iflag_coupl>=3 .and. iflag_coupl<=5) then
387
286560
      do ig=1,ngrid
388
!valeur integree de alp_bl * 0.5:
389
286560
        if (dp_int(ig)>0.) then
390
137471
        alp_bl(ig)=alp_int(ig)/dp_int(ig)
391
        endif
392
      enddo!
393
      endif
394
395
396
! Facteur multiplicatif sur alp_bl
397
286560
      alp_bl(:)=alp_bl_k*alp_bl(:)
398
399
!------------------------------------------------------------------------
400
401
402
403
288
      return
404
      end
405
END MODULE lmdz_thermcell_alp