Context Navigation

ocean_slab_mod.F90 @ 3300

Last change on this file since 3300 was 3291, checked in by mturbet, 11 months ago
improved parameterization of sea ice albedo vs thickness
Property svn:executable set to ``*
File size: 46.3 KB

Line
1	!Completed
2	MODULE ocean_slab_mod
3	!
4	!==================================================================
5	!
6	! Purpose
7	! -------
8	! The dynamical slab ocean model of the Generic-PCM. It has the following features:
9	! (a) Computes sea ice creation and evolution.
10	! (b) Snow has thermodynamic properties.
11	! (c) Computes oceanic horizontal transport (diffusion & surface-wind driven Ekman transport).
12	! (d) Can be used in parallel mode.
13	!
14	! Authors
15	! -------
16	! S. Bhatnagar and E. Millour (2023)
17	! Adapted from the ocean modules of LMDZ Earth (F. Codron) and the Generic-PCM (B. Charnay, 2013).
18	!
19	! Notes
20	! -----
21	! Compared to the old model, the new model has the following changes (non-exhaustive):
22	! (a) More realistic description of sea ice creation and evolution - simultaneous
23	! surface, side and bottom melting / freezing depending on fluxes.
24	! (b) Snow has an effective heat capacity.
25	! (c) Snow has "weight"; it can sink an ice block if there is too much of it.
26	! (d) Snow can be blown off by wind.
27	! (e) The two-layer ocean allows for convective adjustment.
28	! (f) Diffusion can follow the Gent-McWilliams scheme + Eddy diffusivity.
29	! (g) Can be used in parallel mode.
30	!
31	!==================================================================
32
33	USE dimphy, ONLY: klon
34	USE mod_grid_phy_lmdz, ONLY: klon_glo
35	USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
36
37	IMPLICIT NONE
38	PRIVATE
39	PUBLIC :: ocean_slab_init, ocean_slab_ice, ocean_slab_noice, &
40	ocean_slab_frac, ocean_slab_get_vars, ocean_slab_final
41
42	!***********************************************************************************
43	! Global saved variables
44	!***********************************************************************************
45	! number of slab vertical layers
46	INTEGER, PUBLIC, SAVE :: nslay=2
47	!$OMP THREADPRIVATE(nslay)
48	! number of oceanic grid points
49	INTEGER, PUBLIC, SAVE :: knon
50	!$OMP THREADPRIVATE(knon)
51	! timestep for coupling (update slab temperature) in timesteps
52	INTEGER, PRIVATE, SAVE :: cpl_pas
53	!$OMP THREADPRIVATE(cpl_pas)
54	! cyang = 1/heat capacity of top layer (rho.c.H)
55	REAL, PRIVATE, SAVE :: cyang
56	!$OMP THREADPRIVATE(cyang)
57	! depth of slab layers (1st or 2nd layer)
58	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh
59	!$OMP THREADPRIVATE(slabh)
60	! slab temperature
61	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: tslab
62	!$OMP THREADPRIVATE(tslab)
63	! heat flux convergence due to Ekman
64	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_ekman
65	!$OMP THREADPRIVATE(dt_ekman)
66	! heat flux convergence due to horiz diffusion
67	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_hdiff
68	!$OMP THREADPRIVATE(dt_hdiff)
69	! heat flux convergence due to GM eddy advection
70	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_gm
71	!$OMP THREADPRIVATE(dt_gm)
72	! fraction of ocean covered by sea ice (sic / (oce+sic))
73	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic
74	!$OMP THREADPRIVATE(fsic)
75	! temperature of the sea ice
76	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice
77	!$OMP THREADPRIVATE(tice)
78	! sea ice thickness, in kg/m2
79	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: seaice
80	!$OMP THREADPRIVATE(seaice)
81	! net surface heat flux, weighted by open ocean fraction
82	! slab_bils accumulated over cpl_pas timesteps
83	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum
84	!$OMP THREADPRIVATE(bils_cum)
85	! net heat flux into the ocean below the ice : conduction + solar radiation
86	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bilg
87	!$OMP THREADPRIVATE(slab_bilg)
88	! slab_bilg cululated over cpl_pas timesteps
89	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bilg_cum
90	!$OMP THREADPRIVATE(bilg_cum)
91	! wind stress saved over cpl_pas timesteps
92	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: taux_cum
93	!$OMP THREADPRIVATE(taux_cum)
94	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tauy_cum
95	!$OMP THREADPRIVATE(tauy_cum)
96
97	!***********************************************************************************
98	! Parameters (could be read in def file: move to slab_init)
99	!***********************************************************************************
100	! snow and ice physical characteristics:
101	REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp [in K]
102	REAL, PARAMETER :: t_melt=273.15 ! melting ice temp [in K]
103	REAL, PARAMETER :: sno_den=300. !mean snow density [in kg/m3]
104	REAL, PARAMETER :: ice_den=917. ! ice density [in kg/m3]
105	REAL, PARAMETER :: sea_den=1026. ! sea water density [in kg/m3]
106	REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice [in W/(m.K) or (W.kg)/(K.m4)]
107	REAL, PRIVATE, SAVE :: sno_cond ! conductivity of snow [in W/(m.K) or (W.kg)/(K.m4)]
108	REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice [in J/(kg.K)]
109	REAL, PARAMETER :: sea_cap=3994. ! specific heat capacity, seawater [in J/(kg.K)]
110	REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice [in J/kg]
111	REAL, PARAMETER :: ice_sub=2834000. ! latent heat of sublimation for snow and ice [in J/kg]
112
113	! control of snow and ice cover & freeze / melt (heights in m converted to kg/m2)
114	REAL, PARAMETER, PUBLIC :: snow_min=0.05*sno_den ! critical snow height [in kg/m2]
115	REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean [in kg/m2]
116	REAL, PARAMETER :: ice_frac_min=0.001
117	REAL, PRIVATE, SAVE :: ice_frac_max ! Max ice fraction (leads)
118	REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness [in kg/m2]
119	REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness [in kg/m2]
120	! below ice_thin, priority is to melt lateral / grow height
121	! ice_thin is also height of new ice
122	REAL, PRIVATE, SAVE :: h_ice_thick ! thin ice thickness
123	! above ice_thick, priority is melt height / grow lateral
124	REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice [in kg/m2]
125	REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height [in kg/m2]
126
127	REAL, PARAMETER :: epsfra=1.0E-05 ! minimial grid fraction size below which there is no ice
128
129	REAL, PARAMETER, PUBLIC :: capcalocean=50.*4.228e+06 ! surface heat capacity [J.K-1.m-2] (assuming 50 m slab ocean)
130	REAL, PARAMETER, PUBLIC :: capcalseaice=5.1444e+06*0.15
131	REAL, PARAMETER, PUBLIC :: capcalsno=2.3867e+06*0.15
132
133	REAL, PARAMETER, PUBLIC :: h_alb_ice=0.3*ice_den ! height (in kg/m2) used in the calculation of sea ice albedo vs thickness
134	! (changed from 50cm to 30cm based on comparisons with Brandt et al. 2005) ; more info in the slab ocean wiki page
135	REAL, PARAMETER, PUBLIC :: h_sno_alb=0.02*sno_den ! height (in kg/m2) for control of snow fraction
136
137	REAL, PARAMETER, PUBLIC :: alb_ice_min=0.08 ! minimum sea ice albedo used for calculation of albedo as a function of sea ice thickness (https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Slab_ocean_model)
138
139	! Horizontal transport parameters
140	! flag for horizontal diffusion
141	LOGICAL, PUBLIC, SAVE :: slab_hdiff
142	!$OMP THREADPRIVATE(slab_hdiff)
143	! flag for GM eddy diffusivity
144	LOGICAL, PUBLIC, SAVE :: slab_gm
145	!$OMP THREADPRIVATE(slab_gm)
146	REAL, PRIVATE, SAVE :: coef_hdiff ! coefficient for horizontal diffusion
147	!$OMP THREADPRIVATE(coef_hdiff)
148	! flags for Ekman, conv adjustment
149	LOGICAL, PUBLIC, SAVE :: slab_ekman
150	!$OMP THREADPRIVATE(slab_ekman)
151	INTEGER, PUBLIC, SAVE :: slab_cadj
152	!$OMP THREADPRIVATE(slab_cadj)
153
154	!***********************************************************************************
155
156	CONTAINS
157	!
158	!***********************************************************************************
159	!
160	SUBROUTINE ocean_slab_init(dtime, pctsrf_rst, tslab_rst, tice_rst, seaice_rst, zmasq)
161
162	! This routine
163	! (1) allocates variables initialised from restart fields
164	! (2) allocates some other variables internal to the ocean module
165
166	USE ioipsl_getin_p_mod, ONLY : getin_p
167	USE mod_phys_lmdz_transfert_para, ONLY : gather
168	USE slab_heat_transp_mod, ONLY : ini_slab_transp
169
170	! Input variables
171	!***********************************************************************************
172	REAL, INTENT(IN) :: dtime
173	! Variables read from restart file
174	REAL, DIMENSION(klon), INTENT(IN) :: pctsrf_rst
175	REAL, DIMENSION(klon,nslay), INTENT(IN) :: tslab_rst
176	REAL, DIMENSION(klon), INTENT(IN) :: tice_rst
177	REAL, DIMENSION(klon), INTENT(IN) :: seaice_rst
178	REAL, DIMENSION(klon), INTENT(IN) :: zmasq
179
180	! Local variables
181	!************************************************************************************
182	INTEGER :: error
183	REAL, DIMENSION(klon_glo) :: zmasq_glo
184	CHARACTER (len = 80) :: abort_message
185	CHARACTER (len = 20) :: modname = 'ocean_slab_init'
186
187	!***********************************************************************************
188	! Define some parameters
189	!***********************************************************************************
190	!
191	! cpl_pas coupling period (update of tslab and ice fraction)
192	! for a calculation at each physical timestep, cpl_pas=1
193	cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
194	CALL getin_p('cpl_pas',cpl_pas)
195	print *,'cpl_pas',cpl_pas
196	! Number of slab layers
197	! nslay=2
198	! CALL getin_p('slab_layers',nslay)
199	print *,'number of slab layers : ',nslay
200	! Layer thickness
201	ALLOCATE(slabh(nslay), stat = error)
202	IF (error /= 0) THEN
203	abort_message='Pb allocation slabh'
204	CALL abort_physic(modname,abort_message,1)
205	ENDIF
206	slabh(1)=50. ! Height of first ocean layer (wind-mixed layer)
207	CALL getin_p('slab_depth',slabh(1))
208	IF (nslay.GT.1) THEN
209	slabh(2)=150. ! Height of second ocean layer (deep ocean layer)
210	END IF
211	! cyang = 1/heat capacity of top layer (rho.c.H)
212	cyang=1/(slabh(1)sea_densea_cap)
213
214	! ******** Sea Ice parameters *********
215	ice_frac_max = 0.999999 ! frac = 1 may lead to some problems.
216	CALL getin_p('ice_frac_max',ice_frac_max)
217	h_ice_thick = 1.5
218	CALL getin_p('h_ice_thick',h_ice_thick)
219	h_ice_thick = h_ice_thick * ice_den
220	sno_cond = 0.31
221	CALL getin_p('sno_cond',sno_cond)
222	sno_cond = sno_cond * sno_den
223
224	! ******** Heat Transport parameters **
225	! Ekman transport
226	! slab_ekman=0
227	slab_ekman=.FALSE.
228	CALL getin_p('slab_ekman',slab_ekman)
229	! GM eddy advection (2-layers only)
230	slab_gm=.FALSE.
231	CALL getin_p('slab_gm',slab_gm)
232	! IF (slab_ekman.LT.2) THEN
233	IF (.NOT.slab_ekman) THEN
234	slab_gm=.FALSE.
235	ENDIF
236	! Horizontal diffusion
237	slab_hdiff=.FALSE.
238	CALL getin_p('slab_hdiff',slab_hdiff)
239	IF (slab_gm) THEN
240	coef_hdiff=8000. ! non-dimensional; coef_hdiff should be 25000 if GM is off
241	ELSE
242	coef_hdiff=25000.
243	ENDIF
244	CALL getin_p('coef_hdiff',coef_hdiff)
245
246	! Convective adjustment
247	! IF (nslay.EQ.1) THEN
248	! slab_cadj=0
249	! ELSE
250	slab_cadj=1
251	! END IF
252	CALL getin_p('slab_cadj',slab_cadj)
253
254	!************************************************************************************
255	! Allocate surface fraction read from restart file
256	!************************************************************************************
257	ALLOCATE(fsic(klon), stat = error)
258	IF (error /= 0) THEN
259	abort_message='Pb allocation tmp_pctsrf_slab'
260	CALL abort_physic(modname,abort_message,1)
261	ENDIF
262	fsic(:)=0.
263	!zmasq = continent fraction
264	WHERE (1.-zmasq(:)>EPSFRA)
265	! fsic(:) = MIN(pctsrf_rst(:,is_sic)/(1.-zmasq(:)),ice_frac_max)
266	fsic(:) = MIN(pctsrf_rst(:)/(1.-zmasq(:)),ice_frac_max)
267	END WHERE
268
269	!************************************************************************************
270	! Allocate saved fields
271	!************************************************************************************
272	ALLOCATE(tslab(klon,nslay), stat=error)
273	IF (error /= 0) CALL abort_physic &
274	(modname,'pb allocation tslab', 1)
275	tslab(:,:) = tslab_rst(:,:)
276
277	ALLOCATE(bils_cum(klon), stat = error)
278	IF (error /= 0) THEN
279	abort_message='Pb allocation slab_bils_cum'
280	CALL abort_physic(modname,abort_message,1)
281	ENDIF
282	bils_cum(:) = 0.0
283
284	! IF (version_ocean=='sicINT') THEN ! interactive sea ice
285	ALLOCATE(slab_bilg(klon), stat = error)
286	IF (error /= 0) THEN
287	abort_message='Pb allocation slab_bilg'
288	CALL abort_physic(modname,abort_message,1)
289	ENDIF
290	slab_bilg(:) = 0.0
291	ALLOCATE(bilg_cum(klon), stat = error)
292	IF (error /= 0) THEN
293	abort_message='Pb allocation slab_bilg_cum'
294	CALL abort_physic(modname,abort_message,1)
295	ENDIF
296	bilg_cum(:) = 0.0
297	ALLOCATE(tice(klon), stat = error)
298	IF (error /= 0) THEN
299	abort_message='Pb allocation slab_tice'
300	CALL abort_physic(modname,abort_message,1)
301	ENDIF
302	tice(:) = tice_rst(:)
303	ALLOCATE(seaice(klon), stat = error)
304	IF (error /= 0) THEN
305	abort_message='Pb allocation slab_seaice'
306	CALL abort_physic(modname,abort_message,1)
307	ENDIF
308	seaice(:) = seaice_rst(:)
309	! END IF
310
311	IF (slab_hdiff) THEN !horizontal diffusion
312	ALLOCATE(dt_hdiff(klon,nslay), stat = error)
313	IF (error /= 0) THEN
314	abort_message='Pb allocation dt_hdiff'
315	CALL abort_physic(modname,abort_message,1)
316	ENDIF
317	dt_hdiff(:,:) = 0.0
318	ENDIF
319
320	IF (slab_gm) THEN !GM advection
321	ALLOCATE(dt_gm(klon,nslay), stat = error)
322	IF (error /= 0) THEN
323	abort_message='Pb allocation dt_gm'
324	CALL abort_physic(modname,abort_message,1)
325	ENDIF
326	dt_gm(:,:) = 0.0
327	ENDIF
328
329	! IF (slab_ekman.GT.0) THEN ! ekman transport
330	IF (slab_ekman) THEN ! ekman transport
331	ALLOCATE(dt_ekman(klon,nslay), stat = error)
332	IF (error /= 0) THEN
333	abort_message='Pb allocation dt_ekman'
334	CALL abort_physic(modname,abort_message,1)
335	ENDIF
336	dt_ekman(:,:) = 0.0
337	ALLOCATE(taux_cum(klon), stat = error)
338	IF (error /= 0) THEN
339	abort_message='Pb allocation taux_cum'
340	CALL abort_physic(modname,abort_message,1)
341	ENDIF
342	taux_cum(:) = 0.0
343	ALLOCATE(tauy_cum(klon), stat = error)
344	IF (error /= 0) THEN
345	abort_message='Pb allocation tauy_cum'
346	CALL abort_physic(modname,abort_message,1)
347	ENDIF
348	tauy_cum(:) = 0.0
349	ENDIF
350
351	! Initialize transport
352	IF (slab_hdiff.OR.(slab_ekman)) THEN
353	CALL gather(zmasq,zmasq_glo)
354	! Master thread/process only
355	!$OMP MASTER
356	IF (is_mpi_root) THEN
357	CALL ini_slab_transp(zmasq_glo)
358	END IF
359	!$OMP END MASTER
360	END IF
361
362	END SUBROUTINE ocean_slab_init
363	!
364	!***********************************************************************************
365	!
366	SUBROUTINE ocean_slab_frac(pctsrf_chg, zmasq)
367
368	! This routine sends back the sea ice and ocean fraction to the main physics routine.
369	! Called only with interactive sea ice.
370
371	! Arguments
372	!************************************************************************************
373	REAL, DIMENSION(klon), INTENT(IN) :: zmasq ! proxy of rnat
374	REAL, DIMENSION(klon), INTENT(OUT) :: pctsrf_chg ! sea-ice fraction
375
376	pctsrf_chg(:)=fsic(:)*(1.-zmasq(:))
377
378	END SUBROUTINE ocean_slab_frac
379	!
380	!************************************************************************************
381	!
382	SUBROUTINE ocean_slab_noice(itime, dtime, knon, knindex, &
383	precip_snow, tsurf_in, &
384	radsol, snow, fluxsens, &
385	tsurf_new, flux_u1, flux_v1, zmasq)
386
387	USE slab_heat_transp_mod, ONLY: divgrad_phy,slab_ekman2,slab_gmdiff
388	USE mod_phys_lmdz_para
389
390	! This routine
391	! (1) computes surface turbulent fluxes over points with some open ocean
392	! (2) reads additional Q-flux (everywhere)
393	! (3) computes horizontal transport (diffusion & Ekman)
394	! (4) updates slab temperature every cpl_pas ; creates new ice if needed.
395
396	! Note :
397	! klon total number of points
398	! knon number of points with open ocean (varies with time)
399	! knindex gives position of the knon points within klon.
400	! In general, local saved variables have klon values
401	! variables exchanged with PBL module have knon.
402
403	! Input arguments
404	!***********************************************************************************
405	INTEGER, INTENT(IN) :: itime ! current timestep INTEGER,
406	INTEGER, INTENT(IN) :: knon ! number of points
407	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
408	REAL, INTENT(IN) :: dtime ! timestep (s)
409	REAL, DIMENSION(klon), INTENT(IN) :: precip_snow !, precip_rain
410
411	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! surface temperature
412	REAL, DIMENSION(klon), INTENT(IN) :: radsol ! net surface (radiative) flux
413	REAL, DIMENSION(klon), INTENT(IN) :: flux_u1, flux_v1 ! Comes from turbdiff
414	REAL, DIMENSION(klon), INTENT(IN) :: fluxsens !, sensible heat flux
415	REAL, DIMENSION(klon), INTENT(IN) :: zmasq ! proxy of rnat
416
417	! In/Output arguments
418	!************************************************************************************
419	REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2
420
421	! Output arguments
422	!************************************************************************************
423	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! new surface tempearture
424
425	! Local variables
426	!************************************************************************************
427	INTEGER :: i,ki,k
428	REAL :: t_cadj
429
430	! for new ice creation
431	REAL :: e_freeze
432	REAL, DIMENSION(klon) :: slab_bils ! weighted surface heat flux
433	! horizontal diffusion and Ekman local vars
434	! dimension = global domain (klon_glo) instead of // subdomains
435	REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo,dt_ekman_glo,dt_gm_glo
436	! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
437	REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
438	REAL, DIMENSION(klon_glo,nslay) :: tslab_glo
439	REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo
440
441	!****************************************************************************************
442	! 1) Surface fluxes calculation
443	! Points with some open ocean only
444	!****************************************************************************************
445	! save total cumulated heat fluxes locally
446	! radiative + turbulent + melt of falling snow
447	slab_bils(:)=0.
448	DO i=1,knon
449	ki=knindex(i)
450	slab_bils(ki)=(1.-fsic(ki))*(fluxsens(ki)+radsol(ki) &
451	-precip_snow(ki)ice_lat(1.+snow_wfact*fsic(ki)))
452	bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
453	END DO
454
455	! Compute surface wind stress
456	! CALL calcul_flux_wind(knon, dtime, &
457	! u0, v0, u1, v1, gustiness, cdragm, &
458	! AcoefU, AcoefV, BcoefU, BcoefV, &
459	! p1lay, temp_air, &
460	! flux_u1, flux_v1)
461
462	! save cumulated wind stress
463	IF (slab_ekman) THEN
464	DO i=1,knon
465	ki=knindex(i)
466	taux_cum(ki)=taux_cum(ki)+flux_u1(ki)*(1.-fsic(ki))/cpl_pas
467	tauy_cum(ki)=tauy_cum(ki)+flux_v1(ki)*(1.-fsic(ki))/cpl_pas
468	END DO
469	ENDIF
470
471	!****************************************************************************************
472	! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file
473	! limit_slab.nc
474	!
475	!****************************************************************************************
476	! CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
477	! lmt_bils and diff_sst,siv saved by limit_slab
478	! qflux = total QFlux correction (in W/m2)
479	! IF (qflux_bnd.EQ.2) THEN
480	! dt_qflux(:,1) = lmt_bils(:,1)+diff_sst(:)/cyang/86400.
481	! dt_qflux_sic(:) = -diff_siv(:)ice_denice_lat/86400.
482	! ELSE
483	! dt_qflux(:,1) = lmt_bils(:,1)+diff_sst(:)/cyang/86400. &
484	! - diff_siv(:)ice_denice_lat/86400.
485	! END IF
486	! IF (nslay.GT.1) THEN
487	! dt_qflux(:,2:nslay)=lmt_bils(:,2:nslay)
488	! END IF
489
490	!****************************************************************************************
491	! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport)
492	! Bring to freezing temp and make sea ice if necessary
493	!
494	!*********************************************o***************************************
495	tsurf_new=tsurf_in
496	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
497	! ***********************************
498	! Horizontal transport
499	! ***********************************
500	IF (slab_ekman) THEN
501	! copy wind stress to global var
502	CALL gather(taux_cum,taux_glo)
503	CALL gather(tauy_cum,tauy_glo)
504	END IF
505
506	IF (slab_hdiff.OR.(slab_ekman)) THEN
507	CALL gather(tslab,tslab_glo)
508	! Compute horiz transport on one process only
509	IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus
510	IF (slab_hdiff) THEN
511	dt_hdiff_glo(:,:)=0.
512	END IF
513	IF (slab_ekman) THEN
514	dt_ekman_glo(:,:)=0.
515	END IF
516	IF (slab_gm) THEN
517	dt_gm_glo(:,:)=0.
518	END IF
519	DO i=1,cpl_pas ! time splitting for numerical stability
520	! IF (slab_ekman.GT.0) THEN
521	! SELECT CASE (slab_ekman)
522	! CASE (1) ! 1.5 layer scheme
523	! CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
524	! CASE (2) ! 2-layers
525	! CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm)
526	! CASE DEFAULT
527	! dt_ekman_tmp(:,:)=0.
528	! END SELECT
529	IF (slab_ekman) THEN
530	CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm)
531
532	dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:)
533	! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1
534	DO k=1,nslay
535	dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den)
536	ENDDO
537	tslab_glo=tslab_glo+dt_ekman_tmp*dtime
538	IF (slab_gm) THEN ! Gent-McWilliams eddy advection
539	dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:)
540	! convert dt from K.s-1.(kg.m-2) to K.s-1
541	DO k=1,nslay
542	dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/(slabh(k)*sea_den)
543	END DO
544	tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
545	END IF
546	ENDIF
547	IF (slab_hdiff) THEN ! horizontal diffusion
548	! laplacian of slab T
549	CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp)
550	! multiply by diff coef and normalize to 50m slab equivalent
551	dt_hdiff_tmp=dt_hdiff_tmpcoef_hdiff50./SUM(slabh)
552	dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:)
553	tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
554	END IF
555	END DO ! time splitting
556	IF (slab_hdiff) THEN
557	!dt_hdiff_glo saved in W/m2
558	DO k=1,nslay
559	dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)slabh(k)sea_den*sea_cap/cpl_pas
560	END DO
561	END IF
562	IF (slab_gm) THEN
563	!dt_hdiff_glo saved in W/m2
564	dt_gm_glo(:,:)=dt_gm_glo(:,:)*sea_cap/cpl_pas
565	END IF
566	IF (slab_ekman) THEN
567	! dt_ekman_glo saved in W/m2
568	dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
569	END IF
570	END IF ! master process
571	!$OMP BARRIER
572	! Send new fields back to all processes
573	CALL Scatter(tslab_glo,tslab)
574	IF (slab_hdiff) THEN
575	CALL Scatter(dt_hdiff_glo,dt_hdiff)
576	END IF
577	IF (slab_gm) THEN
578	CALL Scatter(dt_gm_glo,dt_gm)
579	END IF
580	IF (slab_ekman) THEN
581	CALL Scatter(dt_ekman_glo,dt_ekman)
582	! clear wind stress
583	taux_cum(:)=0.
584	tauy_cum(:)=0.
585	END IF
586	ENDIF ! transport
587
588	! ***********************************
589	! Other heat fluxes
590	! ***********************************
591	! Add cumulated ocean surface fluxes
592	tslab(:,1) = tslab(:,1) + bils_cum(:)cyangdtime
593	! Convective adjustment if 2 layers
594	IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN
595	DO i=1,klon
596	IF (tslab(i,2).GT.tslab(i,1)) THEN
597	! mean (mass-weighted) temperature
598	t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:))
599	tslab(i,1)=t_cadj
600	tslab(i,2)=t_cadj
601	END IF
602	END DO
603	END IF
604	! Add read QFlux
605	! IF (qflux_bnd.EQ.1) THEN
606	! ! QFlux from ocean circ. cannot cool tslab below freezing.
607	! dq_freeze = (t_freeze - tslab(:,1)) / (cyangdtimecpl_pas)
608	! dt_qflux(:,1) = MAX(dt_qflux(:,1), MIN(dq_freeze,0.))
609	! ELSE IF (qflux_bnd.EQ.2) THEN
610	! dq_freeze = (t_freeze - tslab(:,1)) / (cyangdtimecpl_pas) &
611	! - dt_qflux_sic(:)
612	! dt_qflux(:,1) = MAX(dt_qflux(:,1), MIN(dq_freeze,0.)) &
613	! + dt_qflux_sic(:)
614	! END IF
615	! DO k=1,nslay
616	! tslab(:,k) = tslab(:,k) + dt_qflux(:,k)cyangdtime*cpl_pas &
617	! * slabh(1)/slabh(k)
618	! END DO
619
620	! ***********************************
621	! Update surface temperature and ice
622	! ***********************************
623	! SELECT CASE(version_ocean)
624	! CASE('sicNO') ! no sea ice even below freezing !
625	! DO i=1,knon
626	! ki=knindex(i)
627	! tsurf_new(i)=tslab(ki,1)
628	! END DO
629	! CASE('sicOBS') ! "realistic" case, for prescribed sea ice
630	! ! tslab cannot be below freezing
631	! DO i=1,knon
632	! ki=knindex(i)
633	! IF (tslab(ki,1).LT.t_freeze) THEN
634	! tslab(ki,1)=t_freeze
635	! END IF
636	! tsurf_new(i)=tslab(ki,1)
637	! END DO
638	! CASE('sicINT') ! interactive sea ice
639	DO i=1,knon
640	ki=knindex(i)
641	! Check if new slab temperature is below freezing
642	IF (tslab(ki,1).LT.t_freeze) THEN
643	! We need to form ice now over ice-free points
644	! Else points not seen by slab_ice
645	IF (fsic(ki)*(1.-zmasq(ki)).LT.epsfra) THEN
646	! No ice present yet.
647	! quantity of new ice formed
648	e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat &
649	+fsic(ki)*seaice(ki)
650	! new ice forms at h_ice_thin
651	tsurf_new(ki)=t_freeze
652	tice(ki)=t_freeze
653	fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
654	IF (fsic(ki).GT.ice_frac_min) THEN
655	seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
656	tslab(ki,1)=t_freeze
657	ELSE
658	fsic(ki)=0.
659	END IF
660	END IF ! sea ice present
661	! if ice already present, new ice formed in slab_ice routine.
662	! ELSE ! temperature above freezing
663	! tsurf_new(i)=tslab(ki,1)
664	END IF
665	END DO
666	! END SELECT
667	bils_cum(:)=0.0! clear cumulated fluxes
668	END IF ! coupling time
669	END SUBROUTINE ocean_slab_noice
670	!
671	!*****************************************************************************
672
673	! SUBROUTINE ocean_slab_ice( &
674	! itime, dtime, jour, knon, knindex, &
675	! tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
676	! AcoefH, AcoefQ, BcoefH, BcoefQ, &
677	! AcoefU, AcoefV, BcoefU, BcoefV, &
678	! ps, u1, v1, gustiness, &
679	! radsol, snow, qsurf, qsol, agesno, &
680	! alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
681	! tsurf_new, dflux_s, dflux_l, swnet)
682
683	SUBROUTINE ocean_slab_ice(itime, dtime, knon, knindex, &
684	precip_snow, tsurf_in, &
685	radsol, snow, fluxsens, &
686	tsurf_new, evap, flux_u1, flux_v1, zmasq)
687
688	! USE calcul_fluxs_mod
689
690	! INCLUDE "YOMCST.h"
691	! INCLUDE "clesphys.h"
692
693	! This routine
694	! (1) computes surface turbulent fluxes over points with some sea ice
695	! (2) computes condutive fluxes in the snow and ice, and ice-ocean flux
696	! (3) computes snow/ice albedo
697	! (4) updates snow/ice temperature, surface melt if needed.
698	! (5) lateral growth if tslab < freezing
699	! (6) bottom & side melt / growth depending on bottom fluxes
700	! (7) updates slab temperature every cpl_pas
701
702	! Note :
703	! klon total number of points
704	! knon number of points with open ocean (varies with time)
705	! knindex gives position of the knon points within klon.
706	! In general, local saved variables have klon values
707	! variables exchanged with PBL module have knon.
708
709	! Input arguments
710	!****************************************************************************************
711	INTEGER, INTENT(IN) :: itime, knon !, jour
712	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
713	REAL, INTENT(IN) :: dtime
714	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
715	! REAL, DIMENSION(klon), INTENT(IN) :: p1lay
716	! REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
717	REAL, DIMENSION(klon), INTENT(IN) :: precip_snow !, precip_rain
718	REAL, DIMENSION(klon), INTENT(IN) :: evap, fluxsens!,fluxlat
719	REAL, DIMENSION(klon), INTENT(IN) :: flux_u1, flux_v1
720	REAL, DIMENSION(klon), INTENT(IN) :: zmasq ! proxy of rnat
721	! REAL, DIMENSION(klon), INTENT(IN) :: spechum, temp_air
722	! REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
723	! REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
724	! REAL, DIMENSION(klon), INTENT(IN) :: ps
725	! REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
726	! REAL, DIMENSION(klon), INTENT(IN) :: swnet
727
728	! In/Output arguments
729	!****************************************************************************************
730	REAL, DIMENSION(klon), INTENT(INOUT) :: snow!, qsol
731	! REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
732	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
733
734	! Output arguments
735	!****************************************************************************************
736	! REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
737	! REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval
738	! REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval
739	! REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens!, fluxlat
740	! REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
741	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
742	! REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
743
744	! Local variables
745	!****************************************************************************************
746	INTEGER :: i,ki
747	! REAL, DIMENSION(klon) :: cal, beta, dif_grnd
748	! REAL, DIMENSION(klon) :: u0, v0
749	! REAL, DIMENSION(klon) :: u1_lay, v1_lay
750	REAL, DIMENSION(klon) :: f_bot ! bottom ocean - ice flux
751
752	! intermediate heat fluxes:
753	! (conduction snow / ice, shortwave penetration, bottom turbulent fluxes)
754	REAL :: f_cond_s, f_cond_i, f_turb
755	! friction velocity, 1/C and 1/tau conduction for ice
756	REAL :: ustar
757	! REAL :: uscap, ustau
758	! for snow/ice albedo:
759	! REAL :: alb_snow, alb_ice, alb_pond
760	! REAL :: frac_snow, frac_ice, frac_pond
761	! REAL :: z1_i, z2_i, z1_s, zlog ! height parameters
762	! for ice melt / freeze
763	REAL :: e_melt, e_freeze, snow_evap, h_test, h_new
764	! dhsic, dfsic change in ice mass, fraction.
765	! frac_mf ratio of lateral / thickness growth / melt
766	REAL :: dhsic, dfsic, frac_mf
767
768	!*******************************************************************************
769	! 1) Update surface , ice and slab temperature
770	!*******************************************************************************
771	! Wind stress
772	! flux_u1, flux_v1 from physics
773	! save cumulated wind stress
774	! Use ocean-ice stress = wind stress * (1.-fsic)
775	! IF (slab_ekman.GT.0) THEN
776	IF (slab_ekman) THEN
777	DO i=1,knon
778	ki=knindex(i)
779	taux_cum(ki)=taux_cum(ki)+flux_u1(ki)fsic(ki)(1.-fsic(ki))/cpl_pas
780	tauy_cum(ki)=tauy_cum(ki)+flux_v1(ki)fsic(ki)(1.-fsic(ki))/cpl_pas
781	END DO
782	ENDIF
783
784	! Initialize ice-ocean flux
785	slab_bilg(:)=0.
786
787	! Old, explicit scheme for snow & ice conductive fluxes
788	! radsol is total surface fluxes (input) : radiative + turbulent
789	DO i=1,knon
790	ki=knindex(i) ! For PCM : you can probably replace ki by i
791	! ocean-ice turbulent heat flux
792	! turbulent velocity = (tau/rho)^1/2
793	ustar = MAX(5e-4, SQRT((1.-fsic(ki)) &
794	* SQRT(flux_u1(ki)2 + flux_v1(ki)2) / sea_den))
795	f_turb = 0.0057 * sea_den * sea_cap * ustar * (tslab(ki,1) - t_freeze)
796	! f_turb >0 and cannot bring tslab below zero
797	f_turb = MAX(0., MIN(f_turb, &
798	(tslab(ki,1)-t_freeze) / (cyangdtimecpl_pas)))
799
800	! ice conductive flux (pos up)
801	IF (seaice(ki).GT.0) THEN
802	f_cond_i = ice_cond*(t_freeze-tice(ki))/seaice(ki)
803	ELSE
804	f_cond_i = 0
805	END IF
806
807	! If snow layer present, tsurf = tsnow
808	! Problem here, as tsurf_in # tsnow ?
809	IF (snow(ki).GT.snow_min) THEN
810	! snow conductive flux (pos up)
811	f_cond_s=sno_cond*(tice(ki)-tsurf_in(ki))/snow(ki)
812	! update ice temperature
813	tice(ki)=tice(ki) + 2./ice_cap/(snow(ki)+seaice(ki)) &
814	(f_cond_i-f_cond_s)dtime
815	! update snow temperature
816	tsurf_new(ki) = tsurf_in(ki) + 2./ice_cap/snow(ki) &
817	(fluxsens(ki)+radsol(ki)+f_cond_s)dtime
818	ELSE IF (seaice(ki).GT.0) THEN ! bare ice. tsurf = tice
819	! update ice temperature
820	tice(ki) = tice(ki) + 2./ice_cap/seaice(ki) &
821	(fluxsens(ki)+radsol(ki)+f_cond_i)dtime
822	tsurf_new(ki) = tice(ki)
823	END IF
824	! bottom flux (used to grow ice from below)
825	f_bot(ki) = f_turb - f_cond_i
826	slab_bilg(ki) = -f_turb
827	END DO
828	!
829	!! Surface turbulent fluxes (sens, lat etc) and update surface temp.
830	! dif_grnd(:)=0.
831	! beta(:) = 1.
832	! CALL calcul_fluxs(knon, is_sic, dtime, &
833	! tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, &
834	! precip_rain, precip_snow, snow, qsurf, &
835	! radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
836	! f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
837	! tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
838	! DO i=1,knon
839	! IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
840	! END DO
841
842	! Surface, snow-ice and ice-ocean fluxes.
843	! Prepare call to calcul_fluxs (cal, beta, radsol, dif_grnd)
844
845	! Surface turbulent fluxes (sens, lat etc) and update surface temp.
846	! beta(:) = 1.
847	! CALL calcul_fluxs(knon, is_sic, dtime, &
848	! tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, &
849	! precip_rain, precip_snow, snow, qsurf, &
850	! radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
851	! f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
852	! tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
853
854	!! Update remaining temperature and fluxes
855	! DO i=1,knon
856	! ki=knindex(i)
857	! ! ocean-ice turbulent heat flux
858	! ! turbulent velocity = (tau/rho)^1/2 for low ice cover
859	! ! min = 5e-4 for 1cm/s current
860	! ustar = MAX(5e-4, SQRT((1.-fsic(ki)) &
861	! * SQRT(flux_u1(i)2 + flux_v1(i)2) / sea_den))
862	! f_turb = 0.0057 * sea_den * sea_cap * ustar * (tslab(ki,1) - t_freeze)
863	! ! ice_ocean fluxes cannot bring tslab below freezing
864	! f_turb = MAX(0., MIN(f_turb, slab_bilg(ki) + (tslab(ki,1)-t_freeze) &
865	! / (fsic(ki)cyangdtime*cpl_pas)))
866	! IF (snow(i).GT.snow_min) THEN
867	! ! snow conductive flux after calcul_fluxs
868	! f_cond_s = sno_cond * (tice(ki)-tsurf_new(i)) / snow(i)
869	! ! 1 / heat capacity and conductive timescale
870	! uscap = 2. / ice_cap / (snow(i)+seaice(ki))
871	! ustau = uscap * ice_cond / seaice(ki)
872	! ! update ice temp
873	! tice(ki) = (tice(ki) + dtime(ustaut_freeze - uscap*f_cond_s)) &
874	! / (1 + dtime*ustau)
875	! ELSE ! bare ice
876	! tice(ki)=tsurf_new(i)
877	! END IF
878	! ! ice conductive flux (pos up)
879	! f_cond_i = ice_cond * (t_freeze-tice(ki)) / seaice(ki)
880	! f_bot(i) = f_turb - f_cond_i
881	! slab_bilg(ki) = slab_bilg(ki)-f_turb
882	! END DO
883
884	! weight fluxes to ocean by sea ice fraction
885	slab_bilg(:)=slab_bilg(:)*fsic(:)
886
887	!****************************************************************************************
888	! 2) Update snow and ice surface : thickness and fraction
889	!****************************************************************************************
890	DO i=1,knon
891	ki=knindex(i)
892	! snow precip (could be before fluxes above ?)
893	IF (precip_snow(ki) > 0.) THEN
894	snow(ki) = snow(ki)+precip_snow(ki)dtime(1.-snow_wfact*(1.-fsic(ki)))
895	END IF
896	! snow and ice sublimation
897	IF (evap(ki) > 0.) THEN
898	snow_evap = MIN (snow(ki) / dtime, evap(ki))
899	snow(ki) = snow(ki) - snow_evap * dtime
900	snow(ki) = MAX(0.0, snow(ki))
901	seaice(ki) = MAX(0.0,seaice(ki)-(evap(ki)-snow_evap)*dtime)
902	ENDIF
903	! Melt / Freeze snow from above if Tsurf>0
904	IF (tsurf_new(ki).GT.t_melt) THEN
905	! energy available for melting snow (in kg of melted snow /m2)
906	e_melt = MIN(MAX(snow(ki)(tsurf_new(ki)-t_melt)ice_cap/2. &
907	/(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(ki))
908	! remove snow
909	IF (snow(ki).GT.e_melt) THEN
910	snow(ki)=snow(ki)-e_melt
911	tsurf_new(ki)=t_melt
912	ELSE ! all snow is melted
913	! add remaining heat flux to ice
914	e_melt=e_melt-snow(ki)
915	snow(ki)=0.0
916	tice(ki)=tice(ki)+e_meltice_lat2./(ice_cap*seaice(ki))
917	tsurf_new(ki)=tice(ki)
918	END IF
919	END IF
920	! Bottom melt / grow
921	! bottom freeze if bottom flux (cond + oce-ice) <0
922	IF (f_bot(ki).LT.0) THEN
923	! larger fraction of bottom growth
924	frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thick) &
925	/ (h_ice_max-h_ice_thick)))
926	! quantity of new ice (formed at mean ice temp)
927	e_melt= -f_bot(ki) * dtime * fsic(ki) &
928	/ (ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
929	! first increase height to h_thick
930	dhsic=MAX(0.,MIN(h_ice_thick-seaice(ki),e_melt/fsic(ki)))
931	seaice(ki)=dhsic+seaice(ki)
932	e_melt=e_melt-fsic(ki)*dhsic
933	IF (e_melt.GT.0.) THEN
934	! frac_mf fraction used for lateral increase
935	dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
936	fsic(ki)=fsic(ki)+dfsic
937	e_melt=e_melt-dfsic*seaice(ki)
938	! rest used to increase height
939	seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
940	END IF
941	ELSE
942	! melt from below if bottom flux >0
943	! larger fraction of lateral melt from warm ocean
944	frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin) &
945	/ (h_ice_thick-h_ice_thin)))
946	! bring ice to freezing and melt from below
947	! quantity of melted ice
948	e_melt= f_bot(ki) * dtime * fsic(ki) &
949	/ (ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
950	! first decrease height to h_thick
951	IF (fsic(ki).GT.0) THEN
952	dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
953	ELSE
954	dhsic=0
955	ENDIF
956	seaice(ki)=seaice(ki)-dhsic
957	e_melt=e_melt-fsic(ki)*dhsic
958	IF (e_melt.GT.0) THEN
959	! frac_mf fraction used for height decrease
960	dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
961	seaice(ki)=seaice(ki)-dhsic
962	e_melt=e_melt-fsic(ki)*dhsic
963	! rest used to decrease fraction (up to 0!)
964	dfsic=MIN(fsic(ki),e_melt/seaice(ki))
965	! keep remaining in ocean if everything melted
966	e_melt=e_melt-dfsic*seaice(ki)
967	slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime
968	ELSE
969	dfsic=0
970	END IF
971	fsic(ki)=fsic(ki)-dfsic
972	END IF
973	! melt ice from above if Tice>0
974	IF (tice(ki).GT.t_melt) THEN
975	! quantity of ice melted (kg/m2)
976	e_melt=MAX(seaice(ki)(tice(ki)-t_melt)ice_cap/2. &
977	/(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
978	! melt from above, height only
979	dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
980	e_melt=e_melt-dhsic
981	IF (e_melt.GT.0) THEN
982	! lateral melt if ice too thin
983	dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
984	! if all melted add remaining heat to ocean
985	e_melt=MAX(0.,e_meltfsic(ki)-dfsich_ice_min)
986	slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime
987	! update height and fraction
988	fsic(ki)=fsic(ki)-dfsic
989	END IF
990	seaice(ki)=seaice(ki)-dhsic
991	! surface temperature at melting point
992	tice(ki)=t_melt
993	tsurf_new(ki)=t_melt
994	END IF
995	! convert snow to ice if below floating line
996	h_test=(seaice(ki)+snow(ki))ice_den-seaice(ki)sea_den
997	IF ((h_test.GT.0.).AND.(seaice(ki).GT.h_ice_min)) THEN !snow under water
998	! extra snow converted to ice (with added frozen sea water)
999	dhsic=h_test/(sea_den-ice_den+sno_den)
1000	seaice(ki)=seaice(ki)+dhsic
1001	snow(ki)=snow(ki)-dhsic*sno_den/ice_den
1002	! available energy (freeze sea water + bring to tice)
1003	e_melt=dhsic(1.-sno_den/ice_den)(ice_lat+ &
1004	ice_cap/2.*(t_freeze-tice(ki)))
1005	! update ice temperature
1006	tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(ki)+seaice(ki))
1007	END IF
1008	END DO
1009
1010	!*******************************************************************************
1011	! 3) cumulate ice-ocean fluxes, update tslab, lateral grow
1012	!*********************************************o*****************************
1013	!cumul fluxes.
1014	bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
1015	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab
1016	! Add cumulated surface fluxes
1017	tslab(:,1)=tslab(:,1)+bilg_cum(:)cyangdtime
1018	bilg_cum(:)=0.
1019	! If slab temperature below freezing, new lateral growth
1020	DO i=1,knon
1021	ki=knindex(i)
1022	IF (tslab(ki,1).LT.t_freeze) THEN ! create more ice
1023	! quantity of new ice formed over open ocean
1024	! (formed at mean ice temperature)
1025	e_freeze=(t_freeze-tslab(ki,1))/cyang &
1026	/(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
1027	! new ice height and fraction
1028	h_new=MAX(h_ice_thin,MIN(h_ice_new,seaice(ki))) ! new height
1029	! h_new=MIN(h_ice_new,seaice(ki))
1030	dfsic=MIN(ice_frac_max-fsic(ki) &
1031	,MAX(ice_frac_min,e_freeze/h_new))
1032	! update average sea ice height
1033	seaice(ki)=(seaice(ki)*fsic(ki)+e_freeze) &
1034	/ (dfsic+fsic(ki))
1035	! update snow thickness
1036	snow(ki) = snow(ki) * fsic(ki) / (dfsic+fsic(ki))
1037	! update tslab to freezing
1038	tslab(ki,1)=t_freeze
1039	! update sea ice fraction
1040	fsic(ki)=fsic(ki)+dfsic
1041	END IF ! tslab below freezing
1042	END DO
1043	END IF ! coupling time
1044
1045	!****************************************************************************************
1046	! 4) Compute sea-ice and snow albedo
1047	!****************************************************************************************
1048	! Removed all albedo stuff as it is computed through hydrol in the Generic model
1049	! ------ End Albedo ----------
1050
1051	!tests remaining ice fraction
1052	WHERE ((fsic.LT.ice_frac_min).OR.(seaice.LT.h_ice_min))
1053	tslab(:,1)=tslab(:,1)-fsicseaiceice_lat*cyang
1054	tice=t_melt
1055	fsic=0.
1056	seaice=0.
1057	END WHERE
1058
1059	END SUBROUTINE ocean_slab_ice
1060	!
1061	!****************************************************************************************
1062	!
1063	SUBROUTINE ocean_slab_final
1064
1065	!****************************************************************************************
1066	! Deallocate module variables
1067	!****************************************************************************************
1068	IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
1069	IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
1070	IF (ALLOCATED(tice)) DEALLOCATE(tice)
1071	IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
1072	IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
1073	IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
1074	IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
1075	IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum)
1076	IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum)
1077	IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
1078	IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
1079	IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm)
1080	! IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux)
1081	! IF (ALLOCATED(dt_qflux_sic)) DEALLOCATE(dt_qflux_sic)
1082
1083	END SUBROUTINE ocean_slab_final
1084	!
1085	!****************************************************************************************
1086	!
1087	SUBROUTINE ocean_slab_get_vars(ngrid,tslab_loc, seaice_loc, flux_g_loc, &
1088	dt_hdiff_loc,dt_ekman_loc)
1089
1090	! "Get some variables from module ocean_slab_mod"
1091
1092	INTEGER, INTENT(IN) :: ngrid
1093	REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: tslab_loc
1094	REAL, DIMENSION(ngrid), INTENT(OUT) :: seaice_loc
1095	REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_g_loc
1096	REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: dt_hdiff_loc ! [in W/m2]
1097	REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: dt_ekman_loc ! [in W/m2]
1098	INTEGER :: i
1099
1100
1101	! Set the output variables
1102	tslab_loc(:,:) = 0.
1103	dt_hdiff_loc(:,:)=0.
1104	dt_ekman_loc(:,:)=0.
1105	tslab_loc(:,:) = tslab(:,:)
1106	seaice_loc(:) = seaice(:)
1107	flux_g_loc(:) = slab_bilg(:)
1108
1109	!! dt_hdiff_loc(:,i) = dt_hdiff(:,i)slabh(i)1000.*4228. !Convert en W/m2
1110	!! dt_ekman_loc(:,i) = dt_ekman(:,i)slabh(i)1000.*4228.
1111
1112	IF (slab_hdiff) THEN
1113	DO i=1,nslay
1114	dt_hdiff_loc(:,i) = dt_hdiff(:,i)
1115	ENDDO
1116	ENDIF
1117	IF (slab_ekman) THEN
1118	DO i=1,nslay
1119	dt_ekman_loc(:,i) = dt_ekman(:,i)
1120	ENDDO
1121	ENDIF
1122
1123
1124
1125	END SUBROUTINE ocean_slab_get_vars
1126	!
1127	!****************************************************************************************
1128	!
1129
1130	END MODULE ocean_slab_mod

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

Download in other formats:

Original Format