Context Navigation

ocean_slab_mod.F90 @ 3580

Last change on this file since 3580 was 3397, checked in by bhatnags, 6 months ago

GENERIC-PCM: Fix/clarification of the usage of tice (sea ice surface temperature) vs. tsea_ice (temperature of the top layer (ice or snow)).

Property svn:executable set to *

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

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

Download in other formats:

Original Format