Context Navigation

pbl_surface_mod.F90 @ 3324

Last change on this file since 3324 was 3200, checked in by Laurent Fairhead, 7 years ago

Inclusion of r3198 from trunk
Retour vers l'insensibilite au decoupage en sous domaine.
Les routines gwd_rando incluait le calcul de niveaux de reference
sur la base d'un profile pris au milieu du domaine (en klon/2).
Rempace par un test en presnivs.

Une autre intercation entre routines concernant la tke a fait apparaitre
que la tke n'?\195?\169tait pas passee correctement au niveau klev+1 au moment
du regroupement des mailles sous les sous surface.

Ces changements garantissent la convergence numerique si
addtkeoro=0
iflag_pbl<12
et
ok_gwd_rando=n
La convergence n'est pas garantie pour les dernieres versions des physiq.def
mais les differences devraient etre mineures.

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:keywords set to Author Date Id Revision

File size: 129.9 KB

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

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

Download in other formats:

Original Format