Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

pbl_surface_mod.F90

Last change on this file was 5310, checked in by abarral, 7 weeks ago
unify abort_gcm rename wxios -> wxios_mod
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 174.3 KB

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

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

Download in other formats:

Original Format