Context Navigation

pbl_surface_mod.F90 @ 5014

Last change on this file since 5014 was 5009, checked in by musat, 12 months ago
Correction bug : remplace wake_s sur la grille klon par ywake_s sur la grille knon.
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: 153.9 KB

Line
1	!
2	! $Id: pbl_surface_mod.F90 5009 2024-07-04 11:28:51Z abarral $
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	zgeo1(1:knon) = ywake_s(1:knon)zgeo1_w(1:knon) + (1.-ywake_s(1:knon))zgeo1_x(1:knon)
1667
1668	! --- special Dice. JYG+MPL 25112013 puis BOMEX
1669	IF (ok_prescr_ust) THEN
1670	DO i = 1, knon
1671	! print *,'ycdragm_w avant=',ycdragm_w(i)
1672	vent= sqrt(yu_w(i,1)yu_w(i,1)+yv_w(i,1)yv_w(i,1))
1673	ycdragm_w(i) = ust*ust/(1.+vent)/vent
1674	! print *,'ycdragm_w ust yu yv apres=',ycdragm_w(i),ust,yu_w(i,1),yv_w(i,1)
1675	ENDDO
1676	ENDIF
1677	IF (prt_level >=10) print *,'clcdrag -> ycdragh_w ', ycdragh_w(1:knon)
1678	!!!
1679	ENDIF ! (iflag_split .eq.0)
1680	!!!
1681
1682
1683	!****************************************************************************************
1684	! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm.
1685	!
1686	!****************************************************************************************
1687
1688	!!! jyg le 07/02/2012
1689	IF (iflag_split .eq.0) THEN
1690	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1691	IF (prt_level >=10) THEN
1692	print *,' args coef_diff_turb: yu ', yu(1:knon,:)
1693	print *,' args coef_diff_turb: yv ', yv(1:knon,:)
1694	print *,' args coef_diff_turb: yq ', yq(1:knon,:)
1695	print *,' args coef_diff_turb: yt ', yt(1:knon,:)
1696	print *,' args coef_diff_turb: yts ', yts(1:knon)
1697	print *,' args coef_diff_turb: yz0m ', yz0m(1:knon)
1698	print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon)
1699	print *,' args coef_diff_turb: ycdragm ', ycdragm(1:knon)
1700	print *,' args coef_diff_turb: ycdragh ', ycdragh(1:knon)
1701	print *,' args coef_diff_turb: ytke ', ytke(1:knon,:)
1702	ENDIF
1703
1704	IF (iflag_pbl>=50) THEN
1705	CALL call_atke(dtime,knon,klev,ycdragm(1:knon), ycdragh(1:knon),yus0(1:knon),yvs0(1:knon),yts(1:knon), &
1706	yu(1:knon,:),yv(1:knon,:),yt(1:knon,:),yq(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), &
1707	ytke(1:knon,:),yeps(1:knon,:), ycoefm(1:knon,:), ycoefh(1:knon,:))
1708
1709	ELSE
1710
1711	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
1712	ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
1713	ycoefm, ycoefh, ytke, yeps, y_treedrg)
1714	! ycoefm, ycoefh, ytke)
1715	!FC y_treedrg ajoute
1716	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1717	! In this case, coef_diff_turb is called for the Cd only
1718	DO k = 2, klev
1719	DO j = 1, knon
1720	i = ni(j)
1721	ycoefh(j,k) = zcoefh(i,k,nsrf)
1722	ycoefm(j,k) = zcoefm(i,k,nsrf)
1723	ENDDO
1724	ENDDO
1725	ENDIF
1726
1727	ENDIF ! iflag_pbl >= 50
1728
1729	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh(1:knon,:)
1730
1731
1732	ELSE !(iflag_split .eq.0)
1733
1734
1735	IF (prt_level >=10) THEN
1736	print *,' args coef_diff_turb: yu_x ', yu_x(1:knon,:)
1737	print *,' args coef_diff_turb: yv_x ', yv_x(1:knon,:)
1738	print *,' args coef_diff_turb: yq_x ', yq_x(1:knon,:)
1739	print *,' args coef_diff_turb: yt_x ', yt_x(1:knon,:)
1740	print *,' args coef_diff_turb: yts_x ', yts_x(1:knon)
1741	print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon)
1742	print *,' args coef_diff_turb: ycdragm_x ', ycdragm_x(1:knon)
1743	print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x(1:knon)
1744	print *,' args coef_diff_turb: ytke_x ', ytke_x(1:knon,:)
1745	ENDIF
1746
1747
1748	IF (iflag_pbl>=50) THEN
1749
1750	CALL call_atke(dtime,knon,klev,ycdragm_x(1:knon),ycdragh_x(1:knon),yus0(1:knon),yvs0(1:knon),yts_x(1:knon), &
1751	yu_x(1:knon,:),yv_x(1:knon,:),yt_x(1:knon,:),yq_x(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), &
1752	ytke_x(1:knon,:),yeps_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:))
1753
1754	ELSE
1755
1756	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
1757	ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
1758	ycoefm_x, ycoefh_x, ytke_x,yeps_x,y_treedrg)
1759	! ycoefm_x, ycoefh_x, ytke_x)
1760	!FC doit on le mettre ( on ne l utilise pas si il y a du spliting)
1761	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1762	! In this case, coef_diff_turb is called for the Cd only
1763	DO k = 2, klev
1764	DO j = 1, knon
1765	i = ni(j)
1766	ycoefh_x(j,k) = zcoefh(i,k,nsrf)
1767	ycoefm_x(j,k) = zcoefm(i,k,nsrf)
1768	ENDDO
1769	ENDDO
1770	ENDIF
1771
1772	ENDIF ! iflag_pbl >= 50
1773
1774	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x(1:knon,:)
1775	!
1776	IF (prt_level >=10) THEN
1777	print *,' args coef_diff_turb: yu_w ', yu_w(1:knon,:)
1778	print *,' args coef_diff_turb: yv_w ', yv_w(1:knon,:)
1779	print *,' args coef_diff_turb: yq_w ', yq_w(1:knon,:)
1780	print *,' args coef_diff_turb: yt_w ', yt_w(1:knon,:)
1781	print *,' args coef_diff_turb: yts_w ', yts_w(1:knon)
1782	print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon)
1783	print *,' args coef_diff_turb: ycdragm_w ', ycdragm_w(1:knon)
1784	print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w(1:knon)
1785	print *,' args coef_diff_turb: ytke_w ', ytke_w(1:knon,:)
1786	ENDIF
1787
1788	IF (iflag_pbl>=50) THEN
1789
1790	CALL call_atke(dtime,knon,klev,ycdragm_w(1:knon),ycdragh_w(1:knon),yus0(1:knon),yvs0(1:knon),yts_w(1:knon), &
1791	yu_w(1:knon,:),yv_w(1:knon,:),yt_w(1:knon,:),yq_w(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), &
1792	ytke_w(1:knon,:),yeps_w(1:knon,:),ycoefm_w(1:knon,:),ycoefh_w(1:knon,:))
1793
1794	ELSE
1795
1796	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
1797	ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
1798	ycoefm_w, ycoefh_w, ytke_w,yeps_w,y_treedrg)
1799	! ycoefm_w, ycoefh_w, ytke_w)
1800	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1801	! In this case, coef_diff_turb is called for the Cd only
1802	DO k = 2, klev
1803	DO j = 1, knon
1804	i = ni(j)
1805	ycoefh_w(j,k) = zcoefh(i,k,nsrf)
1806	ycoefm_w(j,k) = zcoefm(i,k,nsrf)
1807	ENDDO
1808	ENDDO
1809	ENDIF
1810
1811	ENDIF ! iflag_pbl >= 50
1812
1813
1814	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w(1:knon,:)
1815
1816	!!!jyg le 10/04/2013
1817	!! En attendant de traiter le transport des traceurs dans les poches froides, formule
1818	!! arbitraire pour ycoefh et ycoefm
1819	DO k = 2,klev
1820	DO j = 1,knon
1821	ycoefh(j,k) = ycoefh_x(j,k) + ywake_s(j)*(ycoefh_w(j,k) - ycoefh_x(j,k))
1822	ycoefm(j,k) = ycoefm_x(j,k) + ywake_s(j)*(ycoefm_w(j,k) - ycoefm_x(j,k))
1823	ENDDO
1824	ENDDO
1825
1826
1827	ENDIF ! (iflag_split .eq.0)
1828
1829
1830	!****************************************************************************************
1831	!
1832	! 8) "La descente" - "The downhill"
1833	!
1834	! climb_hq_down and climb_wind_down calculate the coefficients
1835	! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
1836	! Only the coefficients at surface for H and Q are returned.
1837	!
1838	!****************************************************************************************
1839
1840	! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
1841	!!! jyg le 07/02/2012
1842	IF (iflag_split .eq.0) THEN
1843	!!!
1844	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1845	CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
1846	ydelp, yt, yq, dtime, &
1847	!!! jyg le 09/05/2011
1848	CcoefH, CcoefQ, DcoefH, DcoefQ, &
1849	Kcoef_hq, gama_q, gama_h, &
1850	!!!
1851	AcoefH, AcoefQ, BcoefH, BcoefQ)
1852	ELSE !(iflag_split .eq.0)
1853	CALL climb_hq_down(knon, ycoefh_x, ypaprs, ypplay, &
1854	ydelp, yt_x, yq_x, dtime, &
1855	!!! nrlmd le 02/05/2011
1856	CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
1857	Kcoef_hq_x, gama_q_x, gama_h_x, &
1858	!!!
1859	AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x)
1860	!!!
1861	IF (prt_level >=10) THEN
1862	PRINT *,'pbl_surface (climb_hq_down.x->) AcoefH_x ',AcoefH_x
1863	PRINT *,'pbl_surface (climb_hq_down.x->) AcoefQ_x ',AcoefQ_x
1864	PRINT *,'pbl_surface (climb_hq_down.x->) BcoefH_x ',BcoefH_x
1865	PRINT *,'pbl_surface (climb_hq_down.x->) BcoefQ_x ',BcoefQ_x
1866	ENDIF
1867	!
1868	CALL climb_hq_down(knon, ycoefh_w, ypaprs, ypplay, &
1869	ydelp, yt_w, yq_w, dtime, &
1870	!!! nrlmd le 02/05/2011
1871	CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
1872	Kcoef_hq_w, gama_q_w, gama_h_w, &
1873	!!!
1874	AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w)
1875	!!!
1876	IF (prt_level >=10) THEN
1877	PRINT *,'pbl_surface (climb_hq_down.w->) AcoefH_w ',AcoefH_w
1878	PRINT *,'pbl_surface (climb_hq_down.w->) AcoefQ_w ',AcoefQ_w
1879	PRINT *,'pbl_surface (climb_hq_down.w->) BcoefH_w ',BcoefH_w
1880	PRINT *,'pbl_surface (climb_hq_down.w->) BcoefQ_w ',BcoefQ_w
1881	ENDIF
1882	!!!
1883	ENDIF ! (iflag_split .eq.0)
1884	!!!
1885
1886	! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
1887	!!! jyg le 07/02/2012
1888	IF (iflag_split .eq.0) THEN
1889	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1890	CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
1891	!!! jyg le 09/05/2011
1892	CcoefU, CcoefV, DcoefU, DcoefV, &
1893	Kcoef_m, alf_1, alf_2, &
1894	!!!
1895	AcoefU, AcoefV, BcoefU, BcoefV)
1896	ELSE ! (iflag_split .eq.0)
1897	CALL climb_wind_down(knon, dtime, ycoefm_x, ypplay, ypaprs, yt_x, ydelp, yu_x, yv_x, &
1898	!!! nrlmd le 02/05/2011
1899	CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
1900	Kcoef_m_x, alf_1_x, alf_2_x, &
1901	!!!
1902	AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x)
1903	!
1904	CALL climb_wind_down(knon, dtime, ycoefm_w, ypplay, ypaprs, yt_w, ydelp, yu_w, yv_w, &
1905	!!! nrlmd le 02/05/2011
1906	CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
1907	Kcoef_m_w, alf_1_w, alf_2_w, &
1908	!!!
1909	AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w)
1910	!!!
1911	ENDIF ! (iflag_split .eq.0)
1912	!!!
1913
1914	! For blowing snow:
1915	IF (ok_bs) THEN
1916	! following Bintanja et al 2000, part II and Vionnet V PhD thesis
1917	! we assume that the eddy diffsivity coefficient for
1918	! suspended particles is a fraction of Kh
1919	do k=1,klev
1920	do j=1,knon
1921	ycoefqbs(j,k)=ycoefh(j,k)*zeta_bs
1922	enddo
1923	enddo
1924	CALL climb_qbs_down(knon, ycoefqbs, ypaprs, ypplay, &
1925	ydelp, yt, yqbs, dtime, &
1926	CcoefQBS, DcoefQBS, &
1927	Kcoef_qbs, gama_qbs, &
1928	AcoefQBS, BcoefQBS)
1929	ENDIF
1930
1931	!****************************************************************************************
1932	! 9) Small calculations
1933	!
1934	!****************************************************************************************
1935
1936	! - Reference pressure is given the values at surface level
1937	ypsref(:) = ypaprs(:,1)
1938
1939	! - CO2 field on 2D grid to be sent to ORCHIDEE
1940	! Transform to compressed field
1941	IF (carbon_cycle_cpl) THEN
1942	DO i=1,knon
1943	r_co2_ppm(i) = co2_send(ni(i))
1944	ENDDO
1945	ELSE
1946	r_co2_ppm(:) = co2_ppm ! Constant field
1947	ENDIF
1948
1949	!!! nrlmd le 02/05/2011 -----------------------On raccorde les 2 colonnes dans la couche 1
1950	!----------------------------------------------------------------------------------------
1951	!!! jyg le 07/02/2012
1952	!!! jyg le 01/02/2017
1953	IF (iflag_split .eq. 0) THEN
1954	yt1(:) = yt(:,1)
1955	yq1(:) = yq(:,1)
1956	ELSE IF (iflag_split .ge. 1) THEN
1957	!
1958	! Cdragq computation
1959	! ------------------
1960	!******************************************************************************
1961	! Cdragq computed from cdrag
1962	! The difference comes only from a factor (f_z0qh_oce) on z0, so that
1963	! it can be computed inside wx_pbl0_merge
1964	! More complicated appraches may require the propagation through
1965	! pbl_surface of an independant cdragq variable.
1966	!******************************************************************************
1967	!
1968	IF ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) THEN
1969	! Si on suit les formulations par exemple de Tessel, on
1970	! a z0h=0.4nu/u, z0q=0.62nu/u, d'ou f_z0qh_oce=0.62/0.4=1.55
1971	!! ycdragq_x(1:knon)=ycdragh_x(1:knon)* &
1972	!! log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
1973	!! ycdragq_w(1:knon)=ycdragh_w(1:knon)* &
1974	!! log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
1975	!
1976	DO j = 1,knon
1977	z1lay = zgeo1(j)/RG
1978	fact_cdrag = log(z1lay/yz0h(j))/log(z1lay/(f_z0qh_oce*yz0h(j)))
1979	ycdragq_x(j)=ycdragh_x(j)*fact_cdrag
1980	ycdragq_w(j)=ycdragh_w(j)*fact_cdrag
1981	!! Print *,'YYYYpbl0: fact_cdrag ', fact_cdrag
1982	ENDDO ! j = 1,knon
1983	!
1984	!! Print *,'YYYYpbl0: z1lay, yz0h, f_z0qh_oce, ycdragh_w, ycdragq_w ', &
1985	!! z1lay, yz0h(1:knon), f_z0qh_oce, ycdragh_w(1:knon), ycdragq_w(1:knon)
1986	ELSE
1987	ycdragq_x(1:knon)=ycdragh_x(1:knon)
1988	ycdragq_w(1:knon)=ycdragh_w(1:knon)
1989	ENDIF ! ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce)
1990	!
1991	CALL wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, ywake_s, &
1992	yts, y_delta_tsurf, ygustiness, &
1993	yt_x, yt_w, yq_x, yq_w, &
1994	yu_x, yu_w, yv_x, yv_w, &
1995	ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
1996	ycdragm_x, ycdragm_w, &
1997	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
1998	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
1999	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
2000	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
2001	Kech_h_x, Kech_h_w, Kech_h &
2002	)
2003	CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, &
2004	BcoefQ_x, BcoefQ_w &
2005	)
2006	CALL wx_pbl0_merge(knon, ypplay, ypaprs, &
2007	ywake_s, ydTs0, ydqs0, &
2008	yt_x, yt_w, yq_x, yq_w, &
2009	yu_x, yu_w, yv_x, yv_w, &
2010	ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
2011	ycdragm_x, ycdragm_w, &
2012	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
2013	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
2014	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
2015	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
2016	AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
2017	BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
2018	ycdragh, ycdragq, ycdragm, &
2019	yt1, yq1, yu1, yv1 &
2020	)
2021	IF (iflag_split .eq. 2 .AND. nsrf .ne. is_oce) THEN
2022	CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
2023	ywake_s, ybeta, ywake_cstar, ywake_dens, &
2024	AcoefH_x, AcoefH_w, &
2025	BcoefH_x, BcoefH_w, &
2026	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2027	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2028	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2029	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2030	yg_T, yg_Q, &
2031	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2032	ydTs_ins, ydqs_ins &
2033	)
2034	ELSE !
2035	AcoefH(:) = AcoefH_0(:)
2036	AcoefQ(:) = AcoefQ_0(:)
2037	BcoefH(:) = BcoefH_0(:)
2038	BcoefQ(:) = BcoefQ_0(:)
2039	yg_T(:) = 0.
2040	yg_Q(:) = 0.
2041	yGamma_dTs_phiT(:) = 0.
2042	yGamma_dQs_phiQ(:) = 0.
2043	ydTs_ins(:) = 0.
2044	ydqs_ins(:) = 0.
2045	ENDIF ! (iflag_split .eq. 2)
2046	ENDIF ! (iflag_split .eq.0)
2047	!!!
2048	IF (prt_level >=10) THEN
2049	DO i = 1, min(1,knon)
2050	PRINT *,'pbl_surface (merge->): yt(1,:) ',yt(i,:)
2051	PRINT *,'pbl_surface (merge->): yq(1,:) ',yq(i,:)
2052	PRINT *,'pbl_surface (merge->): yu(1,:) ',yu(i,:)
2053	PRINT *,'pbl_surface (merge->): yv(1,:) ',yv(i,:)
2054	PRINT *,'pbl_surface (merge->): AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1) ', &
2055	AcoefH(i), AcoefQ(i), AcoefU(i), AcoefV(i)
2056	PRINT *,'pbl_surface (merge->): BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1) ', &
2057	BcoefH(i), BcoefQ(i), BcoefU(i), BcoefV(i)
2058	ENDDO
2059
2060	ENDIF
2061
2062	! Save initial value of z0h for use in evappot (z0h wiil be computed again in the surface models)
2063	yz0h_old(1:knon) = yz0h(1:knon)
2064	!
2065	!****************************************************************************************
2066	!
2067	! Calulate t2m and q2m for the case of calculation at land grid points
2068	! t2m and q2m are needed as input to ORCHIDEE
2069	!
2070	!****************************************************************************************
2071	IF (nsrf == is_ter) THEN
2072
2073	DO i = 1, knon
2074	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
2075	* (ypaprs(i,1)-ypplay(i,1))
2076	ENDDO
2077
2078	! Calculate the temperature et relative humidity at 2m and the wind at 10m
2079	IF (iflag_new_t2mq2m==1) THEN
2080	CALL stdlevvarn(klon, knon, is_ter, zxli, &
2081	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
2082	yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
2083	yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
2084	yn2mout(:, nsrf, :))
2085	ELSE
2086	CALL stdlevvar(klon, knon, is_ter, zxli, &
2087	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
2088	yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
2089	yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
2090	ENDIF
2091
2092	ENDIF
2093
2094	!****************************************************************************************
2095	!
2096	! 10) Switch according to current surface
2097	! It is necessary to start with the continental surfaces because the ocean
2098	! needs their run-off.
2099	!
2100	!****************************************************************************************
2101	SELECT CASE(nsrf)
2102
2103	CASE(is_ter)
2104	! print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)
2105	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
2106	rlon, rlat, yrmu0, &
2107	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
2108	!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
2109	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
2110	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2111	AcoefU, AcoefV, BcoefU, BcoefV, &
2112	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
2113	ylwdown, yq2m, yt2m, &
2114	ysnow, yqsol, yagesno, ytsoil, &
2115	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,yfluxbs,&
2116	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
2117	y_flux_u1, y_flux_v1, &
2118	yveget,ylai,yheight )
2119
2120	!FC quid qd yveget ylai yheight ne sont pas definit
2121	!FC yveget,ylai,yheight, &
2122	IF (ifl_pbltree .ge. 1) THEN
2123	CALL freinage(knon, yu, yv, yt, &
2124	! yveget,ylai, yheight,ypaprs,ypplay,y_d_u_frein,y_d_v_frein)
2125	yveget,ylai, yheight,ypaprs,ypplay,y_treedrg, y_d_u_frein,y_d_v_frein)
2126	ENDIF
2127
2128
2129	! Special DICE MPL 05082013 puis BOMEX
2130	IF (ok_prescr_ust) THEN
2131	DO j=1,knon
2132	! ysnow(:)=0.
2133	! yqsol(:)=0.
2134	! yagesno(:)=50.
2135	! ytsoil(:,:)=300.
2136	! yz0_new(:)=0.001
2137	! yevap(:)=flat/RLVTT
2138	! yfluxlat(:)=-flat
2139	! yfluxsens(:)=-fsens
2140	! yqsurf(:)=0.
2141	! ytsurf_new(:)=tg
2142	! y_dflux_t(:)=0.
2143	! y_dflux_q(:)=0.
2144	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)
2145	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)
2146	ENDDO
2147	ENDIF
2148
2149	CASE(is_lic)
2150	! Martin
2151
2152	IF (landice_opt .LT. 2) THEN
2153	! Land ice is treated by LMDZ and not by ORCHIDEE
2154	CALL surf_landice(itap, dtime, knon, ni, &
2155	rlon, rlat, debut, lafin, &
2156	yrmu0, ylwdown, yalb, zgeo1, &
2157	ysolsw, ysollw, yts, ypplay(:,1), &
2158	ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
2159	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2160	AcoefU, AcoefV, BcoefU, BcoefV, &
2161	AcoefQBS, BcoefQBS, &
2162	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
2163	ysnow, yqsurf, yqsol,yqbs1, yagesno, &
2164	ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
2165	yfluxbs, ytsurf_new, y_dflux_t, y_dflux_q, &
2166	yzmea, yzsig, ycldt, &
2167	ysnowhgt, yqsnow, ytoice, ysissnow, &
2168	yalb3_new, yrunoff, &
2169	y_flux_u1, y_flux_v1)
2170
2171	!jyg<
2172	!! alb3_lic(:)=0.
2173	!>jyg
2174	DO j = 1, knon
2175	i = ni(j)
2176	alb3_lic(i) = yalb3_new(j)
2177	snowhgt(i) = ysnowhgt(j)
2178	qsnow(i) = yqsnow(j)
2179	to_ice(i) = ytoice(j)
2180	sissnow(i) = ysissnow(j)
2181	runoff(i) = yrunoff(j)
2182	ENDDO
2183	! Martin
2184	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2185	IF (ok_prescr_ust) THEN
2186	DO j=1,knon
2187	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)
2188	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)
2189	ENDDO
2190	ENDIF
2191
2192	END IF
2193
2194	CASE(is_oce)
2195	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
2196	ywindsp, rmu0, yfder, yts, &
2197	itap, dtime, jour, knon, ni, &
2198	!!jyg ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
2199	ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt(:,1), yq(:,1),& ! ym missing init
2200	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2201	AcoefU, AcoefV, BcoefU, BcoefV, &
2202	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
2203	ysnow, yqsurf, yagesno, &
2204	yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
2205	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
2206	y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
2207	yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
2208	ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss)
2209	IF (prt_level >=10) THEN
2210	print *,'arg de surf_ocean: ycdragh ',ycdragh(1:knon)
2211	print *,'arg de surf_ocean: ycdragm ',ycdragm(1:knon)
2212	print *,'arg de surf_ocean: yt ', yt(1:knon,:)
2213	print *,'arg de surf_ocean: yq ', yq(1:knon,:)
2214	print *,'arg de surf_ocean: yts ', yts(1:knon)
2215	print *,'arg de surf_ocean: AcoefH ',AcoefH(1:knon)
2216	print *,'arg de surf_ocean: AcoefQ ',AcoefQ(1:knon)
2217	print *,'arg de surf_ocean: BcoefH ',BcoefH(1:knon)
2218	print *,'arg de surf_ocean: BcoefQ ',BcoefQ(1:knon)
2219	print *,'arg de surf_ocean: yevap ',yevap(1:knon)
2220	print *,'arg de surf_ocean: yfluxsens ',yfluxsens(1:knon)
2221	print *,'arg de surf_ocean: yfluxlat ',yfluxlat(1:knon)
2222	print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new(1:knon)
2223	ENDIF
2224	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2225	IF (ok_prescr_ust) THEN
2226	DO j=1,knon
2227	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)
2228	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)
2229	ENDDO
2230	ENDIF
2231
2232	CASE(is_sic)
2233	CALL surf_seaice( &
2234	!albedo SB >>>
2235	rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, &
2236	!albedo SB <<<
2237	itap, dtime, jour, knon, ni, &
2238	lafin, &
2239	!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
2240	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
2241	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2242	AcoefU, AcoefV, BcoefU, BcoefV, &
2243	ypsref, yu1, yv1, ygustiness, pctsrf, &
2244	ysnow, yqsurf, yqsol, yagesno, ytsoil, &
2245	!albedo SB >>>
2246	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
2247	!albedo SB <<<
2248	ytsurf_new, y_dflux_t, y_dflux_q, &
2249	y_flux_u1, y_flux_v1)
2250
2251	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2252	IF (ok_prescr_ust) THEN
2253	DO j=1,knon
2254	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)
2255	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)
2256	ENDDO
2257	ENDIF
2258
2259	CASE DEFAULT
2260	WRITE(lunout,*) 'Surface index = ', nsrf
2261	abort_message = 'Surface index not valid'
2262	CALL abort_physic(modname,abort_message,1)
2263	END SELECT
2264
2265
2266	!****************************************************************************************
2267	! 11) - Calcul the increment of surface temperature
2268	!
2269	!****************************************************************************************
2270
2271	IF (evap0>=0.) THEN
2272	yevap(1:knon)=evap0
2273	yevap(1:knon)=RLVTT*evap0
2274	ENDIF
2275
2276	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
2277
2278	!****************************************************************************************
2279	!
2280	! 12) "La remontee" - "The uphill"
2281	!
2282	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
2283	! for X=H, Q, U and V, for all vertical levels.
2284	!
2285	!****************************************************************************************
2286	!!
2287	!!!
2288	!!! jyg le 10/04/2013 et EV 10/2020
2289
2290	IF (ok_forc_tsurf) THEN
2291	DO j=1,knon
2292	ytsurf_new(j)=tg
2293	y_d_ts(j) = ytsurf_new(j) - yts(j)
2294	ENDDO
2295	ENDIF ! ok_forc_tsurf
2296
2297	!!!
2298	IF (ok_flux_surf) THEN
2299	IF (prt_level >=10) THEN
2300	PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
2301	ENDIF
2302	y_flux_t1(:) = fsens
2303	y_flux_q1(:) = flat/RLVTT
2304	yfluxlat(:) = flat
2305	!
2306	!! Test sur iflag_split retire le 2/02/2018, sans vraiment comprendre la raison de ce test. (jyg)
2307	!! IF (iflag_split .eq.0) THEN
2308	DO j=1,knon
2309	Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)2+yv(j,1)2)) * &
2310	ypplay(j,1)/(RD*yt(j,1))
2311	ENDDO
2312	!! ENDIF ! (iflag_split .eq.0)
2313
2314	DO j = 1, knon
2315	yt1_new=(1./RCPD)(AcoefH(j)+BcoefH(j)y_flux_t1(j)*dtime)
2316	ytsurf_new(j)=yt1_new-y_flux_t1(j)/(Kech_h(j)*RCPD)
2317	ENDDO
2318
2319	DO j=1,knon
2320	y_d_ts(j) = ytsurf_new(j) - yts(j)
2321	ENDDO
2322
2323	ELSE ! (ok_flux_surf)
2324	DO j=1,knon
2325	y_flux_t1(j) = yfluxsens(j)
2326	y_flux_q1(j) = -yevap(j)
2327	ENDDO
2328	ENDIF ! (ok_flux_surf)
2329
2330	! flux of blowing snow at the first level
2331	IF (ok_bs) THEN
2332	DO j=1,knon
2333	y_flux_bs(j)=yfluxbs(j)
2334	ENDDO
2335	ENDIF
2336	!
2337	! ------------------------------------------------------------------------------
2338	! 12a) Splitting
2339	! ------------------------------------------------------------------------------
2340
2341	IF (iflag_split .GE. 1) THEN
2342	!
2343	IF (nsrf .ne. is_oce) THEN
2344	!
2345	! Compute potential evaporation and aridity factor (jyg, 20200328)
2346	ybeta_prev(:) = ybeta(:)
2347	DO j = 1, knon
2348	yqa(j) = AcoefQ(j) - BcoefQ(j)yevap(j)dtime
2349	ENDDO
2350	!
2351	CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot)
2352	!
2353	ybeta(1:knon) = min(yevap(1:knon)/yevap_pot(1:knon), 1.)
2354
2355	IF (prt_level >=10) THEN
2356	DO j=1,knon
2357	print*,'y_flux_t1,yfluxlat,wakes' &
2358	& , y_flux_t1(j), yfluxlat(j), ywake_s(j)
2359	print*,'beta_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j)
2360	print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
2361	ENDDO
2362	ENDIF ! (prt_level >=10)
2363	!
2364	! Second call to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account
2365	! the update of the aridity coeficient beta.
2366	!
2367	CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, &
2368	BcoefQ_x, BcoefQ_w &
2369	)
2370	CALL wx_pbl0_merge(knon, ypplay, ypaprs, &
2371	ywake_s, ydTs0, ydqs0, &
2372	yt_x, yt_w, yq_x, yq_w, &
2373	yu_x, yu_w, yv_x, yv_w, &
2374	ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
2375	ycdragm_x, ycdragm_w, &
2376	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
2377	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
2378	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
2379	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
2380	AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
2381	BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
2382	ycdragh, ycdragq, ycdragm, &
2383	yt1, yq1, yu1, yv1 &
2384	)
2385	IF (iflag_split .eq. 2) THEN
2386	CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
2387	ywake_s, ybeta, ywake_cstar, ywake_dens, &
2388	AcoefH_x, AcoefH_w, &
2389	BcoefH_x, BcoefH_w, &
2390	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2391	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2392	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2393	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2394	yg_T, yg_Q, &
2395	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2396	ydTs_ins, ydqs_ins &
2397	)
2398	ELSE !
2399	AcoefH(:) = AcoefH_0(:)
2400	AcoefQ(:) = AcoefQ_0(:)
2401	BcoefH(:) = BcoefH_0(:)
2402	BcoefQ(:) = BcoefQ_0(:)
2403	yg_T(:) = 0.
2404	yg_Q(:) = 0.
2405	yGamma_dTs_phiT(:) = 0.
2406	yGamma_dQs_phiQ(:) = 0.
2407	ydTs_ins(:) = 0.
2408	ydqs_ins(:) = 0.
2409	ENDIF ! (iflag_split .eq. 2)
2410	!
2411	ELSE ! (nsrf .ne. is_oce)
2412	ybeta(1:knon) = 1.
2413	yevap_pot(1:knon) = yevap(1:knon)
2414	AcoefH(:) = AcoefH_0(:)
2415	AcoefQ(:) = AcoefQ_0(:)
2416	BcoefH(:) = BcoefH_0(:)
2417	BcoefQ(:) = BcoefQ_0(:)
2418	yg_T(:) = 0.
2419	yg_Q(:) = 0.
2420	yGamma_dTs_phiT(:) = 0.
2421	yGamma_dQs_phiQ(:) = 0.
2422	ydTs_ins(:) = 0.
2423	ydqs_ins(:) = 0.
2424	ENDIF ! (nsrf .ne. is_oce)
2425	!
2426	CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, &
2427	yg_T, yg_Q, &
2428	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2429	ydTs_ins, ydqs_ins, &
2430	y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
2431	!!!! HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
2432	phiQ0_b, phiT0_b, &
2433	y_flux_t1_x, y_flux_t1_w, &
2434	y_flux_q1_x, y_flux_q1_w, &
2435	y_flux_u1_x, y_flux_u1_w, &
2436	y_flux_v1_x, y_flux_v1_w, &
2437	yfluxlat_x, yfluxlat_w, &
2438	y_delta_qsats, &
2439	y_delta_tsurf_new, y_delta_qsurf &
2440	)
2441	!
2442	CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
2443	yTs, y_delta_tsurf, &
2444	yqsurf, yTsurf_new, &
2445	y_delta_tsurf_new, y_delta_qsats, &
2446	AcoefH_x, AcoefH_w, &
2447	BcoefH_x, BcoefH_w, &
2448	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2449	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2450	y_flux_t1, y_flux_q1, &
2451	y_flux_t1_x, y_flux_t1_w, &
2452	y_flux_q1_x, y_flux_q1_w)
2453	!
2454	IF (nsrf .ne. is_oce) THEN
2455	CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
2456	yTs, y_delta_tsurf, &
2457	yqsurf, yTsurf_new, &
2458	y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf, &
2459	AcoefH_x, AcoefH_w, &
2460	BcoefH_x, BcoefH_w, &
2461	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2462	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2463	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2464	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2465	yg_T, yg_Q, &
2466	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2467	ydTs_ins, ydqs_ins, &
2468	y_flux_t1, y_flux_q1, &
2469	y_flux_t1_x, y_flux_t1_w, &
2470	y_flux_q1_x, y_flux_q1_w )
2471	ENDIF ! (nsrf .ne. is_oce)
2472	!
2473	ELSE ! (iflag_split .ge. 1)
2474	ybeta(1:knon) = 1.
2475	yevap_pot(1:knon) = yevap(1:knon)
2476	ENDIF ! (iflag_split .ge. 1)
2477	!
2478	IF (prt_level >= 10) THEN
2479	print *,'pbl_surface, ybeta , yevap, yevap_pot ', &
2480	ybeta(1:knon) , yevap(1:knon), yevap_pot(1:knon)
2481	ENDIF ! (prt_level >= 10)
2482	!
2483	!>jyg
2484	!
2485
2486	!!jyg!! A reprendre apres reflexion ===============================================
2487	!!jyg!!
2488	!!jyg!! DO j=1,knon
2489	!!jyg!!!!! nrlmd le 13/06/2011
2490	!!jyg!!
2491	!!jyg!!!----Diffusion dans le sol dans le cas continental seulement
2492	!!jyg!! IF (nsrf.eq.is_ter) THEN
2493	!!jyg!!!----Calcul du coefficient delta_coeff
2494	!!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)))
2495	!!jyg!!
2496	!!jyg!!! delta_coef(j)=dtime/(inertia*sqrt(tau_eq(j)))
2497	!!jyg!! delta_coef(j)=facteur*sqrt(tau_eq(j))/inertia
2498	!!jyg!!! delta_coef(j)=0.
2499	!!jyg!! ELSE
2500	!!jyg!! delta_coef(j)=0.
2501	!!jyg!! ENDIF
2502	!!jyg!!
2503	!!jyg!!!----Calcul de delta_tsurf
2504	!!jyg!! y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j)
2505	!!jyg!!
2506	!!jyg!!!----Si il n'y a pas des poches...
2507	!!jyg!! IF (wake_cstar(j).le.0.01) THEN
2508	!!jyg!! y_delta_tsurf(j)=0.
2509	!!jyg!! y_delta_flux_t1(j)=0.
2510	!!jyg!! ENDIF
2511	!!jyg!!
2512	!!jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle)
2513	!!jyg!!!!!!! jyg le 23/02/2012
2514	!!jyg!!!!!!!
2515	!!jyg!!!! ybeta(j)=y_flux_q1(j) / &
2516	!!jyg!!!! & (Kech_h(j)*(yq(j,1)-yqsatsurf(j)))
2517	!!jyg!!!!!! ybeta(j)=-1.*yevap(j) / &
2518	!!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)))
2519	!!jyg!!!!!!! fin jyg
2520	!!jyg!!!!!
2521	!!jyg!!
2522	!!jyg!! ENDDO
2523	!!jyg!!
2524	!!jyg!!!!! fin nrlmd le 13/06/2011
2525	!!jyg!!
2526	IF (iflag_split .ge. 1) THEN
2527	IF (prt_level >=10) THEN
2528	DO j = 1, knon
2529	print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
2530	print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
2531	print*,'t1x, t1w, t1, t1_ancien', &
2532	& yt_x(j,1), yt_w(j,1), yt(j,1), t(j,1)
2533	print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j)
2534	ENDDO
2535
2536	DO j=1,knon
2537	print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
2538	& , 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)
2539	print*,'beta, ytsurf_new ', ybeta(j), ytsurf_new(j)
2540	print*,'inertia, facteur, cstar', inertia, facteur,wake_cstar(j)
2541	ENDDO
2542	ENDIF ! (prt_level >=10)
2543
2544	!!! jyg le 07/02/2012
2545	ENDIF ! (iflag_split .ge.1)
2546	!!!
2547
2548	!!! jyg le 07/02/2012
2549	IF (iflag_split .eq.0) THEN
2550	!!!
2551	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2552	CALL climb_hq_up(knon, dtime, yt, yq, &
2553	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
2554	!!! jyg le 07/02/2012
2555	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2556	CcoefH, CcoefQ, DcoefH, DcoefQ, &
2557	Kcoef_hq, gama_q, gama_h, &
2558	!!!
2559	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
2560	ELSE !(iflag_split .eq.0)
2561	CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
2562	y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &
2563	!!! nrlmd le 02/05/2011
2564	AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &
2565	CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
2566	Kcoef_hq_x, gama_q_x, gama_h_x, &
2567	!!!
2568	y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:))
2569	!
2570	CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
2571	y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &
2572	!!! nrlmd le 02/05/2011
2573	AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &
2574	CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
2575	Kcoef_hq_w, gama_q_w, gama_h_w, &
2576	!!!
2577	y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:))
2578	!!!
2579	ENDIF ! (iflag_split .eq.0)
2580	!!!
2581
2582	!!! jyg le 07/02/2012
2583	IF (iflag_split .eq.0) THEN
2584	!!!
2585	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2586	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
2587	!!! jyg le 07/02/2012
2588	AcoefU, AcoefV, BcoefU, BcoefV, &
2589	CcoefU, CcoefV, DcoefU, DcoefV, &
2590	Kcoef_m, &
2591	!!!
2592	y_flux_u, y_flux_v, y_d_u, y_d_v)
2593	y_d_t_diss(:,:)=0.
2594	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2595	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2596	& ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
2597	& ,iflag_pbl)
2598	ENDIF
2599	! print*,'yamada_c OK'
2600
2601	ELSE !(iflag_split .eq.0)
2602	CALL climb_wind_up(knon, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &
2603	!!! nrlmd le 02/05/2011
2604	AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &
2605	CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
2606	Kcoef_m_x, &
2607	!!!
2608	y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)
2609	!
2610	y_d_t_diss_x(:,:)=0.
2611	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2612	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2613	& ,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 &
2614	,ycoefq_x,y_d_t_diss_x,yustar_x &
2615	& ,iflag_pbl)
2616	ENDIF
2617	! print*,'yamada_c OK'
2618
2619	CALL climb_wind_up(knon, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &
2620	!!! nrlmd le 02/05/2011
2621	AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &
2622	CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
2623	Kcoef_m_w, &
2624	!!!
2625	y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)
2626	!!!
2627	y_d_t_diss_w(:,:)=0.
2628	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2629	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2630	& ,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 &
2631	,ycoefq_w,y_d_t_diss_w,yustar_w &
2632	& ,iflag_pbl)
2633	ENDIF
2634	! print*,'yamada_c OK'
2635	!
2636	IF (prt_level >=10) THEN
2637	print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &
2638	yfluxlat_x(1:knon), yfluxlat_w(1:knon)
2639	ENDIF
2640	!
2641	ENDIF ! (iflag_split .eq.0)
2642
2643	IF (ok_bs) THEN
2644	CALL climb_qbs_up(knon, dtime, yqbs, &
2645	y_flux_bs, ypaprs, ypplay, &
2646	AcoefQBS, BcoefQBS, &
2647	CcoefQBS, DcoefQBS, &
2648	Kcoef_qbs, gama_qbs, &
2649	y_flux_qbs(:,:), y_d_qbs(:,:))
2650	ENDIF
2651
2652	!!!
2653	!!
2654	!! DO j = 1, knon
2655	!! y_dflux_t(j) = y_dflux_t(j) * ypct(j)
2656	!! y_dflux_q(j) = y_dflux_q(j) * ypct(j)
2657	!! ENDDO
2658	!!
2659	!****************************************************************************************
2660	! 13) Transform variables for output format :
2661	! - Decompress
2662	! - Multiply with pourcentage of current surface
2663	! - Cumulate in global variable
2664	!
2665	!****************************************************************************************
2666
2667
2668	!!! jyg le 07/02/2012
2669	IF (iflag_split.EQ.0) THEN
2670	!!!
2671	DO k = 1, klev
2672	DO j = 1, knon
2673	i = ni(j)
2674	y_d_t_diss(j,k) = y_d_t_diss(j,k) * ypct(j)
2675	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
2676	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
2677	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
2678	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
2679	!FC
2680	IF (nsrf .EQ. is_ter .and. ifl_pbltree .GE. 1) THEN
2681	! if (y_d_u_frein(j,k).ne.0. ) then
2682	! print, nsrf,'IS_TER ++', y_d_u_frein(j,k)ypct(j),y_d_u(j,k),j,k
2683	! ENDIF
2684	y_d_u(j,k) =y_d_u(j,k) + y_d_u_frein(j,k)*ypct(j)
2685	y_d_v(j,k) =y_d_v(j,k) + y_d_v_frein(j,k)*ypct(j)
2686	treedrg(i,k,nsrf)=y_treedrg(j,k)
2687	ELSE
2688	treedrg(i,k,nsrf)=0.
2689	ENDIF
2690	!FC
2691	flux_t(i,k,nsrf) = y_flux_t(j,k)
2692	flux_q(i,k,nsrf) = y_flux_q(j,k)
2693	flux_u(i,k,nsrf) = y_flux_u(j,k)
2694	flux_v(i,k,nsrf) = y_flux_v(j,k)
2695	ENDDO
2696	ENDDO
2697
2698	ELSE !(iflag_split .eq.0)
2699
2700	! Tendances hors poches
2701	DO k = 1, klev
2702	DO j = 1, knon
2703	i = ni(j)
2704	y_d_t_diss_x(j,k) = y_d_t_diss_x(j,k) * ypct(j)
2705	y_d_t_x(j,k) = y_d_t_x(j,k) * ypct(j)
2706	y_d_q_x(j,k) = y_d_q_x(j,k) * ypct(j)
2707	y_d_u_x(j,k) = y_d_u_x(j,k) * ypct(j)
2708	y_d_v_x(j,k) = y_d_v_x(j,k) * ypct(j)
2709
2710	flux_t_x(i,k,nsrf) = y_flux_t_x(j,k)
2711	flux_q_x(i,k,nsrf) = y_flux_q_x(j,k)
2712	flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
2713	flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)
2714	ENDDO
2715	ENDDO
2716
2717	! Tendances dans les poches
2718	DO k = 1, klev
2719	DO j = 1, knon
2720	i = ni(j)
2721	y_d_t_diss_w(j,k) = y_d_t_diss_w(j,k) * ypct(j)
2722	y_d_t_w(j,k) = y_d_t_w(j,k) * ypct(j)
2723	y_d_q_w(j,k) = y_d_q_w(j,k) * ypct(j)
2724	y_d_u_w(j,k) = y_d_u_w(j,k) * ypct(j)
2725	y_d_v_w(j,k) = y_d_v_w(j,k) * ypct(j)
2726
2727	flux_t_w(i,k,nsrf) = y_flux_t_w(j,k)
2728	flux_q_w(i,k,nsrf) = y_flux_q_w(j,k)
2729	flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
2730	flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)
2731	ENDDO
2732	ENDDO
2733
2734	! Flux, tendances et Tke moyenne dans la maille
2735	DO k = 1, klev
2736	DO j = 1, knon
2737	i = ni(j)
2738	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))
2739	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))
2740	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))
2741	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))
2742	ENDDO
2743	ENDDO
2744	DO j=1,knon
2745	yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j))
2746	ENDDO
2747	IF (prt_level >=10) THEN
2748	print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &
2749	nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf)
2750	ENDIF
2751
2752	DO k = 1, klev
2753	DO j = 1, knon
2754	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))
2755	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))
2756	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))
2757	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))
2758	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))
2759	ENDDO
2760	ENDDO
2761
2762	ENDIF ! (iflag_split .eq.0)
2763	!!!
2764
2765	! tendencies of blowing snow
2766	IF (ok_bs) THEN
2767	DO k = 1, klev
2768	DO j = 1, knon
2769	i = ni(j)
2770	y_d_qbs(j,k)=y_d_qbs(j,k) * ypct(j)
2771	flux_qbs(i,k,nsrf) = y_flux_qbs(j,k)
2772	ENDDO
2773	ENDDO
2774	ENDIF
2775
2776
2777	DO j = 1, knon
2778	i = ni(j)
2779	evap(i,nsrf) = - flux_q(i,1,nsrf) !jyg
2780	if (ok_bs) then ; snowerosion(i,nsrf)=flux_qbs(i,1,nsrf); endif
2781	beta(i,nsrf) = ybeta(j) !jyg
2782	d_ts(i,nsrf) = y_d_ts(j)
2783	!albedo SB >>>
2784	DO k=1,nsw
2785	alb_dir(i,k,nsrf) = yalb_dir_new(j,k)
2786	alb_dif(i,k,nsrf) = yalb_dif_new(j,k)
2787	ENDDO
2788	!albedo SB <<<
2789	snow(i,nsrf) = ysnow(j)
2790	qsurf(i,nsrf) = yqsurf(j)
2791	z0m(i,nsrf) = yz0m(j)
2792	z0h(i,nsrf) = yz0h(j)
2793	fluxlat(i,nsrf) = yfluxlat(j)
2794	agesno(i,nsrf) = yagesno(j)
2795	cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
2796	cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
2797	dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)
2798	dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)
2799	ENDDO
2800
2801	! print*,'Dans pbl OK2'
2802
2803	!!! jyg le 07/02/2012
2804	IF (iflag_split .ge.1) THEN
2805	!!!
2806	!!! nrlmd le 02/05/2011
2807	DO j = 1, knon
2808	i = ni(j)
2809	fluxlat_x(i,nsrf) = yfluxlat_x(j)
2810	fluxlat_w(i,nsrf) = yfluxlat_w(j)
2811	!!!
2812	!!! nrlmd le 13/06/2011
2813	!!jyg20170131 delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
2814	!!jyg20210118 delta_tsurf(i,nsrf)=y_delta_tsurf(j)
2815	delta_tsurf(i,nsrf)=y_delta_tsurf_new(j)
2816	!
2817	delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j)*ypct(j)
2818	!
2819	cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
2820	cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)
2821	cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j)
2822	cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j)
2823	kh(i) = kh(i) + Kech_h(j)*ypct(j)
2824	kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j)
2825	kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j)
2826	!!!
2827	ENDDO
2828	!!!
2829	ENDIF ! (iflag_split .ge.1)
2830	!!!
2831	!!! nrlmd le 02/05/2011
2832	!!jyg le 20/02/2011
2833	!! tke_x(:,:,nsrf)=0.
2834	!! tke_w(:,:,nsrf)=0.
2835	!!jyg le 20/02/2011
2836	!! DO k = 1, klev+1
2837	!! DO j = 1, knon
2838	!! i = ni(j)
2839	!! wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2840	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2841	!! ENDDO
2842	!! ENDDO
2843	!!jyg le 20/02/2011
2844	!! DO k = 1, klev+1
2845	!! DO j = 1, knon
2846	!! i = ni(j)
2847	!! tke(i,k,nsrf)=(1.-ywake_s(j))tke_x(i,k,nsrf)+ywake_s(j)tke_w(i,k,nsrf)
2848	!! ENDDO
2849	!! ENDDO
2850	!!!
2851	IF (iflag_split .eq.0) THEN
2852	wake_dltke(:,:,nsrf) = 0.
2853	DO k = 1, klev+1
2854	DO j = 1, knon
2855	i = ni(j)
2856	!jyg<
2857	!! tke(i,k,nsrf) = ytke(j,k)
2858	!! tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
2859	tke_x(i,k,nsrf) = ytke(j,k)
2860	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
2861	eps_x(i,k,nsrf) = yeps(j,k)
2862	eps_x(i,k,is_ave) = eps_x(i,k,is_ave) + yeps(j,k)*ypct(j)
2863	!>jyg
2864	ENDDO
2865	ENDDO
2866
2867	ELSE ! (iflag_split .eq.0)
2868	DO k = 1, klev+1
2869	DO j = 1, knon
2870	i = ni(j)
2871	wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2872	!jyg<
2873	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2874	!! tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
2875	tke_x(i,k,nsrf) = ytke_x(j,k)
2876	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)
2877	eps_x(i,k,nsrf) = yeps_x(j,k)
2878	eps_x(i,k,is_ave) = eps_x(i,k,is_ave) + eps_x(i,k,nsrf)*ypct(j)
2879	wake_dltke(i,k,is_ave) = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
2880
2881
2882	!>jyg
2883	ENDDO
2884	ENDDO
2885	ENDIF ! (iflag_split .eq.0)
2886	!!!
2887	DO k = 2, klev
2888	DO j = 1, knon
2889	i = ni(j)
2890	zcoefh(i,k,nsrf) = ycoefh(j,k)
2891	zcoefm(i,k,nsrf) = ycoefm(j,k)
2892	zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
2893	zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
2894	ENDDO
2895	ENDDO
2896
2897	! print*,'Dans pbl OK3'
2898
2899	IF ( nsrf .EQ. is_ter ) THEN
2900	DO j = 1, knon
2901	i = ni(j)
2902	qsol(i) = yqsol(j)
2903	ENDDO
2904	ENDIF
2905
2906	!jyg<
2907	!! ftsoil(:,:,nsrf) = 0.
2908	!>jyg
2909	DO k = 1, nsoilmx
2910	DO j = 1, knon
2911	i = ni(j)
2912	ftsoil(i, k, nsrf) = ytsoil(j,k)
2913	ENDDO
2914	ENDDO
2915
2916	!!! jyg le 07/02/2012
2917	IF (iflag_split .ge.1) THEN
2918	!!!
2919	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
2920	DO k = 1, klev
2921	DO j = 1, knon
2922	i = ni(j)
2923	d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)
2924	d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)
2925	d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)
2926	d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)
2927	d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)
2928	!
2929	d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)
2930	d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)
2931	d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)
2932	d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
2933	d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
2934	!
2935	!! d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
2936	!! d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)
2937	ENDDO
2938	ENDDO
2939	!!!
2940	ENDIF ! (iflag_split .ge.1)
2941	!!!
2942
2943	DO k = 1, klev
2944	DO j = 1, knon
2945	i = ni(j)
2946	d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
2947	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
2948	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
2949	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
2950	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
2951	ENDDO
2952	ENDDO
2953
2954
2955	IF (ok_bs) THEN
2956	DO k = 1, klev
2957	DO j = 1, knon
2958	i = ni(j)
2959	d_qbs(i,k) = d_qbs(i,k) + y_d_qbs(j,k)
2960	ENDDO
2961	ENDDO
2962	ENDIF
2963
2964	! print*,'Dans pbl OK4'
2965
2966	IF (prt_level >=10) THEN
2967	print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
2968	d_t_w(1:knon,1), d_t_x(1:knon,1), d_t(1:knon,1)
2969	ENDIF
2970
2971	if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
2972	delta_sal = missing_val
2973	ds_ns = missing_val
2974	dt_ns = missing_val
2975	delta_sst = missing_val
2976	dter = missing_val
2977	dser = missing_val
2978	tkt = missing_val
2979	tks = missing_val
2980	taur = missing_val
2981	sss = missing_val
2982
2983	delta_sal(ni(:knon)) = ydelta_sal(:knon)
2984	ds_ns(ni(:knon)) = yds_ns(:knon)
2985	dt_ns(ni(:knon)) = ydt_ns(:knon)
2986	delta_sst(ni(:knon)) = ydelta_sst(:knon)
2987	dter(ni(:knon)) = ydter(:knon)
2988	dser(ni(:knon)) = ydser(:knon)
2989	tkt(ni(:knon)) = ytkt(:knon)
2990	tks(ni(:knon)) = ytks(:knon)
2991	taur(ni(:knon)) = ytaur(:knon)
2992	sss(ni(:knon)) = ysss(:knon)
2993
2994	if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
2995	dt_ds = missing_val
2996	dt_ds(ni(:knon)) = ydt_ds(:knon)
2997	end if
2998	end if
2999
3000
3001	!****************************************************************************************
3002	! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m
3003	! Call HBTM
3004	!
3005	!****************************************************************************************
3006	!!!
3007	!
3008	#undef T2m
3009	#define T2m
3010	#ifdef T2m
3011	! Calculations of diagnostic t,q at 2m and u, v at 10m
3012
3013	! print*,'Dans pbl OK41'
3014	! print*,'tair1,yt(:,1),y_d_t(:,1)'
3015	! print*, tair1,yt(:,1),y_d_t(:,1)
3016	!!! jyg le 07/02/2012
3017	IF (iflag_split .eq.0) THEN
3018	DO j=1, knon
3019	uzon(j) = yu(j,1) + y_d_u(j,1)
3020	vmer(j) = yv(j,1) + y_d_v(j,1)
3021	tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
3022	qair1(j) = yq(j,1) + y_d_q(j,1)
3023	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
3024	* (ypaprs(j,1)-ypplay(j,1))
3025	tairsol(j) = yts(j) + y_d_ts(j)
3026	qairsol(j) = yqsurf(j)
3027	ENDDO
3028	ELSE ! (iflag_split .eq.0)
3029	DO j=1, knon
3030	uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
3031	vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
3032	tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
3033	qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
3034	zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
3035	* (ypaprs(j,1)-ypplay(j,1))
3036	tairsol(j) = yts(j) + y_d_ts(j)
3037	!! tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
3038	tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf_new(j)
3039	qairsol(j) = yqsurf(j)
3040	ENDDO
3041	DO j=1, knon
3042	uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
3043	vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
3044	tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
3045	qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
3046	zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
3047	* (ypaprs(j,1)-ypplay(j,1))
3048	tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
3049	qairsol(j) = yqsurf(j)
3050	ENDDO
3051	!!!
3052	ENDIF ! (iflag_split .eq.0)
3053	!!!
3054	DO j=1, knon
3055	! i = ni(j)
3056	! yz0h_oupas(j) = yz0m(j)
3057	! IF(nsrf.EQ.is_oce) THEN
3058	! yz0h_oupas(j) = z0m(i,nsrf)
3059	! ENDIF
3060	psfce(j)=ypaprs(j,1)
3061	patm(j)=ypplay(j,1)
3062	ENDDO
3063
3064	IF (iflag_pbl_surface_t2m_bug==1) THEN
3065	yz0h_oupas(1:knon)=yz0m(1:knon)
3066	ELSE
3067	yz0h_oupas(1:knon)=yz0h(1:knon)
3068	ENDIF
3069
3070	! print*,'Dans pbl OK42A'
3071	! print*,'tair1,yt(:,1),y_d_t(:,1)'
3072	! print*, tair1,yt(:,1),y_d_t(:,1)
3073
3074	! Calculate the temperature and relative humidity at 2m and the wind at 10m
3075	!!! jyg le 07/02/2012
3076	IF (iflag_split .eq.0) THEN
3077	IF (iflag_new_t2mq2m==1) THEN
3078	CALL stdlevvarn(klon, knon, nsrf, zxli, &
3079	uzon, vmer, tair1, qair1, zgeo1, &
3080	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3081	yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
3082	yn2mout(:, nsrf, :))
3083	ELSE
3084	CALL stdlevvar(klon, knon, nsrf, zxli, &
3085	uzon, vmer, tair1, qair1, zgeo1, &
3086	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3087	yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
3088	ENDIF
3089	ELSE !(iflag_split .eq.0)
3090	IF (iflag_new_t2mq2m==1) THEN
3091	CALL stdlevvarn(klon, knon, nsrf, zxli, &
3092	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
3093	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3094	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, &
3095	yn2mout_x(:, nsrf, :))
3096	CALL stdlevvarn(klon, knon, nsrf, zxli, &
3097	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
3098	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3099	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, &
3100	yn2mout_w(:, nsrf, :))
3101	ELSE
3102	CALL stdlevvar(klon, knon, nsrf, zxli, &
3103	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
3104	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3105	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, ypblh_x, rain_f, zxtsol)
3106	CALL stdlevvar(klon, knon, nsrf, zxli, &
3107	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
3108	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3109	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, ypblh_w, rain_f, zxtsol)
3110	ENDIF
3111	!!!
3112	ENDIF ! (iflag_split .eq.0)
3113	!!!
3114	!!! jyg le 07/02/2012
3115	IF (iflag_split .eq.0) THEN
3116	DO j=1, knon
3117	i = ni(j)
3118	t2m(i,nsrf)=yt2m(j)
3119	q2m(i,nsrf)=yq2m(j)
3120	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3121	ustar(i,nsrf)=yustar(j)
3122	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
3123	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
3124	!
3125	DO k = 1, 6
3126	n2mout(i,nsrf,k) = yn2mout(j,nsrf,k)
3127	END DO
3128	!
3129	ENDDO
3130	ELSE !(iflag_split .eq.0)
3131	DO j=1, knon
3132	i = ni(j)
3133	t2m_x(i,nsrf)=yt2m_x(j)
3134	q2m_x(i,nsrf)=yq2m_x(j)
3135	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3136	ustar_x(i,nsrf)=yustar_x(j)
3137	u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
3138	v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
3139	!
3140	DO k = 1, 6
3141	n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k)
3142	END DO
3143	!
3144	ENDDO
3145	DO j=1, knon
3146	i = ni(j)
3147	t2m_w(i,nsrf)=yt2m_w(j)
3148	q2m_w(i,nsrf)=yq2m_w(j)
3149	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3150	ustar_w(i,nsrf)=yustar_w(j)
3151	u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
3152	v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
3153	!
3154	ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
3155	u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
3156	v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
3157	!
3158	DO k = 1, 6
3159	n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k)
3160	END DO
3161	!
3162	ENDDO
3163	!!!
3164	ENDIF ! (iflag_split .eq.0)
3165	!!!
3166
3167	! print*,'Dans pbl OK43'
3168	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
3169	!IM Ajoute dependance type surface
3170	IF (thermcep) THEN
3171	!!! jyg le 07/02/2012
3172	IF (iflag_split .eq.0) THEN
3173	DO j = 1, knon
3174	i=ni(j)
3175	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
3176	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
3177	zx_qs1 = MIN(0.5,zx_qs1)
3178	zcor1 = 1./(1.-RETV*zx_qs1)
3179	zx_qs1 = zx_qs1*zcor1
3180
3181	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
3182	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
3183	ENDDO
3184	ELSE ! (iflag_split .eq.0)
3185	DO j = 1, knon
3186	i=ni(j)
3187	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
3188	zx_qs1 = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)
3189	zx_qs1 = MIN(0.5,zx_qs1)
3190	zcor1 = 1./(1.-RETV*zx_qs1)
3191	zx_qs1 = zx_qs1*zcor1
3192
3193	rh2m_x(i) = rh2m_x(i) + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
3194	qsat2m_x(i) = qsat2m_x(i) + zx_qs1 * pctsrf(i,nsrf)
3195	ENDDO
3196	DO j = 1, knon
3197	i=ni(j)
3198	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
3199	zx_qs1 = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
3200	zx_qs1 = MIN(0.5,zx_qs1)
3201	zcor1 = 1./(1.-RETV*zx_qs1)
3202	zx_qs1 = zx_qs1*zcor1
3203
3204	rh2m_w(i) = rh2m_w(i) + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
3205	qsat2m_w(i) = qsat2m_w(i) + zx_qs1 * pctsrf(i,nsrf)
3206	ENDDO
3207	!!!
3208	ENDIF ! (iflag_split .eq.0)
3209	!!!
3210	ENDIF
3211	!
3212	IF (prt_level >=10) THEN
3213	print *, 'T2m, q2m, RH2m ', &
3214	t2m(1:knon,:), q2m(1:knon,:), rh2m(1:knon)
3215	ENDIF
3216
3217	! print*,'OK pbl 5'
3218	!
3219	!!! jyg le 07/02/2012
3220	IF (iflag_split .eq.0) THEN
3221	CALL hbtm(knon, ypaprs, ypplay, &
3222	yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
3223	y_flux_t,y_flux_q,yu,yv,yt,yq, &
3224	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
3225	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
3226	IF (prt_level >=10) THEN
3227	print *,' Arg. de HBTM: yt2m ',yt2m(1:knon)
3228	print *,' Arg. de HBTM: yt10m ',yt10m(1:knon)
3229	print *,' Arg. de HBTM: yq2m ',yq2m(1:knon)
3230	print *,' Arg. de HBTM: yq10m ',yq10m(1:knon)
3231	print *,' Arg. de HBTM: yustar ',yustar(1:knon)
3232	print *,' Arg. de HBTM: y_flux_t ',y_flux_t(1:knon,:)
3233	print *,' Arg. de HBTM: y_flux_q ',y_flux_q(1:knon,:)
3234	print *,' Arg. de HBTM: yu ',yu(1:knon,:)
3235	print *,' Arg. de HBTM: yv ',yv(1:knon,:)
3236	print *,' Arg. de HBTM: yt ',yt(1:knon,:)
3237	print *,' Arg. de HBTM: yq ',yq(1:knon,:)
3238	ENDIF
3239	ELSE ! (iflag_split .eq.0)
3240	CALL HBTM(knon, ypaprs, ypplay, &
3241	yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
3242	y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
3243	ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
3244	ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
3245	IF (prt_level >=10) THEN
3246	print *,' Arg. de HBTM: yt2m_x ',yt2m_x(1:knon)
3247	print *,' Arg. de HBTM: yt10m_x ',yt10m_x(1:knon)
3248	print *,' Arg. de HBTM: yq2m_x ',yq2m_x(1:knon)
3249	print *,' Arg. de HBTM: yq10m_x ',yq10m_x(1:knon)
3250	print *,' Arg. de HBTM: yustar_x ',yustar_x(1:knon)
3251	print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x(1:knon,:)
3252	print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x(1:knon,:)
3253	print *,' Arg. de HBTM: yu_x ',yu_x(1:knon,:)
3254	print *,' Arg. de HBTM: yv_x ',yv_x(1:knon,:)
3255	print *,' Arg. de HBTM: yt_x ',yt_x(1:knon,:)
3256	print *,' Arg. de HBTM: yq_x ',yq_x(1:knon,:)
3257	ENDIF
3258	CALL HBTM(knon, ypaprs, ypplay, &
3259	yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
3260	y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
3261	ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
3262	ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
3263	!!!
3264	ENDIF ! (iflag_split .eq.0)
3265	!!!
3266
3267	!!! jyg le 07/02/2012
3268	IF (iflag_split .eq.0) THEN
3269	!!!
3270	DO j=1, knon
3271	i = ni(j)
3272	pblh(i,nsrf) = ypblh(j)
3273	wstar(i,nsrf) = ywstar(j)
3274	plcl(i,nsrf) = ylcl(j)
3275	capCL(i,nsrf) = ycapCL(j)
3276	oliqCL(i,nsrf) = yoliqCL(j)
3277	cteiCL(i,nsrf) = ycteiCL(j)
3278	pblT(i,nsrf) = ypblT(j)
3279	therm(i,nsrf) = ytherm(j)
3280	trmb1(i,nsrf) = ytrmb1(j)
3281	trmb2(i,nsrf) = ytrmb2(j)
3282	trmb3(i,nsrf) = ytrmb3(j)
3283	ENDDO
3284	IF (prt_level >=10) THEN
3285	print *, 'After HBTM: pblh ', pblh(1:knon,:)
3286	print *, 'After HBTM: plcl ', plcl(1:knon,:)
3287	print *, 'After HBTM: cteiCL ', cteiCL(1:knon,:)
3288	ENDIF
3289	ELSE !(iflag_split .eq.0)
3290	DO j=1, knon
3291	i = ni(j)
3292	pblh_x(i,nsrf) = ypblh_x(j)
3293	wstar_x(i,nsrf) = ywstar_x(j)
3294	plcl_x(i,nsrf) = ylcl_x(j)
3295	capCL_x(i,nsrf) = ycapCL_x(j)
3296	oliqCL_x(i,nsrf) = yoliqCL_x(j)
3297	cteiCL_x(i,nsrf) = ycteiCL_x(j)
3298	pblT_x(i,nsrf) = ypblT_x(j)
3299	therm_x(i,nsrf) = ytherm_x(j)
3300	trmb1_x(i,nsrf) = ytrmb1_x(j)
3301	trmb2_x(i,nsrf) = ytrmb2_x(j)
3302	trmb3_x(i,nsrf) = ytrmb3_x(j)
3303	ENDDO
3304	IF (prt_level >=10) THEN
3305	print *, 'After HBTM: pblh_x ', pblh_x(1:knon,:)
3306	print *, 'After HBTM: plcl_x ', plcl_x(1:knon,:)
3307	print *, 'After HBTM: cteiCL_x ', cteiCL_x(1:knon,:)
3308	ENDIF
3309	DO j=1, knon
3310	i = ni(j)
3311	pblh_w(i,nsrf) = ypblh_w(j)
3312	wstar_w(i,nsrf) = ywstar_w(j)
3313	plcl_w(i,nsrf) = ylcl_w(j)
3314	capCL_w(i,nsrf) = ycapCL_w(j)
3315	oliqCL_w(i,nsrf) = yoliqCL_w(j)
3316	cteiCL_w(i,nsrf) = ycteiCL_w(j)
3317	pblT_w(i,nsrf) = ypblT_w(j)
3318	therm_w(i,nsrf) = ytherm_w(j)
3319	trmb1_w(i,nsrf) = ytrmb1_w(j)
3320	trmb2_w(i,nsrf) = ytrmb2_w(j)
3321	trmb3_w(i,nsrf) = ytrmb3_w(j)
3322	ENDDO
3323	IF (prt_level >=10) THEN
3324	print *, 'After HBTM: pblh_w ', pblh_w(1:knon,:)
3325	print *, 'After HBTM: plcl_w ', plcl_w(1:knon,:)
3326	print *, 'After HBTM: cteiCL_w ', cteiCL_w(1:knon,:)
3327	ENDIF
3328	!!!
3329	ENDIF ! (iflag_split .eq.0)
3330	!!!
3331
3332	! print*,'OK pbl 6'
3333	#else
3334	! T2m not defined
3335	! No calculation
3336	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
3337	#endif
3338
3339	!****************************************************************************************
3340	! 15) End of loop over different surfaces
3341	!
3342	!****************************************************************************************
3343	ENDDO loop_nbsrf
3344	!
3345	!----------------------------------------------------------------------------------------
3346	! Reset iflag_split
3347	!
3348	iflag_split=iflag_split_ref
3349
3350	!****************************************************************************************
3351	! 16) Calculate the mean value over all sub-surfaces for some variables
3352	!
3353	!****************************************************************************************
3354
3355	z0m(:,nbsrf+1) = 0.0
3356	z0h(:,nbsrf+1) = 0.0
3357	DO nsrf = 1, nbsrf
3358	DO i = 1, klon
3359	z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)
3360	z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)
3361	ENDDO
3362	ENDDO
3363
3364	! print*,'OK pbl 7'
3365	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
3366	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
3367	zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
3368	zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
3369	zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
3370	zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
3371
3372	!!! jyg le 07/02/2012
3373	IF (iflag_split .ge.1) THEN
3374	!!!
3375	!!! nrlmd & jyg les 02/05/2011, 05/02/2012
3376
3377	DO nsrf = 1, nbsrf
3378	DO k = 1, klev
3379	DO i = 1, klon
3380	zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
3381	zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
3382	zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
3383	zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
3384	!
3385	zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
3386	zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
3387	zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
3388	zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
3389	ENDDO
3390	ENDDO
3391	ENDDO
3392
3393	DO i = 1, klon
3394	zxsens_x(i) = - zxfluxt_x(i,1)
3395	zxsens_w(i) = - zxfluxt_w(i,1)
3396	ENDDO
3397	!!!
3398	ENDIF ! (iflag_split .ge.1)
3399	!!!
3400
3401	DO nsrf = 1, nbsrf
3402	DO k = 1, klev
3403	DO i = 1, klon
3404	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
3405	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
3406	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
3407	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
3408	ENDDO
3409	ENDDO
3410	ENDDO
3411
3412	DO i = 1, klon
3413	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
3414	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
3415	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
3416	ENDDO
3417
3418	! if blowing snow
3419	if (ok_bs) then
3420	DO nsrf = 1, nbsrf
3421	DO k = 1, klev
3422	DO i = 1, klon
3423	zxfluxqbs(i,k) = zxfluxqbs(i,k) + flux_qbs(i,k,nsrf) * pctsrf(i,nsrf)
3424	ENDDO
3425	ENDDO
3426	ENDDO
3427
3428	DO i = 1, klon
3429	zxsnowerosion(i) = zxfluxqbs(i,1) ! blowings snow flux at the surface
3430	END DO
3431	endif
3432
3433	!!!
3434
3435	!
3436	! Incrementer la temperature du sol
3437	!
3438	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
3439	zt2m(:) = 0.0 ; zq2m(:) = 0.0 ; zn2mout(:,:) = 0
3440	zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
3441	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
3442	!!! jyg le 07/02/2012
3443	s_pblh_x(:) = 0.0 ; s_plcl_x(:) = 0.0
3444	s_pblh_w(:) = 0.0 ; s_plcl_w(:) = 0.0
3445	!!!
3446	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
3447	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
3448	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
3449	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
3450	wstar(:,is_ave)=0.
3451
3452	! print*,'OK pbl 9'
3453
3454	!!! nrlmd le 02/05/2011
3455	zxfluxlat_x(:) = 0.0 ; zxfluxlat_w(:) = 0.0
3456	!!!
3457
3458	DO nsrf = 1, nbsrf
3459	DO i = 1, klon
3460	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
3461
3462	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
3463	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
3464
3465	wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)
3466
3467	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
3468	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
3469	ENDDO
3470	ENDDO
3471	!
3472	!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
3473	IF (iflag_order2_sollw == 1) THEN
3474	meansqT(:) = 0. ! as working buffer
3475	DO nsrf = 1, nbsrf
3476	DO i = 1, klon
3477	meansqT(i) = meansqT(i)+(ts(i,nsrf)-zxtsol(i))*2 pctsrf(i,nsrf)
3478	ENDDO
3479	ENDDO
3480	zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)
3481	ENDIF ! iflag_order2_sollw == 1
3482	!>al1
3483
3484	!!! jyg le 07/02/2012
3485	IF (iflag_split .eq.0) THEN
3486	DO nsrf = 1, nbsrf
3487	DO i = 1, klon
3488	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
3489	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
3490	!
3491	DO k = 1, 6
3492	zn2mout(i,k) = zn2mout(i,k) + n2mout(i,nsrf,k) * pctsrf(i,nsrf)
3493	ENDDO
3494	!
3495	zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
3496	wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
3497	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
3498	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
3499
3500	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
3501	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
3502	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
3503	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
3504	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
3505	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
3506	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
3507	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
3508	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
3509	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
3510	ENDDO
3511	ENDDO
3512	ELSE !(iflag_split .eq.0)
3513	DO nsrf = 1, nbsrf
3514	DO i = 1, klon
3515	!!! nrlmd le 02/05/2011
3516	zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
3517	zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
3518	!!!
3519	!!! jyg le 08/02/2012
3520	!! Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;
3521	!! pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
3522	!! pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
3523	!! pour les autres variables, on sort les valeurs de la region (x).
3524	zt2m(i) = zt2m(i) + (t2m_x(i,nsrf)+wake_s(i)(t2m_w(i,nsrf)-t2m_x(i,nsrf))) pctsrf(i,nsrf)
3525	zq2m(i) = zq2m(i) + q2m_x(i,nsrf) * pctsrf(i,nsrf)
3526	!
3527	DO k = 1, 6
3528	zn2mout(i,k) = zn2mout(i,k) + n2mout_x(i,nsrf,k) * pctsrf(i,nsrf)
3529	ENDDO
3530	!
3531	zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
3532	wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
3533	zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
3534	zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
3535	!
3536	s_pblh(i) = s_pblh(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
3537	s_pblh_x(i) = s_pblh_x(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
3538	s_pblh_w(i) = s_pblh_w(i) + pblh_w(i,nsrf) * pctsrf(i,nsrf)
3539	!
3540	s_plcl(i) = s_plcl(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
3541	s_plcl_x(i) = s_plcl_x(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
3542	s_plcl_w(i) = s_plcl_w(i) + plcl_w(i,nsrf) * pctsrf(i,nsrf)
3543	!
3544	s_capCL(i) = s_capCL(i) + capCL_x(i,nsrf) * pctsrf(i,nsrf)
3545	s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
3546	s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
3547	s_pblT(i) = s_pblT(i) + pblT_x(i,nsrf) * pctsrf(i,nsrf)
3548	s_therm(i) = s_therm(i) + therm_x(i,nsrf) * pctsrf(i,nsrf)
3549	s_trmb1(i) = s_trmb1(i) + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
3550	s_trmb2(i) = s_trmb2(i) + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
3551	s_trmb3(i) = s_trmb3(i) + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
3552	ENDDO
3553	ENDDO
3554	DO i = 1, klon
3555	qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
3556	ENDDO
3557	!!!
3558	ENDIF ! (iflag_split .eq.0)
3559	!!!
3560
3561	IF (check) THEN
3562	amn=MIN(ts(1,is_ter),1000.)
3563	amx=MAX(ts(1,is_ter),-1000.)
3564	DO i=2, klon
3565	amn=MIN(ts(i,is_ter),amn)
3566	amx=MAX(ts(i,is_ter),amx)
3567	ENDDO
3568	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
3569	ENDIF
3570
3571	!jg ?
3572	!!$!
3573	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
3574	!!$! sub-surfaces is distributed.
3575	!!$!
3576	!!$ DO nsrf = 1, nbsrf
3577	!!$ DO i = 1, klon
3578	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
3579	!!$ ts(i,nsrf) = zxtsol(i)
3580	!!$ t2m(i,nsrf) = zt2m(i)
3581	!!$ q2m(i,nsrf) = zq2m(i)
3582	!!$ u10m(i,nsrf) = zu10m(i)
3583	!!$ v10m(i,nsrf) = zv10m(i)
3584	!!$
3585	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
3586	!!$ pblh(i,nsrf) = s_pblh(i)
3587	!!$ plcl(i,nsrf) = s_plcl(i)
3588	!!$ capCL(i,nsrf) = s_capCL(i)
3589	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
3590	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
3591	!!$ pblT(i,nsrf) = s_pblT(i)
3592	!!$ therm(i,nsrf) = s_therm(i)
3593	!!$ trmb1(i,nsrf) = s_trmb1(i)
3594	!!$ trmb2(i,nsrf) = s_trmb2(i)
3595	!!$ trmb3(i,nsrf) = s_trmb3(i)
3596	!!$ ENDIF
3597	!!$ ENDDO
3598	!!$ ENDDO
3599
3600
3601	DO i = 1, klon
3602	fder(i) = - 4.0RSIGMAzxtsol(i)**3
3603	ENDDO
3604
3605	zxqsurf(:) = 0.0
3606	zxsnow(:) = 0.0
3607	DO nsrf = 1, nbsrf
3608	DO i = 1, klon
3609	zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
3610	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
3611	ENDDO
3612	ENDDO
3613
3614	! Premier niveau de vent sortie dans physiq.F
3615	zu1(:) = u(:,1)
3616	zv1(:) = v(:,1)
3617
3618	END SUBROUTINE pbl_surface
3619	!
3620	!****************************************************************************************
3621	!
3622	SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
3623
3624	USE indice_sol_mod
3625
3626	INCLUDE "dimsoil.h"
3627
3628	! Ouput variables
3629	!****************************************************************************************
3630	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
3631	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
3632	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
3633	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
3634
3635
3636	!****************************************************************************************
3637	! Return module variables for writing to restart file
3638	!
3639	!****************************************************************************************
3640	fder_rst(:) = fder(:)
3641	snow_rst(:,:) = snow(:,:)
3642	qsurf_rst(:,:) = qsurf(:,:)
3643	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
3644
3645	!****************************************************************************************
3646	! Deallocate module variables
3647	!
3648	!****************************************************************************************
3649	! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
3650	IF (ALLOCATED(fder)) DEALLOCATE(fder)
3651	IF (ALLOCATED(snow)) DEALLOCATE(snow)
3652	IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
3653	IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
3654	IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
3655	IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
3656
3657	!jyg<
3658	!****************************************************************************************
3659	! Deallocate variables for pbl splitting
3660	!
3661	!****************************************************************************************
3662
3663	CALL wx_pbl_final
3664	!>jyg
3665
3666	END SUBROUTINE pbl_surface_final
3667	!
3668	!****************************************************************************************
3669	!
3670
3671	!albedo SB >>>
3672	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
3673	evap, z0m, z0h, agesno, &
3674	tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)
3675	!albedo SB <<<
3676	! Give default values where new fraction has appread
3677
3678	USE indice_sol_mod
3679	use phys_state_var_mod, only: delta_sal, ds_ns, dt_ns, delta_sst, dter, &
3680	dser, dt_ds
3681	use config_ocean_skin_m, only: activate_ocean_skin
3682
3683	INCLUDE "dimsoil.h"
3684	INCLUDE "clesphys.h"
3685	INCLUDE "compbl.h"
3686
3687	! Input variables
3688	!****************************************************************************************
3689	INTEGER, INTENT(IN) :: itime
3690	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
3691
3692	! InOutput variables
3693	!****************************************************************************************
3694	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
3695	!albedo SB >>>
3696	REAL, DIMENSION(klon,nsw,nbsrf), INTENT(INOUT) :: alb_dir, alb_dif
3697	INTEGER :: k
3698	!albedo SB <<<
3699	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m
3700	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: evap, agesno
3701	REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT) :: z0m,z0h
3702	REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
3703
3704	! Local variables
3705	!****************************************************************************************
3706	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
3707	CHARACTER(len=80) :: abort_message
3708	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
3709	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
3710	!
3711	! All at once !!
3712	!****************************************************************************************
3713
3714	DO nsrf = 1, nbsrf
3715	! First decide complement sub-surfaces
3716	SELECT CASE (nsrf)
3717	CASE(is_oce)
3718	nsrf_comp1=is_sic
3719	nsrf_comp2=is_ter
3720	nsrf_comp3=is_lic
3721	CASE(is_sic)
3722	nsrf_comp1=is_oce
3723	nsrf_comp2=is_ter
3724	nsrf_comp3=is_lic
3725	CASE(is_ter)
3726	nsrf_comp1=is_lic
3727	nsrf_comp2=is_oce
3728	nsrf_comp3=is_sic
3729	CASE(is_lic)
3730	nsrf_comp1=is_ter
3731	nsrf_comp2=is_oce
3732	nsrf_comp3=is_sic
3733	END SELECT
3734
3735	! Initialize all new fractions
3736	DO i=1, klon
3737	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
3738
3739	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
3740	! Use the complement sub-surface, keeping the continents unchanged
3741	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
3742	evap(i,nsrf) = evap(i,nsrf_comp1)
3743	z0m(i,nsrf) = z0m(i,nsrf_comp1)
3744	z0h(i,nsrf) = z0h(i,nsrf_comp1)
3745	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
3746	!albedo SB >>>
3747	DO k=1,nsw
3748	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp1)
3749	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp1)
3750	ENDDO
3751	!albedo SB <<<
3752	ustar(i,nsrf) = ustar(i,nsrf_comp1)
3753	u10m(i,nsrf) = u10m(i,nsrf_comp1)
3754	v10m(i,nsrf) = v10m(i,nsrf_comp1)
3755	IF (iflag_pbl > 1) THEN
3756	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
3757	ENDIF
3758	mfois(nsrf) = mfois(nsrf) + 1
3759	! F. Codron sensible default values for ocean and sea ice
3760	IF (nsrf.EQ.is_oce) THEN
3761	tsurf(i,nsrf) = 271.35
3762	! (temperature of sea water under sea ice, so that
3763	! is also the temperature of appearing sea water)
3764	DO k=1,nsw
3765	alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
3766	alb_dif(i,k,nsrf) = 0.06
3767	ENDDO
3768	if (activate_ocean_skin >= 1) then
3769	if (activate_ocean_skin == 2 &
3770	.and. type_ocean == "couple") then
3771	delta_sal(i) = 0.
3772	delta_sst(i) = 0.
3773	dter(i) = 0.
3774	dser(i) = 0.
3775	dt_ds(i) = 0.
3776	end if
3777
3778	ds_ns(i) = 0.
3779	dt_ns(i) = 0.
3780	end if
3781	ELSE IF (nsrf.EQ.is_sic) THEN
3782	tsurf(i,nsrf) = 271.35
3783	! (Temperature at base of sea ice. Surface
3784	! temperature could be higher, up to 0 Celsius
3785	! degrees. We set it to -1.8 Celsius degrees for
3786	! consistency with the ocean slab model.)
3787	DO k=1,nsw
3788	alb_dir(i,k,nsrf) = 0.3 ! thin ice
3789	alb_dif(i,k,nsrf) = 0.3
3790	ENDDO
3791	ENDIF
3792	ELSE
3793	! The continents have changed. The new fraction receives the mean sum of the existent fractions
3794	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3795	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3796	z0m(i,nsrf) = z0m(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0m(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3797	z0h(i,nsrf) = z0h(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0h(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3798	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3799	!albedo SB >>>
3800	DO k=1,nsw
3801	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3802	alb_dir(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3803	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3804	alb_dif(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3805	ENDDO
3806	!albedo SB <<<
3807	ustar(i,nsrf) = ustar(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3808	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3809	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3810	IF (iflag_pbl > 1) THEN
3811	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3812	ENDIF
3813
3814	! Security abort. This option has never been tested. To test, comment the following line.
3815	! abort_message='The fraction of the continents have changed!'
3816	! CALL abort_physic(modname,abort_message,1)
3817	nfois(nsrf) = nfois(nsrf) + 1
3818	ENDIF
3819	snow(i,nsrf) = 0.
3820	agesno(i,nsrf) = 0.
3821	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
3822	ELSE
3823	pfois(nsrf) = pfois(nsrf)+ 1
3824	ENDIF
3825	ENDDO
3826
3827	ENDDO
3828
3829	END SUBROUTINE pbl_surface_newfrac
3830	!
3831	!****************************************************************************************
3832	!
3833	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format