Context Navigation

pbl_surface_mod.F90 @ 3823

Last change on this file since 3823 was 3823, checked in by musat, 3 years ago
Nouveaux calculs a 2m et 10m
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: 134.3 KB

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

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

Download in other formats:

Original Format