Context Navigation

ocean_forced_mod.F90 @ 5802

Last change on this file since 5802 was 5662, checked in by Laurent Fairhead, 6 months ago
Ajout du modèle thermodynamique de glace de mer interactive améliorant les flux échangés à la surface de la banquise (Doctorat de Nicolas Michalezyk, Contact : Nicolas Michaleyk, Guillaume Gastineau)
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: 39.7 KB

Line
1	!
2	! $Id: ocean_forced_mod.F90 5662 2025-05-20 14:24:41Z snguyen $
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, tsurf_in, &
22	radsol, snow, agesno, &
23	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
24	tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa, &
25	dthetadz300,pctsrf,Ampl &
26	#ifdef ISO
27	,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
28	xtsnow,xtevap,h1 &
29	#endif
30	)
31	!
32	! This subroutine treats the "open ocean", all grid points that are not entierly covered
33	! by ice.
34	! The routine receives data from climatologie file limit.nc and does some calculations at the
35	! surface.
36	!
37	USE dimphy
38	USE calcul_fluxs_mod
39	USE limit_read_mod
40	USE mod_grid_phy_lmdz
41	USE indice_sol_mod
42	USE surface_data, ONLY : iflag_leads
43	USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
44	use config_ocean_skin_m, only: activate_ocean_skin
45	#ifdef ISO
46	USE infotrac_phy, ONLY: ntiso,niso
47	USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall
48	#ifdef ISOVERIF
49	USE isotopes_mod, ONLY: iso_eau,ridicule
50	!USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
51	USE isotopes_verif_mod
52	#endif
53	#endif
54	USE flux_arp_mod_h
55	USE clesphys_mod_h
56	USE yomcst_mod_h
57
58	! Input arguments
59	!****************************************************************************************
60	INTEGER, INTENT(IN) :: itime, jour, knon
61	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
62	REAL, INTENT(IN) :: dtime
63	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
64	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm
65	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
66	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
67	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
68	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
69	REAL, DIMENSION(klon), INTENT(IN) :: ps
70	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
71	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
72	real, intent(in):: rhoa(:) ! (knon) density of moist air (kg / m3)
73	!GG
74	REAL, DIMENSION(klon), INTENT(IN) :: dthetadz300
75	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
76	!
77
78	#ifdef ISO
79	REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
80	REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum
81	REAL, DIMENSION(klon), INTENT(IN) :: rlat
82	#endif
83
84	! In/Output arguments
85	!****************************************************************************************
86	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
87	REAL, DIMENSION(klon), INTENT(INOUT) :: snow
88	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno !? put to 0 in ocean
89	#ifdef ISO
90	REAL, DIMENSION(niso,klon), INTENT(IN) :: xtsnow
91	REAL, DIMENSION(niso,klon), INTENT(INOUT):: Roce
92	#endif
93
94	! Output arguments
95	!****************************************************************************************
96	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
97	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
98	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
99	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
100	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
101	REAL, INTENT(out):: sens_prec_liq(:) ! (knon)
102	!GG
103	REAL, DIMENSION(klon), INTENT(OUT) :: Ampl
104	!
105
106	#ifdef ISO
107	REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap ! isotopes in evaporation flux
108	REAL, DIMENSION(klon), INTENT(OUT) :: h1 ! just a diagnostic, not useful for the simulation
109	#endif
110
111	! Local variables
112	!****************************************************************************************
113	INTEGER :: i, j
114	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
115	REAL, DIMENSION(klon) :: alb_neig, tsurf_lim, zx_sl
116	REAL, DIMENSION(klon) :: u0, v0
117	REAL, DIMENSION(klon) :: u1_lay, v1_lay
118	LOGICAL :: check=.FALSE.
119	REAL, DIMENSION(knon) :: sens_prec_sol
120	REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol
121	! GG
122	REAL, DIMENSION(klon) :: l_CBL, sicfra
123	!
124	#ifdef ISO
125	REAL, PARAMETER :: t_coup = 273.15
126	#endif
127
128
129	!****************************************************************************************
130	! Start calculation
131	!****************************************************************************************
132	IF (check) WRITE(,)' Entering ocean_forced_noice'
133
134	#ifdef ISO
135	#ifdef ISOVERIF
136	DO i = 1, knon
137	IF (iso_eau > 0) THEN
138	CALL iso_verif_egalite_choix(xtspechum(iso_eau,i), &
139	& spechum(i),'ocean_forced_mod 111', &
140	& errmax,errmaxrel)
141	CALL iso_verif_egalite_choix(snow(i), &
142	& xtsnow(iso_eau,i),'ocean_forced_mod 117', &
143	& errmax,errmaxrel)
144	ENDIF !IF (iso_eau > 0) THEN
145	ENDDO !DO i=1,knon
146	#endif
147	#endif
148
149	!****************************************************************************************
150	! 1)
151	! Read sea-surface temperature from file limit.nc
152	!
153	!****************************************************************************************
154	!--sb:
155	!!jyg if (knon.eq.1) then ! single-column model
156	if (klon_glo.eq.1) then ! single-column model
157	! EV: now surface Tin flux_arp.h
158	!CALL read_tsurf1d(knon,tsurf_lim) ! new
159	DO i = 1, knon
160	tsurf_lim(i) = tg
161	ENDDO
162
163	else ! GCM
164	CALL limit_read_sst(knon,knindex,tsurf_lim &
165	#ifdef ISO
166	& ,Roce,rlat &
167	#endif
168	& )
169	endif ! knon
170	!sb--
171
172	!****************************************************************************************
173	! 2)
174	! Flux calculation
175	!
176	!****************************************************************************************
177	! Set some variables for calcul_fluxs
178	!cal = 0.
179	!beta = 1.
180	!dif_grnd = 0.
181
182
183	! EV: use calbeta to calculate beta
184	! Need to initialize qsurf for calbeta but it is not modified by this routine
185	qsurf(:)=0.
186	CALL calbeta(dtime, is_oce, knon, snow, qsurf, beta, cal, dif_grnd)
187
188
189	alb_neig(:) = 0.
190	agesno(:) = 0.
191	lat_prec_liq = 0.; lat_prec_sol = 0.
192
193	! Suppose zero surface speed
194	u0(:)=0.0
195	v0(:)=0.0
196	u1_lay(:) = u1(:) - u0(:)
197	v1_lay(:) = v1(:) - v0(:)
198
199	! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf
200	CALL calcul_fluxs(knon, is_oce, dtime, &
201	merge(tsurf_in, tsurf_lim, activate_ocean_skin == 2), p1lay, cal, &
202	beta, cdragh, cdragq, ps, &
203	precip_rain, precip_snow, snow, qsurf, &
204	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
205	f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
206	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
207	sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
208	if (activate_ocean_skin == 2) tsurf_new = tsurf_lim
209
210	do j = 1, knon
211	i = knindex(j)
212	sens_prec_liq_o(i,1) = sens_prec_liq(j)
213	sens_prec_sol_o(i,1) = sens_prec_sol(j)
214	lat_prec_liq_o(i,1) = lat_prec_liq(j)
215	lat_prec_sol_o(i,1) = lat_prec_sol(j)
216	enddo
217
218	!GG
219	if (iflag_leads == 1) then
220	l_CBL = -52381.*dthetadz300 + 3008.1
221	Ampl = 6.012e-08l_CBL2 - 4.036e-04l_CBL + 1.4979
222	WHERE(Ampl(:)>1.2) Ampl(:)=1.2
223	sicfra(:)=pctsrf(:,is_sic)/(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
224	WHERE(pctsrf(:,is_sic)+pctsrf(:,is_oce)<EPSFRA) sicfra(:)=0.
225	WHERE(sicfra<0.7) Ampl(:)=1.
226	WHERE((sicfra>0.7).and.(sicfra<0.9)) Ampl=((sicfra-0.7)/0.2)*Ampl+((0.9-sicfra)/0.2)
227	fluxsens=Ampl*fluxsens
228	dflux_s=Ampl*dflux_s
229	endif
230
231
232	! - Flux calculation at first modele level for U and V
233	CALL calcul_flux_wind(knon, dtime, &
234	u0, v0, u1, v1, gustiness, cdragm, &
235	AcoefU, AcoefV, BcoefU, BcoefV, &
236	p1lay, temp_air, &
237	flux_u1, flux_v1)
238
239	#ifdef ISO
240	CALL calcul_iso_surf_oce_vectall(klon, knon,t_coup, &
241	& ps,tsurf_new,spechum,u1_lay, v1_lay, xtspechum, &
242	& evap, Roce,xtevap,h1 &
243	#ifdef ISOTRAC
244	& ,knindex &
245	#endif
246	& )
247	#endif
248
249	#ifdef ISO
250	#ifdef ISOVERIF
251	! write(,) 'ocean_forced_mod 176: sortie de ocean_forced_noice'
252	IF (iso_eau > 0) THEN
253	DO i = 1, knon
254	CALL iso_verif_egalite_choix(snow(i), &
255	& xtsnow(iso_eau,i),'ocean_forced_mod 180', &
256	& errmax,errmaxrel)
257	ENDDO ! DO j=1,knon
258	ENDIF !IF (iso_eau > 0) THEN
259	#endif
260	#endif
261
262	END SUBROUTINE ocean_forced_noice
263	!
264	!***************************************************************************************
265	!
266	SUBROUTINE ocean_forced_ice( &
267	itime, dtime, jour, knon, knindex, &
268	tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air,spechum, &
269	AcoefH, AcoefQ, BcoefH, BcoefQ, &
270	AcoefU, AcoefV, BcoefU, BcoefV, &
271	!GG ps, u1, v1, gustiness, &
272	ps, u1, v1, gustiness, pctsrf, &
273	!GG
274	radsol, snow, qsol, agesno, tsoil, &
275	qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
276	!GG tsurf_new, dflux_s, dflux_l, rhoa)
277	tsurf_new, dflux_s, dflux_l, rhoa, swnet, hice, tice, bilg_cumul, &
278	fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, &
279	dtice_melt, dtice_snow2sic &
280	!GG
281	#ifdef ISO
282	,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &
283	xtsnow, xtsol,xtevap,Rland_ice &
284	#endif
285	)
286	!
287	! This subroutine treats the ocean where there is ice.
288	! The routine reads data from climatologie file and does flux calculations at the
289	! surface.
290	!
291	USE dimphy
292	USE geometry_mod, ONLY: longitude,latitude
293	USE calcul_fluxs_mod
294	!GG USE surface_data, ONLY : calice, calsno
295	USE surface_data, ONLY : calice, calsno, iflag_seaice, iflag_seaice_alb, &
296	sice_cond, sisno_cond, sisno_den, sisno_min, sithick_min, sisno_wfact, &
297	amax_s,amax_n, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, &
298	si_pen_frac, si_pen_ext, fseaN, fseaS, iflag_leads
299
300	USE geometry_mod, ONLY: longitude,latitude,latitude_deg
301	!GG
302	USE limit_read_mod
303	USE fonte_neige_mod, ONLY : fonte_neige
304	USE indice_sol_mod
305	USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
306	#ifdef ISO
307	USE infotrac_phy, ONLY: niso, ntiso
308	USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall
309	#ifdef ISOVERIF
310	USE isotopes_mod, ONLY: iso_eau,ridicule
311	!USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
312	USE isotopes_verif_mod
313	#endif
314	#endif
315	USE flux_arp_mod_h
316	USE clesphys_mod_h
317	USE yomcst_mod_h
318	USE dimsoil_mod_h, ONLY: nsoilmx
319
320	! INCLUDE "indicesol.h"
321
322
323	! Input arguments
324	!****************************************************************************************
325	INTEGER, INTENT(IN) :: itime, jour, knon
326	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
327	REAL, INTENT(IN) :: dtime
328	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
329	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
330	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
331	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
332	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
333	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
334	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
335	REAL, DIMENSION(klon), INTENT(IN) :: ps
336	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
337	real, intent(in):: rhoa(:) ! (knon) density of moist air (kg / m3)
338	!GG
339	REAL, DIMENSION(klon), INTENT(IN) :: swnet
340	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
341	!GG
342	#ifdef ISO
343	REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
344	REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum
345	REAL, DIMENSION(niso,klon), INTENT(IN) :: Roce
346	REAL, DIMENSION(niso,klon), INTENT(IN) :: Rland_ice
347	#endif
348
349	! In/Output arguments
350	!****************************************************************************************
351	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
352	REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
353	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
354	REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
355	!GG
356	REAL, DIMENSION(klon), INTENT(INOUT) :: hice
357	REAL, DIMENSION(klon), INTENT(INOUT) :: tice
358	REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul
359	REAL, DIMENSION(klon), INTENT(INOUT) :: fcds
360	REAL, DIMENSION(klon), INTENT(INOUT) :: fcdi
361	REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_growth
362	REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_melt
363	REAL, DIMENSION(klon), INTENT(INOUT) :: dh_top_melt
364	REAL, DIMENSION(klon), INTENT(INOUT) :: dh_snow2sic
365	REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_melt
366	REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_snow2sic
367	!GG
368	#ifdef ISO
369	REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtsnow
370	REAL, DIMENSION(niso,klon), INTENT(IN) :: xtsol
371	#endif
372
373	! Output arguments
374	!****************************************************************************************
375	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
376	REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval
377	REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval
378	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
379	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
380	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
381	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
382	#ifdef ISO
383	REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap
384	#endif
385
386	! Local variables
387	!****************************************************************************************
388	LOGICAL :: check=.FALSE.
389	INTEGER :: i, j
390	REAL :: zfra
391	REAL, PARAMETER :: t_grnd=271.35
392	REAL, DIMENSION(klon) :: cal, beta, dif_grnd, capsol, icesub
393	REAL, DIMENSION(klon) :: alb_neig, tsurf_tmp
394	REAL, DIMENSION(klon) :: soilcap, soilflux
395	REAL, DIMENSION(klon) :: u0, v0
396	REAL, DIMENSION(klon) :: u1_lay, v1_lay
397	REAL, DIMENSION(knon) :: sens_prec_liq, sens_prec_sol
398	REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol
399	!GG
400	INTEGER :: ki
401	INTEGER :: cpl_pas
402	REAL, DIMENSION(klon) :: bilg, fsic, f_bot
403	REAL, PARAMETER :: latent_ice = 334.0e3
404	REAL, PARAMETER :: rau_ice = 917.0
405	REAL, PARAMETER :: kice=2.2
406	REAL :: f_cond, f_swpen, f_cond_s, f_cond_i
407	REAL :: ustar, uscap, ustau
408	! for snow/ice albedo:
409	REAL :: alb_snow, alb_ice, alb_pond
410	REAL :: frac_snow, frac_ice, frac_pond
411	REAL :: z1_i, z2_i, z1_s, zlog ! height parameters
412	! for ice melt / freeze
413	REAL :: e_melt, snow_evap, h_test
414	! dhsic, dfsic change in ice mass, fraction.
415	REAL :: dhsic, dfsic, frac_mf
416	REAL :: fsea, amax
417	REAL :: hice_i, tice_i, fsic_new
418	! snow and ice physical characteristics:
419	REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
420	REAL, PARAMETER :: t_melt=273.15 ! melting ice temp
421	REAL :: sno_den!=sisno_den !mean snow density, kg/m3
422	REAL, PARAMETER :: ice_den=917. ! ice density
423	REAL, PARAMETER :: sea_den=1025. ! sea water density
424	REAL :: ice_cond!=sice_cond*ice_den !conductivity of ice
425	REAL :: sno_cond!=sisno_cond*sno_den ! conductivity of snow
426	REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice
427	REAL, PARAMETER :: sea_cap=3995. ! specific heat capacity, water
428	REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
429
430	! control of snow and ice cover & freeze / melt (heights converted to kg/m2)
431	REAL :: snow_min!=sisno_min*sno_den !critical snow height 5 cm
432	REAL :: snow_wfact!=sisno_wfact ! max fraction of falling snow blown into ocean
433	REAL, PARAMETER :: ice_frac_min=0.005
434	REAL :: h_ice_min!=sithick_min ! min ice thickness
435	! below ice_thin, priority is melt lateral / grow height
436	! ice_thin is also height of new ice
437	REAL, PARAMETER :: h_ice_max=7 ! max ice height
438	! Ice thickness parameter for lateral growth
439	REAL, PARAMETER :: h_ice_thick=1.5
440	REAL, PARAMETER :: h_ice_thin=0.15
441
442	! albedo and radiation parameters
443	INTEGER, SAVE :: iflag_sic_albedo
444	! albedo old or NEMO
445	REAL :: alb_sno_dry!=rn_alb_sdry !dry snow albedo
446	REAL :: alb_sno_wet!=rn_alb_smlt !wet snow albedo
447	REAL :: alb_ice_dry!=rn_alb_idry !dry thick ice
448	REAL :: alb_ice_wet!=rn_alb_imlt !melting thick ice
449	! new (Toyoda 2020) albedo
450	! Values for snow / ice, dry / melting, visible / near IR
451	REAL, PARAMETER :: alb_sdry_vis=0.98
452	REAL, PARAMETER :: alb_smlt_vis=0.88
453	REAL, PARAMETER :: alb_sdry_nir=0.7
454	REAL, PARAMETER :: alb_smlt_nir=0.55
455	REAL, PARAMETER :: alb_idry_vis=0.78
456	REAL, PARAMETER :: alb_imlt_vis=0.705
457	REAL, PARAMETER :: alb_idry_nir=0.36
458	REAL, PARAMETER :: alb_imlt_nir=0.285
459	REAL, PARAMETER :: h_ice_alb=0.5*ice_den ! height for full ice albedo
460	REAL, PARAMETER :: h_sno_alb=0.02*300 ! height for control of snow fraction
461
462	REAL :: pen_frac !=si_pen_frac !fraction of shortwave penetrating into the
463	! ice (not snow). Should be visible only, not NIR
464	REAL :: pen_ext !=si_pen_ext !extinction length of penetrating shortwave (m-1)
465
466	! HF from ocean below ice
467	! REAL, PARAMETER :: fseaN=2.0 ! NH
468	! REAL, PARAMETER :: fseaS=4.0 ! SH
469	!GG
470
471	#ifdef ISO
472	REAL, PARAMETER :: t_coup = 273.15
473	REAL, DIMENSION(klon) :: fq_fonte_diag
474	REAL, DIMENSION(klon) :: fqfonte_diag
475	REAL, DIMENSION(klon) :: snow_evap_diag
476	REAL, DIMENSION(klon) :: fqcalving_diag
477	REAL, DIMENSION(klon) :: run_off_lic_diag
478	REAL :: coeff_rel_diag
479	REAL :: max_eau_sol_diag
480	REAL, DIMENSION(klon) :: runoff_diag
481	INTEGER IXT
482	REAL, DIMENSION(niso,klon) :: xtsnow_prec, xtsol_prec
483	REAL, DIMENSION(klon) :: snow_prec, qsol_prec
484	#endif
485
486	!****************************************************************************************
487	! Start calculation
488	!****************************************************************************************
489	IF (check) WRITE(,)'Entering surface_seaice, knon=',knon
490
491	!****************************************************************************************
492	! 1)
493	! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, flux_u1, flux_v1
494	! dflux_s, dflux_l and qsurf
495	!****************************************************************************************
496
497	tsurf_tmp(:) = tsurf_in(:)
498
499	!GG
500	IF (iflag_seaice==0) THEN ! Old LMDZ sea ice surface
501	!GG
502	! calculate the parameters cal, beta, capsol and dif_grnd and then recalculate cal
503	CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
504
505
506	IF (soil_model) THEN
507	! update tsoil and calculate soilcap and soilflux
508	CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, qsol, &
509	& longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil,soilcap, soilflux)
510	cal(1:knon) = RCPD / soilcap(1:knon)
511	radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
512	dif_grnd = 1.0 / tau_gl
513	ELSE
514	dif_grnd = 1.0 / tau_gl
515	cal = RCPD * calice
516	WHERE (snow > 0.0) cal = RCPD * calsno
517	ENDIF
518
519	!GG
520	ELSEIF (iflag_seaice==2) THEN
521
522	sno_den=sisno_den !mean snow density, kg/m3
523	ice_cond=sice_cond*ice_den !conductivity of ice
524	sno_cond=sisno_cond*sno_den ! conductivity of snow
525	snow_min=sisno_min*sno_den !critical snow height 5 cm
526	snow_wfact=sisno_wfact ! max fraction of falling snow blown into ocean
527	h_ice_min=sithick_min ! min ice thickness
528	alb_sno_dry=rn_alb_sdry !dry snow albedo
529	alb_sno_wet=rn_alb_smlt !wet snow albedo
530	alb_ice_dry=rn_alb_idry !dry thick ice
531	alb_ice_wet=rn_alb_imlt !melting thick ice
532	pen_frac=si_pen_frac !fraction of shortwave penetrating into the
533	pen_ext=si_pen_ext !extinction length of penetrating shortwave (m-1)
534
535	bilg(:)=0.
536	dif_grnd(:)=0.
537	beta(:) = 1.
538	fsic(:) = pctsrf(:,is_sic)
539	cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
540
541	! Surface, snow-ice and ice-ocean fluxes.
542	! Prepare call to calcul_fluxs (cal, beta, radsol, dif_grnd)
543
544	!write(,) 'radsol 1',radsol(1:100)
545	DO i=1,knon
546	ki=knindex(i)
547	IF (snow(i).GT.snow_min) THEN
548	! 1 / snow-layer heat capacity
549	cal(i)=2.RCPD/(snow(i)ice_cap)
550	! adjustment time-scale of conductive flux
551	dif_grnd(i) = cal(i) * sno_cond / snow(i) / RCPD
552	! for conductive flux
553	f_cond_s = sno_cond * (tice(ki)-t_freeze) / snow(i)
554	radsol(i) = radsol(i)+f_cond_s
555	! all shortwave flux absorbed
556	f_swpen=0.
557	ELSE ! bare ice.
558	f_cond_s = 0.
559	! 1 / ice-layer heat capacity
560	cal(i) = 2.RCPD/(hice(ki)ice_den*ice_cap)
561	! adjustment time-scale of conductive flux
562	dif_grnd(i) = cal(i) * ice_cond / (hice(ki)*ice_den) / RCPD
563	! penetrative shortwave flux...
564	f_swpen=swnet(i)pen_fracexp(-pen_ext*hice(ki))
565	radsol(i) = radsol(i)-f_swpen
566	! GG no conductive flux in this case?
567	END IF
568	bilg(ki)=f_swpen-f_cond_s
569	END DO
570
571	endif
572
573	! beta = 1.0
574	lat_prec_liq = 0.; lat_prec_sol = 0.
575
576	! Suppose zero surface speed
577	u0(:)=0.0
578	v0(:)=0.0
579	u1_lay(:) = u1(:) - u0(:)
580	v1_lay(:) = v1(:) - v0(:)
581	CALL calcul_fluxs(knon, is_sic, dtime, &
582	tsurf_tmp, p1lay, cal, beta, cdragh, cdragh, ps, &
583	precip_rain, precip_snow, snow, qsurf, &
584	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
585	f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
586	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
587	sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
588	do j = 1, knon
589	i = knindex(j)
590	sens_prec_liq_o(i,2) = sens_prec_liq(j)
591	sens_prec_sol_o(i,2) = sens_prec_sol(j)
592	lat_prec_liq_o(i,2) = lat_prec_liq(j)
593	lat_prec_sol_o(i,2) = lat_prec_sol(j)
594	enddo
595
596	! - Flux calculation at first modele level for U and V
597	CALL calcul_flux_wind(knon, dtime, &
598	u0, v0, u1, v1, gustiness, cdragm, &
599	AcoefU, AcoefV, BcoefU, BcoefV, &
600	p1lay, temp_air, &
601	flux_u1, flux_v1)
602
603	!****************************************************************************************
604	! 2)
605	! Calculations due to snow and runoff
606	!
607	!****************************************************************************************
608	!GG
609	if (iflag_seaice==0) then
610	!GG
611	#ifdef ISO
612	! verif
613	#ifdef ISOVERIF
614	DO i = 1, knon
615	IF (iso_eau > 0) THEN
616	IF (snow(i) > ridicule) THEN
617	CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
618	& 'interfsurf 964',errmax,errmaxrel)
619	ENDIF !IF ((snow(i) > ridicule)) THEN
620	ENDIF !IF (iso_eau > 0) THEN
621	ENDDO !DO i=1,knon
622	#endif
623	! end verif
624
625	DO i = 1, knon
626	snow_prec(i) = snow(i)
627	DO ixt = 1, niso
628	xtsnow_prec(ixt,i) = xtsnow(ixt,i)
629	ENDDO !DO ixt=1,niso
630	! initialisation:
631	fq_fonte_diag(i) = 0.0
632	fqfonte_diag(i) = 0.0
633	snow_evap_diag(i)= 0.0
634	ENDDO !DO i=1,knon
635	#endif
636
637
638	CALL fonte_neige( knon, is_sic, knindex, dtime, &
639	tsurf_tmp, precip_rain, precip_snow, &
640	snow, qsol, tsurf_new, evap, icesub &
641	#ifdef ISO
642	& ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag &
643	& ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag &
644	#endif
645	& )
646
647
648	#ifdef ISO
649	! isotopes: tout est externalisé
650	!#ifdef ISOVERIF
651	! write(,) 'ocean_forced_mod 377: call calcul_iso_surf_sic_vectall'
652	! write(,) 'klon,knon=',klon,knon
653	!#endif
654	CALL calcul_iso_surf_sic_vectall(klon,knon, &
655	& evap,snow_evap_diag,Tsurf_new,Roce,snow, &
656	& fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
657	& precip_snow,xtprecip_snow,xtprecip_rain, snow_prec,xtsnow_prec, &
658	& xtspechum,spechum,ps, &
659	& xtevap,xtsnow,fqcalving_diag, &
660	& knindex,is_sic,run_off_lic_diag,coeff_rel_diag,Rland_ice &
661	& )
662	#ifdef ISOVERIF
663	!write(,) 'ocean_forced_mod 391: sortie calcul_iso_surf_sic_vectall'
664	IF (iso_eau > 0) THEN
665	DO i = 1, knon
666	CALL iso_verif_egalite_choix(snow(i), &
667	& xtsnow(iso_eau,i),'ocean_forced_mod 396', &
668	& errmax,errmaxrel)
669	ENDDO ! DO j=1,knon
670	ENDIF !IF (iso_eau > 0) then
671	#endif
672	!#ifdef ISOVERIF
673	#endif
674	!#ifdef ISO
675
676	! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
677	!
678	CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:))
679
680	WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.
681
682	alb1_new(:) = 0.0
683	DO i=1, knon
684	zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
685	alb1_new(i) = alb_neig(i) * zfra + 0.6 * (1.0-zfra)
686	ENDDO
687
688	alb2_new(:) = alb1_new(:)
689
690	!GG
691	else
692
693	DO i=1,knon
694	ki=knindex(i)
695
696	! ocean-ice heat flux
697	fsea=fseaS
698	amax=amax_s
699	if (latitude(ki)>0) THEN
700	fsea=fseaN
701	amax=amax_n
702	ENDIF
703
704	IF (snow(i).GT.snow_min) THEN
705	! snow conductive flux after calcul_fluxs (pos up)
706	f_cond_s = sno_cond * (tice(ki)-tsurf_new(i)) / snow(i)
707	! 1 / heat capacity and conductive timescale
708	uscap = 2. / ice_cap / (snow(i)+hice(ki)*ice_den)
709	ustau = uscap * ice_cond / (hice(ki)*ice_den)
710	! update ice temp
711	tice(ki) = (tice(ki) + dtime(ustaut_freeze - uscap*f_cond_s)) &
712	/ (1 + dtime*ustau)
713	ELSE ! bare ice
714	tice(ki)=tsurf_new(i)
715	ENDIF
716	! ice conductive flux (pos up)
717	f_cond_i = ice_cond * (t_freeze-tice(ki)) / (hice(ki)*ice_den)
718	f_bot(i) = fsea - f_cond_i
719	fcdi(ki) = f_cond_i - fsea
720	fcds(i) = f_cond_s
721	!bilg(ki) = bilg(ki)+f_cond_i
722	END DO
723
724	!****************************************************************************************
725	! 2) Update snow and ice surface : thickness
726	!****************************************************************************************
727
728	IF (iflag_seaice==1) THEN
729	! Read from limit
730	CALL limit_read_hice(knon,knindex,hice)
731	ENDIF
732	! Formula Krinner et al. 1997 : h = (0.2 + 3.8(f_min**2))(1 + 2(f- f_min))
733
734
735
736	DO i=1,knon
737	ki=knindex(i)
738	IF (precip_snow(i) > 0.) THEN
739	snow(i) = snow(i)+precip_snow(i)dtime(1.-snow_wfact*(1.-fsic(ki)))
740	END IF
741	! snow and ice sublimation
742	IF (evap(i) > 0.) THEN
743	snow_evap = MIN (snow(i) / dtime, evap(i))
744	snow(i) = snow(i) - snow_evap * dtime
745	snow(i) = MAX(0.0, snow(i))
746	IF (iflag_seaice==2) THEN
747	hice(ki) = MAX(0.0,hice(ki)-(evap(i)-snow_evap)*dtime/ice_den)
748	ENDIF
749	ENDIF
750	! Melt / Freeze snow from above if Tsurf>0
751	IF (tsurf_new(i).GT.t_melt) THEN
752	! energy available for melting snow (in kg of melted snow /m2)
753	e_melt = MIN(MAX(snow(i)(tsurf_new(i)-t_melt)ice_cap/2. &
754	/(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
755	! remove snow
756	tice_i=tice(ki)
757	IF (snow(i).GT.e_melt) THEN
758	snow(i)=snow(i)-e_melt
759	tsurf_new(i)=t_melt
760	ELSE ! all snow is melted
761	! add remaining heat flux to ice
762	e_melt=e_melt-snow(i)
763	tice(ki)=tice(ki)+e_meltice_lat2./(ice_caphice(ki)ice_den)
764	tsurf_new(i)=tice(ki)
765	END IF
766	dtice_melt(ki)=(tice(ki)-tice_i)/dtime
767	END IF
768	! Bottom melt / grow
769	! bottom freeze if bottom flux (cond + oce-ice) <0
770	IF (iflag_seaice==2) THEN
771	IF (f_bot(i).LT.0) THEN
772	! larger fraction of bottom growth
773	frac_mf=MIN(1.,MAX(0.,(hice(ki)-h_ice_thick) &
774	/ (h_ice_max-h_ice_thick)))
775	! quantity of new ice (formed at mean ice temp)
776	e_melt= -f_bot(i) * dtime * fsic(ki) &
777	/ (ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
778	! first increase height to h_thick
779	dhsic=MAX(0.,MIN(h_ice_thick-hice(ki),e_melt/(fsic(ki)*ice_den)))
780	hice_i=hice(ki)
781	hice(ki)=dhsic+hice(ki)
782	e_melt=e_melt-fsic(ki)*dhsic
783	IF (e_melt.GT.0.) THEN
784	! frac_mf fraction used for lateral increase
785	dfsic=MIN(amax-fsic(ki),e_meltfrac_mf/ (hice(ki)ice_den) )
786	! No lateral growth -> forced ocean
787	!fsic(ki)=fsic(ki)+dfsic
788	e_melt=e_melt-dfsichice(ki)ice_den
789	! rest used to increase height
790	hice(ki)=MIN(h_ice_max,hice(ki)+e_melt/( fsic(ki) * ice_den ) )
791	END IF
792	dh_basal_growth(ki)=(hice(ki)-hice_i)/dtime
793
794	! melt from below if bottom flux >0
795	ELSE
796	! larger fraction of lateral melt from warm ocean
797	frac_mf=MIN(1.,MAX(0.,(hice(ki)-h_ice_thin) &
798	/ (h_ice_thick-h_ice_thin)))
799	! bring ice to freezing and melt from below
800	! quantity of melted ice
801	e_melt= f_bot(i) * dtime * fsic(ki) &
802	/ (ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
803	! first decrease height to h_thick
804	hice_i=hice(ki)
805	dhsic=MAX(0.,MIN(hice(ki)-h_ice_thick,e_melt/(fsic(ki)*ice_den)))
806	hice(ki)=hice(ki)-dhsic
807	e_melt=e_melt-fsic(ki)dhsicice_den
808
809	IF (e_melt.GT.0) THEN
810	! frac_mf fraction used for height decrease
811	dhsic=MAX(0.,MIN(hice(ki)-h_ice_min,e_melt/ice_den*frac_mf/fsic(ki)))
812	hice(ki)=hice(ki)-dhsic
813	e_melt=e_melt-fsic(ki)dhsicice_den
814	! rest used to decrease fraction (up to 0!)
815	dfsic=MIN(fsic(ki),e_melt/(hice(ki)*ice_den))
816	! Remaining heat not used if everything melted
817	e_melt=e_melt-dfsichice(ki)ice_den
818	! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime
819	END IF
820	dh_basal_melt(ki)=(hice(ki)-hice_i)/dtime
821	END IF
822	END IF
823
824	! melt ice from above if Tice>0
825	tice_i=tice(ki)
826	IF (tice(ki).GT.t_melt) THEN
827	IF (iflag_seaice==2) THEN
828	! quantity of ice melted (kg/m2)
829	e_melt=MAX(hice(ki)ice_den(tice(ki)-t_melt)*ice_cap/2. &
830	/(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
831	! melt from above, height only
832	hice_i=hice(ki)
833	dhsic=MIN(hice(ki)-h_ice_min,e_melt/ice_den)
834	dh_top_melt(i)=dhsic
835	e_melt=e_melt-dhsic
836	IF (e_melt.GT.0) THEN
837	! lateral melt if ice too thin
838	dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/(h_ice_minice_den)fsic(ki))
839	! if all melted do nothing with remaining heat
840	e_melt=MAX(0.,e_meltfsic(ki)-dfsich_ice_min*ice_den)
841	! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime
842	END IF
843	hice(ki)=hice(ki)-dhsic
844	dh_top_melt(ki)=(hice(ki)-hice_i)/dtime
845	! surface temperature at melting point
846	END IF
847	tice(ki)=t_melt
848	tsurf_new(i)=t_melt
849	END IF
850	dtice_melt(ki)=dtice_melt(ki)+tice(ki)-tice_i
851
852	! convert snow to ice if below floating line
853	h_test=(hice(ki)ice_den+snow(i))-hice(ki)sea_den
854	IF ((h_test.GT.0.).AND.(hice(ki).GT.h_ice_min)) THEN !snow under water
855	! extra snow converted to ice (with added frozen sea water)
856	IF (iflag_seaice==2) THEN
857	dhsic=h_test/(sea_den-ice_den+sno_den)
858	hice(ki)=hice(ki)+dhsic
859	ENDIF
860	snow(i)=snow(i)-dhsic*sno_den
861	! available energy (freeze sea water + bring to tice)
862	e_melt=dhsicice_den(1.-sno_den/ice_den)*(ice_lat+ &
863	ice_cap/2.*(t_freeze-tice(ki)))
864	! update ice temperature
865	tice_i=tice(ki)
866	tice(ki)=tice(ki)+2.e_melt/ice_cap/(snow(i)+hice(ki)ice_den)
867	IF (iflag_seaice==2) THEN
868	dh_snow2sic(ki)=dhsic/dtime
869	END IF
870	dtice_snow2sic(ki)=(tice(ki)-tice_i)/dtime
871	END IF
872	END DO
873
874	!write(,) 'hice 2',hice(1:100)
875	!write(,) 'tice 2',tice(1:100)
876
877	iflag_sic_albedo=iflag_seaice_alb
878
879	!*******************************************************************************
880	! 3) cumulate ice-ocean fluxes, update tslab, lateral grow
881	!*********************************************o*****************************
882	!cumul fluxes.
883	bilg_cumul(:)=bilg_cumul(:)+bilg(:)/float(cpl_pas)
884	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab
885	bilg_cumul(:)=0.
886	END IF ! coupling time
887
888	! write(,) 'hice 3',hice(1:100)
889	! write(,) 'tice 3',tice(1:100)
890	!tests ice fraction
891	WHERE (fsic.LT.ice_frac_min)
892	tice=t_melt
893	hice=h_ice_min
894	END WHERE
895
896	!write(,) 'hice 4',hice(1:100)
897	!write(,) 'tice 4',tice(1:100)
898
899	endif
900
901	!****************************************************************************************
902	! 4) Compute sea-ice and snow albedo
903	!****************************************************************************************
904	IF (iflag_seaice > 0) THEN
905	SELECT CASE (iflag_sic_albedo)
906	CASE(0)
907	! old slab albedo : single value. age of snow, melt ponds.
908	DO i=1,knon
909	ki=knindex(i)
910	! snow albedo: update snow age
911	IF (snow(i).GT.0.0001) THEN
912	agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
913	* EXP(-1.MAX(0.0,precip_snow(i))dtime/5.)
914	ELSE
915	agesno(i)=0.0
916	END IF
917	! snow albedo
918	alb_snow=alb_sno_wet+(alb_sno_dry-alb_sno_wet)*EXP(-agesno(i)/50.)
919	! ice albedo (varies with ice tkickness and temp)
920	alb_ice=MAX(0.0,0.13LOG(100.hice(ki))+0.1)
921	!alb_ice=MAX(0.0,0.13LOG(100.seaice(ki)/ice_den)+0.1)
922	IF (tice(ki).GT.t_freeze-0.01) THEN
923	alb_ice=MIN(alb_ice,alb_ice_wet)
924	ELSE
925	alb_ice=MIN(alb_ice,alb_ice_dry)
926	END IF
927	! pond albedo
928	alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
929	! pond fraction
930	frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
931	! snow fraction
932	frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
933	! ice fraction
934	frac_ice=MAX(0.0,1.-frac_pond-frac_snow)
935	! total albedo
936	alb1_new(i)=alb_snowfrac_snow+alb_icefrac_ice+alb_pond*frac_pond
937	END DO
938	alb2_new(:) = alb1_new(:)
939
940	CASE(1)
941	! New visible and IR albedos, dry / melting snow
942	! based on Toyoda et al, 2020
943	DO i=1,knon
944	ki=knindex(i)
945	! snow fraction
946	frac_snow = snow(i) / (snow(i) + h_sno_alb)
947	! dependence of ice albedo with ice thickness
948	frac_ice = MIN(1.,ATAN(4.hice(ki)ice_den) / ATAN(4.*h_ice_alb))
949	! Total (for ice, min = 0.066 = alb_ocean)
950	IF (tice(ki).GT.t_melt) THEN
951	alb_ice = 0.066 + (alb_imlt_vis - 0.066)*frac_ice
952	alb1_new(i)=alb_smlt_visfrac_snow + alb_ice(1.-frac_snow)
953	alb_ice = 0.066 + (alb_imlt_nir - 0.066)*frac_ice
954	alb2_new(i)=alb_smlt_nirfrac_snow + alb_ice(1.-frac_snow)
955	ELSEIF (tice(ki).GT.t_melt - 1.) THEN
956	frac_pond = tice(ki) - t_freeze
957	alb_snow = alb_smlt_visfrac_pond + alb_sdry_vis(1.-frac_pond)
958	alb_ice = alb_imlt_visfrac_pond + alb_idry_vis(1.-frac_pond)
959	alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice
960	alb1_new(i)= alb_snowfrac_snow + alb_ice(1.-frac_snow)
961	alb_snow = alb_smlt_nirfrac_pond + alb_sdry_nir(1.-frac_pond)
962	alb_ice = alb_imlt_nirfrac_pond + alb_idry_nir(1.-frac_pond)
963	alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice
964	alb2_new(i)= alb_snowfrac_snow + alb_ice(1.-frac_snow)
965	ELSE
966	alb_ice = 0.066 + (alb_idry_vis - 0.066)*frac_ice
967	alb1_new(i)=alb_sdry_visfrac_snow + alb_ice(1.-frac_snow)
968	alb_ice = 0.066 + (alb_idry_nir - 0.066)*frac_ice
969	alb2_new(i)=alb_sdry_nirfrac_snow + alb_ice(1.-frac_snow)
970	ENDIF
971	END DO
972
973	CASE(2)
974	! LIM3 scheme. Uses clear sky / overcast value, with 50% clear sky
975	z1_i = 1.5 * ice_den
976	z2_i = 0.05 * ice_den
977	zlog = 1. / (LOG(1.5) - LOG(0.05))
978	z1_s = 1. / (0.025 * sno_den)
979	DO i=1,knon
980	ki=knindex(i)
981	! temperature above / below 0
982	IF (tice(ki).GE.t_melt) THEN
983	alb_ice = alb_ice_wet
984	alb_snow = alb_sno_wet
985	ELSE
986	alb_ice = alb_ice_dry
987	alb_snow = alb_sno_dry
988	ENDIF
989	! ice thickness
990	IF (hice(ki)*ice_den.LT.z2_i) THEN
991	alb_ice = 0.066 + 0.114 * hice(ki)*ice_den / z2_i
992	ELSEIF (hice(ki)*ice_den.LT.z1_i) THEN
993	alb_ice = alb_ice + (0.18 - alb_ice) &
994	* (LOG(z1_i) - LOG(hice(ki)ice_den)) zlog
995	ENDIF
996	! ice or snow depending on snow thickness
997	alb_snow = alb_snow - (alb_snow -alb_ice) * EXP(- snow(i) * z1_s)
998	! Effect of clouds (polynomial fit with 50% clouds)
999	alb1_new(i) = alb_snow - 0.5 * (-0.1010 * alb_snow*alb_snow &
1000	+ 0.1933*alb_snow - 0.0148)
1001	alb2_new(i) = alb1_new(i)
1002	END DO
1003
1004	CASE(3)
1005	CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:))
1006	WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.
1007	alb1_new(:) = 0.0
1008	DO i=1, knon
1009	zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
1010	alb1_new(i) = alb_neig(i) * zfra + 0.6 * (1.0-zfra)
1011	ENDDO
1012	alb2_new(:) = alb1_new(:)
1013
1014	print*,'alb_neig=',alb_neig
1015	print*,'zfra=',zfra
1016	print*,'snow=',snow
1017	print*,'alb1_new=',alb1_new
1018	print*,'alb2_new=',alb2_new
1019	END SELECT
1020	END IF
1021	! ------ End Albedo ----------
1022
1023	!GG
1024	END SUBROUTINE ocean_forced_ice
1025
1026	!************************************************************************
1027	! 1D case
1028	!************************************************************************
1029	! SUBROUTINE read_tsurf1d(knon,sst_out)
1030	!
1031	! This subroutine specifies the surface temperature to be used in 1D simulations
1032	!
1033	! USE dimphy, ONLY : klon
1034	!
1035	! INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
1036	! REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model
1037	!
1038	! INTEGER :: i
1039	! COMMON defined in lmdz1d.F:
1040	! real ts_cur
1041	! common /sst_forcing/ts_cur
1042	!
1043	! DO i = 1, knon
1044	! sst_out(i) = ts_cur
1045	! ENDDO
1046	!
1047	! END SUBROUTINE read_tsurf1d
1048	!
1049	!
1050	!************************************************************************
1051	END MODULE ocean_forced_mod
1052
1053
1054
1055
1056
1057

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

Download in other formats:

Original Format