Context Navigation

pbl_surface_mod.F90 @ 2435

Last change on this file since 2435 was 2435, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes r2396:2434 into testing branch
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: 125.0 KB

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

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

Download in other formats:

Original Format