Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

pbl_surface_mod.F90

Last change on this file was 4368, checked in by lguez, 19 months ago
Sync latest trunk changes to Ocean_skin
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: 148.1 KB

Rev	Line
[781]	1	!
[1279]	2	! $Id: pbl_surface_mod.F90 4368 2022-12-05 23:01:16Z abarral $
	3	!
[781]	4	MODULE pbl_surface_mod
	5	!
	6	! Planetary Boundary Layer and Surface module
	7	!
[2852]	8	! This module manages the calculation of turbulent diffusion in the boundary layer
[781]	9	! and all interactions towards the differents sub-surfaces.
	10	!
	11	!
	12	USE dimphy
	13	USE mod_phys_lmdz_para, ONLY : mpi_size
[2181]	14	USE mod_grid_phy_lmdz, ONLY : klon_glo
[781]	15	USE ioipsl
[4368]	16	USE surface_data, ONLY : type_ocean, ok_veget, landice_opt
[781]	17	USE surf_land_mod, ONLY : surf_land
	18	USE surf_landice_mod, ONLY : surf_landice
	19	USE surf_ocean_mod, ONLY : surf_ocean
	20	USE surf_seaice_mod, ONLY : surf_seaice
	21	USE cpl_mod, ONLY : gath2cpl
	22	USE climb_hq_mod, ONLY : climb_hq_down, climb_hq_up
	23	USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up
	24	USE coef_diff_turb_mod, ONLY : coef_diff_turb
[4013]	25	USE ioipsl_getin_p_mod, ONLY : getin_p
	26	USE cdrag_mod
	27	USE stdlevvar_mod
	28	USE wx_pbl_var_mod, ONLY : wx_pbl_init, wx_pbl_final, &
	29	wx_pbl_prelim_0, wx_pbl_prelim_beta
	30	USE wx_pbl_mod, ONLY : wx_pbl0_merge, wx_pbl_split, wx_pbl_dts_merge, &
	31	wx_pbl_check, wx_pbl_dts_check, wx_evappot
[3458]	32	use config_ocean_skin_m, only: activate_ocean_skin
[781]	33
	34	IMPLICIT NONE
	35
	36	! Declaration of variables saved in restart file
[888]	37	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: fder ! flux drift
[781]	38	!$OMP THREADPRIVATE(fder)
[4013]	39	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: snow ! snow at surface
[781]	40	!$OMP THREADPRIVATE(snow)
[888]	41	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qsurf ! humidity at surface
[781]	42	!$OMP THREADPRIVATE(qsurf)
[4013]	43	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ftsoil ! soil temperature
[781]	44	!$OMP THREADPRIVATE(ftsoil)
[4013]	45	REAL, ALLOCATABLE, DIMENSION(:), SAVE :: ydTs0, ydqs0
	46	! nul forced temperature and humidity differences
	47	!$OMP THREADPRIVATE(ydTs0, ydqs0)
[781]	48
[3391]	49	INTEGER, SAVE :: iflag_pbl_surface_t2m_bug
[2896]	50	!$OMP THREADPRIVATE(iflag_pbl_surface_t2m_bug)
[4013]	51	INTEGER, SAVE :: iflag_new_t2mq2m
	52	!$OMP THREADPRIVATE(iflag_new_t2mq2m)
	53
[2952]	54	!FC
	55	! integer, save :: iflag_frein
	56	! !$OMP THREADPRIVATE(iflag_frein)
[2896]	57
[781]	58	CONTAINS
	59	!
	60	!****************************************************************************************
	61	!
[2243]	62	SUBROUTINE pbl_surface_init(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
[781]	63
	64	! This routine should be called after the restart file has been read.
	65	! This routine initialize the restart variables and does some validation tests
	66	! for the index of the different surfaces and tests the choice of type of ocean.
	67
[1785]	68	USE indice_sol_mod
[2311]	69	USE print_control_mod, ONLY: lunout
[2898]	70	USE ioipsl_getin_p_mod, ONLY : getin_p
[3605]	71	IMPLICIT NONE
[1785]	72
[781]	73	INCLUDE "dimsoil.h"
	74
	75	! Input variables
	76	!****************************************************************************************
	77	REAL, DIMENSION(klon), INTENT(IN) :: fder_rst
	78	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: snow_rst
	79	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: qsurf_rst
	80	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
	81
	82	! Local variables
	83	!****************************************************************************************
	84	INTEGER :: ierr
	85	CHARACTER(len=80) :: abort_message
	86	CHARACTER(len = 20) :: modname = 'pbl_surface_init'
	87
	88	!****************************************************************************************
	89	! Allocate and initialize module variables with fields read from restart file.
	90	!
	91	!****************************************************************************************
	92	ALLOCATE(fder(klon), stat=ierr)
[2311]	93	IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
[781]	94
	95	ALLOCATE(snow(klon,nbsrf), stat=ierr)
[2311]	96	IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
[781]	97
	98	ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
[2311]	99	IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
[781]	100
	101	ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
[2311]	102	IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
[781]	103
[4013]	104	ALLOCATE(ydTs0(klon), stat=ierr)
	105	IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
	106
	107	ALLOCATE(ydqs0(klon), stat=ierr)
	108	IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
	109
[781]	110	fder(:) = fder_rst(:)
	111	snow(:,:) = snow_rst(:,:)
	112	qsurf(:,:) = qsurf_rst(:,:)
	113	ftsoil(:,:,:) = ftsoil_rst(:,:,:)
[4013]	114	ydTs0(:) = 0.
	115	ydqs0(:) = 0.
[781]	116
	117	!****************************************************************************************
	118	! Test for sub-surface indices
	119	!
	120	!****************************************************************************************
	121	IF (is_ter /= 1) THEN
	122	WRITE(lunout,)" Warning *"
	123	WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
	124	WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
	125	abort_message="voir ci-dessus"
[2311]	126	CALL abort_physic(modname,abort_message,1)
[781]	127	ENDIF
	128
	129	IF ( is_oce > is_sic ) THEN
	130	WRITE(lunout,)' Warning *'
	131	WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
	132	WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
	133	WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
	134	abort_message='voir ci-dessus'
[2311]	135	CALL abort_physic(modname,abort_message,1)
[781]	136	ENDIF
	137
	138	IF ( is_lic > is_sic ) THEN
	139	WRITE(lunout,)' Warning *'
	140	WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
	141	WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
	142	WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
	143	abort_message='voir ci-dessus'
[2311]	144	CALL abort_physic(modname,abort_message,1)
[781]	145	ENDIF
	146
	147	!****************************************************************************************
	148	! Validation of ocean mode
	149	!
	150	!****************************************************************************************
	151
[996]	152	IF (type_ocean /= 'slab ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN
[1064]	153	WRITE(lunout,)' Warning *'
	154	WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean
	155	abort_message='option pour l''ocean non valable'
[2311]	156	CALL abort_physic(modname,abort_message,1)
[781]	157	ENDIF
	158
[2896]	159	iflag_pbl_surface_t2m_bug=0
	160	CALL getin_p('iflag_pbl_surface_t2m_bug',iflag_pbl_surface_t2m_bug)
[4013]	161	WRITE(lunout,*) 'iflag_pbl_surface_t2m_bug=',iflag_pbl_surface_t2m_bug
[2952]	162	!FC
	163	! iflag_frein = 0
	164	! CALL getin_p('iflag_frein',iflag_frein)
	165	!
[3179]	166	!jyg<
	167	!****************************************************************************************
	168	! Allocate variables for pbl splitting
	169	!
	170	!****************************************************************************************
	171
	172	CALL wx_pbl_init
	173	!>jyg
	174
[781]	175	END SUBROUTINE pbl_surface_init
	176	!
	177	!****************************************************************************************
	178	!
	179
	180	SUBROUTINE pbl_surface( &
	181	dtime, date0, itap, jour, &
	182	debut, lafin, &
	183	rlon, rlat, rugoro, rmu0, &
[4013]	184	lwdown_m, cldt, &
[3798]	185	rain_f, snow_f, solsw_m, solswfdiff_m, sollw_m, &
[2240]	186	gustiness, &
[781]	187	t, q, u, v, &
[2159]	188	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
	189	!! t_x, q_x, t_w, q_w, &
	190	wake_dlt, wake_dlq, &
	191	wake_cstar, wake_s, &
	192	!!!
[781]	193	pplay, paprs, pctsrf, &
[2227]	194	ts,SFRWL, alb_dir, alb_dif,ustar, u10m, v10m,wstar, &
[2188]	195	cdragh, cdragm, zu1, zv1, &
[4013]	196	!jyg< (26/09/2019)
	197	beta, &
	198	!>jyg
[2227]	199	alb_dir_m, alb_dif_m, zxsens, zxevap, &
[1865]	200	alb3_lic, runoff, snowhgt, qsnow, to_ice, sissnow, &
[4013]	201	zxtsol, zxfluxlat, zt2m, qsat2m, zn2mout, &
[1761]	202	d_t, d_q, d_u, d_v, d_t_diss, &
[2159]	203	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
	204	d_t_w, d_q_w, &
	205	d_t_x, d_q_x, &
	206	!! d_wake_dlt,d_wake_dlq, &
	207	zxsens_x, zxfluxlat_x,zxsens_w,zxfluxlat_w, &
	208	!!!
	209	!!! nrlmd le 13/06/2011
	210	delta_tsurf,wake_dens,cdragh_x,cdragh_w, &
	211	cdragm_x,cdragm_w,kh,kh_x,kh_w, &
	212	!!!
[1539]	213	zcoefh, zcoefm, slab_wfbils, &
[2243]	214	qsol, zq2m, s_pblh, s_plcl, &
[2159]	215	!!!
	216	!!! jyg le 08/02/2012
	217	s_pblh_x, s_plcl_x, s_pblh_w, s_plcl_w, &
	218	!!!
[781]	219	s_capCL, s_oliqCL, s_cteiCL, s_pblT, &
	220	s_therm, s_trmb1, s_trmb2, s_trmb3, &
[2243]	221	zustar,zu10m, zv10m, fder_print, &
[4013]	222	zxqsurf, delta_qsurf, &
	223	rh2m, zxfluxu, zxfluxv, &
[2243]	224	z0m, z0h, agesno, sollw, solsw, &
	225	d_ts, evap, fluxlat, t2m, &
[2670]	226	wfbils, wfbilo, wfevap, wfrain, wfsnow, &
	227	flux_t, flux_u, flux_v, &
[781]	228	dflux_t, dflux_q, zxsnow, &
[2181]	229	!jyg<
	230	!! zxfluxt, zxfluxq, q2m, flux_q, tke, &
[3798]	231	zxfluxt, zxfluxq, q2m, flux_q, tke_x, &
[2181]	232	!>jyg
[2159]	233	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
	234	!! tke_x, tke_w &
[3179]	235	wake_dltke, &
	236	treedrg &
[2952]	237	!FC
[2159]	238	!!!
	239	)
[781]	240	!****************************************************************************************
	241	! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	242	! Objet: interface de "couche limite" (diffusion verticale)
	243	!
	244	!AA REM:
	245	!AA-----
	246	!AA Tout ce qui a trait au traceurs est dans phytrac maintenant
	247	!AA pour l'instant le calcul de la couche limite pour les traceurs
	248	!AA se fait avec cltrac et ne tient pas compte de la differentiation
	249	!AA des sous-fraction de sol.
	250	!AA REM bis :
	251	!AA----------
	252	!AA Pour pouvoir extraire les coefficient d'echanges et le vent
	253	!AA dans la premiere couche, 3 champs supplementaires ont ete crees
	254	!AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
	255	!AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
	256	!AA si les informations des subsurfaces doivent etre prises en compte
	257	!AA il faudra sortir ces memes champs en leur ajoutant une dimension,
	258	!AA c'est a dire nbsrf (nbre de subsurface).
	259	!
	260	! Arguments:
	261	!
	262	! dtime----input-R- interval du temps (secondes)
	263	! itap-----input-I- numero du pas de temps
	264	! date0----input-R- jour initial
	265	! t--------input-R- temperature (K)
	266	! q--------input-R- vapeur d'eau (kg/kg)
	267	! u--------input-R- vitesse u
	268	! v--------input-R- vitesse v
[2159]	269	! wake_dlt-input-R- temperatre difference between (w) and (x) (K)
	270	! wake_dlq-input-R- humidity difference between (w) and (x) (kg/kg)
	271	!wake_cstar-input-R- wake gust front speed (m/s)
	272	! wake_s---input-R- wake fractionnal area
[781]	273	! ts-------input-R- temperature du sol (en Kelvin)
	274	! paprs----input-R- pression a intercouche (Pa)
	275	! pplay----input-R- pression au milieu de couche (Pa)
	276	! rlat-----input-R- latitude en degree
[2243]	277	! z0m, z0h ----input-R- longeur de rugosite (en m)
[1865]	278	! Martin
	279	! cldt-----input-R- total cloud fraction
	280	! Martin
[781]	281	!
	282	! d_t------output-R- le changement pour "t"
	283	! d_q------output-R- le changement pour "q"
	284	! d_u------output-R- le changement pour "u"
	285	! d_v------output-R- le changement pour "v"
	286	! d_ts-----output-R- le changement pour "ts"
	287	! flux_t---output-R- flux de chaleur sensible (CpT) J/m2/s (W/m2)
	288	! (orientation positive vers le bas)
[2181]	289	! tke_x---input/output-R- tke in the (x) region (kg/m**2/s)
[2159]	290	! wake_dltke-input/output-R- tke difference between (w) and (x) (kg/m**2/s)
[781]	291	! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
	292	! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
	293	! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
	294	! dflux_t--output-R- derive du flux sensible
	295	! dflux_q--output-R- derive du flux latent
	296	! zu1------output-R- le vent dans la premiere couche
	297	! zv1------output-R- le vent dans la premiere couche
	298	! trmb1----output-R- deep_cape
	299	! trmb2----output-R- inhibition
	300	! trmb3----output-R- Point Omega
	301	! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL
	302	! plcl-----output-R- Niveau de condensation
	303	! pblh-----output-R- HCL
	304	! pblT-----output-R- T au nveau HCL
[3179]	305	! treedrg--output-R- tree drag (m)
[781]	306	!
[3391]	307	USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, carbon_cycle_tr, level_coupling_esm
	308	USE carbon_cycle_mod, ONLY : co2_send, nbcf_out, fields_out, yfields_out, cfname_out
[3798]	309	use hbtm_mod, only: hbtm
[1785]	310	USE indice_sol_mod
[3179]	311	USE time_phylmdz_mod, ONLY : day_ini,annee_ref,itau_phy
[3605]	312	USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat, grid1dto2d_glo
[3179]	313	USE print_control_mod, ONLY : prt_level,lunout
	314	USE ioipsl_getin_p_mod, ONLY : getin_p
[4020]	315	use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, dter, &
	316	dser, dt_ds, zsig, zmea
	317	use phys_output_var_mod, only: tkt, tks, taur, sss
[3626]	318	#ifdef CPP_XIOS
	319	USE wxios, ONLY: missing_val
	320	#else
	321	use netcdf, only: missing_val => nf90_fill_real
	322	#endif
[1785]	323
[4013]	324
	325
	326
[1279]	327	IMPLICIT NONE
	328
[781]	329	INCLUDE "dimsoil.h"
[793]	330	INCLUDE "YOMCST.h"
[1932]	331	INCLUDE "YOETHF.h"
[793]	332	INCLUDE "FCTTRE.h"
	333	INCLUDE "clesphys.h"
[781]	334	INCLUDE "compbl.h"
[1887]	335	INCLUDE "flux_arp.h"
[2952]	336	!FC
	337	INCLUDE "dimpft.h"
	338
[781]	339	!****************************************************************************************
[888]	340	REAL, INTENT(IN) :: dtime ! time interval (s)
	341	REAL, INTENT(IN) :: date0 ! initial day
	342	INTEGER, INTENT(IN) :: itap ! time step
	343	INTEGER, INTENT(IN) :: jour ! current day of the year
	344	LOGICAL, INTENT(IN) :: debut ! true if first run step
	345	LOGICAL, INTENT(IN) :: lafin ! true if last run step
	346	REAL, DIMENSION(klon), INTENT(IN) :: rlon ! longitudes in degrees
	347	REAL, DIMENSION(klon), INTENT(IN) :: rlat ! latitudes in degrees
	348	REAL, DIMENSION(klon), INTENT(IN) :: rugoro ! rugosity length
	349	REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cosine of solar zenith angle
	350	REAL, DIMENSION(klon), INTENT(IN) :: rain_f ! rain fall
	351	REAL, DIMENSION(klon), INTENT(IN) :: snow_f ! snow fall
	352	REAL, DIMENSION(klon), INTENT(IN) :: solsw_m ! net shortwave radiation at mean surface
[3798]	353	REAL, DIMENSION(klon), INTENT(IN) :: solswfdiff_m ! diffuse fraction fordownward shortwave radiation at mean surface
[888]	354	REAL, DIMENSION(klon), INTENT(IN) :: sollw_m ! net longwave radiation at mean surface
	355	REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K)
	356	REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! water vapour (kg/kg)
	357	REAL, DIMENSION(klon,klev), INTENT(IN) :: u ! u speed
	358	REAL, DIMENSION(klon,klev), INTENT(IN) :: v ! v speed
	359	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pression (Pa)
	360	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression between layers (Pa)
	361	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! sub-surface fraction
[1865]	362	! Martin
[2188]	363	REAL, DIMENSION(klon), INTENT(IN) :: lwdown_m ! downward longwave radiation at mean s
[2240]	364	REAL, DIMENSION(klon), INTENT(IN) :: gustiness ! gustiness
	365
[1865]	366	REAL, DIMENSION(klon), INTENT(IN) :: cldt ! total cloud fraction
[781]	367
[2159]	368	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
[2450]	369	!! REAL, DIMENSION(klon,klev), INTENT(IN) :: t_x ! Temp\'erature hors poche froide
	370	!! REAL, DIMENSION(klon,klev), INTENT(IN) :: t_w ! Temp\'erature dans la poches froide
[2159]	371	!! REAL, DIMENSION(klon,klev), INTENT(IN) :: q_x !
[2450]	372	!! REAL, DIMENSION(klon,klev), INTENT(IN) :: q_w ! Pareil pour l'humidit\'e
[2159]	373	REAL, DIMENSION(klon,klev), INTENT(IN) :: wake_dlt !temperature difference between (w) and (x) (K)
	374	REAL, DIMENSION(klon,klev), INTENT(IN) :: wake_dlq !humidity difference between (w) and (x) (K)
	375	REAL, DIMENSION(klon), INTENT(IN) :: wake_s ! Fraction de poches froides
	376	REAL, DIMENSION(klon), INTENT(IN) :: wake_cstar! Vitesse d'expansion des poches froides
	377	REAL, DIMENSION(klon), INTENT(IN) :: wake_dens
	378	!!!
	379
[781]	380	! Input/Output variables
	381	!****************************************************************************************
[4013]	382	!jyg<
	383	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: beta ! Aridity factor
	384	!>jyg
[888]	385	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ts ! temperature at surface (K)
[2159]	386	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: delta_tsurf !surface temperature difference between
	387	!wake and off-wake regions
[2227]	388	!albedo SB >>>
	389	REAL, DIMENSIOn(6),intent(in) :: SFRWL
	390	REAL, DIMENSION(klon, nsw, nbsrf), INTENT(INOUT) :: alb_dir,alb_dif
	391	!albedo SB <<<
[2181]	392	!jyg Pourquoi ustar et wstar sont-elles INOUT ?
[1670]	393	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ustar ! u* (m/s)
[1816]	394	REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT) :: wstar ! w* (m/s)
[888]	395	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: u10m ! u speed at 10m
	396	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: v10m ! v speed at 10m
[2181]	397	!jyg<
	398	!! REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke
	399	REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke_x
	400	!>jyg
[2159]	401
	402	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
[2181]	403	REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: wake_dltke ! TKE_w - TKE_x
[2159]	404	!!!
	405
[781]	406	! Output variables
	407	!****************************************************************************************
[888]	408	REAL, DIMENSION(klon), INTENT(OUT) :: cdragh ! drag coefficient for T and Q
	409	REAL, DIMENSION(klon), INTENT(OUT) :: cdragm ! drag coefficient for wind
	410	REAL, DIMENSION(klon), INTENT(OUT) :: zu1 ! u wind speed in first layer
	411	REAL, DIMENSION(klon), INTENT(OUT) :: zv1 ! v wind speed in first layer
[2227]	412	!albedo SB >>>
[3391]	413	REAL, DIMENSION(klon, nsw), INTENT(OUT) :: alb_dir_m,alb_dif_m
[2227]	414	!albedo SB <<<
[1865]	415	! Martin
[2126]	416	REAL, DIMENSION(klon), INTENT(OUT) :: alb3_lic
[1865]	417	! Martin
[888]	418	REAL, DIMENSION(klon), INTENT(OUT) :: zxsens ! sensible heat flux at surface with inversed sign
	419	! (=> positive sign upwards)
	420	REAL, DIMENSION(klon), INTENT(OUT) :: zxevap ! water vapour flux at surface, positiv upwards
	421	REAL, DIMENSION(klon), INTENT(OUT) :: zxtsol ! temperature at surface, mean for each grid point
[2159]	422	!!! jyg le ???
	423	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t_w ! !
	424	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q_w ! ! Tendances dans les poches
	425	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t_x ! !
	426	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q_x ! ! Tendances hors des poches
	427	!!! jyg
[888]	428	REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat ! latent flux, mean for each grid point
	429	REAL, DIMENSION(klon), INTENT(OUT) :: zt2m ! temperature at 2m, mean for each grid point
[4013]	430	INTEGER, DIMENSION(klon, 6), INTENT(OUT) :: zn2mout ! number of times the 2m temperature is out of the [tsol,temp]
[781]	431	REAL, DIMENSION(klon), INTENT(OUT) :: qsat2m
[888]	432	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! change in temperature
[1761]	433	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t_diss ! change in temperature
[888]	434	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_q ! change in water vapour
	435	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u ! change in u speed
	436	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v ! change in v speed
[781]	437
[1919]	438	REAL, INTENT(OUT):: zcoefh(:, :, :) ! (klon, klev, nbsrf + 1)
	439	! coef for turbulent diffusion of T and Q, mean for each grid point
	440
	441	REAL, INTENT(OUT):: zcoefm(:, :, :) ! (klon, klev, nbsrf + 1)
	442	! coef for turbulent diffusion of U and V (?), mean for each grid point
	443
