Context Navigation

pbl_surface_mod.F90 @ 5273

Last change on this file since 5273 was 5273, checked in by abarral, 3 hours ago
Turn dimsoil.h into a module
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 174.3 KB

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

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

Download in other formats:

Original Format