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

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

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