Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

calcul_fluxs_mod.f90

Last change on this file was 5274, checked in by abarral, 9 hours ago
Replace yomcst.h by existing module
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 14.6 KB

Line
1	!
2	! $Id: calcul_fluxs_mod.f90 5274 2024-10-25 13:41:23Z abarral $
3	!
4	MODULE calcul_fluxs_mod
5
6	IMPLICIT NONE
7
8	CONTAINS
9	SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
10	tsurf, p1lay, cal, beta, cdragh, cdragq, ps, &
11	precip_rain, precip_snow, snow, qsurf, &
12	radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, gustiness, &
13	fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, &
14	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
15	sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
16
17
18	USE dimphy, ONLY : klon
19	USE indice_sol_mod
20	use sens_heat_rain_m, only: sens_heat_rain
21	USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO &
22	, RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA &
23	, R_ecc, R_peri, R_incl &
24	, RA, RG, R1SA &
25	, RSIGMA &
26	, R, RMD, RMV, RD, RV, RCPD &
27	, RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12 &
28	, RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w &
29	, RCW, RCS &
30	, RLVTT, RLSTT, RLMLT, RTT, RATM &
31	, RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS &
32	, RALPD, RBETD, RGAMD
33
34	INCLUDE "clesphys.h"
35
36	! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
37	! une temperature de surface (au cas ou ok_veget = false)
38	!
39	! L. Fairhead 4/2000
40	!
41	! input:
42	! knon nombre de points a traiter
43	! nisurf surface a traiter
44	! tsurf temperature de surface
45	! p1lay pression 1er niveau (milieu de couche)
46	! cal capacite calorifique du sol
47	! beta evap reelle
48	! cdragh coefficient d'echange temperature
49	! cdragq coefficient d'echange evaporation
50	! ps pression au sol
51	! precip_rain precipitations liquides
52	! precip_snow precipitations solides
53	! snow champs hauteur de neige
54	! runoff runoff en cas de trop plein
55	! petAcoef coeff. A de la resolution de la CL pour t
56	! peqAcoef coeff. A de la resolution de la CL pour q
57	! petBcoef coeff. B de la resolution de la CL pour t
58	! peqBcoef coeff. B de la resolution de la CL pour q
59	! radsol rayonnement net aus sol (LW + SW)
60	! dif_grnd coeff. diffusion vers le sol profond
61	!
62	! output:
63	! tsurf_new temperature au sol
64	! qsurf humidite de l'air au dessus du sol
65	! fluxsens flux de chaleur sensible
66	! fluxlat flux de chaleur latente
67	! dflux_s derivee du flux de chaleur sensible / Ts
68	! dflux_l derivee du flux de chaleur latente / Ts
69	! sens_prec_liq flux sensible li� aux echanges de precipitations liquides
70	! sens_prec_sol precipitations solides
71	! lat_prec_liq flux latent li� aux echanges de precipitations liquides
72	! lat_prec_sol precipitations solides
73
74	INCLUDE "YOETHF.h"
75	INCLUDE "FCTTRE.h"
76
77	! Parametres d'entree
78	!****************************************************************************************
79	INTEGER, INTENT(IN) :: knon, nisurf
80	REAL , INTENT(IN) :: dtime
81	REAL, DIMENSION(klon), INTENT(IN) :: petAcoef, peqAcoef
82	REAL, DIMENSION(klon), INTENT(IN) :: petBcoef, peqBcoef
83	REAL, DIMENSION(klon), INTENT(IN) :: ps, q1lay
84	REAL, DIMENSION(klon), INTENT(IN) :: tsurf, p1lay, cal, beta, cdragh,cdragq
85	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
86	REAL, DIMENSION(klon), INTENT(IN) :: radsol, dif_grnd
87	REAL, DIMENSION(klon), INTENT(IN) :: t1lay, u1lay, v1lay,gustiness
88	REAL, INTENT(IN) :: fqsat ! correction factor on qsat (generally 0.98 over salty water, 1 everywhere else)
89
90	real, intent(in), optional:: rhoa(:) ! (knon)
91	! density of moist air (kg / m3)
92
93	! Parametres entree-sorties
94	!****************************************************************************************
95	REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! snow pas utile
96
97	! Parametres sorties
98	!****************************************************************************************
99	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
100	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new, evap, fluxsens, fluxlat
101	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
102	REAL, intent(out), OPTIONAL:: sens_prec_liq(:), sens_prec_sol(:) ! (knon)
103	REAL, DIMENSION(klon), OPTIONAL :: lat_prec_liq, lat_prec_sol
104
105	! Variables locales
106	!****************************************************************************************
107	INTEGER :: i
108	REAL, DIMENSION(klon) :: zx_mh, zx_nh, zx_oh
109	REAL, DIMENSION(klon) :: zx_mq, zx_nq, zx_oq
110	REAL, DIMENSION(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat
111	REAL, DIMENSION(klon) :: zx_sl, zx_coefh, zx_coefq, zx_wind
112	REAL, DIMENSION(klon) :: d_ts
113	REAL :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
114	REAL :: qsat_new, q1_new
115	REAL, PARAMETER :: t_grnd = 271.35, t_coup = 273.15
116	REAL, PARAMETER :: max_eau_sol = 150.0
117	CHARACTER (len = 20) :: modname = 'calcul_fluxs'
118	LOGICAL :: fonte_neige
119	LOGICAL, SAVE :: check = .FALSE.
120	!$OMP THREADPRIVATE(check)
121
122	! End definition
123	!****************************************************************************************
124
125	IF (check) WRITE(,)'Entree ', modname,' surface = ',nisurf
126
127	IF (check) THEN
128	WRITE(,)' radsol (min, max)', &
129	MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
130	ENDIF
131
132	! Traitement neige et humidite du sol
133	!****************************************************************************************
134	!
135	!!$ WRITE(,)'test calcul_flux, surface ', nisurf
136	!!PB test
137	!!$ if (nisurf == is_oce) then
138	!!$ snow = 0.
139	!!$ qsol = max_eau_sol
140	!!$ else
141	!!$ where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
142	!!$ where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
143	!!$! snow = max(0.0, snow + (precip_snow - evap) * dtime)
144	!!$ where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
145	!!$ endif
146	!!$ IF (nisurf /= is_ter) qsol = max_eau_sol
147
148
149	!
150	! Initialisation
151	!****************************************************************************************
152	evap = 0.
153	fluxsens=0.
154	fluxlat=0.
155	dflux_s = 0.
156	dflux_l = 0.
157	if (PRESENT(lat_prec_liq)) lat_prec_liq = 0.
158	if (PRESENT(lat_prec_sol)) lat_prec_sol = 0.
159	!
160	! zx_qs = qsat en kg/kg
161	!****************************************************************************************
162	DO i = 1, knon
163	zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
164	IF (thermcep) THEN
165	zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
166	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
167	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
168	zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
169	zx_qs=MIN(0.5,zx_qs)
170	zcor=1./(1.-retv*zx_qs)
171	zx_qs=zx_qs*zcor
172	zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
173	/RLVTT / zx_pkh(i)
174	ELSE
175	IF (tsurf(i).LT.t_coup) THEN
176	zx_qs = qsats(tsurf(i)) / ps(i)
177	zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
178	/ zx_pkh(i)
179	ELSE
180	zx_qs = qsatl(tsurf(i)) / ps(i)
181	zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
182	/ zx_pkh(i)
183	ENDIF
184	ENDIF
185	zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
186	zx_qsat(i) = zx_qs
187	zx_wind(i)=min_wind_speed+SQRT(gustiness(i)+u1lay(i)2+v1lay(i)2)
188	zx_coefh(i) = cdragh(i) * zx_wind(i) * p1lay(i)/(RD*t1lay(i))
189	zx_coefq(i) = cdragq(i) * zx_wind(i) * p1lay(i)/(RD*t1lay(i))
190	! zx_wind(i)=min_wind_speed+SQRT(gustiness(i)+u1lay(i)2+v1lay(i)2) &
191	! * p1lay(i)/(RD*t1lay(i))
192	! zx_coefh(i) = cdragh(i) * zx_wind(i)
193	! zx_coefq(i) = cdragq(i) * zx_wind(i)
194	ENDDO
195
196
197	! === Calcul de la temperature de surface ===
198	! zx_sl = chaleur latente d'evaporation ou de sublimation
199	!****************************************************************************************
200
201	DO i = 1, knon
202	zx_sl(i) = RLVTT
203	IF (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
204	ENDDO
205
206
207	DO i = 1, knon
208	! Q
209	zx_oq(i) = 1. - (beta(i) * zx_coefq(i) * peqBcoef(i) * dtime)
210	zx_mq(i) = beta(i) * zx_coefq(i) * &
211	(peqAcoef(i) - &
212	! conv num avec precedente version
213	fqsat * zx_qsat(i) + fqsat * zx_dq_s_dt(i) * tsurf(i)) &
214	! fqsat * ( zx_qsat(i) - zx_dq_s_dt(i) * tsurf(i)) ) &
215	/ zx_oq(i)
216	zx_nq(i) = beta(i) * zx_coefq(i) * (- fqsat * zx_dq_s_dt(i)) &
217	/ zx_oq(i)
218
219	! H
220	zx_oh(i) = 1. - (zx_coefh(i) * petBcoef(i) * dtime)
221	zx_mh(i) = zx_coefh(i) * petAcoef(i) / zx_oh(i)
222	zx_nh(i) = - (zx_coefh(i) * RCPD * zx_pkh(i))/ zx_oh(i)
223
224	! Tsurface
225	tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
226	(radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
227	+ dif_grnd(i) * t_grnd * dtime)/ &
228	( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
229	zx_nh(i) + zx_sl(i) * zx_nq(i)) &
230	+ dtime * dif_grnd(i))
231
232	!
233	! Y'a-t-il fonte de neige?
234	!
235	! fonte_neige = (nisurf /= is_oce) .AND. &
236	! & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
237	! & .AND. (tsurf_new(i) >= RTT)
238	! if (fonte_neige) tsurf_new(i) = RTT
239	d_ts(i) = tsurf_new(i) - tsurf(i)
240	! zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
241	! zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
242
243	!== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
244	!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
245	evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
246	fluxlat(i) = - evap(i) * zx_sl(i)
247	fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
248
249	! Derives des flux dF/dTs (W m-2 K-1):
250	dflux_s(i) = zx_nh(i)
251	dflux_l(i) = (zx_sl(i) * zx_nq(i))
252
253	! Nouvelle valeure de l'humidite au dessus du sol
254	qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
255	q1_new = peqAcoef(i) - peqBcoef(i)evap(i)dtime
256	qsurf(i)=q1_new(1.-beta(i)) + beta(i)qsat_new
257	!
258	! en cas de fonte de neige
259	!
260	! if (fonte_neige) then
261	! bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
262	! & dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
263	! & RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
264	! bilan_f = max(0., bilan_f)
265	! fq_fonte = bilan_f / zx_sl(i)
266	! snow(i) = max(0., snow(i) - fq_fonte * dtime)
267	! qsol(i) = qsol(i) + (fq_fonte * dtime)
268	! endif
269	!!$ if (nisurf == is_ter) &
270	!!$ & run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
271	!!$ qsol(i) = min(qsol(i), max_eau_sol)
272	!
273	! calcul de l'enthalpie des precipitations liquides et solides
274	if (PRESENT(sens_prec_liq)) sens_prec_liq(i) &
275	= - sens_heat_rain(precip_rain(i) + precip_snow(i), t1lay(i), &
276	q1lay(i), rhoa(i), rlvtt, tsurf_new(i), ps(i))
277	if (PRESENT(sens_prec_sol)) sens_prec_sol(i) = 0.
278	! On calcule par rapport a T=0
279	!! sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * precip_rain(i)
280	!! sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * precip_snow(i)
281
282	if (PRESENT(lat_prec_liq)) &
283	lat_prec_liq(i) = precip_rain(i) * (RLVTT - RLVTT)
284	if (PRESENT(lat_prec_sol)) &
285	lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
286	ENDDO
287
288
289	!**************************************************************************
290	!
291	END SUBROUTINE calcul_fluxs
292	!
293	!****************************************************************************************
294	!
295	SUBROUTINE calcul_flux_wind(knon, dtime, &
296	u0, v0, u1, v1, gustiness, cdrag_m, &
297	AcoefU, AcoefV, BcoefU, BcoefV, &
298	p1lay, t1lay, &
299	flux_u1, flux_v1)
300
301	USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO &
302	, RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA &
303	, R_ecc, R_peri, R_incl &
304	, RA, RG, R1SA &
305	, RSIGMA &
306	, R, RMD, RMV, RD, RV, RCPD &
307	, RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12 &
308	, RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w &
309	, RCW, RCS &
310	, RLVTT, RLSTT, RLMLT, RTT, RATM &
311	, RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS &
312	, RALPD, RBETD, RGAMD
313	USE dimphy
314
315	INCLUDE "clesphys.h"
316
317	! Input arguments
318	!****************************************************************************************
319	INTEGER, INTENT(IN) :: knon
320	REAL, INTENT(IN) :: dtime
321	REAL, DIMENSION(klon), INTENT(IN) :: u0, v0 ! u and v at niveau 0
322	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness ! u and v at niveau 1
323	REAL, DIMENSION(klon), INTENT(IN) :: cdrag_m ! cdrag pour momentum
324	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
325	REAL, DIMENSION(klon), INTENT(IN) :: p1lay ! pression 1er niveau (milieu de couche)
326	REAL, DIMENSION(klon), INTENT(IN) :: t1lay ! temperature
327	! Output arguments
328	!****************************************************************************************
329	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1
330	REAL, DIMENSION(klon), INTENT(OUT) :: flux_v1
331
332	! Local variables
333	!****************************************************************************************
334	INTEGER :: i
335	REAL :: mod_wind, buf
336
337	!****************************************************************************************
338	! Calculate the surface flux
339	!
340	!****************************************************************************************
341	DO i=1,knon
342	mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))2 + (v1(i)-v0(i))2)
343	buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
344	flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
345	flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
346	END DO
347
348	END SUBROUTINE calcul_flux_wind
349	!
350	!****************************************************************************************
351	!
352	END MODULE calcul_fluxs_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format