Context Navigation

pbl_surface_mod.F90 @ 4648

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

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

Download in other formats:

Original Format