Context Navigation

ocean_slab_mod.F90 @ 2254

Last change on this file since 2254 was 2254, checked in by fhourdin, 9 years ago

Modification du calcul des flux air/mer
1) Introduction d'un facteur f_qsat_oce=0.98 devant qsat dans le calcul
de l'évaporation sur océan pour tenir compte de la moindre évaporation
de l'eau salée.
2) Introduction d'une différentiation entre z0 pour le sensible, z0h,
et le latent, z0q, imposé constant z0q=f_z0qh_oce*z0h

Modification of air/sea fluxes computation
1) Introduction of a correcting factor f_qsat_oce=0.98 on qsat
to account for the weaker evaporation of salty water.
2) Introduction of z0q=f_z0qh_oce*z0h

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:executable set to *
Property svn:keywords set to Author Date Id Revision

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

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

Download in other formats:

Original Format