Context Navigation

pbl_surface_mod.F90 @ 4370

Last change on this file since 4370 was 4370, checked in by lguez, 17 months ago
Merge branch Ocean_skin back into trunk
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: 148.1 KB

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

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

Download in other formats:

Original Format