Context Navigation

pbl_surface_mod.F90 @ 5144

Last change on this file since 5144 was 5144, checked in by abarral, 8 weeks ago
Put YOMCST.h into modules
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: 172.7 KB

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

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

Download in other formats:

Original Format