Context Navigation

pbl_surface_mod.F90 @ 5669

Last change on this file since 5669 was 5669, checked in by Laurent Fairhead, 3 weeks ago
Initialisation problem for 1D
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: 182.9 KB

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

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

Download in other formats:

Original Format