Context Navigation

ocean_slab_mod.F90 @ 5360

Last change on this file since 5360 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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	! 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, cdragq, 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	INCLUDE "clesphys.h"
228
229	! Input arguments
230	!****************************************************************************************
231	INTEGER, INTENT(IN) :: itime
232	INTEGER, INTENT(IN) :: jour
233	INTEGER, INTENT(IN) :: knon
234	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
235	REAL, INTENT(IN) :: dtime
236	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
237	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm
238	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
239	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
240	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
241	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
242	REAL, DIMENSION(klon), INTENT(IN) :: ps
243	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
244	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
245	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
246
247	! In/Output arguments
248	!****************************************************************************************
249	REAL, DIMENSION(klon), INTENT(INOUT) :: snow
250
251	! Output arguments
252	!****************************************************************************************
253	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
254	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
255	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
256	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
257	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
258	REAL, DIMENSION(klon), INTENT(OUT) :: qflux
259
260	! Local variables
261	!****************************************************************************************
262	INTEGER :: i,ki
263	! surface fluxes
264	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
265	! flux correction
266	REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
267	! surface wind stress
268	REAL, DIMENSION(klon) :: u0, v0
269	REAL, DIMENSION(klon) :: u1_lay, v1_lay
270	! ice creation
271	REAL :: e_freeze, h_new, dfsic
272
273	!****************************************************************************************
274	! 1) Flux calculation
275	!
276	!****************************************************************************************
277	cal(:) = 0.
278	beta(:) = 1.
279	dif_grnd(:) = 0.
280
281	! Suppose zero surface speed
282	u0(:)=0.0
283	v0(:)=0.0
284	u1_lay(:) = u1(:) - u0(:)
285	v1_lay(:) = v1(:) - v0(:)
286
287	CALL calcul_fluxs(knon, is_oce, dtime, &
288	tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, &
289	precip_rain, precip_snow, snow, qsurf, &
290	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
291	f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
292	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
293
294	! - Flux calculation at first modele level for U and V
295	CALL calcul_flux_wind(knon, dtime, &
296	u0, v0, u1, v1, gustiness, cdragm, &
297	AcoefU, AcoefV, BcoefU, BcoefV, &
298	p1lay, temp_air, &
299	flux_u1, flux_v1)
300
301	! Accumulate total fluxes locally
302	slab_bils(:)=0.
303	DO i=1,knon
304	ki=knindex(i)
305	slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
306	-precip_snow(i)ice_lat(1.+snow_wfact*fsic(ki)))
307	bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
308	! Also taux, tauy, saved vars...
309	END DO
310
311	!****************************************************************************************
312	! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
313	!
314	!****************************************************************************************
315	lmt_bils(:)=0.
316	CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
317	! lmt_bils and diff_sst,siv saved by limit_slab
318	qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)ice_denice_lat/86400.
319	! qflux = total QFlux correction (in W/m2)
320
321	!****************************************************************************************
322	! 3) Recalculate new temperature
323	!
324	!*********************************************o***************************************
325	tsurf_new=tsurf_in
326	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
327	! Compute transport
328	! Add read QFlux and SST tendency
329	tslab(:,1)=tslab(:,1)+qflux(:)cyangdtime*cpl_pas
330	! Add cumulated surface fluxes
331	tslab(:,1)=tslab(:,1)+bils_cum(:)cyangdtime
332	! Update surface temperature
333	SELECT CASE(version_ocean)
334	CASE('sicNO')
335	DO i=1,knon
336	ki=knindex(i)
337	tsurf_new(i)=tslab(ki,1)
338	END DO
339	CASE('sicOBS') ! check for sea ice or tslab below freezing
340	DO i=1,knon
341	ki=knindex(i)
342	IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
343	tslab(ki,1)=t_freeze
344	END IF
345	tsurf_new(i)=tslab(ki,1)
346	END DO
347	CASE('sicINT')
348	DO i=1,knon
349	ki=knindex(i)
350	IF (fsic(ki).LT.epsfra) THEN ! Free of ice
351	IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
352	! quantity of new ice formed
353	e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
354	! new ice
355	tice(ki)=t_freeze
356	fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
357	IF (fsic(ki).GT.ice_frac_min) THEN
358	seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
359	tslab(ki,1)=t_freeze
360	ELSE
361	fsic(ki)=0.
362	END IF
363	tsurf_new(i)=t_freeze
364	ELSE
365	tsurf_new(i)=tslab(ki,1)
366	END IF
367	ELSE ! ice present
368	tsurf_new(i)=t_freeze
369	IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
370	! quantity of new ice formed over open ocean
371	e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
372	/(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
373	! new ice height and fraction
374	h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
375	dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
376	h_new=MIN(e_freeze/dfsic,h_ice_max)
377	! update tslab to freezing over open ocean only
378	tslab(ki,1)=tslab(ki,1)fsic(ki)+t_freeze(1.-fsic(ki))
379	! update sea ice
380	seaice(ki)=(h_newdfsic+seaice(ki)fsic(ki)) &
381	/(dfsic+fsic(ki))
382	fsic(ki)=fsic(ki)+dfsic
383	! update snow?
384	END IF !freezing
385	END IF ! sea ice present
386	END DO
387	END SELECT
388	bils_cum(:)=0.0! clear cumulated fluxes
389	END IF ! coupling time
390	END SUBROUTINE ocean_slab_noice
391	!
392	!****************************************************************************************
393
394	SUBROUTINE ocean_slab_ice( &
395	itime, dtime, jour, knon, knindex, &
396	tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
397	AcoefH, AcoefQ, BcoefH, BcoefQ, &
398	AcoefU, AcoefV, BcoefU, BcoefV, &
399	ps, u1, v1, gustiness, &
400	radsol, snow, qsurf, qsol, agesno, &
401	alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
402	tsurf_new, dflux_s, dflux_l, swnet)
403
404	USE calcul_fluxs_mod
405
406	INCLUDE "YOMCST.h"
407	INCLUDE "clesphys.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, cdragh, ps, &
501	precip_rain, precip_snow, snow, qsurf, &
502	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
503	f_qsat_oce,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