Context Navigation

pbl_surface_mod.F90 @ 4755

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

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

Download in other formats:

Original Format