Context Navigation

pbl_surface_mod.F90 @ 2896

Last change on this file since 2896 was 2896, checked in by fhourdin, 7 years ago
Correction d'une erreur sur le calcul de T2m diagnostic. On passait deux fois la rugozité pour la quantité de mouvement z0m,z0m au lieu de z0m,z0h. L'ancienne version peut être réactivée avec iflag_pbl_surface_t2m_bug=1 Par défaut on met iflag_pbl_surface_t2m_bug=0 qui ne devrait faire perdre la convergence numérique qu'avec des versions récentes d'Orchidee différenciant les deux z0.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 127.8 KB

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

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

Download in other formats:

Original Format