Context Navigation

pbl_surface_mod.F90 @ 4881

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

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

Download in other formats:

Original Format