Context Navigation

ocean_slab_mod.F90 @ 3983

Last change on this file since 3983 was 3817, checked in by millour, 10 years ago
Further cleanup and removal of references to iniprint.h. Also added bench testcase 48x36x19. EM
Property svn:executable set to ``*
File size: 28.6 KB

Line
1	!
2	MODULE ocean_slab_mod
3	!
4	! This module is used for both surface ocean and sea-ice when using the slab ocean,
5	! "ocean=slab".
6	!
7
8	USE dimphy
9	USE indice_sol_mod
10	USE surface_data
11
12	IMPLICIT NONE
13	PRIVATE
14	PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
15
16	!****************************************************************************************
17	! Global saved variables
18	!****************************************************************************************
19
20	INTEGER, PRIVATE, SAVE :: cpl_pas
21	!$OMP THREADPRIVATE(cpl_pas)
22	REAL, PRIVATE, SAVE :: cyang
23	!$OMP THREADPRIVATE(cyang)
24	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh
25	!$OMP THREADPRIVATE(slabh)
26	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: tslab
27	!$OMP THREADPRIVATE(tslab)
28	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic
29	!$OMP THREADPRIVATE(fsic)
30	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice
31	!$OMP THREADPRIVATE(tice)
32	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: seaice
33	!$OMP THREADPRIVATE(seaice)
34	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bils
35	!$OMP THREADPRIVATE(slab_bils)
36	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum
37	!$OMP THREADPRIVATE(bils_cum)
38	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bilg
39	!$OMP THREADPRIVATE(slab_bilg)
40	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bilg_cum
41	!$OMP THREADPRIVATE(bilg_cum)
42
43	!****************************************************************************************
44	! Parameters (could be read in def file: move to slab_init)
45	!****************************************************************************************
46	! snow and ice physical characteristics:
47	REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
48	REAL, PARAMETER :: t_melt=273.15 ! melting ice temp
49	REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
50	REAL, PARAMETER :: ice_den=917.
51	REAL, PARAMETER :: sea_den=1025.
52	REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity
53	REAL, PARAMETER :: sno_cond=0.31*sno_den
54	REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice
55	REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
56
57	! control of snow and ice cover & freeze / melt (heights converted to kg/m2)
58	REAL, PARAMETER :: snow_min=0.05*sno_den !critical snow height 5 cm
59	REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean
60	REAL, PARAMETER :: ice_frac_min=0.001
61	REAL, PARAMETER :: ice_frac_max=1. ! less than 1. if min leads fraction
62	REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness
63	REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness
64	! below ice_thin, priority is melt lateral / grow height
65	! ice_thin is also height of new ice
66	REAL, PARAMETER :: h_ice_thick=2.5*ice_den ! thin ice thickness
67	! above ice_thick, priority is melt height / grow lateral
68	REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice
69	REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height
70
71	! albedo and radiation parameters
72	REAL, PARAMETER :: alb_sno_min=0.55 !min snow albedo
73	REAL, PARAMETER :: alb_sno_del=0.3 !max snow albedo = min + del
74	REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
75	REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
76	REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow)
77	REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
78
79	!****************************************************************************************
80
81	CONTAINS
82	!
83	!****************************************************************************************
84	!
85	SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
86	!, seaice_rst etc
87
88	use IOIPSL
89
90	! For ok_xxx vars (Ekman...)
91	INCLUDE "clesphys.h"
92
93	! Input variables
94	!****************************************************************************************
95	REAL, INTENT(IN) :: dtime
96	! Variables read from restart file
97	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
98
99	! Local variables
100	!****************************************************************************************
101	INTEGER :: error
102	CHARACTER (len = 80) :: abort_message
103	CHARACTER (len = 20) :: modname = 'ocean_slab_intit'
104
105	!****************************************************************************************
106	! Allocate surface fraction read from restart file
107	!****************************************************************************************
108	ALLOCATE(fsic(klon), stat = error)
109	IF (error /= 0) THEN
110	abort_message='Pb allocation tmp_pctsrf_slab'
111	CALL abort_physic(modname,abort_message,1)
112	ENDIF
113	fsic(:)=0.
114	WHERE (1.-zmasq(:)>EPSFRA)
115	fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
116	END WHERE
117
118	!****************************************************************************************
119	! Allocate saved variables
120	!****************************************************************************************
121	ALLOCATE(tslab(klon,nslay), stat=error)
122	IF (error /= 0) CALL abort_physic &
123	(modname,'pb allocation tslab', 1)
124
125	ALLOCATE(slab_bils(klon), stat = error)
126	IF (error /= 0) THEN
127	abort_message='Pb allocation slab_bils'
128	CALL abort_physic(modname,abort_message,1)
129	ENDIF
130	slab_bils(:) = 0.0
131	ALLOCATE(bils_cum(klon), stat = error)
132	IF (error /= 0) THEN
133	abort_message='Pb allocation slab_bils_cum'
134	CALL abort_physic(modname,abort_message,1)
135	ENDIF
136	bils_cum(:) = 0.0
137
138	IF (version_ocean=='sicINT') THEN
139	ALLOCATE(slab_bilg(klon), stat = error)
140	IF (error /= 0) THEN
141	abort_message='Pb allocation slab_bilg'
142	CALL abort_physic(modname,abort_message,1)
143	ENDIF
144	slab_bilg(:) = 0.0
145	ALLOCATE(bilg_cum(klon), stat = error)
146	IF (error /= 0) THEN
147	abort_message='Pb allocation slab_bilg_cum'
148	CALL abort_physic(modname,abort_message,1)
149	ENDIF
150	bilg_cum(:) = 0.0
151	ALLOCATE(tice(klon), stat = error)
152	IF (error /= 0) THEN
153	abort_message='Pb allocation slab_tice'
154	CALL abort_physic(modname,abort_message,1)
155	ENDIF
156	ALLOCATE(seaice(klon), stat = error)
157	IF (error /= 0) THEN
158	abort_message='Pb allocation slab_seaice'
159	CALL abort_physic(modname,abort_message,1)
160	ENDIF
161	END IF
162
163	!****************************************************************************************
164	! Define some parameters
165	!****************************************************************************************
166	! Layer thickness
167	ALLOCATE(slabh(nslay), stat = error)
168	IF (error /= 0) THEN
169	abort_message='Pb allocation slabh'
170	CALL abort_physic(modname,abort_message,1)
171	ENDIF
172	slabh(1)=50.
173	! cyang = 1/heat capacity of top layer (rho.c.H)
174	cyang=1/(slabh(1)*4.228e+06)
175
176	! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
177	! pour un calcul à chaque pas de temps, cpl_pas=1
178	cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
179	CALL getin('cpl_pas',cpl_pas)
180	print *,'cpl_pas',cpl_pas
181
182	END SUBROUTINE ocean_slab_init
183	!
184	!****************************************************************************************
185	!
186	SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
187
188	USE limit_read_mod
189
190	! INCLUDE "clesphys.h"
191
192	! Arguments
193	!****************************************************************************************
194	INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant
195	INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee
196	REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s)
197	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg ! sub-surface fraction
198	LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step
199
200	! Local variables
201	!****************************************************************************************
202
203	IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN
204	CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
205	ELSE
206	pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
207	pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
208	is_modified=.TRUE.
209	END IF
210
211	END SUBROUTINE ocean_slab_frac
212	!
213	!****************************************************************************************
214	!
215	SUBROUTINE ocean_slab_noice( &
216	itime, dtime, jour, knon, knindex, &
217	p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
218	AcoefH, AcoefQ, BcoefH, BcoefQ, &
219	AcoefU, AcoefV, BcoefU, BcoefV, &
220	ps, u1, v1, gustiness, tsurf_in, &
221	radsol, snow, &
222	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
223	tsurf_new, dflux_s, dflux_l, qflux)
224
225	USE calcul_fluxs_mod
226
227	! Input arguments
228	!****************************************************************************************
229	INTEGER, INTENT(IN) :: itime
230	INTEGER, INTENT(IN) :: jour
231	INTEGER, INTENT(IN) :: knon
232	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
233	REAL, INTENT(IN) :: dtime
234	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
235	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
236	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
237	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
238	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
239	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
240	REAL, DIMENSION(klon), INTENT(IN) :: ps
241	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
242	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
243	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
244
245	! In/Output arguments
246	!****************************************************************************************
247	REAL, DIMENSION(klon), INTENT(INOUT) :: snow
248
249	! Output arguments
250	!****************************************************************************************
251	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
252	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
253	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
254	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
255	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
256	REAL, DIMENSION(klon), INTENT(OUT) :: qflux
257
258	! Local variables
259	!****************************************************************************************
260	INTEGER :: i,ki
261	! surface fluxes
262	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
263	! flux correction
264	REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
265	! surface wind stress
266	REAL, DIMENSION(klon) :: u0, v0
267	REAL, DIMENSION(klon) :: u1_lay, v1_lay
268	! ice creation
269	REAL :: e_freeze, h_new, dfsic
270
271	!****************************************************************************************
272	! 1) Flux calculation
273	!
274	!****************************************************************************************
275	cal(:) = 0.
276	beta(:) = 1.
277	dif_grnd(:) = 0.
278
279	! Suppose zero surface speed
280	u0(:)=0.0
281	v0(:)=0.0
282	u1_lay(:) = u1(:) - u0(:)
283	v1_lay(:) = v1(:) - v0(:)
284
285	CALL calcul_fluxs(knon, is_oce, dtime, &
286	tsurf_in, p1lay, cal, beta, cdragh, ps, &
287	precip_rain, precip_snow, snow, qsurf, &
288	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
289	AcoefH, AcoefQ, BcoefH, BcoefQ, &
290	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
291
292	! - Flux calculation at first modele level for U and V
293	CALL calcul_flux_wind(knon, dtime, &
294	u0, v0, u1, v1, gustiness, cdragm, &
295	AcoefU, AcoefV, BcoefU, BcoefV, &
296	p1lay, temp_air, &
297	flux_u1, flux_v1)
298
299	! Accumulate total fluxes locally
300	slab_bils(:)=0.
301	DO i=1,knon
302	ki=knindex(i)
303	slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
304	-precip_snow(i)ice_lat(1.+snow_wfact*fsic(ki)))
305	bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
306	! Also taux, tauy, saved vars...
307	END DO
308
309	!****************************************************************************************
310	! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
311	!
312	!****************************************************************************************
313	lmt_bils(:)=0.
314	CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
315	! lmt_bils and diff_sst,siv saved by limit_slab
316	qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)ice_denice_lat/86400.
317	! qflux = total QFlux correction (in W/m2)
318
319	!****************************************************************************************
320	! 3) Recalculate new temperature
321	!
322	!*********************************************o***************************************
323	tsurf_new=tsurf_in
324	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
325	! Compute transport
326	! Add read QFlux and SST tendency
327	tslab(:,1)=tslab(:,1)+qflux(:)cyangdtime*cpl_pas
328	! Add cumulated surface fluxes
329	tslab(:,1)=tslab(:,1)+bils_cum(:)cyangdtime
330	! Update surface temperature
331	SELECT CASE(version_ocean)
332	CASE('sicNO')
333	DO i=1,knon
334	ki=knindex(i)
335	tsurf_new(i)=tslab(ki,1)
336	END DO
337	CASE('sicOBS') ! check for sea ice or tslab below freezing
338	DO i=1,knon
339	ki=knindex(i)
340	IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
341	tslab(ki,1)=t_freeze
342	END IF
343	tsurf_new(i)=tslab(ki,1)
344	END DO
345	CASE('sicINT')
346	DO i=1,knon
347	ki=knindex(i)
348	IF (fsic(ki).LT.epsfra) THEN ! Free of ice
349	IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
350	! quantity of new ice formed
351	e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
352	! new ice
353	tice(ki)=t_freeze
354	fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
355	IF (fsic(ki).GT.ice_frac_min) THEN
356	seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
357	tslab(ki,1)=t_freeze
358	ELSE
359	fsic(ki)=0.
360	END IF
361	tsurf_new(i)=t_freeze
362	ELSE
363	tsurf_new(i)=tslab(ki,1)
364	END IF
365	ELSE ! ice present
366	tsurf_new(i)=t_freeze
367	IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
368	! quantity of new ice formed over open ocean
369	e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
370	/(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
371	! new ice height and fraction
372	h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
373	dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
374	h_new=MIN(e_freeze/dfsic,h_ice_max)
375	! update tslab to freezing over open ocean only
376	tslab(ki,1)=tslab(ki,1)fsic(ki)+t_freeze(1.-fsic(ki))
377	! update sea ice
378	seaice(ki)=(h_newdfsic+seaice(ki)fsic(ki)) &
379	/(dfsic+fsic(ki))
380	fsic(ki)=fsic(ki)+dfsic
381	! update snow?
382	END IF !freezing
383	END IF ! sea ice present
384	END DO
385	END SELECT
386	bils_cum(:)=0.0! clear cumulated fluxes
387	END IF ! coupling time
388	END SUBROUTINE ocean_slab_noice
389	!
390	!****************************************************************************************
391
392	SUBROUTINE ocean_slab_ice( &
393	itime, dtime, jour, knon, knindex, &
394	tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
395	AcoefH, AcoefQ, BcoefH, BcoefQ, &
396	AcoefU, AcoefV, BcoefU, BcoefV, &
397	ps, u1, v1, gustiness, &
398	radsol, snow, qsurf, qsol, agesno, &
399	alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
400	tsurf_new, dflux_s, dflux_l, swnet)
401
402	USE calcul_fluxs_mod
403
404	INCLUDE "YOMCST.h"
405
406	! Input arguments
407	!****************************************************************************************
408	INTEGER, INTENT(IN) :: itime, jour, knon
409	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
410	REAL, INTENT(IN) :: dtime
411	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
412	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
413	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
414	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
415	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
416	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
417	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
418	REAL, DIMENSION(klon), INTENT(IN) :: ps
419	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
420	REAL, DIMENSION(klon), INTENT(IN) :: swnet
421
422	! In/Output arguments
423	!****************************************************************************************
424	REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
425	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
426	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
427
428	! Output arguments
429	!****************************************************************************************
430	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
431	REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval
432	REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval
433	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
434	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
435	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
436	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
437
438	! Local variables
439	!****************************************************************************************
440	INTEGER :: i,ki
441	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
442	REAL, DIMENSION(klon) :: u0, v0
443	REAL, DIMENSION(klon) :: u1_lay, v1_lay
444	! intermediate heat fluxes:
445	REAL :: f_cond, f_swpen
446	! for snow/ice albedo:
447	REAL :: alb_snow, alb_ice, alb_pond
448	REAL :: frac_snow, frac_ice, frac_pond
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
454	!****************************************************************************************
455	! 1) Flux calculation
456	!****************************************************************************************
457	! Suppose zero surface speed
458	u0(:)=0.0
459	v0(:)=0.0
460	u1_lay(:) = u1(:) - u0(:)
461	v1_lay(:) = v1(:) - v0(:)
462
463	! set beta, cal, compute conduction fluxes inside ice/snow
464	slab_bilg(:)=0.
465	dif_grnd(:)=0.
466	beta(:) = 1.
467	DO i=1,knon
468	ki=knindex(i)
469	IF (snow(i).GT.snow_min) THEN
470	! snow-layer heat capacity
471	cal(i)=2.RCPD/(snow(i)ice_cap)
472	! snow conductive flux
473	f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
474	! all shortwave flux absorbed
475	f_swpen=0.
476	! bottom flux (ice conduction)
477	slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
478	! update ice temperature
479	tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
480	(slab_bilg(ki)+f_cond)dtime
481	ELSE ! bare ice
482	! ice-layer heat capacity
483	cal(i)=2.RCPD/(seaice(ki)ice_cap)
484	! conductive flux
485	f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
486	! penetrative shortwave flux...
487	f_swpen=swnet(i)pen_fracexp(-pen_ext*seaice(ki)/ice_den)
488	slab_bilg(ki)=f_swpen-f_cond
489	END IF
490	radsol(i)=radsol(i)+f_cond-f_swpen
491	END DO
492	! weight fluxes to ocean by sea ice fraction
493	slab_bilg(:)=slab_bilg(:)*fsic(:)
494
495	! calcul_fluxs (sens, lat etc)
496	CALL calcul_fluxs(knon, is_sic, dtime, &
497	tsurf_in, p1lay, cal, beta, cdragh, ps, &
498	precip_rain, precip_snow, snow, qsurf, &
499	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
500	AcoefH, AcoefQ, BcoefH, BcoefQ, &
501	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
502	DO i=1,knon
503	IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
504	END DO
505
506	! calcul_flux_wind
507	CALL calcul_flux_wind(knon, dtime, &
508	u0, v0, u1, v1, gustiness, cdragm, &
509	AcoefU, AcoefV, BcoefU, BcoefV, &
510	p1lay, temp_air, &
511	flux_u1, flux_v1)
512
513	!****************************************************************************************
514	! 2) Update snow and ice surface
515	!****************************************************************************************
516	! snow precip
517	DO i=1,knon
518	ki=knindex(i)
519	IF (precip_snow(i) > 0.) THEN
520	snow(i) = snow(i)+precip_snow(i)dtime(1.-snow_wfact*(1.-fsic(ki)))
521	END IF
522	! snow and ice sublimation
523	IF (evap(i) > 0.) THEN
524	snow_evap = MIN (snow(i) / dtime, evap(i))
525	snow(i) = snow(i) - snow_evap * dtime
526	snow(i) = MAX(0.0, snow(i))
527	seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
528	ENDIF
529	! Melt / Freeze from above if Tsurf>0
530	IF (tsurf_new(i).GT.t_melt) THEN
531	! energy available for melting snow (in kg/m2 of snow)
532	e_melt = MIN(MAX(snow(i)(tsurf_new(i)-t_melt)ice_cap/2. &
533	/(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
534	! remove snow
535	IF (snow(i).GT.e_melt) THEN
536	snow(i)=snow(i)-e_melt
537	tsurf_new(i)=t_melt
538	ELSE ! all snow is melted
539	! add remaining heat flux to ice
540	e_melt=e_melt-snow(i)
541	tice(ki)=tice(ki)+e_meltice_lat2./(ice_cap*seaice(ki))
542	tsurf_new(i)=tice(ki)
543	END IF
544	END IF
545	! melt ice from above if Tice>0
546	IF (tice(ki).GT.t_melt) THEN
547	! quantity of ice melted (kg/m2)
548	e_melt=MAX(seaice(ki)(tice(ki)-t_melt)ice_cap/2. &
549	/(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
550	! melt from above, height only
551	dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
552	e_melt=e_melt-dhsic
553	IF (e_melt.GT.0) THEN
554	! lateral melt if ice too thin
555	dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
556	! if all melted add remaining heat to ocean
557	e_melt=MAX(0.,e_meltfsic(ki)-dfsich_ice_min)
558	slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime
559	! update height and fraction
560	fsic(ki)=fsic(ki)-dfsic
561	END IF
562	seaice(ki)=seaice(ki)-dhsic
563	! surface temperature at melting point
564	tice(ki)=t_melt
565	tsurf_new(i)=t_melt
566	END IF
567	! convert snow to ice if below floating line
568	h_test=(seaice(ki)+snow(i))ice_den-seaice(ki)sea_den
569	IF (h_test.GT.0.) THEN !snow under water
570	! extra snow converted to ice (with added frozen sea water)
571	dhsic=h_test/(sea_den-ice_den+sno_den)
572	seaice(ki)=seaice(ki)+dhsic
573	snow(i)=snow(i)-dhsic*sno_den/ice_den
574	! available energy (freeze sea water + bring to tice)
575	e_melt=dhsic(1.-sno_den/ice_den)(ice_lat+ &
576	ice_cap/2.*(t_freeze-tice(ki)))
577	! update ice temperature
578	tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
579	END IF
580	END DO
581
582	! New albedo
583	DO i=1,knon
584	ki=knindex(i)
585	! snow albedo: update snow age
586	IF (snow(i).GT.0.0001) THEN
587	agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
588	* EXP(-1.MAX(0.0,precip_snow(i))dtime/5.)
589	ELSE
590	agesno(i)=0.0
591	END IF
592	! snow albedo
593	alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.)
594	! ice albedo (varies with ice tkickness and temp)
595	alb_ice=MAX(0.0,0.13LOG(100.seaice(ki)/ice_den)+0.1)
596	IF (tice(ki).GT.t_freeze-0.01) THEN
597	alb_ice=MIN(alb_ice,alb_ice_wet)
598	ELSE
599	alb_ice=MIN(alb_ice,alb_ice_dry)
600	END IF
601	! pond albedo
602	alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
603	! pond fraction
604	frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
605	! snow fraction
606	frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
607	! ice fraction
608	frac_ice=MAX(0.0,1.-frac_pond-frac_snow)
609	! total albedo
610	alb1_new(i)=alb_snowfrac_snow+alb_icefrac_ice+alb_pond*frac_pond
611	END DO
612	alb2_new(:) = alb1_new(:)
613
614	!****************************************************************************************
615	! 3) Recalculate new ocean temperature (add fluxes below ice)
616	! Melt / freeze from below
617	!*********************************************o***************************************
618	!cumul fluxes
619	bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
620	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
621	! Add cumulated surface fluxes
622	tslab(:,1)=tslab(:,1)+bilg_cum(:)cyangdtime
623	DO i=1,knon
624	ki=knindex(i)
625	! split lateral/top melt-freeze
626	frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin)))
627	IF (tslab(ki,1).LE.t_freeze) THEN
628	! **** Form new ice from below *****
629	! quantity of new ice
630	e_melt=(t_freeze-tslab(ki,1))/cyang &
631	/(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
632	! first increase height to h_thin
633	dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
634	seaice(ki)=dhsic+seaice(ki)
635	e_melt=e_melt-fsic(ki)*dhsic
636	IF (e_melt.GT.0.) THEN
637	! frac_mf fraction used for lateral increase
638	dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
639	fsic(ki)=fsic(ki)+dfsic
640	e_melt=e_melt-dfsic*seaice(ki)
641	! rest used to increase height
642	seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
643	END IF
644	tslab(ki,1)=t_freeze
645	ELSE ! slab temperature above freezing
646	! **** melt ice from below *****
647	! quantity of melted ice
648	e_melt=(tslab(ki,1)-t_freeze)/cyang &
649	/(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
650	! first decrease height to h_thick
651	dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
652	seaice(ki)=seaice(ki)-dhsic
653	e_melt=e_melt-fsic(ki)*dhsic
654	IF (e_melt.GT.0) THEN
655	! frac_mf fraction used for height decrease
656	dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
657	seaice(ki)=seaice(ki)-dhsic
658	e_melt=e_melt-fsic(ki)*dhsic
659	! rest used to decrease fraction (up to 0!)
660	dfsic=MIN(fsic(ki),e_melt/seaice(ki))
661	! keep remaining in ocean
662	e_melt=e_melt-dfsic*seaice(ki)
663	END IF
664	tslab(ki,1)=t_freeze+e_meltice_latcyang
665	fsic(ki)=fsic(ki)-dfsic
666	END IF
667	END DO
668	bilg_cum(:)=0.
669	END IF ! coupling time
670
671	!tests fraction glace
672	WHERE (fsic.LT.ice_frac_min)
673	tslab(:,1)=tslab(:,1)-fsicseaiceice_lat*cyang
674	tice=t_melt
675	fsic=0.
676	seaice=0.
677	END WHERE
678
679	END SUBROUTINE ocean_slab_ice
680	!
681	!****************************************************************************************
682	!
683	SUBROUTINE ocean_slab_final
684
685	!****************************************************************************************
686	! Deallocate module variables
687	!****************************************************************************************
688	IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
689	IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
690	IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
691	IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
692	IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
693	IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
694	IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
695
696	END SUBROUTINE ocean_slab_final
697
698	END MODULE ocean_slab_mod

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

Download in other formats:

Original Format