Context Navigation

pbl_surface_mod.F90 @ 5209

Last change on this file since 5209 was 5171, checked in by evignon, 2 months ago
adding a condition to prescribe surface temperature in 1D cases forced in flux this is needed for the new MOSAI 1D case to avoid spurious crashes
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 174.3 KB

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

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

Download in other formats:

Original Format