Context Navigation

pbl_surface_mod.F90 @ 3713

Last change on this file since 3713 was 3687, checked in by lguez, 5 years ago

Modify sensible heat due to rain sent to the ocean

Modify the sensible heat flux due to rain which is sent to the
ocean. Replace the computation of sens_prec_liq in procedure
calcul_fluxs by a call to sens_heat_rain. Set sens_prec_sol in
procedure calcul_fluxs to 0 because, for now, sens_heat_rain is
supposed to account for both rain and snow.

For the call to sens_heat_rain in procedure calcul_fluxs, we need
an additional dummy argument rhoa of calcul_fluxs. Add dummy
argument rhoa to ocean_cpl_noice, ocean_forced_noice,
ocean_forced_ice and ocean_cpl_ice because we need to pass it down
to calcul_fluxs.

Change the dimension of sens_prec_liq and sens_prec_sol in
procedures calcul_fluxs, ocean_cpl_noice, ocean_cpl_ice,
ocean_forced_noice, ocean_forced_ice, cpl_send_ocean_fields and
cpl_send_seaice_fields to knon.

In procedures ocean_forced_noice and ocean_cpl_noice, promote
sens_prec_liq from local variable to dummy argument because we need
it in surf_ocean. Remove useless initialization of sens_prec_liq
and sens_prec_sol in ocean_cpl_noice, ocean_cpl_ice,
ocean_forced_ice and ocean_forced_noice: they are intent out in
calcul_fluxs.

Remove variable rf of module phys_output_var_mod, we use
sens_prec_liq instead. Remove local variable yrf of procedure
pbl_surface. rf and yrf appeared in pbl_surface only to be output.
Remove variable o_rf of module phys_output_ctrlout_mod. Remove
dummy argument rf of procedure surf_ocean.

Do not call sens_heat_rain in surf_ocean since we now call it
from calcul_fluxs. Move the computation of rhoa in surf_ocean
before the calls to ocean_cpl_noice and ocean_forced_noice.

Add the computation of rhoa in surf_seaice, to pass it down to
ocean_cpl_ice and ocean_forced_ice.

If activate_ocean_skin == 1 then the results are changed because the
call to sens_heat_rain in calcul_fluxs now uses the SST from the
current time step of physics. On this point, the present revision
reverses revision [3463].

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 134.1 KB

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

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

Download in other formats:

Original Format