Context Navigation

pbl_surface_mod.F90 @ 2934

Last change on this file since 2934 was 2924, checked in by fcheruy, 7 years ago

Creation of LMDZ branch to incorporate tree drag from ORCHIDEE.
Should merge in LMDZ trunk quickly

Cheruy with assistance from LF

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: 129.7 KB

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

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

Download in other formats:

Original Format