Context Navigation

pbl_surface_mod.F90 @ 4660

Last change on this file since 4660 was 4653, checked in by evignon, 15 months ago
prise en compte de l'humidite pour le calcul du flux de flottabilite dans atke + petits ajustements
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.3 KB

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

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

Download in other formats:

Original Format