Context Navigation

pbl_surface_mod.F90 @ 2791

Last change on this file since 2791 was 2720, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2664:2719 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.3 KB

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

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

Download in other formats:

Original Format