Context Navigation

pbl_surface_mod.F90 @ 4143

Last change on this file since 4143 was 4143, checked in by dcugnet, 2 years ago

Some variables are renamed or replaced by direct equivalents:
- iso_indnum -> tracers(:)%iso_iName
- niso_possibles -> niso
- iqiso -> iqIsoPha ; index_trac -> itZonIso
- ok_iso_verif -> isoCheck
- ntraceurs_zone -> nzone ; ntraciso -> ntiso
- qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio
Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times.
Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18).
The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).

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

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

Download in other formats:

Original Format