source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/alpale_th.F90 @ 4613

Last change on this file since 4613 was 3331, checked in by acozic, 6 years ago

Add modification for isotopes

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