source: LMDZ6/trunk/libf/phylmd/lmdz_ratqs_multi.f90 @ 5308

Last change on this file since 5308 was 5284, checked in by abarral, 4 days ago

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