Context Navigation

pbl_surface_mod.F90 @ 5209

Last change on this file since 5209 was 4531, checked in by evignon, 18 months ago
petite reecriture pour neige soufflee dans pbl_surface
File size: 172.2 KB

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

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

Download in other formats:

Original Format