Context Navigation

ocean_slab_mod.F90 @ 5101

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

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

Download in other formats:

Original Format