Context Navigation

pbl_surface_mod.F90 @ 4662

Last change on this file since 4662 was 4662, checked in by Laurent Fairhead, 10 months ago
Debugging stuff
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: 149.1 KB

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

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

Download in other formats:

Original Format