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

Last change on this file since 5225 was 5160, checked in by abarral, 3 months ago

Put .h into modules

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