Context Navigation

pbl_surface_mod.F90 @ 4283

Last change on this file since 4283 was 4283, checked in by jghattas, 19 months ago

Added landice_opt=2 : Treat continental land ice fractions in ORCHIDEE => pctsrf(:,is_lic) = 0.0 in LMDZ.

For this option, some more variables are needed from ORCHIDEE. Therfor change in the interface LMDZ-ORCHIDEE in surf_land_orchidee_mod is done. Previous interface is moved to surf_land_orchidee_nolic_mod.f90. To compile with previous interface, cpp key ORCHIDEE_NOLIC is added. Previous interface is compiled with argument orchidee2.1 in makelmdz and makelmdz_fcm.

At the same time, when the interface was changed, the variable yrmu0(coszang) was added in the call to intersurf_initialize_gathered. This is needed in ORCHIDEE to better initialize the model.

Modifications done by Etienne Vignon and Josefine Ghattas

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: 147.6 KB

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

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

Download in other formats:

Original Format