Context Navigation

pbl_surface_mod.F90 @ 5506

Last change on this file since 5506 was 5486, checked in by evignon, 7 days ago
inclusion d'un diagnostique de la sublimation de la glace sur les landice pour la conservation de l'eau
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.5 KB

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

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

Download in other formats:

Original Format