Context Navigation

pbl_surface_mod.F90 @ 4523

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

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

Download in other formats:

Original Format