Context Navigation

pbl_surface_mod.F90 @ 5028

Last change on this file since 5028 was 4368, checked in by lguez, 2 years ago
Sync latest trunk changes to Ocean_skin
File size: 167.6 KB

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

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

Download in other formats:

Original Format