Context Navigation

pbl_surface_mod.F90 @ 5128

Last change on this file since 5128 was 5128, checked in by abarral, 4 months ago
Correct bug in vlspltqs_loc.f90 from r2270 where we call SSUM with incorrect arguments. Merge the three different versions of abort_gcm into one Fix seq, para 3D compilation broken from r5107 onwards (lint) usual + Remove uneeded fixed-form continuations
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.4 KB

Line
1	! $Id: pbl_surface_mod.F90 5128 2024-07-25 15:47:25Z 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
417	IMPLICIT NONE
418
419	INCLUDE "dimsoil.h"
420	INCLUDE "YOMCST.h"
421	INCLUDE "YOETHF.h"
422	INCLUDE "FCTTRE.h"
423	INCLUDE "clesphys.h"
424	INCLUDE "compbl.h"
425	INCLUDE "flux_arp.h"
426	!FC
427	INCLUDE "dimpft.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	ENDDO
2704
2705	DO j = 1, knon
2706	y_d_ts(j) = ytsurf_new(j) - yts(j)
2707	ENDDO
2708
2709	ELSE ! (ok_flux_surf)
2710	DO j = 1, knon
2711	y_flux_t1(j) = yfluxsens(j)
2712	y_flux_q1(j) = -yevap(j)
2713	#ifdef ISO
2714	y_flux_xt1(:,:) = -yxtevap(:,:)
2715	#endif
2716	ENDDO
2717	ENDIF ! (ok_flux_surf)
2718
2719	! flux of blowing snow at the first level
2720	IF (ok_bs) THEN
2721	DO j = 1, knon
2722	y_flux_bs(j) = yfluxbs(j)
2723	ENDDO
2724	ENDIF
2725
2726	! ------------------------------------------------------------------------------
2727	! 12a) Splitting
2728	! ------------------------------------------------------------------------------
2729
2730	IF (iflag_split >= 1) THEN
2731	#ifdef ISO
2732	CALL abort_gcm('pbl_surface_mod 2607','isos pas encore dans iflag_split=1',1)
2733	#endif
2734
2735
2736	IF (nsrf /= is_oce) THEN
2737
2738	! Compute potential evaporation and aridity factor (jyg, 20200328)
2739	ybeta_prev(:) = ybeta(:)
2740	DO j = 1, knon
2741	yqa(j) = AcoefQ(j) - BcoefQ(j) * yevap(j) * dtime
2742	ENDDO
2743
2744	CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot)
2745
2746	ybeta(1:knon) = min(yevap(1:knon) / yevap_pot(1:knon), 1.)
2747
2748	IF (prt_level >=10) THEN
2749	DO j = 1, knon
2750	PRINT*, 'y_flux_t1,yfluxlat,wakes' &
2751	, y_flux_t1(j), yfluxlat(j), ywake_s(j)
2752	PRINT*, 'beta_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j)
2753	PRINT*, 'inertia,facteur,cstar', inertia, facteur, wake_cstar(j)
2754	ENDDO
2755	ENDIF ! (prt_level >=10)
2756
2757	! Second CALL to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account
2758	! the update of the aridity coeficient beta.
2759
2760	CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, &
2761	BcoefQ_x, BcoefQ_w &
2762	)
2763	CALL wx_pbl0_merge(knon, ypplay, ypaprs, &
2764	ywake_s, ydTs0, ydqs0, &
2765	yt_x, yt_w, yq_x, yq_w, &
2766	yu_x, yu_w, yv_x, yv_w, &
2767	ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
2768	ycdragm_x, ycdragm_w, &
2769	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
2770	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
2771	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
2772	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
2773	AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
2774	BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
2775	ycdragh, ycdragq, ycdragm, &
2776	yt1, yq1, yu1, yv1 &
2777	)
2778	IF (iflag_split == 2) THEN
2779	CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
2780	ywake_s, ybeta, ywake_cstar, ywake_dens, &
2781	AcoefH_x, AcoefH_w, &
2782	BcoefH_x, BcoefH_w, &
2783	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2784	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2785	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2786	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2787	yg_T, yg_Q, &
2788	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2789	ydTs_ins, ydqs_ins &
2790	)
2791	ELSE !
2792	AcoefH(:) = AcoefH_0(:)
2793	AcoefQ(:) = AcoefQ_0(:)
2794	BcoefH(:) = BcoefH_0(:)
2795	BcoefQ(:) = BcoefQ_0(:)
2796	yg_T(:) = 0.
2797	yg_Q(:) = 0.
2798	yGamma_dTs_phiT(:) = 0.
2799	yGamma_dQs_phiQ(:) = 0.
2800	ydTs_ins(:) = 0.
2801	ydqs_ins(:) = 0.
2802	ENDIF ! (iflag_split .EQ. 2)
2803
2804	ELSE ! (nsrf .NE. is_oce)
2805	ybeta(1:knon) = 1.
2806	yevap_pot(1:knon) = yevap(1:knon)
2807	AcoefH(:) = AcoefH_0(:)
2808	AcoefQ(:) = AcoefQ_0(:)
2809	BcoefH(:) = BcoefH_0(:)
2810	BcoefQ(:) = BcoefQ_0(:)
2811	yg_T(:) = 0.
2812	yg_Q(:) = 0.
2813	yGamma_dTs_phiT(:) = 0.
2814	yGamma_dQs_phiQ(:) = 0.
2815	ydTs_ins(:) = 0.
2816	ydqs_ins(:) = 0.
2817	ENDIF ! (nsrf .NE. is_oce)
2818
2819	CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, &
2820	yg_T, yg_Q, &
2821	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2822	ydTs_ins, ydqs_ins, &
2823	y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
2824	!!!! HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
2825	phiQ0_b, phiT0_b, &
2826	y_flux_t1_x, y_flux_t1_w, &
2827	y_flux_q1_x, y_flux_q1_w, &
2828	y_flux_u1_x, y_flux_u1_w, &
2829	y_flux_v1_x, y_flux_v1_w, &
2830	yfluxlat_x, yfluxlat_w, &
2831	y_delta_qsats, &
2832	y_delta_tsurf_new, y_delta_qsurf &
2833	)
2834
2835	CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
2836	yTs, y_delta_tsurf, &
2837	yqsurf, yTsurf_new, &
2838	y_delta_tsurf_new, y_delta_qsats, &
2839	AcoefH_x, AcoefH_w, &
2840	BcoefH_x, BcoefH_w, &
2841	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2842	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2843	y_flux_t1, y_flux_q1, &
2844	y_flux_t1_x, y_flux_t1_w, &
2845	y_flux_q1_x, y_flux_q1_w)
2846
2847	IF (nsrf /= is_oce) THEN
2848	CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
2849	yTs, y_delta_tsurf, &
2850	yqsurf, yTsurf_new, &
2851	y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf, &
2852	AcoefH_x, AcoefH_w, &
2853	BcoefH_x, BcoefH_w, &
2854	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2855	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2856	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2857	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2858	yg_T, yg_Q, &
2859	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2860	ydTs_ins, ydqs_ins, &
2861	y_flux_t1, y_flux_q1, &
2862	y_flux_t1_x, y_flux_t1_w, &
2863	y_flux_q1_x, y_flux_q1_w)
2864	ENDIF ! (nsrf .NE. is_oce)
2865
2866	ELSE ! (iflag_split .ge. 1)
2867	ybeta(1:knon) = 1.
2868	yevap_pot(1:knon) = yevap(1:knon)
2869	ENDIF ! (iflag_split .ge. 1)
2870
2871	IF (prt_level >= 10) THEN
2872	print *, 'pbl_surface, ybeta , yevap, yevap_pot ', &
2873	ybeta(1:knon), yevap(1:knon), yevap_pot(1:knon)
2874	ENDIF ! (prt_level >= 10)
2875
2876	!>jyg
2877
2878	!!jyg!! A reprendre apres reflexion ===============================================
2879	!!jyg!!
2880	!!jyg!! DO j=1,knon
2881	!!jyg!!!!! nrlmd le 13/06/2011
2882	!!jyg!!
2883	!!jyg!!!----Diffusion dans le sol dans le cas continental seulement
2884	!!jyg!! IF (nsrf.EQ.is_ter) THEN
2885	!!jyg!!!----Calcul du coefficient delta_coeff
2886	!!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)))
2887	!!jyg!!
2888	!!jyg!!! delta_coef(j)=dtime/(inertia*sqrt(tau_eq(j)))
2889	!!jyg!! delta_coef(j)=facteur*sqrt(tau_eq(j))/inertia
2890	!!jyg!!! delta_coef(j)=0.
2891	!!jyg!! ELSE
2892	!!jyg!! delta_coef(j)=0.
2893	!!jyg!! ENDIF
2894	!!jyg!!
2895	!!jyg!!!----Calcul de delta_tsurf
2896	!!jyg!! y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j)
2897	!!jyg!!
2898	!!jyg!!!----Si il n'y a pas des poches...
2899	!!jyg!! IF (wake_cstar(j).le.0.01) THEN
2900	!!jyg!! y_delta_tsurf(j)=0.
2901	!!jyg!! y_delta_flux_t1(j)=0.
2902	!!jyg!! ENDIF
2903	!!jyg!!
2904	!!jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle)
2905	!!jyg!!!!!!! jyg le 23/02/2012
2906	!!jyg!!!!!!!
2907	!!jyg!!!! ybeta(j)=y_flux_q1(j) / &
2908	!!jyg!!!! & (Kech_h(j)*(yq(j,1)-yqsatsurf(j)))
2909	!!jyg!!!!!! ybeta(j)=-1.*yevap(j) / &
2910	!!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)))
2911	!!jyg!!!!!!! fin jyg
2912	!!jyg!!!!!
2913	!!jyg!!
2914	!!jyg!! ENDDO
2915	!!jyg!!
2916	!!jyg!!!!! fin nrlmd le 13/06/2011
2917	!!jyg!!
2918	IF (iflag_split >= 1) THEN
2919	IF (prt_level >=10) THEN
2920	DO j = 1, knon
2921	PRINT*, 'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
2922	PRINT*, 'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
2923	PRINT*, 't1x, t1w, t1, t1_ancien', &
2924	yt_x(j, 1), yt_w(j, 1), yt(j, 1), t(j, 1)
2925	PRINT*, 'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j)
2926	ENDDO
2927
2928	DO j = 1, knon
2929	PRINT*, 'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
2930	, 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)
2931	PRINT*, 'beta, ytsurf_new ', ybeta(j), ytsurf_new(j)
2932	PRINT*, 'inertia, facteur, cstar', inertia, facteur, wake_cstar(j)
2933	ENDDO
2934	ENDIF ! (prt_level >=10)
2935
2936	!!! jyg le 07/02/2012
2937	ENDIF ! (iflag_split .ge.1)
2938	!!!
2939
2940	!!! jyg le 07/02/2012
2941	IF (iflag_split ==0) THEN
2942	!!!
2943	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2944	CALL climb_hq_up(knon, dtime, yt, yq, &
2945	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
2946	!!! jyg le 07/02/2012
2947	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2948	CcoefH, CcoefQ, DcoefH, DcoefQ, &
2949	Kcoef_hq, gama_q, gama_h, &
2950	!!!
2951	y_flux_q(:, :), y_flux_t(:, :), y_d_q(:, :), y_d_t(:, :) &
2952	#ifdef ISO
2953	,yxt,y_flux_xt1 &
2954	,AcoefXT,BcoefXT,CcoefXT,DcoefXT,gama_xt &
2955	,y_flux_xt(:,:,:),y_d_xt(:,:,:) &
2956	#endif
2957	)
2958	ELSE !(iflag_split .EQ.0)
2959	CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
2960	y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &
2961	!!! nrlmd le 02/05/2011
2962	AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &
2963	CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
2964	Kcoef_hq_x, gama_q_x, gama_h_x, &
2965	!!!
2966	y_flux_q_x(:, :), y_flux_t_x(:, :), y_d_q_x(:, :), y_d_t_x(:, :) &
2967	#ifdef ISO
2968	,yxt_x,y_flux_xt1_x &
2969	,AcoefXT_x,BcoefXT_x,CcoefXT_x,DcoefXT_x,gama_xt_x &
2970	,y_flux_xt_x(:,:,:),y_d_xt_x(:,:,:) &
2971	#endif
2972	)
2973
2974	CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
2975	y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &
2976	!!! nrlmd le 02/05/2011
2977	AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &
2978	CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
2979	Kcoef_hq_w, gama_q_w, gama_h_w, &
2980	!!!
2981	y_flux_q_w(:, :), y_flux_t_w(:, :), y_d_q_w(:, :), y_d_t_w(:, :) &
2982	#ifdef ISO
2983	,yxt_w,y_flux_xt1_w &
2984	,AcoefXT_w,BcoefXT_w,CcoefXT_w,DcoefXT_w,gama_xt_w &
2985	,y_flux_xt_w(:,:,:),y_d_xt_w(:,:,:) &
2986	#endif
2987	)
2988	!!!
2989	ENDIF ! (iflag_split .EQ.0)
2990	!!!
2991
2992	!!! jyg le 07/02/2012
2993	IF (iflag_split ==0) THEN
2994	!!!
2995	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2996	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
2997	!!! jyg le 07/02/2012
2998	AcoefU, AcoefV, BcoefU, BcoefV, &
2999	CcoefU, CcoefV, DcoefU, DcoefV, &
3000	Kcoef_m, &
3001	!!!
3002	y_flux_u, y_flux_v, y_d_u, y_d_v)
3003	y_d_t_diss(:, :) = 0.
3004	IF (iflag_pbl>=20 .AND. iflag_pbl<30) THEN
3005	CALL yamada_c(knon, dtime, ypaprs, ypplay &
3006	, yu, yv, yt, y_d_u, y_d_v, y_d_t, ycdragm, ytke, ycoefm, ycoefh, ycoefq, y_d_t_diss, yustar &
3007	, iflag_pbl)
3008	ENDIF
3009	! PRINT*,'yamada_c OK'
3010
3011	ELSE !(iflag_split .EQ.0)
3012	CALL climb_wind_up(knon, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &
3013	!!! nrlmd le 02/05/2011
3014	AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &
3015	CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
3016	Kcoef_m_x, &
3017	!!!
3018	y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)
3019
3020	y_d_t_diss_x(:, :) = 0.
3021	IF (iflag_pbl>=20 .AND. iflag_pbl<30) THEN
3022	CALL yamada_c(knon, dtime, ypaprs, ypplay &
3023	, 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 &
3024	, ycoefq_x, y_d_t_diss_x, yustar_x &
3025	, iflag_pbl)
3026	ENDIF
3027	! PRINT*,'yamada_c OK'
3028
3029	CALL climb_wind_up(knon, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &
3030	!!! nrlmd le 02/05/2011
3031	AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &
3032	CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
3033	Kcoef_m_w, &
3034	!!!
3035	y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)
3036	!!!
3037	y_d_t_diss_w(:, :) = 0.
3038	IF (iflag_pbl>=20 .AND. iflag_pbl<30) THEN
3039	CALL yamada_c(knon, dtime, ypaprs, ypplay &
3040	, 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 &
3041	, ycoefq_w, y_d_t_diss_w, yustar_w &
3042	, iflag_pbl)
3043	ENDIF
3044	! PRINT*,'yamada_c OK'
3045
3046	IF (prt_level >=10) THEN
3047	print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &
3048	yfluxlat_x(1:knon), yfluxlat_w(1:knon)
3049	ENDIF
3050
3051	ENDIF ! (iflag_split .EQ.0)
3052
3053	IF (ok_bs) THEN
3054	CALL climb_qbs_up(knon, dtime, yqbs, &
3055	y_flux_bs, ypaprs, ypplay, &
3056	AcoefQBS, BcoefQBS, &
3057	CcoefQBS, DcoefQBS, &
3058	Kcoef_qbs, gama_qbs, &
3059	y_flux_qbs(:, :), y_d_qbs(:, :))
3060	ENDIF
3061
3062	!!!
3063	!!
3064	!! DO j = 1, knon
3065	!! y_dflux_t(j) = y_dflux_t(j) * ypct(j)
3066	!! y_dflux_q(j) = y_dflux_q(j) * ypct(j)
3067	!! ENDDO
3068	!!
3069	!****************************************************************************************
3070	! 13) Transform variables for output format :
3071	! - Decompress
3072	! - Multiply with pourcentage of current surface
3073	! - Cumulate in global variable
3074
3075	!****************************************************************************************
3076
3077
3078	!!! jyg le 07/02/2012
3079	IF (iflag_split==0) THEN
3080	!!!
3081	DO k = 1, klev
3082	DO j = 1, knon
3083	i = ni(j)
3084	y_d_t_diss(j, k) = y_d_t_diss(j, k) * ypct(j)
3085	y_d_t(j, k) = y_d_t(j, k) * ypct(j)
3086	y_d_q(j, k) = y_d_q(j, k) * ypct(j)
3087	y_d_u(j, k) = y_d_u(j, k) * ypct(j)
3088	y_d_v(j, k) = y_d_v(j, k) * ypct(j)
3089	!FC
3090	IF (nsrf == is_ter .AND. ifl_pbltree >= 1) THEN
3091	! if (y_d_u_frein(j,k).NE.0. ) THEN
3092	! PRINT, nsrf,'IS_TER ++', y_d_u_frein(j,k)ypct(j),y_d_u(j,k),j,k
3093	! ENDIF
3094	y_d_u(j, k) = y_d_u(j, k) + y_d_u_frein(j, k) * ypct(j)
3095	y_d_v(j, k) = y_d_v(j, k) + y_d_v_frein(j, k) * ypct(j)
3096	treedrg(i, k, nsrf) = y_treedrg(j, k)
3097	ELSE
3098	treedrg(i, k, nsrf) = 0.
3099	ENDIF
3100	!FC
3101	flux_t(i, k, nsrf) = y_flux_t(j, k)
3102	flux_q(i, k, nsrf) = y_flux_q(j, k)
3103	flux_u(i, k, nsrf) = y_flux_u(j, k)
3104	flux_v(i, k, nsrf) = y_flux_v(j, k)
3105
3106	#ifdef ISO
3107	DO ixt=1,ntraciso
3108	y_d_xt(ixt,j,k) = y_d_xt(ixt,j,k) * ypct(j)
3109	flux_xt(ixt,i,k,nsrf) = y_flux_xt(ixt,j,k)
3110	ENDDO ! DO ixt=1,ntraciso
3111	h1_diag(i)=h1(j)
3112	#endif
3113
3114	ENDDO
3115	ENDDO
3116
3117	#ifdef ISO
3118	#ifdef ISOVERIF
3119	IF (iso_eau.gt.0) THEN
3120	CALL iso_verif_egalite_vect2D( &
3121	y_d_xt,y_d_q, &
3122	'pbl_surface_mod 2600',ntraciso,klon,klev)
3123	endif
3124	#endif
3125	#endif
3126
3127	ELSE !(iflag_split .EQ.0)
3128
3129	! Tendances hors poches
3130	DO k = 1, klev
3131	DO j = 1, knon
3132	i = ni(j)
3133	y_d_t_diss_x(j, k) = y_d_t_diss_x(j, k) * ypct(j)
3134	y_d_t_x(j, k) = y_d_t_x(j, k) * ypct(j)
3135	y_d_q_x(j, k) = y_d_q_x(j, k) * ypct(j)
3136	y_d_u_x(j, k) = y_d_u_x(j, k) * ypct(j)
3137	y_d_v_x(j, k) = y_d_v_x(j, k) * ypct(j)
3138
3139	flux_t_x(i, k, nsrf) = y_flux_t_x(j, k)
3140	flux_q_x(i, k, nsrf) = y_flux_q_x(j, k)
3141	flux_u_x(i, k, nsrf) = y_flux_u_x(j, k)
3142	flux_v_x(i, k, nsrf) = y_flux_v_x(j, k)
3143
3144	#ifdef ISO
3145	DO ixt=1,ntraciso
3146	y_d_xt_x(ixt,j,k) = y_d_xt_x(ixt,j,k) * ypct(j)
3147	flux_xt_x(ixt,i,k,nsrf) = y_flux_xt_x(ixt,j,k)
3148	ENDDO ! DO ixt=1,ntraciso
3149	#endif
3150	ENDDO
3151	ENDDO
3152
3153	! Tendances dans les poches
3154	DO k = 1, klev
3155	DO j = 1, knon
3156	i = ni(j)
3157	y_d_t_diss_w(j, k) = y_d_t_diss_w(j, k) * ypct(j)
3158	y_d_t_w(j, k) = y_d_t_w(j, k) * ypct(j)
3159	y_d_q_w(j, k) = y_d_q_w(j, k) * ypct(j)
3160	y_d_u_w(j, k) = y_d_u_w(j, k) * ypct(j)
3161	y_d_v_w(j, k) = y_d_v_w(j, k) * ypct(j)
3162
3163	flux_t_w(i, k, nsrf) = y_flux_t_w(j, k)
3164	flux_q_w(i, k, nsrf) = y_flux_q_w(j, k)
3165	flux_u_w(i, k, nsrf) = y_flux_u_w(j, k)
3166	flux_v_w(i, k, nsrf) = y_flux_v_w(j, k)
3167
3168	#ifdef ISO
3169	DO ixt=1,ntraciso
3170	y_d_xt_w(ixt,j,k) = y_d_xt_w(ixt,j,k) * ypct(j)
3171	flux_xt_w(ixt,i,k,nsrf) = y_flux_xt_w(ixt,j,k)
3172	ENDDO ! do ixt=1,ntraciso
3173	#endif
3174
3175	ENDDO
3176	ENDDO
3177
3178	! Flux, tendances et Tke moyenne dans la maille
3179	DO k = 1, klev
3180	DO j = 1, knon
3181	i = ni(j)
3182	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))
3183	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))
3184	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))
3185	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))
3186	#ifdef ISO
3187	DO ixt=1,ntraciso
3188	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))
3189	ENDDO ! do ixt=1,ntraciso
3190	#endif
3191	ENDDO
3192	ENDDO
3193	DO j = 1, knon
3194	yfluxlat(j) = yfluxlat_x(j) + ywake_s(j) * (yfluxlat_w(j) - yfluxlat_x(j))
3195	ENDDO
3196	IF (prt_level >=10) THEN
3197	print *, ' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &
3198	nsrf, flux_t(:, 1, nsrf), flux_t_x(:, 1, nsrf), flux_t_w(:, 1, nsrf)
3199	ENDIF
3200
3201	DO k = 1, klev
3202	DO j = 1, knon
3203	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))
3204	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))
3205	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))
3206	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))
3207	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))
3208	ENDDO
3209	ENDDO
3210
3211	ENDIF ! (iflag_split .EQ.0)
3212	!!!
3213
3214	! tendencies of blowing snow
3215	IF (ok_bs) THEN
3216	DO k = 1, klev
3217	DO j = 1, knon
3218	i = ni(j)
3219	y_d_qbs(j, k) = y_d_qbs(j, k) * ypct(j)
3220	flux_qbs(i, k, nsrf) = y_flux_qbs(j, k)
3221	ENDDO
3222	ENDDO
3223	ENDIF
3224
3225	DO j = 1, knon
3226	i = ni(j)
3227	evap(i, nsrf) = - flux_q(i, 1, nsrf) !jyg
3228	IF (ok_bs) then ; snowerosion(i, nsrf) = flux_qbs(i, 1, nsrf);
3229	endif
3230	beta(i, nsrf) = ybeta(j) !jyg
3231	d_ts(i, nsrf) = y_d_ts(j)
3232	!albedo SB >>>
3233	DO k = 1, nsw
3234	alb_dir(i, k, nsrf) = yalb_dir_new(j, k)
3235	alb_dif(i, k, nsrf) = yalb_dif_new(j, k)
3236	ENDDO
3237	!albedo SB <<<
3238	snow(i, nsrf) = ysnow(j)
3239	qsurf(i, nsrf) = yqsurf(j)
3240	z0m(i, nsrf) = yz0m(j)
3241	z0h(i, nsrf) = yz0h(j)
3242	fluxlat(i, nsrf) = yfluxlat(j)
3243	agesno(i, nsrf) = yagesno(j)
3244	cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
3245	cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
3246	dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypct(j)
3247	dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypct(j)
3248	#ifdef ISO
3249	DO ixt=1,niso
3250	xtsnow(ixt,i,nsrf) = yxtsnow(ixt,j)
3251	ENDDO
3252	DO ixt=1,ntraciso
3253	xtevap(ixt,i,nsrf) = - flux_xt(ixt,i,1,nsrf)
3254	dflux_xt(ixt,i) = dflux_xt(ixt,i) + y_dflux_xt(ixt,j)*ypct(j)
3255	ENDDO
3256	IF (nsrf == is_lic) THEN
3257	DO ixt=1,niso
3258	Rland_ice(ixt,i) = yRland_ice(ixt,j)
3259	ENDDO
3260	ENDIF !IF (nsrf == is_lic) THEN
3261	#ifdef ISOVERIF
3262	IF (iso_eau.gt.0) THEN
3263	CALL iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
3264	'pbl_surf_mod 1230',errmax,errmaxrel)
3265	ENDIF !if (iso_eau.gt.0) THEN
3266	#endif
3267	#endif
3268	ENDDO
3269
3270	! PRINT*,'Dans pbl OK2'
3271
3272	!!! jyg le 07/02/2012
3273	IF (iflag_split >=1) THEN
3274	!!!
3275	!!! nrlmd le 02/05/2011
3276	DO j = 1, knon
3277	i = ni(j)
3278	fluxlat_x(i, nsrf) = yfluxlat_x(j)
3279	fluxlat_w(i, nsrf) = yfluxlat_w(j)
3280	!!!
3281	!!! nrlmd le 13/06/2011
3282	!!jyg20170131 delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
3283	!!jyg20210118 delta_tsurf(i,nsrf)=y_delta_tsurf(j)
3284	delta_tsurf(i, nsrf) = y_delta_tsurf_new(j)
3285
3286	delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j) * ypct(j)
3287
3288	cdragh_x(i) = cdragh_x(i) + ycdragh_x(j) * ypct(j)
3289	cdragh_w(i) = cdragh_w(i) + ycdragh_w(j) * ypct(j)
3290	cdragm_x(i) = cdragm_x(i) + ycdragm_x(j) * ypct(j)
3291	cdragm_w(i) = cdragm_w(i) + ycdragm_w(j) * ypct(j)
3292	kh(i) = kh(i) + Kech_h(j) * ypct(j)
3293	kh_x(i) = kh_x(i) + Kech_h_x(j) * ypct(j)
3294	kh_w(i) = kh_w(i) + Kech_h_w(j) * ypct(j)
3295	!!!
3296	ENDDO
3297	!!!
3298	ENDIF ! (iflag_split .ge.1)
3299	!!!
3300	!!! nrlmd le 02/05/2011
3301	!!jyg le 20/02/2011
3302	!! tke_x(:,:,nsrf)=0.
3303	!! tke_w(:,:,nsrf)=0.
3304	!!jyg le 20/02/2011
3305	!! DO k = 1, klev+1
3306	!! DO j = 1, knon
3307	!! i = ni(j)
3308	!! wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
3309	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
3310	!! ENDDO
3311	!! ENDDO
3312	!!jyg le 20/02/2011
3313	!! DO k = 1, klev+1
3314	!! DO j = 1, knon
3315	!! i = ni(j)
3316	!! tke(i,k,nsrf)=(1.-ywake_s(j))tke_x(i,k,nsrf)+ywake_s(j)tke_w(i,k,nsrf)
3317	!! ENDDO
3318	!! ENDDO
3319	!!!
3320	IF (iflag_split ==0) THEN
3321	wake_dltke(:, :, nsrf) = 0.
3322	DO k = 1, klev + 1
3323	DO j = 1, knon
3324	i = ni(j)
3325	!jyg<
3326	!! tke(i,k,nsrf) = ytke(j,k)
3327	!! tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
3328	tke_x(i, k, nsrf) = ytke(j, k)
3329	tke_x(i, k, is_ave) = tke_x(i, k, is_ave) + ytke(j, k) * ypct(j)
3330	eps_x(i, k, nsrf) = yeps(j, k)
3331	eps_x(i, k, is_ave) = eps_x(i, k, is_ave) + yeps(j, k) * ypct(j)
3332	!>jyg
3333	ENDDO
3334	ENDDO
3335
3336	ELSE ! (iflag_split .EQ.0)
3337	DO k = 1, klev + 1
3338	DO j = 1, knon
3339	i = ni(j)
3340	wake_dltke(i, k, nsrf) = ytke_w(j, k) - ytke_x(j, k)
3341	!jyg<
3342	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
3343	!! tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
3344	tke_x(i, k, nsrf) = ytke_x(j, k)
3345	tke_x(i, k, is_ave) = tke_x(i, k, is_ave) + tke_x(i, k, nsrf) * ypct(j)
3346	eps_x(i, k, nsrf) = yeps_x(j, k)
3347	eps_x(i, k, is_ave) = eps_x(i, k, is_ave) + eps_x(i, k, nsrf) * ypct(j)
3348	wake_dltke(i, k, is_ave) = wake_dltke(i, k, is_ave) + wake_dltke(i, k, nsrf) * ypct(j)
3349
3350
3351	!>jyg
3352	ENDDO
3353	ENDDO
3354	ENDIF ! (iflag_split .EQ.0)
3355	!!!
3356	DO k = 2, klev
3357	DO j = 1, knon
3358	i = ni(j)
3359	zcoefh(i, k, nsrf) = ycoefh(j, k)
3360	zcoefm(i, k, nsrf) = ycoefm(j, k)
3361	zcoefh(i, k, is_ave) = zcoefh(i, k, is_ave) + ycoefh(j, k) * ypct(j)
3362	zcoefm(i, k, is_ave) = zcoefm(i, k, is_ave) + ycoefm(j, k) * ypct(j)
3363	ENDDO
3364	ENDDO
3365
3366	! PRINT*,'Dans pbl OK3'
3367
3368	IF (nsrf == is_ter) THEN
3369	DO j = 1, knon
3370	i = ni(j)
3371	qsol(i) = yqsol(j)
3372	#ifdef ISO
3373	runoff_diag(i)=yrunoff_diag(j)
3374	DO ixt=1,niso
3375	xtsol(ixt,i) = yxtsol(ixt,j)
3376	xtrunoff_diag(ixt,i)=yxtrunoff_diag(ixt,j)
3377	ENDDO
3378	#endif
3379	ENDDO
3380	ENDIF
3381
3382	!jyg<
3383	!! ftsoil(:,:,nsrf) = 0.
3384	!>jyg
3385	DO k = 1, nsoilmx
3386	DO j = 1, knon
3387	i = ni(j)
3388	ftsoil(i, k, nsrf) = ytsoil(j, k)
3389	ENDDO
3390	ENDDO
3391
3392	#ifdef ISO
3393	#ifdef ISOVERIF
3394	!WRITE(,) 'pbl_surface 2858'
3395	DO i = 1, klon
3396	DO ixt=1,niso
3397	CALL iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 1405')
3398	ENDDO
3399	ENDDO
3400	#endif
3401	#ifdef ISOVERIF
3402	IF (iso_eau.gt.0) THEN
3403	CALL iso_verif_egalite_vect2D( &
3404	y_d_xt,y_d_q, &
3405	'pbl_surface_mod 1261',ntraciso,klon,klev)
3406	ENDIF !if (iso_eau.gt.0) THEN
3407	#endif
3408	#endif
3409	!!! jyg le 07/02/2012
3410	IF (iflag_split >=1) THEN
3411	!!!
3412	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
3413	DO k = 1, klev
3414	DO j = 1, knon
3415	i = ni(j)
3416	d_t_diss_x(i, k) = d_t_diss_x(i, k) + y_d_t_diss_x(j, k)
3417	d_t_x(i, k) = d_t_x(i, k) + y_d_t_x(j, k)
3418	d_q_x(i, k) = d_q_x(i, k) + y_d_q_x(j, k)
3419	d_u_x(i, k) = d_u_x(i, k) + y_d_u_x(j, k)
3420	d_v_x(i, k) = d_v_x(i, k) + y_d_v_x(j, k)
3421
3422	d_t_diss_w(i, k) = d_t_diss_w(i, k) + y_d_t_diss_w(j, k)
3423	d_t_w(i, k) = d_t_w(i, k) + y_d_t_w(j, k)
3424	d_q_w(i, k) = d_q_w(i, k) + y_d_q_w(j, k)
3425	d_u_w(i, k) = d_u_w(i, k) + y_d_u_w(j, k)
3426	d_v_w(i, k) = d_v_w(i, k) + y_d_v_w(j, k)
3427	#ifdef ISO
3428	DO ixt=1,ntraciso
3429	d_xt_x(ixt,i,k) = d_xt_x(ixt,i,k) + y_d_xt_x(ixt,j,k)
3430	d_xt_w(ixt,i,k) = d_xt_w(ixt,i,k) + y_d_xt_w(ixt,j,k)
3431	ENDDO ! DO ixt=1,ntraciso
3432	#endif
3433
3434	!! d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
3435	!! d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)
3436	ENDDO
3437	ENDDO
3438	!!!
3439	ENDIF ! (iflag_split .ge.1)
3440	!!!
3441
3442	DO k = 1, klev
3443	DO j = 1, knon
3444	i = ni(j)
3445	d_t_diss(i, k) = d_t_diss(i, k) + y_d_t_diss(j, k)
3446	d_t(i, k) = d_t(i, k) + y_d_t(j, k)
3447	d_q(i, k) = d_q(i, k) + y_d_q(j, k)
3448	#ifdef ISO
3449	DO ixt=1,ntraciso
3450	d_xt(ixt,i,k) = d_xt(ixt,i,k) + y_d_xt(ixt,j,k)
3451	ENDDO !DO ixt=1,ntraciso
3452	#endif
3453	d_u(i, k) = d_u(i, k) + y_d_u(j, k)
3454	d_v(i, k) = d_v(i, k) + y_d_v(j, k)
3455	ENDDO
3456	ENDDO
3457
3458	IF (ok_bs) THEN
3459	DO k = 1, klev
3460	DO j = 1, knon
3461	i = ni(j)
3462	d_qbs(i, k) = d_qbs(i, k) + y_d_qbs(j, k)
3463	ENDDO
3464	ENDDO
3465	ENDIF
3466
3467	#ifdef ISO
3468	#ifdef ISOVERIF
3469	! WRITE(,) 'd_q,d_xt(iso_eau,554,19)=',d_q(554,19),d_xt(iso_eau,554,19)
3470	! WRITE(,) 'pbl_surface 2929: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1)
3471	! WRITE(,) 'y_d_q,y_d_xt(iso_eau,2,1)=',y_d_q(2,1),y_d_xt(iso_eau,2,1)
3472	! WRITE(,) 'iso_eau.gt.0=',iso_eau.gt.0
3473	CALL iso_verif_noNaN_vect2D( &
3474	d_xt, &
3475	'pbl_surface 1385',ntraciso,klon,klev)
3476	IF (iso_eau >= 0) THEN
3477	CALL iso_verif_egalite_vect2D( &
3478	y_d_xt,y_d_q, &
3479	'pbl_surface_mod 2945',ntraciso,klon,klev)
3480	CALL iso_verif_egalite_vect2D( &
3481	d_xt,d_q, &
3482	'pbl_surface_mod 1276',ntraciso,klon,klev)
3483	ENDIF !IF (iso_eau >= 0) THEN
3484	#endif
3485	#endif
3486
3487	! PRINT*,'Dans pbl OK4'
3488
3489	IF (prt_level >=10) THEN
3490	print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
3491	d_t_w(1:knon, 1), d_t_x(1:knon, 1), d_t(1:knon, 1)
3492	ENDIF
3493
3494	IF (nsrf == is_oce .AND. activate_ocean_skin >= 1) THEN
3495	delta_sal = missing_val
3496	ds_ns = missing_val
3497	dt_ns = missing_val
3498	delta_sst = missing_val
3499	dter = missing_val
3500	dser = missing_val
3501	tkt = missing_val
3502	tks = missing_val
3503	taur = missing_val
3504	sss = missing_val
3505
3506	delta_sal(ni(:knon)) = ydelta_sal(:knon)
3507	ds_ns(ni(:knon)) = yds_ns(:knon)
3508	dt_ns(ni(:knon)) = ydt_ns(:knon)
3509	delta_sst(ni(:knon)) = ydelta_sst(:knon)
3510	dter(ni(:knon)) = ydter(:knon)
3511	dser(ni(:knon)) = ydser(:knon)
3512	tkt(ni(:knon)) = ytkt(:knon)
3513	tks(ni(:knon)) = ytks(:knon)
3514	taur(ni(:knon)) = ytaur(:knon)
3515	sss(ni(:knon)) = ysss(:knon)
3516
3517	IF (activate_ocean_skin == 2 .AND. type_ocean == "couple") THEN
3518	dt_ds = missing_val
3519	dt_ds(ni(:knon)) = ydt_ds(:knon)
3520	end if
3521	end if
3522
3523
3524	!****************************************************************************************
3525	! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m
3526	! Call HBTM
3527
3528	!****************************************************************************************
3529	!!!
3530
3531	#undef T2m
3532	#define T2m
3533	#ifdef T2m
3534	! Calculations of diagnostic t,q at 2m and u, v at 10m
3535
3536	! PRINT*,'Dans pbl OK41'
3537	! PRINT*,'tair1,yt(:,1),y_d_t(:,1)'
3538	! PRINT*, tair1,yt(:,1),y_d_t(:,1)
3539	!!! jyg le 07/02/2012
3540	IF (iflag_split ==0) THEN
3541	DO j = 1, knon
3542	uzon(j) = yu(j, 1) + y_d_u(j, 1)
3543	vmer(j) = yv(j, 1) + y_d_v(j, 1)
3544	tair1(j) = yt(j, 1) + y_d_t(j, 1) + y_d_t_diss(j, 1)
3545	qair1(j) = yq(j, 1) + y_d_q(j, 1)
3546	zgeo1(j) = RD * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, 1))) &
3547	* (ypaprs(j, 1) - ypplay(j, 1))
3548	tairsol(j) = yts(j) + y_d_ts(j)
3549	qairsol(j) = yqsurf(j)
3550	ENDDO
3551	ELSE ! (iflag_split .EQ.0)
3552	DO j = 1, knon
3553	uzon_x(j) = yu_x(j, 1) + y_d_u_x(j, 1)
3554	vmer_x(j) = yv_x(j, 1) + y_d_v_x(j, 1)
3555	tair1_x(j) = yt_x(j, 1) + y_d_t_x(j, 1) + y_d_t_diss_x(j, 1)
3556	qair1_x(j) = yq_x(j, 1) + y_d_q_x(j, 1)
3557	zgeo1_x(j) = RD * tair1_x(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, 1))) &
3558	* (ypaprs(j, 1) - ypplay(j, 1))
3559	tairsol(j) = yts(j) + y_d_ts(j)
3560	!! tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
3561	tairsol_x(j) = tairsol(j) - ywake_s(j) * y_delta_tsurf_new(j)
3562	qairsol(j) = yqsurf(j)
3563	ENDDO
3564	DO j = 1, knon
3565	uzon_w(j) = yu_w(j, 1) + y_d_u_w(j, 1)
3566	vmer_w(j) = yv_w(j, 1) + y_d_v_w(j, 1)
3567	tair1_w(j) = yt_w(j, 1) + y_d_t_w(j, 1) + y_d_t_diss_w(j, 1)
3568	qair1_w(j) = yq_w(j, 1) + y_d_q_w(j, 1)
3569	zgeo1_w(j) = RD * tair1_w(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, 1))) &
3570	* (ypaprs(j, 1) - ypplay(j, 1))
3571	tairsol_w(j) = tairsol(j) + (1. - ywake_s(j)) * y_delta_tsurf(j)
3572	qairsol(j) = yqsurf(j)
3573	ENDDO
3574	!!!
3575	ENDIF ! (iflag_split .EQ.0)
3576	!!!
3577	DO j = 1, knon
3578	! i = ni(j)
3579	! yz0h_oupas(j) = yz0m(j)
3580	! IF(nsrf.EQ.is_oce) THEN
3581	! yz0h_oupas(j) = z0m(i,nsrf)
3582	! ENDIF
3583	psfce(j) = ypaprs(j, 1)
3584	patm(j) = ypplay(j, 1)
3585	ENDDO
3586
3587	IF (iflag_pbl_surface_t2m_bug==1) THEN
3588	yz0h_oupas(1:knon) = yz0m(1:knon)
3589	ELSE
3590	yz0h_oupas(1:knon) = yz0h(1:knon)
3591	ENDIF
3592
3593	! PRINT*,'Dans pbl OK42A'
3594	! PRINT*,'tair1,yt(:,1),y_d_t(:,1)'
3595	! PRINT*, tair1,yt(:,1),y_d_t(:,1)
3596
3597	! Calculate the temperature and relative humidity at 2m and the wind at 10m
3598	!!! jyg le 07/02/2012
3599	IF (iflag_split ==0) THEN
3600	IF (iflag_new_t2mq2m==1) THEN
3601	CALL stdlevvarn(klon, knon, nsrf, zxli, &
3602	uzon, vmer, tair1, qair1, zgeo1, &
3603	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3604	yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
3605	yn2mout(:, nsrf, :))
3606	ELSE
3607	CALL stdlevvar(klon, knon, nsrf, zxli, &
3608	uzon, vmer, tair1, qair1, zgeo1, &
3609	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3610	yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
3611	ENDIF
3612	ELSE !(iflag_split .EQ.0)
3613	IF (iflag_new_t2mq2m==1) THEN
3614	CALL stdlevvarn(klon, knon, nsrf, zxli, &
3615	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
3616	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3617	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, &
3618	yn2mout_x(:, nsrf, :))
3619	CALL stdlevvarn(klon, knon, nsrf, zxli, &
3620	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
3621	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3622	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, &
3623	yn2mout_w(:, nsrf, :))
3624	ELSE
3625	CALL stdlevvar(klon, knon, nsrf, zxli, &
3626	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
3627	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3628	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, ypblh_x, rain_f, zxtsol)
3629	CALL stdlevvar(klon, knon, nsrf, zxli, &
3630	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
3631	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3632	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, ypblh_w, rain_f, zxtsol)
3633	ENDIF
3634	!!!
3635	ENDIF ! (iflag_split .EQ.0)
3636	!!!
3637	!!! jyg le 07/02/2012
3638	IF (iflag_split ==0) THEN
3639	DO j = 1, knon
3640	i = ni(j)
3641	t2m(i, nsrf) = yt2m(j)
3642	q2m(i, nsrf) = yq2m(j)
3643	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3644	ustar(i, nsrf) = yustar(j)
3645	u10m(i, nsrf) = (yu10m(j) * uzon(j)) / SQRT(uzon(j)2 + vmer(j)2)
3646	v10m(i, nsrf) = (yu10m(j) * vmer(j)) / SQRT(uzon(j)2 + vmer(j)2)
3647
3648	DO k = 1, 6
3649	n2mout(i, nsrf, k) = yn2mout(j, nsrf, k)
3650	END DO
3651
3652	ENDDO
3653	ELSE !(iflag_split .EQ.0)
3654	DO j = 1, knon
3655	i = ni(j)
3656	t2m_x(i, nsrf) = yt2m_x(j)
3657	q2m_x(i, nsrf) = yq2m_x(j)
3658	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3659	ustar_x(i, nsrf) = yustar_x(j)
3660	u10m_x(i, nsrf) = (yu10m_x(j) * uzon_x(j)) / SQRT(uzon_x(j)2 + vmer_x(j)2)
3661	v10m_x(i, nsrf) = (yu10m_x(j) * vmer_x(j)) / SQRT(uzon_x(j)2 + vmer_x(j)2)
3662
3663	DO k = 1, 6
3664	n2mout_x(i, nsrf, k) = yn2mout_x(j, nsrf, k)
3665	END DO
3666
3667	ENDDO
3668	DO j = 1, knon
3669	i = ni(j)
3670	t2m_w(i, nsrf) = yt2m_w(j)
3671	q2m_w(i, nsrf) = yq2m_w(j)
3672	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3673	ustar_w(i, nsrf) = yustar_w(j)
3674	u10m_w(i, nsrf) = (yu10m_w(j) * uzon_w(j)) / SQRT(uzon_w(j)2 + vmer_w(j)2)
3675	v10m_w(i, nsrf) = (yu10m_w(j) * vmer_w(j)) / SQRT(uzon_w(j)2 + vmer_w(j)2)
3676
3677	ustar(i, nsrf) = ustar_x(i, nsrf) + wake_s(i) * (ustar_w(i, nsrf) - ustar_x(i, nsrf))
3678	u10m(i, nsrf) = u10m_x(i, nsrf) + wake_s(i) * (u10m_w(i, nsrf) - u10m_x(i, nsrf))
3679	v10m(i, nsrf) = v10m_x(i, nsrf) + wake_s(i) * (v10m_w(i, nsrf) - v10m_x(i, nsrf))
3680
3681	DO k = 1, 6
3682	n2mout_w(i, nsrf, k) = yn2mout_w(j, nsrf, k)
3683	END DO
3684
3685	ENDDO
3686	!!!
3687	ENDIF ! (iflag_split .EQ.0)
3688	!!!
3689
3690	! PRINT*,'Dans pbl OK43'
3691	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
3692	!IM Ajoute dependance type surface
3693	IF (thermcep) THEN
3694	!!! jyg le 07/02/2012
3695	IF (iflag_split ==0) THEN
3696	DO j = 1, knon
3697	i = ni(j)
3698	zdelta1 = MAX(0., SIGN(1., rtt - yt2m(j)))
3699	zx_qs1 = r2es * FOEEW(yt2m(j), zdelta1) / paprs(i, 1)
3700	zx_qs1 = MIN(0.5, zx_qs1)
3701	zcor1 = 1. / (1. - RETV * zx_qs1)
3702	zx_qs1 = zx_qs1 * zcor1
3703
3704	rh2m(i) = rh2m(i) + yq2m(j) / zx_qs1 * pctsrf(i, nsrf)
3705	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i, nsrf)
3706	ENDDO
3707	ELSE ! (iflag_split .EQ.0)
3708	DO j = 1, knon
3709	i = ni(j)
3710	zdelta1 = MAX(0., SIGN(1., rtt - yt2m_x(j)))
3711	zx_qs1 = r2es * FOEEW(yt2m_x(j), zdelta1) / paprs(i, 1)
3712	zx_qs1 = MIN(0.5, zx_qs1)
3713	zcor1 = 1. / (1. - RETV * zx_qs1)
3714	zx_qs1 = zx_qs1 * zcor1
3715
3716	rh2m_x(i) = rh2m_x(i) + yq2m_x(j) / zx_qs1 * pctsrf(i, nsrf)
3717	qsat2m_x(i) = qsat2m_x(i) + zx_qs1 * pctsrf(i, nsrf)
3718	ENDDO
3719	DO j = 1, knon
3720	i = ni(j)
3721	zdelta1 = MAX(0., SIGN(1., rtt - yt2m_w(j)))
3722	zx_qs1 = r2es * FOEEW(yt2m_w(j), zdelta1) / paprs(i, 1)
3723	zx_qs1 = MIN(0.5, zx_qs1)
3724	zcor1 = 1. / (1. - RETV * zx_qs1)
3725	zx_qs1 = zx_qs1 * zcor1
3726
3727	rh2m_w(i) = rh2m_w(i) + yq2m_w(j) / zx_qs1 * pctsrf(i, nsrf)
3728	qsat2m_w(i) = qsat2m_w(i) + zx_qs1 * pctsrf(i, nsrf)
3729	ENDDO
3730	!!!
3731	ENDIF ! (iflag_split .EQ.0)
3732	!!!
3733	ENDIF
3734
3735	IF (prt_level >=10) THEN
3736	print *, 'T2m, q2m, RH2m ', &
3737	t2m(1:knon, :), q2m(1:knon, :), rh2m(1:knon)
3738	ENDIF
3739
3740	! PRINT*,'OK pbl 5'
3741
3742	!!! jyg le 07/02/2012
3743	IF (iflag_split ==0) THEN
3744	CALL hbtm(knon, ypaprs, ypplay, &
3745	yt2m, yt10m, yq2m, yq10m, yustar, ywstar, &
3746	y_flux_t, y_flux_q, yu, yv, yt, yq, &
3747	ypblh, ycapCL, yoliqCL, ycteiCL, ypblT, &
3748	ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
3749	IF (prt_level >=10) THEN
3750	print *, ' Arg. de HBTM: yt2m ', yt2m(1:knon)
3751	print *, ' Arg. de HBTM: yt10m ', yt10m(1:knon)
3752	print *, ' Arg. de HBTM: yq2m ', yq2m(1:knon)
3753	print *, ' Arg. de HBTM: yq10m ', yq10m(1:knon)
3754	print *, ' Arg. de HBTM: yustar ', yustar(1:knon)
3755	print *, ' Arg. de HBTM: y_flux_t ', y_flux_t(1:knon, :)
3756	print *, ' Arg. de HBTM: y_flux_q ', y_flux_q(1:knon, :)
3757	print *, ' Arg. de HBTM: yu ', yu(1:knon, :)
3758	print *, ' Arg. de HBTM: yv ', yv(1:knon, :)
3759	print *, ' Arg. de HBTM: yt ', yt(1:knon, :)
3760	print *, ' Arg. de HBTM: yq ', yq(1:knon, :)
3761	ENDIF
3762	ELSE ! (iflag_split .EQ.0)
3763	CALL HBTM(knon, ypaprs, ypplay, &
3764	yt2m_x, yt10m_x, yq2m_x, yq10m_x, yustar_x, ywstar_x, &
3765	y_flux_t_x, y_flux_q_x, yu_x, yv_x, yt_x, yq_x, &
3766	ypblh_x, ycapCL_x, yoliqCL_x, ycteiCL_x, ypblT_x, &
3767	ytherm_x, ytrmb1_x, ytrmb2_x, ytrmb3_x, ylcl_x)
3768	IF (prt_level >=10) THEN
3769	print *, ' Arg. de HBTM: yt2m_x ', yt2m_x(1:knon)
3770	print *, ' Arg. de HBTM: yt10m_x ', yt10m_x(1:knon)
3771	print *, ' Arg. de HBTM: yq2m_x ', yq2m_x(1:knon)
3772	print *, ' Arg. de HBTM: yq10m_x ', yq10m_x(1:knon)
3773	print *, ' Arg. de HBTM: yustar_x ', yustar_x(1:knon)
3774	print *, ' Arg. de HBTM: y_flux_t_x ', y_flux_t_x(1:knon, :)
3775	print *, ' Arg. de HBTM: y_flux_q_x ', y_flux_q_x(1:knon, :)
3776	print *, ' Arg. de HBTM: yu_x ', yu_x(1:knon, :)
3777	print *, ' Arg. de HBTM: yv_x ', yv_x(1:knon, :)
3778	print *, ' Arg. de HBTM: yt_x ', yt_x(1:knon, :)
3779	print *, ' Arg. de HBTM: yq_x ', yq_x(1:knon, :)
3780	ENDIF
3781	CALL HBTM(knon, ypaprs, ypplay, &
3782	yt2m_w, yt10m_w, yq2m_w, yq10m_w, yustar_w, ywstar_w, &
3783	y_flux_t_w, y_flux_q_w, yu_w, yv_w, yt_w, yq_w, &
3784	ypblh_w, ycapCL_w, yoliqCL_w, ycteiCL_w, ypblT_w, &
3785	ytherm_w, ytrmb1_w, ytrmb2_w, ytrmb3_w, ylcl_w)
3786	!!!
3787	ENDIF ! (iflag_split .EQ.0)
3788	!!!
3789
3790	!!! jyg le 07/02/2012
3791	IF (iflag_split ==0) THEN
3792	!!!
3793	DO j = 1, knon
3794	i = ni(j)
3795	pblh(i, nsrf) = ypblh(j)
3796	wstar(i, nsrf) = ywstar(j)
3797	plcl(i, nsrf) = ylcl(j)
3798	capCL(i, nsrf) = ycapCL(j)
3799	oliqCL(i, nsrf) = yoliqCL(j)
3800	cteiCL(i, nsrf) = ycteiCL(j)
3801	pblT(i, nsrf) = ypblT(j)
3802	therm(i, nsrf) = ytherm(j)
3803	trmb1(i, nsrf) = ytrmb1(j)
3804	trmb2(i, nsrf) = ytrmb2(j)
3805	trmb3(i, nsrf) = ytrmb3(j)
3806	ENDDO
3807	IF (prt_level >=10) THEN
3808	print *, 'After HBTM: pblh ', pblh(1:knon, :)
3809	print *, 'After HBTM: plcl ', plcl(1:knon, :)
3810	print *, 'After HBTM: cteiCL ', cteiCL(1:knon, :)
3811	ENDIF
3812	ELSE !(iflag_split .EQ.0)
3813	DO j = 1, knon
3814	i = ni(j)
3815	pblh_x(i, nsrf) = ypblh_x(j)
3816	wstar_x(i, nsrf) = ywstar_x(j)
3817	plcl_x(i, nsrf) = ylcl_x(j)
3818	capCL_x(i, nsrf) = ycapCL_x(j)
3819	oliqCL_x(i, nsrf) = yoliqCL_x(j)
3820	cteiCL_x(i, nsrf) = ycteiCL_x(j)
3821	pblT_x(i, nsrf) = ypblT_x(j)
3822	therm_x(i, nsrf) = ytherm_x(j)
3823	trmb1_x(i, nsrf) = ytrmb1_x(j)
3824	trmb2_x(i, nsrf) = ytrmb2_x(j)
3825	trmb3_x(i, nsrf) = ytrmb3_x(j)
3826	ENDDO
3827	IF (prt_level >=10) THEN
3828	print *, 'After HBTM: pblh_x ', pblh_x(1:knon, :)
3829	print *, 'After HBTM: plcl_x ', plcl_x(1:knon, :)
3830	print *, 'After HBTM: cteiCL_x ', cteiCL_x(1:knon, :)
3831	ENDIF
3832	DO j = 1, knon
3833	i = ni(j)
3834	pblh_w(i, nsrf) = ypblh_w(j)
3835	wstar_w(i, nsrf) = ywstar_w(j)
3836	plcl_w(i, nsrf) = ylcl_w(j)
3837	capCL_w(i, nsrf) = ycapCL_w(j)
3838	oliqCL_w(i, nsrf) = yoliqCL_w(j)
3839	cteiCL_w(i, nsrf) = ycteiCL_w(j)
3840	pblT_w(i, nsrf) = ypblT_w(j)
3841	therm_w(i, nsrf) = ytherm_w(j)
3842	trmb1_w(i, nsrf) = ytrmb1_w(j)
3843	trmb2_w(i, nsrf) = ytrmb2_w(j)
3844	trmb3_w(i, nsrf) = ytrmb3_w(j)
3845	ENDDO
3846	IF (prt_level >=10) THEN
3847	print *, 'After HBTM: pblh_w ', pblh_w(1:knon, :)
3848	print *, 'After HBTM: plcl_w ', plcl_w(1:knon, :)
3849	print *, 'After HBTM: cteiCL_w ', cteiCL_w(1:knon, :)
3850	ENDIF
3851	!!!
3852	ENDIF ! (iflag_split .EQ.0)
3853	!!!
3854
3855	! PRINT*,'OK pbl 6'
3856	#else
3857	! T2m not defined
3858	! No calculation
3859	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
3860	#endif
3861
3862	!****************************************************************************************
3863	! 15) End of loop over different surfaces
3864
3865	!****************************************************************************************
3866	ENDDO loop_nbsrf
3867
3868	!----------------------------------------------------------------------------------------
3869	! Reset iflag_split
3870
3871	iflag_split = iflag_split_ref
3872
3873	#ifdef ISO
3874	#ifdef ISOVERIF
3875	! WRITE(,) 'pbl_surface tmp 3249: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1)
3876	IF (iso_eau >= 0) THEN
3877	CALL iso_verif_egalite_vect2D( &
3878	d_xt,d_q, &
3879	'pbl_surface_mod 1276',ntraciso,klon,klev)
3880	ENDIF !IF (iso_eau >= 0) THEN
3881	#endif
3882	#endif
3883
3884	!****************************************************************************************
3885	! 16) Calculate the mean value over all sub-surfaces for some variables
3886
3887	!****************************************************************************************
3888
3889	z0m(:, nbsrf + 1) = 0.0
3890	z0h(:, nbsrf + 1) = 0.0
3891	DO nsrf = 1, nbsrf
3892	DO i = 1, klon
3893	z0m(i, nbsrf + 1) = z0m(i, nbsrf + 1) + z0m(i, nsrf) * pctsrf(i, nsrf)
3894	z0h(i, nbsrf + 1) = z0h(i, nbsrf + 1) + z0h(i, nsrf) * pctsrf(i, nsrf)
3895	ENDDO
3896	ENDDO
3897
3898	! PRINT*,'OK pbl 7'
3899	zxfluxt(:, :) = 0.0 ; zxfluxq(:, :) = 0.0
3900	zxfluxu(:, :) = 0.0 ; zxfluxv(:, :) = 0.0
3901	zxfluxt_x(:, :) = 0.0 ; zxfluxq_x(:, :) = 0.0
3902	zxfluxu_x(:, :) = 0.0 ; zxfluxv_x(:, :) = 0.0
3903	zxfluxt_w(:, :) = 0.0 ; zxfluxq_w(:, :) = 0.0
3904	zxfluxu_w(:, :) = 0.0 ; zxfluxv_w(:, :) = 0.0
3905	#ifdef ISO
3906	zxfluxxt(:,:,:) = 0.0
3907	zxfluxxt_x(:,:,:) = 0.0
3908	zxfluxxt_w(:,:,:) = 0.0
3909	#endif
3910
3911
3912	!!! jyg le 07/02/2012
3913	IF (iflag_split >=1) THEN
3914	!!!
3915	!!! nrlmd & jyg les 02/05/2011, 05/02/2012
3916
3917	DO nsrf = 1, nbsrf
3918	DO k = 1, klev
3919	DO i = 1, klon
3920	zxfluxt_x(i, k) = zxfluxt_x(i, k) + flux_t_x(i, k, nsrf) * pctsrf(i, nsrf)
3921	zxfluxq_x(i, k) = zxfluxq_x(i, k) + flux_q_x(i, k, nsrf) * pctsrf(i, nsrf)
3922	zxfluxu_x(i, k) = zxfluxu_x(i, k) + flux_u_x(i, k, nsrf) * pctsrf(i, nsrf)
3923	zxfluxv_x(i, k) = zxfluxv_x(i, k) + flux_v_x(i, k, nsrf) * pctsrf(i, nsrf)
3924
3925	zxfluxt_w(i, k) = zxfluxt_w(i, k) + flux_t_w(i, k, nsrf) * pctsrf(i, nsrf)
3926	zxfluxq_w(i, k) = zxfluxq_w(i, k) + flux_q_w(i, k, nsrf) * pctsrf(i, nsrf)
3927	zxfluxu_w(i, k) = zxfluxu_w(i, k) + flux_u_w(i, k, nsrf) * pctsrf(i, nsrf)
3928	zxfluxv_w(i, k) = zxfluxv_w(i, k) + flux_v_w(i, k, nsrf) * pctsrf(i, nsrf)
3929	#ifdef ISO
3930	DO ixt=1,ntraciso
3931	zxfluxxt_x(ixt,i,k) = zxfluxxt_x(ixt,i,k) + flux_xt_x(ixt,i,k,nsrf) * pctsrf(i,nsrf)
3932	zxfluxxt_w(ixt,i,k) = zxfluxxt_w(ixt,i,k) + flux_xt_w(ixt,i,k,nsrf) * pctsrf(i,nsrf)
3933	ENDDO ! DO ixt=1,ntraciso
3934	#endif
3935	ENDDO
3936	ENDDO
3937	ENDDO
3938
3939	DO i = 1, klon
3940	zxsens_x(i) = - zxfluxt_x(i, 1)
3941	zxsens_w(i) = - zxfluxt_w(i, 1)
3942	ENDDO
3943	!!!
3944	ENDIF ! (iflag_split .ge.1)
3945	!!!
3946
3947	DO nsrf = 1, nbsrf
3948	DO k = 1, klev
3949	DO i = 1, klon
3950	zxfluxt(i, k) = zxfluxt(i, k) + flux_t(i, k, nsrf) * pctsrf(i, nsrf)
3951	zxfluxq(i, k) = zxfluxq(i, k) + flux_q(i, k, nsrf) * pctsrf(i, nsrf)
3952	zxfluxu(i, k) = zxfluxu(i, k) + flux_u(i, k, nsrf) * pctsrf(i, nsrf)
3953	zxfluxv(i, k) = zxfluxv(i, k) + flux_v(i, k, nsrf) * pctsrf(i, nsrf)
3954	#ifdef ISO
3955	DO ixt=1,niso
3956	zxfluxxt(ixt,i,k) = zxfluxxt(ixt,i,k) + flux_xt(ixt,i,k,nsrf) * pctsrf(i,nsrf)
3957	ENDDO ! DO ixt=1,niso
3958	#endif
3959	ENDDO
3960	ENDDO
3961	ENDDO
3962
3963	DO i = 1, klon
3964	zxsens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
3965	zxevap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol
3966	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
3967	ENDDO
3968
3969	! if blowing snow
3970	IF (ok_bs) THEN
3971	DO nsrf = 1, nbsrf
3972	DO k = 1, klev
3973	DO i = 1, klon
3974	zxfluxqbs(i, k) = zxfluxqbs(i, k) + flux_qbs(i, k, nsrf) * pctsrf(i, nsrf)
3975	ENDDO
3976	ENDDO
3977	ENDDO
3978
3979	DO i = 1, klon
3980	zxsnowerosion(i) = zxfluxqbs(i, 1) ! blowings snow flux at the surface
3981	END DO
3982	endif
3983
3984	#ifdef ISO
3985	DO i = 1, klon
3986	DO ixt=1,ntraciso
3987	zxxtevap(ixt,i) = - zxfluxxt(ixt,i,1)
3988	ENDDO
3989	ENDDO
3990	#endif
3991
3992	!!!
3993
3994	! Incrementer la temperature du sol
3995
3996	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
3997	zt2m(:) = 0.0 ; zq2m(:) = 0.0 ; zn2mout(:, :) = 0
3998	zustar(:) = 0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
3999	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
4000	!!! jyg le 07/02/2012
4001	s_pblh_x(:) = 0.0 ; s_plcl_x(:) = 0.0
4002	s_pblh_w(:) = 0.0 ; s_plcl_w(:) = 0.0
4003	!!!
4004	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
4005	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
4006	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
4007	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
4008	wstar(:, is_ave) = 0.
4009
4010	! PRINT*,'OK pbl 9'
4011
4012	!!! nrlmd le 02/05/2011
4013	zxfluxlat_x(:) = 0.0 ; zxfluxlat_w(:) = 0.0
4014	!!!
4015
4016	DO nsrf = 1, nbsrf
4017	DO i = 1, klon
4018	ts(i, nsrf) = ts(i, nsrf) + d_ts(i, nsrf)
4019
4020	wfbils(i, nsrf) = (solsw(i, nsrf) + sollw(i, nsrf) &
4021	+ flux_t(i, 1, nsrf) + fluxlat(i, nsrf)) * pctsrf(i, nsrf)
4022
4023	wfevap(i, nsrf) = evap(i, nsrf) * pctsrf(i, nsrf)
4024
4025	zxtsol(i) = zxtsol(i) + ts(i, nsrf) * pctsrf(i, nsrf)
4026	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf) * pctsrf(i, nsrf)
4027	ENDDO
4028	ENDDO
4029
4030	!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
4031	IF (iflag_order2_sollw == 1) THEN
4032	meansqT(:) = 0. ! as working buffer
4033	DO nsrf = 1, nbsrf
4034	DO i = 1, klon
4035	meansqT(i) = meansqT(i) + (ts(i, nsrf) - zxtsol(i))*2 pctsrf(i, nsrf)
4036	ENDDO
4037	ENDDO
4038	zxtsol(:) = zxtsol(:) + 1.5 * meansqT(:) / zxtsol(:)
4039	ENDIF ! iflag_order2_sollw == 1
4040	!>al1
4041
4042	!!! jyg le 07/02/2012
4043	IF (iflag_split ==0) THEN
4044	DO nsrf = 1, nbsrf
4045	DO i = 1, klon
4046	zt2m(i) = zt2m(i) + t2m(i, nsrf) * pctsrf(i, nsrf)
4047	zq2m(i) = zq2m(i) + q2m(i, nsrf) * pctsrf(i, nsrf)
4048
4049	DO k = 1, 6
4050	zn2mout(i, k) = zn2mout(i, k) + n2mout(i, nsrf, k) * pctsrf(i, nsrf)
4051	ENDDO
4052
4053	zustar(i) = zustar(i) + ustar(i, nsrf) * pctsrf(i, nsrf)
4054	wstar(i, is_ave) = wstar(i, is_ave) + wstar(i, nsrf) * pctsrf(i, nsrf)
4055	zu10m(i) = zu10m(i) + u10m(i, nsrf) * pctsrf(i, nsrf)
4056	zv10m(i) = zv10m(i) + v10m(i, nsrf) * pctsrf(i, nsrf)
4057
4058	s_pblh(i) = s_pblh(i) + pblh(i, nsrf) * pctsrf(i, nsrf)
4059	s_plcl(i) = s_plcl(i) + plcl(i, nsrf) * pctsrf(i, nsrf)
4060	s_capCL(i) = s_capCL(i) + capCL(i, nsrf) * pctsrf(i, nsrf)
4061	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) * pctsrf(i, nsrf)
4062	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) * pctsrf(i, nsrf)
4063	s_pblT(i) = s_pblT(i) + pblT(i, nsrf) * pctsrf(i, nsrf)
4064	s_therm(i) = s_therm(i) + therm(i, nsrf) * pctsrf(i, nsrf)
4065	s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) * pctsrf(i, nsrf)
4066	s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) * pctsrf(i, nsrf)
4067	s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) * pctsrf(i, nsrf)
4068	ENDDO
4069	ENDDO
4070	ELSE !(iflag_split .EQ.0)
4071	DO nsrf = 1, nbsrf
4072	DO i = 1, klon
4073	!!! nrlmd le 02/05/2011
4074	zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i, nsrf) * pctsrf(i, nsrf)
4075	zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i, nsrf) * pctsrf(i, nsrf)
4076	!!!
4077	!!! jyg le 08/02/2012
4078	!! Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;
4079	!! pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
4080	!! pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
4081	!! pour les autres variables, on sort les valeurs de la region (x).
4082	zt2m(i) = zt2m(i) + (t2m_x(i, nsrf) + wake_s(i) * (t2m_w(i, nsrf) - t2m_x(i, nsrf))) * pctsrf(i, nsrf)
4083	zq2m(i) = zq2m(i) + q2m_x(i, nsrf) * pctsrf(i, nsrf)
4084
4085	DO k = 1, 6
4086	zn2mout(i, k) = zn2mout(i, k) + n2mout_x(i, nsrf, k) * pctsrf(i, nsrf)
4087	ENDDO
4088
4089	zustar(i) = zustar(i) + ustar_x(i, nsrf) * pctsrf(i, nsrf)
4090	wstar(i, is_ave) = wstar(i, is_ave) + wstar_x(i, nsrf) * pctsrf(i, nsrf)
4091	zu10m(i) = zu10m(i) + u10m_x(i, nsrf) * pctsrf(i, nsrf)
4092	zv10m(i) = zv10m(i) + v10m_x(i, nsrf) * pctsrf(i, nsrf)
4093
4094	s_pblh(i) = s_pblh(i) + pblh_x(i, nsrf) * pctsrf(i, nsrf)
4095	s_pblh_x(i) = s_pblh_x(i) + pblh_x(i, nsrf) * pctsrf(i, nsrf)
4096	s_pblh_w(i) = s_pblh_w(i) + pblh_w(i, nsrf) * pctsrf(i, nsrf)
4097
4098	s_plcl(i) = s_plcl(i) + plcl_x(i, nsrf) * pctsrf(i, nsrf)
4099	s_plcl_x(i) = s_plcl_x(i) + plcl_x(i, nsrf) * pctsrf(i, nsrf)
4100	s_plcl_w(i) = s_plcl_w(i) + plcl_w(i, nsrf) * pctsrf(i, nsrf)
4101
4102	s_capCL(i) = s_capCL(i) + capCL_x(i, nsrf) * pctsrf(i, nsrf)
4103	s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i, nsrf) * pctsrf(i, nsrf)
4104	s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i, nsrf) * pctsrf(i, nsrf)
4105	s_pblT(i) = s_pblT(i) + pblT_x(i, nsrf) * pctsrf(i, nsrf)
4106	s_therm(i) = s_therm(i) + therm_x(i, nsrf) * pctsrf(i, nsrf)
4107	s_trmb1(i) = s_trmb1(i) + trmb1_x(i, nsrf) * pctsrf(i, nsrf)
4108	s_trmb2(i) = s_trmb2(i) + trmb2_x(i, nsrf) * pctsrf(i, nsrf)
4109	s_trmb3(i) = s_trmb3(i) + trmb3_x(i, nsrf) * pctsrf(i, nsrf)
4110	ENDDO
4111	ENDDO
4112	DO i = 1, klon
4113	qsat2m(i) = qsat2m_x(i) + wake_s(i) * (qsat2m_x(i) - qsat2m_w(i))
4114	ENDDO
4115	!!!
4116	ENDIF ! (iflag_split .EQ.0)
4117	!!!
4118
4119	IF (check) THEN
4120	amn = MIN(ts(1, is_ter), 1000.)
4121	amx = MAX(ts(1, is_ter), -1000.)
4122	DO i = 2, klon
4123	amn = MIN(ts(i, is_ter), amn)
4124	amx = MAX(ts(i, is_ter), amx)
4125	ENDDO
4126	PRINT*, ' debut apres d_ts min max ftsol(ts)', itap, amn, amx
4127	ENDIF
4128
4129	!jg ?
4130	!!$!
4131	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
4132	!!$! sub-surfaces is distributed.
4133	!!$!
4134	!!$ DO nsrf = 1, nbsrf
4135	!!$ DO i = 1, klon
4136	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
4137	!!$ ts(i,nsrf) = zxtsol(i)
4138	!!$ t2m(i,nsrf) = zt2m(i)
4139	!!$ q2m(i,nsrf) = zq2m(i)
4140	!!$ u10m(i,nsrf) = zu10m(i)
4141	!!$ v10m(i,nsrf) = zv10m(i)
4142	!!$
4143	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
4144	!!$ pblh(i,nsrf) = s_pblh(i)
4145	!!$ plcl(i,nsrf) = s_plcl(i)
4146	!!$ capCL(i,nsrf) = s_capCL(i)
4147	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
4148	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
4149	!!$ pblT(i,nsrf) = s_pblT(i)
4150	!!$ therm(i,nsrf) = s_therm(i)
4151	!!$ trmb1(i,nsrf) = s_trmb1(i)
4152	!!$ trmb2(i,nsrf) = s_trmb2(i)
4153	!!$ trmb3(i,nsrf) = s_trmb3(i)
4154	!!$ ENDIF
4155	!!$ ENDDO
4156	!!$ ENDDO
4157
4158	DO i = 1, klon
4159	fder(i) = - 4.0 * RSIGMA * zxtsol(i)**3
4160	ENDDO
4161
4162	zxqsurf(:) = 0.0
4163	zxsnow(:) = 0.0
4164	#ifdef ISO
4165	zxxtsnow(:,:) = 0.0
4166	#endif
4167
4168	DO nsrf = 1, nbsrf
4169	DO i = 1, klon
4170	zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i, nsrf), 0.0) * pctsrf(i, nsrf)
4171	zxsnow(i) = zxsnow(i) + snow(i, nsrf) * pctsrf(i, nsrf)
4172	#ifdef ISO
4173	DO ixt=1,niso
4174	zxxtsnow(ixt,i) = zxxtsnow(ixt,i) + xtsnow(ixt,i,nsrf) * pctsrf(i,nsrf)
4175	ENDDO ! DO ixt=1,niso
4176	#endif
4177	ENDDO
4178	ENDDO
4179
4180	! Premier niveau de vent sortie dans physiq.F
4181	zu1(:) = u(:, 1)
4182	zv1(:) = v(:, 1)
4183
4184	END SUBROUTINE pbl_surface
4185
4186	!****************************************************************************************
4187
4188	SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst &
4189	#ifdef ISO
4190	,xtsnow_rst,Rland_ice_rst &
4191	#endif
4192	)
4193
4194	USE indice_sol_mod
4195	#ifdef ISO
4196	#ifdef ISOVERIF
4197	USE isotopes_mod, ONLY: iso_eau,ridicule
4198	USE isotopes_verif_mod, ONLY: errmax,errmaxrel
4199	#endif
4200	#endif
4201
4202	INCLUDE "dimsoil.h"
4203
4204	! Ouput variables
4205	!****************************************************************************************
4206	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
4207	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
4208	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
4209	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
4210	#ifdef ISO
4211	REAL, DIMENSION(niso,klon, nbsrf), INTENT(OUT) :: xtsnow_rst
4212	REAL, DIMENSION(niso,klon), INTENT(OUT) :: Rland_ice_rst
4213	#endif
4214
4215
4216	!****************************************************************************************
4217	! Return module variables for writing to restart file
4218
4219	!****************************************************************************************
4220	fder_rst(:) = fder(:)
4221	snow_rst(:, :) = snow(:, :)
4222	qsurf_rst(:, :) = qsurf(:, :)
4223	ftsoil_rst(:, :, :) = ftsoil(:, :, :)
4224	#ifdef ISO
4225	xtsnow_rst(:,:,:) = xtsnow(:,:,:)
4226	Rland_ice_rst(:,:) = Rland_ice(:,:)
4227	#endif
4228
4229	!****************************************************************************************
4230	! Deallocate module variables
4231
4232	!****************************************************************************************
4233	! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
4234	IF (ALLOCATED(fder)) DEALLOCATE(fder)
4235	IF (ALLOCATED(snow)) DEALLOCATE(snow)
4236	IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
4237	IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
4238	IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
4239	IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
4240	#ifdef ISO
4241	IF (ALLOCATED(xtsnow)) DEALLOCATE(xtsnow)
4242	IF (ALLOCATED(Rland_ice)) DEALLOCATE(Rland_ice)
4243	IF (ALLOCATED(Roce)) DEALLOCATE(Roce)
4244	#endif
4245
4246	!jyg<
4247	!****************************************************************************************
4248	! Deallocate variables for pbl splitting
4249
4250	!****************************************************************************************
4251
4252	CALL wx_pbl_final
4253	!>jyg
4254
4255	END SUBROUTINE pbl_surface_final
4256
4257	!****************************************************************************************
4258
4259	!albedo SB >>>
4260	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
4261	evap, z0m, z0h, agesno, &
4262	tsurf, alb_dir, alb_dif, ustar, u10m, v10m, tke &
4263	#ifdef ISO
4264	,xtevap &
4265	#endif
4266	&)
4267	!albedo SB <<<
4268	! Give default values where new fraction has appread
4269
4270	USE indice_sol_mod
4271	USE phys_state_var_mod, ONLY: delta_sal, ds_ns, dt_ns, delta_sst, dter, &
4272	dser, dt_ds
4273	USE config_ocean_skin_m, ONLY: activate_ocean_skin
4274
4275	INCLUDE "dimsoil.h"
4276	INCLUDE "clesphys.h"
4277	INCLUDE "compbl.h"
4278
4279	! Input variables
4280	!****************************************************************************************
4281	INTEGER, INTENT(IN) :: itime
4282	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
4283
4284	! InOutput variables
4285	!****************************************************************************************
4286	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: tsurf
4287	!albedo SB >>>
4288	REAL, DIMENSION(klon, nsw, nbsrf), INTENT(INOUT) :: alb_dir, alb_dif
4289	INTEGER :: k
4290	!albedo SB <<<
4291	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ustar, u10m, v10m
4292	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: evap, agesno
4293	REAL, DIMENSION(klon, nbsrf + 1), INTENT(INOUT) :: z0m, z0h
4294	REAL, DIMENSION(klon, klev + 1, nbsrf + 1), INTENT(INOUT) :: tke
4295	#ifdef ISO
4296	REAL, DIMENSION(ntraciso,klon,nbsrf), INTENT(INOUT) :: xtevap
4297	#endif
4298
4299	! Local variables
4300	!****************************************************************************************
4301	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
4302	CHARACTER(len = 80) :: abort_message
4303	CHARACTER(len = 20) :: modname = 'pbl_surface_newfrac'
4304	INTEGER, DIMENSION(nbsrf) :: nfois = 0, mfois = 0, pfois = 0
4305	#ifdef ISO
4306	INTEGER :: ixt
4307	#endif
4308
4309	! All at once !!
4310	!****************************************************************************************
4311
4312	DO nsrf = 1, nbsrf
4313	! First decide complement sub-surfaces
4314	SELECT CASE (nsrf)
4315	CASE(is_oce)
4316	nsrf_comp1 = is_sic
4317	nsrf_comp2 = is_ter
4318	nsrf_comp3 = is_lic
4319	CASE(is_sic)
4320	nsrf_comp1 = is_oce
4321	nsrf_comp2 = is_ter
4322	nsrf_comp3 = is_lic
4323	CASE(is_ter)
4324	nsrf_comp1 = is_lic
4325	nsrf_comp2 = is_oce
4326	nsrf_comp3 = is_sic
4327	CASE(is_lic)
4328	nsrf_comp1 = is_ter
4329	nsrf_comp2 = is_oce
4330	nsrf_comp3 = is_sic
4331	END SELECT
4332
4333	! Initialize all new fractions
4334	DO i = 1, klon
4335	IF (pctsrf_new(i, nsrf) > 0. .AND. pctsrf_old(i, nsrf) == 0.) THEN
4336
4337	IF (pctsrf_old(i, nsrf_comp1) > 0.) THEN
4338	! Use the complement sub-surface, keeping the continents unchanged
4339	qsurf(i, nsrf) = qsurf(i, nsrf_comp1)
4340	evap(i, nsrf) = evap(i, nsrf_comp1)
4341	z0m(i, nsrf) = z0m(i, nsrf_comp1)
4342	z0h(i, nsrf) = z0h(i, nsrf_comp1)
4343	tsurf(i, nsrf) = tsurf(i, nsrf_comp1)
4344	!albedo SB >>>
4345	DO k = 1, nsw
4346	alb_dir(i, k, nsrf) = alb_dir(i, k, nsrf_comp1)
4347	alb_dif(i, k, nsrf) = alb_dif(i, k, nsrf_comp1)
4348	ENDDO
4349	!albedo SB <<<
4350	ustar(i, nsrf) = ustar(i, nsrf_comp1)
4351	u10m(i, nsrf) = u10m(i, nsrf_comp1)
4352	v10m(i, nsrf) = v10m(i, nsrf_comp1)
4353	#ifdef ISO
4354	DO ixt=1,ntraciso
4355	xtevap(ixt,i,nsrf) = xtevap(ixt,i,nsrf_comp1)
4356	ENDDO
4357	#endif
4358	IF (iflag_pbl > 1) THEN
4359	tke(i, :, nsrf) = tke(i, :, nsrf_comp1)
4360	ENDIF
4361	mfois(nsrf) = mfois(nsrf) + 1
4362	! F. Codron sensible default values for ocean and sea ice
4363	IF (nsrf==is_oce) THEN
4364	tsurf(i, nsrf) = 271.35
4365	! (temperature of sea water under sea ice, so that
4366	! is also the temperature of appearing sea water)
4367	DO k = 1, nsw
4368	alb_dir(i, k, nsrf) = 0.06 ! typical Ocean albedo
4369	alb_dif(i, k, nsrf) = 0.06
4370	ENDDO
4371	IF (activate_ocean_skin >= 1) THEN
4372	IF (activate_ocean_skin == 2 &
4373	.AND. type_ocean == "couple") THEN
4374	delta_sal(i) = 0.
4375	delta_sst(i) = 0.
4376	dter(i) = 0.
4377	dser(i) = 0.
4378	dt_ds(i) = 0.
4379	end if
4380
4381	ds_ns(i) = 0.
4382	dt_ns(i) = 0.
4383	end if
4384	ELSE IF (nsrf==is_sic) THEN
4385	tsurf(i, nsrf) = 271.35
4386	! (Temperature at base of sea ice. Surface
4387	! temperature could be higher, up to 0 Celsius
4388	! degrees. We set it to -1.8 Celsius degrees for
4389	! consistency with the ocean slab model.)
4390	DO k = 1, nsw
4391	alb_dir(i, k, nsrf) = 0.3 ! thin ice
4392	alb_dif(i, k, nsrf) = 0.3
4393	ENDDO
4394	ENDIF
4395	ELSE
4396	! The continents have changed. The new fraction receives the mean sum of the existent fractions
4397	qsurf(i, nsrf) = qsurf(i, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + qsurf(i, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4398	evap(i, nsrf) = evap(i, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + evap(i, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4399	z0m(i, nsrf) = z0m(i, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + z0m(i, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4400	z0h(i, nsrf) = z0h(i, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + z0h(i, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4401	tsurf(i, nsrf) = tsurf(i, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + tsurf(i, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4402	!albedo SB >>>
4403	DO k = 1, nsw
4404	alb_dir(i, k, nsrf) = alb_dir(i, k, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + &
4405	alb_dir(i, k, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4406	alb_dif(i, k, nsrf) = alb_dif(i, k, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + &
4407	alb_dif(i, k, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4408	ENDDO
4409	!albedo SB <<<
4410	ustar(i, nsrf) = ustar(i, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + ustar(i, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4411	u10m(i, nsrf) = u10m(i, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + u10m(i, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4412	v10m(i, nsrf) = v10m(i, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + v10m(i, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4413	#ifdef ISO
4414	DO ixt=1,ntraciso
4415	xtevap(ixt,i,nsrf) = xtevap(ixt,i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) &
4416	+ xtevap(ixt,i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
4417	ENDDO
4418	#endif
4419	IF (iflag_pbl > 1) THEN
4420	tke(i, :, nsrf) = tke(i, :, nsrf_comp2) * pctsrf_old(i, nsrf_comp2) + tke(i, :, nsrf_comp3) * pctsrf_old(i, nsrf_comp3)
4421	ENDIF
4422
4423	! Security abort. This option has never been tested. To test, comment the following line.
4424	! abort_message='The fraction of the continents have changed!'
4425	! CALL abort_physic(modname,abort_message,1)
4426	nfois(nsrf) = nfois(nsrf) + 1
4427	ENDIF
4428	snow(i, nsrf) = 0.
4429	agesno(i, nsrf) = 0.
4430	ftsoil(i, :, nsrf) = tsurf(i, nsrf)
4431	#ifdef ISO
4432	xtsnow(:,i,nsrf) = 0.
4433	#endif
4434	ELSE
4435	pfois(nsrf) = pfois(nsrf) + 1
4436	ENDIF
4437	ENDDO
4438
4439	ENDDO
4440
4441	END SUBROUTINE pbl_surface_newfrac
4442
4443	!****************************************************************************************
4444
4445	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format