Context Navigation

pbl_surface_mod.F90 @ 5017

Last change on this file since 5017 was 4743, checked in by Laurent Fairhead, 10 months ago
Merge of ACC branch with 4740 revision from trunk
File size: 172.2 KB

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

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

Download in other formats:

Original Format