source: LMDZ6/trunk/libf/phylmd/alpale_th.f90 @ 5308

Last change on this file since 5308 was 5284, checked in by abarral, 4 days ago

Turn alpale.h alpale.f90 YOETHF.h into modules

  • Property svn:keywords set to Id
File size: 13.2 KB
Line 
1!
2! $Id: alpale_th.f90 5284 2024-10-28 13:12:55Z fairhead $
3!
4SUBROUTINE alpale_th ( dtime, lmax_th, t_seri, cell_area,  &
5                       cin, s2, n2, strig,  &
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  USE alpale_mod
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, strig
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! Local variables
50!----------------
51  INTEGER                                                    :: i
52  LOGICAL, SAVE                                              :: first = .TRUE.
53  LOGICAL, SAVE                                              :: multiply_proba_notrig = .FALSE.
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(multiply_proba_notrig)
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!
68     CHARACTER (LEN=20) :: modname='alpale_th'
69     CHARACTER (LEN=80) :: abort_message
70     
71 umexp(x) = max(sign(1.,x-1.e-3),0.)*(1.-exp(-x))/max(x,1.e-3) + &
72            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! correct formula            (jyg)
73!!!            (1.-max(sign(1.,x-1.e-3),0.))*(-0.5*x*(1.-x/3.*(1.-0.25*x))) !!! bug introduced by mistake  (jyg)
74!!!            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! initial correct formula (jyg)
75!
76!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
77!  JYG, 20160513 : Introduction of the Effective Lifting Power (ELP), which
78! takes into account the area (cv_feed_area) covered by thermals contributing
79! to each cumulonimbus.
80!   The use of ELP prevents singularities when the trigger probability tends to
81! zero. It is activated by iflag_clos_bl = 3.
82!   The ELP values are stored in the ALP_bl variable.
83!   
84!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
85!
86
87    IF (first) THEN
88      CALL getin_p('multiply_proba_notrig',multiply_proba_notrig)
89      IF (iflag_clos_bl .LT. 3) THEN
90         random_notrig_max=1.
91         CALL getin_p('random_notrig_max',random_notrig_max)
92      ELSEIF (iflag_clos_bl .EQ. 3) THEN  ! (iflag_clos_bl .LT. 3)
93         cv_feed_area = 1.e10   ! m2
94         CALL getin_p('cv_feed_area', cv_feed_area)
95      ENDIF  !! (iflag_clos_bl .LT. 3)
96      first=.FALSE.
97    ENDIF
98
99!!
100!!   Control of the multiplication of no-trigger probabilities between calls
101!! to the convection scheme. If multiply_proba_notrig is .false., then
102!! proba_notrig is set to 1 at each call to alpale_th, so that only the last call
103!! plays a role in the triggering of convection. If it is .true., then propa_notrig
104!! is saved between calls to convection and is reset to 1 only after calling the
105!! convection scheme.
106!!    For instance, if the probability of no_trigger is 0.9 at each call, and if
107!! there are 3 calls to alpale_th between calls to the convection scheme, then the
108!! probability of triggering convection will be 0.1 (= 1.-0.9) if
109!! multiply_proba_notrig is .false. and 0.271 (= 1.-0.9^3) if multiply_proba_notrig
110!! is .true.
111!!
112    IF (.NOT.multiply_proba_notrig) THEN
113             DO i=1,klon
114                proba_notrig(i)=1.
115             ENDDO
116    ENDIF  !! (.NOT.multiply_proba_notrig)
117!!
118!!
119!---------------------------------------
120  IF (iflag_clos_bl .LT. 3) THEN
121!---------------------------------------
122!
123!      Original code (Nicolas Rochetin)
124!     --------------------------------
125          !cc nrlmd le 10/04/2012
126          !-----------Stochastic triggering-----------
127          if (iflag_trig_bl.ge.1) then
128             !
129             IF (prt_level .GE. 10) THEN
130                WRITE(lunout,*)'cin, ale_bl_stat, alp_bl, alp_bl_stat ', &
131                     cin, ale_bl_stat, alp_bl, alp_bl_stat
132             ENDIF
133
134
135             !----Initialisations
136             do i=1,klon
137!!jyg                proba_notrig(i)=1.
138                random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
139                if ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0.
140                if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
141                   tau_trig(i)=tau_trig_shallow
142                else
143                   tau_trig(i)=tau_trig_deep
144                endif
145             enddo
146             !
147             IF (prt_level .GE. 10) THEN
148                WRITE(lunout,*)'random_notrig, tau_trig ', &
149                     random_notrig, tau_trig
150                WRITE(lunout,*)'s_trig,s2,n2 ', &
151                     s_trig,s2,n2
152             ENDIF
153
154             !Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
155             IF (iflag_trig_bl.eq.1) then
156
157                !----Tirage al\'eatoire et calcul de ale_bl_trig
158                do i=1,klon
159                   if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
160                      proba_notrig(i)=proba_notrig(i)* &
161                         (1.-exp(-strig(i)/s2(i)))**(n2(i)*dtime/tau_trig(i))
162                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
163                      if (random_notrig(i) .ge. proba_notrig(i)) then
164                         ale_bl_trig(i)=ale_bl_stat(i)
165                      else
166                         ale_bl_trig(i)=0.
167                      endif
168                      birth_rate(i) = n2(i)*exp(-strig(i)/s2(i))/(tau_trig(i)*cell_area(i))
169!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
170                   else
171!!jyg                      proba_notrig(i)=1.
172                      birth_rate(i) = 0.
173                      random_notrig(i)=0.
174                      ale_bl_trig(i)=0.
175                   endif
176                enddo
177
178             ELSE IF (iflag_trig_bl.ge.2) then
179
180                do i=1,klon
181                   if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) )  then
182                      proba_notrig(i)=proba_notrig(i)* &
183                         (1.-exp(-strig(i)/s2(i)))**(n2(i)*dtime/tau_trig(i))
184                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
185                      if (random_notrig(i) .ge. proba_notrig(i)) then
186                         ale_bl_trig(i)=Ale_bl(i)
187                      else
188                         ale_bl_trig(i)=0.
189                      endif
190                      birth_rate(i) = n2(i)*exp(-strig(i)/s2(i))/(tau_trig(i)*cell_area(i))
191!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
192                   else
193!!jyg                      proba_notrig(i)=1.
194                      birth_rate(i) = 0.
195                      random_notrig(i)=0.
196                      ale_bl_trig(i)=0.
197                   endif
198                enddo
199
200             ENDIF
201
202             !
203             IF (prt_level .GE. 10) THEN
204                WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
205                     proba_notrig, ale_bl_trig
206             ENDIF
207
208          endif !(iflag_trig_bl)
209
210          !-----------Statistical closure-----------
211          if (iflag_clos_bl.eq.1) then
212
213             do i=1,klon
214                !CR: alp probabiliste
215                if (ale_bl_trig(i).gt.0.) then
216                   alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
217                endif
218             enddo
219
220          else if (iflag_clos_bl.eq.2) then
221
222             !CR: alp calculee dans thermcell_main
223             do i=1,klon
224                alp_bl(i)=alp_bl_stat(i)
225             enddo
226
227          else
228
229             alp_bl_stat(:)=0.
230
231          endif !(iflag_clos_bl)
232
233!
234!---------------------------------------
235  ELSEIF (iflag_clos_bl .EQ. 3) THEN  ! (iflag_clos_bl .LT. 3)
236!---------------------------------------
237!
238!      New code with Effective Lifting Power
239!     -------------------------------------
240
241          !-----------Stochastic triggering-----------
242     if (iflag_trig_bl.ge.1) then
243        !
244        IF (prt_level .GE. 10) THEN
245           WRITE(lunout,*)'cin, ale_bl_stat, alp_bl_stat ', &
246                cin, ale_bl_stat, alp_bl_stat
247        ENDIF
248
249        ! Use ale_bl_stat (Rochetin's code) or ale_bl (old code) according to
250        ! iflag_trig_bl value.
251        IF (iflag_trig_bl.eq.1) then         ! use ale_bl_stat (Rochetin computation)
252         do i=1,klon
253              ale_bl_ref(i)=ale_bl_stat(i)
254         enddo
255        ELSE IF (iflag_trig_bl.ge.2) then    ! use ale_bl (old computation)
256         do i=1,klon
257              ale_bl_ref(i)=Ale_bl(i)
258         enddo
259        ENDIF ! (iflag_trig_bl.eq.1)
260
261
262        !----Initializations and random number generation
263        do i=1,klon
264!!jyg           proba_notrig(i)=1.
265           random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
266           if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
267              tau_trig(i)=tau_trig_shallow
268           else
269              tau_trig(i)=tau_trig_deep
270           endif
271        enddo
272        !
273        IF (prt_level .GE. 10) THEN
274           WRITE(lunout,*)'random_notrig, tau_trig ', &
275                random_notrig, tau_trig
276           WRITE(lunout,*)'s_trig,s2,n2 ', &
277                s_trig,s2,n2
278        ENDIF
279
280        !----alp_bl computation
281        do i=1,klon
282           if ( (ale_bl_ref(i) .gt. abs(cin(i))+1.e-10) )  then
283              birth_number = n2(i)*exp(-strig(i)/s2(i))
284              birth_rate(i) = birth_number/(tau_trig(i)*cell_area(i))
285!!!              birth_rate(i) = max(birth_rate(i),1.e-18)
286              proba_notrig(i)=proba_notrig(i)*exp(-birth_number*dtime/tau_trig(i))
287              Alp_bl(i) = Alp_bl(i)* &
288                          umexp(-birth_number*cv_feed_area/cell_area(i))/ &
289                          umexp(-birth_number*dtime/tau_trig(i))*  &
290                          tau_trig(i)*cv_feed_area/(dtime*cell_area(i))
291          else
292!!jyg              proba_notrig(i)=1.
293              birth_rate(i)=0.
294              random_notrig(i)=0.
295              alp_bl(i)=0.
296           endif
297        enddo
298
299        !----ale_bl_trig computation
300         do i=1,klon
301           if (random_notrig(i) .ge. proba_notrig(i)) then
302              ale_bl_trig(i)=ale_bl_ref(i)
303           else
304              ale_bl_trig(i)=0.
305           endif
306         enddo
307
308        !
309        IF (prt_level .GE. 10) THEN
310           WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
311                proba_notrig, ale_bl_trig
312        ENDIF
313
314     endif !(iflag_trig_bl .ge. 1)
315
316!---------------------------------------
317  ENDIF ! (iflag_clos_bl .LT. 3)
318!---------------------------------------
319
320          IF (prt_level .GE. 10) THEN
321             WRITE(lunout,*)'alpale_th: ale_bl_trig, alp_bl_stat, birth_rate ', &
322                      ale_bl_trig(1), alp_bl_stat(1), birth_rate(1)
323          ENDIF
324
325          !cc fin nrlmd le 10/04/2012
326!
327          !IM/FH: 2011/02/23
328          ! Couplage Thermiques/Emanuel seulement si T<0
329          if (iflag_coupl==2) then
330             IF (prt_level .GE. 10) THEN
331                WRITE(lunout,*)'Couplage Thermiques/Emanuel seulement si T<0'
332             ENDIF
333             do i=1,klon
334                if (t_seri(i,lmax_th(i))>273.) then
335                   Ale_bl(i)=0.
336                endif
337             enddo
338!    print *,'In order to run with iflag_coupl=2, you have to comment out the following stop'
339!             STOP
340             abort_message='In order to run with iflag_coupl=2, you have to comment out the following abort'
341             CALL abort_physic(modname,abort_message,1)
342          endif
343   RETURN
344   END
345
Note: See TracBrowser for help on using the repository browser.