Context Navigation

pbl_surface_mod.F90 @ 4747

Last change on this file since 4747 was 4747, checked in by jyg, 6 months ago
smalle bug fix in pbl_surface_mod.F90
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: 153.5 KB

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

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

Download in other formats:

Original Format