source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_ratqs_multi.F90 @ 5116

Last change on this file since 5116 was 5116, checked in by abarral, 4 months ago

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