Context Navigation

pbl_surface_mod.F90 @ 5006

Last change on this file since 5006 was 3406, checked in by jghattas, 6 years ago
Added all modifications in the model code that were used for the simulations with DYANMICO during the Grand Challeng 2018. Modifications done by Y. Meurdesoif, L. Fairhead and A.K. Traore
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 129.8 KB

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

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

Download in other formats:

Original Format