Context Navigation

pbl_surface_mod.F90 @ 5408

Last change on this file since 5408 was 3838, checked in by musat, 4 years ago
Ajout ustar pour diagnostics pbl
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: 134.4 KB

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

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

Download in other formats:

Original Format