Context Navigation

pbl_surface_mod.F90 @ 5006

Last change on this file since 5006 was 4916, checked in by evignon, 9 months ago
correction formulation de l'erosion de la neige soufflee suite aux investigations de Nicolas Chiabrando, prise en compte d'un temps d'erosion de la neige fraiche nouvellement accumulee
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: 154.0 KB

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

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

Download in other formats:

Original Format