Context Navigation

ocean_slab_mod.f90 @ 5662

Last change on this file since 5662 was 5662, checked in by Laurent Fairhead, 3 weeks ago
Ajout du modèle thermodynamique de glace de mer interactive améliorant les flux échangés à la surface de la banquise (Doctorat de Nicolas Michalezyk, Contact : Nicolas Michaleyk, Guillaume Gastineau)
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 41.3 KB

Line
1	!Completed
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	USE mod_grid_phy_lmdz, ONLY: klon_glo
12	USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
13
14	IMPLICIT NONE
15	PRIVATE
16	PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
17
18	!***********************************************************************************
19	! Global saved variables
20	!***********************************************************************************
21	! number of slab vertical layers
22	INTEGER, PUBLIC, SAVE :: nslay
23	!$OMP THREADPRIVATE(nslay)
24	! timestep for coupling (update slab temperature) in timesteps
25	INTEGER, PRIVATE, SAVE :: cpl_pas
26	!$OMP THREADPRIVATE(cpl_pas)
27	! cyang = 1/heat capacity of top layer (rho.c.H)
28	REAL, PRIVATE, SAVE :: cyang
29	!$OMP THREADPRIVATE(cyang)
30	! depth of slab layers (1 or 2)
31	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh
32	!$OMP THREADPRIVATE(slabh)
33	! slab temperature
34	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: tslab
35	!$OMP THREADPRIVATE(tslab)
36	! heat flux convergence due to Ekman
37	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_ekman
38	!$OMP THREADPRIVATE(dt_ekman)
39	! heat flux convergence due to horiz diffusion
40	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_hdiff
41	!$OMP THREADPRIVATE(dt_hdiff)
42	! heat flux convergence due to GM eddy advection
43	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_gm
44	!$OMP THREADPRIVATE(dt_gm)
45	! Heat Flux correction
46	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_qflux
47	!$OMP THREADPRIVATE(dt_qflux)
48	! fraction of ocean covered by sea ice (sic / (oce+sic))
49	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic
50	!$OMP THREADPRIVATE(fsic)
51	! temperature of the sea ice
52	!GG
53	!REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice_slab
54	!!$OMP THREADPRIVATE(tice_slab)
55	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice_slab
56	!$OMP THREADPRIVATE(tice_slab)
57	!GG
58	! sea ice thickness, in kg/m2
59	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: seaice
60	!$OMP THREADPRIVATE(seaice)
61	! net surface heat flux, weighted by open ocean fraction
62	! slab_bils accumulated over cpl_pas timesteps
63	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum
64	!$OMP THREADPRIVATE(bils_cum)
65	! net heat flux into the ocean below the ice : conduction + solar radiation
66	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bilg
67	!$OMP THREADPRIVATE(slab_bilg)
68	! slab_bilg over cpl_pas timesteps
69	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bilg_cum
70	!$OMP THREADPRIVATE(bilg_cum)
71	! wind stress saved over cpl_pas timesteps
72	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: taux_cum
73	!$OMP THREADPRIVATE(taux_cum)
74	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tauy_cum
75	!$OMP THREADPRIVATE(tauy_cum)
76
77	!***********************************************************************************
78	! Parameters (could be read in def file: move to slab_init)
79	!***********************************************************************************
80	! snow and ice physical characteristics:
81	REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
82	REAL, PARAMETER :: t_melt=273.15 ! melting ice temp
83	REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
84	REAL, PARAMETER :: ice_den=917. ! ice density
85	REAL, PARAMETER :: sea_den=1025. ! sea water density
86	REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice
87	REAL, PARAMETER :: sno_cond=0.31*sno_den ! conductivity of snow
88	REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice
89	REAL, PARAMETER :: sea_cap=3995. ! specific heat capacity, snow and ice
90	REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
91
92	! control of snow and ice cover & freeze / melt (heights converted to kg/m2)
93	REAL, PARAMETER :: snow_min=0.05*sno_den !critical snow height 5 cm
94	REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean
95	REAL, PARAMETER :: ice_frac_min=0.001
96	REAL, PARAMETER :: ice_frac_max=1. ! less than 1. if min leads fraction
97	REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness
98	REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness
99	! below ice_thin, priority is melt lateral / grow height
100	! ice_thin is also height of new ice
101	REAL, PARAMETER :: h_ice_thick=2.5*ice_den ! thin ice thickness
102	! above ice_thick, priority is melt height / grow lateral
103	REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice
104	REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height
105
106	! albedo and radiation parameters
107	REAL, PARAMETER :: alb_sno_min=0.55 !min snow albedo
108	REAL, PARAMETER :: alb_sno_del=0.3 !max snow albedo = min + del
109	REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
110	REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
111	REAL, PARAMETER :: pen_frac=0.3 !fraction of shortwave penetrating into the
112	! ice (no snow)
113	REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
114
115	! horizontal transport
116	LOGICAL, PUBLIC, SAVE :: slab_hdiff
117	!$OMP THREADPRIVATE(slab_hdiff)
118	LOGICAL, PUBLIC, SAVE :: slab_gm
119	!$OMP THREADPRIVATE(slab_gm)
120	REAL, PRIVATE, SAVE :: coef_hdiff ! coefficient for horizontal diffusion
121	!$OMP THREADPRIVATE(coef_hdiff)
122	INTEGER, PUBLIC, SAVE :: slab_ekman, slab_cadj ! Ekman, conv adjustment
123	!$OMP THREADPRIVATE(slab_ekman)
124	!$OMP THREADPRIVATE(slab_cadj)
125
126	!***********************************************************************************
127
128	CONTAINS
129	!
130	!***********************************************************************************
131	!
132	SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
133	!, seaice_rst etc
134
135	USE ioipsl_getin_p_mod, ONLY : getin_p
136	USE mod_phys_lmdz_transfert_para, ONLY : gather
137	USE slab_heat_transp_mod, ONLY : ini_slab_transp
138
139	! Input variables
140	!***********************************************************************************
141	REAL, INTENT(IN) :: dtime
142	! Variables read from restart file
143	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
144	! surface fractions from start file
145
146	! Local variables
147	!************************************************************************************
148	INTEGER :: error
149	REAL, DIMENSION(klon_glo) :: zmasq_glo
150	CHARACTER (len = 80) :: abort_message
151	CHARACTER (len = 20) :: modname = 'ocean_slab_intit'
152
153	!***********************************************************************************
154	! Define some parameters
155	!***********************************************************************************
156	! Number of slab layers
157	nslay=2
158	CALL getin_p('slab_layers',nslay)
159	print *,'number of slab layers : ',nslay
160	! Layer thickness
161	ALLOCATE(slabh(nslay), stat = error)
162	IF (error /= 0) THEN
163	abort_message='Pb allocation slabh'
164	CALL abort_physic(modname,abort_message,1)
165	ENDIF
166	slabh(1)=50.
167	CALL getin_p('slab_depth',slabh(1))
168	IF (nslay.GT.1) THEN
169	slabh(2)=150.
170	END IF
171
172	! cyang = 1/heat capacity of top layer (rho.c.H)
173	cyang=1/(slabh(1)sea_densea_cap)
174
175	! cpl_pas coupling period (update of tslab and ice fraction)
176	! pour un calcul a chaque pas de temps, cpl_pas=1
177	cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
178	CALL getin_p('cpl_pas',cpl_pas)
179	print *,'cpl_pas',cpl_pas
180
181	! Horizontal diffusion
182	slab_hdiff=.FALSE.
183	CALL getin_p('slab_hdiff',slab_hdiff)
184	coef_hdiff=25000.
185	CALL getin_p('coef_hdiff',coef_hdiff)
186	! Ekman transport
187	slab_ekman=0
188	CALL getin_p('slab_ekman',slab_ekman)
189	! GM eddy advection (2-layers only)
190	slab_gm=.FALSE.
191	CALL getin_p('slab_gm',slab_gm)
192	IF (slab_ekman.LT.2) THEN
193	slab_gm=.FALSE.
194	ENDIF
195	! Convective adjustment
196	IF (nslay.EQ.1) THEN
197	slab_cadj=0
198	ELSE
199	slab_cadj=1
200	END IF
201	CALL getin_p('slab_cadj',slab_cadj)
202
203	!************************************************************************************
204	! Allocate surface fraction read from restart file
205	!************************************************************************************
206	ALLOCATE(fsic(klon), stat = error)
207	IF (error /= 0) THEN
208	abort_message='Pb allocation tmp_pctsrf_slab'
209	CALL abort_physic(modname,abort_message,1)
210	ENDIF
211	fsic(:)=0.
212	!zmasq = continent fraction
213	WHERE (1.-zmasq(:)>EPSFRA)
214	fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
215	END WHERE
216
217	!************************************************************************************
218	! Allocate saved fields
219	!************************************************************************************
220	ALLOCATE(tslab(klon,nslay), stat=error)
221	IF (error /= 0) CALL abort_physic &
222	(modname,'pb allocation tslab', 1)
223
224	ALLOCATE(bils_cum(klon), stat = error)
225	IF (error /= 0) THEN
226	abort_message='Pb allocation slab_bils_cum'
227	CALL abort_physic(modname,abort_message,1)
228	ENDIF
229	bils_cum(:) = 0.0
230
231	IF (version_ocean=='sicINT') THEN ! interactive sea ice
232	ALLOCATE(slab_bilg(klon), stat = error)
233	IF (error /= 0) THEN
234	abort_message='Pb allocation slab_bilg'
235	CALL abort_physic(modname,abort_message,1)
236	ENDIF
237	slab_bilg(:) = 0.0
238	ALLOCATE(bilg_cum(klon), stat = error)
239	IF (error /= 0) THEN
240	abort_message='Pb allocation slab_bilg_cum'
241	CALL abort_physic(modname,abort_message,1)
242	ENDIF
243	bilg_cum(:) = 0.0
244	!GG ALLOCATE(tice(klon), stat = error)
245	ALLOCATE(tice_slab(klon), stat = error)
246	IF (error /= 0) THEN
247	abort_message='Pb allocation slab_tice'
248	CALL abort_physic(modname,abort_message,1)
249	ENDIF
250	ALLOCATE(seaice(klon), stat = error)
251	IF (error /= 0) THEN
252	abort_message='Pb allocation slab_seaice'
253	CALL abort_physic(modname,abort_message,1)
254	ENDIF
255	END IF
256
257	IF (slab_hdiff) THEN !horizontal diffusion
258	ALLOCATE(dt_hdiff(klon,nslay), stat = error)
259	IF (error /= 0) THEN
260	abort_message='Pb allocation dt_hdiff'
261	CALL abort_physic(modname,abort_message,1)
262	ENDIF
263	dt_hdiff(:,:) = 0.0
264	ENDIF
265
266	ALLOCATE(dt_qflux(klon,nslay), stat = error)
267	IF (error /= 0) THEN
268	abort_message='Pb allocation dt_qflux'
269	CALL abort_physic(modname,abort_message,1)
270	ENDIF
271	dt_qflux(:,:) = 0.0
272
273	IF (slab_gm) THEN !GM advection
274	ALLOCATE(dt_gm(klon,nslay), stat = error)
275	IF (error /= 0) THEN
276	abort_message='Pb allocation dt_gm'
277	CALL abort_physic(modname,abort_message,1)
278	ENDIF
279	dt_gm(:,:) = 0.0
280	ENDIF
281
282	IF (slab_ekman.GT.0) THEN ! ekman transport
283	ALLOCATE(dt_ekman(klon,nslay), stat = error)
284	IF (error /= 0) THEN
285	abort_message='Pb allocation dt_ekman'
286	CALL abort_physic(modname,abort_message,1)
287	ENDIF
288	dt_ekman(:,:) = 0.0
289	ALLOCATE(taux_cum(klon), stat = error)
290	IF (error /= 0) THEN
291	abort_message='Pb allocation taux_cum'
292	CALL abort_physic(modname,abort_message,1)
293	ENDIF
294	taux_cum(:) = 0.0
295	ALLOCATE(tauy_cum(klon), stat = error)
296	IF (error /= 0) THEN
297	abort_message='Pb allocation tauy_cum'
298	CALL abort_physic(modname,abort_message,1)
299	ENDIF
300	tauy_cum(:) = 0.0
301	ENDIF
302
303	! Initialize transport
304	IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
305	CALL gather(zmasq,zmasq_glo)
306	! Master thread/process only
307	!$OMP MASTER
308	IF (is_mpi_root) THEN
309	CALL ini_slab_transp(zmasq_glo)
310	END IF
311	!$OMP END MASTER
312	END IF
313
314	END SUBROUTINE ocean_slab_init
315	!
316	!***********************************************************************************
317	!
318	SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
319
320	! this routine sends back the sea ice and ocean fraction to the main physics
321	! routine. Called only with interactive sea ice
322
323	! Arguments
324	!************************************************************************************
325	INTEGER, INTENT(IN) :: itime ! current timestep
326	INTEGER, INTENT(IN) :: jour ! day in year (not
327	REAL , INTENT(IN) :: dtime ! physics timestep (s)
328	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg ! sub-surface fraction
329	LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is
330	! modified at this time step
331
332	pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
333	pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
334	is_modified=.TRUE.
335
336	END SUBROUTINE ocean_slab_frac
337	!
338	!************************************************************************************
339	!
340	SUBROUTINE ocean_slab_noice( &
341	itime, dtime, jour, knon, knindex, &
342	p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum, &
343	AcoefH, AcoefQ, BcoefH, BcoefQ, &
344	AcoefU, AcoefV, BcoefU, BcoefV, &
345	ps, u1, v1, gustiness, tsurf_in, &
346	radsol, snow, &
347	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
348	tsurf_new, dflux_s, dflux_l, slab_bils)
349
350	USE clesphys_mod_h
351	USE calcul_fluxs_mod
352	USE slab_heat_transp_mod, ONLY: divgrad_phy,slab_ekman1,slab_ekman2,slab_gmdiff
353	USE mod_phys_lmdz_para
354
355
356	! This routine
357	! (1) computes surface turbulent fluxes over points with some open ocean
358	! (2) reads additional Q-flux (everywhere)
359	! (3) computes horizontal transport (diffusion & Ekman)
360	! (4) updates slab temperature every cpl_pas ; creates new ice if needed.
361
362	! Note :
363	! klon total number of points
364	! knon number of points with open ocean (varies with time)
365	! knindex gives position of the knon points within klon.
366	! In general, local saved variables have klon values
367	! variables exchanged with PBL module have knon.
368
369	! Input arguments
370	!***********************************************************************************
371	INTEGER, INTENT(IN) :: itime ! current timestep INTEGER,
372	INTEGER, INTENT(IN) :: jour ! day in year (for Q-Flux)
373	INTEGER, INTENT(IN) :: knon ! number of points
374	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
375	REAL, INTENT(IN) :: dtime ! timestep (s)
376	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
377	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm
378	! drag coefficients
379	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
380	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum ! near surface T, q
381	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
382	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
383	! exchange coefficients for boundary layer scheme
384	REAL, DIMENSION(klon), INTENT(IN) :: ps ! surface pressure
385	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness ! surface wind
386	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! surface temperature
387	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! net surface radiative flux
388
389	! In/Output arguments
390	!************************************************************************************
391	REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2
392
393	! Output arguments
394	!************************************************************************************
395	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
396	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
397	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
398	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! new surface tempearture
399	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
400	REAL, DIMENSION(klon), INTENT(OUT) :: slab_bils
401
402	! Local variables
403	!************************************************************************************
404	INTEGER :: i,ki,k
405	REAL :: t_cadj
406	! for surface heat fluxes
407	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
408	! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes
409	REAL, DIMENSION(klon) :: diff_sst, diff_siv
410	REAL, DIMENSION(klon,nslay) :: lmt_bils
411	! for surface wind stress
412	REAL, DIMENSION(klon) :: u0, v0
413	REAL, DIMENSION(klon) :: u1_lay, v1_lay
414	! for new ice creation
415	REAL :: e_freeze, h_new, dfsic
416	! horizontal diffusion and Ekman local vars
417	! dimension = global domain (klon_glo) instead of // subdomains
418	REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo,dt_ekman_glo,dt_gm_glo
419	! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
420	REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
421	REAL, DIMENSION(klon_glo,nslay) :: tslab_glo
422	REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo
423
424	!****************************************************************************************
425	! 1) Surface fluxes calculation
426	!
427	!****************************************************************************************
428	!cal(:) = 0. ! infinite thermal inertia
429	!beta(:) = 1. ! wet surface
430	!dif_grnd(:) = 0. ! no diffusion into ground
431	! EV: use calbeta
432	CALL calbeta(dtime, is_oce, knon, snow,qsurf, beta, cal, dif_grnd)
433
434
435
436	! Suppose zero surface speed
437	u0(:)=0.0
438	v0(:)=0.0
439	u1_lay(:) = u1(:) - u0(:)
440	v1_lay(:) = v1(:) - v0(:)
441
442	! Compute latent & sensible fluxes
443	CALL calcul_fluxs(knon, is_oce, dtime, &
444	tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, &
445	precip_rain, precip_snow, snow, qsurf, &
446	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
447	f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
448	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
449
450	! save total cumulated heat fluxes locally
451	! radiative + turbulent + melt of falling snow
452	slab_bils(:)=0.
453	DO i=1,knon
454	ki=knindex(i)
455	slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
456	-precip_snow(i)ice_lat(1.+snow_wfact*fsic(ki)))
457	bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
458	END DO
459
460	! Compute surface wind stress
461	CALL calcul_flux_wind(knon, dtime, &
462	u0, v0, u1, v1, gustiness, cdragm, &
463	AcoefU, AcoefV, BcoefU, BcoefV, &
464	p1lay, temp_air, &
465	flux_u1, flux_v1)
466
467	! save cumulated wind stress
468	IF (slab_ekman.GT.0) THEN
469	DO i=1,knon
470	ki=knindex(i)
471	taux_cum(ki)=taux_cum(ki)+flux_u1(i)*(1.-fsic(ki))/cpl_pas
472	tauy_cum(ki)=tauy_cum(ki)+flux_v1(i)*(1.-fsic(ki))/cpl_pas
473	END DO
474	ENDIF
475
476	!****************************************************************************************
477	! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file limit_slab.nc
478	!
479	!****************************************************************************************
480	CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
481	! lmt_bils and diff_sst,siv saved by limit_slab
482	! qflux = total QFlux correction (in W/m2)
483	dt_qflux(:,1)=lmt_bils(:,1)+diff_sst(:)/cyang/86400.-diff_siv(:)ice_denice_lat/86400.
484	IF (nslay.GT.1) THEN
485	dt_qflux(:,2:nslay)=lmt_bils(:,2:nslay)
486	END IF
487
488	!****************************************************************************************
489	! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport)
490	! Bring to freezing temp and make sea ice if necessary
491	!
492	!*********************************************o***************************************
493	tsurf_new=tsurf_in
494	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
495	! ***********************************
496	! Horizontal transport
497	! ***********************************
498	IF (slab_ekman.GT.0) THEN
499	! copy wind stress to global var
500	CALL gather(taux_cum,taux_glo)
501	CALL gather(tauy_cum,tauy_glo)
502	END IF
503
504	IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
505	CALL gather(tslab,tslab_glo)
506	! Compute horiz transport on one process only
507	IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus
508	IF (slab_hdiff) THEN
509	dt_hdiff_glo(:,:)=0.
510	END IF
511	IF (slab_ekman.GT.0) THEN
512	dt_ekman_glo(:,:)=0.
513	END IF
514	IF (slab_gm) THEN
515	dt_gm_glo(:,:)=0.
516	END IF
517	DO i=1,cpl_pas ! time splitting for numerical stability
518	IF (slab_ekman.GT.0) THEN
519	SELECT CASE (slab_ekman)
520	CASE (1)
521	CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
522	CASE (2)
523	CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm)
524	CASE DEFAULT
525	dt_ekman_tmp(:,:)=0.
526	END SELECT
527	dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:)
528	! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1
529	DO k=1,nslay
530	dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den)
531	ENDDO
532	tslab_glo=tslab_glo+dt_ekman_tmp*dtime
533	IF (slab_gm) THEN ! Gent-McWilliams eddy advection
534	dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:)
535	! convert dt from K.s-1.(kg.m-2) to K.s-1
536	DO k=1,nslay
537	dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/(slabh(k)*sea_den)
538	END DO
539	tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
540	END IF
541	ENDIF
542	! GM included in Ekman_2
543	! IF (slab_gm) THEN ! Gent-McWilliams eddy advection
544	! CALL slab_gmdiff(tslab_glo,dt_hdiff_tmp)
545	! ! convert dt_gm from K.m.s-1 to K.s-1
546	! DO k=1,nslay
547	! dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/slabh(k)
548	! END DO
549	! tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
550	! dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:)
551	! END IF
552	IF (slab_hdiff) THEN ! horizontal diffusion
553	! laplacian of slab T
554	CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp)
555	! multiply by diff coef and normalize to 50m slab equivalent
556	dt_hdiff_tmp=dt_hdiff_tmpcoef_hdiff50./SUM(slabh)
557	dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:)
558	tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
559	END IF
560	END DO ! time splitting
561	IF (slab_hdiff) THEN
562	!dt_hdiff_glo saved in W/m2
563	DO k=1,nslay
564	dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)slabh(k)sea_den*sea_cap/cpl_pas
565	END DO
566	END IF
567	IF (slab_gm) THEN
568	!dt_hdiff_glo saved in W/m2
569	dt_gm_glo(:,:)=dt_gm_glo(:,:)*sea_cap/cpl_pas
570	END IF
571	IF (slab_ekman.GT.0) THEN
572	! dt_ekman_glo saved in W/m2
573	dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
574	END IF
575	END IF ! master process
576	!$OMP BARRIER
577	! Send new fields back to all processes
578	CALL Scatter(tslab_glo,tslab)
579	IF (slab_hdiff) THEN
580	CALL Scatter(dt_hdiff_glo,dt_hdiff)
581	END IF
582	IF (slab_gm) THEN
583	CALL Scatter(dt_gm_glo,dt_gm)
584	END IF
585	IF (slab_ekman.GT.0) THEN
586	CALL Scatter(dt_ekman_glo,dt_ekman)
587	! clear wind stress
588	taux_cum(:)=0.
589	tauy_cum(:)=0.
590	END IF
591	ENDIF ! transport
592
593	! ***********************************
594	! Other heat fluxes
595	! ***********************************
596	! Add read QFlux
597	DO k=1,nslay
598	tslab(:,k)=tslab(:,k)+dt_qflux(:,k)cyangdtime*cpl_pas &
599	*slabh(1)/slabh(k)
600	END DO
601	! Add cumulated surface fluxes
602	tslab(:,1)=tslab(:,1)+bils_cum(:)cyangdtime
603	! Convective adjustment if 2 layers
604	IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN
605	DO i=1,klon
606	IF (tslab(i,2).GT.tslab(i,1)) THEN
607	! mean (mass-weighted) temperature
608	t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:))
609	tslab(i,1)=t_cadj
610	tslab(i,2)=t_cadj
611	END IF
612	END DO
613	END IF
614	! ***********************************
615	! Update surface temperature and ice
616	! ***********************************
617	SELECT CASE(version_ocean)
618	CASE('sicNO') ! no sea ice even below freezing !
619	DO i=1,knon
620	ki=knindex(i)
621	tsurf_new(i)=tslab(ki,1)
622	END DO
623	CASE('sicOBS') ! "realistic" case, for prescribed sea ice
624	! tslab cannot be below freezing, or above it if there is sea ice
625	DO i=1,knon
626	ki=knindex(i)
627	IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
628	tslab(ki,1)=t_freeze
629	END IF
630	tsurf_new(i)=tslab(ki,1)
631	END DO
632	CASE('sicINT') ! interactive sea ice
633	DO i=1,knon
634	ki=knindex(i)
635	IF (fsic(ki).LT.epsfra) THEN ! Free of ice
636	IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
637	! quantity of new ice formed
638	e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
639	! new ice
640	!GG tice(ki)=t_freeze
641	tice_slab(ki)=t_freeze
642	fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
643	IF (fsic(ki).GT.ice_frac_min) THEN
644	seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
645	tslab(ki,1)=t_freeze
646	ELSE
647	fsic(ki)=0.
648	END IF
649	tsurf_new(i)=t_freeze
650	ELSE
651	tsurf_new(i)=tslab(ki,1)
652	END IF
653	ELSE ! ice present
654	tsurf_new(i)=t_freeze
655	IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
656	! quantity of new ice formed over open ocean
657	e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
658	/(ice_lat+ice_cap/2.*(t_freeze-tice_slab(ki)))
659	!GG /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
660	! new ice height and fraction
661	h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
662	dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
663	h_new=MIN(e_freeze/dfsic,h_ice_max)
664	! update tslab to freezing over open ocean only
665	tslab(ki,1)=tslab(ki,1)fsic(ki)+t_freeze(1.-fsic(ki))
666	! update sea ice
667	seaice(ki)=(h_newdfsic+seaice(ki)fsic(ki)) &
668	/(dfsic+fsic(ki))
669	fsic(ki)=fsic(ki)+dfsic
670	! update snow?
671	END IF ! tslab below freezing
672	END IF ! sea ice present
673	END DO
674	END SELECT
675	bils_cum(:)=0.0! clear cumulated fluxes
676	END IF ! coupling time
677	END SUBROUTINE ocean_slab_noice
678	!
679	!****************************************************************************************
680
681	SUBROUTINE ocean_slab_ice( &
682	itime, dtime, jour, knon, knindex, &
683	tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
684	AcoefH, AcoefQ, BcoefH, BcoefQ, &
685	AcoefU, AcoefV, BcoefU, BcoefV, &
686	ps, u1, v1, gustiness, &
687	radsol, snow, qsurf, qsol, agesno, &
688	alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
689	tsurf_new, dflux_s, dflux_l, swnet)
690
691	USE clesphys_mod_h
692	USE yomcst_mod_h
693	USE calcul_fluxs_mod
694
695
696
697	! Input arguments
698	!****************************************************************************************
699	INTEGER, INTENT(IN) :: itime, jour, knon
700	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
701	REAL, INTENT(IN) :: dtime
702	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
703	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
704	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
705	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
706	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
707	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
708	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
709	REAL, DIMENSION(klon), INTENT(IN) :: ps
710	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
711	REAL, DIMENSION(klon), INTENT(IN) :: swnet
712
713	! In/Output arguments
714	!****************************************************************************************
715	REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
716	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
717	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
718
719	! Output arguments
720	!****************************************************************************************
721	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
722	REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval
723	REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval
724	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
725	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
726	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
727	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
728
729	! Local variables
730	!****************************************************************************************
731	INTEGER :: i,ki
732	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
733	REAL, DIMENSION(klon) :: u0, v0
734	REAL, DIMENSION(klon) :: u1_lay, v1_lay
735	! intermediate heat fluxes:
736	REAL :: f_cond, f_swpen
737	! for snow/ice albedo:
738	REAL :: alb_snow, alb_ice, alb_pond
739	REAL :: frac_snow, frac_ice, frac_pond
740	! for ice melt / freeze
741	REAL :: e_melt, snow_evap, h_test
742	! dhsic, dfsic change in ice mass, fraction.
743	REAL :: dhsic, dfsic, frac_mf
744
745	!****************************************************************************************
746	! 1) Flux calculation
747	!****************************************************************************************
748	! Suppose zero surface speed
749	u0(:)=0.0
750	v0(:)=0.0
751	u1_lay(:) = u1(:) - u0(:)
752	v1_lay(:) = v1(:) - v0(:)
753
754	! set beta, cal, compute conduction fluxes inside ice/snow
755	slab_bilg(:)=0.
756	!dif_grnd(:)=0.
757	!beta(:) = 1.
758	! EV: use calbeta to calculate beta and then recalculate properly cal
759	CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, cal, dif_grnd)
760
761
762	DO i=1,knon
763	ki=knindex(i)
764	IF (snow(i).GT.snow_min) THEN
765	! snow-layer heat capacity
766	cal(i)=2.RCPD/(snow(i)ice_cap)
767	! snow conductive flux
768	f_cond=sno_cond*(tice_slab(ki)-tsurf_in(i))/snow(i)
769	!GG f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
770	! all shortwave flux absorbed
771	f_swpen=0.
772	! bottom flux (ice conduction)
773	slab_bilg(ki)=ice_cond*(tice_slab(ki)-t_freeze)/seaice(ki)
774	!GG slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
775	! update ice temperature
776	!GG tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
777	tice_slab(ki)=tice_slab(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
778	(slab_bilg(ki)+f_cond)dtime
779	ELSE ! bare ice
780	! ice-layer heat capacity
781	cal(i)=2.RCPD/(seaice(ki)ice_cap)
782	! conductive flux
783	!GG f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
784	f_cond=ice_cond*(t_freeze-tice_slab(ki))/seaice(ki)
785	! penetrative shortwave flux...
786	f_swpen=swnet(i)pen_fracexp(-pen_ext*seaice(ki)/ice_den)
787	slab_bilg(ki)=f_swpen-f_cond
788	END IF
789	radsol(i)=radsol(i)+f_cond-f_swpen
790	END DO
791	! weight fluxes to ocean by sea ice fraction
792	slab_bilg(:)=slab_bilg(:)*fsic(:)
793
794	! calcul_fluxs (sens, lat etc)
795	CALL calcul_fluxs(knon, is_sic, dtime, &
796	tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, &
797	precip_rain, precip_snow, snow, qsurf, &
798	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
799	f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
800	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
801	DO i=1,knon
802	!GG IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
803	IF (snow(i).LT.snow_min) tice_slab(knindex(i))=tsurf_new(i)
804	END DO
805
806	! calcul_flux_wind
807	CALL calcul_flux_wind(knon, dtime, &
808	u0, v0, u1, v1, gustiness, cdragm, &
809	AcoefU, AcoefV, BcoefU, BcoefV, &
810	p1lay, temp_air, &
811	flux_u1, flux_v1)
812
813	!****************************************************************************************
814	! 2) Update snow and ice surface
815	!****************************************************************************************
816	! snow precip
817	DO i=1,knon
818	ki=knindex(i)
819	IF (precip_snow(i) > 0.) THEN
820	snow(i) = snow(i)+precip_snow(i)dtime(1.-snow_wfact*(1.-fsic(ki)))
821	END IF
822	! snow and ice sublimation
823	IF (evap(i) > 0.) THEN
824	snow_evap = MIN (snow(i) / dtime, evap(i))
825	snow(i) = snow(i) - snow_evap * dtime
826	snow(i) = MAX(0.0, snow(i))
827	seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
828	ENDIF
829	! Melt / Freeze snow from above if Tsurf>0
830	IF (tsurf_new(i).GT.t_melt) THEN
831	! energy available for melting snow (in kg of melted snow /m2)
832	e_melt = MIN(MAX(snow(i)(tsurf_new(i)-t_melt)ice_cap/2. &
833	/(ice_lat+ice_cap/2.*(t_melt-tice_slab(ki))),0.0),snow(i))
834	!GG /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
835	! remove snow
836	IF (snow(i).GT.e_melt) THEN
837	snow(i)=snow(i)-e_melt
838	tsurf_new(i)=t_melt
839	ELSE ! all snow is melted
840	! add remaining heat flux to ice
841	e_melt=e_melt-snow(i)
842	tice_slab(ki)=tice_slab(ki)+e_meltice_lat2./(ice_cap*seaice(ki))
843	!GG tice(ki)=tice(ki)+e_meltice_lat2./(ice_cap*seaice(ki))
844	tsurf_new(i)=tice_slab(ki)
845	!GG tsurf_new(i)=tice(ki)
846	END IF
847	END IF
848	! melt ice from above if Tice>0
849	!GG IF (tice(ki).GT.t_melt) THEN
850	IF (tice_slab(ki).GT.t_melt) THEN
851	! quantity of ice melted (kg/m2)
852	!GG e_melt=MAX(seaice(ki)(tice(ki)-t_melt)ice_cap/2. &
853	e_melt=MAX(seaice(ki)(tice_slab(ki)-t_melt)ice_cap/2. &
854	/(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
855	! melt from above, height only
856	dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
857	e_melt=e_melt-dhsic
858	IF (e_melt.GT.0) THEN
859	! lateral melt if ice too thin
860	dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
861	! if all melted add remaining heat to ocean
862	e_melt=MAX(0.,e_meltfsic(ki)-dfsich_ice_min)
863	slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime
864	! update height and fraction
865	fsic(ki)=fsic(ki)-dfsic
866	END IF
867	seaice(ki)=seaice(ki)-dhsic
868	! surface temperature at melting point
869	!GG tice(ki)=t_melt
870	tice_slab(ki)=t_melt
871	tsurf_new(i)=t_melt
872	END IF
873	! convert snow to ice if below floating line
874	h_test=(seaice(ki)+snow(i))ice_den-seaice(ki)sea_den
875	IF (h_test.GT.0.) THEN !snow under water
876	! extra snow converted to ice (with added frozen sea water)
877	dhsic=h_test/(sea_den-ice_den+sno_den)
878	seaice(ki)=seaice(ki)+dhsic
879	snow(i)=snow(i)-dhsic*sno_den/ice_den
880	! available energy (freeze sea water + bring to tice_slab)
881	e_melt=dhsic(1.-sno_den/ice_den)(ice_lat+ &
882	ice_cap/2.*(t_freeze-tice_slab(ki)))
883	!GG ice_cap/2.*(t_freeze-tice(ki)))
884	! update ice temperature
885	!GG tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
886	tice_slab(ki)=tice_slab(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
887	END IF
888	END DO
889
890	! New albedo
891	DO i=1,knon
892	ki=knindex(i)
893	! snow albedo: update snow age
894	IF (snow(i).GT.0.0001) THEN
895	agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
896	* EXP(-1.MAX(0.0,precip_snow(i))dtime/5.)
897	ELSE
898	agesno(i)=0.0
899	END IF
900	! snow albedo
901	alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.)
902	! ice albedo (varies with ice tkickness and temp)
903	alb_ice=MAX(0.0,0.13LOG(100.seaice(ki)/ice_den)+0.1)
904	!GG IF (tice(ki).GT.t_freeze-0.01) THEN
905	IF (tice_slab(ki).GT.t_freeze-0.01) THEN
906	alb_ice=MIN(alb_ice,alb_ice_wet)
907	ELSE
908	alb_ice=MIN(alb_ice,alb_ice_dry)
909	END IF
910	! pond albedo
911	alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice_slab(ki)-t_melt,-2.0)))
912	!GG alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
913	! pond fraction
914	frac_pond=0.2*(2.0+MIN(0.0,MAX(tice_slab(ki)-t_melt,-2.0)))
915	!GG frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
916	! snow fraction
917	frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
918	! ice fraction
919	frac_ice=MAX(0.0,1.-frac_pond-frac_snow)
920	! total albedo
921	alb1_new(i)=alb_snowfrac_snow+alb_icefrac_ice+alb_pond*frac_pond
922	END DO
923	alb2_new(:) = alb1_new(:)
924
925	!****************************************************************************************
926	! 3) Recalculate new ocean temperature (add fluxes below ice)
927	! Melt / freeze from below
928	!*********************************************o***************************************
929	!cumul fluxes
930	bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
931	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
932	! Add cumulated surface fluxes
933	tslab(:,1)=tslab(:,1)+bilg_cum(:)cyangdtime
934	DO i=1,knon
935	ki=knindex(i)
936	! split lateral/top melt-freeze
937	frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin)))
938	IF (tslab(ki,1).LE.t_freeze) THEN
939	! **** Form new ice from below *****
940	! quantity of new ice
941	e_melt=(t_freeze-tslab(ki,1))/cyang &
942	/(ice_lat+ice_cap/2.*(t_freeze-tice_slab(ki)))
943	!GG /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
944	! first increase height to h_thin
945	dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
946	seaice(ki)=dhsic+seaice(ki)
947	e_melt=e_melt-fsic(ki)*dhsic
948	IF (e_melt.GT.0.) THEN
949	! frac_mf fraction used for lateral increase
950	dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
951	fsic(ki)=fsic(ki)+dfsic
952	e_melt=e_melt-dfsic*seaice(ki)
953	! rest used to increase height
954	seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
955	END IF
956	tslab(ki,1)=t_freeze
957	ELSE ! slab temperature above freezing
958	! **** melt ice from below *****
959	! quantity of melted ice
960	e_melt=(tslab(ki,1)-t_freeze)/cyang &
961	/(ice_lat+ice_cap/2.*(tice_slab(ki)-t_freeze))
962	!GG /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
963	! first decrease height to h_thick
964	dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
965	seaice(ki)=seaice(ki)-dhsic
966	e_melt=e_melt-fsic(ki)*dhsic
967	IF (e_melt.GT.0) THEN
968	! frac_mf fraction used for height decrease
969	dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
970	seaice(ki)=seaice(ki)-dhsic
971	e_melt=e_melt-fsic(ki)*dhsic
972	! rest used to decrease fraction (up to 0!)
973	dfsic=MIN(fsic(ki),e_melt/seaice(ki))
974	! keep remaining in ocean
975	e_melt=e_melt-dfsic*seaice(ki)
976	END IF
977	tslab(ki,1)=t_freeze+e_meltice_latcyang
978	fsic(ki)=fsic(ki)-dfsic
979	END IF
980	END DO
981	bilg_cum(:)=0.
982	END IF ! coupling time
983
984	!tests ice fraction
985	WHERE (fsic.LT.ice_frac_min)
986	tslab(:,1)=tslab(:,1)-fsicseaiceice_lat*cyang
987	tice_slab=t_melt
988	!GG tice=t_melt
989	fsic=0.
990	seaice=0.
991	END WHERE
992
993	END SUBROUTINE ocean_slab_ice
994	!
995	!****************************************************************************************
996	!
997	SUBROUTINE ocean_slab_final
998
999	!****************************************************************************************
1000	! Deallocate module variables
1001	!****************************************************************************************
1002	IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
1003	IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
1004	IF (ALLOCATED(tice_slab)) DEALLOCATE(tice_slab)
1005	IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
1006	IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
1007	IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
1008	IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
1009	IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum)
1010	IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum)
1011	IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
1012	IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
1013	IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm)
1014	IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux)
1015
1016	END SUBROUTINE ocean_slab_final
1017
1018	END MODULE ocean_slab_mod

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

Download in other formats:

Original Format