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

Last change on this file since 3208 was 3208, checked in by jyg, 7 years ago

Implementation of a first crude model of the
dynamic of wake population.

File size: 12.3 KB
Line 
1SUBROUTINE alpale_th ( dtime, lmax_th, t_seri, cell_area,  &
2                       cin, s2, n2,  &
3                       ale_bl_trig, ale_bl_stat, ale_bl,  &
4                       alp_bl, alp_bl_stat, &
5                       proba_notrig, random_notrig, birth_rate,  &
6                       q_alp)
7
8! **************************************************************
9! *
10! ALPALE_TH                                                    *
11! *
12! *
13! written by   : Jean-Yves Grandpeix, 11/05/2016              *
14! modified by :                                               *
15! **************************************************************
16
17  USE dimphy
18  USE ioipsl_getin_p_mod, ONLY : getin_p
19  USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
20!
21  IMPLICIT NONE
22
23!================================================================
24! Auteur(s)   : Jean-Yves Grandpeix, 11/05/2016
25! Objet : Contribution of the thermal scheme to Ale and Alp
26!================================================================
27
28! Input arguments
29!----------------
30  REAL, INTENT(IN)                                           :: dtime
31  REAL, DIMENSION(klon), INTENT(IN)                          :: cell_area
32  INTEGER, DIMENSION(klon), INTENT(IN)                       :: lmax_th
33  REAL, DIMENSION(klon,klev), INTENT(IN)                     :: t_seri
34  REAL, DIMENSION(klon), INTENT(IN)                          :: ale_bl_stat
35  REAL, DIMENSION(klon), INTENT(IN)                          :: cin
36  REAL, DIMENSION(klon), INTENT(IN)                          :: s2, n2
37                                                               
38  REAL, DIMENSION(klon), INTENT(INOUT)                       :: ale_bl_trig, ale_bl
39  REAL, DIMENSION(klon), INTENT(INOUT)                       :: alp_bl
40  REAL, DIMENSION(klon), INTENT(INOUT)                       :: alp_bl_stat
41  REAL, DIMENSION(klon), INTENT(INOUT)                       :: proba_notrig
42
43  REAL, DIMENSION(klon), INTENT(OUT)                         :: random_notrig
44
45  REAL, DIMENSION(klon), INTENT(OUT)                         :: birth_rate
46  REAL, DIMENSION(klon), INTENT(OUT)                         :: q_alp
47
48  include "thermcell.h"
49
50! Local variables
51!----------------
52  INTEGER                                                    :: i
53  LOGICAL, SAVE                                              :: first = .TRUE.
54  REAL, SAVE                                                 :: random_notrig_max=1.
55  REAL, SAVE                                                 :: cv_feed_area
56  REAL                                                       :: birth_number
57  REAL, DIMENSION(klon)                                      :: ale_bl_ref
58  REAL, DIMENSION(klon)                                      :: tau_trig
59!
60    !$OMP THREADPRIVATE(random_notrig_max)
61    !$OMP THREADPRIVATE(cv_feed_area)
62    !$OMP THREADPRIVATE(first)
63!
64 REAL umexp  ! expression of (1.-exp(-x))/x valid for all x, especially when x->0
65 REAL x
66 umexp(x) = max(sign(1.,x-1.e-3),0.)*(1.-exp(-x))/max(x,1.e-3) + &
67            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! correct formula            (jyg)
68!!!            (1.-max(sign(1.,x-1.e-3),0.))*(-0.5*x*(1.-x/3.*(1.-0.25*x))) !!! bug introduced by mistake  (jyg)
69!!!            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! initial correct formula (jyg)
70!
71!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
72!  JYG, 20160513 : Introduction of the Effective Lifting Power (ELP), which
73! takes into account the area (cv_feed_area) covered by thermals contributing
74! to each cumulonimbus.
75!   The use of ELP prevents singularities when the trigger probability tends to
76! zero. It is activated by iflag_clos_bl = 3.
77!   The ELP values are stored in the ALP_bl variable.
78!   
79!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
80!
81!!
82!!  The following 3 lines should be commented if one wants to activate the
83!! multiplication of no-trigger probabilities between calls to the convection
84!! scheme.
85!!
86             do i=1,klon
87                proba_notrig(i)=1.
88             enddo
89!!
90!!
91!---------------------------------------
92  IF (iflag_clos_bl .LT. 3) THEN
93!---------------------------------------
94!
95!      Original code (Nicolas Rochetin)
96!     --------------------------------
97
98    IF (first) THEN
99       random_notrig_max=1.
100       CALL getin_p('random_notrig_max',random_notrig_max)
101       first=.FALSE.
102    ENDIF
103          !cc nrlmd le 10/04/2012
104          !-----------Stochastic triggering-----------
105          if (iflag_trig_bl.ge.1) then
106             !
107             IF (prt_level .GE. 10) THEN
108                print *,'cin, ale_bl_stat, alp_bl, alp_bl_stat ', &
109                     cin, ale_bl_stat, alp_bl, alp_bl_stat
110             ENDIF
111
112
113             !----Initialisations
114             do i=1,klon
115!!jyg                proba_notrig(i)=1.
116                random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
117                if ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0.
118                if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
119                   tau_trig(i)=tau_trig_shallow
120                else
121                   tau_trig(i)=tau_trig_deep
122                endif
123             enddo
124             !
125             IF (prt_level .GE. 10) THEN
126                print *,'random_notrig, tau_trig ', &
127                     random_notrig, tau_trig
128                print *,'s_trig,s2,n2 ', &
129                     s_trig,s2,n2
130             ENDIF
131
132             !Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
133             IF (iflag_trig_bl.eq.1) then
134
135                !----Tirage al\'eatoire et calcul de ale_bl_trig
136                do i=1,klon
137                   if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
138                      proba_notrig(i)=proba_notrig(i)* &
139                         (1.-exp(-s_trig/s2(i)))**(n2(i)*dtime/tau_trig(i))
140                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
141                      if (random_notrig(i) .ge. proba_notrig(i)) then
142                         ale_bl_trig(i)=ale_bl_stat(i)
143                      else
144                         ale_bl_trig(i)=0.
145                      endif
146                      birth_rate(i) = n2(i)*exp(-s_trig/s2(i))/(tau_trig(i)*cell_area(i))
147!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
148                      q_alp(i) = alp_bl(i)/max(birth_rate(i),1.e-18)
149                   else
150!!jyg                      proba_notrig(i)=1.
151                      birth_rate(i) = 0.
152                      random_notrig(i)=0.
153                      ale_bl_trig(i)=0.
154                      q_alp(i) = 0.
155                   endif
156                enddo
157
158             ELSE IF (iflag_trig_bl.ge.2) then
159
160                do i=1,klon
161                   if ( (Ale_bl(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(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                      q_alp(i) = alp_bl(i)/max(birth_rate(i),1.e-18)
173                   else
174!!jyg                      proba_notrig(i)=1.
175                      birth_rate(i) = 0.
176                      random_notrig(i)=0.
177                      ale_bl_trig(i)=0.
178                      q_alp(i) = 0.
179                   endif
180                enddo
181
182             ENDIF
183
184             !
185             IF (prt_level .GE. 10) THEN
186                print *,'proba_notrig, ale_bl_trig ', &
187                     proba_notrig, ale_bl_trig
188             ENDIF
189
190          endif !(iflag_trig_bl)
191
192          !-----------Statistical closure-----------
193          if (iflag_clos_bl.eq.1) then
194
195             do i=1,klon
196                !CR: alp probabiliste
197                if (ale_bl_trig(i).gt.0.) then
198                   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           print *,'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           print *,'random_notrig, tau_trig ', &
262                random_notrig, tau_trig
263           print *,'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              q_alp(i) = alp_bl(i)/max(birth_rate(i),1.e-18)
275              Alp_bl(i) = Alp_bl(i)* &
276                          umexp(-birth_number*cv_feed_area/cell_area(i))/ &
277                          umexp(-birth_number*dtime/tau_trig(i))*  &
278                          tau_trig(i)*cv_feed_area/(dtime*cell_area(i))
279          else
280!!jyg              proba_notrig(i)=1.
281              birth_rate(i)=0.
282              random_notrig(i)=0.
283              alp_bl(i)=0.
284              q_alp(i) = 0.
285           endif
286        enddo
287
288        !----ale_bl_trig computation
289         do i=1,klon
290           if (random_notrig(i) .ge. proba_notrig(i)) then
291              ale_bl_trig(i)=ale_bl_ref(i)
292           else
293              ale_bl_trig(i)=0.
294           endif
295         enddo
296
297        !
298        IF (prt_level .GE. 10) THEN
299           print *,'proba_notrig, ale_bl_trig ', &
300                proba_notrig, ale_bl_trig
301        ENDIF
302
303     endif !(iflag_trig_bl .ge. 1)
304
305!---------------------------------------
306  ENDIF ! (iflag_clos_bl .LT. 3)
307!---------------------------------------
308
309          IF (prt_level .GE. 10) THEN
310             print *,'alpale_th: ale_bl_trig, alp_bl_stat, birth_rate, q_alp ', &
311                      ale_bl_trig(1), alp_bl_stat(1), birth_rate(1), q_alp(1)
312          ENDIF
313
314          !cc fin nrlmd le 10/04/2012
315!
316          !IM/FH: 2011/02/23
317          ! Couplage Thermiques/Emanuel seulement si T<0
318          if (iflag_coupl==2) then
319             IF (prt_level .GE. 10) THEN
320                print*,'Couplage Thermiques/Emanuel seulement si T<0'
321             ENDIF
322             do i=1,klon
323                if (t_seri(i,lmax_th(i))>273.) then
324                   Ale_bl(i)=0.
325                endif
326             enddo
327    print *,'In order to run with iflag_coupl=2, you have to comment out the following stop'
328             STOP
329          endif
330   RETURN
331   END
332
Note: See TracBrowser for help on using the repository browser.