Context Navigation

pbl_surface_mod.F90 @ 3774

Last change on this file since 3774 was 3774, checked in by lguez, 4 years ago

Bug fix: encapsulate hbtm in a module

Fixes bug introduced in commit [3772]: therm was declared as an
assumed-shape dummy argument, which requires an explicit interface.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 132.4 KB

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

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

Download in other formats:

Original Format