source: LMDZ5/trunk/libf/phylmd/alpale_th.F90 @ 4273

Last change on this file since 4273 was 2730, checked in by jyg, 8 years ago

New input parameter: nbapp_wk = number of calls to
wake routines per day. If =0, call at every physic
time step.

File size: 11.1 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)))
65!
66!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
67!  JYG, 20160513 : Introduction of the Effective Lifting Power (ELP), which
68! takes into account the area (cv_feed_area) covered by thermals contributing
69! to each cumulonimbus.
70!   The use of ELP prevents singularities when the trigger probability tends to
71! zero. It is activated by iflag_clos_bl = 3.
72!   The ELP values are stored in the ALP_bl variable.
73!   
74!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
75!
76!!
77!!  The following 3 lines should be commented if one wants to activate the
78!! multiplication of no-trigger probabilities between calls to the convection
79!! scheme.
80!!
81             do i=1,klon
82                proba_notrig(i)=1.
83             enddo
84!!
85!!
86!---------------------------------------
87  IF (iflag_clos_bl .LT. 3) THEN
88!---------------------------------------
89!
90!      Original code (Nicolas Rochetin)
91!     --------------------------------
92
93    IF (first) THEN
94       random_notrig_max=1.
95       CALL getin_p('random_notrig_max',random_notrig_max)
96       first=.FALSE.
97    ENDIF
98          !cc nrlmd le 10/04/2012
99          !-----------Stochastic triggering-----------
100          if (iflag_trig_bl.ge.1) then
101             !
102             IF (prt_level .GE. 10) THEN
103                print *,'cin, ale_bl_stat, alp_bl_stat ', &
104                     cin, ale_bl_stat, alp_bl_stat
105             ENDIF
106
107
108             !----Initialisations
109             do i=1,klon
110!!jyg                proba_notrig(i)=1.
111                random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
112                if ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0.
113                if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
114                   tau_trig(i)=tau_trig_shallow
115                else
116                   tau_trig(i)=tau_trig_deep
117                endif
118             enddo
119             !
120             IF (prt_level .GE. 10) THEN
121                print *,'random_notrig, tau_trig ', &
122                     random_notrig, tau_trig
123                print *,'s_trig,s2,n2 ', &
124                     s_trig,s2,n2
125             ENDIF
126
127             !Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
128             IF (iflag_trig_bl.eq.1) then
129
130                !----Tirage al\'eatoire et calcul de ale_bl_trig
131                do i=1,klon
132                   if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
133                      proba_notrig(i)=proba_notrig(i)* &
134                         (1.-exp(-s_trig/s2(i)))**(n2(i)*dtime/tau_trig(i))
135                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
136                      if (random_notrig(i) .ge. proba_notrig(i)) then
137                         ale_bl_trig(i)=ale_bl_stat(i)
138                      else
139                         ale_bl_trig(i)=0.
140                      endif
141                   else
142!!jyg                      proba_notrig(i)=1.
143                      random_notrig(i)=0.
144                      ale_bl_trig(i)=0.
145                   endif
146                enddo
147
148             ELSE IF (iflag_trig_bl.ge.2) then
149
150                do i=1,klon
151                   if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) )  then
152                      proba_notrig(i)=proba_notrig(i)* &
153                         (1.-exp(-s_trig/s2(i)))**(n2(i)*dtime/tau_trig(i))
154                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
155                      if (random_notrig(i) .ge. proba_notrig(i)) then
156                         ale_bl_trig(i)=Ale_bl(i)
157                      else
158                         ale_bl_trig(i)=0.
159                      endif
160                   else
161!!jyg                      proba_notrig(i)=1.
162                      random_notrig(i)=0.
163                      ale_bl_trig(i)=0.
164                   endif
165                enddo
166
167             ENDIF
168
169             !
170             IF (prt_level .GE. 10) THEN
171                print *,'proba_notrig, ale_bl_trig ', &
172                     proba_notrig, ale_bl_trig
173             ENDIF
174
175          endif !(iflag_trig_bl)
176
177          !-----------Statistical closure-----------
178          if (iflag_clos_bl.eq.1) then
179
180             do i=1,klon
181                !CR: alp probabiliste
182                if (ale_bl_trig(i).gt.0.) then
183                   alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
184                endif
185             enddo
186
187          else if (iflag_clos_bl.eq.2) then
188
189             !CR: alp calculee dans thermcell_main
190             do i=1,klon
191                alp_bl(i)=alp_bl_stat(i)
192             enddo
193
194          else
195
196             alp_bl_stat(:)=0.
197
198          endif !(iflag_clos_bl)
199
200!
201!---------------------------------------
202  ELSEIF (iflag_clos_bl .EQ. 3) THEN  ! (iflag_clos_bl .LT. 3)
203!---------------------------------------
204!
205!      New code with Effective Lifting Power
206!     -------------------------------------
207    IF (first) THEN
208       cv_feed_area = 1.e10   ! m2
209       CALL getin_p('cv_feed_area', cv_feed_area)
210       first=.FALSE.
211    ENDIF
212
213          !-----------Stochastic triggering-----------
214     if (iflag_trig_bl.ge.1) then
215        !
216        IF (prt_level .GE. 10) THEN
217           print *,'cin, ale_bl_stat, alp_bl_stat ', &
218                cin, ale_bl_stat, alp_bl_stat
219        ENDIF
220
221        ! Use ale_bl_stat (Rochetin's code) or ale_bl (old code) according to
222        ! iflag_trig_bl value.
223        IF (iflag_trig_bl.eq.1) then         ! use ale_bl_stat (Rochetin computation)
224         do i=1,klon
225              ale_bl_ref(i)=ale_bl_stat(i)
226         enddo
227        ELSE IF (iflag_trig_bl.ge.2) then    ! use ale_bl (old computation)
228         do i=1,klon
229              ale_bl_ref(i)=Ale_bl(i)
230         enddo
231        ENDIF ! (iflag_trig_bl.eq.1)
232
233
234        !----Initializations and random number generation
235        do i=1,klon
236!!jyg           proba_notrig(i)=1.
237           random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
238           if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
239              tau_trig(i)=tau_trig_shallow
240           else
241              tau_trig(i)=tau_trig_deep
242           endif
243        enddo
244        !
245        IF (prt_level .GE. 10) THEN
246           print *,'random_notrig, tau_trig ', &
247                random_notrig, tau_trig
248           print *,'s_trig,s2,n2 ', &
249                s_trig,s2,n2
250        ENDIF
251
252        !----alp_bl computation
253        do i=1,klon
254           if ( (ale_bl_ref(i) .gt. abs(cin(i))+1.e-10) )  then
255              birth_number = n2(i)*exp(-s_trig/s2(i))
256              birth_rate(i) = birth_number/(tau_trig(i)*cell_area(i))
257              proba_notrig(i)=proba_notrig(i)*exp(-birth_number*dtime/tau_trig(i))
258              Alp_bl(i) = Alp_bl(i)* &
259                          umexp(-birth_number*cv_feed_area/cell_area(i))/ &
260                          umexp(-birth_number*dtime/tau_trig(i))*  &
261                          tau_trig(i)*cv_feed_area/(dtime*cell_area(i))
262          else
263!!jyg              proba_notrig(i)=1.
264              random_notrig(i)=0.
265              alp_bl(i)=0.
266           endif
267        enddo
268
269        !----ale_bl_trig computation
270         do i=1,klon
271           if (random_notrig(i) .ge. proba_notrig(i)) then
272              ale_bl_trig(i)=ale_bl_ref(i)
273           else
274              ale_bl_trig(i)=0.
275           endif
276         enddo
277
278        !
279        IF (prt_level .GE. 10) THEN
280           print *,'proba_notrig, ale_bl_trig ', &
281                proba_notrig, ale_bl_trig
282        ENDIF
283
284     endif !(iflag_trig_bl .ge. 1)
285
286!---------------------------------------
287  ENDIF ! (iflag_clos_bl .LT. 3)
288!---------------------------------------
289
290          IF (prt_level .GE. 10) THEN
291             print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
292          ENDIF
293
294          !cc fin nrlmd le 10/04/2012
295!
296          !IM/FH: 2011/02/23
297          ! Couplage Thermiques/Emanuel seulement si T<0
298          if (iflag_coupl==2) then
299             IF (prt_level .GE. 10) THEN
300                print*,'Couplage Thermiques/Emanuel seulement si T<0'
301             ENDIF
302             do i=1,klon
303                if (t_seri(i,lmax_th(i))>273.) then
304                   Ale_bl(i)=0.
305                endif
306             enddo
307    print *,'In order to run with iflag_coupl=2, you have to comment out the following stop'
308             STOP
309          endif
310   RETURN
311   END
312
Note: See TracBrowser for help on using the repository browser.