Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

pbl_surface_mod.F90 @ 3770

Last change on this file since 3770 was 3331, checked in by acozic, 6 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 148.1 KB

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

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

Download in other formats:

Original Format