source: LMDZ6/branches/Amaury_dev/libf/phylmd/alpale_th.F90 @ 5135

Last change on this file since 5135 was 5134, checked in by abarral, 16 months ago

Replace academic.h, alpale.h, comdissip.h, comdissipn.h, comdissnew.h by modules
Remove unused clesph0.h

  • Property svn:keywords set to Id
File size: 13.1 KB
Line 
1
2! $Id: alpale_th.F90 5134 2024-07-26 15:56:37Z abarral $
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 lmdz_ioipsl_getin_p, ONLY: getin_p
21  USE lmdz_print_control, ONLY: mydebug=>debug , lunout, prt_level
22  USE lmdz_abort_physic, ONLY: abort_physic
23  USE lmdz_alpale
24
25  IMPLICIT NONE
26
27!================================================================
28! Auteur(s)   : Jean-Yves Grandpeix, 11/05/2016
29! Objet : Contribution of the thermal scheme to Ale and Alp
30!================================================================
31
32! Input arguments
33!----------------
34  REAL, INTENT(IN)                                           :: dtime
35  REAL, DIMENSION(klon), INTENT(IN)                          :: cell_area
36  INTEGER, DIMENSION(klon), INTENT(IN)                       :: lmax_th
37  REAL, DIMENSION(klon,klev), INTENT(IN)                     :: t_seri
38  REAL, DIMENSION(klon), INTENT(IN)                          :: ale_bl_stat
39  REAL, DIMENSION(klon), INTENT(IN)                          :: cin
40  REAL, DIMENSION(klon), INTENT(IN)                          :: s2, n2, strig
41                                                               
42  REAL, DIMENSION(klon), INTENT(INOUT)                       :: ale_bl_trig, ale_bl
43  REAL, DIMENSION(klon), INTENT(INOUT)                       :: alp_bl
44  REAL, DIMENSION(klon), INTENT(INOUT)                       :: alp_bl_stat
45  REAL, DIMENSION(klon), INTENT(INOUT)                       :: proba_notrig
46
47  REAL, DIMENSION(klon), INTENT(OUT)                         :: random_notrig
48
49  REAL, DIMENSION(klon), INTENT(OUT)                         :: birth_rate
50
51! Local variables
52!----------------
53  INTEGER                                                    :: i
54  LOGICAL, SAVE                                              :: first = .TRUE.
55  LOGICAL, SAVE                                              :: multiply_proba_notrig = .FALSE.
56  REAL, SAVE                                                 :: random_notrig_max=1.
57  REAL, SAVE                                                 :: cv_feed_area
58  REAL                                                       :: birth_number
59  REAL, DIMENSION(klon)                                      :: ale_bl_ref
60  REAL, DIMENSION(klon)                                      :: tau_trig
61
62    !$OMP THREADPRIVATE(multiply_proba_notrig)
63    !$OMP THREADPRIVATE(random_notrig_max)
64    !$OMP THREADPRIVATE(cv_feed_area)
65    !$OMP THREADPRIVATE(first)
66
67 REAL umexp  ! expression of (1.-exp(-x))/x valid for all x, especially when x->0
68 REAL x
69
70     CHARACTER (LEN=20) :: modname='alpale_th'
71     CHARACTER (LEN=80) :: abort_message
72     
73 umexp(x) = max(sign(1.,x-1.e-3),0.)*(1.-exp(-x))/max(x,1.e-3) + &
74            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! correct formula            (jyg)
75!!!            (1.-max(sign(1.,x-1.e-3),0.))*(-0.5*x*(1.-x/3.*(1.-0.25*x))) !!! bug introduced by mistake  (jyg)
76!!!            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))  !!! initial correct formula (jyg)
77
78!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
79!  JYG, 20160513 : Introduction of the Effective Lifting Power (ELP), which
80! takes into account the area (cv_feed_area) covered by thermals contributing
81! to each cumulonimbus.
82!   The use of ELP prevents singularities when the trigger probability tends to
83! zero. It is activated by iflag_clos_bl = 3.
84!   The ELP values are stored in the ALP_bl variable.
85
86!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87
88    IF (first) THEN
89      CALL getin_p('multiply_proba_notrig',multiply_proba_notrig)
90      IF (iflag_clos_bl < 3) THEN
91         random_notrig_max=1.
92         CALL getin_p('random_notrig_max',random_notrig_max)
93      ELSEIF (iflag_clos_bl == 3) THEN  ! (iflag_clos_bl .LT. 3)
94         cv_feed_area = 1.e10   ! m2
95         CALL getin_p('cv_feed_area', cv_feed_area)
96      ENDIF  !! (iflag_clos_bl .LT. 3)
97      first=.FALSE.
98    ENDIF
99
100!!
101!!   Control of the multiplication of no-trigger probabilities between calls
102!! to the convection scheme. If multiply_proba_notrig is .FALSE., THEN
103!! proba_notrig is set to 1 at each CALL to alpale_th, so that only the last CALL
104!! plays a role in the triggering of convection. If it is .TRUE., then propa_notrig
105!! is saved between calls to convection and is reset to 1 only after calling the
106!! convection scheme.
107!!    For instance, if the probability of no_trigger is 0.9 at each call, and if
108!! there are 3 calls to alpale_th between calls to the convection scheme, then the
109!! probability of triggering convection will be 0.1 (= 1.-0.9) if
110!! multiply_proba_notrig is .FALSE. and 0.271 (= 1.-0.9^3) if multiply_proba_notrig
111!! is .TRUE.
112!!
113    IF (.NOT.multiply_proba_notrig) THEN
114             DO i=1,klon
115                proba_notrig(i)=1.
116             ENDDO
117    ENDIF  !! (.NOT.multiply_proba_notrig)
118!!
119!!
120!---------------------------------------
121  IF (iflag_clos_bl < 3) THEN
122!---------------------------------------
123
124!      Original code (Nicolas Rochetin)
125!     --------------------------------
126          !cc nrlmd le 10/04/2012
127          !-----------Stochastic triggering-----------
128          IF (iflag_trig_bl>=1) THEN
129             IF (prt_level >= 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) < 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 >= 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==1) THEN
156                !----Tirage al\'eatoire et calcul de ale_bl_trig
157                do i=1,klon
158                   IF ( (ale_bl_stat(i) > abs(cin(i))+1.e-10) )  THEN
159                      proba_notrig(i)=proba_notrig(i)* &
160                         (1.-exp(-strig(i)/s2(i)))**(n2(i)*dtime/tau_trig(i))
161                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
162                      IF (random_notrig(i) >= proba_notrig(i)) THEN
163                         ale_bl_trig(i)=ale_bl_stat(i)
164                      else
165                         ale_bl_trig(i)=0.
166                      endif
167                      birth_rate(i) = n2(i)*exp(-strig(i)/s2(i))/(tau_trig(i)*cell_area(i))
168!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
169                   else
170!!jyg                      proba_notrig(i)=1.
171                      birth_rate(i) = 0.
172                      random_notrig(i)=0.
173                      ale_bl_trig(i)=0.
174                   endif
175                enddo
176
177             ELSE IF (iflag_trig_bl>=2) THEN
178                do i=1,klon
179                   IF ( (Ale_bl(i) > abs(cin(i))+1.e-10) )  THEN
180                      proba_notrig(i)=proba_notrig(i)* &
181                         (1.-exp(-strig(i)/s2(i)))**(n2(i)*dtime/tau_trig(i))
182                      !        print *, 'proba_notrig(i) ',proba_notrig(i)
183                      IF (random_notrig(i) >= proba_notrig(i)) THEN
184                         ale_bl_trig(i)=Ale_bl(i)
185                      else
186                         ale_bl_trig(i)=0.
187                      endif
188                      birth_rate(i) = n2(i)*exp(-strig(i)/s2(i))/(tau_trig(i)*cell_area(i))
189!!!                      birth_rate(i) = max(birth_rate(i),1.e-18)
190                   else
191!!jyg                      proba_notrig(i)=1.
192                      birth_rate(i) = 0.
193                      random_notrig(i)=0.
194                      ale_bl_trig(i)=0.
195                   endif
196                enddo
197
198             ENDIF
199
200             IF (prt_level >= 10) THEN
201                WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
202                     proba_notrig, ale_bl_trig
203             ENDIF
204
205          endif !(iflag_trig_bl)
206
207          !-----------Statistical closure-----------
208          IF (iflag_clos_bl==1) THEN
209             do i=1,klon
210                !CR: alp probabiliste
211                IF (ale_bl_trig(i)>0.) THEN
212                   alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
213                endif
214             enddo
215
216          ELSE IF (iflag_clos_bl==2) THEN
217             !CR: alp calculee dans thermcell_main
218             do i=1,klon
219                alp_bl(i)=alp_bl_stat(i)
220             enddo
221
222          else
223
224             alp_bl_stat(:)=0.
225
226          endif !(iflag_clos_bl)
227
228!---------------------------------------
229  ELSEIF (iflag_clos_bl == 3) THEN  ! (iflag_clos_bl .LT. 3)
230!---------------------------------------
231
232!      New code with Effective Lifting Power
233!     -------------------------------------
234
235          !-----------Stochastic triggering-----------
236     IF (iflag_trig_bl>=1) THEN
237        IF (prt_level >= 10) THEN
238           WRITE(lunout,*)'cin, ale_bl_stat, alp_bl_stat ', &
239                cin, ale_bl_stat, alp_bl_stat
240        ENDIF
241
242        ! Use ale_bl_stat (Rochetin's code) or ale_bl (old code) according to
243        ! iflag_trig_bl value.
244        IF (iflag_trig_bl==1) then         ! use ale_bl_stat (Rochetin computation)
245         do i=1,klon
246              ale_bl_ref(i)=ale_bl_stat(i)
247         enddo
248        ELSE IF (iflag_trig_bl>=2) then    ! use ale_bl (old computation)
249         do i=1,klon
250              ale_bl_ref(i)=Ale_bl(i)
251         enddo
252        ENDIF ! (iflag_trig_bl.EQ.1)
253
254
255        !----Initializations and random number generation
256        do i=1,klon
257!!jyg           proba_notrig(i)=1.
258           random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
259           IF ( ale_bl_trig(i) < abs(cin(i))+1.e-10 ) THEN
260              tau_trig(i)=tau_trig_shallow
261           else
262              tau_trig(i)=tau_trig_deep
263           endif
264        enddo
265
266        IF (prt_level >= 10) THEN
267           WRITE(lunout,*)'random_notrig, tau_trig ', &
268                random_notrig, tau_trig
269           WRITE(lunout,*)'s_trig,s2,n2 ', &
270                s_trig,s2,n2
271        ENDIF
272
273        !----alp_bl computation
274        do i=1,klon
275           IF ( (ale_bl_ref(i) > abs(cin(i))+1.e-10) )  THEN
276              birth_number = n2(i)*exp(-strig(i)/s2(i))
277              birth_rate(i) = birth_number/(tau_trig(i)*cell_area(i))
278!!!              birth_rate(i) = max(birth_rate(i),1.e-18)
279              proba_notrig(i)=proba_notrig(i)*exp(-birth_number*dtime/tau_trig(i))
280              Alp_bl(i) = Alp_bl(i)* &
281                          umexp(-birth_number*cv_feed_area/cell_area(i))/ &
282                          umexp(-birth_number*dtime/tau_trig(i))*  &
283                          tau_trig(i)*cv_feed_area/(dtime*cell_area(i))
284          else
285!!jyg              proba_notrig(i)=1.
286              birth_rate(i)=0.
287              random_notrig(i)=0.
288              alp_bl(i)=0.
289           endif
290        enddo
291
292        !----ale_bl_trig computation
293         do i=1,klon
294           IF (random_notrig(i) >= proba_notrig(i)) THEN
295              ale_bl_trig(i)=ale_bl_ref(i)
296           else
297              ale_bl_trig(i)=0.
298           endif
299         enddo
300
301        IF (prt_level >= 10) THEN
302           WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
303                proba_notrig, ale_bl_trig
304        ENDIF
305
306     endif !(iflag_trig_bl .ge. 1)
307
308!---------------------------------------
309  ENDIF ! (iflag_clos_bl .LT. 3)
310!---------------------------------------
311
312          IF (prt_level >= 10) THEN
313             WRITE(lunout,*)'alpale_th: ale_bl_trig, alp_bl_stat, birth_rate ', &
314                      ale_bl_trig(1), alp_bl_stat(1), birth_rate(1)
315          ENDIF
316
317          !cc fin nrlmd le 10/04/2012
318
319          !IM/FH: 2011/02/23
320          ! Couplage Thermiques/Emanuel seulement si T<0
321          IF (iflag_coupl==2) THEN
322             IF (prt_level >= 10) THEN
323                WRITE(lunout,*)'Couplage Thermiques/Emanuel seulement si T<0'
324             ENDIF
325             do i=1,klon
326                IF (t_seri(i,lmax_th(i))>273.) THEN
327                   Ale_bl(i)=0.
328                endif
329             enddo
330!    print *,'In order to run with iflag_coupl=2, you have to comment out the following stop '
331!             STOP
332             abort_message='In order to run with iflag_coupl=2, you have to comment out the following abort'
333             CALL abort_physic(modname,abort_message,1)
334          endif
335   RETURN
336   END
337
Note: See TracBrowser for help on using the repository browser.