Context Navigation

pbl_surface_mod.F90 @ 2824

Last change on this file since 2824 was 2680, checked in by fhourdin, 8 years ago
Rustine sur ocean_albedo pour le 1D Correction sur la version conservative de Mellor et Yamada
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: 127.3 KB

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

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

Download in other formats:

Original Format