Context Navigation

pbl_surface_mod.F90 @ 4500

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

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

Download in other formats:

Original Format