Context Navigation

pbl_surface_mod.F90 @ 3957

Last change on this file since 3957 was 3851, checked in by dcugnet, 4 years ago
Update the branch to the current trunk.
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.4 KB

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

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

Download in other formats:

Original Format