Context Navigation

pbl_surface_mod.F90 @ 3953

Last change on this file since 3953 was 3953, checked in by jyg, 3 years ago
Bug fix: output arguments are added to subroutine wx_pbl_prelim0 so that exchange coefficients Kh, Kh_w and Kh_x are no longer undefined.
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: 147.0 KB

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

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

Download in other formats:

Original Format