Context Navigation

pbl_surface_mod.F90 @ 4009

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

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

Download in other formats:

Original Format