Context Navigation

ocean_slab_mod.F90 @ 3567

Last change on this file since 3567 was 3397, checked in by bhatnags, 5 months ago

GENERIC-PCM: Fix/clarification of the usage of tice (sea ice surface temperature) vs. tsea_ice (temperature of the top layer (ice or snow)).

Property svn:executable set to *

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

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

Download in other formats:

Original Format