Context Navigation

pbl_surface_mod.F90 @ 4657

Last change on this file since 4657 was 4482, checked in by lguez, 2 years ago
Sync latest trunk changes to branch LMDZ_ECRad
File size: 168.8 KB

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

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

Download in other formats:

Original Format