Context Navigation

pbl_surface_mod.F90 @ 5394

Last change on this file since 5394 was 2455, checked in by jyg, 9 years ago

Implementation of a second order distribution on

sub-surfaces of longwave net radiance (Alain
Lahellec).

Introduction of the flag iflag_order2_sollw:

=1 ==> second order distribution.
=0 ==> linear distribution (Default).

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: 126.1 KB

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

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

Download in other formats:

Original Format