GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/alpale_th.F90 Lines: 39 112 34.8 %
Date: 2023-06-30 12:51:15 Branches: 29 90 32.2 %

Line Branch Exec Source
1
!
2
! $Id: alpale_th.F90 4089 2022-03-10 18:23:47Z fhourdin $
3
!
4
344299
SUBROUTINE alpale_th ( dtime, lmax_th, t_seri, cell_area,  &
5
                       cin, s2, n2,  &
6
                       ale_bl_trig, ale_bl_stat, ale_bl,  &
7
                       alp_bl, alp_bl_stat, &
8
288
                       proba_notrig, random_notrig, birth_rate)
9
10
! **************************************************************
11
! *
12
! ALPALE_TH                                                    *
13
! *
14
! *
15
! written by   : Jean-Yves Grandpeix, 11/05/2016              *
16
! modified by :                                               *
17
! **************************************************************
18
19
  USE dimphy
20
  USE ioipsl_getin_p_mod, ONLY : getin_p
21
  USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
22
!
23
  IMPLICIT NONE
24
25
!================================================================
26
! Auteur(s)   : Jean-Yves Grandpeix, 11/05/2016
27
! Objet : Contribution of the thermal scheme to Ale and Alp
28
!================================================================
29
30
! Input arguments
31
!----------------
32
  REAL, INTENT(IN)                                           :: dtime
33
  REAL, DIMENSION(klon), INTENT(IN)                          :: cell_area
34
  INTEGER, DIMENSION(klon), INTENT(IN)                       :: lmax_th
35
  REAL, DIMENSION(klon,klev), INTENT(IN)                     :: t_seri
36
  REAL, DIMENSION(klon), INTENT(IN)                          :: ale_bl_stat
37
  REAL, DIMENSION(klon), INTENT(IN)                          :: cin
38
  REAL, DIMENSION(klon), INTENT(IN)                          :: s2, n2
39
40
  REAL, DIMENSION(klon), INTENT(INOUT)                       :: ale_bl_trig, ale_bl
41
  REAL, DIMENSION(klon), INTENT(INOUT)                       :: alp_bl
42
  REAL, DIMENSION(klon), INTENT(INOUT)                       :: alp_bl_stat
43
  REAL, DIMENSION(klon), INTENT(INOUT)                       :: proba_notrig
44
45
  REAL, DIMENSION(klon), INTENT(OUT)                         :: random_notrig
46
47
  REAL, DIMENSION(klon), INTENT(OUT)                         :: birth_rate
48
49
  include "alpale.h"
50
51
! Local variables
52
!----------------
53
  INTEGER                                                    :: i
54
  LOGICAL, SAVE                                              :: first = .TRUE.
55
  REAL, SAVE                                                 :: random_notrig_max=1.
56
  REAL, SAVE                                                 :: cv_feed_area
57
  REAL                                                       :: birth_number
58
576
  REAL, DIMENSION(klon)                                      :: ale_bl_ref
59
288
  REAL, DIMENSION(klon)                                      :: tau_trig
60
!
61
    !$OMP THREADPRIVATE(random_notrig_max)
62
    !$OMP THREADPRIVATE(cv_feed_area)
63
    !$OMP THREADPRIVATE(first)
64
!
65
 REAL umexp  ! expression of (1.-exp(-x))/x valid for all x, especially when x->0
66
 REAL x
67
     CHARACTER (LEN=20) :: modname='alpale_th'
68
     CHARACTER (LEN=80) :: abort_message
69
70
 umexp(x) = max(sign(1.,x-1.e-3),0.)*(1.-exp(-x))/max(x,1.e-3) + &
71
            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! correct formula            (jyg)
72
!!!            (1.-max(sign(1.,x-1.e-3),0.))*(-0.5*x*(1.-x/3.*(1.-0.25*x))) !!! bug introduced by mistake  (jyg)
73
!!!            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! initial correct formula (jyg)
74
!
75
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
76
!  JYG, 20160513 : Introduction of the Effective Lifting Power (ELP), which
77
! takes into account the area (cv_feed_area) covered by thermals contributing
78
! to each cumulonimbus.
79
!   The use of ELP prevents singularities when the trigger probability tends to
80
! zero. It is activated by iflag_clos_bl = 3.
81
!   The ELP values are stored in the ALP_bl variable.
82
!
83
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
84
!
85
!!
86
!!  The following 3 lines should be commented if one wants to activate the
87
!! multiplication of no-trigger probabilities between calls to the convection
88
!! scheme.
89
!!
90
286560
             do i=1,klon
91
286560
                proba_notrig(i)=1.
92
             enddo
93
!!
94
!!
95
!---------------------------------------
96
288
  IF (iflag_clos_bl .LT. 3) THEN
