Context Navigation

pbl_surface_mod.F90 @ 5761

Last change on this file since 5761 was 5717, checked in by aborella, 4 weeks ago
Merge with trunk r5653
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: 178.6 KB

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

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

Download in other formats:

Original Format