Context Navigation

pbl_surface_mod.F90 @ 3184

Last change on this file since 3184 was 3179, checked in by jyg, 7 years ago

New structure for the representation of vdf
splitting in pbl_surface. Computations are
gathered in two subroutines included in the module
wx_pbl_mod: wx_pbl0_fuse, called before the
sub-surfaces, determines the single column
equivalent to the (w) and (x) columns.
wx_pbl0_split, called after the subsurfaces,
determines the distinct (w) and (x) surface
fluxes.

This is a first version with no surface
temperature difference between (w) and (x) (hence
the index 0).

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 130.3 KB

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

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

Download in other formats:

Original Format