Context Navigation

pbl_surface_mod.F90 @ 2344

Last change on this file since 2344 was 2344, checked in by Ehouarn Millour, 9 years ago
Physics/dynamics separation: get rid of all the 'include "temps.h"' in the physics; variables in module time_phylmdz_mod must be used instead. Also added JD_cur, JH_cur and JD_ref in module phys_cal_mod, in preparation for having physics handle its calendar internally. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 125.0 KB

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

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

Download in other formats:

Original Format