Context Navigation

pbl_surface_mod.F90 @ 4462

Last change on this file since 4462 was 4449, checked in by evignon, 17 months ago
commission du nouveau schema de turbulence developpe dans le cadre de l'atelier tke
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: 148.7 KB

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

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

Download in other formats:

Original Format