Context Navigation

pbl_surface_mod.F90 @ 3828

Last change on this file since 3828 was 3826, checked in by musat, 3 years ago
Nouveaux calculs 2m/10m activables par iflag_new_t2mq2m=1 Valeur par defaut iflag_new_t2mq2m=0 (<=> anciens calculs CMIP6)
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:keywords set to `Author Date Id Revision`
File size: 132.1 KB

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

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

Download in other formats:

Original Format