source: LMDZ6/trunk/libf/phylmd/calcratqs_multi_mod.F90 @ 4041

Last change on this file since 4041 was 4010, checked in by Ehouarn Millour, 3 years ago

Fix OpenMP bug in previous revision.
While at it removed dead link in phylmdiso.
EM
File size: 12.1 KB
Line 
1MODULE calcratqs_multi_mod
2
3!=============================================
4! module containing subroutines that take
5! into account the effect of convection, orography,
6! surface heterogeneities and subgrid-scale
7! turbulence on ratqs, i.e. on the width of the
8! total water subgrid distribution.
9!=============================================
10
11IMPLICIT NONE
12
13! Include
14!=============================================
15    INCLUDE "YOETHF.h"
16    INCLUDE "YOMCST.h"
17
18
19      CONTAINS
20
21
22!========================================================================   
23SUBROUTINE calcratqs_inter(klon,klev,iflag_ratqs,pdtphys, &
24           ratqsbas, wake_deltaq, wake_s, q_seri,qtc_cv, sigt_cv,     &
25           ratqs_inter)
26USE ioipsl_getin_p_mod, ONLY : getin_p
27implicit none
28
29!========================================================================
30! L. d'Alençon, 25/02/2021
31! Cette subroutine calcule une valeur de ratqsbas interactive dépendant de la présence de poches froides dans l'environnement.
32! Elle est appelée par la subroutine calcratqs lorsque iflag_ratqs = 10.
33!========================================================================
34
35! Declarations
36
37
38LOGICAL, SAVE :: first = .TRUE.
39!$OMP THREADPRIVATE(first)
40! Input
41integer,intent(in) :: klon,klev,iflag_ratqs
42real,intent(in) :: pdtphys,ratqsbas
43real, dimension(klon,klev),intent(in) :: wake_deltaq, q_seri,qtc_cv, sigt_cv
44real, dimension(klon),intent(in) :: wake_s
45
46! Output
47real, dimension(klon,klev),intent(inout) :: ratqs_inter
48
49! local
50integer i,k
51real, dimension(klon,klev) :: wake_dq
52REAL, SAVE             :: a_ratqs_cv
53!$OMP THREADPRIVATE(a_ratqs_cv)
54REAL, SAVE             :: tau_ratqs_wake
55!$OMP THREADPRIVATE(tau_ratqs_wake)
56REAL, SAVE             :: a_ratqs_wake
57!$OMP THREADPRIVATE(a_ratqs_wake)
58real, dimension(klon) :: max_wake_dq, max_dqconv,max_sigt
59!-------------------------------------------------------------------------
60!  Caclul de ratqs_inter
61!-------------------------------------------------------------------------
62
63! 
64      if (first) then
65         tau_ratqs_wake = 3600. ! temps de relaxation de la variabilité
66         a_ratqs_wake = 3.    ! paramètre pilotant l'importance du terme dépendant des poches froides
67         a_ratqs_cv = 0.5
68         CALL getin_p('tau_ratqs_wake', tau_ratqs_wake)
69         CALL getin_p('a_ratqs_wake', a_ratqs_wake)
70         CALL getin_p('a_ratqs_cv', a_ratqs_cv)         
71         first=.false.
72      endif       
73      max_wake_dq(:) = 0.
74      max_dqconv (:) = 0
75      max_sigt(:)    = 0.
76      if (iflag_ratqs.eq.10) then
77        do k=1,klev
78          do i=1,klon
79           max_wake_dq(i) = max(abs(wake_deltaq(i,k)),max_wake_dq(i))
80           max_sigt(i) = max(abs(sigt_cv(i,k)),max_sigt(i))
81           max_dqconv(i) = max(abs(q_seri(i,k) - qtc_cv(i,k)),max_dqconv(i))
82          enddo
83        enddo
84     
85        do k=1,klev
86          do i=1,klon
87           ratqs_inter(i,k)= ratqs_inter(i,k)*exp(-pdtphys/tau_ratqs_wake) +   &
88           a_ratqs_wake*(max_wake_dq(i)*(wake_s(i)**0.5/(1.-wake_s(i))))*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1)   
89           if (ratqs_inter(i,k)<ratqsbas) then
90              ratqs_inter(i,k) = ratqsbas
91           endif       
92          enddo
93        enddo
94      endif
95
96      if (iflag_ratqs.eq.11) then
97        do k=1,klev
98          do i=1,klon
99           max_wake_dq(i) = max(abs(wake_deltaq(i,k)),max_wake_dq(i))
100           max_sigt(i) = max(abs(sigt_cv(i,k)),max_sigt(i))
101           max_dqconv(i) = max(abs(q_seri(i,k) - qtc_cv(i,k)),max_dqconv(i))
102          enddo
103        enddo
104     
105        do k=1,klev
106          do i=1,klon
107           ratqs_inter(i,k)= ratqs_inter(i,k)*exp(-pdtphys/tau_ratqs_wake) +   &
108           a_ratqs_wake*(max_wake_dq(i)*(wake_s(i)**0.5/(1.-wake_s(i))))*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1)  +   &
109           a_ratqs_cv*max_dqconv(i)*max_sigt(i)*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1)
110!           if (ratqs_inter(i,k)>0) then
111!              ratqs_inter(i,k) = abs(q_seri(i,k) - qtc_cv(i,k))
112!           endif       
113          enddo
114        enddo
115      endif
116return
117end         
118
119
120!------------------------------------------------------------------
121SUBROUTINE calcratqs_oro(klon,klev,qsat,temp,pplay,paprs,ratqs_oro)
122
123! Etienne Vignon, November 2021: effect of subgrid orography on ratqs
124
125USE phys_state_var_mod, ONLY: zstd
126USE phys_state_var_mod, ONLY: pctsrf
127USE indice_sol_mod, only: nbsrf, is_lic, is_ter
128
129IMPLICIT NONE
130
131! Declarations
132!--------------
133
134! INPUTS
135
136INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
137INTEGER, INTENT(IN) :: klev                       ! number of vertical layers
138REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat    ! saturation specific humidity [kg/kg]
139REAL, DIMENSION(klon,klev), INTENT(IN) :: temp    ! air temperature [K]
140REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay    ! air pressure, layer's center [Pa]
141REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs    ! air pressure, lower inteface [Pa]
142
143! OUTPUTS
144
145REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_oro ! ratqs profile due to subgrid orography
146
147
148! LOCAL
149
150INTEGER :: i,k
151REAL, DIMENSION(klon) :: orogradT,xsi0
152REAL, DIMENSION (klon,klev) :: zlay
153REAL :: Lvs, temp0
154
155
156! Calculation of the near-surface temperature gradient along the topography
157!--------------------------------------------------------------------------
158
159! at the moment, we fix it at a constant value (moist adiab. lapse rate)
160
161orogradT(:)=-6.5/1000. ! K/m
162
163! Calculation of near-surface surface ratqs
164!-------------------------------------------
165
166DO i=1,klon
167    temp0=temp(i,1)
168    IF (temp0 .LT. RTT) THEN
169        Lvs=RLSTT
170    ELSE
171        Lvs=RLVTT
172    ENDIF
173    xsi0(i)=zstd(i)*ABS(orogradT(i))*Lvs/temp0/temp0/RV
174    ratqs_oro(i,1)=xsi0(i)
175END DO
176
177! Vertical profile of ratqs assuming an exponential decrease with height
178!------------------------------------------------------------------------
179     
180! calculation of geop. height AGL       
181zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) &
182           *(paprs(:,1)-pplay(:,1))/RG
183
184DO k=2,klev
185   DO i = 1, klon
186      zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) &
187               /paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG
188               
189      ratqs_oro(i,k)=MAX(0.0,pctsrf(i,is_ter)*xsi0(i)*exp(-zlay(i,k)/MAX(zstd(i),1.)))   
190    END DO
191END DO
192
193
194
195
196END SUBROUTINE calcratqs_oro
197
198!=============================================
199
200SUBROUTINE calcratqs_hetero(klon,klev,t2m,q2m,temp,q,pplay,paprs,ratqs_hetero)
201
202! Etienne Vignon, November 2021
203! Effect of subgrid surface heterogeneities on ratqs
204
205USE phys_local_var_mod, ONLY: s_pblh
206USE phys_state_var_mod, ONLY: pctsrf
207USE indice_sol_mod
208USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF
209
210IMPLICIT NONE
211
212include "YOMCST.h"
213
214! INPUTS
215
216
217INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
218INTEGER, INTENT(IN) :: klev                       ! number of vertical layers
219REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: t2m    ! 2m temperature for each tile [K]
220REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: q2m    ! 2m specific humidity for each tile [kg/kg]
221REAL, DIMENSION(klon,klev), INTENT(IN) :: temp    ! air temperature [K]
222REAL, DIMENSION(klon,klev), INTENT(IN) :: q       ! specific humidity [kg/kg]
223REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay   ! air pressure, layer's center [Pa]
224REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! air pressure, lower inteface [Pa]
225
226! OUTPUTS
227
228REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_hetero ! ratsq profile due to surface heterogeneities
229
230
231INTEGER :: i,k,nsrf
232REAL :: ratiom, qsat2m, dqsatdT
233REAL, DIMENSION(klon) :: xsi0
234REAL, DIMENSION (klon,klev) :: zlay
235
236
237
238! Calculation of near-surface surface ratqs
239!-------------------------------------------
240
241DO i=1,klon
242   
243    ratiom=0.
244    xsi0(i)=0.
245   
246    DO nsrf=1,nbsrf
247    CALL CALC_QSAT_ECMWF(t2m(i,nsrf),q2m(i,nsrf),paprs(i,1),RTT,0,.false.,qsat2m,dqsatdT)
248    ratiom=ratiom+pctsrf(i,nsrf)*(q2m(i,nsrf)/qsat2m)
249    xsi0(i)=xsi0(i)+pctsrf(i,nsrf)*((q2m(i,nsrf)/qsat2m-ratiom)**2)
250    END DO
251   
252    xsi0(i)=sqrt(xsi0(i))/(ratiom+1E-6)
253END DO
254
255
256
257! Vertical profile of ratqs assuming an exponential decrease with height
258!------------------------------------------------------------------------
259       
260! calculation of geop. height AGL
261
262zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) &
263           *(paprs(:,1)-pplay(:,1))/RG
264ratqs_hetero(:,1)=xsi0(:)
265
266DO k=2,klev
267   DO i = 1, klon
268      zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) &
269               /paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG
270               
271      ratqs_hetero(i,k)=MAX(xsi0(i)*exp(-zlay(i,k)/(s_pblh(i)+1.0)),0.0)   
272    END DO
273END DO
274
275END SUBROUTINE calcratqs_hetero
276
277!=============================================
278
279SUBROUTINE calcratqs_tke(klon,klev,pdtphys,temp,q,qsat,pplay,paprs,tke,tke_dissip,lmix,wprime,ratqs_tke)
280
281! References:
282!
283! Etienne Vignon: effect of subgrid turbulence on ratqs
284!
285! Field, P.R., Hill, A., Furtado, K., Korolev, A., 2014b. Mixed-phase clouds in a turbulent environment. Part
286! 2: analytic treatment. Q. J. R. Meteorol. Soc. 21, 2651–2663. https://doi.org/10.1002/qj.2175.
287!
288! Furtado, K., Field, P.R., Boutle, I.A., Morcrette, C.R., Wilkinson, J., 2016. A physically-based, subgrid
289! parametrization for the production and maintenance of mixed-phase clouds in a general circulation
290! model. J. Atmos. Sci. 73, 279–291. https://doi.org/10.1175/JAS-D-15-0021.
291
292USE phys_local_var_mod, ONLY: omega
293
294IMPLICIT NONE
295
296! INPUTS
297
298INTEGER, INTENT(IN) :: klon                             ! number of horizontal grid points
299INTEGER, INTENT(IN) :: klev                             ! number of vertical layers
300REAL, INTENT(IN) :: pdtphys                             ! physics time step [s]
301REAL, DIMENSION(klon,klev), INTENT(IN) :: temp          ! air temperature [K]
302REAL, DIMENSION(klon,klev), INTENT(IN) :: q             ! specific humidity [kg/kg]
303REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat          ! saturation specific humidity [kg/kg]
304REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay         ! air pressure, layer's center [Pa]
305REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs       ! air pressure, lower inteface [Pa]
306REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke         ! Turbulent Kinetic Energy [m2/s2]
307REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke_dissip  ! Turbulent Kinetic Energy Dissipation rate [m2/s3]
308REAL, DIMENSION(klon,klev+1), INTENT(IN) :: lmix  ! Turbulent mixing length
309REAL, DIMENSION(klon,klev+1), INTENT(IN) :: wprime      ! Turbulent vertical velocity scale [m/s]
310
311! OUTPUTS
312
313REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_tke  ! ratsq profile due to subgrid TKE
314
315! LOCAL
316INTEGER :: i, k
317REAL :: AA, DD, NW, AAprime, VARLOG,rho,Lvs,taue,lhomo,dissmin,maxvarlog
318REAL, DIMENSION(klon,klev) :: sigmaw,w
319REAL, PARAMETER :: C0=10.0
320REAL, PARAMETER :: lmin=0.001
321REAL, PARAMETER :: ratqsmin=1E-6
322REAL, PARAMETER :: ratqsmax=0.5
323
324
325! Calculation of large scale and turbulent vertical velocities
326!---------------------------------------------------------------
327
328DO k=1,klev
329    DO i=1,klon
330        rho=pplay(i,k)/temp(i,k)/RD
331        w(i,k)=-rho*RG*omega(i,k)
332        sigmaw(i,k)=0.5*(wprime(i,k+1)+wprime(i,k)) ! turbulent vertical velocity at the middle of model layers.
333    END DO
334END DO
335
336! Calculation of ratqs
337!---------------------------------------------------------------
338ratqs_tke(:,1)=ratqsmin ! set to a very low value to avoid division by 0 in order parts
339                        ! of the code
340DO k=2,klev ! we start from second model level since TKE is not defined at k=1
341    DO i=1,klon
342
343       IF (temp(i,k) .LT. RTT) THEN
344           Lvs=RLSTT
345       ELSE
346           Lvs=RLVTT
347       ENDIF
348       dissmin=0.01*(0.5*(tke(i,k)+tke(i,k+1))/pdtphys)
349       maxvarlog=LOG(1.0+ratqsmax**2)! to prevent ratqs from exceeding an arbitrary threshold value
350       AA=RG*(Lvs/(RCPD*temp(i,k)*temp(i,k)*RV) - 1./(RD*temp(i,k)))
351       lhomo=MAX(0.5*(lmix(i,k)+lmix(i,k+1)),lmin)
352       taue=(lhomo*lhomo/MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin))**(1./3) ! Fields et al. 2014
353       DD=1.0/taue
354       NW=(sigmaw(i,k)**2)*SQRT(2./(C0*MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin)))
355       AAprime=AA*NW
356       VARLOG=AAprime/2./DD
357       VARLOG=MIN(VARLOG,maxvarlog)
358       ratqs_tke(i,k)=SQRT(MAX(EXP(VARLOG)-1.0,ratqsmin))
359       END DO
360END DO
361END SUBROUTINE calcratqs_tke
362
363
364
365
366
367END MODULE calcratqs_multi_mod
Note: See TracBrowser for help on using the repository browser.