Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

pbl_surface_mod.F90

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

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

Download in other formats:

Original Format