Context Navigation

pbl_surface_mod.F90 @ 5495

Last change on this file since 5495 was 5486, checked in by evignon, 3 weeks ago
inclusion d'un diagnostique de la sublimation de la glace sur les landice pour la conservation de l'eau
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 174.5 KB

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

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

Download in other formats:

Original Format