Context Navigation

pbl_surface_mod.F90 @ 3876

Last change on this file since 3876 was 3876, checked in by oboucher, 3 years ago
various minor tweaks
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: 137.5 KB

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

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

Download in other formats:

Original Format