Context Navigation

ocean_slab_mod.F90 @ 3233

Last change on this file since 3233 was 3100, checked in by bhatnags, 20 months ago

Generic-PCM

This commit updates the slab ocean module to a parallelisable dynamic slab ocean module. This is particularly relevant if you want to be able to use oceanic heat transport in parallel mode.

It has the following features:

(a) Computes sea ice creation and evolution.
(b) Snow has thermodynamic properties.
(c) Computes oceanic horizontal transport (diffusion & surface-wind driven Ekman transport).
(d) Can be used in parallel mode.

Required callphys.def flags:
The slab ocean and its dependencies can be activated with the following flags (already added to deftank):
## Ocean options
## ~
# Model slab-ocean (Main flag for slab ocean)
ok_slab_ocean = .true.
# The following flags can only be set to true if ok_slab_ocean is true
# Ekman transport
slab_ekman = .true.
# Ekman zonal advection
slab_ekman_zonadv = .true.
# Horizontal diffusion (default coef_hdiff=25000., can be changed)
slab_hdiff = .true.
# Slab-ocean timestep (in physics timesteps)
cpl_pas = 1
# Gent-McWilliams? Scheme (can only be true if slab_ekman is true)
slab_gm = .true.

Notes:
In the current state, the model crashes if moistadjustment = .true. Unsure whether this is due to the slab or is an inherent issue with moistadj (under investigation).

SB and EM

Property svn:executable set to *

File size: 46.0 KB

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

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

Download in other formats:

Original Format