Context Navigation

pbl_surface_mod.F90 @ 4881

Last change on this file since 4881 was 4881, checked in by evignon, 7 weeks ago
extraction plus propre de la dissipation de TKE
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 154.0 KB

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

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

Download in other formats:

Original Format