Context Navigation

pbl_surface_mod.F90 @ 3821

Last change on this file since 3821 was 3819, checked in by musat, 4 years ago
Derniere commission nvlle formulation t2m/q2m
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 3819 2021-02-03 20:10:38Z musat $
3	!
4	MODULE pbl_surface_mod
5	!
6	! Planetary Boundary Layer and Surface module
7	!
8	! This module manages the calculation of turbulent diffusion in the boundary layer
9	! and all interactions towards the differents sub-surfaces.
10	!
11	!
12	USE dimphy
13	USE mod_phys_lmdz_para, ONLY : mpi_size
14	USE mod_grid_phy_lmdz, ONLY : klon_glo
15	USE ioipsl
16	USE surface_data, ONLY : type_ocean, ok_veget
17	USE surf_land_mod, ONLY : surf_land
18	USE surf_landice_mod, ONLY : surf_landice
19	USE surf_ocean_mod, ONLY : surf_ocean
20	USE surf_seaice_mod, ONLY : surf_seaice
21	USE cpl_mod, ONLY : gath2cpl
22	USE climb_hq_mod, ONLY : climb_hq_down, climb_hq_up
23	USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up
24	USE coef_diff_turb_mod, ONLY : coef_diff_turb
25	USE ioipsl_getin_p_mod, ONLY : getin_p
26	USE cdrag_mod
27	USE stdlevvar_mod
28	USE wx_pbl_mod, ONLY : wx_pbl_init, wx_pbl_final, &
29	!! wx_pbl_fuse_no_dts, wx_pbl_split_no_dts, &
30	!! wx_pbl_fuse, wx_pbl_split
31	wx_pbl0_fuse, wx_pbl0_split
32	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, yn2mout(:, nsrf, :))
1879	ELSE
1880	CALL stdlevvar(klon, knon, is_ter, zxli, &
1881	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
1882	yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
1883	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
1884	ENDIF
1885
1886	ENDIF
1887
1888	!****************************************************************************************
1889	!
1890	! 10) Switch according to current surface
1891	! It is necessary to start with the continental surfaces because the ocean
1892	! needs their run-off.
1893	!
1894	!****************************************************************************************
1895	SELECT CASE(nsrf)
1896
1897	CASE(is_ter)
1898	! print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)
1899	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
1900	rlon, rlat, yrmu0, &
1901	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
1902	!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1903	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
1904	AcoefH, AcoefQ, BcoefH, BcoefQ, &
1905	AcoefU, AcoefV, BcoefU, BcoefV, &
1906	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1907	ylwdown, yq2m, yt2m, &
1908	ysnow, yqsol, yagesno, ytsoil, &
1909	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
1910	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
1911	y_flux_u1, y_flux_v1, &
1912	yveget,ylai,yheight )
1913
1914	!FC quid qd yveget ylai yheight ne sont pas definit
1915	!FC yveget,ylai,yheight, &
1916	IF (ifl_pbltree .ge. 1) THEN
1917	CALL freinage(knon, yu, yv, yt, &
1918	! yveget,ylai, yheight,ypaprs,ypplay,y_d_u_frein,y_d_v_frein)
1919	yveget,ylai, yheight,ypaprs,ypplay,y_treedrg, y_d_u_frein,y_d_v_frein)
1920	ENDIF
1921
1922
1923	! Special DICE MPL 05082013 puis BOMEX
1924	IF (ok_prescr_ust) THEN
1925	DO j=1,knon
1926	! ysnow(:)=0.
1927	! yqsol(:)=0.
1928	! yagesno(:)=50.
1929	! ytsoil(:,:)=300.
1930	! yz0_new(:)=0.001
1931	! yevap(:)=flat/RLVTT
1932	! yfluxlat(:)=-flat
1933	! yfluxsens(:)=-fsens
1934	! yqsurf(:)=0.
1935	! ytsurf_new(:)=tg
1936	! y_dflux_t(:)=0.
1937	! y_dflux_q(:)=0.
1938	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)
1939	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)
1940	ENDDO
1941	ENDIF
1942
1943	CASE(is_lic)
1944	! Martin
1945	CALL surf_landice(itap, dtime, knon, ni, &
1946	rlon, rlat, debut, lafin, &
1947	yrmu0, ylwdown, yalb, ypphi(:,1), &
1948	ysolsw, ysollw, yts, ypplay(:,1), &
1949	!!jyg ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1950	ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
1951	AcoefH, AcoefQ, BcoefH, BcoefQ, &
1952	AcoefU, AcoefV, BcoefU, BcoefV, &
1953	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1954	ysnow, yqsurf, yqsol, yagesno, &
1955	ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
1956	ytsurf_new, y_dflux_t, y_dflux_q, &
1957	yzsig, ycldt, &
1958	ysnowhgt, yqsnow, ytoice, ysissnow, &
1959	yalb3_new, yrunoff, &
1960	y_flux_u1, y_flux_v1)
1961
1962	!jyg<
1963	!! alb3_lic(:)=0.
1964	!>jyg
1965	DO j = 1, knon
1966	i = ni(j)
1967	alb3_lic(i) = yalb3_new(j)
1968	snowhgt(i) = ysnowhgt(j)
1969	qsnow(i) = yqsnow(j)
1970	to_ice(i) = ytoice(j)
1971	sissnow(i) = ysissnow(j)
1972	runoff(i) = yrunoff(j)
1973	ENDDO
1974	! Martin
1975	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
1976	IF (ok_prescr_ust) THEN
1977	DO j=1,knon
1978	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)
1979	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)
1980	ENDDO
1981	ENDIF
1982
1983	CASE(is_oce)
1984	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
1985	ywindsp, rmu0, yfder, yts, &
1986	itap, dtime, jour, knon, ni, &
1987	!!jyg ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1988	ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& ! ym missing init
1989	AcoefH, AcoefQ, BcoefH, BcoefQ, &
1990	AcoefU, AcoefV, BcoefU, BcoefV, &
1991	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1992	ysnow, yqsurf, yagesno, &
1993	yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
1994	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
1995	y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
1996	yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
1997	ytkt(:knon), ytks(:knon), ytaur(:knon), ysss)
1998	IF (prt_level >=10) THEN
1999	print *,'arg de surf_ocean: ycdragh ',ycdragh
2000	print *,'arg de surf_ocean: ycdragm ',ycdragm
2001	print *,'arg de surf_ocean: yt ', yt
2002	print *,'arg de surf_ocean: yq ', yq
2003	print *,'arg de surf_ocean: yts ', yts
2004	print *,'arg de surf_ocean: AcoefH ',AcoefH
2005	print *,'arg de surf_ocean: AcoefQ ',AcoefQ
2006	print *,'arg de surf_ocean: BcoefH ',BcoefH
2007	print *,'arg de surf_ocean: BcoefQ ',BcoefQ
2008	print *,'arg de surf_ocean: yevap ',yevap
2009	print *,'arg de surf_ocean: yfluxsens ',yfluxsens
2010	print *,'arg de surf_ocean: yfluxlat ',yfluxlat
2011	print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new
2012	ENDIF
2013	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2014	IF (ok_prescr_ust) THEN
2015	DO j=1,knon
2016	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)
2017	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)
2018	ENDDO
2019	ENDIF
2020
2021	CASE(is_sic)
2022	CALL surf_seaice( &
2023	!albedo SB >>>
2024	rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, &
2025	!albedo SB <<<
2026	itap, dtime, jour, knon, ni, &
2027	lafin, &
2028	!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
2029	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
2030	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2031	AcoefU, AcoefV, BcoefU, BcoefV, &
2032	ypsref, yu1, yv1, ygustiness, pctsrf, &
2033	ysnow, yqsurf, yqsol, yagesno, ytsoil, &
2034	!albedo SB >>>
2035	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
2036	!albedo SB <<<
2037	ytsurf_new, y_dflux_t, y_dflux_q, &
2038	y_flux_u1, y_flux_v1)
2039
2040	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2041	IF (ok_prescr_ust) THEN
2042	DO j=1,knon
2043	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)
2044	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)
2045	ENDDO
2046	ENDIF
2047
2048	CASE DEFAULT
2049	WRITE(lunout,*) 'Surface index = ', nsrf
2050	abort_message = 'Surface index not valid'
2051	CALL abort_physic(modname,abort_message,1)
2052	END SELECT
2053
2054
2055	!****************************************************************************************
2056	! 11) - Calcul the increment of surface temperature
2057	!
2058	!****************************************************************************************
2059
2060	IF (evap0>=0.) THEN
2061	yevap(:)=evap0
2062	yevap(:)=RLVTT*evap0
2063	ENDIF
2064
2065	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
2066
2067	!****************************************************************************************
2068	!
2069	! 12) "La remontee" - "The uphill"
2070	!
2071	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
2072	! for X=H, Q, U and V, for all vertical levels.
2073	!
2074	!****************************************************************************************
2075	!!
2076	!!!
2077	!!! jyg le 10/04/2013 et EV 10/2020
2078
2079	IF (ok_forc_tsurf) THEN
2080	DO j=1,knon
2081	ytsurf_new(j)=tg
2082	y_d_ts(j) = ytsurf_new(j) - yts(j)
2083	ENDDO
2084	ENDIF ! ok_forc_tsurf
2085
2086	!!!
2087	IF (ok_flux_surf) THEN
2088	IF (prt_level >=10) THEN
2089	PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
2090	ENDIF
2091	y_flux_t1(:) = fsens
2092	y_flux_q1(:) = flat/RLVTT
2093	yfluxlat(:) = flat
2094	!
2095	!! Test sur iflag_split retire le 2/02/2018, sans vraiment comprendre la raison de ce test. (jyg)
2096	!! IF (iflag_split .eq.0) THEN
2097	DO j=1,knon
2098	Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)2+yv(j,1)2)) * &
2099	ypplay(j,1)/(RD*yt(j,1))
2100	ENDDO
2101	!! ENDIF ! (iflag_split .eq.0)
2102
2103	DO j = 1, knon
2104	yt1_new=(1./RCPD)(AcoefH(j)+BcoefH(j)y_flux_t1(j)*dtime)
2105	ytsurf_new(j)=yt1_new-y_flux_t1(j)/(Kech_h(j)*RCPD)
2106	ENDDO
2107
2108	DO j=1,knon
2109	y_d_ts(j) = ytsurf_new(j) - yts(j)
2110	ENDDO
2111
2112	ELSE ! (ok_flux_surf)
2113	DO j=1,knon
2114	y_flux_t1(j) = yfluxsens(j)
2115	y_flux_q1(j) = -yevap(j)
2116	ENDDO
2117	ENDIF
2118
2119	IF (prt_level >=10) THEN
2120	DO j=1,knon
2121	print*,'y_flux_t1,yfluxlat,wakes' &
2122	& , y_flux_t1(j), yfluxlat(j), ywake_s(j)
2123	print*,'beta,ytsurf_new', ybeta(j), ytsurf_new(j)
2124	print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
2125	ENDDO
2126	ENDIF
2127
2128	!!! jyg le 07/02/2012 puis le 10/04/2013
2129	!! IF (iflag_split .eq.1) THEN
2130	!!!!!
2131	!!!jyg<
2132	!! CALL wx_pbl_split_no_dts(knon, ywake_s, &
2133	!! AcoefH_x, AcoefH_w, &
2134	!! AcoefQ_x, AcoefQ_w, &
2135	!! AcoefU_x, AcoefU_w, &
2136	!! AcoefV_x, AcoefV_w, &
2137	!! y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
2138	!! y_flux_t1_x, y_flux_t1_w, &
2139	!! y_flux_q1_x, y_flux_q1_w, &
2140	!! y_flux_u1_x, y_flux_u1_w, &
2141	!! y_flux_v1_x, y_flux_v1_w, &
2142	!! yfluxlat_x, yfluxlat_w &
2143	!! )
2144	!! ELSE IF (iflag_split .ge. 2) THEN
2145	IF (iflag_split .GE. 1) THEN
2146	CALL wx_pbl0_split(knon, dtime, ywake_s, &
2147	y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
2148	y_flux_t1_x, y_flux_t1_w, &
2149	y_flux_q1_x, y_flux_q1_w, &
2150	y_flux_u1_x, y_flux_u1_w, &
2151	y_flux_v1_x, y_flux_v1_w, &
2152	yfluxlat_x, yfluxlat_w, &
2153	y_delta_tsurf &
2154	)
2155	ENDIF ! (iflag_split .ge. 1)
2156	!>jyg
2157	!
2158
2159	!!jyg!! A reprendre apres reflexion ===============================================
2160	!!jyg!!
2161	!!jyg!! DO j=1,knon
2162	!!jyg!!!!! nrlmd le 13/06/2011
2163	!!jyg!!
2164	!!jyg!!!----Diffusion dans le sol dans le cas continental seulement
2165	!!jyg!! IF (nsrf.eq.is_ter) THEN
2166	!!jyg!!!----Calcul du coefficient delta_coeff
2167	!!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)))
2168	!!jyg!!
2169	!!jyg!!! delta_coef(j)=dtime/(inertia*sqrt(tau_eq(j)))
2170	!!jyg!! delta_coef(j)=facteur*sqrt(tau_eq(j))/inertia
2171	!!jyg!!! delta_coef(j)=0.
2172	!!jyg!! ELSE
2173	!!jyg!! delta_coef(j)=0.
2174	!!jyg!! ENDIF
2175	!!jyg!!
2176	!!jyg!!!----Calcul de delta_tsurf
2177	!!jyg!! y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j)
2178	!!jyg!!
2179	!!jyg!!!----Si il n'y a pas des poches...
2180	!!jyg!! IF (wake_cstar(j).le.0.01) THEN
2181	!!jyg!! y_delta_tsurf(j)=0.
2182	!!jyg!! y_delta_flux_t1(j)=0.
2183	!!jyg!! ENDIF
2184	!!jyg!!
2185	!!jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle)
2186	!!jyg!!!!!!! jyg le 23/02/2012
2187	!!jyg!!!!!!!
2188	!!jyg!!!! ybeta(j)=y_flux_q1(j) / &
2189	!!jyg!!!! & (Kech_h(j)*(yq(j,1)-yqsatsurf(j)))
2190	!!jyg!!!!!! ybeta(j)=-1.*yevap(j) / &
2191	!!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)))
2192	!!jyg!!!!!!! fin jyg
2193	!!jyg!!!!!
2194	!!jyg!!
2195	!!jyg!! ENDDO
2196	!!jyg!!
2197	!!jyg!!!!! fin nrlmd le 13/06/2011
2198	!!jyg!!
2199	IF (iflag_split .ge. 1) THEN
2200	IF (prt_level >=10) THEN
2201	DO j = 1, knon
2202	print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
2203	print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
2204	! print*,'tsurf_x,tsurf_w,tsurf,t1', ytsurf_th_x(j), ytsurf_th_w(j), ytsurf_th(j), yt(j,1)
2205	print*,'tsurf_x,t1x,tsurf_w,t1w,tsurf,t1,t1_ancien', &
2206	& 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)
2207	print*,'qsatsurf,qsatsurf_x,qsatsurf_w', yqsatsurf(j), yqsatsurf_x(j), yqsatsurf_w(j)
2208	print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j)
2209	ENDDO
2210
2211	DO j=1,knon
2212	print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
2213	& , 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)
2214	print*,'beta,ytsurf_new,yqsatsurf', ybeta(j), ytsurf_new(j), yqsatsurf(j)
2215	print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
2216	ENDDO
2217	ENDIF ! (prt_level >=10)
2218
2219	!!! jyg le 07/02/2012
2220	ENDIF ! (iflag_split .ge.1)
2221	!!!
2222
2223	!!! jyg le 07/02/2012
2224	IF (iflag_split .eq.0) THEN
2225	!!!
2226	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2227	CALL climb_hq_up(knon, dtime, yt, yq, &
2228	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
2229	!!! jyg le 07/02/2012
2230	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2231	CcoefH, CcoefQ, DcoefH, DcoefQ, &
2232	Kcoef_hq, gama_q, gama_h, &
2233	!!!
2234	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
2235	ELSE !(iflag_split .eq.0)
2236	CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
2237	y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &
2238	!!! nrlmd le 02/05/2011
2239	AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &
2240	CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
2241	Kcoef_hq_x, gama_q_x, gama_h_x, &
2242	!!!
2243	y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:))
2244	!
2245	CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
2246	y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &
2247	!!! nrlmd le 02/05/2011
2248	AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &
2249	CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
2250	Kcoef_hq_w, gama_q_w, gama_h_w, &
2251	!!!
2252	y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:))
2253	!!!
2254	ENDIF ! (iflag_split .eq.0)
2255	!!!
2256
2257	!!! jyg le 07/02/2012
2258	IF (iflag_split .eq.0) THEN
2259	!!!
2260	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2261	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
2262	!!! jyg le 07/02/2012
2263	AcoefU, AcoefV, BcoefU, BcoefV, &
2264	CcoefU, CcoefV, DcoefU, DcoefV, &
2265	Kcoef_m, &
2266	!!!
2267	y_flux_u, y_flux_v, y_d_u, y_d_v)
2268	y_d_t_diss(:,:)=0.
2269	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2270	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2271	& ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
2272	& ,iflag_pbl)
2273	ENDIF
2274	! print*,'yamada_c OK'
2275
2276	ELSE !(iflag_split .eq.0)
2277	CALL climb_wind_up(knon, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &
2278	!!! nrlmd le 02/05/2011
2279	AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &
2280	CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
2281	Kcoef_m_x, &
2282	!!!
2283	y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)
2284	!
2285	y_d_t_diss_x(:,:)=0.
2286	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2287	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2288	& ,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 &
2289	,ycoefq_x,y_d_t_diss_x,yustar_x &
2290	& ,iflag_pbl)
2291	ENDIF
2292	! print*,'yamada_c OK'
2293
2294	CALL climb_wind_up(knon, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &
2295	!!! nrlmd le 02/05/2011
2296	AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &
2297	CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
2298	Kcoef_m_w, &
2299	!!!
2300	y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)
2301	!!!
2302	y_d_t_diss_w(:,:)=0.
2303	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2304	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2305	& ,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 &
2306	,ycoefq_w,y_d_t_diss_w,yustar_w &
2307	& ,iflag_pbl)
2308	ENDIF
2309	! print*,'yamada_c OK'
2310	!
2311	IF (prt_level >=10) THEN
2312	print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &
2313	yfluxlat_x, yfluxlat_w
2314	ENDIF
2315	!
2316	ENDIF ! (iflag_split .eq.0)
2317	!!!
2318
2319	DO j = 1, knon
2320	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
2321	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
2322	ENDDO
2323
2324	!****************************************************************************************
2325	! 13) Transform variables for output format :
2326	! - Decompress
2327	! - Multiply with pourcentage of current surface
2328	! - Cumulate in global variable
2329	!
2330	!****************************************************************************************
2331
2332
2333	!!! jyg le 07/02/2012
2334	IF (iflag_split.EQ.0) THEN
2335	!!!
2336	DO k = 1, klev
2337	DO j = 1, knon
2338	i = ni(j)
2339	y_d_t_diss(j,k) = y_d_t_diss(j,k) * ypct(j)
2340	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
2341	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
2342	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
2343	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
2344	!FC
2345	IF (nsrf .EQ. is_ter .and. ifl_pbltree .GE. 1) THEN
2346	! if (y_d_u_frein(j,k).ne.0. ) then
2347	! print, nsrf,'IS_TER ++', y_d_u_frein(j,k)ypct(j),y_d_u(j,k),j,k
2348	! ENDIF
2349	y_d_u(j,k) =y_d_u(j,k) + y_d_u_frein(j,k)*ypct(j)
2350	y_d_v(j,k) =y_d_v(j,k) + y_d_v_frein(j,k)*ypct(j)
2351	treedrg(i,k,nsrf)=y_treedrg(j,k)
2352	ELSE
2353	treedrg(i,k,nsrf)=0.
2354	ENDIF
2355	!FC
2356	flux_t(i,k,nsrf) = y_flux_t(j,k)
2357	flux_q(i,k,nsrf) = y_flux_q(j,k)
2358	flux_u(i,k,nsrf) = y_flux_u(j,k)
2359	flux_v(i,k,nsrf) = y_flux_v(j,k)
2360	ENDDO
2361	ENDDO
2362
2363	ELSE !(iflag_split .eq.0)
2364
2365	! Tendances hors poches
2366	DO k = 1, klev
2367	DO j = 1, knon
2368	i = ni(j)
2369	y_d_t_diss_x(j,k) = y_d_t_diss_x(j,k) * ypct(j)
2370	y_d_t_x(j,k) = y_d_t_x(j,k) * ypct(j)
2371	y_d_q_x(j,k) = y_d_q_x(j,k) * ypct(j)
2372	y_d_u_x(j,k) = y_d_u_x(j,k) * ypct(j)
2373	y_d_v_x(j,k) = y_d_v_x(j,k) * ypct(j)
2374
2375	flux_t_x(i,k,nsrf) = y_flux_t_x(j,k)
2376	flux_q_x(i,k,nsrf) = y_flux_q_x(j,k)
2377	flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
2378	flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)
2379	ENDDO
2380	ENDDO
2381
2382	! Tendances dans les poches
2383	DO k = 1, klev
2384	DO j = 1, knon
2385	i = ni(j)
2386	y_d_t_diss_w(j,k) = y_d_t_diss_w(j,k) * ypct(j)
2387	y_d_t_w(j,k) = y_d_t_w(j,k) * ypct(j)
2388	y_d_q_w(j,k) = y_d_q_w(j,k) * ypct(j)
2389	y_d_u_w(j,k) = y_d_u_w(j,k) * ypct(j)
2390	y_d_v_w(j,k) = y_d_v_w(j,k) * ypct(j)
2391
2392	flux_t_w(i,k,nsrf) = y_flux_t_w(j,k)
2393	flux_q_w(i,k,nsrf) = y_flux_q_w(j,k)
2394	flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
2395	flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)
2396	ENDDO
2397	ENDDO
2398
2399	! Flux, tendances et Tke moyenne dans la maille
2400	DO k = 1, klev
2401	DO j = 1, knon
2402	i = ni(j)
2403	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))
2404	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))
2405	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))
2406	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))
2407	ENDDO
2408	ENDDO
2409	DO j=1,knon
2410	yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j))
2411	ENDDO
2412	IF (prt_level >=10) THEN
2413	print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &
2414	nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf)
2415	ENDIF
2416
2417	DO k = 1, klev
2418	DO j = 1, knon
2419	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))
2420	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))
2421	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))
2422	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))
2423	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))
2424	ENDDO
2425	ENDDO
2426
2427	ENDIF ! (iflag_split .eq.0)
2428	!!!
2429
2430	! print*,'Dans pbl OK1'
2431
2432	!jyg<
2433	!! evap(:,nsrf) = - flux_q(:,1,nsrf)
2434	!>jyg
2435	DO j = 1, knon
2436	i = ni(j)
2437	evap(i,nsrf) = - flux_q(i,1,nsrf) !jyg
2438	d_ts(i,nsrf) = y_d_ts(j)
2439	!albedo SB >>>
2440	DO k=1,nsw
2441	alb_dir(i,k,nsrf) = yalb_dir_new(j,k)
2442	alb_dif(i,k,nsrf) = yalb_dif_new(j,k)
2443	ENDDO
2444	!albedo SB <<<
2445	snow(i,nsrf) = ysnow(j)
2446	qsurf(i,nsrf) = yqsurf(j)
2447	z0m(i,nsrf) = yz0m(j)
2448	z0h(i,nsrf) = yz0h(j)
2449	fluxlat(i,nsrf) = yfluxlat(j)
2450	agesno(i,nsrf) = yagesno(j)
2451	cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
2452	cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
2453	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
2454	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
2455	ENDDO
2456
2457	! print*,'Dans pbl OK2'
2458
2459	!!! jyg le 07/02/2012
2460	IF (iflag_split .ge.1) THEN
2461	!!!
2462	!!! nrlmd le 02/05/2011
2463	DO j = 1, knon
2464	i = ni(j)
2465	fluxlat_x(i,nsrf) = yfluxlat_x(j)
2466	fluxlat_w(i,nsrf) = yfluxlat_w(j)
2467	!!!
2468	!!! nrlmd le 13/06/2011
2469	!!jyg20170131 delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
2470	delta_tsurf(i,nsrf)=y_delta_tsurf(j)
2471	!
2472	cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
2473	cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)
2474	cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j)
2475	cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j)
2476	kh(i) = kh(i) + Kech_h(j)*ypct(j)
2477	kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j)
2478	kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j)
2479	!!!
2480	ENDDO
2481	!!!
2482	ENDIF ! (iflag_split .ge.1)
2483	!!!
2484	!!! nrlmd le 02/05/2011
2485	!!jyg le 20/02/2011
2486	!! tke_x(:,:,nsrf)=0.
2487	!! tke_w(:,:,nsrf)=0.
2488	!!jyg le 20/02/2011
2489	!! DO k = 1, klev+1
2490	!! DO j = 1, knon
2491	!! i = ni(j)
2492	!! wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2493	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2494	!! ENDDO
2495	!! ENDDO
2496	!!jyg le 20/02/2011
2497	!! DO k = 1, klev+1
2498	!! DO j = 1, knon
2499	!! i = ni(j)
2500	!! tke(i,k,nsrf)=(1.-ywake_s(j))tke_x(i,k,nsrf)+ywake_s(j)tke_w(i,k,nsrf)
2501	!! ENDDO
2502	!! ENDDO
2503	!!!
2504	IF (iflag_split .eq.0) THEN
2505	wake_dltke(:,:,nsrf) = 0.
2506	DO k = 1, klev+1
2507	DO j = 1, knon
2508	i = ni(j)
2509	!jyg<
2510	!! tke(i,k,nsrf) = ytke(j,k)
2511	!! tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
2512	tke_x(i,k,nsrf) = ytke(j,k)
2513	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
2514
2515	!>jyg
2516	ENDDO
2517	ENDDO
2518
2519	ELSE ! (iflag_split .eq.0)
2520	DO k = 1, klev+1
2521	DO j = 1, knon
2522	i = ni(j)
2523	wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2524	!jyg<
2525	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2526	!! tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
2527	tke_x(i,k,nsrf) = ytke_x(j,k)
2528	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)
2529	wake_dltke(i,k,is_ave) = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
2530
2531
2532	!>jyg
2533	ENDDO
2534	ENDDO
2535	ENDIF ! (iflag_split .eq.0)
2536	!!!
2537	DO k = 2, klev
2538	DO j = 1, knon
2539	i = ni(j)
2540	zcoefh(i,k,nsrf) = ycoefh(j,k)
2541	zcoefm(i,k,nsrf) = ycoefm(j,k)
2542	zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
2543	zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
2544	ENDDO
2545	ENDDO
2546
2547	! print*,'Dans pbl OK3'
2548
2549	IF ( nsrf .EQ. is_ter ) THEN
2550	DO j = 1, knon
2551	i = ni(j)
2552	qsol(i) = yqsol(j)
2553	ENDDO
2554	ENDIF
2555
2556	!jyg<
2557	!! ftsoil(:,:,nsrf) = 0.
2558	!>jyg
2559	DO k = 1, nsoilmx
2560	DO j = 1, knon
2561	i = ni(j)
2562	ftsoil(i, k, nsrf) = ytsoil(j,k)
2563	ENDDO
2564	ENDDO
2565
2566	!!! jyg le 07/02/2012
2567	IF (iflag_split .ge.1) THEN
2568	!!!
2569	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
2570	DO k = 1, klev
2571	DO j = 1, knon
2572	i = ni(j)
2573	d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)
2574	d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)
2575	d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)
2576	d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)
2577	d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)
2578	!
2579	d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)
2580	d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)
2581	d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)
2582	d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
2583	d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
2584	!
2585	!! d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
2586	!! d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)
2587	ENDDO
2588	ENDDO
2589	!!!
2590	ENDIF ! (iflag_split .ge.1)
2591	!!!
2592
2593	DO k = 1, klev
2594	DO j = 1, knon
2595	i = ni(j)
2596	d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
2597	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
2598	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
2599	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
2600	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
2601	ENDDO
2602	ENDDO
2603
2604	! print*,'Dans pbl OK4'
2605
2606	IF (prt_level >=10) THEN
2607	print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
2608	d_t_w(:,1), d_t_x(:,1), d_t(:,1)
2609	ENDIF
2610
2611	if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
2612	delta_sal = missing_val
2613	ds_ns = missing_val
2614	dt_ns = missing_val
2615	delta_sst = missing_val
2616	dter = missing_val
2617	dser = missing_val
2618	tkt = missing_val
2619	tks = missing_val
2620	taur = missing_val
2621	sss = missing_val
2622
2623	delta_sal(ni(:knon)) = ydelta_sal(:knon)
2624	ds_ns(ni(:knon)) = yds_ns(:knon)
2625	dt_ns(ni(:knon)) = ydt_ns(:knon)
2626	delta_sst(ni(:knon)) = ydelta_sst(:knon)
2627	dter(ni(:knon)) = ydter(:knon)
2628	dser(ni(:knon)) = ydser(:knon)
2629	tkt(ni(:knon)) = ytkt(:knon)
2630	tks(ni(:knon)) = ytks(:knon)
2631	taur(ni(:knon)) = ytaur(:knon)
2632	sss(ni(:knon)) = ysss(:knon)
2633	end if
2634
2635	!****************************************************************************************
2636	! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m
2637	! Call HBTM
2638	!
2639	!****************************************************************************************
2640	!!!
2641	!
2642	#undef T2m
2643	#define T2m
2644	#ifdef T2m
2645	! Calculations of diagnostic t,q at 2m and u, v at 10m
2646
2647	! print*,'Dans pbl OK41'
2648	! print*,'tair1,yt(:,1),y_d_t(:,1)'
2649	! print*, tair1,yt(:,1),y_d_t(:,1)
2650	!!! jyg le 07/02/2012
2651	IF (iflag_split .eq.0) THEN
2652	DO j=1, knon
2653	uzon(j) = yu(j,1) + y_d_u(j,1)
2654	vmer(j) = yv(j,1) + y_d_v(j,1)
2655	tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
2656	qair1(j) = yq(j,1) + y_d_q(j,1)
2657	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2658	* (ypaprs(j,1)-ypplay(j,1))
2659	tairsol(j) = yts(j) + y_d_ts(j)
2660	qairsol(j) = yqsurf(j)
2661	ENDDO
2662	ELSE ! (iflag_split .eq.0)
2663	DO j=1, knon
2664	uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
2665	vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
2666	tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
2667	qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
2668	zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2669	* (ypaprs(j,1)-ypplay(j,1))
2670	tairsol(j) = yts(j) + y_d_ts(j)
2671	tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
2672	qairsol(j) = yqsurf(j)
2673	ENDDO
2674	DO j=1, knon
2675	uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
2676	vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
2677	tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
2678	qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
2679	zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2680	* (ypaprs(j,1)-ypplay(j,1))
2681	tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
2682	qairsol(j) = yqsurf(j)
2683	ENDDO
2684	!!!
2685	ENDIF ! (iflag_split .eq.0)
2686	!!!
2687	DO j=1, knon
2688	! i = ni(j)
2689	! yz0h_oupas(j) = yz0m(j)
2690	! IF(nsrf.EQ.is_oce) THEN
2691	! yz0h_oupas(j) = z0m(i,nsrf)
2692	! ENDIF
2693	psfce(j)=ypaprs(j,1)
2694	patm(j)=ypplay(j,1)
2695	ENDDO
2696
2697	IF (iflag_pbl_surface_t2m_bug==1) THEN
2698	yz0h_oupas(1:knon)=yz0m(1:knon)
2699	ELSE
2700	yz0h_oupas(1:knon)=yz0h(1:knon)
2701	ENDIF
2702
2703	! print*,'Dans pbl OK42A'
2704	! print*,'tair1,yt(:,1),y_d_t(:,1)'
2705	! print*, tair1,yt(:,1),y_d_t(:,1)
2706
2707	! Calculate the temperature and relative humidity at 2m and the wind at 10m
2708	!!! jyg le 07/02/2012
2709	IF (iflag_split .eq.0) THEN
2710	IF (iflag_new_t2mq2m==1) THEN
2711	CALL stdlevvarn(klon, knon, nsrf, zxli, &
2712	uzon, vmer, tair1, qair1, zgeo1, &
2713	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
2714	yt2m, yq2m, yt10m, yq10m, yu10m, yn2mout(:, nsrf, :))
2715	ELSE
2716	CALL stdlevvar(klon, knon, nsrf, zxli, &
2717	uzon, vmer, tair1, qair1, zgeo1, &
2718	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
2719	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
2720	ENDIF
2721	ELSE !(iflag_split .eq.0)
2722	IF (iflag_new_t2mq2m==1) THEN
2723	CALL stdlevvarn(klon, knon, nsrf, zxli, &
2724	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
2725	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
2726	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yn2mout_x(:, nsrf, :))
2727	CALL stdlevvarn(klon, knon, nsrf, zxli, &
2728	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
2729	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
2730	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yn2mout_w(:, nsrf, :))
2731	ELSE
2732	CALL stdlevvar(klon, knon, nsrf, zxli, &
2733	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
2734	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
2735	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
2736	CALL stdlevvar(klon, knon, nsrf, zxli, &
2737	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
2738	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
2739	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
2740	ENDIF
2741	!!!
2742	ENDIF ! (iflag_split .eq.0)
2743	!!!
2744	!!! jyg le 07/02/2012
2745	IF (iflag_split .eq.0) THEN
2746	DO j=1, knon
2747	i = ni(j)
2748	t2m(i,nsrf)=yt2m(j)
2749	q2m(i,nsrf)=yq2m(j)
2750	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2751	ustar(i,nsrf)=yustar(j)
2752	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
2753	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
2754	!
2755	DO k = 1, 6
2756	n2mout(i,nsrf,k) = yn2mout(j,nsrf,k)
2757	END DO
2758	!
2759	ENDDO
2760	ELSE !(iflag_split .eq.0)
2761	DO j=1, knon
2762	i = ni(j)
2763	t2m_x(i,nsrf)=yt2m_x(j)
2764	q2m_x(i,nsrf)=yq2m_x(j)
2765	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2766	ustar_x(i,nsrf)=yustar_x(j)
2767	u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
2768	v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
2769	!
2770	DO k = 1, 6
2771	n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k)
2772	END DO
2773	!
2774	ENDDO
2775	DO j=1, knon
2776	i = ni(j)
2777	t2m_w(i,nsrf)=yt2m_w(j)
2778	q2m_w(i,nsrf)=yq2m_w(j)
2779	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2780	ustar_w(i,nsrf)=yustar_w(j)
2781	u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
2782	v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
2783	!
2784	ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
2785	u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
2786	v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
2787	!
2788	DO k = 1, 6
2789	n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k)
2790	END DO
2791	!
2792	ENDDO
2793	!!!
2794	ENDIF ! (iflag_split .eq.0)
2795	!!!
2796
2797	! print*,'Dans pbl OK43'
2798	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
2799	!IM Ajoute dependance type surface
2800	IF (thermcep) THEN
2801	!!! jyg le 07/02/2012
2802	IF (iflag_split .eq.0) THEN
2803	DO j = 1, knon
2804	i=ni(j)
2805	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
2806	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
2807	zx_qs1 = MIN(0.5,zx_qs1)
2808	zcor1 = 1./(1.-RETV*zx_qs1)
2809	zx_qs1 = zx_qs1*zcor1
2810
2811	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
2812	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
2813	ENDDO
2814	ELSE ! (iflag_split .eq.0)
2815	DO j = 1, knon
2816	i=ni(j)
2817	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
2818	zx_qs1 = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)
2819	zx_qs1 = MIN(0.5,zx_qs1)
2820	zcor1 = 1./(1.-RETV*zx_qs1)
2821	zx_qs1 = zx_qs1*zcor1
2822
2823	rh2m_x(i) = rh2m_x(i) + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
2824	qsat2m_x(i) = qsat2m_x(i) + zx_qs1 * pctsrf(i,nsrf)
2825	ENDDO
2826	DO j = 1, knon
2827	i=ni(j)
2828	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
2829	zx_qs1 = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
2830	zx_qs1 = MIN(0.5,zx_qs1)
2831	zcor1 = 1./(1.-RETV*zx_qs1)
2832	zx_qs1 = zx_qs1*zcor1
2833
2834	rh2m_w(i) = rh2m_w(i) + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
2835	qsat2m_w(i) = qsat2m_w(i) + zx_qs1 * pctsrf(i,nsrf)
2836	ENDDO
2837	!!!
2838	ENDIF ! (iflag_split .eq.0)
2839	!!!
2840	ENDIF
2841	!
2842	IF (prt_level >=10) THEN
2843	print *, 'T2m, q2m, RH2m ', &
2844	t2m, q2m, rh2m
2845	ENDIF
2846
2847	! print*,'OK pbl 5'
2848	!
2849	!!! jyg le 07/02/2012
2850	IF (iflag_split .eq.0) THEN
2851	CALL hbtm(knon, ypaprs, ypplay, &
2852	yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
2853	y_flux_t,y_flux_q,yu,yv,yt,yq, &
2854	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
2855	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
2856	IF (prt_level >=10) THEN
2857	print *,' Arg. de HBTM: yt2m ',yt2m
2858	print *,' Arg. de HBTM: yt10m ',yt10m
2859	print *,' Arg. de HBTM: yq2m ',yq2m
2860	print *,' Arg. de HBTM: yq10m ',yq10m
2861	print *,' Arg. de HBTM: yustar ',yustar
2862	print *,' Arg. de HBTM: y_flux_t ',y_flux_t
2863	print *,' Arg. de HBTM: y_flux_q ',y_flux_q
2864	print *,' Arg. de HBTM: yu ',yu
2865	print *,' Arg. de HBTM: yv ',yv
2866	print *,' Arg. de HBTM: yt ',yt
2867	print *,' Arg. de HBTM: yq ',yq
2868	ENDIF
2869	ELSE ! (iflag_split .eq.0)
2870	CALL HBTM(knon, ypaprs, ypplay, &
2871	yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
2872	y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
2873	ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
2874	ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
2875	IF (prt_level >=10) THEN
2876	print *,' Arg. de HBTM: yt2m_x ',yt2m_x
2877	print *,' Arg. de HBTM: yt10m_x ',yt10m_x
2878	print *,' Arg. de HBTM: yq2m_x ',yq2m_x
2879	print *,' Arg. de HBTM: yq10m_x ',yq10m_x
2880	print *,' Arg. de HBTM: yustar_x ',yustar_x
2881	print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x
2882	print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x
2883	print *,' Arg. de HBTM: yu_x ',yu_x
2884	print *,' Arg. de HBTM: yv_x ',yv_x
2885	print *,' Arg. de HBTM: yt_x ',yt_x
2886	print *,' Arg. de HBTM: yq_x ',yq_x
2887	ENDIF
2888	CALL HBTM(knon, ypaprs, ypplay, &
2889	yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
2890	y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
2891	ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
2892	ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
2893	!!!
2894	ENDIF ! (iflag_split .eq.0)
2895	!!!
2896
2897	!!! jyg le 07/02/2012
2898	IF (iflag_split .eq.0) THEN
2899	!!!
2900	DO j=1, knon
2901	i = ni(j)
2902	pblh(i,nsrf) = ypblh(j)
2903	wstar(i,nsrf) = ywstar(j)
2904	plcl(i,nsrf) = ylcl(j)
2905	capCL(i,nsrf) = ycapCL(j)
2906	oliqCL(i,nsrf) = yoliqCL(j)
2907	cteiCL(i,nsrf) = ycteiCL(j)
2908	pblT(i,nsrf) = ypblT(j)
2909	therm(i,nsrf) = ytherm(j)
2910	trmb1(i,nsrf) = ytrmb1(j)
2911	trmb2(i,nsrf) = ytrmb2(j)
2912	trmb3(i,nsrf) = ytrmb3(j)
2913	ENDDO
2914	IF (prt_level >=10) THEN
2915	print *, 'After HBTM: pblh ', pblh
2916	print *, 'After HBTM: plcl ', plcl
2917	print *, 'After HBTM: cteiCL ', cteiCL
2918	ENDIF
2919	ELSE !(iflag_split .eq.0)
2920	DO j=1, knon
2921	i = ni(j)
2922	pblh_x(i,nsrf) = ypblh_x(j)
2923	wstar_x(i,nsrf) = ywstar_x(j)
2924	plcl_x(i,nsrf) = ylcl_x(j)
2925	capCL_x(i,nsrf) = ycapCL_x(j)
2926	oliqCL_x(i,nsrf) = yoliqCL_x(j)
2927	cteiCL_x(i,nsrf) = ycteiCL_x(j)
2928	pblT_x(i,nsrf) = ypblT_x(j)
2929	therm_x(i,nsrf) = ytherm_x(j)
2930	trmb1_x(i,nsrf) = ytrmb1_x(j)
2931	trmb2_x(i,nsrf) = ytrmb2_x(j)
2932	trmb3_x(i,nsrf) = ytrmb3_x(j)
2933	ENDDO
2934	IF (prt_level >=10) THEN
2935	print *, 'After HBTM: pblh_x ', pblh_x
2936	print *, 'After HBTM: plcl_x ', plcl_x
2937	print *, 'After HBTM: cteiCL_x ', cteiCL_x
2938	ENDIF
2939	DO j=1, knon
2940	i = ni(j)
2941	pblh_w(i,nsrf) = ypblh_w(j)
2942	wstar_w(i,nsrf) = ywstar_w(j)
2943	plcl_w(i,nsrf) = ylcl_w(j)
2944	capCL_w(i,nsrf) = ycapCL_w(j)
2945	oliqCL_w(i,nsrf) = yoliqCL_w(j)
2946	cteiCL_w(i,nsrf) = ycteiCL_w(j)
2947	pblT_w(i,nsrf) = ypblT_w(j)
2948	therm_w(i,nsrf) = ytherm_w(j)
2949	trmb1_w(i,nsrf) = ytrmb1_w(j)
2950	trmb2_w(i,nsrf) = ytrmb2_w(j)
2951	trmb3_w(i,nsrf) = ytrmb3_w(j)
2952	ENDDO
2953	IF (prt_level >=10) THEN
2954	print *, 'After HBTM: pblh_w ', pblh_w
2955	print *, 'After HBTM: plcl_w ', plcl_w
2956	print *, 'After HBTM: cteiCL_w ', cteiCL_w
2957	ENDIF
2958	!!!
2959	ENDIF ! (iflag_split .eq.0)
2960	!!!
2961
2962	! print*,'OK pbl 6'
2963	#else
2964	! T2m not defined
2965	! No calculation
2966	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
2967	#endif
2968
2969	!****************************************************************************************
2970	! 15) End of loop over different surfaces
2971	!
2972	!****************************************************************************************
2973	ENDDO loop_nbsrf
2974
2975	!****************************************************************************************
2976	! 16) Calculate the mean value over all sub-surfaces for some variables
2977	!
2978	!****************************************************************************************
2979
2980	z0m(:,nbsrf+1) = 0.0
2981	z0h(:,nbsrf+1) = 0.0
2982	DO nsrf = 1, nbsrf
2983	DO i = 1, klon
2984	z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)
2985	z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)
2986	ENDDO
2987	ENDDO
2988
2989	! print*,'OK pbl 7'
2990	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
2991	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
2992	zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
2993	zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
2994	zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
2995	zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
2996
2997	!!! jyg le 07/02/2012
2998	IF (iflag_split .ge.1) THEN
2999	!!!
3000	!!! nrlmd & jyg les 02/05/2011, 05/02/2012
3001
3002	DO nsrf = 1, nbsrf
3003	DO k = 1, klev
3004	DO i = 1, klon
3005	zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
3006	zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
3007	zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
3008	zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
3009	!
3010	zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
3011	zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
3012	zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
3013	zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
3014	ENDDO
3015	ENDDO
3016	ENDDO
3017
3018	DO i = 1, klon
3019	zxsens_x(i) = - zxfluxt_x(i,1)
3020	zxsens_w(i) = - zxfluxt_w(i,1)
3021	ENDDO
3022	!!!
3023	ENDIF ! (iflag_split .ge.1)
3024	!!!
3025
3026	DO nsrf = 1, nbsrf
3027	DO k = 1, klev
3028	DO i = 1, klon
3029	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
3030	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
3031	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
3032	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
3033	ENDDO
3034	ENDDO
3035	ENDDO
3036
3037	DO i = 1, klon
3038	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
3039	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
3040	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
3041	ENDDO
3042	!!!
3043
3044	!
3045	! Incrementer la temperature du sol
3046	!
3047	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
3048	zt2m(:) = 0.0 ; zq2m(:) = 0.0 ; zn2mout(:,:) = 0
3049	zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
3050	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
3051	!!! jyg le 07/02/2012
3052	s_pblh_x(:) = 0.0 ; s_plcl_x(:) = 0.0
3053	s_pblh_w(:) = 0.0 ; s_plcl_w(:) = 0.0
3054	!!!
3055	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
3056	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
3057	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
3058	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
3059	wstar(:,is_ave)=0.
3060
3061	! print*,'OK pbl 9'
3062
3063	!!! nrlmd le 02/05/2011
3064	zxfluxlat_x(:) = 0.0 ; zxfluxlat_w(:) = 0.0
3065	!!!
3066
3067	DO nsrf = 1, nbsrf
3068	DO i = 1, klon
3069	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
3070
3071	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
3072	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
3073	wfbilo(i,nsrf) = (evap(i,nsrf)-(rain_f(i)+snow_f(i)))*pctsrf(i,nsrf)
3074	wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)
3075	wfrain(i,nsrf) = rain_f(i)*pctsrf(i,nsrf)
3076	wfsnow(i,nsrf) = snow_f(i)*pctsrf(i,nsrf)
3077
3078	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
3079	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
3080	ENDDO
3081	ENDDO
3082	!
3083	!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
3084	IF (iflag_order2_sollw == 1) THEN
3085	meansqT(:) = 0. ! as working buffer
3086	DO nsrf = 1, nbsrf
3087	DO i = 1, klon
3088	meansqT(i) = meansqT(i)+(ts(i,nsrf)-zxtsol(i))*2 pctsrf(i,nsrf)
3089	ENDDO
3090	ENDDO
3091	zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)
3092	ENDIF ! iflag_order2_sollw == 1
3093	!>al1
3094
3095	!!! jyg le 07/02/2012
3096	IF (iflag_split .eq.0) THEN
3097	DO nsrf = 1, nbsrf
3098	DO i = 1, klon
3099	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
3100	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
3101	!
3102	DO k = 1, 6
3103	zn2mout(i,k) = zn2mout(i,k) + n2mout(i,nsrf,k) * pctsrf(i,nsrf)
3104	ENDDO
3105	!
3106	zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
3107	wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
3108	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
3109	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
3110
3111	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
3112	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
3113	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
3114	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
3115	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
3116	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
3117	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
3118	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
3119	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
3120	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
3121	ENDDO
3122	ENDDO
3123	ELSE !(iflag_split .eq.0)
3124	DO nsrf = 1, nbsrf
3125	DO i = 1, klon
3126	!!! nrlmd le 02/05/2011
3127	zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
3128	zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
3129	!!!
3130	!!! jyg le 08/02/2012
3131	!! Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;
3132	!! pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
3133	!! pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
3134	!! pour les autres variables, on sort les valeurs de la region (x).
3135	zt2m(i) = zt2m(i) + (t2m_x(i,nsrf)+wake_s(i)(t2m_w(i,nsrf)-t2m_x(i,nsrf))) pctsrf(i,nsrf)
3136	zq2m(i) = zq2m(i) + q2m_x(i,nsrf) * pctsrf(i,nsrf)
3137	!
3138	DO k = 1, 6
3139	zn2mout(i,k) = zn2mout(i,k) + n2mout_x(i,nsrf,k) * pctsrf(i,nsrf)
3140	ENDDO
3141	!
3142	zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
3143	wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
3144	zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
3145	zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
3146	!
3147	s_pblh(i) = s_pblh(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
3148	s_pblh_x(i) = s_pblh_x(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
3149	s_pblh_w(i) = s_pblh_w(i) + pblh_w(i,nsrf) * pctsrf(i,nsrf)
3150	!
3151	s_plcl(i) = s_plcl(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
3152	s_plcl_x(i) = s_plcl_x(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
3153	s_plcl_w(i) = s_plcl_w(i) + plcl_w(i,nsrf) * pctsrf(i,nsrf)
3154	!
3155	s_capCL(i) = s_capCL(i) + capCL_x(i,nsrf) * pctsrf(i,nsrf)
3156	s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
3157	s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
3158	s_pblT(i) = s_pblT(i) + pblT_x(i,nsrf) * pctsrf(i,nsrf)
3159	s_therm(i) = s_therm(i) + therm_x(i,nsrf) * pctsrf(i,nsrf)
3160	s_trmb1(i) = s_trmb1(i) + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
3161	s_trmb2(i) = s_trmb2(i) + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
3162	s_trmb3(i) = s_trmb3(i) + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
3163	ENDDO
3164	ENDDO
3165	DO i = 1, klon
3166	qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
3167	ENDDO
3168	!!!
3169	ENDIF ! (iflag_split .eq.0)
3170	!!!
3171
3172	IF (check) THEN
3173	amn=MIN(ts(1,is_ter),1000.)
3174	amx=MAX(ts(1,is_ter),-1000.)
3175	DO i=2, klon
3176	amn=MIN(ts(i,is_ter),amn)
3177	amx=MAX(ts(i,is_ter),amx)
3178	ENDDO
3179	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
3180	ENDIF
3181
3182	!jg ?
3183	!!$!
3184	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
3185	!!$! sub-surfaces is distributed.
3186	!!$!
3187	!!$ DO nsrf = 1, nbsrf
3188	!!$ DO i = 1, klon
3189	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
3190	!!$ ts(i,nsrf) = zxtsol(i)
3191	!!$ t2m(i,nsrf) = zt2m(i)
3192	!!$ q2m(i,nsrf) = zq2m(i)
3193	!!$ u10m(i,nsrf) = zu10m(i)
3194	!!$ v10m(i,nsrf) = zv10m(i)
3195	!!$
3196	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
3197	!!$ pblh(i,nsrf) = s_pblh(i)
3198	!!$ plcl(i,nsrf) = s_plcl(i)
3199	!!$ capCL(i,nsrf) = s_capCL(i)
3200	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
3201	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
3202	!!$ pblT(i,nsrf) = s_pblT(i)
3203	!!$ therm(i,nsrf) = s_therm(i)
3204	!!$ trmb1(i,nsrf) = s_trmb1(i)
3205	!!$ trmb2(i,nsrf) = s_trmb2(i)
3206	!!$ trmb3(i,nsrf) = s_trmb3(i)
3207	!!$ ENDIF
3208	!!$ ENDDO
3209	!!$ ENDDO
3210
3211
3212	DO i = 1, klon
3213	fder(i) = - 4.0RSIGMAzxtsol(i)**3
3214	ENDDO
3215
3216	zxqsurf(:) = 0.0
3217	zxsnow(:) = 0.0
3218	DO nsrf = 1, nbsrf
3219	DO i = 1, klon
3220	zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
3221	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
3222	ENDDO
3223	ENDDO
3224
3225	! Premier niveau de vent sortie dans physiq.F
3226	zu1(:) = u(:,1)
3227	zv1(:) = v(:,1)
3228
3229	END SUBROUTINE pbl_surface
3230	!
3231	!****************************************************************************************
3232	!
3233	SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
3234
3235	USE indice_sol_mod
3236
3237	INCLUDE "dimsoil.h"
3238
3239	! Ouput variables
3240	!****************************************************************************************
3241	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
3242	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
3243	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
3244	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
3245
3246
3247	!****************************************************************************************
3248	! Return module variables for writing to restart file
3249	!
3250	!****************************************************************************************
3251	fder_rst(:) = fder(:)
3252	snow_rst(:,:) = snow(:,:)
3253	qsurf_rst(:,:) = qsurf(:,:)
3254	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
3255
3256	!****************************************************************************************
3257	! Deallocate module variables
3258	!
3259	!****************************************************************************************
3260	! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
3261	IF (ALLOCATED(fder)) DEALLOCATE(fder)
3262	IF (ALLOCATED(snow)) DEALLOCATE(snow)
3263	IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
3264	IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
3265
3266	!jyg<
3267	!****************************************************************************************
3268	! Deallocate variables for pbl splitting
3269	!
3270	!****************************************************************************************
3271
3272	CALL wx_pbl_final
3273	!>jyg
3274
3275	END SUBROUTINE pbl_surface_final
3276	!
3277	!****************************************************************************************
3278	!
3279
3280	!albedo SB >>>
3281	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
3282	evap, z0m, z0h, agesno, &
3283	tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)
3284	!albedo SB <<<
3285	! Give default values where new fraction has appread
3286
3287	USE indice_sol_mod
3288	use phys_state_var_mod, only: delta_sal, ds_ns, dt_ns, delta_sst
3289	use config_ocean_skin_m, only: activate_ocean_skin
3290
3291	INCLUDE "dimsoil.h"
3292	INCLUDE "clesphys.h"
3293	INCLUDE "compbl.h"
3294
3295	! Input variables
3296	!****************************************************************************************
3297	INTEGER, INTENT(IN) :: itime
3298	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
3299
3300	! InOutput variables
3301	!****************************************************************************************
3302	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
3303	!albedo SB >>>
3304	REAL, DIMENSION(klon,nsw,nbsrf), INTENT(INOUT) :: alb_dir, alb_dif
3305	INTEGER :: k
3306	!albedo SB <<<
3307	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m
3308	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: evap, agesno
3309	REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT) :: z0m,z0h
3310	REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
3311
3312	! Local variables
3313	!****************************************************************************************
3314	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
3315	CHARACTER(len=80) :: abort_message
3316	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
3317	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
3318	!
3319	! All at once !!
3320	!****************************************************************************************
3321
3322	DO nsrf = 1, nbsrf
3323	! First decide complement sub-surfaces
3324	SELECT CASE (nsrf)
3325	CASE(is_oce)
3326	nsrf_comp1=is_sic
3327	nsrf_comp2=is_ter
3328	nsrf_comp3=is_lic
3329	CASE(is_sic)
3330	nsrf_comp1=is_oce
3331	nsrf_comp2=is_ter
3332	nsrf_comp3=is_lic
3333	CASE(is_ter)
3334	nsrf_comp1=is_lic
3335	nsrf_comp2=is_oce
3336	nsrf_comp3=is_sic
3337	CASE(is_lic)
3338	nsrf_comp1=is_ter
3339	nsrf_comp2=is_oce
3340	nsrf_comp3=is_sic
3341	END SELECT
3342
3343	! Initialize all new fractions
3344	DO i=1, klon
3345	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
3346
3347	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
3348	! Use the complement sub-surface, keeping the continents unchanged
3349	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
3350	evap(i,nsrf) = evap(i,nsrf_comp1)
3351	z0m(i,nsrf) = z0m(i,nsrf_comp1)
3352	z0h(i,nsrf) = z0h(i,nsrf_comp1)
3353	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
3354	!albedo SB >>>
3355	DO k=1,nsw
3356	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp1)
3357	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp1)
3358	ENDDO
3359	!albedo SB <<<
3360	ustar(i,nsrf) = ustar(i,nsrf_comp1)
3361	u10m(i,nsrf) = u10m(i,nsrf_comp1)
3362	v10m(i,nsrf) = v10m(i,nsrf_comp1)
3363	IF (iflag_pbl > 1) THEN
3364	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
3365	ENDIF
3366	mfois(nsrf) = mfois(nsrf) + 1
3367	! F. Codron sensible default values for ocean and sea ice
3368	IF (nsrf.EQ.is_oce) THEN
3369	tsurf(i,nsrf) = 271.35
3370	! (temperature of sea water under sea ice, so that
3371	! is also the temperature of appearing sea water)
3372	DO k=1,nsw
3373	alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
3374	alb_dif(i,k,nsrf) = 0.06
3375	ENDDO
3376	if (activate_ocean_skin >= 1) then
3377	if (activate_ocean_skin == 2 &
3378	.and. type_ocean == "couple") then
3379	delta_sal(i) = 0.
3380	delta_sst(i) = 0.
3381	end if
3382
3383	ds_ns(i) = 0.
3384	dt_ns(i) = 0.
3385	end if
3386	ELSE IF (nsrf.EQ.is_sic) THEN
3387	tsurf(i,nsrf) = 271.35
3388	! (Temperature at base of sea ice. Surface
3389	! temperature could be higher, up to 0 Celsius
3390	! degrees. We set it to -1.8 Celsius degrees for
3391	! consistency with the ocean slab model.)
3392	DO k=1,nsw
3393	alb_dir(i,k,nsrf) = 0.3 ! thin ice
3394	alb_dif(i,k,nsrf) = 0.3
3395	ENDDO
3396	ENDIF
3397	ELSE
3398	! The continents have changed. The new fraction receives the mean sum of the existent fractions
3399	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3400	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3401	z0m(i,nsrf) = z0m(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0m(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3402	z0h(i,nsrf) = z0h(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0h(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3403	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3404	!albedo SB >>>
3405	DO k=1,nsw
3406	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3407	alb_dir(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3408	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3409	alb_dif(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3410	ENDDO
3411	!albedo SB <<<
3412	ustar(i,nsrf) = ustar(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3413	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3414	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3415	IF (iflag_pbl > 1) THEN
3416	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3417	ENDIF
3418
3419	! Security abort. This option has never been tested. To test, comment the following line.
3420	! abort_message='The fraction of the continents have changed!'
3421	! CALL abort_physic(modname,abort_message,1)
3422	nfois(nsrf) = nfois(nsrf) + 1
3423	ENDIF
3424	snow(i,nsrf) = 0.
3425	agesno(i,nsrf) = 0.
3426	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
3427	ELSE
3428	pfois(nsrf) = pfois(nsrf)+ 1
3429	ENDIF
3430	ENDDO
3431
3432	ENDDO
3433
3434	END SUBROUTINE pbl_surface_newfrac
3435	!
3436	!****************************************************************************************
3437	!
3438	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format