Context Navigation

ocean_forced_mod.F90 @ 5932

Last change on this file since 5932 was 5896, checked in by yann meurdesoif, 5 weeks ago

GPU port of surf_ocean

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: 43.1 KB

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

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

Download in other formats:

Original Format