Context Navigation

pbl_surface_mod.F90 @ 3940

Last change on this file since 3940 was 3906, checked in by jyg, 3 years ago

1/ Bug fix in wx_pbl_mod.F90

2/ New pbl_split option:
iflag_pbl_split=3 ==> iflag_split=1 over ocean,

iflag_split=0 everywhere else.

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: 146.8 KB

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

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

Download in other formats:

Original Format