Context Navigation

pbl_surface_mod.F90 @ 2853

Last change on this file since 2853 was 2410, checked in by jghattas, 9 years ago

Added the variable yrmu0(cosine of the solar zenith angle) into ORCHIDEE:

Current version of LMDZ can now be used with ORCHIDEE trunk revisions from rev 2961 and newer. For revision 1078-2960 on ORCHIDEE trunk, a small modification (change coszang into sinang) in surf_land_orchidee_mod.f90 is needed. For older versions than ORCHIDEE trunk revision 1078, the interface in surf_land_orchidee_noopenmp_mod.f90 should be used (add cpp key ORCHIDEE_NOOPENMP).

For details see ticket https://forge.ipsl.jussieu.fr/orchidee/ticket/217

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 125.0 KB

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

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

Download in other formats:

Original Format