Context Navigation

pbl_surface_mod.F90 @ 5183

Last change on this file since 5183 was 2936, checked in by fcheruy, 7 years ago
corrected restart for drgtree Cd_frein is read from physiq.def
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:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 129.8 KB

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

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

Download in other formats:

Original Format