Context Navigation

pbl_surface_mod.F90 @ 5101

Last change on this file since 5101 was 5101, checked in by abarral, 4 months ago
Handle DEBUG_IO in lmdz_cppkeys_wrapper.F90 Transform some files .F -> .[fF]90 [ne compile pas à cause de writefield_u non défini - en attente de réponse Laurent]
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: 173.1 KB

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

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

Download in other formats:

Original Format