Context Navigation

pbl_surface_mod.F90 @ 4179

Last change on this file since 4179 was 3956, checked in by jyg, 3 years ago
Bug fixes concerning various variables ill-initialized, ill-used, ill-printed, or ill-placed. + cv_gen moved from phys_local_var_mod.F90 to phys_state_var_mod.F90; ==> changes in physiq_mod.F90 and phys_output_write.F90 + awake_dens added in phys_state_var_mod.F90 + cv_gen and awake_dens now initialized in phyetat0.F90 and written in phyredem.F90 + cv_gen, awake_dens, and solswfdiff now initialized in old_lmdz1d.F90 and scm.F90 + useless variables suppressed in pbl_surface_mod.F90.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 147.3 KB

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

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

Download in other formats:

Original Format