Context Navigation

pbl_surface_mod.F90 @ 3780

Last change on this file since 3780 was 3780, checked in by evignon, 4 years ago
Premiere comission Etienne: changements pour le 1D (forcage en Ts au dessus des continents) et inclusion drag arbres dans yamada4_num=6
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 132.7 KB

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

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

Download in other formats:

Original Format