Context Navigation

ocean_slab_mod.F90 @ 2241

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

Revisite de la formule des flux de surface
(en priorité sur l'océan) en tenant compte des bourrasques de
vent et de la différence entre les hauteurs de rugosités pour
la quantité de mouvement, l'enthalpie et éventuellement l'humidité.

Etape 1 :
Introduction d'un calcul de gustiness dans la physique
gustiness(:)=f_gust_bl * ale_bl + f_gust_wk * ame_wk
Cette variable est passée ensuite jusqu'au fin fond de la couche limite.
L'étape 1 est prête à commettre, ne nécessite pas de nouvelles
variables dans les startphy et assure la convergence numérique.

Introduction of gustiness in the surface flux computation.
Gustiness is computed from as
gustiness(:)=f_gust_bl * ale_bl + f_gust_wk * ame_wk
and pass through pbl_surface down to the routines that compute
surface fluxes.

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.7 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, 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
230	! Input arguments
231	!****************************************************************************************
232	INTEGER, INTENT(IN) :: itime
233	INTEGER, INTENT(IN) :: jour
234	INTEGER, INTENT(IN) :: knon
235	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
236	REAL, INTENT(IN) :: dtime
237	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
238	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
239	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
240	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
241	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
242	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
243	REAL, DIMENSION(klon), INTENT(IN) :: ps
244	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
245	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
246	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
247
248	! In/Output arguments
249	!****************************************************************************************
250	REAL, DIMENSION(klon), INTENT(INOUT) :: snow
251
252	! Output arguments
253	!****************************************************************************************
254	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
255	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
256	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
257	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
258	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
259	REAL, DIMENSION(klon), INTENT(OUT) :: qflux
260
261	! Local variables
262	!****************************************************************************************
263	INTEGER :: i,ki
264	! surface fluxes
265	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
266	! flux correction
267	REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
268	! surface wind stress
269	REAL, DIMENSION(klon) :: u0, v0
270	REAL, DIMENSION(klon) :: u1_lay, v1_lay
271	! ice creation
272	REAL :: e_freeze, h_new, dfsic
273
274	!****************************************************************************************
275	! 1) Flux calculation
276	!
277	!****************************************************************************************
278	cal(:) = 0.
279	beta(:) = 1.
280	dif_grnd(:) = 0.
281
282	! Suppose zero surface speed
283	u0(:)=0.0
284	v0(:)=0.0
285	u1_lay(:) = u1(:) - u0(:)
286	v1_lay(:) = v1(:) - v0(:)
287
288	CALL calcul_fluxs(knon, is_oce, dtime, &
289	tsurf_in, p1lay, cal, beta, cdragh, ps, &
290	precip_rain, precip_snow, snow, qsurf, &
291	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
292	AcoefH, AcoefQ, BcoefH, BcoefQ, &
293	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
294
295	! - Flux calculation at first modele level for U and V
296	CALL calcul_flux_wind(knon, dtime, &
297	u0, v0, u1, v1, gustiness, cdragm, &
298	AcoefU, AcoefV, BcoefU, BcoefV, &
299	p1lay, temp_air, &
300	flux_u1, flux_v1)
301
302	! Accumulate total fluxes locally
303	slab_bils(:)=0.
304	DO i=1,knon
305	ki=knindex(i)
306	slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
307	-precip_snow(i)ice_lat(1.+snow_wfact*fsic(ki)))
308	bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
309	! Also taux, tauy, saved vars...
310	END DO
311
312	!****************************************************************************************
313	! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
314	!
315	!****************************************************************************************
316	lmt_bils(:)=0.
317	CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
318	! lmt_bils and diff_sst,siv saved by limit_slab
319	qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)ice_denice_lat/86400.
320	! qflux = total QFlux correction (in W/m2)
321
322	!****************************************************************************************
323	! 3) Recalculate new temperature
324	!
325	!*********************************************o***************************************
326	tsurf_new=tsurf_in
327	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
328	! Compute transport
329	! Add read QFlux and SST tendency
330	tslab(:,1)=tslab(:,1)+qflux(:)cyangdtime*cpl_pas
331	! Add cumulated surface fluxes
332	tslab(:,1)=tslab(:,1)+bils_cum(:)cyangdtime
333	! Update surface temperature
334	SELECT CASE(version_ocean)
335	CASE('sicNO')
336	DO i=1,knon
337	ki=knindex(i)
338	tsurf_new(i)=tslab(ki,1)
339	END DO
340	CASE('sicOBS') ! check for sea ice or tslab below freezing
341	DO i=1,knon
342	ki=knindex(i)
343	IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
344	tslab(ki,1)=t_freeze
345	END IF
346	tsurf_new(i)=tslab(ki,1)
347	END DO
348	CASE('sicINT')
349	DO i=1,knon
350	ki=knindex(i)
351	IF (fsic(ki).LT.epsfra) THEN ! Free of ice
352	IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
353	! quantity of new ice formed
354	e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
355	! new ice
356	tice(ki)=t_freeze
357	fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
358	IF (fsic(ki).GT.ice_frac_min) THEN
359	seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
360	tslab(ki,1)=t_freeze
361	ELSE
362	fsic(ki)=0.
363	END IF
364	tsurf_new(i)=t_freeze
365	ELSE
366	tsurf_new(i)=tslab(ki,1)
367	END IF
368	ELSE ! ice present
369	tsurf_new(i)=t_freeze
370	IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
371	! quantity of new ice formed over open ocean
372	e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
373	/(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
374	! new ice height and fraction
375	h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
376	dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
377	h_new=MIN(e_freeze/dfsic,h_ice_max)
378	! update tslab to freezing over open ocean only
379	tslab(ki,1)=tslab(ki,1)fsic(ki)+t_freeze(1.-fsic(ki))
380	! update sea ice
381	seaice(ki)=(h_newdfsic+seaice(ki)fsic(ki)) &
382	/(dfsic+fsic(ki))
383	fsic(ki)=fsic(ki)+dfsic
384	! update snow?
385	END IF !freezing
386	END IF ! sea ice present
387	END DO
388	END SELECT
389	bils_cum(:)=0.0! clear cumulated fluxes
390	END IF ! coupling time
391	END SUBROUTINE ocean_slab_noice
392	!
393	!****************************************************************************************
394
395	SUBROUTINE ocean_slab_ice( &
396	itime, dtime, jour, knon, knindex, &
397	tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
398	AcoefH, AcoefQ, BcoefH, BcoefQ, &
399	AcoefU, AcoefV, BcoefU, BcoefV, &
400	ps, u1, v1, gustiness, &
401	radsol, snow, qsurf, qsol, agesno, &
402	alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
403	tsurf_new, dflux_s, dflux_l, swnet)
404
405	USE calcul_fluxs_mod
406
407	INCLUDE "YOMCST.h"
408
409	! Input arguments
410	!****************************************************************************************
411	INTEGER, INTENT(IN) :: itime, jour, knon
412	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
413	REAL, INTENT(IN) :: dtime
414	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
415	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
416	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
417	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
418	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
419	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
420	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
421	REAL, DIMENSION(klon), INTENT(IN) :: ps
422	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
423	REAL, DIMENSION(klon), INTENT(IN) :: swnet
424
425	! In/Output arguments
426	!****************************************************************************************
427	REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
428	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
429	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
430
431	! Output arguments
432	!****************************************************************************************
433	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
434	REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval
435	REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval
436	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
437	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
438	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
439	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
440
441	! Local variables
442	!****************************************************************************************
443	INTEGER :: i,ki
444	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
445	REAL, DIMENSION(klon) :: u0, v0
446	REAL, DIMENSION(klon) :: u1_lay, v1_lay
447	! intermediate heat fluxes:
448	REAL :: f_cond, f_swpen
449	! for snow/ice albedo:
450	REAL :: alb_snow, alb_ice, alb_pond
451	REAL :: frac_snow, frac_ice, frac_pond
452	! for ice melt / freeze
453	REAL :: e_melt, snow_evap, h_test
454	! dhsic, dfsic change in ice mass, fraction.
455	REAL :: dhsic, dfsic, frac_mf
456
457	!****************************************************************************************
458	! 1) Flux calculation
459	!****************************************************************************************
460	! Suppose zero surface speed
461	u0(:)=0.0
462	v0(:)=0.0
463	u1_lay(:) = u1(:) - u0(:)
464	v1_lay(:) = v1(:) - v0(:)
465
466	! set beta, cal, compute conduction fluxes inside ice/snow
467	slab_bilg(:)=0.
468	dif_grnd(:)=0.
469	beta(:) = 1.
470	DO i=1,knon
471	ki=knindex(i)
472	IF (snow(i).GT.snow_min) THEN
473	! snow-layer heat capacity
474	cal(i)=2.RCPD/(snow(i)ice_cap)
475	! snow conductive flux
476	f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
477	! all shortwave flux absorbed
478	f_swpen=0.
479	! bottom flux (ice conduction)
480	slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
481	! update ice temperature
482	tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
483	(slab_bilg(ki)+f_cond)dtime
484	ELSE ! bare ice
485	! ice-layer heat capacity
486	cal(i)=2.RCPD/(seaice(ki)ice_cap)
487	! conductive flux
488	f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
489	! penetrative shortwave flux...
490	f_swpen=swnet(i)pen_fracexp(-pen_ext*seaice(ki)/ice_den)
491	slab_bilg(ki)=f_swpen-f_cond
492	END IF
493	radsol(i)=radsol(i)+f_cond-f_swpen
494	END DO
495	! weight fluxes to ocean by sea ice fraction
496	slab_bilg(:)=slab_bilg(:)*fsic(:)
497
498	! calcul_fluxs (sens, lat etc)
499	CALL calcul_fluxs(knon, is_sic, dtime, &
500	tsurf_in, p1lay, cal, beta, cdragh, ps, &
501	precip_rain, precip_snow, snow, qsurf, &
502	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
503	AcoefH, AcoefQ, BcoefH, BcoefQ, &
504	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
505	DO i=1,knon
506	IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
507	END DO
508
509	! calcul_flux_wind
510	CALL calcul_flux_wind(knon, dtime, &
511	u0, v0, u1, v1, gustiness, cdragm, &
512	AcoefU, AcoefV, BcoefU, BcoefV, &
513	p1lay, temp_air, &
514	flux_u1, flux_v1)
515
516	!****************************************************************************************
517	! 2) Update snow and ice surface
518	!****************************************************************************************
519	! snow precip
520	DO i=1,knon
521	ki=knindex(i)
522	IF (precip_snow(i) > 0.) THEN
523	snow(i) = snow(i)+precip_snow(i)dtime(1.-snow_wfact*(1.-fsic(ki)))
524	END IF
525	! snow and ice sublimation
526	IF (evap(i) > 0.) THEN
527	snow_evap = MIN (snow(i) / dtime, evap(i))
528	snow(i) = snow(i) - snow_evap * dtime
529	snow(i) = MAX(0.0, snow(i))
530	seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
531	ENDIF
532	! Melt / Freeze from above if Tsurf>0
533	IF (tsurf_new(i).GT.t_melt) THEN
534	! energy available for melting snow (in kg/m2 of snow)
535	e_melt = MIN(MAX(snow(i)(tsurf_new(i)-t_melt)ice_cap/2. &
536	/(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
537	! remove snow
538	IF (snow(i).GT.e_melt) THEN
539	snow(i)=snow(i)-e_melt
540	tsurf_new(i)=t_melt
541	ELSE ! all snow is melted
542	! add remaining heat flux to ice
543	e_melt=e_melt-snow(i)
544	tice(ki)=tice(ki)+e_meltice_lat2./(ice_cap*seaice(ki))
545	tsurf_new(i)=tice(ki)
546	END IF
547	END IF
548	! melt ice from above if Tice>0
549	IF (tice(ki).GT.t_melt) THEN
550	! quantity of ice melted (kg/m2)
551	e_melt=MAX(seaice(ki)(tice(ki)-t_melt)ice_cap/2. &
552	/(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
553	! melt from above, height only
554	dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
555	e_melt=e_melt-dhsic
556	IF (e_melt.GT.0) THEN
557	! lateral melt if ice too thin
558	dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
559	! if all melted add remaining heat to ocean
560	e_melt=MAX(0.,e_meltfsic(ki)-dfsich_ice_min)
561	slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime
562	! update height and fraction
563	fsic(ki)=fsic(ki)-dfsic
564	END IF
565	seaice(ki)=seaice(ki)-dhsic
566	! surface temperature at melting point
567	tice(ki)=t_melt
568	tsurf_new(i)=t_melt
569	END IF
570	! convert snow to ice if below floating line
571	h_test=(seaice(ki)+snow(i))ice_den-seaice(ki)sea_den
572	IF (h_test.GT.0.) THEN !snow under water
573	! extra snow converted to ice (with added frozen sea water)
574	dhsic=h_test/(sea_den-ice_den+sno_den)
575	seaice(ki)=seaice(ki)+dhsic
576	snow(i)=snow(i)-dhsic*sno_den/ice_den
577	! available energy (freeze sea water + bring to tice)
578	e_melt=dhsic(1.-sno_den/ice_den)(ice_lat+ &
579	ice_cap/2.*(t_freeze-tice(ki)))
580	! update ice temperature
581	tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
582	END IF
583	END DO
584
585	! New albedo
586	DO i=1,knon
587	ki=knindex(i)
588	! snow albedo: update snow age
589	IF (snow(i).GT.0.0001) THEN
590	agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
591	* EXP(-1.MAX(0.0,precip_snow(i))dtime/5.)
592	ELSE
593	agesno(i)=0.0
594	END IF
595	! snow albedo
596	alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.)
597	! ice albedo (varies with ice tkickness and temp)
598	alb_ice=MAX(0.0,0.13LOG(100.seaice(ki)/ice_den)+0.1)
599	IF (tice(ki).GT.t_freeze-0.01) THEN
600	alb_ice=MIN(alb_ice,alb_ice_wet)
601	ELSE
602	alb_ice=MIN(alb_ice,alb_ice_dry)
603	END IF
604	! pond albedo
605	alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
606	! pond fraction
607	frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
608	! snow fraction
609	frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
610	! ice fraction
611	frac_ice=MAX(0.0,1.-frac_pond-frac_snow)
612	! total albedo
613	alb1_new(i)=alb_snowfrac_snow+alb_icefrac_ice+alb_pond*frac_pond
614	END DO
615	alb2_new(:) = alb1_new(:)
616
617	!****************************************************************************************
618	! 3) Recalculate new ocean temperature (add fluxes below ice)
619	! Melt / freeze from below
620	!*********************************************o***************************************
621	!cumul fluxes
622	bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
623	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
624	! Add cumulated surface fluxes
625	tslab(:,1)=tslab(:,1)+bilg_cum(:)cyangdtime
626	DO i=1,knon
627	ki=knindex(i)
628	! split lateral/top melt-freeze
629	frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin)))
630	IF (tslab(ki,1).LE.t_freeze) THEN
631	! **** Form new ice from below *****
632	! quantity of new ice
633	e_melt=(t_freeze-tslab(ki,1))/cyang &
634	/(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
635	! first increase height to h_thin
636	dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
637	seaice(ki)=dhsic+seaice(ki)
638	e_melt=e_melt-fsic(ki)*dhsic
639	IF (e_melt.GT.0.) THEN
640	! frac_mf fraction used for lateral increase
641	dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
642	fsic(ki)=fsic(ki)+dfsic
643	e_melt=e_melt-dfsic*seaice(ki)
644	! rest used to increase height
645	seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
646	END IF
647	tslab(ki,1)=t_freeze
648	ELSE ! slab temperature above freezing
649	! **** melt ice from below *****
650	! quantity of melted ice
651	e_melt=(tslab(ki,1)-t_freeze)/cyang &
652	/(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
653	! first decrease height to h_thick
654	dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
655	seaice(ki)=seaice(ki)-dhsic
656	e_melt=e_melt-fsic(ki)*dhsic
657	IF (e_melt.GT.0) THEN
658	! frac_mf fraction used for height decrease
659	dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
660	seaice(ki)=seaice(ki)-dhsic
661	e_melt=e_melt-fsic(ki)*dhsic
662	! rest used to decrease fraction (up to 0!)
663	dfsic=MIN(fsic(ki),e_melt/seaice(ki))
664	! keep remaining in ocean
665	e_melt=e_melt-dfsic*seaice(ki)
666	END IF
667	tslab(ki,1)=t_freeze+e_meltice_latcyang
668	fsic(ki)=fsic(ki)-dfsic
669	END IF
670	END DO
671	bilg_cum(:)=0.
672	END IF ! coupling time
673
674	!tests fraction glace
675	WHERE (fsic.LT.ice_frac_min)
676	tslab(:,1)=tslab(:,1)-fsicseaiceice_lat*cyang
677	tice=t_melt
678	fsic=0.
679	seaice=0.
680	END WHERE
681
682	END SUBROUTINE ocean_slab_ice
683	!
684	!****************************************************************************************
685	!
686	SUBROUTINE ocean_slab_final
687
688	!****************************************************************************************
689	! Deallocate module variables
690	!****************************************************************************************
691	IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
692	IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
693	IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
694	IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
695	IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
696	IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
697	IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
698
699	END SUBROUTINE ocean_slab_final
700
701	END MODULE ocean_slab_mod

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

Download in other formats:

Original Format