Context Navigation

pbl_surface_mod.F90 @ 3780

Last change on this file since 3780 was 3780, checked in by evignon, 4 years ago
Premiere comission Etienne: changements pour le 1D (forcage en Ts au dessus des continents) et inclusion drag arbres dans yamada4_num=6
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:keywords set to `Author Date Id Revision`
File size: 132.7 KB

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

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

Download in other formats:

Original Format