Context Navigation

pbl_surface_mod.F90 @ 3901

Last change on this file since 3901 was 3901, checked in by evignon, 3 years ago

Commission LMDZ-sisvat, deuxieme phase:

suppression du repertoire sisvat
nettoyage de certains scripts de la physique et des makefile

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: 146.5 KB

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

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

Download in other formats:

Original Format