Context Navigation

pbl_surface_mod.F90 @ 3421

Last change on this file since 3421 was 3402, checked in by lguez, 6 years ago
Set temperature of appearing sea ice to -1.8 Celsius degrees instead of 0 Celsius degrees, following advice of F. Codron.
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: 132.0 KB

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

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

Download in other formats:

Original Format