Context Navigation

pbl_surface_mod.F90 @ 4722

Last change on this file since 4722 was 4687, checked in by evignon, 16 months ago
je renomme les routines atke suivant la convention decidee pour la reecriture
File size: 172.5 KB

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

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

Download in other formats:

Original Format