source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/lmdz_ratqs_multi.F90 @ 4987

Last change on this file since 4987 was 4727, checked in by idelkadi, 13 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch
File size: 16.9 KB
Line 
1MODULE lmdz_ratqs_multi
2
3!=============================================
4! A FAIRE :
5! Traiter le probleme de USE lmdz_lscp_tools, ONLY: CALC_QSAT_ECMWF
6!=============================================
7
8!=============================================
9! module containing subroutines that take
10! into account the effect of convection, orography,
11! surface heterogeneities and subgrid-scale
12! turbulence on ratqs, i.e. on the width of the
13! total water subgrid distribution.
14!=============================================
15
16IMPLICIT NONE
17
18! Include
19!=============================================
20    INCLUDE "YOETHF.h"
21
22
23      CONTAINS
24
25
26!========================================================================   
27SUBROUTINE ratqs_inter(klon,klev,iflag_ratqs,pdtphys,paprs, &
28           ratqsbas, wake_deltaq, wake_s, q_seri,qtc_cv, sigt_cv,     &
29           fm_therm,entr_therm,detr_therm,detrain_cv,fm_cv,fqd,fqcomp,sigd, &
30           ratqs_inter_)
31
32USE lmdz_ratqs_ini, ONLY : a_ratqs_cv,tau_var,fac_tau,tau_cumul,a_ratqs_wake, dqimpl
33USE lmdz_ratqs_ini, ONLY : RG
34USE lmdz_ratqs_ini, ONLY : povariance, var_conv
35USE lmdz_thermcell_dq,  ONLY : thermcell_dq
36
37implicit none
38
39!========================================================================
40! L. d'Alençon, 25/02/2021
41! Cette subroutine calcule une valeur de ratqsbas interactive 
42! Elle est appelée par la subroutine ratqs lorsque iflag_ratqs = 11.
43!========================================================================
44
45! Declarations
46! Input
47integer,intent(in) :: klon,klev,iflag_ratqs
48real,intent(in) :: pdtphys,ratqsbas
49real, dimension(klon,klev+1),intent(in) :: paprs
50real, dimension(klon,klev),intent(in) :: wake_deltaq, q_seri,qtc_cv, sigt_cv
51real, dimension(klon),intent(in) :: wake_s
52real, dimension(klon,klev+1),intent(in) :: fm_therm
53real, dimension(klon,klev),intent(in) :: entr_therm,detr_therm,detrain_cv,fm_cv,fqd,fqcomp
54real, dimension(klon),intent(in) :: sigd
55
56! Output
57real, dimension(klon,klev),intent(inout) :: ratqs_inter_
58
59! local
60LOGICAL :: klein = .false.
61LOGICAL :: klein_conv = .true.
62REAL :: taup0 = 70000
63REAL :: taudp = 500
64integer :: lev_out=10
65REAL, DIMENSION (klon,klev) :: zmasse,entr0,detr0,detraincv,dqp,detrain_p,q0,qd0,tau_diss
66REAL, DIMENSION (klon,klev+1) :: fm0
67integer i,k
68real, dimension(klon,klev) :: wake_dq
69
70
71real, dimension(klon) :: max_sigd, max_dqconv,max_sigt
72real, dimension(klon,klev) :: zoa,zocarrea,pdocarreadj,pocarre,po,pdoadj,varq_therm
73
74
75lev_out=0.
76
77print*,'ratqs_inter'
78
79!-----------------------------------------------------------------------
80!   Calcul des masses
81!-----------------------------------------------------------------------
82
83      do k=1,klev
84         zmasse(:,k)=(paprs(:,k)-paprs(:,k+1))/RG
85      enddo
86!-------------------------------------------------------------------------
87!  Caclul du terme de détrainement de la variance pour les thermiques
88!-------------------------------------------------------------------------
89
90! initialisations
91
92
93      do k=1,klev
94         do i=1,klon
95            tau_diss(i,k)=tau_var +0.5*fac_tau*tau_var*(tanh((taup0-paprs(i,k))/taudp) + 1.)
96         enddo
97      enddo
98       
99                   
100     
101      entr0(:,:) = entr_therm(:,:)
102      fm0(:,:) = fm_therm(:,:) 
103      detr0(:,:) = detr_therm(:,:)
104
105! calcul du carré de l'humidité spécifique
106      po(:,:) = q_seri(:,:)
107      call thermcell_dq(klon,klev,dqimpl,pdtphys,fm0,entr0,zmasse,  &
108     &                   po,pdoadj,zoa,lev_out)
109      do k=1,klev
110         do i=1,klon
111            pocarre(i,k)=po(i,k)*po(i,k) + povariance(i,k)
112         enddo
113      enddo
114      call thermcell_dq(klon,klev,dqimpl,pdtphys,fm0,entr0,zmasse,  &
115      &                   pocarre,pdocarreadj,zocarrea,lev_out)
116
117     
118      do k=1,klev
119         do i=1,klon     
120            varq_therm(i,k)=zocarrea(i,k)-zoa(i,k)*zoa(i,k)
121         enddo
122      enddo
123
124     
125      if (klein) then
126         do k=1,klev-1
127            do i=1,klon
128              qd0(:,:) = 0.0
129              if (sigd(i).ne.0) then
130                qd0(i,k) = fqd(i,k)*tau_cumul/sigd(i)
131              endif
132            enddo
133         enddo
134         do k=1,klev-1
135            do i=1,klon     
136               povariance(i,k)= (detr0(i,k)*((zoa(i,k)-po(i,k))**2 + (varq_therm(i,k)-povariance(i,k)))/zmasse(i,k) + fm0(i,k+1)*povariance(i,k+1)/zmasse(i,k) -  &
137               fm0(i,k)*povariance(i,k)/zmasse(i,k) +  &
138               a_ratqs_cv*(detrain_cv(i,k)/zmasse(i,k)) +  sigd(i)*(1-sigd(i))*qd0(i,k)**2/tau_cumul  &   
139                + ((povariance(i,k+1)-povariance(i,k))*(fm_cv(i,k)/zmasse(i,k))))*pdtphys + povariance(i,k)
140               povariance(i,k)= povariance(i,k)*exp(-pdtphys/tau_diss(i,k))
141            enddo
142         enddo
143         povariance(:,klev) = povariance(:,klev-1)
144         
145      else ! calcul direct     
146         qd0(:,:) = 0.0
147         q0(:,:) = 0.0
148         do k=1,klev-1
149            do i=1,klon
150               if (sigd(i).ne.0) then
151                 qd0(i,k) = fqd(i,k)*tau_cumul/sigd(i)
152               endif
153               if (sigt_cv(i,k).ne.0) then
154                 q0(i,k) = fqcomp(i,k)*tau_cumul/sigt_cv(i,k)
155               endif
156            enddo
157         enddo
158         do k=1,klev-1
159            do i=1,klon     
160               povariance(i,k)= (pdocarreadj(i,k)-2.*po(i,k)*pdoadj(i,k) +  &
161               a_ratqs_cv*(sigt_cv(i,k)*(1-sigt_cv(i,k))*q0(i,k)**2/tau_cumul + sigd(i)*(1-sigd(i))*qd0(i,k)**2/tau_cumul) +  &
162               ((povariance(i,k+1)-povariance(i,k))*(fm_cv(i,k)/zmasse(i,k))))*pdtphys + povariance(i,k)
163               povariance(i,k)=povariance(i,k)*exp(-pdtphys/tau_diss(i,k))
164            enddo
165         enddo
166         povariance(:,klev) = povariance(:,klev-1)
167!         fqd(:,:)=sigt_cv(:,:)*(1-sigt_cv(:,:))*q0(:,:)**2/tau_cumul
168      endif
169     
170!-------------------------------------------------------------------------
171!  Caclul du terme de détrainement de la variance pour la convection (version fausse avec deux calculs de variance indépendants)
172!-------------------------------------------------------------------------
173
174!      if (klein_conv) then
175!        detrain_p(:,:) = 0.
176!        detraincv(:,:) = 0.
177!        dqp(:,:)    = 0.
178     
179!        do k=1,klev-1
180!          do i=1,klon
181!           dqp(i,k) = q_seri(i,k) - qp(i,k)
182!           detraincv(i,k) = abs(detrain_cv(i,k))
183         
184!         if ((mp(i,k)-mp(i,k+1)).le.0) then
185!            detrain_p(i,k) = (mp(i,k+1)-mp(i,k))
186!         endif
187!          enddo
188!        enddo
189
190!        do k=1,klev-1
191!          do i=1,klon
192!            qd0(:,:) = 0.0
193!            if (sigd(i).ne.0) then
194!              qd0(i,k) = fqd(i,k)*tau_cumul/sigd(i)
195!            endif
196!          enddo
197!        enddo
198
199!        do k=1,klev-1
200!          do i=1,klon           
201!            var_conv(i,k)= var_conv(i,k)*exp(-pdtphys/tau_conv) +  &
202!                (a_ratqs_cv*(detraincv(i,k)/zmasse(i,k)) +  sigd(i)*(1-sigd(i))*qd0(i,k)**2/tau_cumul  &   
203!                + ((var_conv(i,k+1)-var_conv(i,k))*(fm_cv(i,k)/zmasse(i,k))))* &
204!              ((dqp(i,k)**2)*(detrain_p(i,k)/zmasse(i,k))))* & ! les termes de descentes précipitantes qui seront traités autrement
205!              (((mp(i,k)-mp(i,k+1))/zmasse(i,k))*dqp(i,k)**2) + (2*mp(i,k)*dqp(i,k)*(qp(i,k+1)-qp(i,k))/zmasse(i,k)))* &
206!                (1.-exp(-pdtphys/tau_conv))/(1/tau_conv)     
207!          enddo
208!        enddo
209!        var_conv(:,klev) = var_conv(:,klev-1)
210!        fqd(:,:) = var_conv(:,:)
211!      else
212
213
214!        do k=1,klev-1
215!          do i=1,klon           
216!            var_conv(i,k)= var_conv(i,k)*exp(-pdtphys/tau_conv) +  &
217!                (a_ratqs_cv*(sigt_cv(i,k)*(1-sigt_cv(i,k))*q0(i,k)**2/pdtphys + sigd(i)*(1-sigd(i))*qd0(i,k)**2/tau_cumul) +  &   ! somme des contributions des descentes précipitantes et des flux mélangés
218!                + ((var_conv(i,k+1)-var_conv(i,k))*(fm_cv(i,k)/zmasse(i,k))))* & ! flux compensatoires dans l'environnement
219!                (1.-exp(-pdtphys/tau_conv))/(1/tau_conv)
220!             ((var_conv(i,k+1)-var_conv(i,k))*(fm_cv(i,k)/zmasse(i,k))) 
221!          enddo
222!        enddo
223!        var_conv(:,klev) = var_conv(:,klev-1)
224!      endif
225
226!-------------------------------------------------------------------------
227!  Caclul du ratqs_inter_
228!-------------------------------------------------------------------------
229
230      do k=1,klev
231        do i=1,klon
232           if(q_seri(i,k).ge.1E-7) then
233               ratqs_inter_(i,k) = abs(povariance(i,k))**0.5/q_seri(i,k)
234           else
235               ratqs_inter_(i,k) = 0. 
236           endif
237        enddo
238      enddo
239     
240return
241end             
242
243!------------------------------------------------------------------
244SUBROUTINE ratqs_oro(klon,klev,pctsrf,zstd,qsat,temp,pplay,paprs,ratqs_oro_)
245
246! Etienne Vignon, November 2021: effect of subgrid orography on ratqs
247
248USE lmdz_ratqs_ini, ONLY : RG,RV,RD,RLSTT,RLVTT,RTT,nbsrf,is_lic,is_ter
249
250IMPLICIT NONE
251
252! Declarations
253!--------------
254
255! INPUTS
256
257INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
258INTEGER, INTENT(IN) :: klev                       ! number of vertical layers
259REAL, DIMENSION(klon,nbsrf) :: pctsrf
260REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat    ! saturation specific humidity [kg/kg]
261REAL, DIMENSION(klon), INTENT(IN) :: zstd    ! sub grid orography standard deviation
262REAL, DIMENSION(klon,klev), INTENT(IN) :: temp    ! air temperature [K]
263REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay    ! air pressure, layer's center [Pa]
264REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs    ! air pressure, lower inteface [Pa]
265
266! OUTPUTS
267
268REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_oro_ ! ratqs profile due to subgrid orography
269
270
271! LOCAL
272
273INTEGER :: i,k
274REAL, DIMENSION(klon) :: orogradT,xsi0
275REAL, DIMENSION (klon,klev) :: zlay
276REAL :: Lvs, temp0
277
278
279! Calculation of the near-surface temperature gradient along the topography
280!--------------------------------------------------------------------------
281
282! at the moment, we fix it at a constant value (moist adiab. lapse rate)
283
284orogradT(:)=-6.5/1000. ! K/m
285
286! Calculation of near-surface surface ratqs
287!-------------------------------------------
288
289DO i=1,klon
290    temp0=temp(i,1)
291    IF (temp0 .LT. RTT) THEN
292        Lvs=RLSTT
293    ELSE
294        Lvs=RLVTT
295    ENDIF
296    xsi0(i)=zstd(i)*ABS(orogradT(i))*Lvs/temp0/temp0/RV
297    ratqs_oro_(i,1)=xsi0(i)
298END DO
299
300! Vertical profile of ratqs assuming an exponential decrease with height
301!------------------------------------------------------------------------
302     
303! calculation of geop. height AGL       
304zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) &
305           *(paprs(:,1)-pplay(:,1))/RG
306
307DO k=2,klev
308   DO i = 1, klon
309      zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) &
310               /paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG
311               
312      ratqs_oro_(i,k)=MAX(0.0,pctsrf(i,is_ter)*xsi0(i)*exp(-zlay(i,k)/MAX(zstd(i),1.)))   
313    END DO
314END DO
315
316
317
318
319END SUBROUTINE ratqs_oro
320
321!=============================================
322
323SUBROUTINE ratqs_hetero(klon,klev,pctsrf,s_pblh,t2m,q2m,temp,q,pplay,paprs,ratqs_hetero_)
324
325! Etienne Vignon, November 2021
326! Effect of subgrid surface heterogeneities on ratqs
327
328USE lmdz_lscp_tools, ONLY: CALC_QSAT_ECMWF
329
330USE lmdz_ratqs_ini, ONLY : RG,RD,RTT,nbsrf
331
332IMPLICIT NONE
333
334! INPUTS
335
336
337INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
338INTEGER, INTENT(IN) :: klev                       ! number of vertical layers
339REAL, DIMENSION(klon)                   :: s_pblh ! height of the planetary boundary layer(HPBL)
340REAL, DIMENSION(klon,nbsrf)             :: pctsrf ! Fractional cover of subsurfaces
341REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: t2m    ! 2m temperature for each tile [K]
342REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: q2m    ! 2m specific humidity for each tile [kg/kg]
343REAL, DIMENSION(klon,klev), INTENT(IN) :: temp    ! air temperature [K]
344REAL, DIMENSION(klon,klev), INTENT(IN) :: q       ! specific humidity [kg/kg]
345REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay   ! air pressure, layer's center [Pa]
346REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! air pressure, lower inteface [Pa]
347
348! OUTPUTS
349
350REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_hetero_ ! ratsq profile due to surface heterogeneities
351
352
353INTEGER :: i,k,nsrf
354REAL, DIMENSION(klon) :: xsi0, ratiom, qsat2m, dqsatdT
355REAL, DIMENSION (klon,klev) :: zlay
356
357
358
359! Calculation of near-surface surface ratqs
360!-------------------------------------------
361
362   
363    ratiom(:)=0.
364    xsi0(:)=0.
365   
366    DO nsrf=1,nbsrf
367    CALL CALC_QSAT_ECMWF(klon,t2m(:,nsrf),q2m(:,nsrf),paprs(:,1),RTT,0,.false.,qsat2m,dqsatdT)
368    ratiom(:)=ratiom(:)+pctsrf(:,nsrf)*(q2m(:,nsrf)/qsat2m(:))
369    xsi0(:)=xsi0(:)+pctsrf(:,nsrf)*((q2m(:,nsrf)/qsat2m(:)-ratiom(:))**2)
370    END DO
371   
372    xsi0(:)=sqrt(xsi0(:))/(ratiom(:)+1E-6)
373
374
375
376! Vertical profile of ratqs assuming an exponential decrease with height
377!------------------------------------------------------------------------
378       
379! calculation of geop. height AGL
380
381zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) &
382           *(paprs(:,1)-pplay(:,1))/RG
383ratqs_hetero_(:,1)=xsi0(:)
384
385DO k=2,klev
386   DO i = 1, klon
387      zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) &
388               /paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG
389               
390      ratqs_hetero_(i,k)=MAX(xsi0(i)*exp(-zlay(i,k)/(s_pblh(i)+1.0)),0.0)   
391    END DO
392END DO
393
394END SUBROUTINE ratqs_hetero
395
396!=============================================
397
398SUBROUTINE ratqs_tke(klon,klev,pdtphys,temp,q,qsat,pplay,paprs,omega,tke,tke_dissip,lmix,wprime,ratqs_tke_)
399
400! References:
401!
402! Etienne Vignon: effect of subgrid turbulence on ratqs
403!
404! Field, P.R., Hill, A., Furtado, K., Korolev, A., 2014b. Mixed-phase clouds in a turbulent environment. Part
405! 2: analytic treatment. Q. J. R. Meteorol. Soc. 21, 2651–2663. https://doi.org/10.1002/qj.2175.
406!
407! Furtado, K., Field, P.R., Boutle, I.A., Morcrette, C.R., Wilkinson, J., 2016. A physically-based, subgrid
408! parametrization for the production and maintenance of mixed-phase clouds in a general circulation
409! model. J. Atmos. Sci. 73, 279–291. https://doi.org/10.1175/JAS-D-15-0021.
410
411USE lmdz_ratqs_ini, ONLY : RG,RV,RD,RCPD,RLSTT,RLVTT,RTT
412
413IMPLICIT NONE
414
415! INPUTS
416
417INTEGER, INTENT(IN) :: klon                             ! number of horizontal grid points
418INTEGER, INTENT(IN) :: klev                             ! number of vertical layers
419REAL, INTENT(IN) :: pdtphys                             ! physics time step [s]
420REAL, DIMENSION(klon,klev), INTENT(IN) :: temp          ! air temperature [K]
421REAL, DIMENSION(klon,klev), INTENT(IN) :: q             ! specific humidity [kg/kg]
422REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat          ! saturation specific humidity [kg/kg]
423REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay         ! air pressure, layer's center [Pa]
424REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs       ! air pressure, lower inteface [Pa]
425REAL, DIMENSION(klon,klev), INTENT(IN) :: omega       ! air pressure, lower inteface [Pa]
426REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke         ! Turbulent Kinetic Energy [m2/s2]
427REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke_dissip  ! Turbulent Kinetic Energy Dissipation rate [m2/s3]
428REAL, DIMENSION(klon,klev+1), INTENT(IN) :: lmix  ! Turbulent mixing length
429REAL, DIMENSION(klon,klev+1), INTENT(IN) :: wprime      ! Turbulent vertical velocity scale [m/s]
430
431! OUTPUTS
432
433REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_tke_  ! ratsq profile due to subgrid TKE
434
435! LOCAL
436INTEGER :: i, k
437REAL :: AA, DD, NW, AAprime, VARLOG,rho,Lvs,taue,lhomo,dissmin,maxvarlog
438REAL, DIMENSION(klon,klev) :: sigmaw,w
439REAL, PARAMETER :: C0=10.0
440REAL, PARAMETER :: lmin=0.001
441REAL, PARAMETER :: ratqsmin=1E-6
442REAL, PARAMETER :: ratqsmax=0.5
443
444
445! Calculation of large scale and turbulent vertical velocities
446!---------------------------------------------------------------
447
448DO k=1,klev
449    DO i=1,klon
450        rho=pplay(i,k)/temp(i,k)/RD
451        w(i,k)=-rho*RG*omega(i,k)
452        sigmaw(i,k)=0.5*(wprime(i,k+1)+wprime(i,k)) ! turbulent vertical velocity at the middle of model layers.
453    END DO
454END DO
455
456! Calculation of ratqs
457!---------------------------------------------------------------
458ratqs_tke_(:,1)=ratqsmin ! set to a very low value to avoid division by 0 in order parts
459                        ! of the code
460DO k=2,klev ! we start from second model level since TKE is not defined at k=1
461    DO i=1,klon
462
463       IF (temp(i,k) .LT. RTT) THEN
464           Lvs=RLSTT
465       ELSE
466           Lvs=RLVTT
467       ENDIF
468       dissmin=0.01*(0.5*(tke(i,k)+tke(i,k+1))/pdtphys)
469       maxvarlog=LOG(1.0+ratqsmax**2)! to prevent ratqs from exceeding an arbitrary threshold value
470       AA=RG*(Lvs/(RCPD*temp(i,k)*temp(i,k)*RV) - 1./(RD*temp(i,k)))
471       lhomo=MAX(0.5*(lmix(i,k)+lmix(i,k+1)),lmin)
472       taue=(lhomo*lhomo/MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin))**(1./3) ! Fields et al. 2014
473       DD=1.0/taue
474       NW=(sigmaw(i,k)**2)*SQRT(2./(C0*MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin)))
475       AAprime=AA*NW
476       VARLOG=AAprime/2./DD
477       VARLOG=MIN(VARLOG,maxvarlog)
478       ratqs_tke_(i,k)=SQRT(MAX(EXP(VARLOG)-1.0,ratqsmin))
479       END DO
480END DO
481END SUBROUTINE ratqs_tke
482
483END MODULE lmdz_ratqs_multi
Note: See TracBrowser for help on using the repository browser.