Context Navigation

pbl_surface_mod.F90 @ 4741

Last change on this file since 4741 was 4727, checked in by idelkadi, 9 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
File size: 172.8 KB

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

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

Download in other formats:

Original Format