Context Navigation

pbl_surface_mod.F90 @ 2670

Last change on this file since 2670 was 2670, checked in by oboucher, 8 years ago
adding new diagnostics for evap/rain/snow weighted by surface types this parallels what was done for wbilo (net flux) new diagnostics are wevap_??? wrain_??? wsnow_??? where ??? is ter, oce, lic or sic.
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 2670 2016-10-16 09:17:32Z oboucher $
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	zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
852	zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
853	cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0.
854	kh(:)=0. ; kh_x(:)=0. ; kh_w(:)=0.
855	slab_wfbils(:)=0.
856	s_pblh(:)=0. ; s_pblh_x(:)=0. ; s_pblh_w(:)=0.
857	s_plcl(:)=0. ; s_plcl_x(:)=0. ; s_plcl_w(:)=0.
858	s_capCL(:)=0. ; s_oliqCL(:)=0. ; s_cteiCL(:)=0. ; s_pblT(:)=0.
859	s_therm(:)=0.
860	s_trmb1(:)=0. ; s_trmb2(:)=0. ; s_trmb3(:)=0.
861	zustar(:)=0.
862	zu10m(:)=0. ; zv10m(:)=0.
863	fder_print(:)=0.
864	zxqsurf(:)=0.
865	zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
866	solsw(:,:)=0. ; sollw(:,:)=0.
867	d_ts(:,:)=0.
868	evap(:,:)=0.
869	fluxlat(:,:)=0.
870	wfbils(:,:)=0. ; wfbilo(:,:)=0.
871	wfevap(:,:)=0. ; wfrain(:,:)=0. ; wfsnow(:,:)=0.
872	flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.
873	dflux_t(:)=0. ; dflux_q(:)=0.
874	zxsnow(:)=0.
875	zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.
876	qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
877	runoff(:)=0.
878	IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
879	zcoefh(:,:,:) = 0.0
880	zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used
881	zcoefm(:,:,:) = 0.0
882	zcoefm(:,1,:) = 999999. !
883	ELSE
884	zcoefm(:,:,is_ave)=0.
885	zcoefh(:,:,is_ave)=0.
886	ENDIF
887	!!
888	! The components "is_ave" of tke_x and wake_deltke are "OUT" variables
889	!jyg<
890	!! tke(:,:,is_ave)=0.
891	tke_x(:,:,is_ave)=0.
892	wake_dltke(:,:,is_ave)=0.
893	!>jyg
894	!!! jyg le 23/02/2013
895	t2m(:,:) = 999999. ! t2m and q2m are meaningfull only over sub-surfaces
896	q2m(:,:) = 999999. ! actually present in the grid cell.
897	!!!
898	rh2m(:) = 0. ; qsat2m(:) = 0.
899	!!!
900	!!! jyg le 10/02/2012
901	rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0.
902	!!!
903
904	! 2b) Initialization of all local variables that will be compressed later
905	!****************************************************************************************
906	!! cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0
907	ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0
908	!! zv1 = 0.0 ; yqsurf = 0.0
909	!albedo SB >>>
910	yqsurf = 0.0 ; yalb = 0.0 ; yalb_vis = 0.0
911	!albedo SB <<<
912	yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0
913	ysollw = 0.0 ; yz0m = 0.0 ; yz0h = 0.0 ; yu1 = 0.0
914	yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0
915	ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0
916	yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
917	yrugoro = 0.0 ; ywindsp = 0.0
918	!! d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0
919	yfluxlat=0.0
920	!! flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0
921	!! d_t_diss= 0.0 ;d_u = 0.0 ; d_v = 0.0
922	yqsol = 0.0
923	ytherm = 0.0 ; ytke=0.
924	! Martin
925	ysnowhgt = 0.0; yqsnow = 0.0 ; yrunoff = 0.0 ; ytoice =0.0
926	yalb3_new = 0.0 ; ysissnow = 0.0
927	ypphi = 0.0 ; ycldt = 0.0 ; yrmu0 = 0.0
928	! Martin
929
930	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
931	ytke_x=0. ; ytke_w=0. ; ywake_dltke=0.
932	y_d_t_x=0. ; y_d_t_w=0. ; y_d_q_x=0. ; y_d_q_w=0.
933	!! d_t_w=0. ; d_q_w=0.
934	!! d_t_x=0. ; d_q_x=0.
935	!! d_wake_dlt=0. ; d_wake_dlq=0.
936	yfluxlat_x=0. ; yfluxlat_w=0.
937	ywake_s=0. ; ywake_cstar=0. ;ywake_dens=0.
938	!!!
939	!!! nrlmd le 13/06/2011
940	tau_eq=0. ; delta_coef=0.
941	y_delta_flux_t1=0.
942	ydtsurf_th=0.
943	yts_x=0. ; yts_w=0.
944	y_delta_tsurf=0.
945	!!!
946	ytsoil = 999999.
947
948
949	! 2c) Initialization of all local variables computed within the subsurface loop and used later on
950	!****************************************************************************************
951	d_t_diss_x(:,:) = 0. ; d_t_diss_w(:,:) = 0.
952	d_u_x(:,:)=0. ; d_u_w(:,:)=0.
953	d_v_x(:,:)=0. ; d_v_w(:,:)=0.
954	flux_t_x(:,:,:)=0. ; flux_t_w(:,:,:)=0.
955	flux_q_x(:,:,:)=0. ; flux_q_w(:,:,:)=0.
956	!
957	!jyg<
958	flux_u_x(:,:,:)=0. ; flux_u_w(:,:,:)=0.
959	flux_v_x(:,:,:)=0. ; flux_v_w(:,:,:)=0.
960	fluxlat_x(:,:)=0. ; fluxlat_w(:,:)=0.
961	!>jyg
962	!
963	!jyg<
964	! pblh,plcl,capCL,cteiCL ... are meaningfull only over sub-surfaces
965	! actually present in the grid cell ==> value set to 999999.
966	!
967	!jyg<
968	ustar(:,:) = 999999.
969	wstar(:,:) = 999999.
970	windsp(:,:) = SQRT(u10m(:,:)2 + v10m(:,:)2 )
971	u10m(:,:) = 999999.
972	v10m(:,:) = 999999.
973	!>jyg
974	!
975	pblh(:,:) = 999999. ! Hauteur de couche limite
976	plcl(:,:) = 999999. ! Niveau de condensation de la CLA
977	capCL(:,:) = 999999. ! CAPE de couche limite
978	oliqCL(:,:) = 999999. ! eau_liqu integree de couche limite
979	cteiCL(:,:) = 999999. ! cloud top instab. crit. couche limite
980	pblt(:,:) = 999999. ! T a la Hauteur de couche limite
981	therm(:,:) = 999999.
982	trmb1(:,:) = 999999. ! deep_cape
983	trmb2(:,:) = 999999. ! inhibition
984	trmb3(:,:) = 999999. ! Point Omega
985	!
986	t2m_x(:,:) = 999999.
987	q2m_x(:,:) = 999999.
988	ustar_x(:,:) = 999999.
989	wstar_x(:,:) = 999999.
990	u10m_x(:,:) = 999999.
991	v10m_x(:,:) = 999999.
992	!
993	pblh_x(:,:) = 999999. ! Hauteur de couche limite
994	plcl_x(:,:) = 999999. ! Niveau de condensation de la CLA
995	capCL_x(:,:) = 999999. ! CAPE de couche limite
996	oliqCL_x(:,:) = 999999. ! eau_liqu integree de couche limite
997	cteiCL_x(:,:) = 999999. ! cloud top instab. crit. couche limite
998	pblt_x(:,:) = 999999. ! T a la Hauteur de couche limite
999	therm_x(:,:) = 999999.
1000	trmb1_x(:,:) = 999999. ! deep_cape
1001	trmb2_x(:,:) = 999999. ! inhibition
1002	trmb3_x(:,:) = 999999. ! Point Omega
1003	!
1004	t2m_w(:,:) = 999999.
1005	q2m_w(:,:) = 999999.
1006	ustar_w(:,:) = 999999.
1007	wstar_w(:,:) = 999999.
1008	u10m_w(:,:) = 999999.
1009	v10m_w(:,:) = 999999.
1010
1011	pblh_w(:,:) = 999999. ! Hauteur de couche limite
1012	plcl_w(:,:) = 999999. ! Niveau de condensation de la CLA
1013	capCL_w(:,:) = 999999. ! CAPE de couche limite
1014	oliqCL_w(:,:) = 999999. ! eau_liqu integree de couche limite
1015	cteiCL_w(:,:) = 999999. ! cloud top instab. crit. couche limite
1016	pblt_w(:,:) = 999999. ! T a la Hauteur de couche limite
1017	therm_w(:,:) = 999999.
1018	trmb1_w(:,:) = 999999. ! deep_cape
1019	trmb2_w(:,:) = 999999. ! inhibition
1020	trmb3_w(:,:) = 999999. ! Point Omega
1021	!!!
1022	!
1023	!!!
1024	!****************************************************************************************
1025	! 3) - Calculate pressure thickness of each layer
1026	! - Calculate the wind at first layer
1027	! - Mean calculations of albedo
1028	! - Calculate net radiance at sub-surface
1029	!****************************************************************************************
1030	DO k = 1, klev
1031	DO i = 1, klon
1032	delp(i,k) = paprs(i,k)-paprs(i,k+1)
1033	ENDDO
1034	ENDDO
1035
1036	!****************************************************************************************
1037	! Test for rugos........ from physiq.. A la fin plutot???
1038	!
1039	!****************************************************************************************
1040
1041	DO nsrf = 1, nbsrf
1042	DO i = 1, klon
1043	z0m(i,nsrf) = MAX(z0m(i,nsrf),z0min)
1044	z0h(i,nsrf) = MAX(z0h(i,nsrf),z0min)
1045	ENDDO
1046	ENDDO
1047
1048	! Mean calculations of albedo
1049	!
1050	! * alb : mean albedo for whole SW interval
1051	!
1052	! Mean albedo for grid point
1053	! * alb_m : mean albedo at whole SW interval
1054
1055	alb_dir_m(:,:) = 0.0
1056	alb_dif_m(:,:) = 0.0
1057	DO k = 1, nsw
1058	DO nsrf = 1, nbsrf
1059	DO i = 1, klon
1060	alb_dir_m(i,k) = alb_dir_m(i,k) + alb_dir(i,k,nsrf) * pctsrf(i,nsrf)
1061	alb_dif_m(i,k) = alb_dif_m(i,k) + alb_dif(i,k,nsrf) * pctsrf(i,nsrf)
1062	ENDDO
1063	ENDDO
1064	ENDDO
1065
1066	! We here suppose the fraction f1 of incoming radiance of visible radiance
1067	! as a fraction of all shortwave radiance
1068	f1 = 0.5
1069	! f1 = 1 ! put f1=1 to recreate old calculations
1070
1071	!f1 is already included with SFRWL values in each surf files
1072	alb=0.0
1073	DO k=1,nsw
1074	DO nsrf = 1, nbsrf
1075	DO i = 1, klon
1076	alb(i,nsrf) = alb(i,nsrf) + alb_dir(i,k,nsrf)*SFRWL(k)
1077	ENDDO
1078	ENDDO
1079	ENDDO
1080
1081	alb_m=0.0
1082	DO k = 1,nsw
1083	DO i = 1, klon
1084	alb_m(i) = alb_m(i) + alb_dir_m(i,k)*SFRWL(k)
1085	END DO
1086	ENDDO
1087	!albedo SB <<<
1088
1089
1090
1091	! Calculation of mean temperature at surface grid points
1092	ztsol(:) = 0.0
1093	DO nsrf = 1, nbsrf
1094	DO i = 1, klon
1095	ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
1096	ENDDO
1097	ENDDO
1098
1099	! Linear distrubution on sub-surface of long- and shortwave net radiance
1100	DO nsrf = 1, nbsrf
1101	DO i = 1, klon
1102	sollw(i,nsrf) = sollw_m(i) + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ts(i,nsrf))
1103
1104	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1105	! ! Martin
1106	! Apparently introduced for sisvat but not used
1107	! sollwd(i,nsrf)= sollwd_m(i)
1108	! ! Martin
1109	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1110
1111	solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
1112	ENDDO
1113	ENDDO
1114	!
1115	!<al1: second order corrections
1116	!- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)
1117	IF (iflag_order2_sollw == 1) THEN
1118	meansqT(:) = 0. ! as working buffer
1119	DO nsrf = 1, nbsrf
1120	DO i = 1, klon
1121	meansqT(i) = meansqT(i)+(ts(i,nsrf)-ztsol(i))*2 pctsrf(i,nsrf)
1122	END DO
1123	END DO
1124	DO nsrf = 1, nbsrf
1125	DO i = 1, klon
1126	sollw(i,nsrf) = sollw(i,nsrf) &
1127	+ 6.0RSIGMAztsol(i)*2 (meansqT(i)-(ztsol(i)-ts(i,nsrf))**2)
1128	ENDDO
1129	ENDDO
1130	ENDIF ! iflag_order2_sollw == 1
1131	!>al1
1132
1133	!****************************************************************************************
1134	! 4) Loop over different surfaces
1135	!
1136	! Only points containing a fraction of the sub surface will be treated.
1137	!
1138	!****************************************************************************************
1139
1140	loop_nbsrf: DO nsrf = 1, nbsrf
1141	IF (prt_level >=10) print *,' Loop nsrf ',nsrf
1142
1143	! Search for index(ni) and size(knon) of domaine to treat
1144	ni(:) = 0
1145	knon = 0
1146	DO i = 1, klon
1147	IF (pctsrf(i,nsrf) > 0.) THEN
1148	knon = knon + 1
1149	ni(knon) = i
1150	ENDIF
1151	ENDDO
1152
1153	!!! jyg le 19/08/2012
1154	! IF (knon <= 0) THEN
1155	! IF (prt_level >= 10) print *,' no grid point for nsrf= ',nsrf
1156	! cycle loop_nbsrf
1157	! ENDIF
1158	!!!
1159
1160	! write index, with IOIPSL
1161	IF (debugindex .AND. mpi_size==1) THEN
1162	tabindx(:)=0.
1163	DO i=1,knon
1164	tabindx(i)=REAL(i)
1165	END DO
1166	debugtab(:,:) = 0.
1167	ndexbg(:) = 0
1168	CALL gath2cpl(tabindx,debugtab,knon,ni)
1169	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,nbp_lon*nbp_lat, ndexbg)
1170	ENDIF
1171
1172	!****************************************************************************************
1173	! 5) Compress variables
1174	!
1175	!****************************************************************************************
1176
1177	DO j = 1, knon
1178	i = ni(j)
1179	ypct(j) = pctsrf(i,nsrf)
1180	yts(j) = ts(i,nsrf)
1181	ysnow(j) = snow(i,nsrf)
1182	yqsurf(j) = qsurf(i,nsrf)
1183	yalb(j) = alb(i,nsrf)
1184	!albedo SB >>>
1185	yalb_vis(j) = alb_dir(i,1,nsrf)
1186	if(nsw==6)then
1187	yalb_vis(j)=(alb_dir(i,1,nsrf)SFRWL(1)+alb_dir(i,2,nsrf)SFRWL(2) &
1188	+alb_dir(i,3,nsrf)*SFRWL(3))/(SFRWL(1)+SFRWL(2)+SFRWL(3))
1189	endif
1190	!albedo SB <<<
1191	yrain_f(j) = rain_f(i)
1192	ysnow_f(j) = snow_f(i)
1193	yagesno(j) = agesno(i,nsrf)
1194	yfder(j) = fder(i)
1195	ylwdown(j) = lwdown_m(i)
1196	ygustiness(j) = gustiness(i)
1197	ysolsw(j) = solsw(i,nsrf)
1198	ysollw(j) = sollw(i,nsrf)
1199	yz0m(j) = z0m(i,nsrf)
1200	yz0h(j) = z0h(i,nsrf)
1201	yrugoro(j) = rugoro(i)
1202	yu1(j) = u(i,1)
1203	yv1(j) = v(i,1)
1204	ypaprs(j,klev+1) = paprs(i,klev+1)
1205	!jyg<
1206	!! ywindsp(j) = SQRT(u10m(i,nsrf)2 + v10m(i,nsrf)2 )
1207	ywindsp(j) = windsp(i,nsrf)
1208	!>jyg
1209	! Martin
1210	yzsig(j) = zsig(i)
1211	ycldt(j) = cldt(i)
1212	yrmu0(j) = rmu0(i)
1213	! Martin
1214	!!! nrlmd le 13/06/2011
1215	y_delta_tsurf(j)=delta_tsurf(i,nsrf)
1216	!!!
1217	END DO
1218
1219	DO k = 1, klev
1220	DO j = 1, knon
1221	i = ni(j)
1222	ypaprs(j,k) = paprs(i,k)
1223	ypplay(j,k) = pplay(i,k)
1224	ydelp(j,k) = delp(i,k)
1225	ENDDO
1226	ENDDO
1227	!!! jyg le 07/02/2012 et le 10/04/2013
1228	DO k = 1, klev
1229	DO j = 1, knon
1230	i = ni(j)
1231	!jyg<
1232	!! ytke(j,k) = tke(i,k,nsrf)
1233	ytke(j,k) = tke_x(i,k,nsrf)
1234	!>jyg
1235	yu(j,k) = u(i,k)
1236	yv(j,k) = v(i,k)
1237	yt(j,k) = t(i,k)
1238	yq(j,k) = q(i,k)
1239	ENDDO
1240	ENDDO
1241	!
1242	IF (iflag_split .eq.1) THEN
1243	!!! nrlmd le 02/05/2011
1244	DO k = 1, klev
1245	DO j = 1, knon
1246	i = ni(j)
1247	yu_x(j,k) = u(i,k)
1248	yv_x(j,k) = v(i,k)
1249	yt_x(j,k) = t(i,k)-wake_s(i)*wake_dlt(i,k)
1250	yq_x(j,k) = q(i,k)-wake_s(i)*wake_dlq(i,k)
1251	yu_w(j,k) = u(i,k)
1252	yv_w(j,k) = v(i,k)
1253	yt_w(j,k) = t(i,k)+(1.-wake_s(i))*wake_dlt(i,k)
1254	yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k)
1255	!!!
1256	ENDDO
1257	ENDDO
1258	!!! nrlmd le 02/05/2011
1259	DO k = 1, klev+1
1260	DO j = 1, knon
1261	i = ni(j)
1262	!jyg<
1263	!! ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)*wake_dltke(i,k,nsrf)
1264	!! ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))*wake_dltke(i,k,nsrf)
1265	!! ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
1266	!! ytke(j,k) = tke(i,k,nsrf)
1267	!
1268	ytke_x(j,k) = tke_x(i,k,nsrf)
1269	ytke(j,k) = tke_x(i,k,nsrf)+wake_s(i)*wake_dltke(i,k,nsrf)
1270	ytke_w(j,k) = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf)
1271	ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
1272	!>jyg
1273	ENDDO
1274	ENDDO
1275	!!!
1276	!!! jyg le 07/02/2012
1277	DO j = 1, knon
1278	i = ni(j)
1279	ywake_s(j)=wake_s(i)
1280	ywake_cstar(j)=wake_cstar(i)
1281	ywake_dens(j)=wake_dens(i)
1282	ENDDO
1283	!!!
1284	!!! nrlmd le 13/06/2011
1285	DO j=1,knon
1286	yts_x(j)=yts(j)-ywake_s(j)*y_delta_tsurf(j)
1287	yts_w(j)=yts(j)+(1.-ywake_s(j))*y_delta_tsurf(j)
1288	ENDDO
1289	!!!
1290	ENDIF ! (iflag_split .eq.1)
1291	!!!
1292	DO k = 1, nsoilmx
1293	DO j = 1, knon
1294	i = ni(j)
1295	ytsoil(j,k) = ftsoil(i,k,nsrf)
1296	END DO
1297	END DO
1298
1299	! qsol(water height in soil) only for bucket continental model
1300	IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
1301	DO j = 1, knon
1302	i = ni(j)
1303	yqsol(j) = qsol(i)
1304	END DO
1305	ENDIF
1306
1307	!****************************************************************************************
1308	! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
1309	!
1310	!****************************************************************************************
1311
1312	!!! jyg le 07/02/2012
1313	IF (iflag_split .eq.0) THEN
1314	!!!
1315	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1316	! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1317	! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1318	! yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
1319	! yts, yqsurf, yrugos, &
1320	! ycdragm, ycdragh )
1321	! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1322	DO i = 1, knon
1323	! print*,'PBL ',i,RD
1324	! print*,'PBL ',yt(i,1),ypaprs(i,1),ypplay(i,1)
1325	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1326	* (ypaprs(i,1)-ypplay(i,1))
1327	speed(i) = SQRT(yu(i,1)2+yv(i,1)2)
1328	END DO
1329	CALL cdrag(knon, nsrf, &
1330	speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
1331	yts, yqsurf, yz0m, yz0h, &
1332	ycdragm, ycdragh, zri1, pref )
1333
1334	! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
1335	IF (ok_prescr_ust) then
1336	DO i = 1, knon
1337	print *,'ycdragm avant=',ycdragm(i)
1338	vent= sqrt(yu(i,1)yu(i,1)+yv(i,1)yv(i,1))
1339	! ycdragm(i) = ustust/(1.+(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1)))
1340	! ycdragm(i) = ustust/((1.+sqrt(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1))) &
1341	! sqrt(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1)))
1342	ycdragm(i) = ust*ust/(1.+vent)/vent
1343	! print *,'ycdragm ust yu yv apres=',ycdragm(i),ust,yu(i,1),yv(i,1)
1344	ENDDO
1345	ENDIF
1346	IF (prt_level >=10) print *,'clcdrag -> ycdragh ', ycdragh
1347	ELSE !(iflag_split .eq.0)
1348
1349	! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1350	! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1351	! yu_x(:,1), yv_x(:,1), yt_x(:,1), yq_x(:,1), &
1352	! yts_x, yqsurf, yrugos, &
1353	! ycdragm_x, ycdragh_x )
1354	! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1355	DO i = 1, knon
1356	zgeo1_x(i) = RD * yt_x(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1357	* (ypaprs(i,1)-ypplay(i,1))
1358	speed_x(i) = SQRT(yu_x(i,1)2+yv_x(i,1)2)
1359	END DO
1360	CALL cdrag(knon, nsrf, &
1361	speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
1362	yts_x, yqsurf, yz0m, yz0h, &
1363	ycdragm_x, ycdragh_x, zri1_x, pref_x )
1364
1365	! --- special Dice. JYG+MPL 25112013
1366	IF (ok_prescr_ust) then
1367	DO i = 1, knon
1368	! print *,'ycdragm_x avant=',ycdragm_x(i)
1369	vent= sqrt(yu_x(i,1)yu_x(i,1)+yv_x(i,1)yv_x(i,1))
1370	ycdragm_x(i) = ust*ust/(1.+vent)/vent
1371	! print *,'ycdragm_x ust yu yv apres=',ycdragm_x(i),ust,yu_x(i,1),yv_x(i,1)
1372	ENDDO
1373	ENDIF
1374	IF (prt_level >=10) print *,'clcdrag -> ycdragh_x ', ycdragh_x
1375	!
1376	! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1377	! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1378	! yu_w(:,1), yv_w(:,1), yt_w(:,1), yq_w(:,1), &
1379	! yts_w, yqsurf, yz0m, &
1380	! ycdragm_w, ycdragh_w )
1381	! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1382	DO i = 1, knon
1383	zgeo1_w(i) = RD * yt_w(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1384	* (ypaprs(i,1)-ypplay(i,1))
1385	speed_w(i) = SQRT(yu_w(i,1)2+yv_w(i,1)2)
1386	END DO
1387	CALL cdrag(knon, nsrf, &
1388	speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
1389	yts_w, yqsurf, yz0m, yz0h, &
1390	ycdragm_w, ycdragh_w, zri1_w, pref_w )
1391
1392	! --- special Dice. JYG+MPL 25112013 puis BOMEX
1393	IF (ok_prescr_ust) then
1394	DO i = 1, knon
1395	! print *,'ycdragm_w avant=',ycdragm_w(i)
1396	vent= sqrt(yu_w(i,1)yu_w(i,1)+yv_w(i,1)yv_w(i,1))
1397	ycdragm_w(i) = ust*ust/(1.+vent)/vent
1398	! print *,'ycdragm_w ust yu yv apres=',ycdragm_w(i),ust,yu_w(i,1),yv_w(i,1)
1399	ENDDO
1400	ENDIF
1401	IF (prt_level >=10) print *,'clcdrag -> ycdragh_w ', ycdragh_w
1402	!!!
1403	ENDIF ! (iflag_split .eq.0)
1404	!!!
1405
1406
1407	!****************************************************************************************
1408	! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm.
1409	!
1410	!****************************************************************************************
1411
1412	!!! jyg le 07/02/2012
1413	IF (iflag_split .eq.0) THEN
1414	!!!
1415	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1416	IF (prt_level >=10) THEN
1417	print *,' args coef_diff_turb: yu ', yu
1418	print *,' args coef_diff_turb: yv ', yv
1419	print *,' args coef_diff_turb: yq ', yq
1420	print *,' args coef_diff_turb: yt ', yt
1421	print *,' args coef_diff_turb: yts ', yts
1422	print *,' args coef_diff_turb: yz0m ', yz0m
1423	print *,' args coef_diff_turb: yqsurf ', yqsurf
1424	print *,' args coef_diff_turb: ycdragm ', ycdragm
1425	print *,' args coef_diff_turb: ycdragh ', ycdragh
1426	print *,' args coef_diff_turb: ytke ', ytke
1427	ENDIF
1428	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
1429	ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
1430	ycoefm, ycoefh, ytke)
1431	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1432	! In this case, coef_diff_turb is called for the Cd only
1433	DO k = 2, klev
1434	DO j = 1, knon
1435	i = ni(j)
1436	ycoefh(j,k) = zcoefh(i,k,nsrf)
1437	ycoefm(j,k) = zcoefm(i,k,nsrf)
1438	ENDDO
1439	ENDDO
1440	ENDIF
1441	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh
1442	!
1443	ELSE !(iflag_split .eq.0)
1444	IF (prt_level >=10) THEN
1445	print *,' args coef_diff_turb: yu_x ', yu_x
1446	print *,' args coef_diff_turb: yv_x ', yv_x
1447	print *,' args coef_diff_turb: yq_x ', yq_x
1448	print *,' args coef_diff_turb: yt_x ', yt_x
1449	print *,' args coef_diff_turb: yts_x ', yts_x
1450	print *,' args coef_diff_turb: yqsurf ', yqsurf
1451	print *,' args coef_diff_turb: ycdragm_x ', ycdragm_x
1452	print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x
1453	print *,' args coef_diff_turb: ytke_x ', ytke_x
1454	ENDIF
1455	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
1456	ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf, ycdragm_x, &
1457	ycoefm_x, ycoefh_x, ytke_x)
1458	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1459	! In this case, coef_diff_turb is called for the Cd only
1460	DO k = 2, klev
1461	DO j = 1, knon
1462	i = ni(j)
1463	ycoefh_x(j,k) = zcoefh(i,k,nsrf)
1464	ycoefm_x(j,k) = zcoefm(i,k,nsrf)
1465	ENDDO
1466	ENDDO
1467	ENDIF
1468	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x
1469	!
1470	IF (prt_level >=10) THEN
1471	print *,' args coef_diff_turb: yu_w ', yu_w
1472	print *,' args coef_diff_turb: yv_w ', yv_w
1473	print *,' args coef_diff_turb: yq_w ', yq_w
1474	print *,' args coef_diff_turb: yt_w ', yt_w
1475	print *,' args coef_diff_turb: yts_w ', yts_w
1476	print *,' args coef_diff_turb: yqsurf ', yqsurf
1477	print *,' args coef_diff_turb: ycdragm_w ', ycdragm_w
1478	print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w
1479	print *,' args coef_diff_turb: ytke_w ', ytke_w
1480	ENDIF
1481	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
1482	ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf, ycdragm_w, &
1483	ycoefm_w, ycoefh_w, ytke_w)
1484	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1485	! In this case, coef_diff_turb is called for the Cd only
1486	DO k = 2, klev
1487	DO j = 1, knon
1488	i = ni(j)
1489	ycoefh_w(j,k) = zcoefh(i,k,nsrf)
1490	ycoefm_w(j,k) = zcoefm(i,k,nsrf)
1491	ENDDO
1492	ENDDO
1493	ENDIF
1494	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w
1495	!
1496	!!!jyg le 10/04/2013
1497	!! En attendant de traiter le transport des traceurs dans les poches froides, formule
1498	!! arbitraire pour ycoefh et ycoefm
1499	DO k = 2,klev
1500	DO j = 1,knon
1501	ycoefh(j,k) = ycoefh_x(j,k) + ywake_s(j)*(ycoefh_w(j,k) - ycoefh_x(j,k))
1502	ycoefm(j,k) = ycoefm_x(j,k) + ywake_s(j)*(ycoefm_w(j,k) - ycoefm_x(j,k))
1503	ENDDO
1504	ENDDO
1505	!!!
1506	ENDIF ! (iflag_split .eq.0)
1507	!!!
1508
1509	!****************************************************************************************
1510	!
1511	! 8) "La descente" - "The downhill"
1512	!
1513	! climb_hq_down and climb_wind_down calculate the coefficients
1514	! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
1515	! Only the coefficients at surface for H and Q are returned.
1516	!
1517	!****************************************************************************************
1518
1519	! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
1520	!!! jyg le 07/02/2012
1521	IF (iflag_split .eq.0) THEN
1522	!!!
1523	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1524	CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
1525	ydelp, yt, yq, dtime, &
1526	!!! jyg le 09/05/2011
1527	CcoefH, CcoefQ, DcoefH, DcoefQ, &
1528	Kcoef_hq, gama_q, gama_h, &
1529	!!!
1530	AcoefH, AcoefQ, BcoefH, BcoefQ)
1531	ELSE !(iflag_split .eq.0)
1532	CALL climb_hq_down(knon, ycoefh_x, ypaprs, ypplay, &
1533	ydelp, yt_x, yq_x, dtime, &
1534	!!! nrlmd le 02/05/2011
1535	CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
1536	Kcoef_hq_x, gama_q_x, gama_h_x, &
1537	!!!
1538	AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x)
1539	!
1540	CALL climb_hq_down(knon, ycoefh_w, ypaprs, ypplay, &
1541	ydelp, yt_w, yq_w, dtime, &
1542	!!! nrlmd le 02/05/2011
1543	CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
1544	Kcoef_hq_w, gama_q_w, gama_h_w, &
1545	!!!
1546	AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w)
1547	!!!
1548	ENDIF ! (iflag_split .eq.0)
1549	!!!
1550
1551	! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
1552	!!! jyg le 07/02/2012
1553	IF (iflag_split .eq.0) THEN
1554	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1555	CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
1556	!!! jyg le 09/05/2011
1557	CcoefU, CcoefV, DcoefU, DcoefV, &
1558	Kcoef_m, alf_1, alf_2, &
1559	!!!
1560	AcoefU, AcoefV, BcoefU, BcoefV)
1561	ELSE ! (iflag_split .eq.0)
1562	CALL climb_wind_down(knon, dtime, ycoefm_x, ypplay, ypaprs, yt_x, ydelp, yu_x, yv_x, &
1563	!!! nrlmd le 02/05/2011
1564	CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
1565	Kcoef_m_x, alf_1_x, alf_2_x, &
1566	!!!
1567	AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x)
1568	!
1569	CALL climb_wind_down(knon, dtime, ycoefm_w, ypplay, ypaprs, yt_w, ydelp, yu_w, yv_w, &
1570	!!! nrlmd le 02/05/2011
1571	CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
1572	Kcoef_m_w, alf_1_w, alf_2_w, &
1573	!!!
1574	AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w)
1575	!!!
1576	ENDIF ! (iflag_split .eq.0)
1577	!!!
1578
1579	!****************************************************************************************
1580	! 9) Small calculations
1581	!
1582	!****************************************************************************************
1583
1584	! - Reference pressure is given the values at surface level
1585	ypsref(:) = ypaprs(:,1)
1586
1587	! - CO2 field on 2D grid to be sent to ORCHIDEE
1588	! Transform to compressed field
1589	IF (carbon_cycle_cpl) THEN
1590	DO i=1,knon
1591	r_co2_ppm(i) = co2_send(ni(i))
1592	END DO
1593	ELSE
1594	r_co2_ppm(:) = co2_ppm ! Constant field
1595	END IF
1596
1597	!!! nrlmd le 13/06/2011
1598	!----- 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
1599	! Kech_h_x(j) = ycdragh_x(j) * &
1600	! (1.0+SQRT(yu_x(j,1)2+yv_x(j,1)2)) * &
1601	! ypplay(j,1)/(RD*yt_x(j,1))
1602	! Kech_h_w(j) = ycdragh_w(j) * &
1603	! (1.0+SQRT(yu_w(j,1)2+yv_w(j,1)2)) * &
1604	! ypplay(j,1)/(RD*yt_w(j,1))
1605	! Kech_h(j) = (1.-ywake_s(j))Kech_h_x(j)+ywake_s(j)Kech_h_w(j)
1606	!
1607	! Kech_m_x(j) = ycdragm_x(j) * &
1608	! (1.0+SQRT(yu_x(j,1)2+yv_x(j,1)2)) * &
1609	! ypplay(j,1)/(RD*yt_x(j,1))
1610	! Kech_m_w(j) = ycdragm_w(j) * &
1611	! (1.0+SQRT(yu_w(j,1)2+yv_w(j,1)2)) * &
1612	! ypplay(j,1)/(RD*yt_w(j,1))
1613	! Kech_m(j) = (1.-ywake_s(j))Kech_m_x(j)+ywake_s(j)Kech_m_w(j)
1614	!!!
1615
1616	!!! nrlmd le 02/05/2011 -----------------------On raccorde les 2 colonnes dans la couche 1
1617	!----------------------------------------------------------------------------------------
1618	!!! jyg le 07/02/2012
1619	IF (iflag_split .eq.1) THEN
1620	!!!
1621	!!! jyg le 09/04/2013 ; passage aux nouvelles expressions en differences
1622
1623	DO j=1,knon
1624	!
1625	! Calcul des coefficients d echange
1626	mod_wind_x = 1.0+SQRT(yu_x(j,1)2+yv_x(j,1)2)
1627	mod_wind_w = 1.0+SQRT(yu_w(j,1)2+yv_w(j,1)2)
1628	rho1 = ypplay(j,1)/(RD*yt(j,1))
1629	Kech_h_x(j) = ycdragh_x(j) * mod_wind_x * rho1
1630	Kech_h_w(j) = ycdragh_w(j) * mod_wind_w * rho1
1631	Kech_m_x(j) = ycdragm_x(j) * mod_wind_x * rho1
1632	Kech_m_w(j) = ycdragm_w(j) * mod_wind_w * rho1
1633	!
1634	dd_Kh = Kech_h_w(j) - Kech_h_x(j)
1635	dd_Km = Kech_m_w(j) - Kech_m_x(j)
1636	IF (prt_level >=10) THEN
1637	print *,' mod_wind_x, mod_wind_w ', mod_wind_x, mod_wind_w
1638	print *,' rho1 ',rho1
1639	print *,' ycdragh_x(j),ycdragm_x(j) ',ycdragh_x(j),ycdragm_x(j)
1640	print *,' ycdragh_w(j),ycdragm_w(j) ',ycdragh_w(j),ycdragm_w(j)
1641	print *,' dd_Kh: ',dd_KH
1642	ENDIF
1643	!
1644	Kech_h(j) = Kech_h_x(j) + ywake_s(j)*dd_Kh
1645	Kech_m(j) = Kech_m_x(j) + ywake_s(j)*dd_Km
1646	!
1647	! Calcul des coefficients d echange corriges des retroactions
1648	Kech_H_xp(j) = Kech_h_x(j)/(1.-BcoefH_x(j)Kech_h_x(j)dtime)
1649	Kech_H_wp(j) = Kech_h_w(j)/(1.-BcoefH_w(j)Kech_h_w(j)dtime)
1650	Kech_Q_xp(j) = Kech_h_x(j)/(1.-BcoefQ_x(j)Kech_h_x(j)dtime)
1651	Kech_Q_wp(j) = Kech_h_w(j)/(1.-BcoefQ_w(j)Kech_h_w(j)dtime)
1652	Kech_U_xp(j) = Kech_m_x(j)/(1.-BcoefU_x(j)Kech_m_x(j)dtime)
1653	Kech_U_wp(j) = Kech_m_w(j)/(1.-BcoefU_w(j)Kech_m_w(j)dtime)
1654	Kech_V_xp(j) = Kech_m_x(j)/(1.-BcoefV_x(j)Kech_m_x(j)dtime)
1655	Kech_V_wp(j) = Kech_m_w(j)/(1.-BcoefV_w(j)Kech_m_w(j)dtime)
1656	!
1657	dd_KHp = Kech_H_wp(j) - Kech_H_xp(j)
1658	dd_KQp = Kech_Q_wp(j) - Kech_Q_xp(j)
1659	dd_KUp = Kech_U_wp(j) - Kech_U_xp(j)
1660	dd_KVp = Kech_V_wp(j) - Kech_V_xp(j)
1661	!
1662	Kech_Hp(j) = Kech_H_xp(j) + ywake_s(j)*dd_KHp
1663	Kech_Qp(j) = Kech_Q_xp(j) + ywake_s(j)*dd_KQp
1664	Kech_Up(j) = Kech_U_xp(j) + ywake_s(j)*dd_KUp
1665	Kech_Vp(j) = Kech_V_xp(j) + ywake_s(j)*dd_KVp
1666	!
1667	! Calcul des differences w-x
1668	dd_CM = ycdragm_w(j) - ycdragm_x(j)
1669	dd_CH = ycdragh_w(j) - ycdragh_x(j)
1670	dd_u = yu_w(j,1) - yu_x(j,1)
1671	dd_v = yv_w(j,1) - yv_x(j,1)
1672	dd_t = yt_w(j,1) - yt_x(j,1)
1673	dd_q = yq_w(j,1) - yq_x(j,1)
1674	dd_AH = AcoefH_w(j) - AcoefH_x(j)
1675	dd_AQ = AcoefQ_w(j) - AcoefQ_x(j)
1676	dd_AU = AcoefU_w(j) - AcoefU_x(j)
1677	dd_AV = AcoefV_w(j) - AcoefV_x(j)
1678	dd_BH = BcoefH_w(j) - BcoefH_x(j)
1679	dd_BQ = BcoefQ_w(j) - BcoefQ_x(j)
1680	dd_BU = BcoefU_w(j) - BcoefU_x(j)
1681	dd_BV = BcoefV_w(j) - BcoefV_x(j)
1682	!
1683	IF (prt_level >=10) THEN
1684	print *,'Variables pour la fusion : Kech_H_xp(j)' ,Kech_H_xp(j)
1685	print *,'Variables pour la fusion : Kech_H_wp(j)' ,Kech_H_wp(j)
1686	print *,'Variables pour la fusion : Kech_Hp(j)' ,Kech_Hp(j)
1687	print *,'Variables pour la fusion : Kech_h(j)' ,Kech_h(j)
1688	ENDIF
1689	!
1690	! Calcul des coef A, B \'equivalents dans la couche 1
1691	!
1692	AcoefH(j) = AcoefH_x(j) + ywake_s(j)(Kech_H_wp(j)/Kech_Hp(j))dd_AH
1693	AcoefQ(j) = AcoefQ_x(j) + ywake_s(j)(Kech_Q_wp(j)/Kech_Qp(j))dd_AQ
1694	AcoefU(j) = AcoefU_x(j) + ywake_s(j)(Kech_U_wp(j)/Kech_Up(j))dd_AU
1695	AcoefV(j) = AcoefV_x(j) + ywake_s(j)(Kech_V_wp(j)/Kech_Vp(j))dd_AV
1696	!
1697	BcoefH(j) = BcoefH_x(j) + ywake_s(j)BcoefH_x(j)(dd_Kh/Kech_h(j))*(1.+Kech_H_wp(j)/Kech_Hp(j)) &
1698	+ ywake_s(j)(Kech_H_wp(j)/Kech_Hp(j))(Kech_h_w(j)/Kech_h(j))*dd_BH
1699
1700	BcoefQ(j) = BcoefQ_x(j) + ywake_s(j)BcoefQ_x(j)(dd_Kh/Kech_h(j))*(1.+Kech_Q_wp(j)/Kech_Qp(j)) &
1701	+ ywake_s(j)(Kech_Q_wp(j)/Kech_Qp(j))(Kech_h_w(j)/Kech_h(j))*dd_BQ
1702
1703	BcoefU(j) = BcoefU_x(j) + ywake_s(j)BcoefU_x(j)(dd_Km/Kech_h(j))*(1.+Kech_U_wp(j)/Kech_Up(j)) &
1704	+ ywake_s(j)(Kech_U_wp(j)/Kech_Up(j))(Kech_m_w(j)/Kech_m(j))*dd_BU
1705
1706	BcoefV(j) = BcoefV_x(j) + ywake_s(j)BcoefV_x(j)(dd_Km/Kech_h(j))*(1.+Kech_V_wp(j)/Kech_Vp(j)) &
1707	+ ywake_s(j)(Kech_V_wp(j)/Kech_Vp(j))(Kech_m_w(j)/Kech_m(j))*dd_BV
1708
1709	!
1710	! Calcul des cdrag \'equivalents dans la couche
1711	!
1712	ycdragm(j) = ycdragm_x(j) + ywake_s(j)*dd_CM
1713	ycdragh(j) = ycdragh_x(j) + ywake_s(j)*dd_CH
1714	!
1715	! Calcul de T, q, u et v \'equivalents dans la couche 1
1716	yt(j,1) = yt_x(j,1) + ywake_s(j)(Kech_h_w(j)/Kech_h(j))dd_t
1717	yq(j,1) = yq_x(j,1) + ywake_s(j)(Kech_h_w(j)/Kech_h(j))dd_q
1718	yu(j,1) = yu_x(j,1) + ywake_s(j)(Kech_m_w(j)/Kech_m(j))dd_u
1719	yv(j,1) = yv_x(j,1) + ywake_s(j)(Kech_m_w(j)/Kech_m(j))dd_v
1720
1721
1722	ENDDO
1723	!!!
1724	ENDIF ! (iflag_split .eq.1)
1725	!!!
1726
1727	!****************************************************************************************
1728	!
1729	! Calulate t2m and q2m for the case of calculation at land grid points
1730	! t2m and q2m are needed as input to ORCHIDEE
1731	!
1732	!****************************************************************************************
1733	IF (nsrf == is_ter) THEN
1734
1735	DO i = 1, knon
1736	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1737	* (ypaprs(i,1)-ypplay(i,1))
1738	END DO
1739
1740	! Calculate the temperature et relative humidity at 2m and the wind at 10m
1741	CALL stdlevvar(klon, knon, is_ter, zxli, &
1742	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
1743	yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
1744	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
1745
1746	END IF
1747
1748	!****************************************************************************************
1749	!
1750	! 10) Switch according to current surface
1751	! It is necessary to start with the continental surfaces because the ocean
1752	! needs their run-off.
1753	!
1754	!****************************************************************************************
1755	SELECT CASE(nsrf)
1756
1757	CASE(is_ter)
1758	! print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)
1759	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
1760	rlon, rlat, yrmu0, &
1761	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
1762	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1763	AcoefH, AcoefQ, BcoefH, BcoefQ, &
1764	AcoefU, AcoefV, BcoefU, BcoefV, &
1765	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1766	ylwdown, yq2m, yt2m, &
1767	ysnow, yqsol, yagesno, ytsoil, &
1768	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
1769	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
1770	y_flux_u1, y_flux_v1 )
1771
1772	! Special DICE MPL 05082013 puis BOMEX
1773	IF (ok_prescr_ust) THEN
1774	do j=1,knon
1775	! ysnow(:)=0.
1776	! yqsol(:)=0.
1777	! yagesno(:)=50.
1778	! ytsoil(:,:)=300.
1779	! yz0_new(:)=0.001
1780	! yevap(:)=flat/RLVTT
1781	! yfluxlat(:)=-flat
1782	! yfluxsens(:)=-fsens
1783	! yqsurf(:)=0.
1784	! ytsurf_new(:)=tg
1785	! y_dflux_t(:)=0.
1786	! y_dflux_q(:)=0.
1787	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)
1788	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)
1789	enddo
1790	ENDIF
1791
1792
1793	CASE(is_lic)
1794	! Martin
1795	CALL surf_landice(itap, dtime, knon, ni, &
1796	rlon, rlat, debut, lafin, &
1797	yrmu0, ylwdown, yalb, ypphi(:,1), &
1798	ysolsw, ysollw, yts, ypplay(:,1), &
1799	ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1800	AcoefH, AcoefQ, BcoefH, BcoefQ, &
1801	AcoefU, AcoefV, BcoefU, BcoefV, &
1802	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1803	ysnow, yqsurf, yqsol, yagesno, &
1804	ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
1805	ytsurf_new, y_dflux_t, y_dflux_q, &
1806	yzsig, ycldt, &
1807	ysnowhgt, yqsnow, ytoice, ysissnow, &
1808	yalb3_new, yrunoff, &
1809	y_flux_u1, y_flux_v1)
1810
1811	!jyg<
1812	!! alb3_lic(:)=0.
1813	!>jyg
1814	DO j = 1, knon
1815	i = ni(j)
1816	alb3_lic(i) = yalb3_new(j)
1817	snowhgt(i) = ysnowhgt(j)
1818	qsnow(i) = yqsnow(j)
1819	to_ice(i) = ytoice(j)
1820	sissnow(i) = ysissnow(j)
1821	runoff(i) = yrunoff(j)
1822	END DO
1823	! Martin
1824	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
1825	IF (ok_prescr_ust) THEN
1826	DO j=1,knon
1827	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)
1828	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)
1829	ENDDO
1830	ENDIF
1831
1832	CASE(is_oce)
1833	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
1834	ywindsp, rmu0, yfder, yts, &
1835	itap, dtime, jour, knon, ni, &
1836	ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1837	AcoefH, AcoefQ, BcoefH, BcoefQ, &
1838	AcoefU, AcoefV, BcoefU, BcoefV, &
1839	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1840	ysnow, yqsurf, yagesno, &
1841	yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
1842	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
1843	y_flux_u1, y_flux_v1)
1844	IF (prt_level >=10) THEN
1845	print *,'arg de surf_ocean: ycdragh ',ycdragh
1846	print *,'arg de surf_ocean: ycdragm ',ycdragm
1847	print *,'arg de surf_ocean: yt ', yt
1848	print *,'arg de surf_ocean: yq ', yq
1849	print *,'arg de surf_ocean: yts ', yts
1850	print *,'arg de surf_ocean: AcoefH ',AcoefH
1851	print *,'arg de surf_ocean: AcoefQ ',AcoefQ
1852	print *,'arg de surf_ocean: BcoefH ',BcoefH
1853	print *,'arg de surf_ocean: BcoefQ ',BcoefQ
1854	print *,'arg de surf_ocean: yevap ',yevap
1855	print *,'arg de surf_ocean: yfluxsens ',yfluxsens
1856	print *,'arg de surf_ocean: yfluxlat ',yfluxlat
1857	print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new
1858	ENDIF
1859	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
1860	IF (ok_prescr_ust) THEN
1861	DO j=1,knon
1862	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)
1863	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)
1864	ENDDO
1865	ENDIF
1866
1867	CASE(is_sic)
1868	CALL surf_seaice( &
1869	!albedo SB >>>
1870	rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, &
1871	!albedo SB <<<
1872	itap, dtime, jour, knon, ni, &
1873	lafin, &
1874	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1875	AcoefH, AcoefQ, BcoefH, BcoefQ, &
1876	AcoefU, AcoefV, BcoefU, BcoefV, &
1877	ypsref, yu1, yv1, ygustiness, pctsrf, &
1878	ysnow, yqsurf, yqsol, yagesno, ytsoil, &
1879	!albedo SB >>>
1880	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
1881	!albedo SB <<<
1882	ytsurf_new, y_dflux_t, y_dflux_q, &
1883	y_flux_u1, y_flux_v1)
1884
1885	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
1886	IF (ok_prescr_ust) THEN
1887	DO j=1,knon
1888	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)
1889	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)
1890	ENDDO
1891	ENDIF
1892
1893	CASE DEFAULT
1894	WRITE(lunout,*) 'Surface index = ', nsrf
1895	abort_message = 'Surface index not valid'
1896	CALL abort_physic(modname,abort_message,1)
1897	END SELECT
1898
1899
1900	!****************************************************************************************
1901	! 11) - Calcul the increment of surface temperature
1902	!
1903	!****************************************************************************************
1904
1905	if (evap0>=0.) then
1906	yevap(:)=evap0
1907	yevap(:)=RLVTT*evap0
1908	endif
1909
1910
1911	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
1912
1913	!****************************************************************************************
1914	!
1915	! 12) "La remontee" - "The uphill"
1916	!
1917	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
1918	! for X=H, Q, U and V, for all vertical levels.
1919	!
1920	!****************************************************************************************
1921
1922	!!!
1923	!!! jyg le 10/04/2013
1924	!!!
1925	IF (ok_flux_surf) THEN
1926	IF (prt_level >=10) THEN
1927	PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
1928	ENDIF
1929	y_flux_t1(:) = fsens
1930	y_flux_q1(:) = flat/RLVTT
1931	yfluxlat(:) = flat
1932	!
1933	IF (iflag_split .eq.0) THEN
1934	do j=1,knon
1935	Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)2+yv(j,1)2)) * &
1936	ypplay(j,1)/(RD*yt(j,1))
1937	enddo
1938	ENDIF ! (iflag_split .eq.0)
1939
1940	DO j = 1, knon
1941	yt1_new=(1./RCPD)(AcoefH(j)+BcoefH(j)y_flux_t1(j)*dtime)
1942	ytsurf_new(j)=yt1_new-y_flux_t1(j)/(Kech_h(j)*RCPD)
1943	ENDDO
1944
1945	do j=1,knon
1946	y_d_ts(j) = ytsurf_new(j) - yts(j)
1947	enddo
1948
1949	ELSE ! (ok_flux_surf)
1950	do j=1,knon
1951	y_flux_t1(j) = yfluxsens(j)
1952	y_flux_q1(j) = -yevap(j)
1953	enddo
1954	ENDIF
1955
1956	IF (prt_level >=10) THEN
1957	DO j=1,knon
1958	print*,'y_flux_t1,yfluxlat,wakes' &
1959	& , y_flux_t1(j), yfluxlat(j), ywake_s(j)
1960	print*,'beta,ytsurf_new', ybeta(j), ytsurf_new(j)
1961	print*,'effusivity,facteur,cstar', effusivity, facteur,wake_cstar(j)
1962	ENDDO
1963	ENDIF
1964
1965	!!! jyg le 07/02/2012 puis le 10/04/2013
1966	IF (iflag_split .eq.1) THEN
1967	!!!
1968	DO j=1,knon
1969	y_delta_flux_t1(j) = ( Kech_H_wp(j)Kech_H_xp(j)(AcoefH_w(j)-AcoefH_x(j)) + &
1970	y_flux_t1(j)*(Kech_H_wp(j)-Kech_H_xp(j)) ) / Kech_Hp(j)
1971	y_delta_flux_q1(j) = ( Kech_Q_wp(j)Kech_Q_xp(j)(AcoefQ_w(j)-AcoefQ_x(j)) + &
1972	y_flux_q1(j)*(Kech_Q_wp(j)-Kech_Q_xp(j)) ) / Kech_Qp(j)
1973	y_delta_flux_u1(j) = ( Kech_U_wp(j)Kech_U_xp(j)(AcoefU_w(j)-AcoefU_x(j)) + &
1974	y_flux_u1(j)*(Kech_U_wp(j)-Kech_U_xp(j)) ) / Kech_Up(j)
1975	y_delta_flux_v1(j) = ( Kech_V_wp(j)Kech_V_xp(j)(AcoefV_w(j)-AcoefV_x(j)) + &
1976	y_flux_v1(j)*(Kech_V_wp(j)-Kech_V_xp(j)) ) / Kech_Vp(j)
1977	!
1978	y_flux_t1_x(j)=y_flux_t1(j) - ywake_s(j)*y_delta_flux_t1(j)
1979	y_flux_t1_w(j)=y_flux_t1(j) + (1.-ywake_s(j))*y_delta_flux_t1(j)
1980	y_flux_q1_x(j)=y_flux_q1(j) - ywake_s(j)*y_delta_flux_q1(j)
1981	y_flux_q1_w(j)=y_flux_q1(j) + (1.-ywake_s(j))*y_delta_flux_q1(j)
1982	y_flux_u1_x(j)=y_flux_u1(j) - ywake_s(j)*y_delta_flux_u1(j)
1983	y_flux_u1_w(j)=y_flux_u1(j) + (1.-ywake_s(j))*y_delta_flux_u1(j)
1984	y_flux_v1_x(j)=y_flux_v1(j) - ywake_s(j)*y_delta_flux_v1(j)
1985	y_flux_v1_w(j)=y_flux_v1(j) + (1.-ywake_s(j))*y_delta_flux_v1(j)
1986	!
1987	yfluxlat_x(j)=y_flux_q1_x(j)*RLVTT
1988	yfluxlat_w(j)=y_flux_q1_w(j)*RLVTT
1989
1990	ENDDO
1991	!
1992
1993	!!jyg!! A reprendre apres reflexion ===============================================
1994	!!jyg!!
1995	!!jyg!! DO j=1,knon
1996	!!jyg!!!!! nrlmd le 13/06/2011
1997	!!jyg!!
1998	!!jyg!!!----Diffusion dans le sol dans le cas continental seulement
1999	!!jyg!! IF (nsrf.eq.is_ter) THEN
2000	!!jyg!!!----Calcul du coefficient delta_coeff
2001	!!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)))
2002	!!jyg!!
2003	!!jyg!!! delta_coef(j)=dtime/(effusivity*sqrt(tau_eq(j)))
2004	!!jyg!! delta_coef(j)=facteur*sqrt(tau_eq(j))/effusivity
2005	!!jyg!!! delta_coef(j)=0.
2006	!!jyg!! ELSE
2007	!!jyg!! delta_coef(j)=0.
2008	!!jyg!! ENDIF
2009	!!jyg!!
2010	!!jyg!!!----Calcul de delta_tsurf
2011	!!jyg!! y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j)
2012	!!jyg!!
2013	!!jyg!!!----Si il n'y a pas des poches...
2014	!!jyg!! IF (wake_cstar(j).le.0.01) THEN
2015	!!jyg!! y_delta_tsurf(j)=0.
2016	!!jyg!! y_delta_flux_t1(j)=0.
2017	!!jyg!! ENDIF
2018	!!jyg!!
2019	!!jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle)
2020	!!jyg!!!!!!! jyg le 23/02/2012
2021	!!jyg!!!!!!!
2022	!!jyg!!!! ybeta(j)=y_flux_q1(j) / &
2023	!!jyg!!!! & (Kech_h(j)*(yq(j,1)-yqsatsurf(j)))
2024	!!jyg!!!!!! ybeta(j)=-1.*yevap(j) / &
2025	!!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)))
2026	!!jyg!!!!!!! fin jyg
2027	!!jyg!!!!!
2028	!!jyg!!
2029	!!jyg!! ENDDO
2030	!!jyg!!
2031	!!jyg!!!!! fin nrlmd le 13/06/2011
2032	!!jyg!!
2033	IF (prt_level >=10) THEN
2034	DO j = 1, knon
2035	print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
2036	print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
2037	! print*,'tsurf_x,tsurf_w,tsurf,t1', ytsurf_th_x(j), ytsurf_th_w(j), ytsurf_th(j), yt(j,1)
2038	print*,'tsurf_x,t1x,tsurf_w,t1w,tsurf,t1,t1_ancien', &
2039	& 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)
2040	print*,'qsatsurf,qsatsurf_x,qsatsurf_w', yqsatsurf(j), yqsatsurf_x(j), yqsatsurf_w(j)
2041	print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j)
2042	ENDDO
2043
2044	DO j=1,knon
2045	print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
2046	& , 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)
2047	print*,'beta,ytsurf_new,yqsatsurf', ybeta(j), ytsurf_new(j), yqsatsurf(j)
2048	print*,'effusivity,facteur,cstar', effusivity, facteur,wake_cstar(j)
2049	ENDDO
2050	ENDIF
2051
2052	!!! jyg le 07/02/2012
2053	ENDIF ! (iflag_split .eq.1)
2054	!!!
2055
2056	!!! jyg le 07/02/2012
2057	IF (iflag_split .eq.0) THEN
2058	!!!
2059	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2060	CALL climb_hq_up(knon, dtime, yt, yq, &
2061	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
2062	!!! jyg le 07/02/2012
2063	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2064	CcoefH, CcoefQ, DcoefH, DcoefQ, &
2065	Kcoef_hq, gama_q, gama_h, &
2066	!!!
2067	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
2068	ELSE !(iflag_split .eq.0)
2069	CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
2070	y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &
2071	!!! nrlmd le 02/05/2011
2072	AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &
2073	CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
2074	Kcoef_hq_x, gama_q_x, gama_h_x, &
2075	!!!
2076	y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:))
2077	!
2078	CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
2079	y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &
2080	!!! nrlmd le 02/05/2011
2081	AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &
2082	CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
2083	Kcoef_hq_w, gama_q_w, gama_h_w, &
2084	!!!
2085	y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:))
2086	!!!
2087	ENDIF ! (iflag_split .eq.0)
2088	!!!
2089
2090	!!! jyg le 07/02/2012
2091	IF (iflag_split .eq.0) THEN
2092	!!!
2093	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2094	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
2095	!!! jyg le 07/02/2012
2096	AcoefU, AcoefV, BcoefU, BcoefV, &
2097	CcoefU, CcoefV, DcoefU, DcoefV, &
2098	Kcoef_m, &
2099	!!!
2100	y_flux_u, y_flux_v, y_d_u, y_d_v)
2101	y_d_t_diss(:,:)=0.
2102	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2103	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2104	& ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
2105	& ,iflag_pbl,nsrf)
2106	ENDIF
2107	! print*,'yamada_c OK'
2108
2109	ELSE !(iflag_split .eq.0)
2110	CALL climb_wind_up(knon, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &
2111	!!! nrlmd le 02/05/2011
2112	AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &
2113	CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
2114	Kcoef_m_x, &
2115	!!!
2116	y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)
2117	!
2118	y_d_t_diss_x(:,:)=0.
2119	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2120	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2121	& ,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 &
2122	,ycoefq_x,y_d_t_diss_x,yustar_x &
2123	& ,iflag_pbl,nsrf)
2124	ENDIF
2125	! print*,'yamada_c OK'
2126
2127	CALL climb_wind_up(knon, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &
2128	!!! nrlmd le 02/05/2011
2129	AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &
2130	CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
2131	Kcoef_m_w, &
2132	!!!
2133	y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)
2134	!!!
2135	y_d_t_diss_w(:,:)=0.
2136	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2137	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2138	& ,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 &
2139	,ycoefq_w,y_d_t_diss_w,yustar_w &
2140	& ,iflag_pbl,nsrf)
2141	ENDIF
2142	! print*,'yamada_c OK'
2143	!
2144	IF (prt_level >=10) THEN
2145	print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &
2146	yfluxlat_x, yfluxlat_w
2147	ENDIF
2148	!
2149	ENDIF ! (iflag_split .eq.0)
2150	!!!
2151
2152	DO j = 1, knon
2153	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
2154	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
2155	ENDDO
2156
2157	!****************************************************************************************
2158	! 13) Transform variables for output format :
2159	! - Decompress
2160	! - Multiply with pourcentage of current surface
2161	! - Cumulate in global variable
2162	!
2163	!****************************************************************************************
2164
2165
2166	!!! jyg le 07/02/2012
2167	IF (iflag_split .eq.0) THEN
2168	!!!
2169	DO k = 1, klev
2170	DO j = 1, knon
2171	i = ni(j)
2172	y_d_t_diss(j,k) = y_d_t_diss(j,k) * ypct(j)
2173	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
2174	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
2175	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
2176	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
2177
2178	flux_t(i,k,nsrf) = y_flux_t(j,k)
2179	flux_q(i,k,nsrf) = y_flux_q(j,k)
2180	flux_u(i,k,nsrf) = y_flux_u(j,k)
2181	flux_v(i,k,nsrf) = y_flux_v(j,k)
2182
2183
2184	ENDDO
2185	ENDDO
2186
2187
2188	ELSE !(iflag_split .eq.0)
2189
2190	! Tendances hors poches
2191	DO k = 1, klev
2192	DO j = 1, knon
2193	i = ni(j)
2194	y_d_t_diss_x(j,k) = y_d_t_diss_x(j,k) * ypct(j)
2195	y_d_t_x(j,k) = y_d_t_x(j,k) * ypct(j)
2196	y_d_q_x(j,k) = y_d_q_x(j,k) * ypct(j)
2197	y_d_u_x(j,k) = y_d_u_x(j,k) * ypct(j)
2198	y_d_v_x(j,k) = y_d_v_x(j,k) * ypct(j)
2199
2200	flux_t_x(i,k,nsrf) = y_flux_t_x(j,k)
2201	flux_q_x(i,k,nsrf) = y_flux_q_x(j,k)
2202	flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
2203	flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)
2204	ENDDO
2205	ENDDO
2206
2207	! Tendances dans les poches
2208	DO k = 1, klev
2209	DO j = 1, knon
2210	i = ni(j)
2211	y_d_t_diss_w(j,k) = y_d_t_diss_w(j,k) * ypct(j)
2212	y_d_t_w(j,k) = y_d_t_w(j,k) * ypct(j)
2213	y_d_q_w(j,k) = y_d_q_w(j,k) * ypct(j)
2214	y_d_u_w(j,k) = y_d_u_w(j,k) * ypct(j)
2215	y_d_v_w(j,k) = y_d_v_w(j,k) * ypct(j)
2216
2217	flux_t_w(i,k,nsrf) = y_flux_t_w(j,k)
2218	flux_q_w(i,k,nsrf) = y_flux_q_w(j,k)
2219	flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
2220	flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)
2221	ENDDO
2222	ENDDO
2223
2224	! Flux, tendances et Tke moyenne dans la maille
2225	DO k = 1, klev
2226	DO j = 1, knon
2227	i = ni(j)
2228	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))
2229	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))
2230	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))
2231	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))
2232	ENDDO
2233	ENDDO
2234	DO j=1,knon
2235	yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j))
2236	ENDDO
2237	IF (prt_level >=10) THEN
2238	print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &
2239	nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf)
2240	ENDIF
2241
2242	DO k = 1, klev
2243	DO j = 1, knon
2244	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))
2245	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))
2246	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))
2247	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))
2248	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))
2249	ENDDO
2250	ENDDO
2251
2252	ENDIF ! (iflag_split .eq.0)
2253	!!!
2254
2255	! print*,'Dans pbl OK1'
2256
2257	!jyg<
2258	!! evap(:,nsrf) = - flux_q(:,1,nsrf)
2259	!>jyg
2260	DO j = 1, knon
2261	i = ni(j)
2262	evap(i,nsrf) = - flux_q(i,1,nsrf) !jyg
2263	d_ts(i,nsrf) = y_d_ts(j)
2264	!albedo SB >>>
2265	do k=1,nsw
2266	alb_dir(i,k,nsrf) = yalb_dir_new(j,k)
2267	alb_dif(i,k,nsrf) = yalb_dif_new(j,k)
2268	enddo
2269	!albedo SB <<<
2270	snow(i,nsrf) = ysnow(j)
2271	qsurf(i,nsrf) = yqsurf(j)
2272	z0m(i,nsrf) = yz0m(j)
2273	z0h(i,nsrf) = yz0h(j)
2274	fluxlat(i,nsrf) = yfluxlat(j)
2275	agesno(i,nsrf) = yagesno(j)
2276	cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
2277	cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
2278	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
2279	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
2280	END DO
2281
2282	! print*,'Dans pbl OK2'
2283
2284	!!! jyg le 07/02/2012
2285	IF (iflag_split .eq.1) THEN
2286	!!!
2287	!!! nrlmd le 02/05/2011
2288	DO j = 1, knon
2289	i = ni(j)
2290	fluxlat_x(i,nsrf) = yfluxlat_x(j)
2291	fluxlat_w(i,nsrf) = yfluxlat_w(j)
2292	!!!
2293	!!! nrlmd le 13/06/2011
2294	delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
2295	cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
2296	cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)
2297	cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j)
2298	cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j)
2299	kh(i) = kh(i) + Kech_h(j)*ypct(j)
2300	kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j)
2301	kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j)
2302	!!!
2303	END DO
2304	!!!
2305	ENDIF ! (iflag_split .eq.1)
2306	!!!
2307	!!! nrlmd le 02/05/2011
2308	!!jyg le 20/02/2011
2309	!! tke_x(:,:,nsrf)=0.
2310	!! tke_w(:,:,nsrf)=0.
2311	!!jyg le 20/02/2011
2312	!! DO k = 1, klev+1
2313	!! DO j = 1, knon
2314	!! i = ni(j)
2315	!! wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2316	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2317	!! ENDDO
2318	!! ENDDO
2319	!!jyg le 20/02/2011
2320	!! DO k = 1, klev+1
2321	!! DO j = 1, knon
2322	!! i = ni(j)
2323	!! tke(i,k,nsrf)=(1.-ywake_s(j))tke_x(i,k,nsrf)+ywake_s(j)tke_w(i,k,nsrf)
2324	!! ENDDO
2325	!! ENDDO
2326	!!!
2327	IF (iflag_split .eq.0) THEN
2328	DO k = 2, klev
2329	DO j = 1, knon
2330	i = ni(j)
2331	!jyg<
2332	!! tke(i,k,nsrf) = ytke(j,k)
2333	!! tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
2334	tke_x(i,k,nsrf) = ytke(j,k)
2335	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
2336	!>jyg
2337	END DO
2338	END DO
2339
2340	ELSE
2341	DO k = 2, klev
2342	DO j = 1, knon
2343	i = ni(j)
2344	wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2345	!jyg<
2346	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2347	!! tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
2348	tke_x(i,k,nsrf) = ytke_x(j,k)
2349	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)
2350	wake_dltke(i,k,is_ave) = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
2351
2352	!>jyg
2353	ENDDO
2354	ENDDO
2355	ENDIF ! (iflag_split .eq.0)
2356	!!!
2357	DO k = 2, klev
2358	DO j = 1, knon
2359	i = ni(j)
2360	zcoefh(i,k,nsrf) = ycoefh(j,k)
2361	zcoefm(i,k,nsrf) = ycoefm(j,k)
2362	zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
2363	zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
2364	END DO
2365	END DO
2366
2367	! print*,'Dans pbl OK3'
2368
2369	IF ( nsrf .EQ. is_ter ) THEN
2370	DO j = 1, knon
2371	i = ni(j)
2372	qsol(i) = yqsol(j)
2373	END DO
2374	END IF
2375
2376	!jyg<
2377	!! ftsoil(:,:,nsrf) = 0.
2378	!>jyg
2379	DO k = 1, nsoilmx
2380	DO j = 1, knon
2381	i = ni(j)
2382	ftsoil(i, k, nsrf) = ytsoil(j,k)
2383	END DO
2384	END DO
2385
2386	!!! jyg le 07/02/2012
2387	IF (iflag_split .eq.1) THEN
2388	!!!
2389	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
2390	DO k = 1, klev
2391	DO j = 1, knon
2392	i = ni(j)
2393	d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)
2394	d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)
2395	d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)
2396	d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)
2397	d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)
2398	!
2399	d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)
2400	d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)
2401	d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)
2402	d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
2403	d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
2404	!
2405	!! d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
2406	!! d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)
2407	END DO
2408	END DO
2409	!!!
2410	ENDIF ! (iflag_split .eq.1)
2411	!!!
2412
2413	DO k = 1, klev
2414	DO j = 1, knon
2415	i = ni(j)
2416	d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
2417	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
2418	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
2419	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
2420	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
2421	END DO
2422	END DO
2423
2424	! print*,'Dans pbl OK4'
2425
2426	IF (prt_level >=10) THEN
2427	print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
2428	d_t_w(:,1), d_t_x(:,1), d_t(:,1)
2429	ENDIF
2430
2431	!****************************************************************************************
2432	! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m
2433	! Call HBTM
2434	!
2435	!****************************************************************************************
2436	!!!
2437	!
2438	#undef T2m
2439	#define T2m
2440	#ifdef T2m
2441	! Calculations of diagnostic t,q at 2m and u, v at 10m
2442
2443	! print*,'Dans pbl OK41'
2444	! print*,'tair1,yt(:,1),y_d_t(:,1)'
2445	! print*, tair1,yt(:,1),y_d_t(:,1)
2446	!!! jyg le 07/02/2012
2447	IF (iflag_split .eq.0) THEN
2448	DO j=1, knon
2449	uzon(j) = yu(j,1) + y_d_u(j,1)
2450	vmer(j) = yv(j,1) + y_d_v(j,1)
2451	tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
2452	qair1(j) = yq(j,1) + y_d_q(j,1)
2453	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2454	* (ypaprs(j,1)-ypplay(j,1))
2455	tairsol(j) = yts(j) + y_d_ts(j)
2456	qairsol(j) = yqsurf(j)
2457	END DO
2458	ELSE ! (iflag_split .eq.0)
2459	DO j=1, knon
2460	uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
2461	vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
2462	tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
2463	qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
2464	zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2465	* (ypaprs(j,1)-ypplay(j,1))
2466	tairsol(j) = yts(j) + y_d_ts(j)
2467	tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
2468	qairsol(j) = yqsurf(j)
2469	END DO
2470	DO j=1, knon
2471	uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
2472	vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
2473	tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
2474	qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
2475	zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2476	* (ypaprs(j,1)-ypplay(j,1))
2477	tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
2478	qairsol(j) = yqsurf(j)
2479	END DO
2480	!!!
2481	ENDIF ! (iflag_split .eq.0)
2482	!!!
2483	DO j=1, knon
2484	i = ni(j)
2485	rugo1(j) = yz0m(j)
2486	IF(nsrf.EQ.is_oce) THEN
2487	rugo1(j) = z0m(i,nsrf)
2488	ENDIF
2489	psfce(j)=ypaprs(j,1)
2490	patm(j)=ypplay(j,1)
2491	END DO
2492
2493	! print*,'Dans pbl OK42A'
2494	! print*,'tair1,yt(:,1),y_d_t(:,1)'
2495	! print*, tair1,yt(:,1),y_d_t(:,1)
2496
2497	! Calculate the temperature et relative humidity at 2m and the wind at 10m
2498	!!! jyg le 07/02/2012
2499	IF (iflag_split .eq.0) THEN
2500	CALL stdlevvar(klon, knon, nsrf, zxli, &
2501	uzon, vmer, tair1, qair1, zgeo1, &
2502	tairsol, qairsol, rugo1, rugo1, psfce, patm, &
2503	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
2504	ELSE !(iflag_split .eq.0)
2505	CALL stdlevvar(klon, knon, nsrf, zxli, &
2506	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
2507	tairsol_x, qairsol, rugo1, rugo1, psfce, patm, &
2508	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
2509	CALL stdlevvar(klon, knon, nsrf, zxli, &
2510	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
2511	tairsol_w, qairsol, rugo1, rugo1, psfce, patm, &
2512	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
2513	!!!
2514	ENDIF ! (iflag_split .eq.0)
2515	!!!
2516	!!! jyg le 07/02/2012
2517	IF (iflag_split .eq.0) THEN
2518	DO j=1, knon
2519	i = ni(j)
2520	t2m(i,nsrf)=yt2m(j)
2521	q2m(i,nsrf)=yq2m(j)
2522	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2523	ustar(i,nsrf)=yustar(j)
2524	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
2525	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
2526	END DO
2527	ELSE !(iflag_split .eq.0)
2528	DO j=1, knon
2529	i = ni(j)
2530	t2m_x(i,nsrf)=yt2m_x(j)
2531	q2m_x(i,nsrf)=yq2m_x(j)
2532	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2533	ustar_x(i,nsrf)=yustar_x(j)
2534	u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
2535	v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
2536	END DO
2537	DO j=1, knon
2538	i = ni(j)
2539	t2m_w(i,nsrf)=yt2m_w(j)
2540	q2m_w(i,nsrf)=yq2m_w(j)
2541	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2542	ustar_w(i,nsrf)=yustar_w(j)
2543	u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
2544	v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
2545	!
2546	ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
2547	u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
2548	v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
2549	END DO
2550	!!!
2551	ENDIF ! (iflag_split .eq.0)
2552	!!!
2553
2554	! print*,'Dans pbl OK43'
2555	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
2556	!IM Ajoute dependance type surface
2557	IF (thermcep) THEN
2558	!!! jyg le 07/02/2012
2559	IF (iflag_split .eq.0) THEN
2560	DO j = 1, knon
2561	i=ni(j)
2562	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
2563	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
2564	zx_qs1 = MIN(0.5,zx_qs1)
2565	zcor1 = 1./(1.-RETV*zx_qs1)
2566	zx_qs1 = zx_qs1*zcor1
2567
2568	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
2569	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
2570	END DO
2571	ELSE ! (iflag_split .eq.0)
2572	DO j = 1, knon
2573	i=ni(j)
2574	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
2575	zx_qs1 = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)
2576	zx_qs1 = MIN(0.5,zx_qs1)
2577	zcor1 = 1./(1.-RETV*zx_qs1)
2578	zx_qs1 = zx_qs1*zcor1
2579
2580	rh2m_x(i) = rh2m_x(i) + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
2581	qsat2m_x(i) = qsat2m_x(i) + zx_qs1 * pctsrf(i,nsrf)
2582	END DO
2583	DO j = 1, knon
2584	i=ni(j)
2585	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
2586	zx_qs1 = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
2587	zx_qs1 = MIN(0.5,zx_qs1)
2588	zcor1 = 1./(1.-RETV*zx_qs1)
2589	zx_qs1 = zx_qs1*zcor1
2590
2591	rh2m_w(i) = rh2m_w(i) + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
2592	qsat2m_w(i) = qsat2m_w(i) + zx_qs1 * pctsrf(i,nsrf)
2593	END DO
2594	!!!
2595	ENDIF ! (iflag_split .eq.0)
2596	!!!
2597	END IF
2598	!
2599	IF (prt_level >=10) THEN
2600	print *, 'T2m, q2m, RH2m ', &
2601	t2m, q2m, rh2m
2602	ENDIF
2603
2604	! print*,'OK pbl 5'
2605	!
2606	!!! jyg le 07/02/2012
2607	IF (iflag_split .eq.0) THEN
2608	CALL hbtm(knon, ypaprs, ypplay, &
2609	yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
2610	y_flux_t,y_flux_q,yu,yv,yt,yq, &
2611	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
2612	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
2613	IF (prt_level >=10) THEN
2614	print *,' Arg. de HBTM: yt2m ',yt2m
2615	print *,' Arg. de HBTM: yt10m ',yt10m
2616	print *,' Arg. de HBTM: yq2m ',yq2m
2617	print *,' Arg. de HBTM: yq10m ',yq10m
2618	print *,' Arg. de HBTM: yustar ',yustar
2619	print *,' Arg. de HBTM: y_flux_t ',y_flux_t
2620	print *,' Arg. de HBTM: y_flux_q ',y_flux_q
2621	print *,' Arg. de HBTM: yu ',yu
2622	print *,' Arg. de HBTM: yv ',yv
2623	print *,' Arg. de HBTM: yt ',yt
2624	print *,' Arg. de HBTM: yq ',yq
2625	ENDIF
2626	ELSE ! (iflag_split .eq.0)
2627	CALL HBTM(knon, ypaprs, ypplay, &
2628	yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
2629	y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
2630	ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
2631	ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
2632	IF (prt_level >=10) THEN
2633	print *,' Arg. de HBTM: yt2m_x ',yt2m_x
2634	print *,' Arg. de HBTM: yt10m_x ',yt10m_x
2635	print *,' Arg. de HBTM: yq2m_x ',yq2m_x
2636	print *,' Arg. de HBTM: yq10m_x ',yq10m_x
2637	print *,' Arg. de HBTM: yustar_x ',yustar_x
2638	print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x
2639	print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x
2640	print *,' Arg. de HBTM: yu_x ',yu_x
2641	print *,' Arg. de HBTM: yv_x ',yv_x
2642	print *,' Arg. de HBTM: yt_x ',yt_x
2643	print *,' Arg. de HBTM: yq_x ',yq_x
2644	ENDIF
2645	CALL HBTM(knon, ypaprs, ypplay, &
2646	yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
2647	y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
2648	ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
2649	ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
2650	!!!
2651	ENDIF ! (iflag_split .eq.0)
2652	!!!
2653
2654	!!! jyg le 07/02/2012
2655	IF (iflag_split .eq.0) THEN
2656	!!!
2657	DO j=1, knon
2658	i = ni(j)
2659	pblh(i,nsrf) = ypblh(j)
2660	wstar(i,nsrf) = ywstar(j)
2661	plcl(i,nsrf) = ylcl(j)
2662	capCL(i,nsrf) = ycapCL(j)
2663	oliqCL(i,nsrf) = yoliqCL(j)
2664	cteiCL(i,nsrf) = ycteiCL(j)
2665	pblT(i,nsrf) = ypblT(j)
2666	therm(i,nsrf) = ytherm(j)
2667	trmb1(i,nsrf) = ytrmb1(j)
2668	trmb2(i,nsrf) = ytrmb2(j)
2669	trmb3(i,nsrf) = ytrmb3(j)
2670	END DO
2671	IF (prt_level >=10) THEN
2672	print *, 'After HBTM: pblh ', pblh
2673	print *, 'After HBTM: plcl ', plcl
2674	print *, 'After HBTM: cteiCL ', cteiCL
2675	ENDIF
2676	ELSE !(iflag_split .eq.0)
2677	DO j=1, knon
2678	i = ni(j)
2679	pblh_x(i,nsrf) = ypblh_x(j)
2680	wstar_x(i,nsrf) = ywstar_x(j)
2681	plcl_x(i,nsrf) = ylcl_x(j)
2682	capCL_x(i,nsrf) = ycapCL_x(j)
2683	oliqCL_x(i,nsrf) = yoliqCL_x(j)
2684	cteiCL_x(i,nsrf) = ycteiCL_x(j)
2685	pblT_x(i,nsrf) = ypblT_x(j)
2686	therm_x(i,nsrf) = ytherm_x(j)
2687	trmb1_x(i,nsrf) = ytrmb1_x(j)
2688	trmb2_x(i,nsrf) = ytrmb2_x(j)
2689	trmb3_x(i,nsrf) = ytrmb3_x(j)
2690	END DO
2691	IF (prt_level >=10) THEN
2692	print *, 'After HBTM: pblh_x ', pblh_x
2693	print *, 'After HBTM: plcl_x ', plcl_x
2694	print *, 'After HBTM: cteiCL_x ', cteiCL_x
2695	ENDIF
2696	DO j=1, knon
2697	i = ni(j)
2698	pblh_w(i,nsrf) = ypblh_w(j)
2699	wstar_w(i,nsrf) = ywstar_w(j)
2700	plcl_w(i,nsrf) = ylcl_w(j)
2701	capCL_w(i,nsrf) = ycapCL_w(j)
2702	oliqCL_w(i,nsrf) = yoliqCL_w(j)
2703	cteiCL_w(i,nsrf) = ycteiCL_w(j)
2704	pblT_w(i,nsrf) = ypblT_w(j)
2705	therm_w(i,nsrf) = ytherm_w(j)
2706	trmb1_w(i,nsrf) = ytrmb1_w(j)
2707	trmb2_w(i,nsrf) = ytrmb2_w(j)
2708	trmb3_w(i,nsrf) = ytrmb3_w(j)
2709	END DO
2710	IF (prt_level >=10) THEN
2711	print *, 'After HBTM: pblh_w ', pblh_w
2712	print *, 'After HBTM: plcl_w ', plcl_w
2713	print *, 'After HBTM: cteiCL_w ', cteiCL_w
2714	ENDIF
2715	!!!
2716	ENDIF ! (iflag_split .eq.0)
2717	!!!
2718
2719	! print*,'OK pbl 6'
2720	#else
2721	! T2m not defined
2722	! No calculation
2723	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
2724	#endif
2725
2726	!****************************************************************************************
2727	! 15) End of loop over different surfaces
2728	!
2729	!****************************************************************************************
2730	END DO loop_nbsrf
2731
2732	!****************************************************************************************
2733	! 16) Calculate the mean value over all sub-surfaces for some variables
2734	!
2735	!****************************************************************************************
2736
2737	z0m(:,nbsrf+1) = 0.0
2738	z0h(:,nbsrf+1) = 0.0
2739	DO nsrf = 1, nbsrf
2740	DO i = 1, klon
2741	z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)
2742	z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)
2743	ENDDO
2744	ENDDO
2745
2746	! print*,'OK pbl 7'
2747	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
2748	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
2749	zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
2750	zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
2751	zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
2752	zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
2753
2754	!!! jyg le 07/02/2012
2755	IF (iflag_split .eq.1) THEN
2756	!!!
2757	!!! nrlmd & jyg les 02/05/2011, 05/02/2012
2758
2759	DO nsrf = 1, nbsrf
2760	DO k = 1, klev
2761	DO i = 1, klon
2762	zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
2763	zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
2764	zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
2765	zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
2766	!
2767	zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
2768	zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
2769	zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
2770	zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
2771	END DO
2772	END DO
2773	END DO
2774
2775	DO i = 1, klon
2776	zxsens_x(i) = - zxfluxt_x(i,1)
2777	zxsens_w(i) = - zxfluxt_w(i,1)
2778	END DO
2779	!!!
2780	ENDIF ! (iflag_split .eq.1)
2781	!!!
2782
2783	DO nsrf = 1, nbsrf
2784	DO k = 1, klev
2785	DO i = 1, klon
2786	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
2787	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
2788	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
2789	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
2790	END DO
2791	END DO
2792	END DO
2793
2794	DO i = 1, klon
2795	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
2796	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
2797	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
2798	ENDDO
2799	!!!
2800
2801	!
2802	! Incrementer la temperature du sol
2803	!
2804	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
2805	zt2m(:) = 0.0 ; zq2m(:) = 0.0
2806	zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
2807	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
2808	!!! jyg le 07/02/2012
2809	s_pblh_x(:) = 0.0 ; s_plcl_x(:) = 0.0
2810	s_pblh_w(:) = 0.0 ; s_plcl_w(:) = 0.0
2811	!!!
2812	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
2813	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
2814	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
2815	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
2816	wstar(:,is_ave)=0.
2817
2818	! print*,'OK pbl 9'
2819
2820	!!! nrlmd le 02/05/2011
2821	zxfluxlat_x(:) = 0.0 ; zxfluxlat_w(:) = 0.0
2822	!!!
2823
2824	DO nsrf = 1, nbsrf
2825	DO i = 1, klon
2826	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
2827
2828	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
2829	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
2830	wfbilo(i,nsrf) = (evap(i,nsrf)-(rain_f(i)+snow_f(i)))*pctsrf(i,nsrf)
2831	wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)
2832	wfrain(i,nsrf) = rain_f(i)*pctsrf(i,nsrf)
2833	wfsnow(i,nsrf) = snow_f(i)*pctsrf(i,nsrf)
2834
2835	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
2836	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
2837	END DO
2838	END DO
2839	!
2840	!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
2841	IF (iflag_order2_sollw == 1) THEN
2842	meansqT(:) = 0. ! as working buffer
2843	DO nsrf = 1, nbsrf
2844	DO i = 1, klon
2845	meansqT(i) = meansqT(i)+(ts(i,nsrf)-zxtsol(i))*2 pctsrf(i,nsrf)
2846	END DO
2847	END DO
2848	zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)
2849	ENDIF ! iflag_order2_sollw == 1
2850	!>al1
2851
2852	!!! jyg le 07/02/2012
2853	IF (iflag_split .eq.0) THEN
2854	DO nsrf = 1, nbsrf
2855	DO i = 1, klon
2856	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
2857	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
2858	zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
2859	wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
2860	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
2861	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
2862
2863	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
2864	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
2865	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
2866	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
2867	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
2868	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
2869	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
2870	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
2871	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
2872	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
2873	END DO
2874	END DO
2875	ELSE !(iflag_split .eq.0)
2876	DO nsrf = 1, nbsrf
2877	DO i = 1, klon
2878	!!! nrlmd le 02/05/2011
2879	zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
2880	zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
2881	!!!
2882	!!! jyg le 08/02/2012
2883	!! Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;
2884	!! pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
2885	!! pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
2886	!! pour les autres variables, on sort les valeurs de la region (x).
2887	zt2m(i) = zt2m(i) + (t2m_x(i,nsrf)+wake_s(i)(t2m_w(i,nsrf)-t2m_x(i,nsrf))) pctsrf(i,nsrf)
2888	zq2m(i) = zq2m(i) + q2m_x(i,nsrf) * pctsrf(i,nsrf)
2889	zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
2890	wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
2891	zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
2892	zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
2893	!
2894	s_pblh(i) = s_pblh(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
2895	s_pblh_x(i) = s_pblh_x(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
2896	s_pblh_w(i) = s_pblh_w(i) + pblh_w(i,nsrf) * pctsrf(i,nsrf)
2897	!
2898	s_plcl(i) = s_plcl(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
2899	s_plcl_x(i) = s_plcl_x(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
2900	s_plcl_w(i) = s_plcl_w(i) + plcl_w(i,nsrf) * pctsrf(i,nsrf)
2901	!
2902	s_capCL(i) = s_capCL(i) + capCL_x(i,nsrf) * pctsrf(i,nsrf)
2903	s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
2904	s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
2905	s_pblT(i) = s_pblT(i) + pblT_x(i,nsrf) * pctsrf(i,nsrf)
2906	s_therm(i) = s_therm(i) + therm_x(i,nsrf) * pctsrf(i,nsrf)
2907	s_trmb1(i) = s_trmb1(i) + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
2908	s_trmb2(i) = s_trmb2(i) + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
2909	s_trmb3(i) = s_trmb3(i) + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
2910	END DO
2911	END DO
2912	DO i = 1, klon
2913	qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
2914	END DO
2915	!!!
2916	ENDIF ! (iflag_split .eq.0)
2917	!!!
2918
2919	IF (check) THEN
2920	amn=MIN(ts(1,is_ter),1000.)
2921	amx=MAX(ts(1,is_ter),-1000.)
2922	DO i=2, klon
2923	amn=MIN(ts(i,is_ter),amn)
2924	amx=MAX(ts(i,is_ter),amx)
2925	ENDDO
2926	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
2927	ENDIF
2928
2929	!jg ?
2930	!!$!
2931	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
2932	!!$! sub-surfaces is distributed.
2933	!!$!
2934	!!$ DO nsrf = 1, nbsrf
2935	!!$ DO i = 1, klon
2936	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
2937	!!$ ts(i,nsrf) = zxtsol(i)
2938	!!$ t2m(i,nsrf) = zt2m(i)
2939	!!$ q2m(i,nsrf) = zq2m(i)
2940	!!$ u10m(i,nsrf) = zu10m(i)
2941	!!$ v10m(i,nsrf) = zv10m(i)
2942	!!$
2943	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
2944	!!$ pblh(i,nsrf) = s_pblh(i)
2945	!!$ plcl(i,nsrf) = s_plcl(i)
2946	!!$ capCL(i,nsrf) = s_capCL(i)
2947	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
2948	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
2949	!!$ pblT(i,nsrf) = s_pblT(i)
2950	!!$ therm(i,nsrf) = s_therm(i)
2951	!!$ trmb1(i,nsrf) = s_trmb1(i)
2952	!!$ trmb2(i,nsrf) = s_trmb2(i)
2953	!!$ trmb3(i,nsrf) = s_trmb3(i)
2954	!!$ ENDIF
2955	!!$ ENDDO
2956	!!$ ENDDO
2957
2958
2959	DO i = 1, klon
2960	fder(i) = - 4.0RSIGMAzxtsol(i)**3
2961	ENDDO
2962
2963	zxqsurf(:) = 0.0
2964	zxsnow(:) = 0.0
2965	DO nsrf = 1, nbsrf
2966	DO i = 1, klon
2967	zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
2968	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
2969	END DO
2970	END DO
2971
2972	! Premier niveau de vent sortie dans physiq.F
2973	zu1(:) = u(:,1)
2974	zv1(:) = v(:,1)
2975
2976
2977	END SUBROUTINE pbl_surface
2978	!
2979	!****************************************************************************************
2980	!
2981	SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
2982
2983	USE indice_sol_mod
2984
2985	INCLUDE "dimsoil.h"
2986
2987	! Ouput variables
2988	!****************************************************************************************
2989	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
2990	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
2991	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
2992	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
2993
2994
2995	!****************************************************************************************
2996	! Return module variables for writing to restart file
2997	!
2998	!****************************************************************************************
2999	fder_rst(:) = fder(:)
3000	snow_rst(:,:) = snow(:,:)
3001	qsurf_rst(:,:) = qsurf(:,:)
3002	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
3003
3004	!****************************************************************************************
3005	! Deallocate module variables
3006	!
3007	!****************************************************************************************
3008	! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
3009	IF (ALLOCATED(fder)) DEALLOCATE(fder)
3010	IF (ALLOCATED(snow)) DEALLOCATE(snow)
3011	IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
3012	IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
3013
3014	END SUBROUTINE pbl_surface_final
3015	!
3016	!****************************************************************************************
3017	!
3018
3019	!albedo SB >>>
3020	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
3021	evap, z0m, z0h, agesno, &
3022	tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)
3023	!albedo SB <<<
3024	! Give default values where new fraction has appread
3025
3026	USE indice_sol_mod
3027
3028	INCLUDE "dimsoil.h"
3029	INCLUDE "clesphys.h"
3030	INCLUDE "compbl.h"
3031
3032	! Input variables
3033	!****************************************************************************************
3034	INTEGER, INTENT(IN) :: itime
3035	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
3036
3037	! InOutput variables
3038	!****************************************************************************************
3039	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
3040	!albedo SB >>>
3041	REAL, DIMENSION(klon,nsw,nbsrf), INTENT(INOUT) :: alb_dir, alb_dif
3042	INTEGER :: k
3043	!albedo SB <<<
3044	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m
3045	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: evap, agesno
3046	REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT) :: z0m,z0h
3047	REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
3048
3049	! Local variables
3050	!****************************************************************************************
3051	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
3052	CHARACTER(len=80) :: abort_message
3053	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
3054	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
3055	!
3056	! All at once !!
3057	!****************************************************************************************
3058
3059	DO nsrf = 1, nbsrf
3060	! First decide complement sub-surfaces
3061	SELECT CASE (nsrf)
3062	CASE(is_oce)
3063	nsrf_comp1=is_sic
3064	nsrf_comp2=is_ter
3065	nsrf_comp3=is_lic
3066	CASE(is_sic)
3067	nsrf_comp1=is_oce
3068	nsrf_comp2=is_ter
3069	nsrf_comp3=is_lic
3070	CASE(is_ter)
3071	nsrf_comp1=is_lic
3072	nsrf_comp2=is_oce
3073	nsrf_comp3=is_sic
3074	CASE(is_lic)
3075	nsrf_comp1=is_ter
3076	nsrf_comp2=is_oce
3077	nsrf_comp3=is_sic
3078	END SELECT
3079
3080	! Initialize all new fractions
3081	DO i=1, klon
3082	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
3083
3084	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
3085	! Use the complement sub-surface, keeping the continents unchanged
3086	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
3087	evap(i,nsrf) = evap(i,nsrf_comp1)
3088	z0m(i,nsrf) = z0m(i,nsrf_comp1)
3089	z0h(i,nsrf) = z0h(i,nsrf_comp1)
3090	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
3091	!albedo SB >>>
3092	DO k=1,nsw
3093	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp1)
3094	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp1)
3095	ENDDO
3096	!albedo SB <<<
3097	ustar(i,nsrf) = ustar(i,nsrf_comp1)
3098	u10m(i,nsrf) = u10m(i,nsrf_comp1)
3099	v10m(i,nsrf) = v10m(i,nsrf_comp1)
3100	if (iflag_pbl > 1) then
3101	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
3102	endif
3103	mfois(nsrf) = mfois(nsrf) + 1
3104	! F. Codron sensible default values for ocean and sea ice
3105	IF (nsrf.EQ.is_oce) THEN
3106	tsurf(i,nsrf) = 271.35 ! Freezing sea water
3107	DO k=1,nsw
3108	alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
3109	alb_dif(i,k,nsrf) = 0.06
3110	END DO
3111	ELSE IF (nsrf.EQ.is_sic) THEN
3112	tsurf(i,nsrf) = 273.15 ! Melting ice
3113	DO k=1,nsw
3114	alb_dir(i,k,nsrf) = 0.3 ! thin ice
3115	alb_dif(i,k,nsrf) = 0.3
3116	END DO
3117	END IF
3118	! F. Codron
3119	ELSE
3120	! The continents have changed. The new fraction receives the mean sum of the existent fractions
3121	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3122	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3123	z0m(i,nsrf) = z0m(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0m(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3124	z0h(i,nsrf) = z0h(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0h(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3125	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3126	!albedo SB >>>
3127	DO k=1,nsw
3128	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3129	alb_dir(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3130	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3131	alb_dif(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3132	ENDDO
3133	!albedo SB <<<
3134	ustar(i,nsrf) = ustar(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3135	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3136	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3137	if (iflag_pbl > 1) then
3138	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3139	endif
3140
3141	! Security abort. This option has never been tested. To test, comment the following line.
3142	! abort_message='The fraction of the continents have changed!'
3143	! CALL abort_physic(modname,abort_message,1)
3144	nfois(nsrf) = nfois(nsrf) + 1
3145	END IF
3146	snow(i,nsrf) = 0.
3147	agesno(i,nsrf) = 0.
3148	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
3149	ELSE
3150	pfois(nsrf) = pfois(nsrf)+ 1
3151	END IF
3152	END DO
3153
3154	END DO
3155
3156	END SUBROUTINE pbl_surface_newfrac
3157
3158	!
3159	!****************************************************************************************
3160	!
3161
3162	END MODULE pbl_surface_mod
3163

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

Download in other formats:

Original Format