Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

pbl_surface_mod.F90

Last change on this file was 4916, checked in by evignon, 3 weeks ago
correction formulation de l'erosion de la neige soufflee suite aux investigations de Nicolas Chiabrando, prise en compte d'un temps d'erosion de la neige fraiche nouvellement accumulee
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 154.0 KB

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

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

Download in other formats:

Original Format