Context Navigation

pbl_surface_mod.F90 @ 4553

Last change on this file since 4553 was 4545, checked in by evignon, 14 months ago
modifications suite au dernier atelier TKE. Ajout de la routine call_atke_mod et test d'une prediction explicit du cisaillement pour calculer les Km et Kh
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: 152.2 KB

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

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

Download in other formats:

Original Format