Context Navigation

pbl_surface_mod.F90 @ 4799

Last change on this file since 4799 was 2952, checked in by Laurent Fairhead, 7 years ago

Parametrization of drag by copses
Need version 4465 of ORCHIDEE at least

Cheruy

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:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 129.8 KB

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

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

Download in other formats:

Original Format