Context Navigation

pbl_surface_mod.F90 @ 3761

Last change on this file since 3761 was 3756, checked in by oboucher, 4 years ago
The fraction of SW down radiation that is diffuse is extracted from radlwsw.F90 and passed on to pbl_surface for inclusion in fields_out for transfer to ORCHIDEE if it requested in input file coupling_fields.def. Hence there is nothing automatic here. The fraction of SW down radiation is also added as a diagnostic. To satisfy the case when this new quantity is used in ORCHIDEE, then it is added in the restart (and read in case it is there).
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: 132.4 KB

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

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

Download in other formats:

Original Format