source: LMDZ6/trunk/libf/phylmd/alpale_th.F90 @ 4827

Last change on this file since 4827 was 4827, checked in by jyg, 3 months ago

Improved comment for multiply_proba_notrig

  • Property svn:keywords set to Id
File size: 13.2 KB
Line 
1!
2! $Id: alpale_th.F90 4827 2024-02-16 13:59:33Z jyg $
3!
4SUBROUTINE 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                       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  LOGICAL, SAVE                                              :: multiply_proba_notrig = .FALSE.
56  REAL, SAVE                                                 :: random_notrig_max=1.
57  REAL, SAVE                                                 :: cv_feed_area
58  REAL                                                       :: birth_number
59  REAL, DIMENSION(klon)                                      :: ale_bl_ref
60  REAL, DIMENSION(klon)                                      :: tau_trig
61!
62    !$OMP THREADPRIVATE(multiply_proba_notrig)
63    !$OMP THREADPRIVATE(random_notrig_max)
64    !$OMP THREADPRIVATE(cv_feed_area)
65    !$OMP THREADPRIVATE(first)
66!
67 REAL umexp  ! expression of (1.-exp(-x))/x valid for all x, especially when x->0
68 REAL x
69!
70     CHARACTER (LEN=20) :: modname='alpale_th'
71     CHARACTER (LEN=80) :: abort_message
72     
73 umexp(x) = max(sign(1.,x-1.e-3),0.)*(1.-exp(-x))/max(x,1.e-3) + &
74            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! correct formula            (jyg)
75!!!            (1.-max(sign(1.,x-1.e-3),0.))*(-0.5*x*(1.-x/3.*(1.-0.25*x))) !!! bug introduced by mistake  (jyg)
76!!!            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! initial correct formula (jyg)
77!
78!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
79!  JYG, 20160513 : Introduction of the Effective Lifting Power (ELP), which
80! takes into account the area (cv_feed_area) covered by thermals contributing
81! to each cumulonimbus.
82!   The use of ELP prevents singularities when the trigger probability tends to
83! zero. It is activated by iflag_clos_bl = 3.
84!   The ELP values are stored in the ALP_bl variable.
85!   
86!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87!
88
89    IF (first) THEN
90      CALL getin_p('multiply_proba_notrig',multiply_proba_notrig)
91      IF (iflag_clos_bl .LT. 3) THEN
92         random_notrig_max=1.
93         CALL getin_p('random_notrig_max',random_notrig_max)
94      ELSEIF (iflag_clos_bl .EQ. 3) THEN  ! (iflag_clos_bl .LT. 3)
95         cv_feed_area = 1.e10   ! m2
96         CALL getin_p('cv_feed_area', cv_feed_area)
97      ENDIF  !! (iflag_clos_bl .LT. 3)
98      first=.FALSE.
99    ENDIF
100
101!!
102!!   Control of the multiplication of no-trigger probabilities between calls
103!! to the convection scheme. If multiply_proba_notrig is .false., then
104!! proba_notrig is set to 1 at each call to alpale_th, so that only the last call
105!! plays a role in the triggering of convection. If it is .true., then propa_notrig
106!! is saved between calls to convection and is reset to 1 only after calling the
107!! convection scheme.
108!!    For instance, if the probability of no_trigger is 0.9 at each call, and if
109!! there are 3 calls to alpale_th between calls to the convection scheme, then the
110!! probability of triggering convection will be 0.1 (= 1.-0.9) if
111!! multiply_proba_notrig is .false. and 0.271 (= 1.-0.9^3) if multiply_proba_notrig
112!! is .true.
113!!
114    IF (.NOT.multiply_proba_notrig) THEN
115             DO i=1,klon
116                proba_notrig(i)=1.
117             ENDDO
118    ENDIF  !! (.NOT.multiply_proba_notrig)
119!!
120!!
121!---------------------------------------
122  IF (iflag_clos_bl .LT. 3) THEN
123!---------------------------------------
124!
125!      Original code (Nicolas Rochetin)
126!     --------------------------------
127          !cc nrlmd le 10/04/2012
128          !-----------Stochastic triggering-----------
129          if (iflag_trig_bl.ge.1) then
130             !
131             IF (prt_level .GE. 10) THEN
132                WRITE(lunout,*)'cin, ale_bl_stat, alp_bl, alp_bl_stat ', &
133                     cin, ale_bl_stat, alp_bl, alp_bl_stat
134             ENDIF
135
136
137             !----Initialisations
138             do i=1,klon
139!!jyg                proba_notrig(i)=1.
140                random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
141                if ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0.
142                if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
143                   tau_trig(i)=tau_trig_shallow
144                else
145                   tau_trig(i)=tau_trig_deep
146                endif
147             enddo
148             !
149             IF (prt_level .GE. 10) THEN
150                WRITE(lunout,*)'random_notrig, tau_trig ', &
151                     random_notrig, tau_trig
152                WRITE(lunout,*)'s_trig,s2,n2 ', &
153                     s_trig,s2,n2
154             ENDIF
155
156             !Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
157             IF (iflag_trig_bl.eq.1) then
158
159                !----Tirage al\'eatoire et calcul de ale_bl_trig
160                do i=1,klon
161                   if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
162                      proba_notrig(i)=proba_notrig(i)* &
163                         (1.-exp(-s_trig/s2(i)))**(n2(i)*dtime/tau_trig(i))
164                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
165                      if (random_notrig(i) .ge. proba_notrig(i)) then
166                         ale_bl_trig(i)=ale_bl_stat(i)
167                      else
168                         ale_bl_trig(i)=0.
169                      endif
170                      birth_rate(i) = n2(i)*exp(-s_trig/s2(i))/(tau_trig(i)*cell_area(i))
171!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
172                   else
173!!jyg                      proba_notrig(i)=1.
174                      birth_rate(i) = 0.
175                      random_notrig(i)=0.
176                      ale_bl_trig(i)=0.
177                   endif
178                enddo
179
180             ELSE IF (iflag_trig_bl.ge.2) then
181
182                do i=1,klon
183                   if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) )  then
184                      proba_notrig(i)=proba_notrig(i)* &
185                         (1.-exp(-s_trig/s2(i)))**(n2(i)*dtime/tau_trig(i))
186                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
187                      if (random_notrig(i) .ge. proba_notrig(i)) then
188                         ale_bl_trig(i)=Ale_bl(i)
189                      else
190                         ale_bl_trig(i)=0.
191                      endif
192                      birth_rate(i) = n2(i)*exp(-s_trig/s2(i))/(tau_trig(i)*cell_area(i))
193!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
194                   else
195!!jyg                      proba_notrig(i)=1.
196                      birth_rate(i) = 0.
197                      random_notrig(i)=0.
198                      ale_bl_trig(i)=0.
199                   endif
200                enddo
201
202             ENDIF
203
204             !
205             IF (prt_level .GE. 10) THEN
206                WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
207                     proba_notrig, ale_bl_trig
208             ENDIF
209
210          endif !(iflag_trig_bl)
211
212          !-----------Statistical closure-----------
213          if (iflag_clos_bl.eq.1) then
214
215             do i=1,klon
216                !CR: alp probabiliste
217                if (ale_bl_trig(i).gt.0.) then
218                   alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
219                endif
220             enddo
221
222          else if (iflag_clos_bl.eq.2) then
223
224             !CR: alp calculee dans thermcell_main
225             do i=1,klon
226                alp_bl(i)=alp_bl_stat(i)
227             enddo
228
229          else
230
231             alp_bl_stat(:)=0.
232
233          endif !(iflag_clos_bl)
234
235!
236!---------------------------------------
237  ELSEIF (iflag_clos_bl .EQ. 3) THEN  ! (iflag_clos_bl .LT. 3)
238!---------------------------------------
239!
240!      New code with Effective Lifting Power
241!     -------------------------------------
242
243          !-----------Stochastic triggering-----------
244     if (iflag_trig_bl.ge.1) then
245        !
246        IF (prt_level .GE. 10) THEN
247           WRITE(lunout,*)'cin, ale_bl_stat, alp_bl_stat ', &
248                cin, ale_bl_stat, alp_bl_stat
249        ENDIF
250
251        ! Use ale_bl_stat (Rochetin's code) or ale_bl (old code) according to
252        ! iflag_trig_bl value.
253        IF (iflag_trig_bl.eq.1) then         ! use ale_bl_stat (Rochetin computation)
254         do i=1,klon
255              ale_bl_ref(i)=ale_bl_stat(i)
256         enddo
257        ELSE IF (iflag_trig_bl.ge.2) then    ! use ale_bl (old computation)
258         do i=1,klon
259              ale_bl_ref(i)=Ale_bl(i)
260         enddo
261        ENDIF ! (iflag_trig_bl.eq.1)
262
263
264        !----Initializations and random number generation
265        do i=1,klon
266!!jyg           proba_notrig(i)=1.
267           random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
268           if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
269              tau_trig(i)=tau_trig_shallow
270           else
271              tau_trig(i)=tau_trig_deep
272           endif
273        enddo
274        !
275        IF (prt_level .GE. 10) THEN
276           WRITE(lunout,*)'random_notrig, tau_trig ', &
277                random_notrig, tau_trig
278           WRITE(lunout,*)'s_trig,s2,n2 ', &
279                s_trig,s2,n2
280        ENDIF
281
282        !----alp_bl computation
283        do i=1,klon
284           if ( (ale_bl_ref(i) .gt. abs(cin(i))+1.e-10) )  then
285              birth_number = n2(i)*exp(-s_trig/s2(i))
286              birth_rate(i) = birth_number/(tau_trig(i)*cell_area(i))
287!!!              birth_rate(i) = max(birth_rate(i),1.e-18)
288              proba_notrig(i)=proba_notrig(i)*exp(-birth_number*dtime/tau_trig(i))
289              Alp_bl(i) = Alp_bl(i)* &
290                          umexp(-birth_number*cv_feed_area/cell_area(i))/ &
291                          umexp(-birth_number*dtime/tau_trig(i))*  &
292                          tau_trig(i)*cv_feed_area/(dtime*cell_area(i))
293          else
294!!jyg              proba_notrig(i)=1.
295              birth_rate(i)=0.
296              random_notrig(i)=0.
297              alp_bl(i)=0.
298           endif
299        enddo
300
301        !----ale_bl_trig computation
302         do i=1,klon
303           if (random_notrig(i) .ge. proba_notrig(i)) then
304              ale_bl_trig(i)=ale_bl_ref(i)
305           else
306              ale_bl_trig(i)=0.
307           endif
308         enddo
309
310        !
311        IF (prt_level .GE. 10) THEN
312           WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
313                proba_notrig, ale_bl_trig
314        ENDIF
315
316     endif !(iflag_trig_bl .ge. 1)
317
318!---------------------------------------
319  ENDIF ! (iflag_clos_bl .LT. 3)
320!---------------------------------------
321
322          IF (prt_level .GE. 10) THEN
323             WRITE(lunout,*)'alpale_th: ale_bl_trig, alp_bl_stat, birth_rate ', &
324                      ale_bl_trig(1), alp_bl_stat(1), birth_rate(1)
325          ENDIF
326
327          !cc fin nrlmd le 10/04/2012
328!
329          !IM/FH: 2011/02/23
330          ! Couplage Thermiques/Emanuel seulement si T<0
331          if (iflag_coupl==2) then
332             IF (prt_level .GE. 10) THEN
333                WRITE(lunout,*)'Couplage Thermiques/Emanuel seulement si T<0'
334             ENDIF
335             do i=1,klon
336                if (t_seri(i,lmax_th(i))>273.) then
337                   Ale_bl(i)=0.
338                endif
339             enddo
340!    print *,'In order to run with iflag_coupl=2, you have to comment out the following stop'
341!             STOP
342             abort_message='In order to run with iflag_coupl=2, you have to comment out the following abort'
343             CALL abort_physic(modname,abort_message,1)
344          endif
345   RETURN
346   END
347
Note: See TracBrowser for help on using the repository browser.