Context Navigation

ocean_slab_mod.F90 @ 5130

Last change on this file since 5130 was 3780, checked in by evignon, 4 years ago
Premiere comission Etienne: changements pour le 1D (forcage en Ts au dessus des continents) et inclusion drag arbres dans yamada4_num=6
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 39.8 KB

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

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

Download in other formats:

Original Format