Context Navigation

pbl_surface_mod.F90 @ 5408

Last change on this file since 5408 was 3886, checked in by Laurent Fairhead, 4 years ago
By request, re-introduction of two 'bugs' so that running without the iflag_new_t2mq2m gives exactly the same results as CMIP6 previous simulations
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 132.3 KB

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

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

Download in other formats:

Original Format