Context Navigation

pbl_surface_mod.F90 @ 4745

Last change on this file since 4745 was 4745, checked in by evignon, 6 months ago
nettoyage et corrections dans les routines atke pour utilisation en 3D (terre + mars)
File size: 172.7 KB

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

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

Download in other formats:

Original Format