Context Navigation

pbl_surface_mod.F90 @ 5456

Last change on this file since 5456 was 5310, checked in by abarral, 8 weeks ago
unify abort_gcm rename wxios -> wxios_mod
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 174.3 KB

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

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

Download in other formats:

Original Format