Context Navigation

pbl_surface_mod.F90 @ 5445

Last change on this file since 5445 was 3838, checked in by musat, 4 years ago
Ajout ustar pour diagnostics pbl
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: 134.4 KB

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

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

Download in other formats:

Original Format