Context Navigation

pbl_surface_mod.F90 @ 5761

Last change on this file since 5761 was 5717, checked in by aborella, 4 weeks ago
Merge with trunk r5653
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 178.6 KB

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

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

Download in other formats:

Original Format