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

Last change on this file since 3709 was 3531, checked in by Laurent Fairhead, 5 years ago

Replaced STOP statements by a call to abort_physic in phylmd as per ticket #86
Still some work to be done in phylmd subdirectories

  • Property svn:keywords set to Id
File size: 12.2 KB
Line 
1!
2! $Id: alpale_th.F90 3531 2019-06-06 15:08:45Z adurocher $
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 "thermcell.h"
50
51! Local variables
52!----------------
53  INTEGER                                                    :: i
54  LOGICAL, SAVE                                              :: first = .TRUE.
55  REAL, SAVE                                                 :: random_notrig_max=1.
56  REAL, SAVE                                                 :: cv_feed_area
57  REAL                                                       :: birth_number
58  REAL, DIMENSION(klon)                                      :: ale_bl_ref
59  REAL, DIMENSION(klon)                                      :: tau_trig
60!
61    !$OMP THREADPRIVATE(random_notrig_max)
62    !$OMP THREADPRIVATE(cv_feed_area)
63    !$OMP THREADPRIVATE(first)
64!
65 REAL umexp  ! expression of (1.-exp(-x))/x valid for all x, especially when x->0
66 REAL x
67     CHARACTER (LEN=20) :: modname='alpale_th'
68     CHARACTER (LEN=80) :: abort_message
69     
70 umexp(x) = max(sign(1.,x-1.e-3),0.)*(1.-exp(-x))/max(x,1.e-3) + &
71            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! correct formula            (jyg)
72!!!            (1.-max(sign(1.,x-1.e-3),0.))*(-0.5*x*(1.-x/3.*(1.-0.25*x))) !!! bug introduced by mistake  (jyg)
73!!!            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! initial correct formula (jyg)
74!
75!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
76!  JYG, 20160513 : Introduction of the Effective Lifting Power (ELP), which
77! takes into account the area (cv_feed_area) covered by thermals contributing
78! to each cumulonimbus.
79!   The use of ELP prevents singularities when the trigger probability tends to
80! zero. It is activated by iflag_clos_bl = 3.
81!   The ELP values are stored in the ALP_bl variable.
82!   
83!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
84!
85!!
86!!  The following 3 lines should be commented if one wants to activate the
87!! multiplication of no-trigger probabilities between calls to the convection
88!! scheme.
89!!
90             do i=1,klon
91                proba_notrig(i)=1.
92             enddo
93!!
94!!
95!---------------------------------------
96  IF (iflag_clos_bl .LT. 3) THEN
97!---------------------------------------
98!
99!      Original code (Nicolas Rochetin)
100!     --------------------------------
101
102    IF (first) THEN
103       random_notrig_max=1.
104       CALL getin_p('random_notrig_max',random_notrig_max)
105       first=.FALSE.
106    ENDIF
107          !cc nrlmd le 10/04/2012
108          !-----------Stochastic triggering-----------
109          if (iflag_trig_bl.ge.1) then
110             !
111             IF (prt_level .GE. 10) THEN
112                WRITE(lunout,*)'cin, ale_bl_stat, alp_bl, alp_bl_stat ', &
113                     cin, ale_bl_stat, alp_bl, alp_bl_stat
114             ENDIF
115
116
117             !----Initialisations
118             do i=1,klon
119!!jyg                proba_notrig(i)=1.
120                random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
121                if ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0.
122                if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
123                   tau_trig(i)=tau_trig_shallow
124                else
125                   tau_trig(i)=tau_trig_deep
126                endif
127             enddo
128             !
129             IF (prt_level .GE. 10) THEN
130                WRITE(lunout,*)'random_notrig, tau_trig ', &
131                     random_notrig, tau_trig
132                WRITE(lunout,*)'s_trig,s2,n2 ', &
133                     s_trig,s2,n2
134             ENDIF
135
136             !Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
137             IF (iflag_trig_bl.eq.1) then
138
139                !----Tirage al\'eatoire et calcul de ale_bl_trig
140                do i=1,klon
141                   if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
142                      proba_notrig(i)=proba_notrig(i)* &
143                         (1.-exp(-s_trig/s2(i)))**(n2(i)*dtime/tau_trig(i))
144                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
145                      if (random_notrig(i) .ge. proba_notrig(i)) then
146                         ale_bl_trig(i)=ale_bl_stat(i)
147                      else
148                         ale_bl_trig(i)=0.
149                      endif
150                      birth_rate(i) = n2(i)*exp(-s_trig/s2(i))/(tau_trig(i)*cell_area(i))
151!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
152                   else
153!!jyg                      proba_notrig(i)=1.
154                      birth_rate(i) = 0.
155                      random_notrig(i)=0.
156                      ale_bl_trig(i)=0.
157                   endif
158                enddo
159
160             ELSE IF (iflag_trig_bl.ge.2) then
161
162                do i=1,klon
163                   if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) )  then
164                      proba_notrig(i)=proba_notrig(i)* &
165                         (1.-exp(-s_trig/s2(i)))**(n2(i)*dtime/tau_trig(i))
166                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
167                      if (random_notrig(i) .ge. proba_notrig(i)) then
168                         ale_bl_trig(i)=Ale_bl(i)
169                      else
170                         ale_bl_trig(i)=0.
171                      endif
172                      birth_rate(i) = n2(i)*exp(-s_trig/s2(i))/(tau_trig(i)*cell_area(i))
173!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
174                   else
175!!jyg                      proba_notrig(i)=1.
176                      birth_rate(i) = 0.
177                      random_notrig(i)=0.
178                      ale_bl_trig(i)=0.
179                   endif
180                enddo
181
182             ENDIF
183
184             !
185             IF (prt_level .GE. 10) THEN
186                WRITE(lunout,*)'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           WRITE(lunout,*)'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           WRITE(lunout,*)'random_notrig, tau_trig ', &
262                random_notrig, tau_trig
263           WRITE(lunout,*)'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              Alp_bl(i) = Alp_bl(i)* &
275                          umexp(-birth_number*cv_feed_area/cell_area(i))/ &
276                          umexp(-birth_number*dtime/tau_trig(i))*  &
277                          tau_trig(i)*cv_feed_area/(dtime*cell_area(i))
278          else
279!!jyg              proba_notrig(i)=1.
280              birth_rate(i)=0.
281              random_notrig(i)=0.
282              alp_bl(i)=0.
283           endif
284        enddo
285
286        !----ale_bl_trig computation
287         do i=1,klon
288           if (random_notrig(i) .ge. proba_notrig(i)) then
289              ale_bl_trig(i)=ale_bl_ref(i)
290           else
291              ale_bl_trig(i)=0.
292           endif
293         enddo
294
295        !
296        IF (prt_level .GE. 10) THEN
297           WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
298                proba_notrig, ale_bl_trig
299        ENDIF
300
301     endif !(iflag_trig_bl .ge. 1)
302
303!---------------------------------------
304  ENDIF ! (iflag_clos_bl .LT. 3)
305!---------------------------------------
306
307          IF (prt_level .GE. 10) THEN
308             WRITE(lunout,*)'alpale_th: ale_bl_trig, alp_bl_stat, birth_rate ', &
309                      ale_bl_trig(1), alp_bl_stat(1), birth_rate(1)
310          ENDIF
311
312          !cc fin nrlmd le 10/04/2012
313!
314          !IM/FH: 2011/02/23
315          ! Couplage Thermiques/Emanuel seulement si T<0
316          if (iflag_coupl==2) then
317             IF (prt_level .GE. 10) THEN
318                WRITE(lunout,*)'Couplage Thermiques/Emanuel seulement si T<0'
319             ENDIF
320             do i=1,klon
321                if (t_seri(i,lmax_th(i))>273.) then
322                   Ale_bl(i)=0.
323                endif
324             enddo
325!    print *,'In order to run with iflag_coupl=2, you have to comment out the following stop'
326!             STOP
327             abort_message='In order to run with iflag_coupl=2, you have to comment out the following abort'
328             CALL abort_physic(modname,abort_message,1)
329          endif
330   RETURN
331   END
332
Note: See TracBrowser for help on using the repository browser.