Context Navigation

pbl_surface_mod.F90 @ 5220

Last change on this file since 5220 was 5179, checked in by abarral, 2 months ago
Removed always-on and <2007 T2m CPP key
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: 172.8 KB

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

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

Download in other formats:

Original Format