Context Navigation

pbl_surface_mod.F90 @ 4096

Last change on this file since 4096 was 4020, checked in by lguez, 3 years ago

Send 3 more fields to the ocean

Send 3 more fields to the ocean to compute CO2 flux at
ocean-atmosphere interface. The three fields are dter and dser, which
already existed, and a newly created field: dt_ds. So dter and dser
have to become state variables. The variable dt_ds of module
phys_state_var_mod is only allocated and defined if
activate_ocean_skin == 2 and type_ocean == "couple".

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 147.8 KB

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

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

Download in other formats:

Original Format