Context Navigation

pbl_surface_mod.F90 @ 4446

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

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

Download in other formats:

Original Format