Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

surf_land_bucket_hetero_mod.F90

Last change on this file was 5629, checked in by Laurent Fairhead, 6 weeks ago
Adaptation of heterogeneous continental surfaces to isotopes
File size: 16.7 KB

Line
1	!
2	MODULE surf_land_bucket_hetero_mod
3	!
4	! 2025/04 A. Maison (adapted from surf_land_bucket_mod)
5	! Surface land bucket module with heterogeneous continental sub-surfaces
6	! This module is used when no external land model is choosen and iflag_hetero_surf = 1 or 2.
7	!
8	IMPLICIT NONE
9
10	CONTAINS
11
12	SUBROUTINE surf_land_bucket_hetero(itime, jour, knon, knindex, debut, dtime,&
13	tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, &
14	spechum, petAcoef, peqAcoef, petBcoef, peqBcoef, pref, plev, &
15	u1, v1, gustiness, rugoro, swnet, lwnet, &
16	snow, qsol, agesno, tsoil, &
17	qsurf, z0m, z0h, alb1_new, alb2_new, evap, &
18	fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l, &
19	tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
20	cdragm_tersrf, cdragh_tersrf, &
21	swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf)
22
23	USE clesphys_mod_h
24	USE dimsoil_mod_h, ONLY: nsoilmx
25	USE yomcst_mod_h, ONLY: RD, RG, RCPD, RSIGMA
26	USE compbl_mod_h
27	USE dimpft_mod_h
28	USE limit_read_mod
29	USE surface_data
30	USE fonte_neige_mod
31	USE calcul_fluxs_mod
32	USE cpl_mod
33	USE dimphy
34	USE geometry_mod, ONLY: longitude,latitude
35	USE mod_grid_phy_lmdz
36	USE mod_phys_lmdz_para
37	USE indice_sol_mod
38	USE phys_state_var_mod, ONLY: frac_tersrf, ratio_z0m_z0h_tersrf, z0m_tersrf, &
39	albedo_tersrf, beta_tersrf, inertie_tersrf, &
40	hcond_tersrf
41	USE surf_param_mod, ONLY: eff_surf_param, average_surf_var
42	USE cdrag_mod
43
44	#ifdef ISO
45	use infotrac_phy, ONLY: niso
46	#endif
47
48	!****************************************************************************************
49	! Bucket calculations for surface.
50	!****************************************************************************************
51	!
52	! Input variables
53	!****************************************************************************************
54	INTEGER, INTENT(IN) :: itime, jour, knon
55	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
56	LOGICAL, INTENT(IN) :: debut
57	REAL, INTENT(IN) :: dtime
58	REAL, DIMENSION(klon), INTENT(IN) :: tsurf
59	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
60	REAL, DIMENSION(klon), INTENT(IN) :: tq_cdrag
61	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
62	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
63	REAL, DIMENSION(klon), INTENT(IN) :: petAcoef, peqAcoef
64	REAL, DIMENSION(klon), INTENT(IN) :: petBcoef, peqBcoef
65	REAL, DIMENSION(klon), INTENT(IN) :: pref
66	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
67	REAL, DIMENSION(klon), INTENT(IN) :: rugoro
68	REAL, DIMENSION(klon), INTENT(IN) :: swnet, lwnet
69	REAL, DIMENSION(klon, nbtersrf), INTENT(IN) :: tsurf_tersrf
70
71
72	! In/Output variables
73	!****************************************************************************************
74	REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
75	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
76	REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
77	REAL, DIMENSION(klon), INTENT(INOUT) :: plev
78	REAL, DIMENSION(klon, nsoilmx, nbtersrf), INTENT(INOUT) :: tsoil_tersrf
79
80
81	! Output variables
82	!****************************************************************************************
83	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
84	REAL, DIMENSION(klon), INTENT(OUT) :: z0m, z0h
85	REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new, alb2_new
86	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
87	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
88	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
89	!
90	REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: tsurf_new_tersrf
91	REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: qsurf_tersrf
92	REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: cdragm_tersrf, cdragh_tersrf
93	REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: swnet_tersrf, lwnet_tersrf
94	REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: fluxsens_tersrf, fluxlat_tersrf
95
96	#ifdef ISO
97	! REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap
98	! REAL, DIMENSION(klon), INTENT(OUT) :: h1
99	! REAL, DIMENSION(niso,klon), INTENT(OUT) :: xtrunoff_diag
100	REAL, DIMENSION(klon) :: runoff_diag
101	! REAL, DIMENSION(niso,klon), INTENT(IN) :: Rland_ice
102	#endif
103	! Local variables
104	!****************************************************************************************
105	REAL, DIMENSION(klon) :: soilcap, soilflux
106	REAL, DIMENSION(klon) :: cal, beta, dif_grnd, capsol
107	REAL, DIMENSION(klon) :: alb_neig, alb_lim, icesub
108	REAL, DIMENSION(klon) :: zfra
109	REAL, DIMENSION(klon) :: radsol
110	REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay
111	REAL, DIMENSION(klon) :: dummy_riverflow,dummy_coastalflow
112	INTEGER :: i, j, k
113	#ifdef ISO
114	INTEGER :: ixt
115	REAL, PARAMETER :: t_coup = 273.15
116	REAL, DIMENSION(klon) :: fq_fonte_diag
117	REAL, DIMENSION(klon) :: fqfonte_diag
118	REAL, DIMENSION(klon) :: snow_evap_diag
119	REAL, DIMENSION(klon) :: fqcalving_diag
120	REAL :: max_eau_sol_diag
121	REAL, DIMENSION(klon) :: run_off_lic_diag
122	REAL :: coeff_rel_diag
123	#endif
124	!
125	REAL, DIMENSION(klon) :: zlev, geop1, speed, pblh, ri_in, sst
126	REAL, DIMENSION(klon) :: beta_eff, inertie_eff, conv_ratio_eff
127	REAL, DIMENSION(klon) :: meansqT
128	REAL, DIMENSION(klon, nbtersrf) :: z0h_tersrf, emis_tersrf, conv_ratio_tersrf
129	REAL, DIMENSION(klon, nbtersrf) :: evap_tersrf
130	REAL, DIMENSION(klon, nbtersrf) :: dflux_s_tersrf, dflux_l_tersrf
131	REAL, DIMENSION(klon, nbtersrf) :: radsol_tersrf
132	REAL, DIMENSION(klon, nbtersrf) :: zri_tersrf
133	REAL, PARAMETER :: klon_1D = 1
134
135	!
136	!****************************************************************************************
137
138	! * Calculations common to the two flag values *
139
140	! average albedo
141	alb_lim = eff_surf_param(klon, nbtersrf, albedo_tersrf, frac_tersrf, 'ARI')
142
143	! suppose zero surface speed
144	u0(:)=0.0
145	v0(:)=0.0
146	u1_lay(:) = u1(:) - u0(:)
147	v1_lay(:) = v1(:) - v0(:)
148	speed(:) = (u1_lay(:)2 + v1_lay(:)2)**0.5
149	!
150	geop1(1:knon) = RD * temp_air(1:knon) / (0.5*(pref(1:knon)+p1lay(1:knon))) &
151	* (pref(1:knon)-p1lay(1:knon))
152	zlev(1:knon) = (plev(1:knon)RDtemp_air(1:knon)/(pref(1:knon)*RG))/2.
153	!
154	! compute roughness lengths
155	DO i=1, knon
156	DO j=1, nbtersrf
157	IF (ratio_z0m_z0h_tersrf(i,j) .NE. 0.) THEN
158	z0h_tersrf(i,j) = z0m_tersrf(i,j) / ratio_z0m_z0h_tersrf(i,j)
159	ELSE
160	z0h_tersrf(i,j) = 1.E-12
161	ENDIF
162	ENDDO
163	ENDDO
164
165	z0m = eff_surf_param(klon, nbtersrf, z0m_tersrf, frac_tersrf, 'CDN', zlev)
166	z0h = eff_surf_param(klon, nbtersrf, z0h_tersrf, frac_tersrf, 'CDN', zlev)
167
168	DO i=1, knon
169	z0m(i) = MAX(1.5e-05,SQRT(z0m(i)2 + rugoro(i)2))
170	END DO
171
172	! compute the ratio to convert and print soil depths in meters (conv_ratio = (cond/cap)^0.5 and cap = I^2/cond)
173	DO j=1, nbtersrf
174	conv_ratio_tersrf(:,j) = hcond_tersrf(:,j)/inertie_tersrf(:,j)
175	ENDDO
176
177	!
178	! * Surface parameter aggregation *
179	!
180	IF (iflag_hetero_surf == 1) THEN
181	!* Calcultaion of fluxes
182
183	! calculate total absorbed radiance at surface
184	radsol(:) = 0.0
185	radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
186
187	! calculate constants (needeed for capsol and dif_grnd)
188	CALL calbeta(dtime, is_ter, knon, snow, qsol, beta, capsol, dif_grnd)
189	IF (type_veget=='betaclim') THEN
190	CALL calbeta_clim(knon,jour,latitude(knindex(1:knon)),beta)
191	ENDIF
192
193	! mean evapotranspiration coefficient
194	beta_eff = eff_surf_param(klon, nbtersrf, beta_tersrf, frac_tersrf, 'ARI')
195	beta = beta_eff
196
197	! calculate temperature, heat capacity and conduction flux in soil
198	IF (soil_model) THEN
199	inertie_eff = eff_surf_param(klon, nbtersrf, inertie_tersrf, frac_tersrf, 'ARI')
200	conv_ratio_eff = eff_surf_param(klon, nbtersrf, conv_ratio_tersrf, frac_tersrf, 'ARI')
201	!
202	CALL soil_hetero(dtime, is_ter, knon, snow, tsurf, qsol, &
203	& longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux, &
204	& inertie_eff, conv_ratio_eff)
205	!
206	DO i=1, knon
207	cal(i) = RCPD / soilcap(i)
208	radsol(i) = radsol(i) + soilflux(i)
209	END DO
210	ELSE
211	cal(:) = RCPD * capsol(:)
212	IF (klon_glo .EQ. 1) THEN
213	cal(:) = 0.
214	ENDIF
215	ENDIF
216
217	! calculate fluxes
218	CALL calcul_fluxs(knon, is_ter, dtime, &
219	tsurf, p1lay, cal, beta, tq_cdrag, tq_cdrag, pref, &
220	precip_rain, precip_snow, snow, qsurf, &
221	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
222	1.,petAcoef, peqAcoef, petBcoef, peqBcoef, &
223	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
224
225
226	!* Calculate snow height, run_off, age of snow
227	#ifdef ISO
228	DO i=1,knon
229	! initialisation:
230	fqfonte_diag(i) =0.0
231	fq_fonte_diag(i) =0.0
232	snow_evap_diag(i)=0.0
233	ENDDO !DO i=1,knon
234	#endif
235
236	CALL fonte_neige( knon, is_ter, knindex, dtime, &
237	tsurf, precip_rain, precip_snow, &
238	snow, qsol, tsurf_new, evap, icesub &
239	#ifdef ISO
240	& ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag &
241	& ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag &
242	#endif
243	& )
244
245	! calculate the age of snow
246	CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
247
248	WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
249
250	DO i=1, knon
251	zfra(i) = MAX(0.0,MIN(1.0, snow(i)/(snow(i)+10.0)))
252	alb_lim(i) = alb_neig(i) zfra(i) + alb_lim(i)(1.0-zfra(i))
253	END DO
254
255
256	!* Return albedo :
257	! alb1_new and alb2_new are here given the same values
258
259	alb1_new(:) = 0.0
260	alb2_new(:) = 0.0
261	alb1_new(1:knon) = alb_lim(1:knon)
262	alb2_new(1:knon) = alb_lim(1:knon)
263
264	!* Send to coupler
265	! The run-off from river and coast are not calculated in the bucket modele.
266	! For testing purpose of the coupled modele we put the run-off to zero.
267	IF (type_ocean=='couple') THEN
268	dummy_riverflow(:) = 0.0
269	dummy_coastalflow(:) = 0.0
270	CALL cpl_send_land_fields(itime, knon, knindex, &
271	dummy_riverflow, dummy_coastalflow)
272	ENDIF
273
274	!
275	! * Flux aggregation *
276	!
277	ELSE IF (iflag_hetero_surf == 2) THEN
278	! initialize output tables
279	evap_tersrf(:,:) = 0.
280	fluxsens_tersrf(:,:) = 0.
281	fluxlat_tersrf(:,:) = 0.
282	tsurf_new_tersrf(:,:) = 0.
283	dflux_s_tersrf(:,:) = 0.
284	dflux_l_tersrf(:,:) = 0.
285	radsol_tersrf(:,:) = 0.
286	swnet_tersrf(:,:) = 0.
287	lwnet_tersrf(:,:) = 0.
288	! hyp: surface emissivity = 1
289	emis_tersrf(:,:) = 1.
290
291	! * calculate total absorbed radiance at surface
292	DO j=1, nbtersrf
293	! SW
294	swnet_tersrf(klon_1D,j) = (1. - albedo_tersrf(klon_1D,j)) / (1. - alb_lim(klon_1D)) * swnet(klon_1D)
295	! LW
296	! first order
297	lwnet_tersrf(klon_1D,j) = lwnet(klon_1D) + 4. * emis_tersrf(klon_1D,j) * RSIGMA * tsurf(klon_1D)*3 &
298	(tsurf(klon_1D) - tsurf_tersrf(klon_1D,j))
299	ENDDO
300	! LW second order corrections
301	!- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)
302	IF (iflag_order2_sollw == 1) THEN
303	meansqT(:) = 0. ! as working buffer
304	!
305	DO j=1, nbtersrf
306	meansqT(klon_1D) = meansqT(klon_1D) + (tsurf_tersrf(klon_1D,j) - tsurf(klon_1D))*2 frac_tersrf(klon_1D,j)
307	ENDDO
308	DO j=1, nbtersrf
309	lwnet_tersrf(klon_1D,j) = lwnet_tersrf(klon_1D,j) + 6. * RSIGMA * tsurf(klon_1D)*2 (meansqT(klon_1D) - &
310	(tsurf(klon_1D) - tsurf_tersrf(klon_1D,j))**2)
311	ENDDO
312	ENDIF
313	! net radiation
314	radsol_tersrf(:,:) = swnet_tersrf(:,:) + lwnet_tersrf(:,:)
315
316	! * compute evapotranspiration coefficient
317	capsol(:) = 1.0/(2.5578E+06*0.15)
318	dif_grnd(:) = 0.
319
320	! unused variables in cdrag routine
321	pblh(:) = 0.
322	ri_in(:) = 0.
323	sst(:) = 0.
324
325	! Loop on sub-surfaces
326	DO j=1, nbtersrf
327	! * drag coefficients
328	CALL cdrag(knon, is_ter, speed, temp_air, spechum, geop1, pref, pblh, &
329	tsurf_tersrf(:,j), qsurf_tersrf(:,j), z0m_tersrf(:,j), z0h_tersrf(:,j), ri_in, 0, &
330	cdragm_tersrf(:,j), cdragh_tersrf(:,j), zri_tersrf(:,j), plev, precip_rain, sst, p1lay)
331
332	! * calculate temperature, heat capacity and conduction flux in soil
333	IF (soil_model) THEN
334	!
335	CALL soil_hetero(dtime, is_ter, knon, snow, tsurf_tersrf(:,j), qsol, &
336	longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil_tersrf(:,:,j), soilcap, soilflux, &
337	inertie_tersrf(:,j), conv_ratio_tersrf(:,j))
338	!
339	cal(:) = RCPD / soilcap(:)
340	radsol_tersrf(:,j) = radsol_tersrf(:,j) + soilflux(:)
341	!
342	ELSE
343	cal = RCPD * capsol
344	IF (klon_glo .EQ. 1) THEN
345	cal = 0.
346	ENDIF
347	ENDIF
348
349	! * calcultaion of fluxes
350	CALL calcul_fluxs(knon, is_ter, dtime, &
351	tsurf_tersrf(klon_1D,j), p1lay, cal, beta_tersrf(klon_1D,j), cdragh_tersrf(klon_1D,j), cdragh_tersrf(klon_1D,j), pref, &
352	precip_rain, precip_snow, snow, qsurf_tersrf(klon_1D,j), &
353	radsol_tersrf(klon_1D,j), dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
354	1.,petAcoef, peqAcoef, petBcoef, peqBcoef, &
355	tsurf_new_tersrf(klon_1D,j), evap_tersrf(klon_1D,j), fluxlat_tersrf(klon_1D,j), fluxsens_tersrf(klon_1D,j), &
356	dflux_s_tersrf(klon_1D,j), dflux_l_tersrf(klon_1D,j))
357
358	! if snow > 0
359	! calculate snow height, run_off, age of snow
360	CALL fonte_neige( knon, is_ter, knindex, dtime, &
361	tsurf_tersrf(:,j), precip_rain, precip_snow, &
362	snow, qsol, tsurf_new_tersrf(:,j), evap_tersrf(:,j), icesub &
363	#ifdef ISO
364	& ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag &
365	& ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag &
366	#endif
367	& )
368
369	ENDDO ! loop on sub-surfaces
370
371	! calculate the age of snow
372	CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
373	WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
374
375	DO i=1, knon
376	zfra(i) = MAX(0.0,MIN(1.0, snow(i)/(snow(i)+10.0)))
377	alb_lim(i) = alb_neig(i) zfra(i) + alb_lim(i)(1.0-zfra(i))
378	END DO
379
380	! return albedo :
381	! alb1_new and alb2_new are here given the same values
382	alb1_new(:) = 0.0
383	alb2_new(:) = 0.0
384	alb1_new(1:knon) = alb_lim(1:knon)
385	alb2_new(1:knon) = alb_lim(1:knon)
386
387	! send to coupler
388	! the run-off from river and coast are not calculated in the bucket modele.
389	! for testing purpose of the coupled modele we put the run-off to zero.
390	IF (type_ocean=='couple') THEN
391	dummy_riverflow(:) = 0.0
392	dummy_coastalflow(:) = 0.0
393	CALL cpl_send_land_fields(itime, knon, knindex, &
394	dummy_riverflow, dummy_coastalflow)
395	ENDIF
396
397	! * average of fluxes and surface variables
398	qsurf = average_surf_var(klon, nbtersrf, qsurf_tersrf, frac_tersrf, 'ARI')
399	tsurf_new = average_surf_var(klon, nbtersrf, tsurf_new_tersrf, frac_tersrf, 'ARI')
400	evap = average_surf_var(klon, nbtersrf, evap_tersrf, frac_tersrf, 'ARI')
401	fluxlat = average_surf_var(klon, nbtersrf, fluxlat_tersrf, frac_tersrf, 'ARI')
402	fluxsens = average_surf_var(klon, nbtersrf, fluxsens_tersrf, frac_tersrf, 'ARI')
403	dflux_l = average_surf_var(klon, nbtersrf, dflux_l_tersrf, frac_tersrf, 'ARI')
404	dflux_s = average_surf_var(klon, nbtersrf, dflux_s_tersrf, frac_tersrf, 'ARI')
405	DO k=1, nsoilmx
406	tsoil(:,k) = average_surf_var(klon, nbtersrf, tsoil_tersrf(:,k,:), frac_tersrf, 'ARI')
407	ENDDO
408
409	! order 2 correction to tsurf_new, for radiation computations (main atm effect of Ts)
410	IF (iflag_order2_sollw == 1) THEN
411	meansqT(:) = 0. ! as working buffer
412	DO j=1, nbtersrf
413	meansqT(klon_1D) = meansqT(klon_1D)+(tsurf_tersrf(klon_1D,j)-tsurf_new(klon_1D))*2 frac_tersrf(klon_1D,j)
414	ENDDO
415	tsurf_new(:) = tsurf_new(:)+1.5*meansqT(:)/tsurf_new(:)
416	ENDIF ! iflag_order2_sollw == 1
417
418	ENDIF ! iflag_hetero_surf
419	!
420	!* End
421	!
422	END SUBROUTINE surf_land_bucket_hetero
423	!
424	!****************************************************************************************
425	!
426	END MODULE surf_land_bucket_hetero_mod

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

Download in other formats:

Original Format