Context Navigation

pbl_surface_mod.F90 @ 5191

Last change on this file since 5191 was 2886, checked in by Laurent Fairhead, 7 years ago
Merged trunk changes r2865:2885 into testing branch
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.4 KB

Line
1	!
2	! $Id: pbl_surface_mod.F90 2886 2017-05-20 07:41:16Z evignon $
3	!
4	MODULE pbl_surface_mod
5	!
6	! Planetary Boundary Layer and Surface module
7	!
8	! This module manages the calculation of turbulent diffusion in the boundary layer
9	! and all interactions towards the differents sub-surfaces.
10	!
11	!
12	USE dimphy
13	USE mod_phys_lmdz_para, ONLY : mpi_size
14	USE mod_grid_phy_lmdz, ONLY : klon_glo
15	USE ioipsl
16	USE surface_data, ONLY : type_ocean, ok_veget
17	USE surf_land_mod, ONLY : surf_land
18	USE surf_landice_mod, ONLY : surf_landice
19	USE surf_ocean_mod, ONLY : surf_ocean
20	USE surf_seaice_mod, ONLY : surf_seaice
21	USE cpl_mod, ONLY : gath2cpl
22	USE climb_hq_mod, ONLY : climb_hq_down, climb_hq_up
23	USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up
24	USE coef_diff_turb_mod, ONLY : coef_diff_turb
25
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	!!jyg iflag_split = mod(iflag_pbl_split,2)
781	iflag_split = mod(iflag_pbl_split,10)
782
783	!****************************************************************************************
784	! 1) Initialisation and validation tests
785	! Only done first time entering this subroutine
786	!
787	!****************************************************************************************
788
789	IF (first_call) THEN
790	print*,'PBL SURFACE AVEC GUSTINESS'
791	first_call=.FALSE.
792
793	! Initialize ok_flux_surf (for 1D model)
794	if (klon_glo>1) ok_flux_surf=.FALSE.
795
796	! Initilize debug IO
797	IF (debugindex .AND. mpi_size==1) THEN
798	! initialize IOIPSL output
799	idayref = day_ini
800	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
801	CALL gr_fi_ecrit(1,klon,nbp_lon,nbp_lat,rlon,zx_lon)
802	DO i = 1, nbp_lon
803	zx_lon(i,1) = rlon(i+1)
804	zx_lon(i,nbp_lat) = rlon(i+1)
805	ENDDO
806	CALL gr_fi_ecrit(1,klon,nbp_lon,nbp_lat,rlat,zx_lat)
807	CALL histbeg("sous_index",nbp_lon,zx_lon(:,1),nbp_lat,zx_lat(1,:), &
808	1,nbp_lon,1,nbp_lat, &
809	itau_phy,zjulian,dtime,nhoridbg,nidbg)
810	! no vertical axis
811	cl_surf(1)='ter'
812	cl_surf(2)='lic'
813	cl_surf(3)='oce'
814	cl_surf(4)='sic'
815	DO nsrf=1,nbsrf
816	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",nbp_lon, &
817	nbp_lat,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
818	END DO
819
820	CALL histend(nidbg)
821	CALL histsync(nidbg)
822
823	END IF
824
825	ENDIF
826
827	!****************************************************************************************
828	! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
829	! instead of ORCHIDEE)
830	IF (qsol0>=0.) THEN
831	PRINT*,'WARNING : On impose qsol=',qsol0
832	qsol(:)=qsol0
833	ENDIF
834	!****************************************************************************************
835
836	!****************************************************************************************
837	! 2) Initialization to zero
838	!****************************************************************************************
839	!
840	! 2a) Initialization of all argument variables with INTENT(OUT)
841	!****************************************************************************************
842	cdragh(:)=0. ; cdragm(:)=0.
843	zu1(:)=0. ; zv1(:)=0.
844	!albedo SB >>>
845	alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0.
846	!albedo SB <<<
847	zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0.
848	d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0.
849	zxfluxlat(:)=0.
850	zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
851	d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=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 .ge.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 .ge.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 ! (prt_level >=10)
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)
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)
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)
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 .ge.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 .ge.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	wake_dltke(:,:,nsrf) = 0.
2329	DO k = 1, klev
2330	DO j = 1, knon
2331	i = ni(j)
2332	!jyg<
2333	!! tke(i,k,nsrf) = ytke(j,k)
2334	!! tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
2335	tke_x(i,k,nsrf) = ytke(j,k)
2336	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
2337	!>jyg
2338	END DO
2339	END DO
2340
2341	ELSE ! (iflag_split .eq.0)
2342	DO k = 1, klev
2343	DO j = 1, knon
2344	i = ni(j)
2345	wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2346	!jyg<
2347	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2348	!! tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
2349	tke_x(i,k,nsrf) = ytke_x(j,k)
2350	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)
2351	wake_dltke(i,k,is_ave) = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
2352
2353	!>jyg
2354	ENDDO
2355	ENDDO
2356	ENDIF ! (iflag_split .eq.0)
2357	!!!
2358	DO k = 2, klev
2359	DO j = 1, knon
2360	i = ni(j)
2361	zcoefh(i,k,nsrf) = ycoefh(j,k)
2362	zcoefm(i,k,nsrf) = ycoefm(j,k)
2363	zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
2364	zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
2365	END DO
2366	END DO
2367
2368	! print*,'Dans pbl OK3'
2369
2370	IF ( nsrf .EQ. is_ter ) THEN
2371	DO j = 1, knon
2372	i = ni(j)
2373	qsol(i) = yqsol(j)
2374	END DO
2375	END IF
2376
2377	!jyg<
2378	!! ftsoil(:,:,nsrf) = 0.
2379	!>jyg
2380	DO k = 1, nsoilmx
2381	DO j = 1, knon
2382	i = ni(j)
2383	ftsoil(i, k, nsrf) = ytsoil(j,k)
2384	END DO
2385	END DO
2386
2387	!!! jyg le 07/02/2012
2388	IF (iflag_split .ge.1) THEN
2389	!!!
2390	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
2391	DO k = 1, klev
2392	DO j = 1, knon
2393	i = ni(j)
2394	d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)
2395	d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)
2396	d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)
2397	d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)
2398	d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)
2399	!
2400	d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)
2401	d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)
2402	d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)
2403	d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
2404	d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
2405	!
2406	!! d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
2407	!! d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)
2408	END DO
2409	END DO
2410	!!!
2411	ENDIF ! (iflag_split .ge.1)
2412	!!!
2413
2414	DO k = 1, klev
2415	DO j = 1, knon
2416	i = ni(j)
2417	d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
2418	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
2419	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
2420	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
2421	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
2422	END DO
2423	END DO
2424
2425	! print*,'Dans pbl OK4'
2426
2427	IF (prt_level >=10) THEN
2428	print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
2429	d_t_w(:,1), d_t_x(:,1), d_t(:,1)
2430	ENDIF
2431
2432	!****************************************************************************************
2433	! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m
2434	! Call HBTM
2435	!
2436	!****************************************************************************************
2437	!!!
2438	!
2439	#undef T2m
2440	#define T2m
2441	#ifdef T2m
2442	! Calculations of diagnostic t,q at 2m and u, v at 10m
2443
2444	! print*,'Dans pbl OK41'
2445	! print*,'tair1,yt(:,1),y_d_t(:,1)'
2446	! print*, tair1,yt(:,1),y_d_t(:,1)
2447	!!! jyg le 07/02/2012
2448	IF (iflag_split .eq.0) THEN
2449	DO j=1, knon
2450	uzon(j) = yu(j,1) + y_d_u(j,1)
2451	vmer(j) = yv(j,1) + y_d_v(j,1)
2452	tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
2453	qair1(j) = yq(j,1) + y_d_q(j,1)
2454	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2455	* (ypaprs(j,1)-ypplay(j,1))
2456	tairsol(j) = yts(j) + y_d_ts(j)
2457	qairsol(j) = yqsurf(j)
2458	END DO
2459	ELSE ! (iflag_split .eq.0)
2460	DO j=1, knon
2461	uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
2462	vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
2463	tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
2464	qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
2465	zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2466	* (ypaprs(j,1)-ypplay(j,1))
2467	tairsol(j) = yts(j) + y_d_ts(j)
2468	tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
2469	qairsol(j) = yqsurf(j)
2470	END DO
2471	DO j=1, knon
2472	uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
2473	vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
2474	tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
2475	qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
2476	zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2477	* (ypaprs(j,1)-ypplay(j,1))
2478	tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
2479	qairsol(j) = yqsurf(j)
2480	END DO
2481	!!!
2482	ENDIF ! (iflag_split .eq.0)
2483	!!!
2484	DO j=1, knon
2485	i = ni(j)
2486	rugo1(j) = yz0m(j)
2487	IF(nsrf.EQ.is_oce) THEN
2488	rugo1(j) = z0m(i,nsrf)
2489	ENDIF
2490	psfce(j)=ypaprs(j,1)
2491	patm(j)=ypplay(j,1)
2492	END DO
2493
2494	! print*,'Dans pbl OK42A'
2495	! print*,'tair1,yt(:,1),y_d_t(:,1)'
2496	! print*, tair1,yt(:,1),y_d_t(:,1)
2497
2498	! Calculate the temperature et relative humidity at 2m and the wind at 10m
2499	!!! jyg le 07/02/2012
2500	IF (iflag_split .eq.0) THEN
2501	CALL stdlevvar(klon, knon, nsrf, zxli, &
2502	uzon, vmer, tair1, qair1, zgeo1, &
2503	tairsol, qairsol, rugo1, rugo1, psfce, patm, &
2504	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
2505	ELSE !(iflag_split .eq.0)
2506	CALL stdlevvar(klon, knon, nsrf, zxli, &
2507	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
2508	tairsol_x, qairsol, rugo1, rugo1, psfce, patm, &
2509	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
2510	CALL stdlevvar(klon, knon, nsrf, zxli, &
2511	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
2512	tairsol_w, qairsol, rugo1, rugo1, psfce, patm, &
2513	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
2514	!!!
2515	ENDIF ! (iflag_split .eq.0)
2516	!!!
2517	!!! jyg le 07/02/2012
2518	IF (iflag_split .eq.0) THEN
2519	DO j=1, knon
2520	i = ni(j)
2521	t2m(i,nsrf)=yt2m(j)
2522	q2m(i,nsrf)=yq2m(j)
2523	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2524	ustar(i,nsrf)=yustar(j)
2525	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
2526	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
2527	END DO
2528	ELSE !(iflag_split .eq.0)
2529	DO j=1, knon
2530	i = ni(j)
2531	t2m_x(i,nsrf)=yt2m_x(j)
2532	q2m_x(i,nsrf)=yq2m_x(j)
2533	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2534	ustar_x(i,nsrf)=yustar_x(j)
2535	u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
2536	v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
2537	END DO
2538	DO j=1, knon
2539	i = ni(j)
2540	t2m_w(i,nsrf)=yt2m_w(j)
2541	q2m_w(i,nsrf)=yq2m_w(j)
2542	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2543	ustar_w(i,nsrf)=yustar_w(j)
2544	u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
2545	v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
2546	!
2547	ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
2548	u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
2549	v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
2550	END DO
2551	!!!
2552	ENDIF ! (iflag_split .eq.0)
2553	!!!
2554
2555	! print*,'Dans pbl OK43'
2556	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
2557	!IM Ajoute dependance type surface
2558	IF (thermcep) THEN
2559	!!! jyg le 07/02/2012
2560	IF (iflag_split .eq.0) THEN
2561	DO j = 1, knon
2562	i=ni(j)
2563	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
2564	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
2565	zx_qs1 = MIN(0.5,zx_qs1)
2566	zcor1 = 1./(1.-RETV*zx_qs1)
2567	zx_qs1 = zx_qs1*zcor1
2568
2569	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
2570	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
2571	END DO
2572	ELSE ! (iflag_split .eq.0)
2573	DO j = 1, knon
2574	i=ni(j)
2575	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
2576	zx_qs1 = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)
2577	zx_qs1 = MIN(0.5,zx_qs1)
2578	zcor1 = 1./(1.-RETV*zx_qs1)
2579	zx_qs1 = zx_qs1*zcor1
2580
2581	rh2m_x(i) = rh2m_x(i) + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
2582	qsat2m_x(i) = qsat2m_x(i) + zx_qs1 * pctsrf(i,nsrf)
2583	END DO
2584	DO j = 1, knon
2585	i=ni(j)
2586	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
2587	zx_qs1 = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
2588	zx_qs1 = MIN(0.5,zx_qs1)
2589	zcor1 = 1./(1.-RETV*zx_qs1)
2590	zx_qs1 = zx_qs1*zcor1
2591
2592	rh2m_w(i) = rh2m_w(i) + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
2593	qsat2m_w(i) = qsat2m_w(i) + zx_qs1 * pctsrf(i,nsrf)
2594	END DO
2595	!!!
2596	ENDIF ! (iflag_split .eq.0)
2597	!!!
2598	END IF
2599	!
2600	IF (prt_level >=10) THEN
2601	print *, 'T2m, q2m, RH2m ', &
2602	t2m, q2m, rh2m
2603	ENDIF
2604
2605	! print*,'OK pbl 5'
2606	!
2607	!!! jyg le 07/02/2012
2608	IF (iflag_split .eq.0) THEN
2609	CALL hbtm(knon, ypaprs, ypplay, &
2610	yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
2611	y_flux_t,y_flux_q,yu,yv,yt,yq, &
2612	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
2613	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
2614	IF (prt_level >=10) THEN
2615	print *,' Arg. de HBTM: yt2m ',yt2m
2616	print *,' Arg. de HBTM: yt10m ',yt10m
2617	print *,' Arg. de HBTM: yq2m ',yq2m
2618	print *,' Arg. de HBTM: yq10m ',yq10m
2619	print *,' Arg. de HBTM: yustar ',yustar
2620	print *,' Arg. de HBTM: y_flux_t ',y_flux_t
2621	print *,' Arg. de HBTM: y_flux_q ',y_flux_q
2622	print *,' Arg. de HBTM: yu ',yu
2623	print *,' Arg. de HBTM: yv ',yv
2624	print *,' Arg. de HBTM: yt ',yt
2625	print *,' Arg. de HBTM: yq ',yq
2626	ENDIF
2627	ELSE ! (iflag_split .eq.0)
2628	CALL HBTM(knon, ypaprs, ypplay, &
2629	yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
2630	y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
2631	ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
2632	ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
2633	IF (prt_level >=10) THEN
2634	print *,' Arg. de HBTM: yt2m_x ',yt2m_x
2635	print *,' Arg. de HBTM: yt10m_x ',yt10m_x
2636	print *,' Arg. de HBTM: yq2m_x ',yq2m_x
2637	print *,' Arg. de HBTM: yq10m_x ',yq10m_x
2638	print *,' Arg. de HBTM: yustar_x ',yustar_x
2639	print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x
2640	print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x
2641	print *,' Arg. de HBTM: yu_x ',yu_x
2642	print *,' Arg. de HBTM: yv_x ',yv_x
2643	print *,' Arg. de HBTM: yt_x ',yt_x
2644	print *,' Arg. de HBTM: yq_x ',yq_x
2645	ENDIF
2646	CALL HBTM(knon, ypaprs, ypplay, &
2647	yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
2648	y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
2649	ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
2650	ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
2651	!!!
2652	ENDIF ! (iflag_split .eq.0)
2653	!!!
2654
2655	!!! jyg le 07/02/2012
2656	IF (iflag_split .eq.0) THEN
2657	!!!
2658	DO j=1, knon
2659	i = ni(j)
2660	pblh(i,nsrf) = ypblh(j)
2661	wstar(i,nsrf) = ywstar(j)
2662	plcl(i,nsrf) = ylcl(j)
2663	capCL(i,nsrf) = ycapCL(j)
2664	oliqCL(i,nsrf) = yoliqCL(j)
2665	cteiCL(i,nsrf) = ycteiCL(j)
2666	pblT(i,nsrf) = ypblT(j)
2667	therm(i,nsrf) = ytherm(j)
2668	trmb1(i,nsrf) = ytrmb1(j)
2669	trmb2(i,nsrf) = ytrmb2(j)
2670	trmb3(i,nsrf) = ytrmb3(j)
2671	END DO
2672	IF (prt_level >=10) THEN
2673	print *, 'After HBTM: pblh ', pblh
2674	print *, 'After HBTM: plcl ', plcl
2675	print *, 'After HBTM: cteiCL ', cteiCL
2676	ENDIF
2677	ELSE !(iflag_split .eq.0)
2678	DO j=1, knon
2679	i = ni(j)
2680	pblh_x(i,nsrf) = ypblh_x(j)
2681	wstar_x(i,nsrf) = ywstar_x(j)
2682	plcl_x(i,nsrf) = ylcl_x(j)
2683	capCL_x(i,nsrf) = ycapCL_x(j)
2684	oliqCL_x(i,nsrf) = yoliqCL_x(j)
2685	cteiCL_x(i,nsrf) = ycteiCL_x(j)
2686	pblT_x(i,nsrf) = ypblT_x(j)
2687	therm_x(i,nsrf) = ytherm_x(j)
2688	trmb1_x(i,nsrf) = ytrmb1_x(j)
2689	trmb2_x(i,nsrf) = ytrmb2_x(j)
2690	trmb3_x(i,nsrf) = ytrmb3_x(j)
2691	END DO
2692	IF (prt_level >=10) THEN
2693	print *, 'After HBTM: pblh_x ', pblh_x
2694	print *, 'After HBTM: plcl_x ', plcl_x
2695	print *, 'After HBTM: cteiCL_x ', cteiCL_x
2696	ENDIF
2697	DO j=1, knon
2698	i = ni(j)
2699	pblh_w(i,nsrf) = ypblh_w(j)
2700	wstar_w(i,nsrf) = ywstar_w(j)
2701	plcl_w(i,nsrf) = ylcl_w(j)
2702	capCL_w(i,nsrf) = ycapCL_w(j)
2703	oliqCL_w(i,nsrf) = yoliqCL_w(j)
2704	cteiCL_w(i,nsrf) = ycteiCL_w(j)
2705	pblT_w(i,nsrf) = ypblT_w(j)
2706	therm_w(i,nsrf) = ytherm_w(j)
2707	trmb1_w(i,nsrf) = ytrmb1_w(j)
2708	trmb2_w(i,nsrf) = ytrmb2_w(j)
2709	trmb3_w(i,nsrf) = ytrmb3_w(j)
2710	END DO
2711	IF (prt_level >=10) THEN
2712	print *, 'After HBTM: pblh_w ', pblh_w
2713	print *, 'After HBTM: plcl_w ', plcl_w
2714	print *, 'After HBTM: cteiCL_w ', cteiCL_w
2715	ENDIF
2716	!!!
2717	ENDIF ! (iflag_split .eq.0)
2718	!!!
2719
2720	! print*,'OK pbl 6'
2721	#else
2722	! T2m not defined
2723	! No calculation
2724	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
2725	#endif
2726
2727	!****************************************************************************************
2728	! 15) End of loop over different surfaces
2729	!
2730	!****************************************************************************************
2731	END DO loop_nbsrf
2732
2733	!****************************************************************************************
2734	! 16) Calculate the mean value over all sub-surfaces for some variables
2735	!
2736	!****************************************************************************************
2737
2738	z0m(:,nbsrf+1) = 0.0
2739	z0h(:,nbsrf+1) = 0.0
2740	DO nsrf = 1, nbsrf
2741	DO i = 1, klon
2742	z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)
2743	z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)
2744	ENDDO
2745	ENDDO
2746
2747	! print*,'OK pbl 7'
2748	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
2749	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
2750	zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
2751	zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
2752	zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
2753	zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
2754
2755	!!! jyg le 07/02/2012
2756	IF (iflag_split .ge.1) THEN
2757	!!!
2758	!!! nrlmd & jyg les 02/05/2011, 05/02/2012
2759
2760	DO nsrf = 1, nbsrf
2761	DO k = 1, klev
2762	DO i = 1, klon
2763	zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
2764	zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
2765	zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
2766	zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
2767	!
2768	zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
2769	zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
2770	zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
2771	zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
2772	END DO
2773	END DO
2774	END DO
2775
2776	DO i = 1, klon
2777	zxsens_x(i) = - zxfluxt_x(i,1)
2778	zxsens_w(i) = - zxfluxt_w(i,1)
2779	END DO
2780	!!!
2781	ENDIF ! (iflag_split .ge.1)
2782	!!!
2783
2784	DO nsrf = 1, nbsrf
2785	DO k = 1, klev
2786	DO i = 1, klon
2787	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
2788	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
2789	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
2790	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
2791	END DO
2792	END DO
2793	END DO
2794
2795	DO i = 1, klon
2796	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
2797	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
2798	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
2799	ENDDO
2800	!!!
2801
2802	!
2803	! Incrementer la temperature du sol
2804	!
2805	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
2806	zt2m(:) = 0.0 ; zq2m(:) = 0.0
2807	zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
2808	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
2809	!!! jyg le 07/02/2012
2810	s_pblh_x(:) = 0.0 ; s_plcl_x(:) = 0.0
2811	s_pblh_w(:) = 0.0 ; s_plcl_w(:) = 0.0
2812	!!!
2813	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
2814	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
2815	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
2816	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
2817	wstar(:,is_ave)=0.
2818
2819	! print*,'OK pbl 9'
2820
2821	!!! nrlmd le 02/05/2011
2822	zxfluxlat_x(:) = 0.0 ; zxfluxlat_w(:) = 0.0
2823	!!!
2824
2825	DO nsrf = 1, nbsrf
2826	DO i = 1, klon
2827	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
2828
2829	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
2830	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
2831	wfbilo(i,nsrf) = (evap(i,nsrf)-(rain_f(i)+snow_f(i)))*pctsrf(i,nsrf)
2832	wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)
2833	wfrain(i,nsrf) = rain_f(i)*pctsrf(i,nsrf)
2834	wfsnow(i,nsrf) = snow_f(i)*pctsrf(i,nsrf)
2835
2836	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
2837	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
2838	END DO
2839	END DO
2840	!
2841	!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
2842	IF (iflag_order2_sollw == 1) THEN
2843	meansqT(:) = 0. ! as working buffer
2844	DO nsrf = 1, nbsrf
2845	DO i = 1, klon
2846	meansqT(i) = meansqT(i)+(ts(i,nsrf)-zxtsol(i))*2 pctsrf(i,nsrf)
2847	END DO
2848	END DO
2849	zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)
2850	ENDIF ! iflag_order2_sollw == 1
2851	!>al1
2852
2853	!!! jyg le 07/02/2012
2854	IF (iflag_split .eq.0) THEN
2855	DO nsrf = 1, nbsrf
2856	DO i = 1, klon
2857	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
2858	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
2859	zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
2860	wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
2861	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
2862	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
2863
2864	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
2865	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
2866	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
2867	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
2868	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
2869	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
2870	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
2871	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
2872	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
2873	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
2874	END DO
2875	END DO
2876	ELSE !(iflag_split .eq.0)
2877	DO nsrf = 1, nbsrf
2878	DO i = 1, klon
2879	!!! nrlmd le 02/05/2011
2880	zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
2881	zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
2882	!!!
2883	!!! jyg le 08/02/2012
2884	!! Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;
2885	!! pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
2886	!! pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
2887	!! pour les autres variables, on sort les valeurs de la region (x).
2888	zt2m(i) = zt2m(i) + (t2m_x(i,nsrf)+wake_s(i)(t2m_w(i,nsrf)-t2m_x(i,nsrf))) pctsrf(i,nsrf)
2889	zq2m(i) = zq2m(i) + q2m_x(i,nsrf) * pctsrf(i,nsrf)
2890	zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
2891	wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
2892	zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
2893	zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
2894	!
2895	s_pblh(i) = s_pblh(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
2896	s_pblh_x(i) = s_pblh_x(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
2897	s_pblh_w(i) = s_pblh_w(i) + pblh_w(i,nsrf) * pctsrf(i,nsrf)
2898	!
2899	s_plcl(i) = s_plcl(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
2900	s_plcl_x(i) = s_plcl_x(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
2901	s_plcl_w(i) = s_plcl_w(i) + plcl_w(i,nsrf) * pctsrf(i,nsrf)
2902	!
2903	s_capCL(i) = s_capCL(i) + capCL_x(i,nsrf) * pctsrf(i,nsrf)
2904	s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
2905	s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
2906	s_pblT(i) = s_pblT(i) + pblT_x(i,nsrf) * pctsrf(i,nsrf)
2907	s_therm(i) = s_therm(i) + therm_x(i,nsrf) * pctsrf(i,nsrf)
2908	s_trmb1(i) = s_trmb1(i) + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
2909	s_trmb2(i) = s_trmb2(i) + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
2910	s_trmb3(i) = s_trmb3(i) + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
2911	END DO
2912	END DO
2913	DO i = 1, klon
2914	qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
2915	END DO
2916	!!!
2917	ENDIF ! (iflag_split .eq.0)
2918	!!!
2919
2920	IF (check) THEN
2921	amn=MIN(ts(1,is_ter),1000.)
2922	amx=MAX(ts(1,is_ter),-1000.)
2923	DO i=2, klon
2924	amn=MIN(ts(i,is_ter),amn)
2925	amx=MAX(ts(i,is_ter),amx)
2926	ENDDO
2927	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
2928	ENDIF
2929
2930	!jg ?
2931	!!$!
2932	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
2933	!!$! sub-surfaces is distributed.
2934	!!$!
2935	!!$ DO nsrf = 1, nbsrf
2936	!!$ DO i = 1, klon
2937	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
2938	!!$ ts(i,nsrf) = zxtsol(i)
2939	!!$ t2m(i,nsrf) = zt2m(i)
2940	!!$ q2m(i,nsrf) = zq2m(i)
2941	!!$ u10m(i,nsrf) = zu10m(i)
2942	!!$ v10m(i,nsrf) = zv10m(i)
2943	!!$
2944	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
2945	!!$ pblh(i,nsrf) = s_pblh(i)
2946	!!$ plcl(i,nsrf) = s_plcl(i)
2947	!!$ capCL(i,nsrf) = s_capCL(i)
2948	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
2949	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
2950	!!$ pblT(i,nsrf) = s_pblT(i)
2951	!!$ therm(i,nsrf) = s_therm(i)
2952	!!$ trmb1(i,nsrf) = s_trmb1(i)
2953	!!$ trmb2(i,nsrf) = s_trmb2(i)
2954	!!$ trmb3(i,nsrf) = s_trmb3(i)
2955	!!$ ENDIF
2956	!!$ ENDDO
2957	!!$ ENDDO
2958
2959
2960	DO i = 1, klon
2961	fder(i) = - 4.0RSIGMAzxtsol(i)**3
2962	ENDDO
2963
2964	zxqsurf(:) = 0.0
2965	zxsnow(:) = 0.0
2966	DO nsrf = 1, nbsrf
2967	DO i = 1, klon
2968	zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
2969	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
2970	END DO
2971	END DO
2972
2973	! Premier niveau de vent sortie dans physiq.F
2974	zu1(:) = u(:,1)
2975	zv1(:) = v(:,1)
2976
2977
2978	END SUBROUTINE pbl_surface
2979	!
2980	!****************************************************************************************
2981	!
2982	SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
2983
2984	USE indice_sol_mod
2985
2986	INCLUDE "dimsoil.h"
2987
2988	! Ouput variables
2989	!****************************************************************************************
2990	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
2991	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
2992	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
2993	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
2994
2995
2996	!****************************************************************************************
2997	! Return module variables for writing to restart file
2998	!
2999	!****************************************************************************************
3000	fder_rst(:) = fder(:)
3001	snow_rst(:,:) = snow(:,:)
3002	qsurf_rst(:,:) = qsurf(:,:)
3003	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
3004
3005	!****************************************************************************************
3006	! Deallocate module variables
3007	!
3008	!****************************************************************************************
3009	! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
3010	IF (ALLOCATED(fder)) DEALLOCATE(fder)
3011	IF (ALLOCATED(snow)) DEALLOCATE(snow)
3012	IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
3013	IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
3014
3015	END SUBROUTINE pbl_surface_final
3016	!
3017	!****************************************************************************************
3018	!
3019
3020	!albedo SB >>>
3021	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
3022	evap, z0m, z0h, agesno, &
3023	tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)
3024	!albedo SB <<<
3025	! Give default values where new fraction has appread
3026
3027	USE indice_sol_mod
3028
3029	INCLUDE "dimsoil.h"
3030	INCLUDE "clesphys.h"
3031	INCLUDE "compbl.h"
3032
3033	! Input variables
3034	!****************************************************************************************
3035	INTEGER, INTENT(IN) :: itime
3036	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
3037
3038	! InOutput variables
3039	!****************************************************************************************
3040	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
3041	!albedo SB >>>
3042	REAL, DIMENSION(klon,nsw,nbsrf), INTENT(INOUT) :: alb_dir, alb_dif
3043	INTEGER :: k
3044	!albedo SB <<<
3045	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m
3046	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: evap, agesno
3047	REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT) :: z0m,z0h
3048	REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
3049
3050	! Local variables
3051	!****************************************************************************************
3052	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
3053	CHARACTER(len=80) :: abort_message
3054	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
3055	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
3056	!
3057	! All at once !!
3058	!****************************************************************************************
3059
3060	DO nsrf = 1, nbsrf
3061	! First decide complement sub-surfaces
3062	SELECT CASE (nsrf)
3063	CASE(is_oce)
3064	nsrf_comp1=is_sic
3065	nsrf_comp2=is_ter
3066	nsrf_comp3=is_lic
3067	CASE(is_sic)
3068	nsrf_comp1=is_oce
3069	nsrf_comp2=is_ter
3070	nsrf_comp3=is_lic
3071	CASE(is_ter)
3072	nsrf_comp1=is_lic
3073	nsrf_comp2=is_oce
3074	nsrf_comp3=is_sic
3075	CASE(is_lic)
3076	nsrf_comp1=is_ter
3077	nsrf_comp2=is_oce
3078	nsrf_comp3=is_sic
3079	END SELECT
3080
3081	! Initialize all new fractions
3082	DO i=1, klon
3083	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
3084
3085	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
3086	! Use the complement sub-surface, keeping the continents unchanged
3087	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
3088	evap(i,nsrf) = evap(i,nsrf_comp1)
3089	z0m(i,nsrf) = z0m(i,nsrf_comp1)
3090	z0h(i,nsrf) = z0h(i,nsrf_comp1)
3091	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
3092	!albedo SB >>>
3093	DO k=1,nsw
3094	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp1)
3095	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp1)
3096	ENDDO
3097	!albedo SB <<<
3098	ustar(i,nsrf) = ustar(i,nsrf_comp1)
3099	u10m(i,nsrf) = u10m(i,nsrf_comp1)
3100	v10m(i,nsrf) = v10m(i,nsrf_comp1)
3101	if (iflag_pbl > 1) then
3102	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
3103	endif
3104	mfois(nsrf) = mfois(nsrf) + 1
3105	! F. Codron sensible default values for ocean and sea ice
3106	IF (nsrf.EQ.is_oce) THEN
3107	tsurf(i,nsrf) = 271.35 ! Freezing sea water
3108	DO k=1,nsw
3109	alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
3110	alb_dif(i,k,nsrf) = 0.06
3111	END DO
3112	ELSE IF (nsrf.EQ.is_sic) THEN
3113	tsurf(i,nsrf) = 273.15 ! Melting ice
3114	DO k=1,nsw
3115	alb_dir(i,k,nsrf) = 0.3 ! thin ice
3116	alb_dif(i,k,nsrf) = 0.3
3117	END DO
3118	END IF
3119	! F. Codron
3120	ELSE
3121	! The continents have changed. The new fraction receives the mean sum of the existent fractions
3122	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3123	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3124	z0m(i,nsrf) = z0m(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0m(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3125	z0h(i,nsrf) = z0h(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0h(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3126	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3127	!albedo SB >>>
3128	DO k=1,nsw
3129	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3130	alb_dir(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3131	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3132	alb_dif(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3133	ENDDO
3134	!albedo SB <<<
3135	ustar(i,nsrf) = ustar(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3136	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3137	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3138	if (iflag_pbl > 1) then
3139	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3140	endif
3141
3142	! Security abort. This option has never been tested. To test, comment the following line.
3143	! abort_message='The fraction of the continents have changed!'
3144	! CALL abort_physic(modname,abort_message,1)
3145	nfois(nsrf) = nfois(nsrf) + 1
3146	END IF
3147	snow(i,nsrf) = 0.
3148	agesno(i,nsrf) = 0.
3149	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
3150	ELSE
3151	pfois(nsrf) = pfois(nsrf)+ 1
3152	END IF
3153	END DO
3154
3155	END DO
3156
3157	END SUBROUTINE pbl_surface_newfrac
3158
3159	!
3160	!****************************************************************************************
3161	!
3162
3163	END MODULE pbl_surface_mod
3164

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

Download in other formats:

Original Format