97
!---------------------------------------
98
!
99
!      Original code (Nicolas Rochetin)
100
!     --------------------------------
101
102
288
    IF (first) THEN
103
1
       random_notrig_max=1.
104
1
       CALL getin_p('random_notrig_max',random_notrig_max)
105
1
       first=.FALSE.
106
    ENDIF
107
          !cc nrlmd le 10/04/2012
108
          !-----------Stochastic triggering-----------
109
288
          if (iflag_trig_bl.ge.1) then
110
             !
111
288
             IF (prt_level .GE. 10) THEN
112
                WRITE(lunout,*)'cin, ale_bl_stat, alp_bl, alp_bl_stat ', &
113
                     cin, ale_bl_stat, alp_bl, alp_bl_stat
114
             ENDIF
115
116
117
             !----Initialisations
118
286560
             do i=1,klon
119
!!jyg                proba_notrig(i)=1.
120
286272
                random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
121
286272
                if ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0.
122
286560
                if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
123
283181
                   tau_trig(i)=tau_trig_shallow
124
                else
125
3091
                   tau_trig(i)=tau_trig_deep
126
                endif
127
             enddo
128
             !
129
288
             IF (prt_level .GE. 10) THEN
130
                WRITE(lunout,*)'random_notrig, tau_trig ', &
131
                     random_notrig, tau_trig
132
                WRITE(lunout,*)'s_trig,s2,n2 ', &
133
                     s_trig,s2,n2
134
             ENDIF
135
136
             !Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
137
288
             IF (iflag_trig_bl.eq.1) then
138
139
                !----Tirage al\'eatoire et calcul de ale_bl_trig
140
286560
                do i=1,klon
141
286560
                   if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
142
                      proba_notrig(i)=proba_notrig(i)* &
143
56299
                         (1.-exp(-s_trig/s2(i)))**(n2(i)*dtime/tau_trig(i))
144
                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
145
56299
                      if (random_notrig(i) .ge. proba_notrig(i)) then
146
3115
                         ale_bl_trig(i)=ale_bl_stat(i)
147
                      else
148
53184
                         ale_bl_trig(i)=0.
149
                      endif
150
56299
                      birth_rate(i) = n2(i)*exp(-s_trig/s2(i))/(tau_trig(i)*cell_area(i))
151
!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
152
                   else
153
!!jyg                      proba_notrig(i)=1.
154
229973
                      birth_rate(i) = 0.
155
229973
                      random_notrig(i)=0.
156
229973
                      ale_bl_trig(i)=0.
157
                   endif
158
                enddo
159
160
             ELSE IF (iflag_trig_bl.ge.2) then
161
162
                do i=1,klon
163
                   if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) )  then
164
                      proba_notrig(i)=proba_notrig(i)* &
165
                         (1.-exp(-s_trig/s2(i)))**(n2(i)*dtime/tau_trig(i))
166
                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
167
                      if (random_notrig(i) .ge. proba_notrig(i)) then
168
                         ale_bl_trig(i)=Ale_bl(i)
169
                      else
170
                         ale_bl_trig(i)=0.
171
                      endif
172
                      birth_rate(i) = n2(i)*exp(-s_trig/s2(i))/(tau_trig(i)*cell_area(i))
173
!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
174
                   else
175
!!jyg                      proba_notrig(i)=1.
176
                      birth_rate(i) = 0.
177
                      random_notrig(i)=0.
178
                      ale_bl_trig(i)=0.
179
                   endif
180
                enddo
181
182
             ENDIF
183
184
             !
185
288
             IF (prt_level .GE. 10) THEN
186
                WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
187
                     proba_notrig, ale_bl_trig
188
             ENDIF
189
190
          endif !(iflag_trig_bl)
191
192
          !-----------Statistical closure-----------
193
288
          if (iflag_clos_bl.eq.1) then
194
195
286560
             do i=1,klon
196
                !CR: alp probabiliste
197
286560
                if (ale_bl_trig(i).gt.0.) then
198
3115
                   alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
199
                endif
200
             enddo
201
202
          else if (iflag_clos_bl.eq.2) then
203
204
             !CR: alp calculee dans thermcell_main
205
             do i=1,klon
206
                alp_bl(i)=alp_bl_stat(i)
207
             enddo
208
209
          else
210
211
             alp_bl_stat(:)=0.
212
213
          endif !(iflag_clos_bl)
214
215
!
216
!---------------------------------------
217
  ELSEIF (iflag_clos_bl .EQ. 3) THEN  ! (iflag_clos_bl .LT. 3)
218
!---------------------------------------
219
!
220
!      New code with Effective Lifting Power
221
!     -------------------------------------
222
    IF (first) THEN
