Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

ocean_slab_mod.F90

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

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

Download in other formats:

Original Format