Context Navigation

ocean_slab_mod.F90 @ 3001

Last change on this file since 3001 was 2656, checked in by Ehouarn Millour, 8 years ago

Making the slab work:

added a slab_heat_transp_mod module for horizontal diffusion and Ekman transport
added storage and output of relevent variables in phyredem, phyetat0, phy_output_ctrlout_mod, phys_output_write_mod
moved nslay (number of slab layers) out of dimphy into ocean_slab_mod.

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:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 37.5 KB

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

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

Download in other formats:

Original Format