[2159]	444	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
	445	REAL, DIMENSION(klon), INTENT(OUT) :: zxsens_x ! Flux sensible hors poche
	446	REAL, DIMENSION(klon), INTENT(OUT) :: zxsens_w ! Flux sensible dans la poche
	447	REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat_x! Flux latent hors poche
	448	REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat_w! Flux latent dans la poche
	449	!! REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_wake_dlt
	450	!! REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_wake_dlq
	451
[781]	452	! Output only for diagnostics
[2159]	453	REAL, DIMENSION(klon), INTENT(OUT) :: cdragh_x
	454	REAL, DIMENSION(klon), INTENT(OUT) :: cdragh_w
	455	REAL, DIMENSION(klon), INTENT(OUT) :: cdragm_x
	456	REAL, DIMENSION(klon), INTENT(OUT) :: cdragm_w
	457	REAL, DIMENSION(klon), INTENT(OUT) :: kh
	458	REAL, DIMENSION(klon), INTENT(OUT) :: kh_x
	459	REAL, DIMENSION(klon), INTENT(OUT) :: kh_w
	460	!!!
[996]	461	REAL, DIMENSION(klon), INTENT(OUT) :: slab_wfbils! heat balance at surface only for slab at ocean points
[2243]	462	REAL, DIMENSION(klon), INTENT(OUT) :: qsol ! water height in the soil (mm)
[888]	463	REAL, DIMENSION(klon), INTENT(OUT) :: zq2m ! water vapour at 2m, mean for each grid point
	464	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh ! height of the planetary boundary layer(HPBL)
[2159]	465	!!! jyg le 08/02/2012
	466	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh_x ! height of the PBL in the off-wake region
	467	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh_w ! height of the PBL in the wake region
	468	!!!
[888]	469	REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl ! condensation level
[2159]	470	!!! jyg le 08/02/2012
	471	REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl_x ! condensation level in the off-wake region
	472	REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl_w ! condensation level in the wake region
	473	!!!
[888]	474	REAL, DIMENSION(klon), INTENT(OUT) :: s_capCL ! CAPE of PBL
	475	REAL, DIMENSION(klon), INTENT(OUT) :: s_oliqCL ! liquid water intergral of PBL
	476	REAL, DIMENSION(klon), INTENT(OUT) :: s_cteiCL ! cloud top instab. crit. of PBL
	477	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblT ! temperature at PBLH
	478	REAL, DIMENSION(klon), INTENT(OUT) :: s_therm ! thermal virtual temperature excess
	479	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb1 ! deep cape, mean for each grid point
	480	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb2 ! inhibition, mean for each grid point
	481	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb3 ! point Omega, mean for each grid point
[1670]	482	REAL, DIMENSION(klon), INTENT(OUT) :: zustar ! u*
[888]	483	REAL, DIMENSION(klon), INTENT(OUT) :: zu10m ! u speed at 10m, mean for each grid point
	484	REAL, DIMENSION(klon), INTENT(OUT) :: zv10m ! v speed at 10m, mean for each grid point
	485	REAL, DIMENSION(klon), INTENT(OUT) :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
	486	REAL, DIMENSION(klon), INTENT(OUT) :: zxqsurf ! humidity at surface, mean for each grid point
[4013]	487	REAL, DIMENSION(klon), INTENT(OUT) :: delta_qsurf! humidity difference at surface, mean for each grid point
[888]	488	REAL, DIMENSION(klon), INTENT(OUT) :: rh2m ! relative humidity at 2m
	489	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxu ! u wind tension, mean for each grid point
	490	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxv ! v wind tension, mean for each grid point
[3798]	491	REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT) :: z0m,z0h ! rugosity length (m)
	492	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: agesno ! age of snow at surface
[888]	493	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: solsw ! net shortwave radiation at surface
	494	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: sollw ! net longwave radiation at surface
	495	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: d_ts ! change in temperature at surface
[3798]	496	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: evap ! evaporation at surface
[888]	497	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: fluxlat ! latent flux
	498	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: t2m ! temperature at 2 meter height
	499	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbils ! heat balance at surface
[2670]	500	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbilo ! water balance at surface weighted by srf
	501	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfevap ! water balance (evap) at surface weighted by srf
	502	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfrain ! water balance (rain) at surface weighted by srf
	503	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfsnow ! water balance (snow) at surface weighted by srf
[888]	504	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t ! sensible heat flux (CpT) J/m2/s (W/m2)
	505	! positve orientation downwards
	506	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u ! u wind tension (kg m/s)/(m**2 s) or Pascal
	507	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v ! v wind tension (kg m/s)/(m**2 s) or Pascal
[2952]	508	!FC
	509	REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: treedrg ! tree drag (m)
[781]	510
[2952]	511
[781]	512	! Output not needed
[888]	513	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_t ! change of sensible heat flux
	514	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_q ! change of water vapour flux
	515	REAL, DIMENSION(klon), INTENT(OUT) :: zxsnow ! snow at surface, mean for each grid point
	516	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxt ! sensible heat flux, mean for each grid point
	517	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxq ! water vapour flux, mean for each grid point
	518	REAL, DIMENSION(klon, nbsrf),INTENT(OUT) :: q2m ! water vapour at 2 meter height
	519	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q ! water vapour flux(latent flux) (kg/m**2/s)
[781]	520
[1865]	521	! Martin
[4013]	522	! inlandsis
[1865]	523	REAL, DIMENSION(klon), INTENT(OUT) :: qsnow ! snow water content
	524	REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt ! snow height
	525	REAL, DIMENSION(klon), INTENT(OUT) :: to_ice ! snow passed to ice
	526	REAL, DIMENSION(klon), INTENT(OUT) :: sissnow ! snow in snow model
	527	REAL, DIMENSION(klon), INTENT(OUT) :: runoff ! runoff on land ice
	528	! Martin
[781]	529
	530	! Local variables with attribute SAVE
	531	!****************************************************************************************
[888]	532	INTEGER, SAVE :: nhoridbg, nidbg ! variables for IOIPSL
[781]	533	!$OMP THREADPRIVATE(nhoridbg, nidbg)
	534	LOGICAL, SAVE :: debugindex=.FALSE.
	535	!$OMP THREADPRIVATE(debugindex)
	536	LOGICAL, SAVE :: first_call=.TRUE.
	537	!$OMP THREADPRIVATE(first_call)
	538	CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
	539	!$OMP THREADPRIVATE(cl_surf)
[3179]	540	REAL, SAVE :: beta_land ! beta for wx_dts
	541	!$OMP THREADPRIVATE(beta_land)
[781]	542
	543	! Other local variables
	544	!****************************************************************************************
[3391]	545	! >> PC
	546	INTEGER :: ierr
	547	INTEGER :: n
	548	! << PC
[4013]	549	INTEGER :: iflag_split, iflag_split_ref
[781]	550	INTEGER :: i, k, nsrf
	551	INTEGER :: knon, j
	552	INTEGER :: idayref
	553	INTEGER , DIMENSION(klon) :: ni
[2159]	554	REAL :: yt1_new
[781]	555	REAL :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
	556	REAL :: amn, amx
[888]	557	REAL :: f1 ! fraction de longeurs visibles parmi tout SW intervalle
[781]	558	REAL, DIMENSION(klon) :: r_co2_ppm ! taux CO2 atmosphere
[2243]	559	REAL, DIMENSION(klon) :: yts, yz0m, yz0h, ypct
[4013]	560	REAL, DIMENSION(klon) :: yz0h_old
[2227]	561	!albedo SB >>>
	562	REAL, DIMENSION(klon) :: yalb,yalb_vis
	563	!albedo SB <<<
[3179]	564	REAL, DIMENSION(klon) :: yt1, yq1, yu1, yv1
[4013]	565	REAL, DIMENSION(klon) :: yqa
[781]	566	REAL, DIMENSION(klon) :: ysnow, yqsurf, yagesno, yqsol
	567	REAL, DIMENSION(klon) :: yrain_f, ysnow_f
[888]	568	REAL, DIMENSION(klon) :: ysolsw, ysollw
[781]	569	REAL, DIMENSION(klon) :: yfder
[888]	570	REAL, DIMENSION(klon) :: yrugoro
[781]	571	REAL, DIMENSION(klon) :: yfluxlat
	572	REAL, DIMENSION(klon) :: y_d_ts
	573	REAL, DIMENSION(klon) :: y_flux_t1, y_flux_q1
	574	REAL, DIMENSION(klon) :: y_dflux_t, y_dflux_q
