source: LMDZ6/branches/LMDZ-QUEST/libf/phylmd/alpale_th.F90

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

back to original (correct) formula in alpale_th

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