Context Navigation

pbl_surface_mod.F90 @ 4774

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

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

Download in other formats:

Original Format