Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

pbl_surface_mod.F90

Last change on this file was 4727, checked in by idelkadi, 9 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 153.5 KB

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

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

Download in other formats:

Original Format