Context Navigation

pbl_surface_mod.F90 @ 5212

Last change on this file since 5212 was 2886, checked in by Laurent Fairhead, 7 years ago
Merged trunk changes r2865:2885 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: 127.4 KB

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

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

Download in other formats:

Original Format