source: LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90 @ 5185

Last change on this file since 5185 was 5019, checked in by aborella, 6 months ago

Modification et simplification du calcul de gammasat, modification du role de temp_nowater
File size: 22.0 KB
Line 
1MODULE lmdz_lscp_tools
2
3    IMPLICIT NONE
4
5CONTAINS
6
7!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo)
9
10    ! Ref:
11    ! Stubenrauch, C. J., Bonazzola, M.,
12    ! Protopapadaki, S. E., & Musat, I. (2019).
13    ! New cloud system metrics to assess bulk
14    ! ice cloud schemes in a GCM. Journal of
15    ! Advances in Modeling Earth Systems, 11,
16    ! 3212–3234. https://doi.org/10.1029/2019MS001642
17   
18    use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc
19    use lmdz_lscp_ini, only: cice_velo, dice_velo
20
21    IMPLICIT NONE
22
23    INTEGER, INTENT(IN) :: klon
24    REAL, INTENT(IN), DIMENSION(klon) :: iwc       ! specific ice water content [kg/m3]
25    REAL, INTENT(IN), DIMENSION(klon) :: temp      ! temperature [K]
26    REAL, INTENT(IN), DIMENSION(klon) :: rho       ! dry air density [kg/m3]
27    REAL, INTENT(IN), DIMENSION(klon) :: pres      ! air pressure [Pa]
28    LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv    ! convective point  [-]
29
30    REAL, INTENT(OUT), DIMENSION(klon) :: velo    ! fallspeed velocity of crystals [m/s]
31
32
33    INTEGER i
34    REAL logvm,iwcg,tempc,phpa,fallv_tun
35    REAL m2ice, m2snow, vmice, vmsnow
36    REAL aice, bice, asnow, bsnow
37   
38
39    DO i=1,klon
40
41        IF (ptconv(i)) THEN
42            fallv_tun=ffallv_con
43        ELSE
44            fallv_tun=ffallv_lsc
45        ENDIF
46
47        tempc=temp(i)-273.15 ! celcius temp
48        iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
49        phpa=pres(i)/100.    ! pressure in hPa
50
51    IF (iflag_vice .EQ. 1) THEN
52        ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
53        if (tempc .GE. -60.0) then
54
55            logvm= -0.0000414122*tempc*tempc*log(iwcg)-0.00538922*tempc*log(iwcg) &
56                    -0.0516344*log(iwcg)+0.00216078*tempc + 1.9714   
57            velo(i)=exp(logvm)
58        else
59            velo(i)=65.0*(iwcg**0.2)*(150./phpa)**0.15
60        endif
61       
62        velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
63
64    ELSE IF (iflag_vice .EQ. 2) THEN
65        ! so called  PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
66        aice=0.587
67        bice=2.45
68        asnow=0.0444
69        bsnow=2.1
70       
71        m2ice=((iwcg*0.001/aice)/(exp(13.6-bice*7.76+0.479*bice**2)* &
72                exp((-0.0361+bice*0.0151+0.00149*bice**2)*tempc)))   &
73                **(1./(0.807+bice*0.00581+0.0457*bice**2))
74
75        vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+ &
76                 (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc) &
77                 *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice)
78
79       
80        vmice=vmice*((1000./phpa)**0.2)
81     
82        m2snow=((iwcg*0.001/asnow)/(exp(13.6-bsnow*7.76+0.479*bsnow**2)* &
83               exp((-0.0361+bsnow*0.0151+0.00149*bsnow**2)*tempc)))         &
84               **(1./(0.807+bsnow*0.00581+0.0457*bsnow**2))
85
86
87        vmsnow=100.*14.3*exp(13.6-(bsnow+.416)*7.76+0.479*(bsnow+.416)**2)&
88                  *exp((-0.0361+(bsnow+.416)*0.0151+0.00149*(bsnow+.416)**2)*tempc)&
89                  *(m2snow**(0.807+(bsnow+.416)*0.00581+0.0457*(bsnow+.416)**2))/(iwcg*0.001/asnow)
90
91        vmsnow=vmsnow*((1000./phpa)**0.35)
92        velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
93
94    ELSE
95        ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
96        velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo)
97    ENDIF
98    ENDDO
99
100END SUBROUTINE FALLICE_VELOCITY
101!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102
103!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
104SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)
105!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
106 
107  ! Compute the ice fraction 1-xliq (see e.g.
108  ! Doutriaux-Boucher & Quaas 2004, section 2.2.)
109  ! as a function of temperature
110  ! see also Fig 3 of Madeleine et al. 2020, JAMES
111!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
112
113
114    USE print_control_mod, ONLY: lunout, prt_level
115    USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
116    USE lmdz_lscp_ini, ONLY : RTT, dist_liq, temp_nowater
117
118    IMPLICIT NONE
119
120
121    INTEGER, INTENT(IN)                 :: klon              ! number of horizontal grid points
122    REAL, INTENT(IN), DIMENSION(klon)   :: temp              ! temperature
123    REAL, INTENT(IN), DIMENSION(klon)   :: distcltop         ! distance to cloud top
124    REAL, INTENT(IN), DIMENSION(klon)   :: temp_cltop        ! temperature of cloud top
125    INTEGER, INTENT(IN)                 :: iflag_ice_thermo
126    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
127    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT
128
129
130    INTEGER i
131    REAL    liqfrac_tmp, dicefrac_tmp
132    REAL    Dv, denomdep,beta,qsi,dqsidt
133    LOGICAL ice_thermo
134
135    CHARACTER (len = 20) :: modname = 'lscp_tools'
136    CHARACTER (len = 80) :: abort_message
137
138    IF ((iflag_t_glace.LT.2) .OR. (iflag_t_glace.GT.6)) THEN
139       abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
140       CALL abort_physic(modname,abort_message,1)
141    ENDIF
142
143    IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
144       abort_message = 'lscp cannot be used without ice thermodynamics'
145       CALL abort_physic(modname,abort_message,1)
146    ENDIF
147
148
149    DO i=1,klon
150 
151        ! old function with sole dependence upon temperature
152        IF (iflag_t_glace .EQ. 2) THEN
153            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
154            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
155            icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
156            IF (icefrac(i) .GT.0.) THEN
157                 dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
158                           / (t_glace_min - t_glace_max)
159            ENDIF
160
161            IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
162                 dicefracdT(i)=0.
163            ENDIF
164
165        ENDIF
166
167        ! function of temperature used in CMIP6 physics
168        IF (iflag_t_glace .EQ. 3) THEN
169            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
170            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
171            icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
172            IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN
173                dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
174                          / (t_glace_min - t_glace_max)
175            ELSE
176                dicefracdT(i)=0.
177            ENDIF
178        ENDIF
179
180        ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
181        ! and then decreases with decreasing height
182
183        !with linear function of temperature at cloud top
184        IF (iflag_t_glace .EQ. 4) THEN 
185                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
186                liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
187                icefrac(i)    =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
188                dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min)
189                dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)
190                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
191                        dicefracdT(i) = 0.
192                ENDIF
193        ENDIF
194
195        ! with CMIP6 function of temperature at cloud top
196        IF (iflag_t_glace .EQ. 5) THEN
197                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
198                liqfrac_tmp =  MIN(MAX(liqfrac_tmp,0.0),1.0)
199                liqfrac_tmp = liqfrac_tmp**exposant_glace
200                icefrac(i)  =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
201                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
202                        dicefracdT(i) = 0.
203                ELSE
204                        dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) &
205                                        *exp(-distcltop(i)/dist_liq)
206                ENDIF
207        ENDIF
208
209        ! with modified function of temperature at cloud top
210        ! to get largere values around 260 K, works well with t_glace_min = 241K
211        IF (iflag_t_glace .EQ. 6) THEN
212                IF (temp(i) .GT. t_glace_max) THEN
213                        liqfrac_tmp = 1.
214                ELSE
215                        liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1.
216                ENDIF
217                liqfrac_tmp  = MIN(MAX(liqfrac_tmp,0.0),1.0)
218                icefrac(i)   = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)       
219                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
220                        dicefracdT(i) = 0.
221                ELSE
222                        dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) &
223                                  *exp(-distcltop(i)/dist_liq)
224                ENDIF
225        ENDIF
226
227        ! if temperature of cloud top <-40°C,
228        IF (iflag_t_glace .GE. 4) THEN
229                IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN
230                        icefrac(i) = 1.
231                        dicefracdT(i) = 0.
232                ENDIF
233        ENDIF
234
235
236     ENDDO ! klon
237 
238     RETURN
239 
240END SUBROUTINE ICEFRAC_LSCP
241!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
242
243
244
245SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
246!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
247    ! Calculate qsat following ECMWF method
248!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
249
250
251    IMPLICIT NONE
252
253    include "YOMCST.h"
254    include "YOETHF.h"
255    include "FCTTRE.h"
256
257    INTEGER, INTENT(IN) :: klon  ! number of horizontal grid points
258    REAL, INTENT(IN), DIMENSION(klon) :: temp     ! temperature in K
259    REAL, INTENT(IN), DIMENSION(klon) :: qtot     ! total specific water in kg/kg
260    REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
261    REAL, INTENT(IN)                  :: tref     ! reference temperature in K
262    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
263    INTEGER, INTENT(IN) :: phase
264    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
265    !        1=liquid
266    !        2=solid
267
268    REAL, INTENT(OUT), DIMENSION(klon) :: qs      ! saturation specific humidity [kg/kg]
269    REAL, INTENT(OUT), DIMENSION(klon) :: dqs     ! derivation of saturation specific humidity wrt T
270
271    REAL delta, cor, cvm5
272    INTEGER i
273
274    DO i=1,klon
275
276    IF (phase .EQ. 1) THEN
277        delta=0.
278    ELSEIF (phase .EQ. 2) THEN
279        delta=1.
280    ELSE
281        delta=MAX(0.,SIGN(1.,tref-temp(i)))
282    ENDIF
283
284    IF (flagth) THEN
285    cvm5=R5LES*(1.-delta) + R5IES*delta
286    ELSE
287    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
288    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
289    ENDIF
290
291    qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
292    qs(i)=MIN(0.5,qs(i))
293    cor=1./(1.-RETV*qs(i))
294    qs(i)=qs(i)*cor
295    dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
296
297    END DO
298
299END SUBROUTINE CALC_QSAT_ECMWF
300!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
301
302
303!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
304SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
305
306!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
307    ! programme that calculates the gammasat parameter that determines the
308    ! homogeneous condensation thresholds for cold (<0oC) clouds
309    ! condensation at q>gammasat*qsat
310    ! Etienne Vignon, March 2021
311!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
312
313    use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT
314    use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez
315
316    IMPLICIT NONE
317
318
319    INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
320    REAL, INTENT(IN), DIMENSION(klon) :: temp         ! temperature in K
321    REAL, INTENT(IN), DIMENSION(klon) :: qtot         ! total specific water in kg/kg
322
323    REAL, INTENT(IN), DIMENSION(klon) :: pressure     ! pressure in Pa
324
325    REAL, INTENT(OUT), DIMENSION(klon) :: gammasat    ! coefficient to multiply qsat with to calculate saturation
326    REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
327
328    REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
329    REAL  f_homofreez, fac
330
331    INTEGER i
332   
333        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
334        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
335
336    DO i = 1, klon
337
338        IF ( temp(i) .GE. RTT ) THEN
339            ! warm clouds: condensation at saturation wrt liquid
340            gammasat(i) = 1.
341            dgammasatdt(i) = 0.
342
343        ELSE
344            ! cold clouds: qsi > qsl
345           
346            ! homogeneous freezing of aerosols, according to
347            ! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
348            ! 'Cirrus regime'
349            ! if f_homofreez > qsl / qsi, liquid nucleation
350            ! if f_homofreez < qsl / qsi, homogeneous freezing of aerosols
351            ! Note: f_homofreez = qsl / qsi for temp ~= -38degC
352            f_homofreez = a_homofreez - temp(i) / b_homofreez
353           
354            IF ( iflag_gammasat .GE. 3 ) THEN
355              ! condensation at homogeneous freezing threshold for temp < -38 degC
356              ! condensation at liquid saturation for temp > -38 degC
357              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
358                gammasat(i) = f_homofreez
359                dgammasatdt(i) = - 1. / b_homofreez
360              ELSE
361                gammasat(i) = qsl(i) / qsi(i)
362                dgammasatdt(i) = ( dqsl(i) * qsi(i) - dqsi(i) * qsl(i) ) / qsi(i) / qsi(i)
363              ENDIF
364
365            ELSEIF ( iflag_gammasat .EQ. 2 ) THEN
366              ! condensation at homogeneous freezing threshold for temp < -38 degC
367              ! condensation at a threshold linearly decreasing between homogeneous
368              ! freezing and ice saturation for -38 degC < temp < temp_nowater
369              ! condensation at ice saturation for temp > temp_nowater
370              ! If temp_nowater = 235.15 K, this is equivalent to iflag_gammasat = 1
371              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
372                gammasat(i) = f_homofreez
373                dgammasatdt(i) = - 1. / b_homofreez
374              ELSEIF ( temp(i) .LE. temp_nowater ) THEN
375                ! Here, we assume that f_homofreez = qsl / qsi for temp = -38 degC = 235.15 K
376                dgammasatdt(i) = ( a_homofreez - 235.15 / b_homofreez - 1. ) &
377                               / ( 235.15 - temp_nowater )
378                gammasat(i) = dgammasatdt(i) * ( temp(i) - temp_nowater ) + 1.
379              ELSE
380                gammasat(i) = 1.
381                dgammasatdt(i) = 0.
382              ENDIF
383
384            ELSEIF ( iflag_gammasat .EQ. 1 ) THEN
385              ! condensation at homogeneous freezing threshold for temp < -38 degC
386              ! condensation at ice saturation for temp > -38 degC
387              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
388                gammasat(i) = f_homofreez
389                dgammasatdt(i) = - 1. / b_homofreez
390              ELSE
391                gammasat(i) = 1.
392                dgammasatdt(i) = 0.
393              ENDIF
394
395            ELSE
396              ! condensation at ice saturation for temp < -38 degC
397              ! condensation at ice saturation for temp > -38 degC
398              gammasat(i) = 1.
399              dgammasatdt(i) = 0.
400
401            ENDIF
402
403            ! Note that the delta_hetfreez parameter allows to linearly decrease the
404            ! condensation threshold between the calculated threshold and the ice saturation
405            ! for delta_hetfreez = 1, the threshold is the calculated condensation threshold
406            ! for delta_hetfreez = 0, the threshold is the ice saturation
407            gammasat(i) = ( 1. - delta_hetfreez ) + delta_hetfreez * gammasat(i)
408            dgammasatdt(i) = delta_hetfreez * dgammasatdt(i)
409
410        ENDIF
411   
412    END DO
413
414
415END SUBROUTINE CALC_GAMMASAT
416!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
417
418
419!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
420SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
421!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
422   
423   USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl
424
425   IMPLICIT NONE
426   
427   INTEGER, INTENT(IN) :: klon,klev                !number of horizontal and vertical grid points
428   INTEGER, INTENT(IN) :: k                        ! vertical index
429   REAL, INTENT(IN), DIMENSION(klon,klev) :: temp  ! temperature in K
430   REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
431   REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
432   REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb  ! cloud fraction
433
434   REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D  ! distance from cloud top
435   REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop     ! temperature of cloud top
436   
437   REAL dzlay(klon,klev)
438   REAL zlay(klon,klev)
439   REAL dzinterf
440   INTEGER i,k_top, kvert
441   LOGICAL bool_cl
442
443
444   DO i=1,klon
445         ! Initialization height middle of first layer
446          dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
447          zlay(i,1) = dzlay(i,1)/2
448
449          DO kvert=2,klev
450                 IF (kvert.EQ.klev) THEN
451                       dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
452                 ELSE
453                       dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
454                 ENDIF
455                       dzinterf       = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
456                       zlay(i,kvert)  = zlay(i,kvert-1) + dzinterf
457           ENDDO
458   ENDDO
459   
460   DO i=1,klon
461          k_top = k
462          IF (rneb(i,k) .LE. tresh_cl) THEN
463                 bool_cl = .FALSE.
464          ELSE
465                 bool_cl = .TRUE.
466          ENDIF
467
468          DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
469          ! find cloud top
470                IF (rneb(i,k_top) .GT. tresh_cl) THEN
471                      k_top = k_top + 1
472                ELSE
473                      bool_cl = .FALSE.
474                      k_top   = k_top - 1
475                ENDIF
476          ENDDO
477          k_top=min(k_top,klev)
478
479          !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
480          !interf for layer of cloud top
481          distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
482          temp_cltop(i)  = temp(i,k_top)
483   ENDDO ! klon
484
485END SUBROUTINE DISTANCE_TO_CLOUD_TOP
486!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
487
488!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
489FUNCTION GAMMAINC ( p, x )
490
491!*****************************************************************************80
492!
493!! GAMMAINC computes the regularized lower incomplete Gamma Integral
494!
495!  Modified:
496!
497!    20 January 2008
498!
499!  Author:
500!
501!    Original FORTRAN77 version by B Shea.
502!    FORTRAN90 version by John Burkardt.
503!
504!  Reference:
505!
506!    B Shea,
507!    Algorithm AS 239:
508!    Chi-squared and Incomplete Gamma Integral,
509!    Applied Statistics,
510!    Volume 37, Number 3, 1988, pages 466-473.
511!
512!  Parameters:
513!
514!    Input, real X, P, the parameters of the incomplete
515!    gamma ratio.  0 <= X, and 0 < P.
516!
517!    Output, real GAMMAINC, the value of the incomplete
518!    Gamma integral.
519!
520  IMPLICIT NONE
521
522  REAL A
523  REAL AN
524  REAL ARG
525  REAL B
526  REAL C
527  REAL, PARAMETER :: ELIMIT = - 88.0E+00
528  REAL GAMMAINC
529  REAL, PARAMETER :: OFLO = 1.0E+37
530  REAL P
531  REAL, PARAMETER :: PLIMIT = 1000.0E+00
532  REAL PN1
533  REAL PN2
534  REAL PN3
535  REAL PN4
536  REAL PN5
537  REAL PN6
538  REAL RN
539  REAL, PARAMETER :: TOL = 1.0E-14
540  REAL X
541  REAL, PARAMETER :: XBIG = 1.0E+08
542
543  GAMMAINC = 0.0E+00
544
545  IF ( X == 0.0E+00 ) THEN
546    GAMMAINC = 0.0E+00
547    RETURN
548  END IF
549!
550!  IF P IS LARGE, USE A NORMAL APPROXIMATION.
551!
552  IF ( PLIMIT < P ) THEN
553
554    PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) &
555    + 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 )
556
557    GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) )
558    RETURN
559
560  END IF
561!
562!  IF X IS LARGE SET GAMMAD = 1.
563!
564  IF ( XBIG < X ) THEN
565    GAMMAINC = 1.0E+00
566    RETURN
567  END IF
568!
569!  USE PEARSON'S SERIES EXPANSION.
570!  (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM).
571!
572  IF ( X <= 1.0E+00 .OR. X < P ) THEN
573
574    ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 )
575    C = 1.0E+00
576    GAMMAINC = 1.0E+00
577    A = P
578
579    DO
580
581      A = A + 1.0E+00
582      C = C * X / A
583      GAMMAINC = GAMMAINC + C
584
585      IF ( C <= TOL ) THEN
586        EXIT
587      END IF
588
589    END DO
590
591    ARG = ARG + LOG ( GAMMAINC )
592
593    IF ( ELIMIT <= ARG ) THEN
594      GAMMAINC = EXP ( ARG )
595    ELSE
596      GAMMAINC = 0.0E+00
597    END IF
598!
599!  USE A CONTINUED FRACTION EXPANSION.
600!
601  ELSE
602
603    ARG = P * LOG ( X ) - X - LOG_GAMMA ( P )
604    A = 1.0E+00 - P
605    B = A + X + 1.0E+00
606    C = 0.0E+00
607    PN1 = 1.0E+00
608    PN2 = X
609    PN3 = X + 1.0E+00
610    PN4 = X * B
611    GAMMAINC = PN3 / PN4
612
613    DO
614
615      A = A + 1.0E+00
616      B = B + 2.0E+00
617      C = C + 1.0E+00
618      AN = A * C
619      PN5 = B * PN3 - AN * PN1
620      PN6 = B * PN4 - AN * PN2
621
622      IF ( PN6 /= 0.0E+00 ) THEN
623
624        RN = PN5 / PN6
625
626        IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN
627          EXIT
628        END IF
629
630        GAMMAINC = RN
631
632      END IF
633
634      PN1 = PN3
635      PN2 = PN4
636      PN3 = PN5
637      PN4 = PN6
638!
639!  RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE.
640!
641      IF ( OFLO <= ABS ( PN5 ) ) THEN
642        PN1 = PN1 / OFLO
643        PN2 = PN2 / OFLO
644        PN3 = PN3 / OFLO
645        PN4 = PN4 / OFLO
646      END IF
647
648    END DO
649
650    ARG = ARG + LOG ( GAMMAINC )
651
652    IF ( ELIMIT <= ARG ) THEN
653      GAMMAINC = 1.0E+00 - EXP ( ARG )
654    ELSE
655      GAMMAINC = 1.0E+00
656    END IF
657
658  END IF
659
660  RETURN
661END FUNCTION GAMMAINC
662!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
663
664END MODULE lmdz_lscp_tools
665
666
Note: See TracBrowser for help on using the repository browser.