Context Navigation

pbl_surface_mod.F90 @ 5015

Last change on this file since 5015 was 5015, checked in by musat, 2 months ago
Ajout flag ok_bug_zg_wk_pbl=y par defaut pour garder la convergence numerique.
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.3 KB

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

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

Download in other formats:

Original Format