source: LMDZ6/branches/blowing_snow/libf/phylmd/calcratqs_multi_mod.F90 @ 5322

Last change on this file since 5322 was 4072, checked in by evignon, 3 years ago

Mise à jour de la routine de nuages lscp
les principaux changements consistent en:

  • des corrections de bug (déclaration et division

par zero) pour que la routine compile avec les modifications d'OB

  • des restructurations pour que les fonctions de lscp_tools_mod

travaillent sur des tableaux et non des scalaires

  • une séparation des coefficients d'évaporation de la neige et de la pluie
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, DIMENSION(klon) :: xsi0, ratiom, qsat2m, dqsatdT
233REAL, DIMENSION (klon,klev) :: zlay
234
235
236
237! Calculation of near-surface surface ratqs
238!-------------------------------------------
239
240   
241    ratiom(:)=0.
242    xsi0(:)=0.
243   
244    DO nsrf=1,nbsrf
245    CALL CALC_QSAT_ECMWF(klon,t2m(:,nsrf),q2m(:,nsrf),paprs(:,1),RTT,0,.false.,qsat2m,dqsatdT)
246    ratiom(:)=ratiom(:)+pctsrf(:,nsrf)*(q2m(:,nsrf)/qsat2m(:))
247    xsi0(:)=xsi0(:)+pctsrf(:,nsrf)*((q2m(:,nsrf)/qsat2m(:)-ratiom(:))**2)
248    END DO
249   
250    xsi0(:)=sqrt(xsi0(:))/(ratiom(:)+1E-6)
251
252
253
254! Vertical profile of ratqs assuming an exponential decrease with height
255!------------------------------------------------------------------------
256       
257! calculation of geop. height AGL
258
259zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) &
260           *(paprs(:,1)-pplay(:,1))/RG
261ratqs_hetero(:,1)=xsi0(:)
262
263DO k=2,klev
264   DO i = 1, klon
265      zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) &
266               /paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG
267               
268      ratqs_hetero(i,k)=MAX(xsi0(i)*exp(-zlay(i,k)/(s_pblh(i)+1.0)),0.0)   
269    END DO
270END DO
271
272END SUBROUTINE calcratqs_hetero
273
274!=============================================
275
276SUBROUTINE calcratqs_tke(klon,klev,pdtphys,temp,q,qsat,pplay,paprs,tke,tke_dissip,lmix,wprime,ratqs_tke)
277
278! References:
279!
280! Etienne Vignon: effect of subgrid turbulence on ratqs
281!
282! Field, P.R., Hill, A., Furtado, K., Korolev, A., 2014b. Mixed-phase clouds in a turbulent environment. Part
283! 2: analytic treatment. Q. J. R. Meteorol. Soc. 21, 2651–2663. https://doi.org/10.1002/qj.2175.
284!
285! Furtado, K., Field, P.R., Boutle, I.A., Morcrette, C.R., Wilkinson, J., 2016. A physically-based, subgrid
286! parametrization for the production and maintenance of mixed-phase clouds in a general circulation
287! model. J. Atmos. Sci. 73, 279–291. https://doi.org/10.1175/JAS-D-15-0021.
288
289USE phys_local_var_mod, ONLY: omega
290
291IMPLICIT NONE
292
293! INPUTS
294
295INTEGER, INTENT(IN) :: klon                             ! number of horizontal grid points
296INTEGER, INTENT(IN) :: klev                             ! number of vertical layers
297REAL, INTENT(IN) :: pdtphys                             ! physics time step [s]
298REAL, DIMENSION(klon,klev), INTENT(IN) :: temp          ! air temperature [K]
299REAL, DIMENSION(klon,klev), INTENT(IN) :: q             ! specific humidity [kg/kg]
300REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat          ! saturation specific humidity [kg/kg]
301REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay         ! air pressure, layer's center [Pa]
302REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs       ! air pressure, lower inteface [Pa]
303REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke         ! Turbulent Kinetic Energy [m2/s2]
304REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke_dissip  ! Turbulent Kinetic Energy Dissipation rate [m2/s3]
305REAL, DIMENSION(klon,klev+1), INTENT(IN) :: lmix  ! Turbulent mixing length
306REAL, DIMENSION(klon,klev+1), INTENT(IN) :: wprime      ! Turbulent vertical velocity scale [m/s]
307
308! OUTPUTS
309
310REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_tke  ! ratsq profile due to subgrid TKE
311
312! LOCAL
313INTEGER :: i, k
314REAL :: AA, DD, NW, AAprime, VARLOG,rho,Lvs,taue,lhomo,dissmin,maxvarlog
315REAL, DIMENSION(klon,klev) :: sigmaw,w
316REAL, PARAMETER :: C0=10.0
317REAL, PARAMETER :: lmin=0.001
318REAL, PARAMETER :: ratqsmin=1E-6
319REAL, PARAMETER :: ratqsmax=0.5
320
321
322! Calculation of large scale and turbulent vertical velocities
323!---------------------------------------------------------------
324
325DO k=1,klev
326    DO i=1,klon
327        rho=pplay(i,k)/temp(i,k)/RD
328        w(i,k)=-rho*RG*omega(i,k)
329        sigmaw(i,k)=0.5*(wprime(i,k+1)+wprime(i,k)) ! turbulent vertical velocity at the middle of model layers.
330    END DO
331END DO
332
333! Calculation of ratqs
334!---------------------------------------------------------------
335ratqs_tke(:,1)=ratqsmin ! set to a very low value to avoid division by 0 in order parts
336                        ! of the code
337DO k=2,klev ! we start from second model level since TKE is not defined at k=1
338    DO i=1,klon
339
340       IF (temp(i,k) .LT. RTT) THEN
341           Lvs=RLSTT
342       ELSE
343           Lvs=RLVTT
344       ENDIF
345       dissmin=0.01*(0.5*(tke(i,k)+tke(i,k+1))/pdtphys)
346       maxvarlog=LOG(1.0+ratqsmax**2)! to prevent ratqs from exceeding an arbitrary threshold value
347       AA=RG*(Lvs/(RCPD*temp(i,k)*temp(i,k)*RV) - 1./(RD*temp(i,k)))
348       lhomo=MAX(0.5*(lmix(i,k)+lmix(i,k+1)),lmin)
349       taue=(lhomo*lhomo/MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin))**(1./3) ! Fields et al. 2014
350       DD=1.0/taue
351       NW=(sigmaw(i,k)**2)*SQRT(2./(C0*MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin)))
352       AAprime=AA*NW
353       VARLOG=AAprime/2./DD
354       VARLOG=MIN(VARLOG,maxvarlog)
355       ratqs_tke(i,k)=SQRT(MAX(EXP(VARLOG)-1.0,ratqsmin))
356       END DO
357END DO
358END SUBROUTINE calcratqs_tke
359
360
361
362
363
364END MODULE calcratqs_multi_mod
Note: See TracBrowser for help on using the repository browser.