Context Navigation

ocean_forced_mod.F90 @ 3927

Last change on this file since 3927 was 3927, checked in by Laurent Fairhead, 3 years ago
Initial import of the physics wih isotopes from Camille Risi CR
File size: 18.5 KB

Line
1	!
2	! $Id: ocean_forced_mod.F90 3327 2018-05-16 13:58:59Z musat $
3	!
4	MODULE ocean_forced_mod
5	!
6	! This module is used for both the sub-surfaces ocean and sea-ice for the case of a
7	! forced ocean, "ocean=force".
8	!
9	IMPLICIT NONE
10
11	CONTAINS
12	!
13	!****************************************************************************************
14	!
15	SUBROUTINE ocean_forced_noice( &
16	itime, dtime, jour, knon, knindex, &
17	p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, &
18	temp_air, spechum, &
19	AcoefH, AcoefQ, BcoefH, BcoefQ, &
20	AcoefU, AcoefV, BcoefU, BcoefV, &
21	ps, u1, v1, gustiness, &
22	radsol, snow, agesno, &
23	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
24	tsurf_new, dflux_s, dflux_l &
25	#ifdef ISO
26	,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
27	xtsnow,xtevap,h1 &
28	#endif
29	)
30	!
31	! This subroutine treats the "open ocean", all grid points that are not entierly covered
32	! by ice.
33	! The routine receives data from climatologie file limit.nc and does some calculations at the
34	! surface.
35	!
36	USE dimphy
37	USE calcul_fluxs_mod
38	USE limit_read_mod
39	USE mod_grid_phy_lmdz
40	USE indice_sol_mod
41	USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
42	#ifdef ISO
43	USE infotrac_phy, ONLY: ntraciso,niso
44	USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, &
45	& calcul_iso_surf_sic_vectall
46	#ifdef ISOVERIF
47	USE isotopes_mod, ONLY: iso_eau,ridicule
48	!USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
49	USE isotopes_verif_mod
50	#endif
51	#endif
52
53	INCLUDE "YOMCST.h"
54	INCLUDE "clesphys.h"
55
56
57	! Input arguments
58	!****************************************************************************************
59	INTEGER, INTENT(IN) :: itime, jour, knon
60	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
61	REAL, INTENT(IN) :: dtime
62	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
63	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm
64	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
65	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
66	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
67	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
68	REAL, DIMENSION(klon), INTENT(IN) :: ps
69	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
70	#ifdef ISO
71	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
72	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtspechum
73	real, dimension(klon), intent(IN) :: rlat
74	#endif
75
76	! In/Output arguments
77	!****************************************************************************************
78	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
79	REAL, DIMENSION(klon), INTENT(INOUT) :: snow
80	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno !? put to 0 in ocean
81	#ifdef ISO
82	REAL, DIMENSION(niso,klon), INTENT(IN) :: xtsnow
83	REAL, DIMENSION(niso,klon), INTENT(INOUT) :: Roce
84	#endif
85
86	! Output arguments
87	!****************************************************************************************
88	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
89	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
90	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
91	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
92	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
93	#ifdef ISO
94	REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: xtevap ! isotopes in evaporation flux
95	REAL, DIMENSION(klon), INTENT(out) :: h1 ! just a diagnostic, not useful for the simulation
96	#endif
97
98	! Local variables
99	!****************************************************************************************
100	INTEGER :: i, j
101	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
102	REAL, DIMENSION(klon) :: alb_neig, tsurf_lim, zx_sl
103	REAL, DIMENSION(klon) :: u0, v0
104	REAL, DIMENSION(klon) :: u1_lay, v1_lay
105	LOGICAL :: check=.FALSE.
106	REAL, DIMENSION(klon) :: sens_prec_liq, sens_prec_sol
107	REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol
108
109	#ifdef ISO
110	real, parameter :: t_coup = 273.15
111	#endif
112
113	!****************************************************************************************
114	! Start calculation
115	!****************************************************************************************
116	IF (check) WRITE(,)' Entering ocean_forced_noice'
117
118	#ifdef ISO
119	#ifdef ISOVERIF
120	do i=1,knon
121	if (iso_eau.gt.0) then
122	call iso_verif_egalite_choix(xtspechum(iso_eau,i), &
123	& spechum(i),'ocean_forced_mod 111', &
124	& errmax,errmaxrel)
125	call iso_verif_egalite_choix(snow(i), &
126	& xtsnow(iso_eau,i),'ocean_forced_mod 117', &
127	& errmax,errmaxrel)
128	endif !if (iso_eau.gt.0) then
129	enddo !do i=1,knon
130	#endif
131	#endif
132
133	!****************************************************************************************
134	! 1)
135	! Read sea-surface temperature from file limit.nc
136	!
137	!****************************************************************************************
138	!--sb:
139	!!jyg if (knon.eq.1) then ! single-column model
140	if (klon_glo.eq.1) then ! single-column model
141	CALL read_tsurf1d(knon,tsurf_lim) ! new
142	#ifdef ISO
143	write(,) 'ocean_forced_mod 143: isotopes pas prévus ici'
144	stop
145	#endif
146	else ! GCM
147	CALL limit_read_sst(knon,knindex,tsurf_lim &
148	#ifdef ISO
149	& ,Roce,rlat &
150	#endif
151	& )
152	endif ! knon
153	!sb--
154
155	!****************************************************************************************
156	! 2)
157	! Flux calculation
158	!
159	!****************************************************************************************
160	! Set some variables for calcul_fluxs
161	cal = 0.
162	beta = 1.
163	dif_grnd = 0.
164	alb_neig(:) = 0.
165	agesno(:) = 0.
166	sens_prec_liq = 0.; sens_prec_sol = 0.; lat_prec_liq = 0.; lat_prec_sol = 0.
167
168	! Suppose zero surface speed
169	u0(:)=0.0
170	v0(:)=0.0
171	u1_lay(:) = u1(:) - u0(:)
172	v1_lay(:) = v1(:) - v0(:)
173
174	! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf
175	CALL calcul_fluxs(knon, is_oce, dtime, &
176	tsurf_lim, p1lay, cal, beta, cdragh, cdragq, ps, &
177	precip_rain, precip_snow, snow, qsurf, &
178	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
179	f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
180	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
181	sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
182
183	do j = 1, knon
184	i = knindex(j)
185	sens_prec_liq_o(i,1) = sens_prec_liq(j)
186	sens_prec_sol_o(i,1) = sens_prec_sol(j)
187	lat_prec_liq_o(i,1) = lat_prec_liq(j)
188	lat_prec_sol_o(i,1) = lat_prec_sol(j)
189	enddo
190
191
192	! - Flux calculation at first modele level for U and V
193	CALL calcul_flux_wind(knon, dtime, &
194	u0, v0, u1, v1, gustiness, cdragm, &
195	AcoefU, AcoefV, BcoefU, BcoefV, &
196	p1lay, temp_air, &
197	flux_u1, flux_v1)
198
199
200	#ifdef ISO
201	call calcul_iso_surf_oce_vectall(klon, knon,t_coup, &
202	& ps,tsurf_new,spechum,u1_lay, v1_lay, xtspechum, &
203	& evap, Roce,xtevap,h1 &
204	#ifdef ISOTRAC
205	& ,knindex &
206	#endif
207	& )
208	#endif
209
210	#ifdef ISO
211	#ifdef ISOVERIF
212	! write(,) 'ocean_forced_mod 176: sortie de ocean_forced_noice'
213	if (iso_eau.gt.0) then
214	do i=1,knon
215	call iso_verif_egalite_choix(snow(i), &
216	& xtsnow(iso_eau,i),'ocean_forced_mod 180', &
217	& errmax,errmaxrel)
218	enddo ! do j=1,knon
219	endif !if (iso_eau.gt.0) then
220	#endif
221	#endif
222
223	END SUBROUTINE ocean_forced_noice
224	!
225	!***************************************************************************************
226	!
227	SUBROUTINE ocean_forced_ice( &
228	itime, dtime, jour, knon, knindex, &
229	tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
230	AcoefH, AcoefQ, BcoefH, BcoefQ, &
231	AcoefU, AcoefV, BcoefU, BcoefV, &
232	ps, u1, v1, gustiness, &
233	radsol, snow, qsol, agesno, tsoil, &
234	qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
235	tsurf_new, dflux_s, dflux_l &
236	#ifdef ISO
237	,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &
238	xtsnow, xtsol,xtevap,Rland_ice &
239	#endif
240	)
241	!
242	! This subroutine treats the ocean where there is ice.
243	! The routine reads data from climatologie file and does flux calculations at the
244	! surface.
245	!
246	USE dimphy
247	USE calcul_fluxs_mod
248	USE surface_data, ONLY : calice, calsno
249	USE limit_read_mod
250	USE fonte_neige_mod, ONLY : fonte_neige
251	USE indice_sol_mod
252	USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
253	#ifdef ISO
254	USE infotrac_phy, ONLY: niso,ntraciso
255	USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, &
256	& calcul_iso_surf_sic_vectall
257	#ifdef ISOVERIF
258	USE isotopes_mod, ONLY: iso_eau,ridicule
259	!USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
260	USE isotopes_verif_mod
261	#endif
262	#endif
263
264	! INCLUDE "indicesol.h"
265	INCLUDE "dimsoil.h"
266	INCLUDE "YOMCST.h"
267	INCLUDE "clesphys.h"
268
269	! Input arguments
270	!****************************************************************************************
271	INTEGER, INTENT(IN) :: itime, jour, knon
272	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
273	REAL, INTENT(IN) :: dtime
274	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
275	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
276	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
277	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
278	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
279	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
280	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
281	REAL, DIMENSION(klon), INTENT(IN) :: ps
282	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
283	#ifdef ISO
284	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
285	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtspechum
286	REAL, DIMENSION(niso,klon), INTENT(IN) :: Roce
287	REAL, DIMENSION(niso,klon), INTENT(IN) :: Rland_ice
288	#endif
289
290	! In/Output arguments
291	!****************************************************************************************
292	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
293	REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
294	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
295	REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
296	#ifdef ISO
297	REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtsnow
298	REAL, DIMENSION(niso,klon), INTENT(IN) :: xtsol
299	#endif
300
301	! Output arguments
302	!****************************************************************************************
303	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
304	REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval
305	REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval
306	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
307	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
308	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
309	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
310	#ifdef ISO
311	REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: xtevap
312	#endif
313
314	! Local variables
315	!****************************************************************************************
316	LOGICAL :: check=.FALSE.
317	INTEGER :: i, j
318	REAL :: zfra
319	REAL, PARAMETER :: t_grnd=271.35
320	REAL, DIMENSION(klon) :: cal, beta, dif_grnd, capsol
321	REAL, DIMENSION(klon) :: alb_neig, tsurf_tmp
322	REAL, DIMENSION(klon) :: soilcap, soilflux
323	REAL, DIMENSION(klon) :: u0, v0
324	REAL, DIMENSION(klon) :: u1_lay, v1_lay
325	REAL, DIMENSION(klon) :: sens_prec_liq, sens_prec_sol
326	REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol
327
328	#ifdef ISO
329	real, parameter :: t_coup = 273.15
330	real, dimension(klon) :: fq_fonte_diag
331	real, dimension(klon) :: fqfonte_diag
332	real, dimension(klon) :: snow_evap_diag
333	real, dimension(klon) :: fqcalving_diag
334	real, DIMENSION(klon) :: run_off_lic_diag
335	real :: coeff_rel_diag
336	real :: max_eau_sol_diag
337	real, dimension(klon) :: runoff_diag
338	integer ixt
339	REAL, DIMENSION(niso,klon) :: xtsnow_prec,xtsol_prec
340	REAL, DIMENSION(klon) :: snow_prec,qsol_prec
341	#endif
342
343	!****************************************************************************************
344	! Start calculation
345	!****************************************************************************************
346	IF (check) WRITE(,)'Entering surface_seaice, knon=',knon
347
348	!****************************************************************************************
349	! 1)
350	! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, flux_u1, flux_v1
351	! dflux_s, dflux_l and qsurf
352	!****************************************************************************************
353
354	tsurf_tmp(:) = tsurf_in(:)
355
356	! calculate the parameters cal, beta, capsol and dif_grnd
357	CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
358
359
360	IF (soil_model) THEN
361	! update tsoil and calculate soilcap and soilflux
362	CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, tsoil,soilcap, soilflux)
363	cal(1:knon) = RCPD / soilcap(1:knon)
364	radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
365	dif_grnd = 1.0 / tau_gl
366	ELSE
367	dif_grnd = 1.0 / tau_gl
368	cal = RCPD * calice
369	WHERE (snow > 0.0) cal = RCPD * calsno
370	ENDIF
371
372	beta = 1.0
373	sens_prec_liq = 0.; sens_prec_sol = 0.; lat_prec_liq = 0.; lat_prec_sol = 0.
374
375	! Suppose zero surface speed
376	u0(:)=0.0
377	v0(:)=0.0
378	u1_lay(:) = u1(:) - u0(:)
379	v1_lay(:) = v1(:) - v0(:)
380	CALL calcul_fluxs(knon, is_sic, dtime, &
381	tsurf_tmp, p1lay, cal, beta, cdragh, cdragh, ps, &
382	precip_rain, precip_snow, snow, qsurf, &
383	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
384	f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
385	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
386	sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
387	do j = 1, knon
388	i = knindex(j)
389	sens_prec_liq_o(i,2) = sens_prec_liq(j)
390	sens_prec_sol_o(i,2) = sens_prec_sol(j)
391	lat_prec_liq_o(i,2) = lat_prec_liq(j)
392	lat_prec_sol_o(i,2) = lat_prec_sol(j)
393	enddo
394
395	! - Flux calculation at first modele level for U and V
396	CALL calcul_flux_wind(knon, dtime, &
397	u0, v0, u1, v1, gustiness, cdragm, &
398	AcoefU, AcoefV, BcoefU, BcoefV, &
399	p1lay, temp_air, &
400	flux_u1, flux_v1)
401
402	!****************************************************************************************
403	! 2)
404	! Calculations due to snow and runoff
405	!
406	!****************************************************************************************
407
408	#ifdef ISO
409	! verif
410	#ifdef ISOVERIF
411	do i=1,knon
412	if (iso_eau.gt.0) then
413	if (snow(i).gt.ridicule) then
414	call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
415	& 'interfsurf 964',errmax,errmaxrel)
416	endif !if ((snow(i).gt.ridicule)) then
417	endif !if (iso_eau.gt.0) then
418	enddo !do i=1,knon
419	#endif
420	! end verif
421
422	do i=1,knon
423	snow_prec(i)=snow(i)
424	do ixt=1,niso
425	xtsnow_prec(ixt,i)=xtsnow(ixt,i)
426	enddo !do ixt=1,niso
427	! initialisation:
428	fq_fonte_diag(i)=0.0
429	fqfonte_diag(i)=0.0
430	snow_evap_diag(i)=0.0
431	enddo !do i=1,knon
432	#endif
433
434	CALL fonte_neige( knon, is_sic, knindex, dtime, &
435	tsurf_tmp, precip_rain, precip_snow, &
436	snow, qsol, tsurf_new, evap &
437	#ifdef ISO
438	& ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag &
439	& ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag &
440	#endif
441	& )
442
443
444	#ifdef ISO
445	! isotopes: tout est externalisé
446	#ifdef ISOVERIF
447	write(,) 'ocean_forced_mod 377: call calcul_iso_surf_sic_vectall'
448	write(,) 'klon,knon=',klon,knon
449	#endif
450	call calcul_iso_surf_sic_vectall(klon,knon, &
451	& evap,snow_evap_diag,Tsurf_new,Roce,snow, &
452	& fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
453	& precip_snow,xtprecip_snow,xtprecip_rain, snow_prec,xtsnow_prec, &
454	& xtspechum,spechum,ps, &
455	& xtevap,xtsnow,fqcalving_diag, &
456	& knindex,is_sic,run_off_lic_diag,coeff_rel_diag,Rland_ice &
457	& )
458	#ifdef ISOVERIF
459	write(,) 'ocean_forced_mod 391: sortie calcul_iso_surf_sic_vectall'
460	if (iso_eau.gt.0) then
461	do i=1,knon
462	call iso_verif_egalite_choix(snow(i), &
463	& xtsnow(iso_eau,i),'ocean_forced_mod 396', &
464	& errmax,errmaxrel)
465	enddo ! do j=1,knon
466	endif !if (iso_eau.gt.0) then
467	#endif
468	!#ifdef ISOVERIF
469	#endif
470	!#ifdef ISO
471
472	! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
473	!
474	CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:))
475
476	WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.
477
478	alb1_new(:) = 0.0
479	DO i=1, knon
480	zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
481	alb1_new(i) = alb_neig(i) * zfra + 0.6 * (1.0-zfra)
482	ENDDO
483
484	alb2_new(:) = alb1_new(:)
485
486	END SUBROUTINE ocean_forced_ice
487
488	!************************************************************************
489	! 1D case
490	!************************************************************************
491	SUBROUTINE read_tsurf1d(knon,sst_out)
492
493	! This subroutine specifies the surface temperature to be used in 1D simulations
494
495	USE dimphy, ONLY : klon
496
497	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
498	REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model
499
500	INTEGER :: i
501	! COMMON defined in lmdz1d.F:
502	real ts_cur
503	common /sst_forcing/ts_cur
504
505	DO i = 1, knon
506	sst_out(i) = ts_cur
507	ENDDO
508
509	END SUBROUTINE read_tsurf1d
510
511	!
512	!************************************************************************
513	!
514	END MODULE ocean_forced_mod
515
516
517
518
519
520

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

Download in other formats:

Original Format