Context Navigation

pbl_surface_mod.F90 @ 5134

Last change on this file since 5134 was 5084, checked in by Laurent Fairhead, 5 months ago
Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks AB, LF
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.0 KB

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

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

Download in other formats:

Original Format