Context Navigation

pbl_surface_mod.F90 @ 5279

Last change on this file since 5279 was 5274, checked in by abarral, 9 months ago
Replace yomcst.h by existing module
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: 175.1 KB

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

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

Download in other formats:

Original Format