[1067]	575	REAL, DIMENSION(klon) :: y_flux_u1, y_flux_v1
[781]	576	REAL, DIMENSION(klon) :: yt2m, yq2m, yu10m
[4013]	577	INTEGER, DIMENSION(klon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w
	578	INTEGER, DIMENSION(klon, nbsrf, 6) :: n2mout, n2mout_x, n2mout_w
[781]	579	REAL, DIMENSION(klon) :: yustar
[1816]	580	REAL, DIMENSION(klon) :: ywstar
[781]	581	REAL, DIMENSION(klon) :: ywindsp
	582	REAL, DIMENSION(klon) :: yt10m, yq10m
	583	REAL, DIMENSION(klon) :: ypblh
	584	REAL, DIMENSION(klon) :: ylcl
	585	REAL, DIMENSION(klon) :: ycapCL
	586	REAL, DIMENSION(klon) :: yoliqCL
	587	REAL, DIMENSION(klon) :: ycteiCL
	588	REAL, DIMENSION(klon) :: ypblT
	589	REAL, DIMENSION(klon) :: ytherm
	590	REAL, DIMENSION(klon) :: ytrmb1
	591	REAL, DIMENSION(klon) :: ytrmb2
	592	REAL, DIMENSION(klon) :: ytrmb3
	593	REAL, DIMENSION(klon) :: uzon, vmer
	594	REAL, DIMENSION(klon) :: tair1, qair1, tairsol
	595	REAL, DIMENSION(klon) :: psfce, patm
[2232]	596	REAL, DIMENSION(klon) :: qairsol, zgeo1, speed, zri1, pref !speed, zri1, pref, added by Fuxing WANG, 04/03/2015
[2896]	597	REAL, DIMENSION(klon) :: yz0h_oupas
[888]	598	REAL, DIMENSION(klon) :: yfluxsens
[4013]	599	REAL, DIMENSION(klon) :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0
[1067]	600	REAL, DIMENSION(klon) :: AcoefH, AcoefQ, BcoefH, BcoefQ
	601	REAL, DIMENSION(klon) :: AcoefU, AcoefV, BcoefU, BcoefV
[888]	602	REAL, DIMENSION(klon) :: ypsref
[4013]	603	REAL, DIMENSION(klon) :: yevap, yevap_pot, ytsurf_new, yalb3_new
[2227]	604	!albedo SB >>>
	605	REAL, DIMENSION(klon,nsw) :: yalb_dir_new, yalb_dif_new
	606	!albedo SB <<<
[781]	607	REAL, DIMENSION(klon) :: ztsol
[2455]	608	REAL, DIMENSION(klon) :: meansqT ! mean square deviation of subsurface temperatures
[888]	609	REAL, DIMENSION(klon) :: alb_m ! mean albedo for whole SW interval
[1761]	610	REAL, DIMENSION(klon,klev) :: y_d_t, y_d_q, y_d_t_diss
[781]	611	REAL, DIMENSION(klon,klev) :: y_d_u, y_d_v
	612	REAL, DIMENSION(klon,klev) :: y_flux_t, y_flux_q
	613	REAL, DIMENSION(klon,klev) :: y_flux_u, y_flux_v
[1761]	614	REAL, DIMENSION(klon,klev) :: ycoefh, ycoefm,ycoefq
[4013]	615	REAL, DIMENSION(klon) :: ycdragh, ycdragq, ycdragm
[781]	616	REAL, DIMENSION(klon,klev) :: yu, yv
	617	REAL, DIMENSION(klon,klev) :: yt, yq
	618	REAL, DIMENSION(klon,klev) :: ypplay, ydelp
	619	REAL, DIMENSION(klon,klev) :: delp
	620	REAL, DIMENSION(klon,klev+1) :: ypaprs
[878]	621	REAL, DIMENSION(klon,klev+1) :: ytke
[781]	622	REAL, DIMENSION(klon,nsoilmx) :: ytsoil
[2952]	623	!FC
	624	REAL, DIMENSION(klon,nvm_lmdz) :: yveget
	625	REAL, DIMENSION(klon,nvm_lmdz) :: ylai
	626	REAL, DIMENSION(klon,nvm_lmdz) :: yheight
	627	REAL, DIMENSION(klon,klev) :: y_d_u_frein
	628	REAL, DIMENSION(klon,klev) :: y_d_v_frein
	629	REAL, DIMENSION(klon,klev) :: y_treedrg
	630	!FC
	631
[781]	632	CHARACTER(len=80) :: abort_message
	633	CHARACTER(len=20) :: modname = 'pbl_surface'
	634	LOGICAL, PARAMETER :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
	635	LOGICAL, PARAMETER :: check=.FALSE.
[2159]	636
	637	!!! nrlmd le 02/05/2011
	638	!!! jyg le 07/02/2012
	639	REAL, DIMENSION(klon) :: ywake_s, ywake_cstar, ywake_dens
	640	!!!
	641	REAL, DIMENSION(klon,klev+1) :: ytke_x, ytke_w
	642	REAL, DIMENSION(klon,klev+1) :: ywake_dltke
	643	REAL, DIMENSION(klon,klev) :: yu_x, yv_x, yu_w, yv_w
	644	REAL, DIMENSION(klon,klev) :: yt_x, yq_x, yt_w, yq_w
	645	REAL, DIMENSION(klon,klev) :: ycoefh_x, ycoefm_x, ycoefh_w, ycoefm_w
	646	REAL, DIMENSION(klon,klev) :: ycoefq_x, ycoefq_w
[4013]	647	REAL, DIMENSION(klon) :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w
	648	REAL, DIMENSION(klon) :: ycdragm_x, ycdragm_w
[2159]	649	REAL, DIMENSION(klon) :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x
	650	REAL, DIMENSION(klon) :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w
	651	REAL, DIMENSION(klon) :: AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x
	652	REAL, DIMENSION(klon) :: AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w
	653	REAL, DIMENSION(klon) :: y_flux_t1_x, y_flux_q1_x, y_flux_t1_w, y_flux_q1_w
	654	REAL, DIMENSION(klon) :: y_flux_u1_x, y_flux_v1_x, y_flux_u1_w, y_flux_v1_w
	655	REAL, DIMENSION(klon,klev) :: y_flux_t_x, y_flux_q_x, y_flux_t_w, y_flux_q_w
	656	REAL, DIMENSION(klon,klev) :: y_flux_u_x, y_flux_v_x, y_flux_u_w, y_flux_v_w
	657	REAL, DIMENSION(klon) :: yfluxlat_x, yfluxlat_w
	658	REAL, DIMENSION(klon,klev) :: y_d_t_x, y_d_q_x, y_d_t_w, y_d_q_w
	659	REAL, DIMENSION(klon,klev) :: y_d_t_diss_x, y_d_t_diss_w
	660	REAL, DIMENSION(klon,klev) :: d_t_diss_x, d_t_diss_w
	661	REAL, DIMENSION(klon,klev) :: y_d_u_x, y_d_v_x, y_d_u_w, y_d_v_w
	662	REAL, DIMENSION(klon, klev, nbsrf) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w
	663	REAL, DIMENSION(klon, klev, nbsrf) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w
	664	REAL, DIMENSION(klon, nbsrf) :: fluxlat_x, fluxlat_w
	665	REAL, DIMENSION(klon, klev) :: zxfluxt_x, zxfluxq_x, zxfluxt_w, zxfluxq_w
	666	REAL, DIMENSION(klon, klev) :: zxfluxu_x, zxfluxv_x, zxfluxu_w, zxfluxv_w
	667	REAL :: zx_qs_surf, zcor_surf, zdelta_surf
[4013]	668	!jyg<
[2159]	669	REAL, DIMENSION(klon) :: ybeta
[4013]	670	REAL, DIMENSION(klon) :: ybeta_prev
	671	!>jyg
[2159]	672	REAL, DIMENSION(klon, klev) :: d_u_x
	673	REAL, DIMENSION(klon, klev) :: d_u_w
	674	REAL, DIMENSION(klon, klev) :: d_v_x
	675	REAL, DIMENSION(klon, klev) :: d_v_w
	676
	677	REAL, DIMENSION(klon,klev) :: CcoefH, CcoefQ, DcoefH, DcoefQ
	678	REAL, DIMENSION(klon,klev) :: CcoefU, CcoefV, DcoefU, DcoefV
	679	REAL, DIMENSION(klon,klev) :: CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x
	680	REAL, DIMENSION(klon,klev) :: CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w
	681	REAL, DIMENSION(klon,klev) :: CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x
	682	REAL, DIMENSION(klon,klev) :: CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w
	683	REAL, DIMENSION(klon,klev) :: Kcoef_hq, Kcoef_m, gama_h, gama_q
	684	REAL, DIMENSION(klon,klev) :: Kcoef_hq_x, Kcoef_m_x, gama_h_x, gama_q_x
	685	REAL, DIMENSION(klon,klev) :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w
	686	REAL, DIMENSION(klon) :: alf_1, alf_2, alf_1_x, alf_2_x, alf_1_w, alf_2_w
	687	!!!
	688	!!!jyg le 08/02/2012
[2181]	689	REAL, DIMENSION(klon, nbsrf) :: windsp
	690	!
[2159]	691	REAL, DIMENSION(klon, nbsrf) :: t2m_x
	692	REAL, DIMENSION(klon, nbsrf) :: q2m_x
	693	REAL, DIMENSION(klon) :: rh2m_x
	694	REAL, DIMENSION(klon) :: qsat2m_x
	695	REAL, DIMENSION(klon, nbsrf) :: u10m_x
	696	REAL, DIMENSION(klon, nbsrf) :: v10m_x
	697	REAL, DIMENSION(klon, nbsrf) :: ustar_x
	698	REAL, DIMENSION(klon, nbsrf) :: wstar_x
	699	!
	700	REAL, DIMENSION(klon, nbsrf) :: pblh_x
	701	REAL, DIMENSION(klon, nbsrf) :: plcl_x
	702	REAL, DIMENSION(klon, nbsrf) :: capCL_x
	703	REAL, DIMENSION(klon, nbsrf) :: oliqCL_x
	704	REAL, DIMENSION(klon, nbsrf) :: cteiCL_x
	705	REAL, DIMENSION(klon, nbsrf) :: pblt_x
	706	REAL, DIMENSION(klon, nbsrf) :: therm_x
	707	REAL, DIMENSION(klon, nbsrf) :: trmb1_x
	708	REAL, DIMENSION(klon, nbsrf) :: trmb2_x
	709	REAL, DIMENSION(klon, nbsrf) :: trmb3_x
	710	!
	711	REAL, DIMENSION(klon, nbsrf) :: t2m_w
	712	REAL, DIMENSION(klon, nbsrf) :: q2m_w
	713	REAL, DIMENSION(klon) :: rh2m_w
	714	REAL, DIMENSION(klon) :: qsat2m_w
	715	REAL, DIMENSION(klon, nbsrf) :: u10m_w
	716	REAL, DIMENSION(klon, nbsrf) :: v10m_w
	717	REAL, DIMENSION(klon, nbsrf) :: ustar_w
	718	REAL, DIMENSION(klon, nbsrf) :: wstar_w
	719	!
	720	REAL, DIMENSION(klon, nbsrf) :: pblh_w
	721	REAL, DIMENSION(klon, nbsrf) :: plcl_w
	722	REAL, DIMENSION(klon, nbsrf) :: capCL_w
	723	REAL, DIMENSION(klon, nbsrf) :: oliqCL_w
	724	REAL, DIMENSION(klon, nbsrf) :: cteiCL_w
	725	REAL, DIMENSION(klon, nbsrf) :: pblt_w
	726	REAL, DIMENSION(klon, nbsrf) :: therm_w
	727	REAL, DIMENSION(klon, nbsrf) :: trmb1_w
	728	REAL, DIMENSION(klon, nbsrf) :: trmb2_w
	729	REAL, DIMENSION(klon, nbsrf) :: trmb3_w
	730	!
	731	REAL, DIMENSION(klon) :: yt2m_x
	732	REAL, DIMENSION(klon) :: yq2m_x
	733	REAL, DIMENSION(klon) :: yt10m_x
	734	REAL, DIMENSION(klon) :: yq10m_x
	735	REAL, DIMENSION(klon) :: yu10m_x
	736	REAL, DIMENSION(klon) :: yv10m_x
	737	REAL, DIMENSION(klon) :: yustar_x
	738	REAL, DIMENSION(klon) :: ywstar_x
	739	!
	740	REAL, DIMENSION(klon) :: ypblh_x
	741	REAL, DIMENSION(klon) :: ylcl_x
	742	REAL, DIMENSION(klon) :: ycapCL_x
	743	REAL, DIMENSION(klon) :: yoliqCL_x
	744	REAL, DIMENSION(klon) :: ycteiCL_x
	745	REAL, DIMENSION(klon) :: ypblt_x
	746	REAL, DIMENSION(klon) :: ytherm_x
	747	REAL, DIMENSION(klon) :: ytrmb1_x
	748	REAL, DIMENSION(klon) :: ytrmb2_x
	749	REAL, DIMENSION(klon) :: ytrmb3_x
	750	!
	751	REAL, DIMENSION(klon) :: yt2m_w
	752	REAL, DIMENSION(klon) :: yq2m_w
	753	REAL, DIMENSION(klon) :: yt10m_w
	754	REAL, DIMENSION(klon) :: yq10m_w
	755	REAL, DIMENSION(klon) :: yu10m_w
	756	REAL, DIMENSION(klon) :: yv10m_w
	757	REAL, DIMENSION(klon) :: yustar_w
	758	REAL, DIMENSION(klon) :: ywstar_w
	759	!
	760	REAL, DIMENSION(klon) :: ypblh_w
	761	REAL, DIMENSION(klon) :: ylcl_w
	762	REAL, DIMENSION(klon) :: ycapCL_w
	763	REAL, DIMENSION(klon) :: yoliqCL_w
	764	REAL, DIMENSION(klon) :: ycteiCL_w
	765	REAL, DIMENSION(klon) :: ypblt_w
	766	REAL, DIMENSION(klon) :: ytherm_w
	767	REAL, DIMENSION(klon) :: ytrmb1_w
	768	REAL, DIMENSION(klon) :: ytrmb2_w
	769	REAL, DIMENSION(klon) :: ytrmb3_w
	770	!
[3391]	771	REAL, DIMENSION(klon) :: uzon_x, vmer_x, speed_x, zri1_x, pref_x !speed_x, zri1_x, pref_x, added by Fuxing WANG, 04/03/2015
	772	REAL, DIMENSION(klon) :: zgeo1_x, tair1_x, qair1_x, tairsol_x
[2159]	773	!
[3391]	774	REAL, DIMENSION(klon) :: uzon_w, vmer_w, speed_w, zri1_w, pref_w !speed_w, zri1_w, pref_w, added by Fuxing WANG, 04/03/2015
	775	REAL, DIMENSION(klon) :: zgeo1_w, tair1_w, qair1_w, tairsol_w
[2159]	776
	777	!!! jyg le 25/03/2013
[2450]	778	!! Variables intermediaires pour le raccord des deux colonnes \`a la surface
[3179]	779	!jyg<
	780	!! REAL :: dd_Ch
	781	!! REAL :: dd_Cm
	782	!! REAL :: dd_Kh
	783	!! REAL :: dd_Km
	784	!! REAL :: dd_u
	785	!! REAL :: dd_v
	786	!! REAL :: dd_t
	787	!! REAL :: dd_q
	788	!! REAL :: dd_AH
	789	!! REAL :: dd_AQ
	790	!! REAL :: dd_AU
	791	!! REAL :: dd_AV
	792	!! REAL :: dd_BH
	793	!! REAL :: dd_BQ
	794	!! REAL :: dd_BU
	795	!! REAL :: dd_BV
	796	!!
	797	!! REAL :: dd_KHp
	798	!! REAL :: dd_KQp
	799	!! REAL :: dd_KUp
	800	!! REAL :: dd_KVp
	801	!>jyg
[2159]	802
	803	!!!
	804	!!! nrlmd le 13/06/2011
	805	REAL, DIMENSION(klon) :: y_delta_flux_t1, y_delta_flux_q1, y_delta_flux_u1, y_delta_flux_v1
[4013]	806	REAL, DIMENSION(klon) :: y_delta_tsurf, y_delta_tsurf_new
	807	REAL, DIMENSION(klon) :: delta_coef, tau_eq
	808	REAL, DIMENSION(klon) :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
	809	REAL, DIMENSION(klon) :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0
	810	REAL, DIMENSION(klon) :: y_delta_qsurf
	811	REAL, DIMENSION(klon) :: y_delta_qsats
	812	REAL, DIMENSION(klon) :: yg_T, yg_Q
	813	REAL, DIMENSION(klon) :: yGamma_dTs_phiT, yGamma_dQs_phiQ
	814	REAL, DIMENSION(klon) :: ydTs_ins, ydqs_ins
	815	!
[2159]	816	REAL, PARAMETER :: facteur=2./sqrt(3.14)
[3179]	817	REAL, PARAMETER :: inertia=2000.
[2159]	818	REAL, DIMENSION(klon) :: ydtsurf_th
	819	REAL :: zdelta_surf_x,zdelta_surf_w,zx_qs_surf_x,zx_qs_surf_w
	820	REAL :: zcor_surf_x,zcor_surf_w
	821	REAL :: mod_wind_x, mod_wind_w
	822	REAL :: rho1
	823	REAL, DIMENSION(klon) :: Kech_h ! Coefficient d'echange pour l'energie
	824	REAL, DIMENSION(klon) :: Kech_h_x, Kech_h_w
	825	REAL, DIMENSION(klon) :: Kech_m
	826	REAL, DIMENSION(klon) :: Kech_m_x, Kech_m_w
[4013]	827	REAL, DIMENSION(klon) :: yts_x, yts_w
	828	REAL, DIMENSION(klon) :: yqsatsrf0_x, yqsatsrf0_w
	829	REAL, DIMENSION(klon) :: yqsurf_x, yqsurf_w
[3179]	830	!jyg<
	831	!! REAL, DIMENSION(klon) :: Kech_Hp, Kech_H_xp, Kech_H_wp
	832	!! REAL, DIMENSION(klon) :: Kech_Qp, Kech_Q_xp, Kech_Q_wp
	833	!! REAL, DIMENSION(klon) :: Kech_Up, Kech_U_xp, Kech_U_wp
	834	!! REAL, DIMENSION(klon) :: Kech_Vp, Kech_V_xp, Kech_V_wp
	835	!>jyg
[2159]	836
[4013]	837	REAL :: fact_cdrag
	838	REAL :: z1lay
	839
[2126]	840	REAL :: vent
[3391]	841	!
[781]	842	! For debugging with IOIPSL
[2344]	843	INTEGER, DIMENSION(nbp_lon*nbp_lat) :: ndexbg
[781]	844	REAL :: zjulian
	845	REAL, DIMENSION(klon) :: tabindx
[2344]	846	REAL, DIMENSION(nbp_lon,nbp_lat) :: zx_lon, zx_lat
	847	REAL, DIMENSION(nbp_lon,nbp_lat) :: debugtab
[781]	848
	849
[888]	850	REAL, DIMENSION(klon,nbsrf) :: pblh ! height of the planetary boundary layer
	851	REAL, DIMENSION(klon,nbsrf) :: plcl ! condensation level
[781]	852	REAL, DIMENSION(klon,nbsrf) :: capCL
	853	REAL, DIMENSION(klon,nbsrf) :: oliqCL
	854	REAL, DIMENSION(klon,nbsrf) :: cteiCL
	855	REAL, DIMENSION(klon,nbsrf) :: pblT
	856	REAL, DIMENSION(klon,nbsrf) :: therm
[888]	857	REAL, DIMENSION(klon,nbsrf) :: trmb1 ! deep cape
	858	REAL, DIMENSION(klon,nbsrf) :: trmb2 ! inhibition
	859	REAL, DIMENSION(klon,nbsrf) :: trmb3 ! point Omega
[1067]	860	REAL, DIMENSION(klon,nbsrf) :: zx_rh2m, zx_qsat2m
[996]	861	REAL, DIMENSION(klon,nbsrf) :: zx_t1
[888]	862	REAL, DIMENSION(klon, nbsrf) :: alb ! mean albedo for whole SW interval
	863	REAL, DIMENSION(klon) :: ylwdown ! jg : temporary (ysollwdown)
[2240]	864	REAL, DIMENSION(klon) :: ygustiness ! jg : temporary (ysollwdown)
[781]	865
[996]	866	REAL :: zx_qs1, zcor1, zdelta1
[781]	867
[1865]	868	! Martin
	869	REAL, DIMENSION(klon, nbsrf) :: sollwd ! net longwave radiation at surface
	870	REAL, DIMENSION(klon) :: ytoice
	871	REAL, DIMENSION(klon) :: ysnowhgt, yqsnow, ysissnow, yrunoff
[4013]	872	REAL, DIMENSION(klon) :: yzmea
[1865]	873	REAL, DIMENSION(klon) :: yzsig
	874	REAL, DIMENSION(klon) :: ycldt
	875	REAL, DIMENSION(klon) :: yrmu0
	876	! Martin
	877
[4020]	878	REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, &
	879	ydser, ydt_ds, ytkt, ytks, ytaur, ysss
	880	! compression of delta_sst, delta_sal, ds_ns, dt_ns, dter, dser,
	881	! dt_ds, tkt, tks, taur, sss on ocean points
[3429]	882
[781]	883	!****************************************************************************************
	884	! End of declarations
	885	!****************************************************************************************
	886
[2159]	887	IF (prt_level >=10) print *,' -> pbl_surface, itap ',itap
	888	!
[2852]	889	!!jyg iflag_split = mod(iflag_pbl_split,2)
[4013]	890	!!jyg iflag_split = mod(iflag_pbl_split,10)
	891	!
	892	! Flags controlling the splitting of the turbulent boundary layer:
	893	! iflag_split_ref = 0 ==> no splitting
	894	! = 1 ==> splitting without coupling with surface temperature
	895	! = 2 ==> splitting with coupling with surface temperature over land
	896	! = 3 ==> splitting over ocean; no splitting over land
	897	! iflag_split: actual flag controlling the splitting.
	898	! iflag_split = iflag_split_ref outside the sub-surface loop
	899	! = iflag_split_ref if iflag_split_ref = 0, 1, or 2
	900	! = 0 over land if iflga_split_ref = 3
	901	! = 1 over ocean if iflga_split_ref = 3
[781]	902
[4013]	903	iflag_split_ref = mod(iflag_pbl_split,10)
	904	iflag_split = iflag_split_ref
	905
[781]	906	!****************************************************************************************
	907	! 1) Initialisation and validation tests
	908	! Only done first time entering this subroutine
	909	!
	910	!****************************************************************************************
	911
	912	IF (first_call) THEN
[4013]	913
	914	iflag_new_t2mq2m=1
	915	CALL getin_p('iflag_new_t2mq2m',iflag_new_t2mq2m)
	916	WRITE(lunout,*) 'pbl_iflag_new_t2mq2m=',iflag_new_t2mq2m
	917
[2240]	918	print*,'PBL SURFACE AVEC GUSTINESS'
[781]	919	first_call=.FALSE.
	920
[1282]	921	! Initialize ok_flux_surf (for 1D model)
[4013]	922	IF (klon_glo>1) ok_flux_surf=.FALSE.
	923	IF (klon_glo>1) ok_forc_tsurf=.FALSE.
[3179]	924
	925	! intialize beta_land
	926	beta_land = 0.5
	927	call getin_p('beta_land', beta_land)
[1282]	928
[781]	929	! Initilize debug IO
	930	IF (debugindex .AND. mpi_size==1) THEN
	931	! initialize IOIPSL output
	932	idayref = day_ini
	933	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
[3605]	934	CALL grid1dTo2d_glo(rlon,zx_lon)
[2344]	935	DO i = 1, nbp_lon
[781]	936	zx_lon(i,1) = rlon(i+1)
[2344]	937	zx_lon(i,nbp_lat) = rlon(i+1)
[781]	938	ENDDO
[3605]	939	CALL grid1dTo2d_glo(rlat,zx_lat)
[2344]	940	CALL histbeg("sous_index",nbp_lon,zx_lon(:,1),nbp_lat,zx_lat(1,:), &
	941	1,nbp_lon,1,nbp_lat, &
[781]	942	itau_phy,zjulian,dtime,nhoridbg,nidbg)
	943	! no vertical axis
	944	cl_surf(1)='ter'
	945	cl_surf(2)='lic'
	946	cl_surf(3)='oce'
	947	cl_surf(4)='sic'
	948	DO nsrf=1,nbsrf
[2344]	949	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",nbp_lon, &
	950	nbp_lat,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
[3391]	951	ENDDO
[781]	952
	953	CALL histend(nidbg)
	954	CALL histsync(nidbg)
	955
[3391]	956	ENDIF
[781]	957
	958	ENDIF
	959
	960	!****************************************************************************************
[889]	961	! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
	962	! instead of ORCHIDEE)
[1894]	963	IF (qsol0>=0.) THEN
[1067]	964	PRINT*,'WARNING : On impose qsol=',qsol0
[889]	965	qsol(:)=qsol0
[1067]	966	ENDIF
[889]	967	!****************************************************************************************
	968
	969	!****************************************************************************************
[781]	970	! 2) Initialization to zero
	971	!****************************************************************************************
[2181]	972	!
	973	! 2a) Initialization of all argument variables with INTENT(OUT)
	974	!****************************************************************************************
	975	cdragh(:)=0. ; cdragm(:)=0.
	976	zu1(:)=0. ; zv1(:)=0.
[2227]	977	!albedo SB >>>
	978	alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0.
	979	!albedo SB <<<
[2181]	980	zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0.
	981	d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0.
	982	zxfluxlat(:)=0.
	983	zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
[4013]	984	zn2mout(:,:)=0 ;
[2181]	985	d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
[3179]	986	zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
[2181]	987	zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
	988	cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0.
	989	kh(:)=0. ; kh_x(:)=0. ; kh_w(:)=0.
	990	slab_wfbils(:)=0.
	991	s_pblh(:)=0. ; s_pblh_x(:)=0. ; s_pblh_w(:)=0.
	992	s_plcl(:)=0. ; s_plcl_x(:)=0. ; s_plcl_w(:)=0.
	993	s_capCL(:)=0. ; s_oliqCL(:)=0. ; s_cteiCL(:)=0. ; s_pblT(:)=0.
	994	s_therm(:)=0.
	995	s_trmb1(:)=0. ; s_trmb2(:)=0. ; s_trmb3(:)=0.
[2243]	996	zustar(:)=0.
[2181]	997	zu10m(:)=0. ; zv10m(:)=0.
	998	fder_print(:)=0.
	999	zxqsurf(:)=0.
[4013]	1000	delta_qsurf(:) = 0.
[2181]	1001	zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
	1002	solsw(:,:)=0. ; sollw(:,:)=0.
	1003	d_ts(:,:)=0.
[2243]	1004	evap(:,:)=0.
[2181]	1005	fluxlat(:,:)=0.
	1006	wfbils(:,:)=0. ; wfbilo(:,:)=0.
[2670]	1007	wfevap(:,:)=0. ; wfrain(:,:)=0. ; wfsnow(:,:)=0.
[2181]	1008	flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.
	1009	dflux_t(:)=0. ; dflux_q(:)=0.
	1010	zxsnow(:)=0.
	1011	zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.
	1012	qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
	1013	runoff(:)=0.
	1014	IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
	1015	zcoefh(:,:,:) = 0.0
	1016	zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used
	1017	zcoefm(:,:,:) = 0.0
	1018	zcoefm(:,1,:) = 999999. !
	1019	ELSE
	1020	zcoefm(:,:,is_ave)=0.
	1021	zcoefh(:,:,is_ave)=0.
	1022	ENDIF
	1023	!!
	1024	! The components "is_ave" of tke_x and wake_deltke are "OUT" variables
	1025	!jyg<
	1026	!! tke(:,:,is_ave)=0.
	1027	tke_x(:,:,is_ave)=0.
[3798]	1028
[2181]	1029	wake_dltke(:,:,is_ave)=0.
	1030	!>jyg
	1031	!!! jyg le 23/02/2013
	1032	t2m(:,:) = 999999. ! t2m and q2m are meaningfull only over sub-surfaces
	1033	q2m(:,:) = 999999. ! actually present in the grid cell.
	1034	!!!
	1035	rh2m(:) = 0. ; qsat2m(:) = 0.
	1036	!!!
	1037	!!! jyg le 10/02/2012
	1038	rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0.
	1039
	1040	! 2b) Initialization of all local variables that will be compressed later
	1041	!****************************************************************************************
	1042	!! cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0
[1067]	1043	ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0
[2241]	1044	!! zv1 = 0.0 ; yqsurf = 0.0
[2227]	1045	!albedo SB >>>
	1046	yqsurf = 0.0 ; yalb = 0.0 ; yalb_vis = 0.0
	1047	!albedo SB <<<
[3798]	1048	yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0
[2243]	1049	ysollw = 0.0 ; yz0m = 0.0 ; yz0h = 0.0 ; yu1 = 0.0
[888]	1050	yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0
[781]	1051	ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0
[996]	1052	yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
[1067]	1053	yrugoro = 0.0 ; ywindsp = 0.0
[2181]	1054	!! d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0
	1055	yfluxlat=0.0
	1056	!! flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0
	1057	!! d_t_diss= 0.0 ;d_u = 0.0 ; d_v = 0.0
	1058	yqsol = 0.0
[3798]	1059
	1060	ytke=0.
[2952]	1061	!FC
	1062	y_treedrg=0.
	1063
[1865]	1064	! Martin
	1065	ysnowhgt = 0.0; yqsnow = 0.0 ; yrunoff = 0.0 ; ytoice =0.0
[2188]	1066	yalb3_new = 0.0 ; ysissnow = 0.0
[4013]	1067	ycldt = 0.0 ; yrmu0 = 0.0
[1865]	1068	! Martin
[2159]	1069
	1070	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
	1071	ytke_x=0. ; ytke_w=0. ; ywake_dltke=0.
	1072	y_d_t_x=0. ; y_d_t_w=0. ; y_d_q_x=0. ; y_d_q_w=0.
[2181]	1073	!! d_t_w=0. ; d_q_w=0.
	1074	!! d_t_x=0. ; d_q_x=0.
[2159]	1075	!! d_wake_dlt=0. ; d_wake_dlq=0.
	1076	yfluxlat_x=0. ; yfluxlat_w=0.
	1077	ywake_s=0. ; ywake_cstar=0. ;ywake_dens=0.
	1078	!!!
	1079	!!! nrlmd le 13/06/2011
	1080	tau_eq=0. ; delta_coef=0.
	1081	y_delta_flux_t1=0.
	1082	ydtsurf_th=0.
[4013]	1083	yts_x(:)=0. ; yts_w(:)=0.
	1084	y_delta_tsurf(:)=0. ; y_delta_qsurf(:)=0.
	1085	yqsurf_x(:)=0. ; yqsurf_w(:)=0.
	1086	yg_T(:) = 0. ; yg_Q(:) = 0.
	1087	yGamma_dTs_phiT(:) = 0. ; yGamma_dQs_phiQ(:) = 0.
	1088	ydTs_ins(:) = 0. ; ydqs_ins(:) = 0.
	1089
[2159]	1090	!!!
[781]	1091	ytsoil = 999999.
[2952]	1092	!FC
	1093	y_d_u_frein(:,:)=0.
	1094	y_d_v_frein(:,:)=0.
	1095	!FC
[781]	1096
[3391]	1097	! >> PC
	1098	!the yfields_out variable is defined in (klon,nbcf_out) even if it is used on
	1099	!the ORCHIDEE grid and as such should be defined in yfields_out(knon,nbcf_out) but
	1100	!the knon variable is not known at that level of pbl_surface_mod
[2181]	1101
[3391]	1102	!the yfields_in variable is defined in (klon,nbcf_in) even if it is used on the
	1103	!ORCHIDEE grid and as such should be defined in yfields_in(knon,nbcf_in) but the
	1104	!knon variable is not known at that level of pbl_surface_mod
[2952]	1105
[3391]	1106	yfields_out(:,:) = 0.
	1107	! << PC
[2952]	1108
[3391]	1109
[2181]	1110	! 2c) Initialization of all local variables computed within the subsurface loop and used later on
	1111	!****************************************************************************************
	1112	d_t_diss_x(:,:) = 0. ; d_t_diss_w(:,:) = 0.
	1113	d_u_x(:,:)=0. ; d_u_w(:,:)=0.
	1114	d_v_x(:,:)=0. ; d_v_w(:,:)=0.
	1115	flux_t_x(:,:,:)=0. ; flux_t_w(:,:,:)=0.
	1116	flux_q_x(:,:,:)=0. ; flux_q_w(:,:,:)=0.
[2159]	1117	!
[2181]	1118	!jyg<
	1119	flux_u_x(:,:,:)=0. ; flux_u_w(:,:,:)=0.
	1120	flux_v_x(:,:,:)=0. ; flux_v_w(:,:,:)=0.
	1121	fluxlat_x(:,:)=0. ; fluxlat_w(:,:)=0.
	1122	!>jyg
	1123	!
	1124	!jyg<
	1125	! pblh,plcl,capCL,cteiCL ... are meaningfull only over sub-surfaces
	1126	! actually present in the grid cell ==> value set to 999999.
	1127	!
	1128	!jyg<
	1129	ustar(:,:) = 999999.
	1130	wstar(:,:) = 999999.
	1131	windsp(:,:) = SQRT(u10m(:,:)2 + v10m(:,:)2 )
	1132	u10m(:,:) = 999999.
	1133	v10m(:,:) = 999999.
	1134	!>jyg
	1135	!
	1136	pblh(:,:) = 999999. ! Hauteur de couche limite
	1137	plcl(:,:) = 999999. ! Niveau de condensation de la CLA
	1138	capCL(:,:) = 999999. ! CAPE de couche limite
	1139	oliqCL(:,:) = 999999. ! eau_liqu integree de couche limite
	1140	cteiCL(:,:) = 999999. ! cloud top instab. crit. couche limite
	1141	pblt(:,:) = 999999. ! T a la Hauteur de couche limite
	1142	therm(:,:) = 999999.
	1143	trmb1(:,:) = 999999. ! deep_cape
	1144	trmb2(:,:) = 999999. ! inhibition
	1145	trmb3(:,:) = 999999. ! Point Omega
	1146	!
	1147	t2m_x(:,:) = 999999.
	1148	q2m_x(:,:) = 999999.
	1149	ustar_x(:,:) = 999999.
	1150	wstar_x(:,:) = 999999.
	1151	u10m_x(:,:) = 999999.
	1152	v10m_x(:,:) = 999999.
	1153	!
	1154	pblh_x(:,:) = 999999. ! Hauteur de couche limite
	1155	plcl_x(:,:) = 999999. ! Niveau de condensation de la CLA
	1156	capCL_x(:,:) = 999999. ! CAPE de couche limite
	1157	oliqCL_x(:,:) = 999999. ! eau_liqu integree de couche limite
	1158	cteiCL_x(:,:) = 999999. ! cloud top instab. crit. couche limite
	1159	pblt_x(:,:) = 999999. ! T a la Hauteur de couche limite
	1160	therm_x(:,:) = 999999.
	1161	trmb1_x(:,:) = 999999. ! deep_cape
	1162	trmb2_x(:,:) = 999999. ! inhibition
	1163	trmb3_x(:,:) = 999999. ! Point Omega
	1164	!
	1165	t2m_w(:,:) = 999999.
	1166	q2m_w(:,:) = 999999.
	1167	ustar_w(:,:) = 999999.
	1168	wstar_w(:,:) = 999999.
	1169	u10m_w(:,:) = 999999.
	1170	v10m_w(:,:) = 999999.
	1171
	1172	pblh_w(:,:) = 999999. ! Hauteur de couche limite
	1173	plcl_w(:,:) = 999999. ! Niveau de condensation de la CLA
	1174	capCL_w(:,:) = 999999. ! CAPE de couche limite
	1175	oliqCL_w(:,:) = 999999. ! eau_liqu integree de couche limite
	1176	cteiCL_w(:,:) = 999999. ! cloud top instab. crit. couche limite
	1177	pblt_w(:,:) = 999999. ! T a la Hauteur de couche limite
	1178	therm_w(:,:) = 999999.
	1179	trmb1_w(:,:) = 999999. ! deep_cape
	1180	trmb2_w(:,:) = 999999. ! inhibition
	1181	trmb3_w(:,:) = 999999. ! Point Omega
	1182	!!!
	1183	!
[2159]	1184	!!!
[781]	1185	!****************************************************************************************
	1186	! 3) - Calculate pressure thickness of each layer
	1187	! - Calculate the wind at first layer
[888]	1188	! - Mean calculations of albedo
	1189	! - Calculate net radiance at sub-surface
[781]	1190	!****************************************************************************************
	1191	DO k = 1, klev
	1192	DO i = 1, klon
	1193	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	1194	ENDDO
	1195	ENDDO
	1196
	1197	!****************************************************************************************
	1198	! Test for rugos........ from physiq.. A la fin plutot???
[888]	1199	!
[781]	1200	!****************************************************************************************
	1201
	1202	DO nsrf = 1, nbsrf
	1203	DO i = 1, klon
[2243]	1204	z0m(i,nsrf) = MAX(z0m(i,nsrf),z0min)
	1205	z0h(i,nsrf) = MAX(z0h(i,nsrf),z0min)
[781]	1206	ENDDO
	1207	ENDDO
	1208
[888]	1209	! Mean calculations of albedo
	1210	!
	1211	! * alb : mean albedo for whole SW interval
	1212	!
	1213	! Mean albedo for grid point
	1214	! * alb_m : mean albedo at whole SW interval
	1215
[2227]	1216	alb_dir_m(:,:) = 0.0
	1217	alb_dif_m(:,:) = 0.0
	1218	DO k = 1, nsw
	1219	DO nsrf = 1, nbsrf
[781]	1220	DO i = 1, klon
[2227]	1221	alb_dir_m(i,k) = alb_dir_m(i,k) + alb_dir(i,k,nsrf) * pctsrf(i,nsrf)
	1222	alb_dif_m(i,k) = alb_dif_m(i,k) + alb_dif(i,k,nsrf) * pctsrf(i,nsrf)
[781]	1223	ENDDO
[2227]	1224	ENDDO
[781]	1225	ENDDO
	1226
[888]	1227	! We here suppose the fraction f1 of incoming radiance of visible radiance
	1228	! as a fraction of all shortwave radiance
[2227]	1229	f1 = 0.5
[1069]	1230	! f1 = 1 ! put f1=1 to recreate old calculations
[781]	1231
[2227]	1232	!f1 is already included with SFRWL values in each surf files
	1233	alb=0.0
	1234	DO k=1,nsw
	1235	DO nsrf = 1, nbsrf
	1236	DO i = 1, klon
	1237	alb(i,nsrf) = alb(i,nsrf) + alb_dir(i,k,nsrf)*SFRWL(k)
	1238	ENDDO
	1239	ENDDO
[888]	1240	ENDDO
[781]	1241
[2227]	1242	alb_m=0.0
	1243	DO k = 1,nsw
	1244	DO i = 1, klon
	1245	alb_m(i) = alb_m(i) + alb_dir_m(i,k)*SFRWL(k)
[3391]	1246	ENDDO
[2227]	1247	ENDDO
	1248	!albedo SB <<<
[888]	1249
[2227]	1250
	1251
[888]	1252	! Calculation of mean temperature at surface grid points
[781]	1253	ztsol(:) = 0.0
	1254	DO nsrf = 1, nbsrf
	1255	DO i = 1, klon
	1256	ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
	1257	ENDDO
	1258	ENDDO
	1259
[888]	1260	! Linear distrubution on sub-surface of long- and shortwave net radiance
[781]	1261	DO nsrf = 1, nbsrf
	1262	DO i = 1, klon
	1263	sollw(i,nsrf) = sollw_m(i) + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ts(i,nsrf))
[3798]	1264	!--OB this line is not satisfactory because alb is the direct albedo not total albedo
[888]	1265	solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
[781]	1266	ENDDO
	1267	ENDDO
[2455]	1268	!
	1269	!<al1: second order corrections
	1270	!- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)
	1271	IF (iflag_order2_sollw == 1) THEN
	1272	meansqT(:) = 0. ! as working buffer
	1273	DO nsrf = 1, nbsrf
	1274	DO i = 1, klon
	1275	meansqT(i) = meansqT(i)+(ts(i,nsrf)-ztsol(i))*2 pctsrf(i,nsrf)
[3391]	1276	ENDDO
	1277	ENDDO
[2455]	1278	DO nsrf = 1, nbsrf
	1279	DO i = 1, klon
	1280	sollw(i,nsrf) = sollw(i,nsrf) &
	1281	+ 6.0RSIGMAztsol(i)*2 (meansqT(i)-(ztsol(i)-ts(i,nsrf))**2)
	1282	ENDDO
	1283	ENDDO
	1284	ENDIF ! iflag_order2_sollw == 1
	1285	!>al1
[781]	1286
[3798]	1287	!--OB add diffuse fraction of SW down
	1288	DO n=1,nbcf_out
	1289	IF (cfname_out(n) == "swdownfdiff" ) fields_out(:,n) = solswfdiff_m(:)
	1290	ENDDO
[3391]	1291	! >> PC
	1292	IF (carbon_cycle_cpl .AND. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN
	1293	r_co2_ppm(:) = co2_send(:)
	1294	DO n=1,nbcf_out
	1295	IF (cfname_out(n) == "atmco2" ) fields_out(:,n) = co2_send(:)
	1296	ENDDO
	1297	ENDIF
	1298	IF ( .NOT. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN
	1299	r_co2_ppm(:) = co2_ppm ! Constant field
	1300	DO n=1,nbcf_out
	1301	IF (cfname_out(n) == "atmco2" ) fields_out(:,n) = co2_ppm
	1302	ENDDO
	1303	ENDIF
	1304	! << PC
	1305
[781]	1306	!****************************************************************************************
	1307	! 4) Loop over different surfaces
	1308	!
[2159]	1309	! Only points containing a fraction of the sub surface will be treated.
[781]	1310	!
	1311	!****************************************************************************************
[4013]	1312	!<<<<<<<<<<<<<
	1313	loop_nbsrf: DO nsrf = 1, nbsrf !<<<<<<<<<<<<<
	1314	!<<<<<<<<<<<<<
[2159]	1315	IF (prt_level >=10) print *,' Loop nsrf ',nsrf
[4013]	1316	!
	1317	IF (iflag_split_ref == 3) THEN
	1318	IF (nsrf == is_oce) THEN
	1319	iflag_split = 1
	1320	ELSE
	1321	iflag_split=0
	1322	ENDIF !! (nsrf == is_oce)
	1323	ELSE
	1324	iflag_split = iflag_split_ref
	1325	ENDIF !! (iflag_split_ref == 3)
[781]	1326
	1327	! Search for index(ni) and size(knon) of domaine to treat
	1328	ni(:) = 0
	1329	knon = 0
	1330	DO i = 1, klon
[996]	1331	IF (pctsrf(i,nsrf) > 0.) THEN
[781]	1332	knon = knon + 1
	1333	ni(knon) = i
	1334	ENDIF
	1335	ENDDO
	1336
[2159]	1337	!!! jyg le 19/08/2012
[2186]	1338	! IF (knon <= 0) THEN
	1339	! IF (prt_level >= 10) print *,' no grid point for nsrf= ',nsrf
	1340	! cycle loop_nbsrf
	1341	! ENDIF
[2159]	1342	!!!
	1343
[781]	1344	! write index, with IOIPSL
	1345	IF (debugindex .AND. mpi_size==1) THEN
	1346	tabindx(:)=0.
	1347	DO i=1,knon
[1403]	1348	tabindx(i)=REAL(i)
[3391]	1349	ENDDO
[781]	1350	debugtab(:,:) = 0.
	1351	ndexbg(:) = 0
	1352	CALL gath2cpl(tabindx,debugtab,knon,ni)
[2344]	1353	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,nbp_lon*nbp_lat, ndexbg)
[781]	1354	ENDIF
	1355
	1356	!****************************************************************************************
	1357	! 5) Compress variables
	1358	!
	1359	!****************************************************************************************
	1360
[4013]	1361	!
	1362	!jyg< (20190926)
	1363	! Provisional : set ybeta to standard values
	1364	IF (nsrf .NE. is_ter) THEN
	1365	ybeta(:) = 1.
	1366	ELSE
	1367	IF (iflag_split .EQ. 0) THEN
	1368	ybeta(:) = 1.
	1369	ELSE
	1370	DO j = 1, knon
	1371	i = ni(j)
	1372	ybeta(j) = beta(i,nsrf)
	1373	ENDDO
	1374	ENDIF ! (iflag_split .LE.1)
	1375	ENDIF ! (nsrf .NE. is_ter)
	1376	!>jyg
	1377	!
[781]	1378	DO j = 1, knon
	1379	i = ni(j)
[888]	1380	ypct(j) = pctsrf(i,nsrf)
	1381	yts(j) = ts(i,nsrf)
	1382	ysnow(j) = snow(i,nsrf)
	1383	yqsurf(j) = qsurf(i,nsrf)
	1384	yalb(j) = alb(i,nsrf)
[2227]	1385	!albedo SB >>>
	1386	yalb_vis(j) = alb_dir(i,1,nsrf)
[3391]	1387	IF (nsw==6) THEN
[2227]	1388	yalb_vis(j)=(alb_dir(i,1,nsrf)SFRWL(1)+alb_dir(i,2,nsrf)SFRWL(2) &
	1389	+alb_dir(i,3,nsrf)*SFRWL(3))/(SFRWL(1)+SFRWL(2)+SFRWL(3))
[3391]	1390	ENDIF
[2227]	1391	!albedo SB <<<
[781]	1392	yrain_f(j) = rain_f(i)
	1393	ysnow_f(j) = snow_f(i)
	1394	yagesno(j) = agesno(i,nsrf)
[888]	1395	yfder(j) = fder(i)
[2188]	1396	ylwdown(j) = lwdown_m(i)
[2240]	1397	ygustiness(j) = gustiness(i)
[888]	1398	ysolsw(j) = solsw(i,nsrf)
	1399	ysollw(j) = sollw(i,nsrf)
[2243]	1400	yz0m(j) = z0m(i,nsrf)
	1401	yz0h(j) = z0h(i,nsrf)
[781]	1402	yrugoro(j) = rugoro(i)
[1067]	1403	yu1(j) = u(i,1)
	1404	yv1(j) = v(i,1)
[781]	1405	ypaprs(j,klev+1) = paprs(i,klev+1)
[2181]	1406	!jyg<
	1407	!! ywindsp(j) = SQRT(u10m(i,nsrf)2 + v10m(i,nsrf)2 )
	1408	ywindsp(j) = windsp(i,nsrf)
	1409	!>jyg
[4013]	1410	! Martin and Etienne
	1411	yzmea(j) = zmea(i)
[1865]	1412	yzsig(j) = zsig(i)
	1413	ycldt(j) = cldt(i)
	1414	yrmu0(j) = rmu0(i)
	1415	! Martin
[2159]	1416	!!! nrlmd le 13/06/2011
	1417	y_delta_tsurf(j)=delta_tsurf(i,nsrf)
	1418	!!!
[3391]	1419	ENDDO
	1420	! >> PC
	1421	!--compressing fields_out onto ORCHIDEE grid
	1422	!--these fields are shared and used directly surf_land_orchidee_mod
	1423	DO n = 1, nbcf_out
	1424	DO j = 1, knon
	1425	i = ni(j)
	1426	yfields_out(j,n) = fields_out(i,n)
	1427	ENDDO
	1428	ENDDO
	1429	! << PC
[781]	1430	DO k = 1, klev
	1431	DO j = 1, knon
	1432	i = ni(j)
	1433	ypaprs(j,k) = paprs(i,k)
	1434	ypplay(j,k) = pplay(i,k)
[996]	1435	ydelp(j,k) = delp(i,k)
[2159]	1436	ENDDO
	1437	ENDDO
[3391]	1438	!
[2159]	1439	!!! jyg le 07/02/2012 et le 10/04/2013
[3198]	1440	DO k = 1, klev+1
[2159]	1441	DO j = 1, knon
	1442	i = ni(j)
[2181]	1443	!jyg<
	1444	!! ytke(j,k) = tke(i,k,nsrf)
	1445	ytke(j,k) = tke_x(i,k,nsrf)
[3198]	1446	ENDDO
	1447	ENDDO
[2181]	1448	!>jyg
[3198]	1449	DO k = 1, klev
	1450	DO j = 1, knon
	1451	i = ni(j)
[2952]	1452	y_treedrg(j,k) = treedrg(i,k,nsrf)
[781]	1453	yu(j,k) = u(i,k)
	1454	yv(j,k) = v(i,k)
	1455	yt(j,k) = t(i,k)
	1456	yq(j,k) = q(i,k)
	1457	ENDDO
[2159]	1458	ENDDO
	1459	!
[3391]	1460	IF (iflag_split.GE.1) THEN
[2159]	1461	!!! nrlmd le 02/05/2011
	1462	DO k = 1, klev
	1463	DO j = 1, knon
	1464	i = ni(j)
	1465	yu_x(j,k) = u(i,k)
	1466	yv_x(j,k) = v(i,k)
	1467	yt_x(j,k) = t(i,k)-wake_s(i)*wake_dlt(i,k)
	1468	yq_x(j,k) = q(i,k)-wake_s(i)*wake_dlq(i,k)
	1469	yu_w(j,k) = u(i,k)
	1470	yv_w(j,k) = v(i,k)
	1471	yt_w(j,k) = t(i,k)+(1.-wake_s(i))*wake_dlt(i,k)
	1472	yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k)
	1473	!!!
	1474	ENDDO
	1475	ENDDO
[3391]	1476
[3179]	1477	IF (prt_level .ge. 10) THEN
	1478	print *,'pbl_surface, wake_s(1), wake_dlt(1,:) ', wake_s(1), wake_dlt(1,:)
	1479	print *,'pbl_surface, wake_s(1), wake_dlq(1,:) ', wake_s(1), wake_dlq(1,:)
	1480	ENDIF
[3391]	1481
[2159]	1482	!!! nrlmd le 02/05/2011
	1483	DO k = 1, klev+1
	1484	DO j = 1, knon
	1485	i = ni(j)
[2181]	1486	!jyg<
	1487	!! ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)*wake_dltke(i,k,nsrf)
	1488	!! ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))*wake_dltke(i,k,nsrf)
	1489	!! ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
	1490	!! ytke(j,k) = tke(i,k,nsrf)
	1491	!
	1492	ytke_x(j,k) = tke_x(i,k,nsrf)
	1493	ytke(j,k) = tke_x(i,k,nsrf)+wake_s(i)*wake_dltke(i,k,nsrf)
	1494	ytke_w(j,k) = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf)
[2159]	1495	ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
[3798]	1496
[2181]	1497	!>jyg
[2159]	1498	ENDDO
	1499	ENDDO
	1500	!!!
	1501	!!! jyg le 07/02/2012
	1502	DO j = 1, knon
	1503	i = ni(j)
	1504	ywake_s(j)=wake_s(i)
	1505	ywake_cstar(j)=wake_cstar(i)
	1506	ywake_dens(j)=wake_dens(i)
	1507	ENDDO
	1508	!!!
	1509	!!! nrlmd le 13/06/2011
	1510	DO j=1,knon
	1511	yts_x(j)=yts(j)-ywake_s(j)*y_delta_tsurf(j)
	1512	yts_w(j)=yts(j)+(1.-ywake_s(j))*y_delta_tsurf(j)
	1513	ENDDO
	1514	!!!
[2852]	1515	ENDIF ! (iflag_split .ge.1)
[2159]	1516	!!!
[781]	1517	DO k = 1, nsoilmx
	1518	DO j = 1, knon
	1519	i = ni(j)
	1520	ytsoil(j,k) = ftsoil(i,k,nsrf)
[3391]	1521	ENDDO
	1522	ENDDO
[781]	1523
	1524	! qsol(water height in soil) only for bucket continental model
	1525	IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
	1526	DO j = 1, knon
	1527	i = ni(j)
	1528	yqsol(j) = qsol(i)
[3391]	1529	ENDDO
[781]	1530	ENDIF
[3429]	1531
[3458]	1532	if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
[3740]	1533	if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
[3767]	1534	ydelta_sal(:knon) = delta_sal(ni(:knon))
[3744]	1535	ydelta_sst(:knon) = delta_sst(ni(:knon))
[4020]	1536	ydter(:knon) = dter(ni(:knon))
	1537	ydser(:knon) = dser(ni(:knon))
	1538	ydt_ds(:knon) = dt_ds(ni(:knon))
[3740]	1539	end if
	1540
[3429]	1541	yds_ns(:knon) = ds_ns(ni(:knon))
	1542	ydt_ns(:knon) = dt_ns(ni(:knon))
	1543	end if
[781]	1544
	1545	!****************************************************************************************
[1067]	1546	! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
[781]	1547	!
	1548	!****************************************************************************************
	1549
[4013]	1550
[2159]	1551	!!! jyg le 07/02/2012
	1552	IF (iflag_split .eq.0) THEN
	1553	!!!
	1554	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
[2232]	1555	! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
	1556	! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
	1557	! yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
	1558	! yts, yqsurf, yrugos, &
	1559	! ycdragm, ycdragh )
	1560	! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
[2234]	1561	DO i = 1, knon
	1562	! print*,'PBL ',i,RD
	1563	! print*,'PBL ',yt(i,1),ypaprs(i,1),ypplay(i,1)
[2232]	1564	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
	1565	* (ypaprs(i,1)-ypplay(i,1))
	1566	speed(i) = SQRT(yu(i,1)2+yv(i,1)2)
[3391]	1567	ENDDO
[2232]	1568	CALL cdrag(knon, nsrf, &
	1569	speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
[2243]	1570	yts, yqsurf, yz0m, yz0h, &
[2232]	1571	ycdragm, ycdragh, zri1, pref )
	1572
[2126]	1573	! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
[3392]	1574	IF (ok_prescr_ust) THEN
[2126]	1575	DO i = 1, knon
	1576	print *,'ycdragm avant=',ycdragm(i)
	1577	vent= sqrt(yu(i,1)yu(i,1)+yv(i,1)yv(i,1))
	1578	! ycdragm(i) = ustust/(1.+(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1)))
	1579	! ycdragm(i) = ustust/((1.+sqrt(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1))) &
	1580	! sqrt(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1)))
	1581	ycdragm(i) = ust*ust/(1.+vent)/vent
[2307]	1582	! print *,'ycdragm ust yu yv apres=',ycdragm(i),ust,yu(i,1),yv(i,1)
[2126]	1583	ENDDO
	1584	ENDIF
[3798]	1585
[2159]	1586	IF (prt_level >=10) print *,'clcdrag -> ycdragh ', ycdragh
	1587	ELSE !(iflag_split .eq.0)
[2232]	1588
	1589	! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
	1590	! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
	1591	! yu_x(:,1), yv_x(:,1), yt_x(:,1), yq_x(:,1), &
	1592	! yts_x, yqsurf, yrugos, &
	1593	! ycdragm_x, ycdragh_x )
	1594	! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
	1595	DO i = 1, knon
	1596	zgeo1_x(i) = RD * yt_x(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
	1597	* (ypaprs(i,1)-ypplay(i,1))
	1598	speed_x(i) = SQRT(yu_x(i,1)2+yv_x(i,1)2)
[3391]	1599	ENDDO
[4013]	1600
	1601
	1602	CALL cdrag(knon, nsrf, &
[2232]	1603	speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
[4013]	1604	yts_x, yqsurf_x, yz0m, yz0h, &
[2232]	1605	ycdragm_x, ycdragh_x, zri1_x, pref_x )
	1606
[2159]	1607	! --- special Dice. JYG+MPL 25112013
[3392]	1608	IF (ok_prescr_ust) THEN
[2159]	1609	DO i = 1, knon
[2307]	1610	! print *,'ycdragm_x avant=',ycdragm_x(i)
[2159]	1611	vent= sqrt(yu_x(i,1)yu_x(i,1)+yv_x(i,1)yv_x(i,1))
	1612	ycdragm_x(i) = ust*ust/(1.+vent)/vent
[2307]	1613	! print *,'ycdragm_x ust yu yv apres=',ycdragm_x(i),ust,yu_x(i,1),yv_x(i,1)
[2159]	1614	ENDDO
	1615	ENDIF
	1616	IF (prt_level >=10) print *,'clcdrag -> ycdragh_x ', ycdragh_x
	1617	!
[2255]	1618	! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
	1619	! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
	1620	! yu_w(:,1), yv_w(:,1), yt_w(:,1), yq_w(:,1), &
	1621	! yts_w, yqsurf, yz0m, &
	1622	! ycdragm_w, ycdragh_w )
	1623	! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
	1624	DO i = 1, knon
	1625	zgeo1_w(i) = RD * yt_w(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
	1626	* (ypaprs(i,1)-ypplay(i,1))
	1627	speed_w(i) = SQRT(yu_w(i,1)2+yv_w(i,1)2)
[3391]	1628	ENDDO
[2255]	1629	CALL cdrag(knon, nsrf, &
	1630	speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
[4013]	1631	yts_w, yqsurf_w, yz0m, yz0h, &
[2255]	1632	ycdragm_w, ycdragh_w, zri1_w, pref_w )
[3179]	1633	!
	1634	zgeo1(:) = wake_s(:)zgeo1_w(:) + (1.-wake_s(:))zgeo1_x(:)
[2243]	1635
[2307]	1636	! --- special Dice. JYG+MPL 25112013 puis BOMEX
[3392]	1637	IF (ok_prescr_ust) THEN
[2159]	1638	DO i = 1, knon
[2307]	1639	! print *,'ycdragm_w avant=',ycdragm_w(i)
[2159]	1640	vent= sqrt(yu_w(i,1)yu_w(i,1)+yv_w(i,1)yv_w(i,1))
	1641	ycdragm_w(i) = ust*ust/(1.+vent)/vent
[2307]	1642	! print *,'ycdragm_w ust yu yv apres=',ycdragm_w(i),ust,yu_w(i,1),yv_w(i,1)
[2159]	1643	ENDDO
	1644	ENDIF
	1645	IF (prt_level >=10) print *,'clcdrag -> ycdragh_w ', ycdragh_w
	1646	!!!
	1647	ENDIF ! (iflag_split .eq.0)
	1648	!!!
	1649
[1067]	1650
	1651	!****************************************************************************************
[2159]	1652	! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm.
[1067]	1653	!
	1654	!****************************************************************************************
	1655
[2159]	1656	!!! jyg le 07/02/2012
	1657	IF (iflag_split .eq.0) THEN
	1658	!!!
	1659	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
	1660	IF (prt_level >=10) THEN
	1661	print *,' args coef_diff_turb: yu ', yu
	1662	print *,' args coef_diff_turb: yv ', yv
	1663	print *,' args coef_diff_turb: yq ', yq
	1664	print *,' args coef_diff_turb: yt ', yt
	1665	print *,' args coef_diff_turb: yts ', yts
[2243]	1666	print *,' args coef_diff_turb: yz0m ', yz0m
[2159]	1667	print *,' args coef_diff_turb: yqsurf ', yqsurf
	1668	print *,' args coef_diff_turb: ycdragm ', ycdragm
	1669	print *,' args coef_diff_turb: ycdragh ', ycdragh
	1670	print *,' args coef_diff_turb: ytke ', ytke
[3798]	1671
[2159]	1672	ENDIF
	1673	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
[2243]	1674	ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
[2952]	1675	ycoefm, ycoefh, ytke, y_treedrg)
	1676	! ycoefm, ycoefh, ytke)
[3179]	1677	!FC y_treedrg ajoute
[1761]	1678	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
	1679	! In this case, coef_diff_turb is called for the Cd only
	1680	DO k = 2, klev
	1681	DO j = 1, knon
	1682	i = ni(j)
	1683	ycoefh(j,k) = zcoefh(i,k,nsrf)
	1684	ycoefm(j,k) = zcoefm(i,k,nsrf)
	1685	ENDDO
	1686	ENDDO
	1687	ENDIF
[2159]	1688	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh
	1689	!
	1690	ELSE !(iflag_split .eq.0)
	1691	IF (prt_level >=10) THEN
	1692	print *,' args coef_diff_turb: yu_x ', yu_x
	1693	print *,' args coef_diff_turb: yv_x ', yv_x
	1694	print *,' args coef_diff_turb: yq_x ', yq_x
	1695	print *,' args coef_diff_turb: yt_x ', yt_x
	1696	print *,' args coef_diff_turb: yts_x ', yts_x
	1697	print *,' args coef_diff_turb: yqsurf ', yqsurf
	1698	print *,' args coef_diff_turb: ycdragm_x ', ycdragm_x
	1699	print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x
	1700	print *,' args coef_diff_turb: ytke_x ', ytke_x
[3798]	1701
[2159]	1702	ENDIF
	1703	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
[4013]	1704	ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
[2952]	1705	ycoefm_x, ycoefh_x, ytke_x,y_treedrg)
	1706	! ycoefm_x, ycoefh_x, ytke_x)
	1707	!FC doit on le mettre ( on ne l utilise pas si il y a du spliting)
[2159]	1708	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
	1709	! In this case, coef_diff_turb is called for the Cd only
	1710	DO k = 2, klev
	1711	DO j = 1, knon
	1712	i = ni(j)
	1713	ycoefh_x(j,k) = zcoefh(i,k,nsrf)
	1714	ycoefm_x(j,k) = zcoefm(i,k,nsrf)
	1715	ENDDO
	1716	ENDDO
	1717	ENDIF
	1718	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x
	1719	!
	1720	IF (prt_level >=10) THEN
	1721	print *,' args coef_diff_turb: yu_w ', yu_w
	1722	print *,' args coef_diff_turb: yv_w ', yv_w
	1723	print *,' args coef_diff_turb: yq_w ', yq_w
	1724	print *,' args coef_diff_turb: yt_w ', yt_w
	1725	print *,' args coef_diff_turb: yts_w ', yts_w
	1726	print *,' args coef_diff_turb: yqsurf ', yqsurf
	1727	print *,' args coef_diff_turb: ycdragm_w ', ycdragm_w
	1728	print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w
	1729	print *,' args coef_diff_turb: ytke_w ', ytke_w
	1730	ENDIF
	1731	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
[4013]	1732	ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
[2952]	1733	ycoefm_w, ycoefh_w, ytke_w,y_treedrg)
	1734	! ycoefm_w, ycoefh_w, ytke_w)
[2159]	1735	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
	1736	! In this case, coef_diff_turb is called for the Cd only
	1737	DO k = 2, klev
	1738	DO j = 1, knon
	1739	i = ni(j)
	1740	ycoefh_w(j,k) = zcoefh(i,k,nsrf)
	1741	ycoefm_w(j,k) = zcoefm(i,k,nsrf)
	1742	ENDDO
	1743	ENDDO
	1744	ENDIF
	1745	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w
	1746	!
	1747	!!!jyg le 10/04/2013
	1748	!! En attendant de traiter le transport des traceurs dans les poches froides, formule
	1749	!! arbitraire pour ycoefh et ycoefm
	1750	DO k = 2,klev
	1751	DO j = 1,knon
	1752	ycoefh(j,k) = ycoefh_x(j,k) + ywake_s(j)*(ycoefh_w(j,k) - ycoefh_x(j,k))
	1753	ycoefm(j,k) = ycoefm_x(j,k) + ywake_s(j)*(ycoefm_w(j,k) - ycoefm_x(j,k))
	1754	ENDDO
	1755	ENDDO
	1756	!!!
	1757	ENDIF ! (iflag_split .eq.0)
	1758	!!!
[781]	1759
	1760	!****************************************************************************************
	1761	!
	1762	! 8) "La descente" - "The downhill"
	1763	!
	1764	! climb_hq_down and climb_wind_down calculate the coefficients
	1765	! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
	1766	! Only the coefficients at surface for H and Q are returned.
	1767	!
	1768	!****************************************************************************************
	1769
	1770	! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
[2159]	1771	!!! jyg le 07/02/2012
	1772	IF (iflag_split .eq.0) THEN
	1773	!!!
	1774	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
	1775	CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
[781]	1776	ydelp, yt, yq, dtime, &
[2159]	1777	!!! jyg le 09/05/2011
	1778	CcoefH, CcoefQ, DcoefH, DcoefQ, &
	1779	Kcoef_hq, gama_q, gama_h, &
	1780	!!!
[1067]	1781	AcoefH, AcoefQ, BcoefH, BcoefQ)
[2159]	1782	ELSE !(iflag_split .eq.0)
	1783	CALL climb_hq_down(knon, ycoefh_x, ypaprs, ypplay, &
	1784	ydelp, yt_x, yq_x, dtime, &
	1785	!!! nrlmd le 02/05/2011
	1786	CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
	1787	Kcoef_hq_x, gama_q_x, gama_h_x, &
	1788	!!!
	1789	AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x)
[3179]	1790	!!!
	1791	IF (prt_level >=10) THEN
	1792	PRINT *,'pbl_surface (climb_hq_down.x->) AcoefH_x ',AcoefH_x
	1793	PRINT *,'pbl_surface (climb_hq_down.x->) AcoefQ_x ',AcoefQ_x
	1794	PRINT *,'pbl_surface (climb_hq_down.x->) BcoefH_x ',BcoefH_x
	1795	PRINT *,'pbl_surface (climb_hq_down.x->) BcoefQ_x ',BcoefQ_x
	1796	ENDIF
[2159]	1797	!
	1798	CALL climb_hq_down(knon, ycoefh_w, ypaprs, ypplay, &
	1799	ydelp, yt_w, yq_w, dtime, &
	1800	!!! nrlmd le 02/05/2011
	1801	CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
	1802	Kcoef_hq_w, gama_q_w, gama_h_w, &
	1803	!!!
	1804	AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w)
	1805	!!!
[3179]	1806	IF (prt_level >=10) THEN
	1807	PRINT *,'pbl_surface (climb_hq_down.w->) AcoefH_w ',AcoefH_w
	1808	PRINT *,'pbl_surface (climb_hq_down.w->) AcoefQ_w ',AcoefQ_w
	1809	PRINT *,'pbl_surface (climb_hq_down.w->) BcoefH_w ',BcoefH_w
	1810	PRINT *,'pbl_surface (climb_hq_down.w->) BcoefQ_w ',BcoefQ_w
	1811	ENDIF
	1812	!!!
[2159]	1813	ENDIF ! (iflag_split .eq.0)
	1814	!!!
[781]	1815
	1816	! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
[2159]	1817	!!! jyg le 07/02/2012
	1818	IF (iflag_split .eq.0) THEN
	1819	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
	1820	CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
	1821	!!! jyg le 09/05/2011
	1822	CcoefU, CcoefV, DcoefU, DcoefV, &
	1823	Kcoef_m, alf_1, alf_2, &
	1824	!!!
[1067]	1825	AcoefU, AcoefV, BcoefU, BcoefV)
[2159]	1826	ELSE ! (iflag_split .eq.0)
	1827	CALL climb_wind_down(knon, dtime, ycoefm_x, ypplay, ypaprs, yt_x, ydelp, yu_x, yv_x, &
	1828	!!! nrlmd le 02/05/2011
	1829	CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
	1830	Kcoef_m_x, alf_1_x, alf_2_x, &
	1831	!!!
	1832	AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x)
	1833	!
	1834	CALL climb_wind_down(knon, dtime, ycoefm_w, ypplay, ypaprs, yt_w, ydelp, yu_w, yv_w, &
	1835	!!! nrlmd le 02/05/2011
	1836	CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
	1837	Kcoef_m_w, alf_1_w, alf_2_w, &
	1838	!!!
	1839	AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w)
	1840	!!!
	1841	ENDIF ! (iflag_split .eq.0)
	1842	!!!
[781]	1843
	1844	!****************************************************************************************
	1845	! 9) Small calculations
	1846	!
	1847	!****************************************************************************************
[888]	1848
	1849	! - Reference pressure is given the values at surface level
[781]	1850	ypsref(:) = ypaprs(:,1)
	1851
[1279]	1852	! - CO2 field on 2D grid to be sent to ORCHIDEE
	1853	! Transform to compressed field
	1854	IF (carbon_cycle_cpl) THEN
	1855	DO i=1,knon
	1856	r_co2_ppm(i) = co2_send(ni(i))
[3391]	1857	ENDDO
[1279]	1858	ELSE
	1859	r_co2_ppm(:) = co2_ppm ! Constant field
[3391]	1860	ENDIF
[781]	1861
[2159]	1862	!!! nrlmd le 02/05/2011 -----------------------On raccorde les 2 colonnes dans la couche 1
	1863	!----------------------------------------------------------------------------------------
	1864	!!! jyg le 07/02/2012
[3179]	1865	!!! jyg le 01/02/2017
	1866	IF (iflag_split .eq. 0) THEN
	1867	yt1(:) = yt(:,1)
	1868	yq1(:) = yq(:,1)
	1869	ELSE IF (iflag_split .ge. 1) THEN
[4013]	1870	!
	1871	! Cdragq computation
	1872	! ------------------
	1873	!******************************************************************************
	1874	! Cdragq computed from cdrag
	1875	! The difference comes only from a factor (f_z0qh_oce) on z0, so that
	1876	! it can be computed inside wx_pbl0_merge
	1877	! More complicated appraches may require the propagation through
	1878	! pbl_surface of an independant cdragq variable.
	1879	!******************************************************************************
	1880	!
	1881	IF ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) THEN
	1882	! Si on suit les formulations par exemple de Tessel, on
	1883	! a z0h=0.4nu/u, z0q=0.62nu/u, d'ou f_z0qh_oce=0.62/0.4=1.55
	1884	!! ycdragq_x(1:knon)=ycdragh_x(1:knon)* &
	1885	!! log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
	1886	!! ycdragq_w(1:knon)=ycdragh_w(1:knon)* &
	1887	!! log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
	1888	!
	1889	DO j = 1,knon
	1890	z1lay = zgeo1(j)/RG
	1891	fact_cdrag = log(z1lay/yz0h(j))/log(z1lay/(f_z0qh_oce*yz0h(j)))
	1892	ycdragq_x(j)=ycdragh_x(j)*fact_cdrag
	1893	ycdragq_w(j)=ycdragh_w(j)*fact_cdrag
	1894	!! Print *,'YYYYpbl0: fact_cdrag ', fact_cdrag
	1895	ENDDO ! j = 1,knon
	1896	!
	1897	!! Print *,'YYYYpbl0: z1lay, yz0h, f_z0qh_oce, ycdragh_w, ycdragq_w ', &
	1898	!! z1lay, yz0h(1:knon), f_z0qh_oce, ycdragh_w(1:knon), ycdragq_w(1:knon)
	1899	ELSE
	1900	ycdragq_x(1:knon)=ycdragh_x(1:knon)
	1901	ycdragq_w(1:knon)=ycdragh_w(1:knon)
	1902	ENDIF ! ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce)
	1903	!
	1904	CALL wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, ywake_s, &
	1905	yts, y_delta_tsurf, ygustiness, &
[3179]	1906	yt_x, yt_w, yq_x, yq_w, &
	1907	yu_x, yu_w, yv_x, yv_w, &
[4013]	1908	ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
	1909	ycdragm_x, ycdragm_w, &
[3179]	1910	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
	1911	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
	1912	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
	1913	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
[4013]	1914	Kech_h_x, Kech_h_w, Kech_h &
	1915	)
	1916	CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, &
	1917	BcoefQ_x, BcoefQ_w &
	1918	)
	1919	CALL wx_pbl0_merge(knon, ypplay, ypaprs, &
	1920	ywake_s, ydTs0, ydqs0, &
	1921	yt_x, yt_w, yq_x, yq_w, &
	1922	yu_x, yu_w, yv_x, yv_w, &
	1923	ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
	1924	ycdragm_x, ycdragm_w, &
	1925	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
	1926	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
	1927	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
	1928	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
	1929	AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
	1930	BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
	1931	ycdragh, ycdragq, ycdragm, &
[3179]	1932	yt1, yq1, yu1, yv1 &
	1933	)
[4013]	1934	IF (iflag_split .eq. 2 .AND. nsrf .ne. is_oce) THEN
	1935	CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
	1936	ywake_s, ybeta, ywake_cstar, ywake_dens, &
	1937	AcoefH_x, AcoefH_w, &
	1938	BcoefH_x, BcoefH_w, &
	1939	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
	1940	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	1941	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
	1942	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
	1943	yg_T, yg_Q, &
	1944	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
	1945	ydTs_ins, ydqs_ins &
	1946	)
	1947	ELSE !
	1948	AcoefH(:) = AcoefH_0(:)
	1949	AcoefQ(:) = AcoefQ_0(:)
	1950	BcoefH(:) = BcoefH_0(:)
	1951	BcoefQ(:) = BcoefQ_0(:)
	1952	yg_T(:) = 0.
	1953	yg_Q(:) = 0.
	1954	yGamma_dTs_phiT(:) = 0.
	1955	yGamma_dQs_phiQ(:) = 0.
	1956	ydTs_ins(:) = 0.
	1957	ydqs_ins(:) = 0.
	1958	ENDIF ! (iflag_split .eq. 2)
	1959	ENDIF ! (iflag_split .eq.0)
[3179]	1960	!!!
[2159]	1961	IF (prt_level >=10) THEN
[4013]	1962	PRINT *,'pbl_surface (merge->): yt(1,:) ',yt(1,:)
	1963	PRINT *,'pbl_surface (merge->): yq(1,:) ',yq(1,:)
	1964	PRINT *,'pbl_surface (merge->): yu(1,:) ',yu(1,:)
	1965	PRINT *,'pbl_surface (merge->): yv(1,:) ',yv(1,:)
	1966	PRINT *,'pbl_surface (merge->): AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1) ', &
	1967	AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1)
	1968	PRINT *,'pbl_surface (merge->): BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1) ', &
	1969	BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1)
	1970
[2159]	1971	ENDIF
	1972
[4013]	1973	! Save initial value of z0h for use in evappot (z0h wiil be computed again in the surface models)
	1974	yz0h_old(1:knon) = yz0h(1:knon)
	1975	!
[781]	1976	!****************************************************************************************
	1977	!
[1146]	1978	! Calulate t2m and q2m for the case of calculation at land grid points
	1979	! t2m and q2m are needed as input to ORCHIDEE
	1980	!
	1981	!****************************************************************************************
	1982	IF (nsrf == is_ter) THEN
	1983
	1984	DO i = 1, knon
	1985	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
	1986	* (ypaprs(i,1)-ypplay(i,1))
[3391]	1987	ENDDO
[1146]	1988
	1989	! Calculate the temperature et relative humidity at 2m and the wind at 10m
[4013]	1990	IF (iflag_new_t2mq2m==1) THEN
	1991	CALL stdlevvarn(klon, knon, is_ter, zxli, &
	1992	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
	1993	yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
	1994	yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
	1995	yn2mout(:, nsrf, :))
	1996	ELSE
[1146]	1997	CALL stdlevvar(klon, knon, is_ter, zxli, &
	1998	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
[2243]	1999	yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
[1146]	2000	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
[4013]	2001	ENDIF
[1146]	2002
[3391]	2003	ENDIF
[1146]	2004
	2005	!****************************************************************************************
	2006	!
[2159]	2007	! 10) Switch according to current surface
[781]	2008	! It is necessary to start with the continental surfaces because the ocean
	2009	! needs their run-off.
	2010	!
	2011	!****************************************************************************************
	2012	SELECT CASE(nsrf)
	2013
	2014	CASE(is_ter)
[2189]	2015	! print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)
[781]	2016	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
[2410]	2017	rlon, rlat, yrmu0, &
[888]	2018	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
[3179]	2019	!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
	2020	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
[1067]	2021	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	2022	AcoefU, AcoefV, BcoefU, BcoefV, &
[2240]	2023	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
[1146]	2024	ylwdown, yq2m, yt2m, &
[888]	2025	ysnow, yqsol, yagesno, ytsoil, &
[2243]	2026	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
[996]	2027	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
[2952]	2028	y_flux_u1, y_flux_v1, &
	2029	yveget,ylai,yheight )
[3391]	2030
[2952]	2031	!FC quid qd yveget ylai yheight ne sont pas definit
	2032	!FC yveget,ylai,yheight, &
[3391]	2033	IF (ifl_pbltree .ge. 1) THEN
	2034	CALL freinage(knon, yu, yv, yt, &
	2035	! yveget,ylai, yheight,ypaprs,ypplay,y_d_u_frein,y_d_v_frein)
	2036	yveget,ylai, yheight,ypaprs,ypplay,y_treedrg, y_d_u_frein,y_d_v_frein)
	2037	ENDIF
[2952]	2038
[1146]	2039
[2307]	2040	! Special DICE MPL 05082013 puis BOMEX
[2126]	2041	IF (ok_prescr_ust) THEN
[3391]	2042	DO j=1,knon
[2126]	2043	! ysnow(:)=0.
	2044	! yqsol(:)=0.
	2045	! yagesno(:)=50.
	2046	! ytsoil(:,:)=300.
	2047	! yz0_new(:)=0.001
	2048	! yevap(:)=flat/RLVTT
	2049	! yfluxlat(:)=-flat
	2050	! yfluxsens(:)=-fsens
	2051	! yqsurf(:)=0.
	2052	! ytsurf_new(:)=tg
	2053	! y_dflux_t(:)=0.
	2054	! y_dflux_q(:)=0.
[2307]	2055	y_flux_u1(j)=ycdragm(j)(1.+sqrt(yu(j,1)yu(j,1)+yv(j,1)yv(j,1)))yu(j,1)*ypplay(j,1)/RD/yt(j,1)
	2056	y_flux_v1(j)=ycdragm(j)(1.+sqrt(yu(j,1)yu(j,1)+yv(j,1)yv(j,1)))yv(j,1)*ypplay(j,1)/RD/yt(j,1)
[3391]	2057	ENDDO
[2126]	2058	ENDIF
[781]	2059
	2060	CASE(is_lic)
[1865]	2061	! Martin
	2062
[4368]	2063	IF (landice_opt .LT. 2) THEN
	2064	! Land ice is treated by LMDZ and not by ORCHIDEE
	2065	CALL surf_landice(itap, dtime, knon, ni, &
	2066	rlon, rlat, debut, lafin, &
	2067	yrmu0, ylwdown, yalb, zgeo1, &
	2068	ysolsw, ysollw, yts, ypplay(:,1), &
	2069	!!jyg ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
	2070	ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
	2071	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	2072	AcoefU, AcoefV, BcoefU, BcoefV, &
	2073	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
	2074	ysnow, yqsurf, yqsol, yagesno, &
	2075	ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
	2076	ytsurf_new, y_dflux_t, y_dflux_q, &
	2077	yzmea, yzsig, ycldt, &
	2078	ysnowhgt, yqsnow, ytoice, ysissnow, &
	2079	yalb3_new, yrunoff, &
	2080	y_flux_u1, y_flux_v1)
	2081
	2082	!jyg<
	2083	!! alb3_lic(:)=0.
	2084	!>jyg
	2085	DO j = 1, knon
	2086	i = ni(j)
	2087	alb3_lic(i) = yalb3_new(j)
	2088	snowhgt(i) = ysnowhgt(j)
	2089	qsnow(i) = yqsnow(j)
	2090	to_ice(i) = ytoice(j)
	2091	sissnow(i) = ysissnow(j)
	2092	runoff(i) = yrunoff(j)
	2093	ENDDO
	2094	! Martin
	2095	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
	2096	IF (ok_prescr_ust) THEN
	2097	DO j=1,knon
	2098	y_flux_u1(j)=ycdragm(j)(1.+sqrt(yu(j,1)yu(j,1)+yv(j,1)yv(j,1)))yu(j,1)*ypplay(j,1)/RD/yt(j,1)
	2099	y_flux_v1(j)=ycdragm(j)(1.+sqrt(yu(j,1)yu(j,1)+yv(j,1)yv(j,1)))yv(j,1)*ypplay(j,1)/RD/yt(j,1)
	2100	ENDDO
	2101	ENDIF
	2102
	2103	END IF
[781]	2104
	2105	CASE(is_oce)
[2227]	2106	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
[2243]	2107	ywindsp, rmu0, yfder, yts, &
[781]	2108	itap, dtime, jour, knon, ni, &
[3179]	2109	!!jyg ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
[3605]	2110	ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& ! ym missing init
[1067]	2111	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	2112	AcoefU, AcoefV, BcoefU, BcoefV, &
[2240]	2113	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
[888]	2114	ysnow, yqsurf, yagesno, &
[2243]	2115	yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
[1067]	2116	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
[3767]	2117	y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
[3429]	2118	yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
[4020]	2119	ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss)
[2159]	2120	IF (prt_level >=10) THEN
	2121	print *,'arg de surf_ocean: ycdragh ',ycdragh
	2122	print *,'arg de surf_ocean: ycdragm ',ycdragm
	2123	print *,'arg de surf_ocean: yt ', yt
	2124	print *,'arg de surf_ocean: yq ', yq
	2125	print *,'arg de surf_ocean: yts ', yts
	2126	print *,'arg de surf_ocean: AcoefH ',AcoefH
	2127	print *,'arg de surf_ocean: AcoefQ ',AcoefQ
	2128	print *,'arg de surf_ocean: BcoefH ',BcoefH
	2129	print *,'arg de surf_ocean: BcoefQ ',BcoefQ
	2130	print *,'arg de surf_ocean: yevap ',yevap
	2131	print *,'arg de surf_ocean: yfluxsens ',yfluxsens
	2132	print *,'arg de surf_ocean: yfluxlat ',yfluxlat
	2133	print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new
	2134	ENDIF
[2307]	2135	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
	2136	IF (ok_prescr_ust) THEN
	2137	DO j=1,knon
	2138	y_flux_u1(j)=ycdragm(j)(1.+sqrt(yu(j,1)yu(j,1)+yv(j,1)yv(j,1)))yu(j,1)*ypplay(j,1)/RD/yt(j,1)
	2139	y_flux_v1(j)=ycdragm(j)(1.+sqrt(yu(j,1)yu(j,1)+yv(j,1)yv(j,1)))yv(j,1)*ypplay(j,1)/RD/yt(j,1)
	2140	ENDDO
	2141	ENDIF
[781]	2142
	2143	CASE(is_sic)
	2144	CALL surf_seaice( &
[2227]	2145	!albedo SB >>>
	2146	rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, &
	2147	!albedo SB <<<
[781]	2148	itap, dtime, jour, knon, ni, &
[1067]	2149	lafin, &
[3179]	2150	!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
	2151	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
[1067]	2152	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	2153	AcoefU, AcoefV, BcoefU, BcoefV, &
[2244]	2154	ypsref, yu1, yv1, ygustiness, pctsrf, &
[888]	2155	ysnow, yqsurf, yqsol, yagesno, ytsoil, &
[2227]	2156	!albedo SB >>>
[2243]	2157	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
[2227]	2158	!albedo SB <<<
[1067]	2159	ytsurf_new, y_dflux_t, y_dflux_q, &
	2160	y_flux_u1, y_flux_v1)
[781]	2161
[2307]	2162	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
	2163	IF (ok_prescr_ust) THEN
	2164	DO j=1,knon
	2165	y_flux_u1(j)=ycdragm(j)(1.+sqrt(yu(j,1)yu(j,1)+yv(j,1)yv(j,1)))yu(j,1)*ypplay(j,1)/RD/yt(j,1)
	2166	y_flux_v1(j)=ycdragm(j)(1.+sqrt(yu(j,1)yu(j,1)+yv(j,1)yv(j,1)))yv(j,1)*ypplay(j,1)/RD/yt(j,1)
	2167	ENDDO
	2168	ENDIF
[781]	2169
	2170	CASE DEFAULT
	2171	WRITE(lunout,*) 'Surface index = ', nsrf
	2172	abort_message = 'Surface index not valid'
[2311]	2173	CALL abort_physic(modname,abort_message,1)
[781]	2174	END SELECT
	2175
	2176
	2177	!****************************************************************************************
	2178	! 11) - Calcul the increment of surface temperature
	2179	!
	2180	!****************************************************************************************
[1894]	2181
[3392]	2182	IF (evap0>=0.) THEN
[1894]	2183	yevap(:)=evap0
	2184	yevap(:)=RLVTT*evap0
[3391]	2185	ENDIF
[1894]	2186
[781]	2187	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
	2188
	2189	!****************************************************************************************
	2190	!
	2191	! 12) "La remontee" - "The uphill"
	2192	!
	2193	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
	2194	! for X=H, Q, U and V, for all vertical levels.
	2195	!
	2196	!****************************************************************************************
[3798]	2197	!!
	2198	!!!
	2199	!!! jyg le 10/04/2013 et EV 10/2020
[2159]	2200
[3798]	2201	IF (ok_forc_tsurf) THEN
	2202	DO j=1,knon
	2203	ytsurf_new(j)=tg
	2204	y_d_ts(j) = ytsurf_new(j) - yts(j)
	2205	ENDDO
	2206	ENDIF ! ok_forc_tsurf
	2207
[2159]	2208	!!!
	2209	IF (ok_flux_surf) THEN
	2210	IF (prt_level >=10) THEN
	2211	PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
	2212	ENDIF
[882]	2213	y_flux_t1(:) = fsens
	2214	y_flux_q1(:) = flat/RLVTT
	2215	yfluxlat(:) = flat
[2159]	2216	!
[3179]	2217	!! Test sur iflag_split retire le 2/02/2018, sans vraiment comprendre la raison de ce test. (jyg)
	2218	!! IF (iflag_split .eq.0) THEN
[3391]	2219	DO j=1,knon
[2307]	2220	Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)2+yv(j,1)2)) * &
	2221	ypplay(j,1)/(RD*yt(j,1))
[3391]	2222	ENDDO
[3179]	2223	!! ENDIF ! (iflag_split .eq.0)
[1555]	2224
[2159]	2225	DO j = 1, knon
[2452]	2226	yt1_new=(1./RCPD)(AcoefH(j)+BcoefH(j)y_flux_t1(j)*dtime)
	2227	ytsurf_new(j)=yt1_new-y_flux_t1(j)/(Kech_h(j)*RCPD)
[2159]	2228	ENDDO
	2229
[3391]	2230	DO j=1,knon
[2307]	2231	y_d_ts(j) = ytsurf_new(j) - yts(j)
[3391]	2232	ENDDO
[1555]	2233
[2159]	2234	ELSE ! (ok_flux_surf)
[3391]	2235	DO j=1,knon
[2307]	2236	y_flux_t1(j) = yfluxsens(j)
	2237	y_flux_q1(j) = -yevap(j)
[3391]	2238	ENDDO
[4013]	2239	ENDIF ! (ok_flux_surf)
	2240	!
	2241	! ------------------------------------------------------------------------------
	2242	! 12a) Splitting
	2243	! ------------------------------------------------------------------------------
[2159]	2244
[3179]	2245	IF (iflag_split .GE. 1) THEN
[4013]	2246	!
	2247	IF (nsrf .ne. is_oce) THEN
	2248	!
	2249	! Compute potential evaporation and aridity factor (jyg, 20200328)
	2250	ybeta_prev(:) = ybeta(:)
	2251	DO j = 1, knon
	2252	yqa(j) = AcoefQ(j) - BcoefQ(j)yevap(j)dtime
	2253	ENDDO
	2254	!
	2255	CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot)
	2256	!
	2257	ybeta(1:knon) = min(yevap(1:knon)/yevap_pot(1:knon), 1.)
	2258
	2259	IF (prt_level >=10) THEN
	2260	DO j=1,knon
	2261	print*,'y_flux_t1,yfluxlat,wakes' &
	2262	& , y_flux_t1(j), yfluxlat(j), ywake_s(j)
	2263	print*,'beta_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j)
	2264	print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
	2265	ENDDO
	2266	ENDIF ! (prt_level >=10)
	2267	!
	2268	! Second call to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account
	2269	! the update of the aridity coeficient beta.
	2270	!
	2271	CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, &
	2272	BcoefQ_x, BcoefQ_w &
	2273	)
	2274	CALL wx_pbl0_merge(knon, ypplay, ypaprs, &
	2275	ywake_s, ydTs0, ydqs0, &
	2276	yt_x, yt_w, yq_x, yq_w, &
	2277	yu_x, yu_w, yv_x, yv_w, &
	2278	ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
	2279	ycdragm_x, ycdragm_w, &
	2280	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
	2281	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
	2282	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
	2283	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
	2284	AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
	2285	BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
	2286	ycdragh, ycdragq, ycdragm, &
	2287	yt1, yq1, yu1, yv1 &
	2288	)
	2289	IF (iflag_split .eq. 2) THEN
	2290	CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
	2291	ywake_s, ybeta, ywake_cstar, ywake_dens, &
	2292	AcoefH_x, AcoefH_w, &
	2293	BcoefH_x, BcoefH_w, &
	2294	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
	2295	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	2296	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
	2297	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
	2298	yg_T, yg_Q, &
	2299	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
	2300	ydTs_ins, ydqs_ins &
	2301	)
	2302	ELSE !
	2303	AcoefH(:) = AcoefH_0(:)
	2304	AcoefQ(:) = AcoefQ_0(:)
	2305	BcoefH(:) = BcoefH_0(:)
	2306	BcoefQ(:) = BcoefQ_0(:)
	2307	yg_T(:) = 0.
	2308	yg_Q(:) = 0.
	2309	yGamma_dTs_phiT(:) = 0.
	2310	yGamma_dQs_phiQ(:) = 0.
	2311	ydTs_ins(:) = 0.
	2312	ydqs_ins(:) = 0.
	2313	ENDIF ! (iflag_split .eq. 2)
	2314	!
	2315	ELSE ! (nsrf .ne. is_oce)
	2316	ybeta(1:knon) = 1.
	2317	yevap_pot(1:knon) = yevap(1:knon)
	2318	AcoefH(:) = AcoefH_0(:)
	2319	AcoefQ(:) = AcoefQ_0(:)
	2320	BcoefH(:) = BcoefH_0(:)
	2321	BcoefQ(:) = BcoefQ_0(:)
	2322	yg_T(:) = 0.
	2323	yg_Q(:) = 0.
	2324	yGamma_dTs_phiT(:) = 0.
	2325	yGamma_dQs_phiQ(:) = 0.
	2326	ydTs_ins(:) = 0.
	2327	ydqs_ins(:) = 0.
	2328	ENDIF ! (nsrf .ne. is_oce)
	2329	!
	2330	CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, &
	2331	yg_T, yg_Q, &
	2332	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
	2333	ydTs_ins, ydqs_ins, &
[3179]	2334	y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
[4013]	2335	!!!! HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
	2336	phiQ0_b, phiT0_b, &
[3179]	2337	y_flux_t1_x, y_flux_t1_w, &
	2338	y_flux_q1_x, y_flux_q1_w, &
	2339	y_flux_u1_x, y_flux_u1_w, &
	2340	y_flux_v1_x, y_flux_v1_w, &
	2341	yfluxlat_x, yfluxlat_w, &
[4013]	2342	y_delta_qsats, &
	2343	y_delta_tsurf_new, y_delta_qsurf &
[3179]	2344	)
[4013]	2345	!
	2346	CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
	2347	yTs, y_delta_tsurf, &
	2348	yqsurf, yTsurf_new, &
	2349	y_delta_tsurf_new, y_delta_qsats, &
	2350	AcoefH_x, AcoefH_w, &
	2351	BcoefH_x, BcoefH_w, &
	2352	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
	2353	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	2354	y_flux_t1, y_flux_q1, &
	2355	y_flux_t1_x, y_flux_t1_w, &
	2356	y_flux_q1_x, y_flux_q1_w)
	2357	!
	2358	IF (nsrf .ne. is_oce) THEN
	2359	CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
	2360	yTs, y_delta_tsurf, &
	2361	yqsurf, yTsurf_new, &
	2362	y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf, &
	2363	AcoefH_x, AcoefH_w, &
	2364	BcoefH_x, BcoefH_w, &
	2365	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
	2366	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	2367	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
	2368	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
	2369	yg_T, yg_Q, &
	2370	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
	2371	ydTs_ins, ydqs_ins, &
	2372	y_flux_t1, y_flux_q1, &
	2373	y_flux_t1_x, y_flux_t1_w, &
	2374	y_flux_q1_x, y_flux_q1_w )
	2375	ENDIF ! (nsrf .ne. is_oce)
	2376	!
	2377	ELSE ! (iflag_split .ge. 1)
	2378	ybeta(1:knon) = 1.
	2379	yevap_pot(1:knon) = yevap(1:knon)
[3179]	2380	ENDIF ! (iflag_split .ge. 1)
[4013]	2381	!
	2382	IF (prt_level >= 10) THEN
	2383	print *,'pbl_surface, ybeta , yevap, yevap_pot ', &
	2384	ybeta , yevap, yevap_pot
	2385	ENDIF ! (prt_level >= 10)
	2386	!
[3179]	2387	!>jyg
[2159]	2388	!
	2389
	2390	!!jyg!! A reprendre apres reflexion ===============================================
	2391	!!jyg!!
	2392	!!jyg!! DO j=1,knon
	2393	!!jyg!!!!! nrlmd le 13/06/2011
	2394	!!jyg!!
	2395	!!jyg!!!----Diffusion dans le sol dans le cas continental seulement
	2396	!!jyg!! IF (nsrf.eq.is_ter) THEN
	2397	!!jyg!!!----Calcul du coefficient delta_coeff
	2398	!!jyg!! tau_eq(j)=(ywake_s(j)/2.)(1./max(wake_cstar(j),0.01))sqrt(0.4/(3.14*max(wake_dens(j),8e-12)))
	2399	!!jyg!!
[3179]	2400	!!jyg!!! delta_coef(j)=dtime/(inertia*sqrt(tau_eq(j)))
	2401	!!jyg!! delta_coef(j)=facteur*sqrt(tau_eq(j))/inertia
[2159]	2402	!!jyg!!! delta_coef(j)=0.
	2403	!!jyg!! ELSE
	2404	!!jyg!! delta_coef(j)=0.
	2405	!!jyg!! ENDIF
	2406	!!jyg!!
	2407	!!jyg!!!----Calcul de delta_tsurf
	2408	!!jyg!! y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j)
	2409	!!jyg!!
	2410	!!jyg!!!----Si il n'y a pas des poches...
	2411	!!jyg!! IF (wake_cstar(j).le.0.01) THEN
	2412	!!jyg!! y_delta_tsurf(j)=0.
	2413	!!jyg!! y_delta_flux_t1(j)=0.
	2414	!!jyg!! ENDIF
	2415	!!jyg!!
[2450]	2416	!!jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle)
[2159]	2417	!!jyg!!!!!!! jyg le 23/02/2012
	2418	!!jyg!!!!!!!
	2419	!!jyg!!!! ybeta(j)=y_flux_q1(j) / &
	2420	!!jyg!!!! & (Kech_h(j)*(yq(j,1)-yqsatsurf(j)))
	2421	!!jyg!!!!!! ybeta(j)=-1.*yevap(j) / &
	2422	!!jyg!!!!!! & (ywake_s(j)Kech_h_w(j)(yq_w(j,1)-yqsatsurf_w(j))+(1.-ywake_s(j))Kech_h_x(j)(yq_x(j,1)-yqsatsurf_x(j)))
	2423	!!jyg!!!!!!! fin jyg
	2424	!!jyg!!!!!
	2425	!!jyg!!
	2426	!!jyg!! ENDDO
	2427	!!jyg!!
	2428	!!jyg!!!!! fin nrlmd le 13/06/2011
	2429	!!jyg!!
[3179]	2430	IF (iflag_split .ge. 1) THEN
[2159]	2431	IF (prt_level >=10) THEN
	2432	DO j = 1, knon
	2433	print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
	2434	print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
[4013]	2435	print*,'t1x, t1w, t1, t1_ancien', &
	2436	& yt_x(j,1), yt_w(j,1), yt(j,1), t(j,1)
[2159]	2437	print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j)
	2438	ENDDO
	2439
	2440	DO j=1,knon
	2441	print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
	2442	& , y_flux_t1_x(j), y_flux_t1_w(j), y_flux_t1(j), y_flux_q1_x(j)RLVTT, y_flux_q1_w(j)RLVTT, yfluxlat(j), ywake_s(j)
[4013]	2443	print*,'beta, ytsurf_new ', ybeta(j), ytsurf_new(j)
	2444	print*,'inertia, facteur, cstar', inertia, facteur,wake_cstar(j)
[2159]	2445	ENDDO
[2852]	2446	ENDIF ! (prt_level >=10)
[2159]	2447
	2448	!!! jyg le 07/02/2012
[3179]	2449	ENDIF ! (iflag_split .ge.1)
[2159]	2450	!!!
	2451
	2452	!!! jyg le 07/02/2012
	2453	IF (iflag_split .eq.0) THEN
	2454	!!!
	2455	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
	2456	CALL climb_hq_up(knon, dtime, yt, yq, &
[781]	2457	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
[2159]	2458	!!! jyg le 07/02/2012
	2459	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	2460	CcoefH, CcoefQ, DcoefH, DcoefQ, &
	2461	Kcoef_hq, gama_q, gama_h, &
	2462	!!!
[781]	2463	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
[2159]	2464	ELSE !(iflag_split .eq.0)
	2465	CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
	2466	y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &
	2467	!!! nrlmd le 02/05/2011
	2468	AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &
	2469	CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
	2470	Kcoef_hq_x, gama_q_x, gama_h_x, &
	2471	!!!
	2472	y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:))
	2473	!
	2474	CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
	2475	y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &
	2476	!!! nrlmd le 02/05/2011
	2477	AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &
	2478	CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
	2479	Kcoef_hq_w, gama_q_w, gama_h_w, &
	2480	!!!
	2481	y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:))
	2482	!!!
	2483	ENDIF ! (iflag_split .eq.0)
	2484	!!!
[1067]	2485
[2159]	2486	!!! jyg le 07/02/2012
	2487	IF (iflag_split .eq.0) THEN
	2488	!!!
	2489	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
	2490	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
	2491	!!! jyg le 07/02/2012
	2492	AcoefU, AcoefV, BcoefU, BcoefV, &
	2493	CcoefU, CcoefV, DcoefU, DcoefV, &
	2494	Kcoef_m, &
	2495	!!!
[781]	2496	y_flux_u, y_flux_v, y_d_u, y_d_v)
[1761]	2497	y_d_t_diss(:,:)=0.
	2498	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
	2499	CALL yamada_c(knon,dtime,ypaprs,ypplay &
	2500	& ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
[2680]	2501	& ,iflag_pbl)
[1761]	2502	ENDIF
	2503	! print*,'yamada_c OK'
	2504
[2159]	2505	ELSE !(iflag_split .eq.0)
	2506	CALL climb_wind_up(knon, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &
	2507	!!! nrlmd le 02/05/2011
	2508	AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &
	2509	CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
	2510	Kcoef_m_x, &
	2511	!!!
	2512	y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)
	2513	!
	2514	y_d_t_diss_x(:,:)=0.
	2515	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
	2516	CALL yamada_c(knon,dtime,ypaprs,ypplay &
	2517	& ,yu_x,yv_x,yt_x,y_d_u_x,y_d_v_x,y_d_t_x,ycdragm_x,ytke_x,ycoefm_x,ycoefh_x &
	2518	,ycoefq_x,y_d_t_diss_x,yustar_x &
[2680]	2519	& ,iflag_pbl)
[2159]	2520	ENDIF
	2521	! print*,'yamada_c OK'
	2522
	2523	CALL climb_wind_up(knon, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &
	2524	!!! nrlmd le 02/05/2011
	2525	AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &
	2526	CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
	2527	Kcoef_m_w, &
	2528	!!!
	2529	y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)
	2530	!!!
	2531	y_d_t_diss_w(:,:)=0.
	2532	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
	2533	CALL yamada_c(knon,dtime,ypaprs,ypplay &
	2534	& ,yu_w,yv_w,yt_w,y_d_u_w,y_d_v_w,y_d_t_w,ycdragm_w,ytke_w,ycoefm_w,ycoefh_w &
	2535	,ycoefq_w,y_d_t_diss_w,yustar_w &
[2680]	2536	& ,iflag_pbl)
[2159]	2537	ENDIF
	2538	! print*,'yamada_c OK'
	2539	!
	2540	IF (prt_level >=10) THEN
	2541	print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &
	2542	yfluxlat_x, yfluxlat_w
	2543	ENDIF
	2544	!
	2545	ENDIF ! (iflag_split .eq.0)
	2546	!!!
[4013]	2547	!!
	2548	!! DO j = 1, knon
	2549	!! y_dflux_t(j) = y_dflux_t(j) * ypct(j)
	2550	!! y_dflux_q(j) = y_dflux_q(j) * ypct(j)
	2551	!! ENDDO
	2552	!!
[781]	2553	!****************************************************************************************
	2554	! 13) Transform variables for output format :
	2555	! - Decompress
	2556	! - Multiply with pourcentage of current surface
	2557	! - Cumulate in global variable
	2558	!
	2559	!****************************************************************************************
	2560
[2159]	2561
	2562	!!! jyg le 07/02/2012
[3392]	2563	IF (iflag_split.EQ.0) THEN
[2159]	2564	!!!
	2565	DO k = 1, klev
	2566	DO j = 1, knon
[781]	2567	i = ni(j)
[1761]	2568	y_d_t_diss(j,k) = y_d_t_diss(j,k) * ypct(j)
[996]	2569	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
	2570	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
	2571	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
	2572	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
[2952]	2573	!FC
[3392]	2574	IF (nsrf .EQ. is_ter .and. ifl_pbltree .GE. 1) THEN
[2952]	2575	! if (y_d_u_frein(j,k).ne.0. ) then
	2576	! print, nsrf,'IS_TER ++', y_d_u_frein(j,k)ypct(j),y_d_u(j,k),j,k
[3391]	2577	! ENDIF
	2578	y_d_u(j,k) =y_d_u(j,k) + y_d_u_frein(j,k)*ypct(j)
	2579	y_d_v(j,k) =y_d_v(j,k) + y_d_v_frein(j,k)*ypct(j)
	2580	treedrg(i,k,nsrf)=y_treedrg(j,k)
	2581	ELSE
	2582	treedrg(i,k,nsrf)=0.
	2583	ENDIF
[2952]	2584	!FC
[781]	2585	flux_t(i,k,nsrf) = y_flux_t(j,k)
	2586	flux_q(i,k,nsrf) = y_flux_q(j,k)
	2587	flux_u(i,k,nsrf) = y_flux_u(j,k)
	2588	flux_v(i,k,nsrf) = y_flux_v(j,k)
[2159]	2589	ENDDO
	2590	ENDDO
	2591
	2592	ELSE !(iflag_split .eq.0)
	2593
	2594	! Tendances hors poches
	2595	DO k = 1, klev
	2596	DO j = 1, knon
	2597	i = ni(j)
	2598	y_d_t_diss_x(j,k) = y_d_t_diss_x(j,k) * ypct(j)
	2599	y_d_t_x(j,k) = y_d_t_x(j,k) * ypct(j)
	2600	y_d_q_x(j,k) = y_d_q_x(j,k) * ypct(j)
	2601	y_d_u_x(j,k) = y_d_u_x(j,k) * ypct(j)
	2602	y_d_v_x(j,k) = y_d_v_x(j,k) * ypct(j)
	2603
	2604	flux_t_x(i,k,nsrf) = y_flux_t_x(j,k)
	2605	flux_q_x(i,k,nsrf) = y_flux_q_x(j,k)
	2606	flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
	2607	flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)
[781]	2608	ENDDO
[2159]	2609	ENDDO
[1067]	2610
[2159]	2611	! Tendances dans les poches
	2612	DO k = 1, klev
	2613	DO j = 1, knon
	2614	i = ni(j)
	2615	y_d_t_diss_w(j,k) = y_d_t_diss_w(j,k) * ypct(j)
	2616	y_d_t_w(j,k) = y_d_t_w(j,k) * ypct(j)
	2617	y_d_q_w(j,k) = y_d_q_w(j,k) * ypct(j)
	2618	y_d_u_w(j,k) = y_d_u_w(j,k) * ypct(j)
	2619	y_d_v_w(j,k) = y_d_v_w(j,k) * ypct(j)
	2620
	2621	flux_t_w(i,k,nsrf) = y_flux_t_w(j,k)
	2622	flux_q_w(i,k,nsrf) = y_flux_q_w(j,k)
	2623	flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
	2624	flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)
	2625	ENDDO
	2626	ENDDO
	2627
	2628	! Flux, tendances et Tke moyenne dans la maille
	2629	DO k = 1, klev
	2630	DO j = 1, knon
	2631	i = ni(j)
	2632	flux_t(i,k,nsrf) = flux_t_x(i,k,nsrf)+ywake_s(j)*(flux_t_w(i,k,nsrf)-flux_t_x(i,k,nsrf))
	2633	flux_q(i,k,nsrf) = flux_q_x(i,k,nsrf)+ywake_s(j)*(flux_q_w(i,k,nsrf)-flux_q_x(i,k,nsrf))
	2634	flux_u(i,k,nsrf) = flux_u_x(i,k,nsrf)+ywake_s(j)*(flux_u_w(i,k,nsrf)-flux_u_x(i,k,nsrf))
	2635	flux_v(i,k,nsrf) = flux_v_x(i,k,nsrf)+ywake_s(j)*(flux_v_w(i,k,nsrf)-flux_v_x(i,k,nsrf))
	2636	ENDDO
	2637	ENDDO
	2638	DO j=1,knon
	2639	yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j))
	2640	ENDDO
	2641	IF (prt_level >=10) THEN
	2642	print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &
	2643	nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf)
	2644	ENDIF
	2645
	2646	DO k = 1, klev
	2647	DO j = 1, knon
	2648	y_d_t_diss(j,k) = y_d_t_diss_x(j,k)+ywake_s(j)*(y_d_t_diss_w(j,k) -y_d_t_diss_x(j,k))
	2649	y_d_t(j,k) = y_d_t_x(j,k)+ywake_s(j)*(y_d_t_w(j,k) -y_d_t_x(j,k))
	2650	y_d_q(j,k) = y_d_q_x(j,k)+ywake_s(j)*(y_d_q_w(j,k) -y_d_q_x(j,k))
	2651	y_d_u(j,k) = y_d_u_x(j,k)+ywake_s(j)*(y_d_u_w(j,k) -y_d_u_x(j,k))
	2652	y_d_v(j,k) = y_d_v_x(j,k)+ywake_s(j)*(y_d_v_w(j,k) -y_d_v_x(j,k))
	2653	ENDDO
	2654	ENDDO
	2655
	2656	ENDIF ! (iflag_split .eq.0)
	2657	!!!
	2658
[1761]	2659	! print*,'Dans pbl OK1'
	2660
[2181]	2661	!jyg<
	2662	!! evap(:,nsrf) = - flux_q(:,1,nsrf)
	2663	!>jyg
[781]	2664	DO j = 1, knon
	2665	i = ni(j)
[2181]	2666	evap(i,nsrf) = - flux_q(i,1,nsrf) !jyg
[4013]	2667	beta(i,nsrf) = ybeta(j) !jyg
[781]	2668	d_ts(i,nsrf) = y_d_ts(j)
[2227]	2669	!albedo SB >>>
[3391]	2670	DO k=1,nsw
	2671	alb_dir(i,k,nsrf) = yalb_dir_new(j,k)
	2672	alb_dif(i,k,nsrf) = yalb_dif_new(j,k)
	2673	ENDDO
[2227]	2674	!albedo SB <<<
[781]	2675	snow(i,nsrf) = ysnow(j)
	2676	qsurf(i,nsrf) = yqsurf(j)
[2243]	2677	z0m(i,nsrf) = yz0m(j)
	2678	z0h(i,nsrf) = yz0h(j)
[781]	2679	fluxlat(i,nsrf) = yfluxlat(j)
	2680	agesno(i,nsrf) = yagesno(j)
[1067]	2681	cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
	2682	cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
[4013]	2683	dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)
	2684	dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)
[3391]	2685	ENDDO
[781]	2686
[1761]	2687	! print*,'Dans pbl OK2'
	2688
[2159]	2689	!!! jyg le 07/02/2012
[2852]	2690	IF (iflag_split .ge.1) THEN
[2159]	2691	!!!
	2692	!!! nrlmd le 02/05/2011
	2693	DO j = 1, knon
	2694	i = ni(j)
	2695	fluxlat_x(i,nsrf) = yfluxlat_x(j)
	2696	fluxlat_w(i,nsrf) = yfluxlat_w(j)
	2697	!!!
	2698	!!! nrlmd le 13/06/2011
[3179]	2699	!!jyg20170131 delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
[4013]	2700	!!jyg20210118 delta_tsurf(i,nsrf)=y_delta_tsurf(j)
	2701	delta_tsurf(i,nsrf)=y_delta_tsurf_new(j)
[3179]	2702	!
[4013]	2703	delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j)*ypct(j)
	2704	!
[2159]	2705	cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
	2706	cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)
	2707	cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j)
	2708	cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j)
	2709	kh(i) = kh(i) + Kech_h(j)*ypct(j)
	2710	kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j)
	2711	kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j)
	2712	!!!
[3391]	2713	ENDDO
[2159]	2714	!!!
[2852]	2715	ENDIF ! (iflag_split .ge.1)
[2159]	2716	!!!
	2717	!!! nrlmd le 02/05/2011
	2718	!!jyg le 20/02/2011
	2719	!! tke_x(:,:,nsrf)=0.
	2720	!! tke_w(:,:,nsrf)=0.
	2721	!!jyg le 20/02/2011
	2722	!! DO k = 1, klev+1
	2723	!! DO j = 1, knon
	2724	!! i = ni(j)
	2725	!! wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
	2726	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
	2727	!! ENDDO
	2728	!! ENDDO
	2729	!!jyg le 20/02/2011
	2730	!! DO k = 1, klev+1
	2731	!! DO j = 1, knon
	2732	!! i = ni(j)
	2733	!! tke(i,k,nsrf)=(1.-ywake_s(j))tke_x(i,k,nsrf)+ywake_s(j)tke_w(i,k,nsrf)
	2734	!! ENDDO
	2735	!! ENDDO
	2736	!!!
	2737	IF (iflag_split .eq.0) THEN
[2884]	2738	wake_dltke(:,:,nsrf) = 0.
[3198]	2739	DO k = 1, klev+1
[2159]	2740	DO j = 1, knon
	2741	i = ni(j)
[2181]	2742	!jyg<
	2743	!! tke(i,k,nsrf) = ytke(j,k)
	2744	!! tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
	2745	tke_x(i,k,nsrf) = ytke(j,k)
	2746	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
[3798]	2747
[2181]	2748	!>jyg
[3391]	2749	ENDDO
	2750	ENDDO
[2159]	2751
[2852]	2752	ELSE ! (iflag_split .eq.0)
[3198]	2753	DO k = 1, klev+1
[2159]	2754	DO j = 1, knon
	2755	i = ni(j)
	2756	wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
[2181]	2757	!jyg<
	2758	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
	2759	!! tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
	2760	tke_x(i,k,nsrf) = ytke_x(j,k)
[3798]	2761	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)
[2181]	2762	wake_dltke(i,k,is_ave) = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
[3798]	2763
[2181]	2764
	2765	!>jyg
[2159]	2766	ENDDO
	2767	ENDDO
	2768	ENDIF ! (iflag_split .eq.0)
	2769	!!!
[1067]	2770	DO k = 2, klev
	2771	DO j = 1, knon
	2772	i = ni(j)
[1761]	2773	zcoefh(i,k,nsrf) = ycoefh(j,k)
	2774	zcoefm(i,k,nsrf) = ycoefm(j,k)
	2775	zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
	2776	zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
[3391]	2777	ENDDO
	2778	ENDDO
[1067]	2779
[1761]	2780	! print*,'Dans pbl OK3'
	2781
[781]	2782	IF ( nsrf .EQ. is_ter ) THEN
	2783	DO j = 1, knon
	2784	i = ni(j)
	2785	qsol(i) = yqsol(j)
[3391]	2786	ENDDO
	2787	ENDIF
[781]	2788
[2181]	2789	!jyg<
	2790	!! ftsoil(:,:,nsrf) = 0.
	2791	!>jyg
[781]	2792	DO k = 1, nsoilmx
	2793	DO j = 1, knon
	2794	i = ni(j)
	2795	ftsoil(i, k, nsrf) = ytsoil(j,k)
[3391]	2796	ENDDO
	2797	ENDDO
[781]	2798
[2159]	2799	!!! jyg le 07/02/2012
[2852]	2800	IF (iflag_split .ge.1) THEN
[2159]	2801	!!!
	2802	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
	2803	DO k = 1, klev
	2804	DO j = 1, knon
	2805	i = ni(j)
	2806	d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)
	2807	d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)
	2808	d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)
	2809	d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)
	2810	d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)
	2811	!
	2812	d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)
	2813	d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)
	2814	d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)
	2815	d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
	2816	d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
	2817	!
	2818	!! d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
	2819	!! d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)
[3391]	2820	ENDDO
	2821	ENDDO
[2159]	2822	!!!
[2852]	2823	ENDIF ! (iflag_split .ge.1)
[2159]	2824	!!!
[781]	2825
	2826	DO k = 1, klev
	2827	DO j = 1, knon
	2828	i = ni(j)
[1761]	2829	d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
[781]	2830	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
	2831	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
	2832	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
	2833	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
[3391]	2834	ENDDO
	2835	ENDDO
[781]	2836
[1761]	2837	! print*,'Dans pbl OK4'
	2838
[2159]	2839	IF (prt_level >=10) THEN
	2840	print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
	2841	d_t_w(:,1), d_t_x(:,1), d_t(:,1)
	2842	ENDIF
	2843
[3458]	2844	if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
[3767]	2845	delta_sal = missing_val
[3628]	2846	ds_ns = missing_val
[3626]	2847	dt_ns = missing_val
[3744]	2848	delta_sst = missing_val
[3628]	2849	dter = missing_val
	2850	dser = missing_val
	2851	tkt = missing_val
	2852	tks = missing_val
	2853	taur = missing_val
[3720]	2854	sss = missing_val
[3628]	2855
[3767]	2856	delta_sal(ni(:knon)) = ydelta_sal(:knon)
[3429]	2857	ds_ns(ni(:knon)) = yds_ns(:knon)
	2858	dt_ns(ni(:knon)) = ydt_ns(:knon)
[3744]	2859	delta_sst(ni(:knon)) = ydelta_sst(:knon)
[3429]	2860	dter(ni(:knon)) = ydter(:knon)
	2861	dser(ni(:knon)) = ydser(:knon)
	2862	tkt(ni(:knon)) = ytkt(:knon)
	2863	tks(ni(:knon)) = ytks(:knon)
	2864	taur(ni(:knon)) = ytaur(:knon)
[3720]	2865	sss(ni(:knon)) = ysss(:knon)
[4020]	2866
	2867	if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
	2868	dt_ds = missing_val
	2869	dt_ds(ni(:knon)) = ydt_ds(:knon)
	2870	end if
[3429]	2871	end if
	2872
[4013]	2873
[781]	2874	!****************************************************************************************
[2159]	2875	! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m
[781]	2876	! Call HBTM
	2877	!
	2878	!****************************************************************************************
[2159]	2879	!!!
	2880	!
[781]	2881	#undef T2m
	2882	#define T2m
	2883	#ifdef T2m
[996]	2884	! Calculations of diagnostic t,q at 2m and u, v at 10m
[781]	2885
[1761]	2886	! print*,'Dans pbl OK41'
	2887	! print*,'tair1,yt(:,1),y_d_t(:,1)'
	2888	! print*, tair1,yt(:,1),y_d_t(:,1)
[2159]	2889	!!! jyg le 07/02/2012
	2890	IF (iflag_split .eq.0) THEN
	2891	DO j=1, knon
[781]	2892	uzon(j) = yu(j,1) + y_d_u(j,1)
	2893	vmer(j) = yv(j,1) + y_d_v(j,1)
[1761]	2894	tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
[781]	2895	qair1(j) = yq(j,1) + y_d_q(j,1)
	2896	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
	2897	* (ypaprs(j,1)-ypplay(j,1))
	2898	tairsol(j) = yts(j) + y_d_ts(j)
[2159]	2899	qairsol(j) = yqsurf(j)
[3391]	2900	ENDDO
[2159]	2901	ELSE ! (iflag_split .eq.0)
	2902	DO j=1, knon
	2903	uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
	2904	vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
	2905	tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
	2906	qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
	2907	zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
	2908	* (ypaprs(j,1)-ypplay(j,1))
	2909	tairsol(j) = yts(j) + y_d_ts(j)
[4013]	2910	!! tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
	2911	tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf_new(j)
[2159]	2912	qairsol(j) = yqsurf(j)
[3391]	2913	ENDDO
[2159]	2914	DO j=1, knon
	2915	uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
	2916	vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
	2917	tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
	2918	qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
	2919	zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
	2920	* (ypaprs(j,1)-ypplay(j,1))
	2921	tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
	2922	qairsol(j) = yqsurf(j)
[3391]	2923	ENDDO
[2159]	2924	!!!
	2925	ENDIF ! (iflag_split .eq.0)
	2926	!!!
	2927	DO j=1, knon
[2896]	2928	! i = ni(j)
	2929	! yz0h_oupas(j) = yz0m(j)
	2930	! IF(nsrf.EQ.is_oce) THEN
	2931	! yz0h_oupas(j) = z0m(i,nsrf)
	2932	! ENDIF
[781]	2933	psfce(j)=ypaprs(j,1)
	2934	patm(j)=ypplay(j,1)
[3391]	2935	ENDDO
[2896]	2936
	2937	IF (iflag_pbl_surface_t2m_bug==1) THEN
	2938	yz0h_oupas(1:knon)=yz0m(1:knon)
	2939	ELSE
	2940	yz0h_oupas(1:knon)=yz0h(1:knon)
	2941	ENDIF
[781]	2942
[1761]	2943	! print*,'Dans pbl OK42A'
	2944	! print*,'tair1,yt(:,1),y_d_t(:,1)'
	2945	! print*, tair1,yt(:,1),y_d_t(:,1)
[781]	2946
[3179]	2947	! Calculate the temperature and relative humidity at 2m and the wind at 10m
[2159]	2948	!!! jyg le 07/02/2012
	2949	IF (iflag_split .eq.0) THEN
[4013]	2950	IF (iflag_new_t2mq2m==1) THEN
	2951	CALL stdlevvarn(klon, knon, nsrf, zxli, &
	2952	uzon, vmer, tair1, qair1, zgeo1, &
	2953	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
	2954	yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
	2955	yn2mout(:, nsrf, :))
	2956	ELSE
[2159]	2957	CALL stdlevvar(klon, knon, nsrf, zxli, &
[781]	2958	uzon, vmer, tair1, qair1, zgeo1, &
[2896]	2959	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
[781]	2960	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
[4013]	2961	ENDIF
[2159]	2962	ELSE !(iflag_split .eq.0)
[4013]	2963	IF (iflag_new_t2mq2m==1) THEN
	2964	CALL stdlevvarn(klon, knon, nsrf, zxli, &
	2965	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
	2966	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
	2967	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, &
	2968	yn2mout_x(:, nsrf, :))
	2969	CALL stdlevvarn(klon, knon, nsrf, zxli, &
	2970	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
	2971	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
	2972	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, &
	2973	yn2mout_w(:, nsrf, :))
	2974	ELSE
[2159]	2975	CALL stdlevvar(klon, knon, nsrf, zxli, &
	2976	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
[2896]	2977	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
[2159]	2978	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
	2979	CALL stdlevvar(klon, knon, nsrf, zxli, &
	2980	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
[2896]	2981	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
[2159]	2982	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
[4013]	2983	ENDIF
[2159]	2984	!!!
	2985	ENDIF ! (iflag_split .eq.0)
	2986	!!!
	2987	!!! jyg le 07/02/2012
	2988	IF (iflag_split .eq.0) THEN
	2989	DO j=1, knon
[781]	2990	i = ni(j)
	2991	t2m(i,nsrf)=yt2m(j)
[996]	2992	q2m(i,nsrf)=yq2m(j)
[2159]	2993	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
[1670]	2994	ustar(i,nsrf)=yustar(j)
[781]	2995	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
	2996	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
[4013]	2997	!
	2998	DO k = 1, 6
	2999	n2mout(i,nsrf,k) = yn2mout(j,nsrf,k)
	3000	END DO
	3001	!
[3391]	3002	ENDDO
[2159]	3003	ELSE !(iflag_split .eq.0)
	3004	DO j=1, knon
	3005	i = ni(j)
	3006	t2m_x(i,nsrf)=yt2m_x(j)
	3007	q2m_x(i,nsrf)=yq2m_x(j)
	3008	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
	3009	ustar_x(i,nsrf)=yustar_x(j)
	3010	u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
	3011	v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
[4013]	3012	!
	3013	DO k = 1, 6
	3014	n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k)
	3015	END DO
	3016	!
[3391]	3017	ENDDO
[2159]	3018	DO j=1, knon
	3019	i = ni(j)
	3020	t2m_w(i,nsrf)=yt2m_w(j)
	3021	q2m_w(i,nsrf)=yq2m_w(j)
	3022	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
	3023	ustar_w(i,nsrf)=yustar_w(j)
	3024	u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
	3025	v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
	3026	!
	3027	ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
	3028	u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
	3029	v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
[4013]	3030	!
	3031	DO k = 1, 6
	3032	n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k)
	3033	END DO
	3034	!
[3391]	3035	ENDDO
[2159]	3036	!!!
	3037	ENDIF ! (iflag_split .eq.0)
	3038	!!!
[1670]	3039
[1761]	3040	! print*,'Dans pbl OK43'
[996]	3041	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
	3042	!IM Ajoute dependance type surface
	3043	IF (thermcep) THEN
[2159]	3044	!!! jyg le 07/02/2012
	3045	IF (iflag_split .eq.0) THEN
[996]	3046	DO j = 1, knon
	3047	i=ni(j)
	3048	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
	3049	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
	3050	zx_qs1 = MIN(0.5,zx_qs1)
	3051	zcor1 = 1./(1.-RETV*zx_qs1)
	3052	zx_qs1 = zx_qs1*zcor1
	3053
	3054	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
	3055	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
[3391]	3056	ENDDO
[2159]	3057	ELSE ! (iflag_split .eq.0)
	3058	DO j = 1, knon
	3059	i=ni(j)
	3060	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
	3061	zx_qs1 = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)
	3062	zx_qs1 = MIN(0.5,zx_qs1)
	3063	zcor1 = 1./(1.-RETV*zx_qs1)
	3064	zx_qs1 = zx_qs1*zcor1
	3065
	3066	rh2m_x(i) = rh2m_x(i) + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
	3067	qsat2m_x(i) = qsat2m_x(i) + zx_qs1 * pctsrf(i,nsrf)
[3391]	3068	ENDDO
[2159]	3069	DO j = 1, knon
	3070	i=ni(j)
	3071	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
	3072	zx_qs1 = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
	3073	zx_qs1 = MIN(0.5,zx_qs1)
	3074	zcor1 = 1./(1.-RETV*zx_qs1)
	3075	zx_qs1 = zx_qs1*zcor1
	3076
	3077	rh2m_w(i) = rh2m_w(i) + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
	3078	qsat2m_w(i) = qsat2m_w(i) + zx_qs1 * pctsrf(i,nsrf)
[3391]	3079	ENDDO
[2159]	3080	!!!
	3081	ENDIF ! (iflag_split .eq.0)
	3082	!!!
[3391]	3083	ENDIF
[2159]	3084	!
	3085	IF (prt_level >=10) THEN
	3086	print *, 'T2m, q2m, RH2m ', &
	3087	t2m, q2m, rh2m
	3088	ENDIF
[781]	3089
[1761]	3090	! print*,'OK pbl 5'
[2159]	3091	!
	3092	!!! jyg le 07/02/2012
	3093	IF (iflag_split .eq.0) THEN
	3094	CALL hbtm(knon, ypaprs, ypplay, &
[1816]	3095	yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
[781]	3096	y_flux_t,y_flux_q,yu,yv,yt,yq, &
	3097	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
	3098	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
[2159]	3099	IF (prt_level >=10) THEN
	3100	print *,' Arg. de HBTM: yt2m ',yt2m
	3101	print *,' Arg. de HBTM: yt10m ',yt10m
	3102	print *,' Arg. de HBTM: yq2m ',yq2m
	3103	print *,' Arg. de HBTM: yq10m ',yq10m
	3104	print *,' Arg. de HBTM: yustar ',yustar
	3105	print *,' Arg. de HBTM: y_flux_t ',y_flux_t
	3106	print *,' Arg. de HBTM: y_flux_q ',y_flux_q
	3107	print *,' Arg. de HBTM: yu ',yu
	3108	print *,' Arg. de HBTM: yv ',yv
	3109	print *,' Arg. de HBTM: yt ',yt
	3110	print *,' Arg. de HBTM: yq ',yq
	3111	ENDIF
	3112	ELSE ! (iflag_split .eq.0)
	3113	CALL HBTM(knon, ypaprs, ypplay, &
	3114	yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
	3115	y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
	3116	ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
	3117	ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
	3118	IF (prt_level >=10) THEN
	3119	print *,' Arg. de HBTM: yt2m_x ',yt2m_x
	3120	print *,' Arg. de HBTM: yt10m_x ',yt10m_x
	3121	print *,' Arg. de HBTM: yq2m_x ',yq2m_x
	3122	print *,' Arg. de HBTM: yq10m_x ',yq10m_x
	3123	print *,' Arg. de HBTM: yustar_x ',yustar_x
	3124	print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x
	3125	print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x
	3126	print *,' Arg. de HBTM: yu_x ',yu_x
	3127	print *,' Arg. de HBTM: yv_x ',yv_x
	3128	print *,' Arg. de HBTM: yt_x ',yt_x
	3129	print *,' Arg. de HBTM: yq_x ',yq_x
	3130	ENDIF
	3131	CALL HBTM(knon, ypaprs, ypplay, &
	3132	yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
	3133	y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
	3134	ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
	3135	ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
	3136	!!!
	3137	ENDIF ! (iflag_split .eq.0)
	3138	!!!
[781]	3139
[2159]	3140	!!! jyg le 07/02/2012
	3141	IF (iflag_split .eq.0) THEN
	3142	!!!
	3143	DO j=1, knon
[781]	3144	i = ni(j)
	3145	pblh(i,nsrf) = ypblh(j)
[1816]	3146	wstar(i,nsrf) = ywstar(j)
[781]	3147	plcl(i,nsrf) = ylcl(j)
	3148	capCL(i,nsrf) = ycapCL(j)
	3149	oliqCL(i,nsrf) = yoliqCL(j)
	3150	cteiCL(i,nsrf) = ycteiCL(j)
	3151	pblT(i,nsrf) = ypblT(j)
	3152	therm(i,nsrf) = ytherm(j)
	3153	trmb1(i,nsrf) = ytrmb1(j)
	3154	trmb2(i,nsrf) = ytrmb2(j)
	3155	trmb3(i,nsrf) = ytrmb3(j)
[3391]	3156	ENDDO
[2159]	3157	IF (prt_level >=10) THEN
	3158	print *, 'After HBTM: pblh ', pblh
	3159	print *, 'After HBTM: plcl ', plcl
	3160	print *, 'After HBTM: cteiCL ', cteiCL
	3161	ENDIF
	3162	ELSE !(iflag_split .eq.0)
	3163	DO j=1, knon
	3164	i = ni(j)
	3165	pblh_x(i,nsrf) = ypblh_x(j)
	3166	wstar_x(i,nsrf) = ywstar_x(j)
	3167	plcl_x(i,nsrf) = ylcl_x(j)
	3168	capCL_x(i,nsrf) = ycapCL_x(j)
	3169	oliqCL_x(i,nsrf) = yoliqCL_x(j)
	3170	cteiCL_x(i,nsrf) = ycteiCL_x(j)
	3171	pblT_x(i,nsrf) = ypblT_x(j)
	3172	therm_x(i,nsrf) = ytherm_x(j)
	3173	trmb1_x(i,nsrf) = ytrmb1_x(j)
	3174	trmb2_x(i,nsrf) = ytrmb2_x(j)
	3175	trmb3_x(i,nsrf) = ytrmb3_x(j)
[3391]	3176	ENDDO
[2159]	3177	IF (prt_level >=10) THEN
	3178	print *, 'After HBTM: pblh_x ', pblh_x
	3179	print *, 'After HBTM: plcl_x ', plcl_x
	3180	print *, 'After HBTM: cteiCL_x ', cteiCL_x
	3181	ENDIF
	3182	DO j=1, knon
	3183	i = ni(j)
	3184	pblh_w(i,nsrf) = ypblh_w(j)
	3185	wstar_w(i,nsrf) = ywstar_w(j)
	3186	plcl_w(i,nsrf) = ylcl_w(j)
	3187	capCL_w(i,nsrf) = ycapCL_w(j)
	3188	oliqCL_w(i,nsrf) = yoliqCL_w(j)
	3189	cteiCL_w(i,nsrf) = ycteiCL_w(j)
	3190	pblT_w(i,nsrf) = ypblT_w(j)
	3191	therm_w(i,nsrf) = ytherm_w(j)
	3192	trmb1_w(i,nsrf) = ytrmb1_w(j)
	3193	trmb2_w(i,nsrf) = ytrmb2_w(j)
	3194	trmb3_w(i,nsrf) = ytrmb3_w(j)
[3391]	3195	ENDDO
[2159]	3196	IF (prt_level >=10) THEN
	3197	print *, 'After HBTM: pblh_w ', pblh_w
	3198	print *, 'After HBTM: plcl_w ', plcl_w
	3199	print *, 'After HBTM: cteiCL_w ', cteiCL_w
	3200	ENDIF
	3201	!!!
	3202	ENDIF ! (iflag_split .eq.0)
	3203	!!!
	3204
[1761]	3205	! print*,'OK pbl 6'
[781]	3206	#else
[996]	3207	! T2m not defined
[781]	3208	! No calculation
[996]	3209	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
[781]	3210	#endif
	3211
	3212	!****************************************************************************************
	3213	! 15) End of loop over different surfaces
	3214	!
	3215	!****************************************************************************************
[3391]	3216	ENDDO loop_nbsrf
[4013]	3217	!
	3218	!----------------------------------------------------------------------------------------
	3219	! Reset iflag_split
	3220	!
	3221	iflag_split=iflag_split_ref
[781]	3222
	3223	!****************************************************************************************
[2159]	3224	! 16) Calculate the mean value over all sub-surfaces for some variables
[781]	3225	!
	3226	!****************************************************************************************
	3227
[2243]	3228	z0m(:,nbsrf+1) = 0.0
	3229	z0h(:,nbsrf+1) = 0.0
	3230	DO nsrf = 1, nbsrf
	3231	DO i = 1, klon
	3232	z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)
	3233	z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)
	3234	ENDDO
	3235	ENDDO
	3236
[1761]	3237	! print*,'OK pbl 7'
[781]	3238	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
	3239	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
[2159]	3240	zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
	3241	zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
	3242	zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
	3243	zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
	3244
	3245	!!! jyg le 07/02/2012
[2852]	3246	IF (iflag_split .ge.1) THEN
[2159]	3247	!!!
	3248	!!! nrlmd & jyg les 02/05/2011, 05/02/2012
	3249
	3250	DO nsrf = 1, nbsrf
	3251	DO k = 1, klev
	3252	DO i = 1, klon
	3253	zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
	3254	zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
	3255	zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
	3256	zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
	3257	!
	3258	zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
	3259	zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
	3260	zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
	3261	zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
[3391]	3262	ENDDO
	3263	ENDDO
	3264	ENDDO
[2159]	3265
	3266	DO i = 1, klon
	3267	zxsens_x(i) = - zxfluxt_x(i,1)
	3268	zxsens_w(i) = - zxfluxt_w(i,1)
[3391]	3269	ENDDO
[2159]	3270	!!!
[2852]	3271	ENDIF ! (iflag_split .ge.1)
[2159]	3272	!!!
	3273
[781]	3274	DO nsrf = 1, nbsrf
	3275	DO k = 1, klev
	3276	DO i = 1, klon
[996]	3277	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
	3278	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
	3279	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
	3280	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
[3391]	3281	ENDDO
	3282	ENDDO
	3283	ENDDO
[781]	3284
	3285	DO i = 1, klon
	3286	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
	3287	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
	3288	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
	3289	ENDDO
[2159]	3290	!!!
[4013]	3291
[781]	3292	!
	3293	! Incrementer la temperature du sol
	3294	!
	3295	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
[4013]	3296	zt2m(:) = 0.0 ; zq2m(:) = 0.0 ; zn2mout(:,:) = 0
[1670]	3297	zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
[781]	3298	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
[2159]	3299	!!! jyg le 07/02/2012
	3300	s_pblh_x(:) = 0.0 ; s_plcl_x(:) = 0.0
	3301	s_pblh_w(:) = 0.0 ; s_plcl_w(:) = 0.0
	3302	!!!
[781]	3303	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
	3304	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
	3305	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
	3306	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
[1816]	3307	wstar(:,is_ave)=0.
[781]	3308
[1761]	3309	! print*,'OK pbl 9'
[781]	3310
[2159]	3311	!!! nrlmd le 02/05/2011
	3312	zxfluxlat_x(:) = 0.0 ; zxfluxlat_w(:) = 0.0
	3313	!!!
	3314
[781]	3315	DO nsrf = 1, nbsrf
	3316	DO i = 1, klon
	3317	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
	3318
	3319	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
[996]	3320	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
[2670]	3321	wfbilo(i,nsrf) = (evap(i,nsrf)-(rain_f(i)+snow_f(i)))*pctsrf(i,nsrf)
	3322	wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)
	3323	wfrain(i,nsrf) = rain_f(i)*pctsrf(i,nsrf)
	3324	wfsnow(i,nsrf) = snow_f(i)*pctsrf(i,nsrf)
[781]	3325
[996]	3326	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
	3327	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
[3391]	3328	ENDDO
	3329	ENDDO
[2455]	3330	!
	3331	!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
	3332	IF (iflag_order2_sollw == 1) THEN
	3333	meansqT(:) = 0. ! as working buffer
	3334	DO nsrf = 1, nbsrf
	3335	DO i = 1, klon
	3336	meansqT(i) = meansqT(i)+(ts(i,nsrf)-zxtsol(i))*2 pctsrf(i,nsrf)
[3391]	3337	ENDDO
	3338	ENDDO
[2455]	3339	zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)
	3340	ENDIF ! iflag_order2_sollw == 1
	3341	!>al1
[781]	3342
[2159]	3343	!!! jyg le 07/02/2012
	3344	IF (iflag_split .eq.0) THEN
	3345	DO nsrf = 1, nbsrf
	3346	DO i = 1, klon
[996]	3347	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
	3348	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
[4013]	3349	!
	3350	DO k = 1, 6
	3351	zn2mout(i,k) = zn2mout(i,k) + n2mout(i,nsrf,k) * pctsrf(i,nsrf)
	3352	ENDDO
	3353	!
[1670]	3354	zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
[1816]	3355	wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
[996]	3356	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
	3357	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
[781]	3358
[996]	3359	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
	3360	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
	3361	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
	3362	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
	3363	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
	3364	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
	3365	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
	3366	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
	3367	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
	3368	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
[3391]	3369	ENDDO
	3370	ENDDO
[2159]	3371	ELSE !(iflag_split .eq.0)
	3372	DO nsrf = 1, nbsrf
	3373	DO i = 1, klon
	3374	!!! nrlmd le 02/05/2011
	3375	zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
	3376	zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
	3377	!!!
	3378	!!! jyg le 08/02/2012
	3379	!! Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;
	3380	!! pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
	3381	!! pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
	3382	!! pour les autres variables, on sort les valeurs de la region (x).
	3383	zt2m(i) = zt2m(i) + (t2m_x(i,nsrf)+wake_s(i)(t2m_w(i,nsrf)-t2m_x(i,nsrf))) pctsrf(i,nsrf)
	3384	zq2m(i) = zq2m(i) + q2m_x(i,nsrf) * pctsrf(i,nsrf)
[4013]	3385	!
	3386	DO k = 1, 6
	3387	zn2mout(i,k) = zn2mout(i,k) + n2mout_x(i,nsrf,k) * pctsrf(i,nsrf)
	3388	ENDDO
	3389	!
[2159]	3390	zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
	3391	wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
	3392	zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
	3393	zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
	3394	!
	3395	s_pblh(i) = s_pblh(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
	3396	s_pblh_x(i) = s_pblh_x(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
	3397	s_pblh_w(i) = s_pblh_w(i) + pblh_w(i,nsrf) * pctsrf(i,nsrf)
	3398	!
	3399	s_plcl(i) = s_plcl(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
	3400	s_plcl_x(i) = s_plcl_x(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
	3401	s_plcl_w(i) = s_plcl_w(i) + plcl_w(i,nsrf) * pctsrf(i,nsrf)
	3402	!
	3403	s_capCL(i) = s_capCL(i) + capCL_x(i,nsrf) * pctsrf(i,nsrf)
	3404	s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
	3405	s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
	3406	s_pblT(i) = s_pblT(i) + pblT_x(i,nsrf) * pctsrf(i,nsrf)
	3407	s_therm(i) = s_therm(i) + therm_x(i,nsrf) * pctsrf(i,nsrf)
	3408	s_trmb1(i) = s_trmb1(i) + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
	3409	s_trmb2(i) = s_trmb2(i) + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
	3410	s_trmb3(i) = s_trmb3(i) + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
[3391]	3411	ENDDO
	3412	ENDDO
[2159]	3413	DO i = 1, klon
	3414	qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
[3391]	3415	ENDDO
[2159]	3416	!!!
	3417	ENDIF ! (iflag_split .eq.0)
	3418	!!!
[781]	3419
	3420	IF (check) THEN
	3421	amn=MIN(ts(1,is_ter),1000.)
	3422	amx=MAX(ts(1,is_ter),-1000.)
	3423	DO i=2, klon
	3424	amn=MIN(ts(i,is_ter),amn)
	3425	amx=MAX(ts(i,is_ter),amx)
	3426	ENDDO
	3427	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
	3428	ENDIF
[1067]	3429
	3430	!jg ?
[996]	3431	!!$!
	3432	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
	3433	!!$! sub-surfaces is distributed.
	3434	!!$!
	3435	!!$ DO nsrf = 1, nbsrf
	3436	!!$ DO i = 1, klon
	3437	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
	3438	!!$ ts(i,nsrf) = zxtsol(i)
	3439	!!$ t2m(i,nsrf) = zt2m(i)
	3440	!!$ q2m(i,nsrf) = zq2m(i)
	3441	!!$ u10m(i,nsrf) = zu10m(i)
	3442	!!$ v10m(i,nsrf) = zv10m(i)
	3443	!!$
	3444	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
	3445	!!$ pblh(i,nsrf) = s_pblh(i)
	3446	!!$ plcl(i,nsrf) = s_plcl(i)
	3447	!!$ capCL(i,nsrf) = s_capCL(i)
	3448	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
	3449	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
	3450	!!$ pblT(i,nsrf) = s_pblT(i)
	3451	!!$ therm(i,nsrf) = s_therm(i)
	3452	!!$ trmb1(i,nsrf) = s_trmb1(i)
	3453	!!$ trmb2(i,nsrf) = s_trmb2(i)
	3454	!!$ trmb3(i,nsrf) = s_trmb3(i)
	3455	!!$ ENDIF
	3456	!!$ ENDDO
	3457	!!$ ENDDO
[781]	3458
	3459
	3460	DO i = 1, klon
	3461	fder(i) = - 4.0RSIGMAzxtsol(i)**3
	3462	ENDDO
	3463
	3464	zxqsurf(:) = 0.0
	3465	zxsnow(:) = 0.0
	3466	DO nsrf = 1, nbsrf
	3467	DO i = 1, klon
[4013]	3468	zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
[996]	3469	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
[3391]	3470	ENDDO
	3471	ENDDO
[781]	3472
[1067]	3473	! Premier niveau de vent sortie dans physiq.F
	3474	zu1(:) = u(:,1)
	3475	zv1(:) = v(:,1)
[781]	3476
	3477	END SUBROUTINE pbl_surface
	3478	!
	3479	!****************************************************************************************
	3480	!
[2243]	3481	SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
[781]	3482
[1785]	3483	USE indice_sol_mod
	3484
[781]	3485	INCLUDE "dimsoil.h"
	3486
	3487	! Ouput variables
	3488	!****************************************************************************************
	3489	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
	3490	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
	3491	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
	3492	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
	3493
	3494
	3495	!****************************************************************************************
	3496	! Return module variables for writing to restart file
	3497	!
	3498	!****************************************************************************************
	3499	fder_rst(:) = fder(:)
	3500	snow_rst(:,:) = snow(:,:)
	3501	qsurf_rst(:,:) = qsurf(:,:)
	3502	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
	3503
	3504	!****************************************************************************************
	3505	! Deallocate module variables
	3506	!
	3507	!****************************************************************************************
[1413]	3508	! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
	3509	IF (ALLOCATED(fder)) DEALLOCATE(fder)
	3510	IF (ALLOCATED(snow)) DEALLOCATE(snow)
	3511	IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
	3512	IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
[4013]	3513	IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
	3514	IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
[781]	3515
[3179]	3516	!jyg<
	3517	!****************************************************************************************
	3518	! Deallocate variables for pbl splitting
	3519	!
	3520	!****************************************************************************************
	3521
	3522	CALL wx_pbl_final
	3523	!>jyg
	3524
[781]	3525	END SUBROUTINE pbl_surface_final
	3526	!
	3527	!****************************************************************************************
[996]	3528	!
	3529
[2227]	3530	!albedo SB >>>
[3429]	3531	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
	3532	evap, z0m, z0h, agesno, &
	3533	tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)
	3534	!albedo SB <<<
[996]	3535	! Give default values where new fraction has appread
	3536
[1785]	3537	USE indice_sol_mod
[4020]	3538	use phys_state_var_mod, only: delta_sal, ds_ns, dt_ns, delta_sst, dter, &
	3539	dser, dt_ds
[3458]	3540	use config_ocean_skin_m, only: activate_ocean_skin
[1785]	3541
[996]	3542	INCLUDE "dimsoil.h"
	3543	INCLUDE "clesphys.h"
[1236]	3544	INCLUDE "compbl.h"
[996]	3545
	3546	! Input variables
	3547	!****************************************************************************************
	3548	INTEGER, INTENT(IN) :: itime
	3549	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
	3550
	3551	! InOutput variables
	3552	!****************************************************************************************
	3553	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
[2227]	3554	!albedo SB >>>
	3555	REAL, DIMENSION(klon,nsw,nbsrf), INTENT(INOUT) :: alb_dir, alb_dif
	3556	INTEGER :: k
	3557	!albedo SB <<<
[1670]	3558	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m
[2243]	3559	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: evap, agesno
	3560	REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT) :: z0m,z0h
[2181]	3561	REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
[996]	3562
	3563	! Local variables
	3564	!****************************************************************************************
	3565	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
	3566	CHARACTER(len=80) :: abort_message
	3567	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
	3568	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
	3569	!
	3570	! All at once !!
	3571	!****************************************************************************************
	3572
	3573	DO nsrf = 1, nbsrf
	3574	! First decide complement sub-surfaces
	3575	SELECT CASE (nsrf)
	3576	CASE(is_oce)
	3577	nsrf_comp1=is_sic
	3578	nsrf_comp2=is_ter
	3579	nsrf_comp3=is_lic
	3580	CASE(is_sic)
	3581	nsrf_comp1=is_oce
	3582	nsrf_comp2=is_ter
	3583	nsrf_comp3=is_lic
	3584	CASE(is_ter)
	3585	nsrf_comp1=is_lic
	3586	nsrf_comp2=is_oce
	3587	nsrf_comp3=is_sic
	3588	CASE(is_lic)
	3589	nsrf_comp1=is_ter
	3590	nsrf_comp2=is_oce
	3591	nsrf_comp3=is_sic
	3592	END SELECT
	3593
	3594	! Initialize all new fractions
	3595	DO i=1, klon
	3596	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
[1067]	3597
[996]	3598	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
	3599	! Use the complement sub-surface, keeping the continents unchanged
	3600	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
	3601	evap(i,nsrf) = evap(i,nsrf_comp1)
[2243]	3602	z0m(i,nsrf) = z0m(i,nsrf_comp1)
	3603	z0h(i,nsrf) = z0h(i,nsrf_comp1)
[996]	3604	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
[2227]	3605	!albedo SB >>>
	3606	DO k=1,nsw
	3607	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp1)
	3608	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp1)
	3609	ENDDO
	3610	!albedo SB <<<
[1670]	3611	ustar(i,nsrf) = ustar(i,nsrf_comp1)
[996]	3612	u10m(i,nsrf) = u10m(i,nsrf_comp1)
	3613	v10m(i,nsrf) = v10m(i,nsrf_comp1)
[3391]	3614	IF (iflag_pbl > 1) THEN
[1236]	3615	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
[3391]	3616	ENDIF
[996]	3617	mfois(nsrf) = mfois(nsrf) + 1
[2655]	3618	! F. Codron sensible default values for ocean and sea ice
	3619	IF (nsrf.EQ.is_oce) THEN
[3402]	3620	tsurf(i,nsrf) = 271.35
	3621	! (temperature of sea water under sea ice, so that
	3622	! is also the temperature of appearing sea water)
	3623	DO k=1,nsw
[2655]	3624	alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
	3625	alb_dif(i,k,nsrf) = 0.06
[3402]	3626	ENDDO
[3458]	3627	if (activate_ocean_skin >= 1) then
[3628]	3628	if (activate_ocean_skin == 2 &
[3740]	3629	.and. type_ocean == "couple") then
[3767]	3630	delta_sal(i) = 0.
[3744]	3631	delta_sst(i) = 0.
[4020]	3632	dter(i) = 0.
	3633	dser(i) = 0.
	3634	dt_ds(i) = 0.
[3740]	3635	end if
	3636
[3458]	3637	ds_ns(i) = 0.
	3638	dt_ns(i) = 0.
	3639	end if
[2655]	3640	ELSE IF (nsrf.EQ.is_sic) THEN
[3402]	3641	tsurf(i,nsrf) = 271.35
	3642	! (Temperature at base of sea ice. Surface
	3643	! temperature could be higher, up to 0 Celsius
	3644	! degrees. We set it to -1.8 Celsius degrees for
	3645	! consistency with the ocean slab model.)
	3646	DO k=1,nsw
[2655]	3647	alb_dir(i,k,nsrf) = 0.3 ! thin ice
	3648	alb_dif(i,k,nsrf) = 0.3
[3402]	3649	ENDDO
[3391]	3650	ENDIF
[996]	3651	ELSE
	3652	! The continents have changed. The new fraction receives the mean sum of the existent fractions
	3653	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
	3654	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
[2243]	3655	z0m(i,nsrf) = z0m(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0m(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
	3656	z0h(i,nsrf) = z0h(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0h(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
[996]	3657	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
[2227]	3658	!albedo SB >>>
	3659	DO k=1,nsw
	3660	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
	3661	alb_dir(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
	3662	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
	3663	alb_dif(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
	3664	ENDDO
	3665	!albedo SB <<<
[1670]	3666	ustar(i,nsrf) = ustar(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
[996]	3667	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
	3668	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
[3391]	3669	IF (iflag_pbl > 1) THEN
[1236]	3670	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
[3391]	3671	ENDIF
[996]	3672
	3673	! Security abort. This option has never been tested. To test, comment the following line.
	3674	! abort_message='The fraction of the continents have changed!'
[2311]	3675	! CALL abort_physic(modname,abort_message,1)
[996]	3676	nfois(nsrf) = nfois(nsrf) + 1
[3391]	3677	ENDIF
[996]	3678	snow(i,nsrf) = 0.
	3679	agesno(i,nsrf) = 0.
	3680	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
	3681	ELSE
	3682	pfois(nsrf) = pfois(nsrf)+ 1
[3391]	3683	ENDIF
	3684	ENDDO
[996]	3685
[3391]	3686	ENDDO
[996]	3687
	3688	END SUBROUTINE pbl_surface_newfrac
[781]	3689	!
[996]	3690	!****************************************************************************************
	3691	!
[781]	3692	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format