223
       cv_feed_area = 1.e10   ! m2
224
       CALL getin_p('cv_feed_area', cv_feed_area)
225
       first=.FALSE.
226
    ENDIF
227
228
          !-----------Stochastic triggering-----------
229
     if (iflag_trig_bl.ge.1) then
230
        !
231
        IF (prt_level .GE. 10) THEN
232
           WRITE(lunout,*)'cin, ale_bl_stat, alp_bl_stat ', &
233
                cin, ale_bl_stat, alp_bl_stat
234
        ENDIF
235
236
        ! Use ale_bl_stat (Rochetin's code) or ale_bl (old code) according to
237
        ! iflag_trig_bl value.
238
        IF (iflag_trig_bl.eq.1) then         ! use ale_bl_stat (Rochetin computation)
239
         do i=1,klon
240
              ale_bl_ref(i)=ale_bl_stat(i)
241
         enddo
242
        ELSE IF (iflag_trig_bl.ge.2) then    ! use ale_bl (old computation)
243
         do i=1,klon
244
              ale_bl_ref(i)=Ale_bl(i)
245
         enddo
246
        ENDIF ! (iflag_trig_bl.eq.1)
247
248
249
        !----Initializations and random number generation
250
        do i=1,klon
251
!!jyg           proba_notrig(i)=1.
252
           random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
253
           if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
254
              tau_trig(i)=tau_trig_shallow
255
           else
256
              tau_trig(i)=tau_trig_deep
257
           endif
258
        enddo
259
        !
260
        IF (prt_level .GE. 10) THEN
261
           WRITE(lunout,*)'random_notrig, tau_trig ', &
262
                random_notrig, tau_trig
263
           WRITE(lunout,*)'s_trig,s2,n2 ', &
264
                s_trig,s2,n2
265
        ENDIF
266
267
        !----alp_bl computation
268
        do i=1,klon
269
           if ( (ale_bl_ref(i) .gt. abs(cin(i))+1.e-10) )  then
270
              birth_number = n2(i)*exp(-s_trig/s2(i))
271
              birth_rate(i) = birth_number/(tau_trig(i)*cell_area(i))
272
!!!              birth_rate(i) = max(birth_rate(i),1.e-18)
273
              proba_notrig(i)=proba_notrig(i)*exp(-birth_number*dtime/tau_trig(i))
274
              Alp_bl(i) = Alp_bl(i)* &
275
                          umexp(-birth_number*cv_feed_area/cell_area(i))/ &
276
                          umexp(-birth_number*dtime/tau_trig(i))*  &
277
                          tau_trig(i)*cv_feed_area/(dtime*cell_area(i))
278
          else
279
!!jyg              proba_notrig(i)=1.
280
              birth_rate(i)=0.
281
              random_notrig(i)=0.
282
              alp_bl(i)=0.
283
           endif
284
        enddo
285
286
        !----ale_bl_trig computation
287
         do i=1,klon
288
           if (random_notrig(i) .ge. proba_notrig(i)) then
289
              ale_bl_trig(i)=ale_bl_ref(i)
290
           else
291
              ale_bl_trig(i)=0.
292
           endif
293
         enddo
294
295
        !
296
        IF (prt_level .GE. 10) THEN
297
           WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
298
                proba_notrig, ale_bl_trig
299
        ENDIF
300
301
     endif !(iflag_trig_bl .ge. 1)
302
303
!---------------------------------------
304
  ENDIF ! (iflag_clos_bl .LT. 3)
305
!---------------------------------------
306
307
288
          IF (prt_level .GE. 10) THEN
308
             WRITE(lunout,*)'alpale_th: ale_bl_trig, alp_bl_stat, birth_rate ', &
309
                      ale_bl_trig(1), alp_bl_stat(1), birth_rate(1)
310
          ENDIF
311
312
          !cc fin nrlmd le 10/04/2012
313
!
314
          !IM/FH: 2011/02/23
315
          ! Couplage Thermiques/Emanuel seulement si T<0
316
288
          if (iflag_coupl==2) then
317
             IF (prt_level .GE. 10) THEN
318
                WRITE(lunout,*)'Couplage Thermiques/Emanuel seulement si T<0'
319
             ENDIF
320
             do i=1,klon
321
                if (t_seri(i,lmax_th(i))>273.) then
322
                   Ale_bl(i)=0.
323
                endif
324
             enddo
325
!    print *,'In order to run with iflag_coupl=2, you have to comment out the following stop'
326
!             STOP
327
             abort_message='In order to run with iflag_coupl=2, you have to comment out the following abort'
328
             CALL abort_physic(modname,abort_message,1)
329
          endif
330
288
   RETURN
331
   END
332