Context Navigation

pbl_surface_mod.F90 @ 2236

Last change on this file since 2236 was 2234, checked in by fhourdin, 9 years ago
Correction d'un bug Bug fixing
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 127.9 KB

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

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

Download in other formats:

Original Format