Context Navigation

pbl_surface_mod.F90 @ 4476

Last change on this file since 4476 was 4449, checked in by evignon, 17 months ago
commission du nouveau schema de turbulence developpe dans le cadre de l'atelier tke
File size: 168.7 KB

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

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

Download in other formats:

Original Format