Context Navigation

pbl_surface_mod.F90 @ 3603

Last change on this file since 3603 was 3406, checked in by jghattas, 6 years ago
Added all modifications in the model code that were used for the simulations with DYANMICO during the Grand Challeng 2018. Modifications done by Y. Meurdesoif, L. Fairhead and A.K. Traore
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 129.8 KB

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

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

Download in other formats:

Original Format