Context Navigation

pbl_surface_mod.F90 @ 5744

Last change on this file since 5744 was 5744, checked in by jyg, 6 days ago
Correction of a bug found by Catherine Rio in pbl_surface_mod.F90: modif = use of smallestreal in order to prevent division by 0.
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: 183.3 KB

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

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

Download in other formats:

Original Format