Context Navigation

pbl_surface_mod.F90 @ 3720

Last change on this file since 3720 was 3720, checked in by lguez, 4 years ago

Output bulk sea surface salinity

Output bulk sea surface salinity when activate_ocean_skin >= 1. So
add variable sss in module phys_output_var_mod. In procedure
surf_ocean, promote sss from local variable to dummy argument.

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.2 KB

Line
1	!
2	! $Id: pbl_surface_mod.F90 3720 2020-06-12 22:23:29Z lguez $
3	!
4	MODULE pbl_surface_mod
5	!
6	! Planetary Boundary Layer and Surface module
7	!
8	! This module manages the calculation of turbulent diffusion in the boundary layer
9	! and all interactions towards the differents sub-surfaces.
10	!
11	!
12	USE dimphy
13	USE mod_phys_lmdz_para, ONLY : mpi_size
14	USE mod_grid_phy_lmdz, ONLY : klon_glo
15	USE ioipsl
16	USE surface_data, ONLY : type_ocean, ok_veget
17	USE surf_land_mod, ONLY : surf_land
18	USE surf_landice_mod, ONLY : surf_landice
19	USE surf_ocean_mod, ONLY : surf_ocean
20	USE surf_seaice_mod, ONLY : surf_seaice
21	USE cpl_mod, ONLY : gath2cpl
22	USE climb_hq_mod, ONLY : climb_hq_down, climb_hq_up
23	USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up
24	USE coef_diff_turb_mod, ONLY : coef_diff_turb
25	USE wx_pbl_mod, ONLY : wx_pbl_init, wx_pbl_final, &
26	!! wx_pbl_fuse_no_dts, wx_pbl_split_no_dts, &
27	!! wx_pbl_fuse, wx_pbl_split
28	wx_pbl0_fuse, wx_pbl0_split
29	use config_ocean_skin_m, only: activate_ocean_skin
30
31	IMPLICIT NONE
32
33	! Declaration of variables saved in restart file
34	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: fder ! flux drift
35	!$OMP THREADPRIVATE(fder)
36	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: snow ! snow at surface
37	!$OMP THREADPRIVATE(snow)
38	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qsurf ! humidity at surface
39	!$OMP THREADPRIVATE(qsurf)
40	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ftsoil ! soil temperature
41	!$OMP THREADPRIVATE(ftsoil)
42
43	INTEGER, SAVE :: iflag_pbl_surface_t2m_bug
44	!$OMP THREADPRIVATE(iflag_pbl_surface_t2m_bug)
45	!FC
46	! integer, save :: iflag_frein
47	! !$OMP THREADPRIVATE(iflag_frein)
48
49	CONTAINS
50	!
51	!****************************************************************************************
52	!
53	SUBROUTINE pbl_surface_init(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
54
55	! This routine should be called after the restart file has been read.
56	! This routine initialize the restart variables and does some validation tests
57	! for the index of the different surfaces and tests the choice of type of ocean.
58
59	USE indice_sol_mod
60	USE print_control_mod, ONLY: lunout
61	USE ioipsl_getin_p_mod, ONLY : getin_p
62	IMPLICIT NONE
63
64	INCLUDE "dimsoil.h"
65
66	! Input variables
67	!****************************************************************************************
68	REAL, DIMENSION(klon), INTENT(IN) :: fder_rst
69	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: snow_rst
70	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: qsurf_rst
71	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
72
73	! Local variables
74	!****************************************************************************************
75	INTEGER :: ierr
76	CHARACTER(len=80) :: abort_message
77	CHARACTER(len = 20) :: modname = 'pbl_surface_init'
78
79	!****************************************************************************************
80	! Allocate and initialize module variables with fields read from restart file.
81	!
82	!****************************************************************************************
83	ALLOCATE(fder(klon), stat=ierr)
84	IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
85
86	ALLOCATE(snow(klon,nbsrf), stat=ierr)
87	IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
88
89	ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
90	IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
91
92	ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
93	IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
94
95	fder(:) = fder_rst(:)
96	snow(:,:) = snow_rst(:,:)
97	qsurf(:,:) = qsurf_rst(:,:)
98	ftsoil(:,:,:) = ftsoil_rst(:,:,:)
99
100	!****************************************************************************************
101	! Test for sub-surface indices
102	!
103	!****************************************************************************************
104	IF (is_ter /= 1) THEN
105	WRITE(lunout,)" Warning *"
106	WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
107	WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
108	abort_message="voir ci-dessus"
109	CALL abort_physic(modname,abort_message,1)
110	ENDIF
111
112	IF ( is_oce > is_sic ) THEN
113	WRITE(lunout,)' Warning *'
114	WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
115	WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
116	WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
117	abort_message='voir ci-dessus'
118	CALL abort_physic(modname,abort_message,1)
119	ENDIF
120
121	IF ( is_lic > is_sic ) THEN
122	WRITE(lunout,)' Warning *'
123	WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
124	WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
125	WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
126	abort_message='voir ci-dessus'
127	CALL abort_physic(modname,abort_message,1)
128	ENDIF
129
130	!****************************************************************************************
131	! Validation of ocean mode
132	!
133	!****************************************************************************************
134
135	IF (type_ocean /= 'slab ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN
136	WRITE(lunout,)' Warning *'
137	WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean
138	abort_message='option pour l''ocean non valable'
139	CALL abort_physic(modname,abort_message,1)
140	ENDIF
141
142	iflag_pbl_surface_t2m_bug=0
143	CALL getin_p('iflag_pbl_surface_t2m_bug',iflag_pbl_surface_t2m_bug)
144	!FC
145	! iflag_frein = 0
146	! CALL getin_p('iflag_frein',iflag_frein)
147	!
148	!jyg<
149	!****************************************************************************************
150	! Allocate variables for pbl splitting
151	!
152	!****************************************************************************************
153
154	CALL wx_pbl_init
155	!>jyg
156
157	END SUBROUTINE pbl_surface_init
158	!
159	!****************************************************************************************
160	!
161
162	SUBROUTINE pbl_surface( &
163	dtime, date0, itap, jour, &
164	debut, lafin, &
165	rlon, rlat, rugoro, rmu0, &
166	zsig, lwdown_m, pphi, cldt, &
167	rain_f, snow_f, solsw_m, sollw_m, &
168	gustiness, &
169	t, q, u, v, &
170	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
171	!! t_x, q_x, t_w, q_w, &
172	wake_dlt, wake_dlq, &
173	wake_cstar, wake_s, &
174	!!!
175	pplay, paprs, pctsrf, &
176	ts,SFRWL, alb_dir, alb_dif,ustar, u10m, v10m,wstar, &
177	cdragh, cdragm, zu1, zv1, &
178	alb_dir_m, alb_dif_m, zxsens, zxevap, &
179	alb3_lic, runoff, snowhgt, qsnow, to_ice, sissnow, &
180	zxtsol, zxfluxlat, zt2m, qsat2m, &
181	d_t, d_q, d_u, d_v, d_t_diss, &
182	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
183	d_t_w, d_q_w, &
184	d_t_x, d_q_x, &
185	!! d_wake_dlt,d_wake_dlq, &
186	zxsens_x, zxfluxlat_x,zxsens_w,zxfluxlat_w, &
187	!!!
188	!!! nrlmd le 13/06/2011
189	delta_tsurf,wake_dens,cdragh_x,cdragh_w, &
190	cdragm_x,cdragm_w,kh,kh_x,kh_w, &
191	!!!
192	zcoefh, zcoefm, slab_wfbils, &
193	qsol, zq2m, s_pblh, s_plcl, &
194	!!!
195	!!! jyg le 08/02/2012
196	s_pblh_x, s_plcl_x, s_pblh_w, s_plcl_w, &
197	!!!
198	s_capCL, s_oliqCL, s_cteiCL, s_pblT, &
199	s_therm, s_trmb1, s_trmb2, s_trmb3, &
200	zustar,zu10m, zv10m, fder_print, &
201	zxqsurf, rh2m, zxfluxu, zxfluxv, &
202	z0m, z0h, agesno, sollw, solsw, &
203	d_ts, evap, fluxlat, t2m, &
204	wfbils, wfbilo, wfevap, wfrain, wfsnow, &
205	flux_t, flux_u, flux_v, &
206	dflux_t, dflux_q, zxsnow, &
207	!jyg<
208	!! zxfluxt, zxfluxq, q2m, flux_q, tke, &
209	zxfluxt, zxfluxq, q2m, flux_q, tke_x, &
210	!>jyg
211	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
212	!! tke_x, tke_w &
213	wake_dltke, &
214	treedrg &
215	!FC
216	!!!
217	)
218	!****************************************************************************************
219	! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
220	! Objet: interface de "couche limite" (diffusion verticale)
221	!
222	!AA REM:
223	!AA-----
224	!AA Tout ce qui a trait au traceurs est dans phytrac maintenant
225	!AA pour l'instant le calcul de la couche limite pour les traceurs
226	!AA se fait avec cltrac et ne tient pas compte de la differentiation
227	!AA des sous-fraction de sol.
228	!AA REM bis :
229	!AA----------
230	!AA Pour pouvoir extraire les coefficient d'echanges et le vent
231	!AA dans la premiere couche, 3 champs supplementaires ont ete crees
232	!AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
233	!AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
234	!AA si les informations des subsurfaces doivent etre prises en compte
235	!AA il faudra sortir ces memes champs en leur ajoutant une dimension,
236	!AA c'est a dire nbsrf (nbre de subsurface).
237	!
238	! Arguments:
239	!
240	! dtime----input-R- interval du temps (secondes)
241	! itap-----input-I- numero du pas de temps
242	! date0----input-R- jour initial
243	! t--------input-R- temperature (K)
244	! q--------input-R- vapeur d'eau (kg/kg)
245	! u--------input-R- vitesse u
246	! v--------input-R- vitesse v
247	! wake_dlt-input-R- temperatre difference between (w) and (x) (K)
248	! wake_dlq-input-R- humidity difference between (w) and (x) (kg/kg)
249	!wake_cstar-input-R- wake gust front speed (m/s)
250	! wake_s---input-R- wake fractionnal area
251	! ts-------input-R- temperature du sol (en Kelvin)
252	! paprs----input-R- pression a intercouche (Pa)
253	! pplay----input-R- pression au milieu de couche (Pa)
254	! rlat-----input-R- latitude en degree
255	! z0m, z0h ----input-R- longeur de rugosite (en m)
256	! Martin
257	! zsig-----input-R- slope
258	! cldt-----input-R- total cloud fraction
259	! pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol)
260	! Martin
261	!
262	! d_t------output-R- le changement pour "t"
263	! d_q------output-R- le changement pour "q"
264	! d_u------output-R- le changement pour "u"
265	! d_v------output-R- le changement pour "v"
266	! d_ts-----output-R- le changement pour "ts"
267	! flux_t---output-R- flux de chaleur sensible (CpT) J/m2/s (W/m2)
268	! (orientation positive vers le bas)
269	! tke_x---input/output-R- tke in the (x) region (kg/m**2/s)
270	! wake_dltke-input/output-R- tke difference between (w) and (x) (kg/m**2/s)
271	! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
272	! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
273	! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
274	! dflux_t--output-R- derive du flux sensible
275	! dflux_q--output-R- derive du flux latent
276	! zu1------output-R- le vent dans la premiere couche
277	! zv1------output-R- le vent dans la premiere couche
278	! trmb1----output-R- deep_cape
279	! trmb2----output-R- inhibition
280	! trmb3----output-R- Point Omega
281	! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL
282	! plcl-----output-R- Niveau de condensation
283	! pblh-----output-R- HCL
284	! pblT-----output-R- T au nveau HCL
285	! treedrg--output-R- tree drag (m)
286	!
287	USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, carbon_cycle_tr, level_coupling_esm
288	USE carbon_cycle_mod, ONLY : co2_send, nbcf_out, fields_out, yfields_out, cfname_out
289	USE indice_sol_mod
290	USE time_phylmdz_mod, ONLY : day_ini,annee_ref,itau_phy
291	USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat, grid1dto2d_glo
292	USE print_control_mod, ONLY : prt_level,lunout
293	USE ioipsl_getin_p_mod, ONLY : getin_p
294	use phys_state_var_mod, only: s_int, ds_ns, dt_ns
295	use phys_output_var_mod, only: sst_nff, dter, dser, tkt, tks, taur, sss
296	#ifdef CPP_XIOS
297	USE wxios, ONLY: missing_val
298	#else
299	use netcdf, only: missing_val => nf90_fill_real
300	#endif
301
302	IMPLICIT NONE
303
304	INCLUDE "dimsoil.h"
305	INCLUDE "YOMCST.h"
306	INCLUDE "YOETHF.h"
307	INCLUDE "FCTTRE.h"
308	INCLUDE "clesphys.h"
309	INCLUDE "compbl.h"
310	INCLUDE "flux_arp.h"
311	!FC
312	INCLUDE "dimpft.h"
313
314	!****************************************************************************************
315	REAL, INTENT(IN) :: dtime ! time interval (s)
316	REAL, INTENT(IN) :: date0 ! initial day
317	INTEGER, INTENT(IN) :: itap ! time step
318	INTEGER, INTENT(IN) :: jour ! current day of the year
319	LOGICAL, INTENT(IN) :: debut ! true if first run step
320	LOGICAL, INTENT(IN) :: lafin ! true if last run step
321	REAL, DIMENSION(klon), INTENT(IN) :: rlon ! longitudes in degrees
322	REAL, DIMENSION(klon), INTENT(IN) :: rlat ! latitudes in degrees
323	REAL, DIMENSION(klon), INTENT(IN) :: rugoro ! rugosity length
324	REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cosine of solar zenith angle
325	REAL, DIMENSION(klon), INTENT(IN) :: rain_f ! rain fall
326	REAL, DIMENSION(klon), INTENT(IN) :: snow_f ! snow fall
327	REAL, DIMENSION(klon), INTENT(IN) :: solsw_m ! net shortwave radiation at mean surface
328	REAL, DIMENSION(klon), INTENT(IN) :: sollw_m ! net longwave radiation at mean surface
329	REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K)
330	REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! water vapour (kg/kg)
331	REAL, DIMENSION(klon,klev), INTENT(IN) :: u ! u speed
332	REAL, DIMENSION(klon,klev), INTENT(IN) :: v ! v speed
333	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pression (Pa)
334	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression between layers (Pa)
335	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! sub-surface fraction
336	! Martin
337	REAL, DIMENSION(klon), INTENT(IN) :: zsig ! slope
338	REAL, DIMENSION(klon), INTENT(IN) :: lwdown_m ! downward longwave radiation at mean s
339	REAL, DIMENSION(klon), INTENT(IN) :: gustiness ! gustiness
340
341	REAL, DIMENSION(klon), INTENT(IN) :: cldt ! total cloud fraction
342	REAL, DIMENSION(klon,klev), INTENT(IN) :: pphi ! geopotential (m2/s2)
343	! Martin
344
345	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
346	!! REAL, DIMENSION(klon,klev), INTENT(IN) :: t_x ! Temp\'erature hors poche froide
347	!! REAL, DIMENSION(klon,klev), INTENT(IN) :: t_w ! Temp\'erature dans la poches froide
348	!! REAL, DIMENSION(klon,klev), INTENT(IN) :: q_x !
349	!! REAL, DIMENSION(klon,klev), INTENT(IN) :: q_w ! Pareil pour l'humidit\'e
350	REAL, DIMENSION(klon,klev), INTENT(IN) :: wake_dlt !temperature difference between (w) and (x) (K)
351	REAL, DIMENSION(klon,klev), INTENT(IN) :: wake_dlq !humidity difference between (w) and (x) (K)
352	REAL, DIMENSION(klon), INTENT(IN) :: wake_s ! Fraction de poches froides
353	REAL, DIMENSION(klon), INTENT(IN) :: wake_cstar! Vitesse d'expansion des poches froides
354	REAL, DIMENSION(klon), INTENT(IN) :: wake_dens
355	!!!
356
357	! Input/Output variables
358	!****************************************************************************************
359	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ts ! temperature at surface (K)
360	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: delta_tsurf !surface temperature difference between
361	!wake and off-wake regions
362	!albedo SB >>>
363	REAL, DIMENSIOn(6),intent(in) :: SFRWL
364	REAL, DIMENSION(klon, nsw, nbsrf), INTENT(INOUT) :: alb_dir,alb_dif
365	!albedo SB <<<
366	!jyg Pourquoi ustar et wstar sont-elles INOUT ?
367	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ustar ! u* (m/s)
368	REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT) :: wstar ! w* (m/s)
369	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: u10m ! u speed at 10m
370	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: v10m ! v speed at 10m
371	!jyg<
372	!! REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke
373	REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke_x
374	!>jyg
375
376	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
377	REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: wake_dltke ! TKE_w - TKE_x
378	!!!
379
380	! Output variables
381	!****************************************************************************************
382	REAL, DIMENSION(klon), INTENT(OUT) :: cdragh ! drag coefficient for T and Q
383	REAL, DIMENSION(klon), INTENT(OUT) :: cdragm ! drag coefficient for wind
384	REAL, DIMENSION(klon), INTENT(OUT) :: zu1 ! u wind speed in first layer
385	REAL, DIMENSION(klon), INTENT(OUT) :: zv1 ! v wind speed in first layer
386	!albedo SB >>>
387	REAL, DIMENSION(klon, nsw), INTENT(OUT) :: alb_dir_m,alb_dif_m
388	!albedo SB <<<
389	! Martin
390	REAL, DIMENSION(klon), INTENT(OUT) :: alb3_lic
391	! Martin
392	REAL, DIMENSION(klon), INTENT(OUT) :: zxsens ! sensible heat flux at surface with inversed sign
393	! (=> positive sign upwards)
394	REAL, DIMENSION(klon), INTENT(OUT) :: zxevap ! water vapour flux at surface, positiv upwards
395	REAL, DIMENSION(klon), INTENT(OUT) :: zxtsol ! temperature at surface, mean for each grid point
396	!!! jyg le ???
397	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t_w ! !
398	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q_w ! ! Tendances dans les poches
399	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t_x ! !
400	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q_x ! ! Tendances hors des poches
401	!!! jyg
402	REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat ! latent flux, mean for each grid point
403	REAL, DIMENSION(klon), INTENT(OUT) :: zt2m ! temperature at 2m, mean for each grid point
404	REAL, DIMENSION(klon), INTENT(OUT) :: qsat2m
405	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! change in temperature
406	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t_diss ! change in temperature
407	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_q ! change in water vapour
408	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u ! change in u speed
409	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v ! change in v speed
410
411	REAL, INTENT(OUT):: zcoefh(:, :, :) ! (klon, klev, nbsrf + 1)
412	! coef for turbulent diffusion of T and Q, mean for each grid point
413
414	REAL, INTENT(OUT):: zcoefm(:, :, :) ! (klon, klev, nbsrf + 1)
415	! coef for turbulent diffusion of U and V (?), mean for each grid point
416
417	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
418	REAL, DIMENSION(klon), INTENT(OUT) :: zxsens_x ! Flux sensible hors poche
419	REAL, DIMENSION(klon), INTENT(OUT) :: zxsens_w ! Flux sensible dans la poche
420	REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat_x! Flux latent hors poche
421	REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat_w! Flux latent dans la poche
422	!! REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_wake_dlt
423	!! REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_wake_dlq
424
425	! Output only for diagnostics
426	REAL, DIMENSION(klon), INTENT(OUT) :: cdragh_x
427	REAL, DIMENSION(klon), INTENT(OUT) :: cdragh_w
428	REAL, DIMENSION(klon), INTENT(OUT) :: cdragm_x
429	REAL, DIMENSION(klon), INTENT(OUT) :: cdragm_w
430	REAL, DIMENSION(klon), INTENT(OUT) :: kh
431	REAL, DIMENSION(klon), INTENT(OUT) :: kh_x
432	REAL, DIMENSION(klon), INTENT(OUT) :: kh_w
433	!!!
434	REAL, DIMENSION(klon), INTENT(OUT) :: slab_wfbils! heat balance at surface only for slab at ocean points
435	REAL, DIMENSION(klon), INTENT(OUT) :: qsol ! water height in the soil (mm)
436	REAL, DIMENSION(klon), INTENT(OUT) :: zq2m ! water vapour at 2m, mean for each grid point
437	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh ! height of the planetary boundary layer(HPBL)
438	!!! jyg le 08/02/2012
439	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh_x ! height of the PBL in the off-wake region
440	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh_w ! height of the PBL in the wake region
441	!!!
442	REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl ! condensation level
443	!!! jyg le 08/02/2012
444	REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl_x ! condensation level in the off-wake region
445	REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl_w ! condensation level in the wake region
446	!!!
447	REAL, DIMENSION(klon), INTENT(OUT) :: s_capCL ! CAPE of PBL
448	REAL, DIMENSION(klon), INTENT(OUT) :: s_oliqCL ! liquid water intergral of PBL
449	REAL, DIMENSION(klon), INTENT(OUT) :: s_cteiCL ! cloud top instab. crit. of PBL
450	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblT ! temperature at PBLH
451	REAL, DIMENSION(klon), INTENT(OUT) :: s_therm ! thermal virtual temperature excess
452	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb1 ! deep cape, mean for each grid point
453	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb2 ! inhibition, mean for each grid point
454	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb3 ! point Omega, mean for each grid point
455	REAL, DIMENSION(klon), INTENT(OUT) :: zustar ! u*
456	REAL, DIMENSION(klon), INTENT(OUT) :: zu10m ! u speed at 10m, mean for each grid point
457	REAL, DIMENSION(klon), INTENT(OUT) :: zv10m ! v speed at 10m, mean for each grid point
458	REAL, DIMENSION(klon), INTENT(OUT) :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
459	REAL, DIMENSION(klon), INTENT(OUT) :: zxqsurf ! humidity at surface, mean for each grid point
460	REAL, DIMENSION(klon), INTENT(OUT) :: rh2m ! relative humidity at 2m
461	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxu ! u wind tension, mean for each grid point
462	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxv ! v wind tension, mean for each grid point
463	REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT) :: z0m,z0h ! rugosity length (m)
464	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: agesno ! age of snow at surface
465	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: solsw ! net shortwave radiation at surface
466	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: sollw ! net longwave radiation at surface
467	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: d_ts ! change in temperature at surface
468	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: evap ! evaporation at surface
469	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: fluxlat ! latent flux
470	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: t2m ! temperature at 2 meter height
471	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbils ! heat balance at surface
472	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbilo ! water balance at surface weighted by srf
473	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfevap ! water balance (evap) at surface weighted by srf
474	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfrain ! water balance (rain) at surface weighted by srf
475	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfsnow ! water balance (snow) at surface weighted by srf
476	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t ! sensible heat flux (CpT) J/m2/s (W/m2)
477	! positve orientation downwards
478	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u ! u wind tension (kg m/s)/(m**2 s) or Pascal
479	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v ! v wind tension (kg m/s)/(m**2 s) or Pascal
480	!FC
481	REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: treedrg ! tree drag (m)
482
483
484	! Output not needed
485	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_t ! change of sensible heat flux
486	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_q ! change of water vapour flux
487	REAL, DIMENSION(klon), INTENT(OUT) :: zxsnow ! snow at surface, mean for each grid point
488	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxt ! sensible heat flux, mean for each grid point
489	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxq ! water vapour flux, mean for each grid point
490	REAL, DIMENSION(klon, nbsrf),INTENT(OUT) :: q2m ! water vapour at 2 meter height
491	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q ! water vapour flux(latent flux) (kg/m**2/s)
492
493	! Martin
494	! sisvat
495	REAL, DIMENSION(klon), INTENT(OUT) :: qsnow ! snow water content
496	REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt ! snow height
497	REAL, DIMENSION(klon), INTENT(OUT) :: to_ice ! snow passed to ice
498	REAL, DIMENSION(klon), INTENT(OUT) :: sissnow ! snow in snow model
499	REAL, DIMENSION(klon), INTENT(OUT) :: runoff ! runoff on land ice
500	! Martin
501
502	! Local variables with attribute SAVE
503	!****************************************************************************************
504	INTEGER, SAVE :: nhoridbg, nidbg ! variables for IOIPSL
505	!$OMP THREADPRIVATE(nhoridbg, nidbg)
506	LOGICAL, SAVE :: debugindex=.FALSE.
507	!$OMP THREADPRIVATE(debugindex)
508	LOGICAL, SAVE :: first_call=.TRUE.
509	!$OMP THREADPRIVATE(first_call)
510	CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
511	!$OMP THREADPRIVATE(cl_surf)
512	REAL, SAVE :: beta_land ! beta for wx_dts
513	!$OMP THREADPRIVATE(beta_land)
514
515	! Other local variables
516	!****************************************************************************************
517	! >> PC
518	INTEGER :: ierr
519	INTEGER :: n
520	! << PC
521	INTEGER :: iflag_split
522	INTEGER :: i, k, nsrf
523	INTEGER :: knon, j
524	INTEGER :: idayref
525	INTEGER , DIMENSION(klon) :: ni
526	REAL :: yt1_new
527	REAL :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
528	REAL :: amn, amx
529	REAL :: f1 ! fraction de longeurs visibles parmi tout SW intervalle
530	REAL, DIMENSION(klon) :: r_co2_ppm ! taux CO2 atmosphere
531	REAL, DIMENSION(klon) :: yts, yz0m, yz0h, ypct
532	!albedo SB >>>
533	REAL, DIMENSION(klon) :: yalb,yalb_vis
534	!albedo SB <<<
535	REAL, DIMENSION(klon) :: yt1, yq1, yu1, yv1
536	REAL, DIMENSION(klon) :: ysnow, yqsurf, yagesno, yqsol
537	REAL, DIMENSION(klon) :: yrain_f, ysnow_f
538	REAL, DIMENSION(klon) :: ysolsw, ysollw
539	REAL, DIMENSION(klon) :: yfder
540	REAL, DIMENSION(klon) :: yrugoro
541	REAL, DIMENSION(klon) :: yfluxlat
542	REAL, DIMENSION(klon) :: y_d_ts
543	REAL, DIMENSION(klon) :: y_flux_t1, y_flux_q1
544	REAL, DIMENSION(klon) :: y_dflux_t, y_dflux_q
545	REAL, DIMENSION(klon) :: y_flux_u1, y_flux_v1
546	REAL, DIMENSION(klon) :: yt2m, yq2m, yu10m
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	real, DIMENSION(klon):: ysst_nff, ys_int, yds_ns, ydt_ns, ydter, ydser, &
833	ytkt, ytks, ytaur, ysss
834	! compression of sst_nff, s_int, ds_ns, dt_ns, dter, dser, tkt, tks,
835	! taur, sss on ocean points
836
837	!****************************************************************************************
838	! End of declarations
839	!****************************************************************************************
840
841	IF (prt_level >=10) print *,' -> pbl_surface, itap ',itap
842	!
843	!!jyg iflag_split = mod(iflag_pbl_split,2)
844	iflag_split = mod(iflag_pbl_split,10)
845
846	!****************************************************************************************
847	! 1) Initialisation and validation tests
848	! Only done first time entering this subroutine
849	!
850	!****************************************************************************************
851
852	IF (first_call) THEN
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	d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
919	zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
920	zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
921	cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0.
922	kh(:)=0. ; kh_x(:)=0. ; kh_w(:)=0.
923	slab_wfbils(:)=0.
924	s_pblh(:)=0. ; s_pblh_x(:)=0. ; s_pblh_w(:)=0.
925	s_plcl(:)=0. ; s_plcl_x(:)=0. ; s_plcl_w(:)=0.
926	s_capCL(:)=0. ; s_oliqCL(:)=0. ; s_cteiCL(:)=0. ; s_pblT(:)=0.
927	s_therm(:)=0.
928	s_trmb1(:)=0. ; s_trmb2(:)=0. ; s_trmb3(:)=0.
929	zustar(:)=0.
930	zu10m(:)=0. ; zv10m(:)=0.
931	fder_print(:)=0.
932	zxqsurf(:)=0.
933	zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
934	solsw(:,:)=0. ; sollw(:,:)=0.
935	d_ts(:,:)=0.
936	evap(:,:)=0.
937	fluxlat(:,:)=0.
938	wfbils(:,:)=0. ; wfbilo(:,:)=0.
939	wfevap(:,:)=0. ; wfrain(:,:)=0. ; wfsnow(:,:)=0.
940	flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.
941	dflux_t(:)=0. ; dflux_q(:)=0.
942	zxsnow(:)=0.
943	zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.
944	qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
945	runoff(:)=0.
946	IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
947	zcoefh(:,:,:) = 0.0
948	zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used
949	zcoefm(:,:,:) = 0.0
950	zcoefm(:,1,:) = 999999. !
951	ELSE
952	zcoefm(:,:,is_ave)=0.
953	zcoefh(:,:,is_ave)=0.
954	ENDIF
955	!!
956	! The components "is_ave" of tke_x and wake_deltke are "OUT" variables
957	!jyg<
958	!! tke(:,:,is_ave)=0.
959	tke_x(:,:,is_ave)=0.
960	wake_dltke(:,:,is_ave)=0.
961	!>jyg
962	!!! jyg le 23/02/2013
963	t2m(:,:) = 999999. ! t2m and q2m are meaningfull only over sub-surfaces
964	q2m(:,:) = 999999. ! actually present in the grid cell.
965	!!!
966	rh2m(:) = 0. ; qsat2m(:) = 0.
967	!!!
968	!!! jyg le 10/02/2012
969	rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0.
970
971	! 2b) Initialization of all local variables that will be compressed later
972	!****************************************************************************************
973	!! cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0
974	ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0
975	!! zv1 = 0.0 ; yqsurf = 0.0
976	!albedo SB >>>
977	yqsurf = 0.0 ; yalb = 0.0 ; yalb_vis = 0.0
978	!albedo SB <<<
979	yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0
980	ysollw = 0.0 ; yz0m = 0.0 ; yz0h = 0.0 ; yu1 = 0.0
981	yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0
982	ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0
983	yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
984	yrugoro = 0.0 ; ywindsp = 0.0
985	!! d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0
986	yfluxlat=0.0
987	!! flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0
988	!! d_t_diss= 0.0 ;d_u = 0.0 ; d_v = 0.0
989	yqsol = 0.0
990	ytherm = 0.0 ; ytke=0.
991	!FC
992	y_treedrg=0.
993
994	! Martin
995	ysnowhgt = 0.0; yqsnow = 0.0 ; yrunoff = 0.0 ; ytoice =0.0
996	yalb3_new = 0.0 ; ysissnow = 0.0
997	ypphi = 0.0 ; ycldt = 0.0 ; yrmu0 = 0.0
998	! Martin
999
1000	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
1001	ytke_x=0. ; ytke_w=0. ; ywake_dltke=0.
1002	y_d_t_x=0. ; y_d_t_w=0. ; y_d_q_x=0. ; y_d_q_w=0.
1003	!! d_t_w=0. ; d_q_w=0.
1004	!! d_t_x=0. ; d_q_x=0.
1005	!! d_wake_dlt=0. ; d_wake_dlq=0.
1006	yfluxlat_x=0. ; yfluxlat_w=0.
1007	ywake_s=0. ; ywake_cstar=0. ;ywake_dens=0.
1008	!!!
1009	!!! nrlmd le 13/06/2011
1010	tau_eq=0. ; delta_coef=0.
1011	y_delta_flux_t1=0.
1012	ydtsurf_th=0.
1013	yts_x=0. ; yts_w=0.
1014	y_delta_tsurf=0.
1015	!!!
1016	ytsoil = 999999.
1017	!FC
1018	y_d_u_frein(:,:)=0.
1019	y_d_v_frein(:,:)=0.
1020	!FC
1021
1022	! >> PC
1023	!the yfields_out variable is defined in (klon,nbcf_out) even if it is used on
1024	!the ORCHIDEE grid and as such should be defined in yfields_out(knon,nbcf_out) but
1025	!the knon variable is not known at that level of pbl_surface_mod
1026
1027	!the yfields_in variable is defined in (klon,nbcf_in) even if it is used on the
1028	!ORCHIDEE grid and as such should be defined in yfields_in(knon,nbcf_in) but the
1029	!knon variable is not known at that level of pbl_surface_mod
1030
1031	yfields_out(:,:) = 0.
1032	! << PC
1033
1034
1035	! 2c) Initialization of all local variables computed within the subsurface loop and used later on
1036	!****************************************************************************************
1037	d_t_diss_x(:,:) = 0. ; d_t_diss_w(:,:) = 0.
1038	d_u_x(:,:)=0. ; d_u_w(:,:)=0.
1039	d_v_x(:,:)=0. ; d_v_w(:,:)=0.
1040	flux_t_x(:,:,:)=0. ; flux_t_w(:,:,:)=0.
1041	flux_q_x(:,:,:)=0. ; flux_q_w(:,:,:)=0.
1042	!
1043	!jyg<
1044	flux_u_x(:,:,:)=0. ; flux_u_w(:,:,:)=0.
1045	flux_v_x(:,:,:)=0. ; flux_v_w(:,:,:)=0.
1046	fluxlat_x(:,:)=0. ; fluxlat_w(:,:)=0.
1047	!>jyg
1048	!
1049	!jyg<
1050	! pblh,plcl,capCL,cteiCL ... are meaningfull only over sub-surfaces
1051	! actually present in the grid cell ==> value set to 999999.
1052	!
1053	!jyg<
1054	ustar(:,:) = 999999.
1055	wstar(:,:) = 999999.
1056	windsp(:,:) = SQRT(u10m(:,:)2 + v10m(:,:)2 )
1057	u10m(:,:) = 999999.
1058	v10m(:,:) = 999999.
1059	!>jyg
1060	!
1061	pblh(:,:) = 999999. ! Hauteur de couche limite
1062	plcl(:,:) = 999999. ! Niveau de condensation de la CLA
1063	capCL(:,:) = 999999. ! CAPE de couche limite
1064	oliqCL(:,:) = 999999. ! eau_liqu integree de couche limite
1065	cteiCL(:,:) = 999999. ! cloud top instab. crit. couche limite
1066	pblt(:,:) = 999999. ! T a la Hauteur de couche limite
1067	therm(:,:) = 999999.
1068	trmb1(:,:) = 999999. ! deep_cape
1069	trmb2(:,:) = 999999. ! inhibition
1070	trmb3(:,:) = 999999. ! Point Omega
1071	!
1072	t2m_x(:,:) = 999999.
1073	q2m_x(:,:) = 999999.
1074	ustar_x(:,:) = 999999.
1075	wstar_x(:,:) = 999999.
1076	u10m_x(:,:) = 999999.
1077	v10m_x(:,:) = 999999.
1078	!
1079	pblh_x(:,:) = 999999. ! Hauteur de couche limite
1080	plcl_x(:,:) = 999999. ! Niveau de condensation de la CLA
1081	capCL_x(:,:) = 999999. ! CAPE de couche limite
1082	oliqCL_x(:,:) = 999999. ! eau_liqu integree de couche limite
1083	cteiCL_x(:,:) = 999999. ! cloud top instab. crit. couche limite
1084	pblt_x(:,:) = 999999. ! T a la Hauteur de couche limite
1085	therm_x(:,:) = 999999.
1086	trmb1_x(:,:) = 999999. ! deep_cape
1087	trmb2_x(:,:) = 999999. ! inhibition
1088	trmb3_x(:,:) = 999999. ! Point Omega
1089	!
1090	t2m_w(:,:) = 999999.
1091	q2m_w(:,:) = 999999.
1092	ustar_w(:,:) = 999999.
1093	wstar_w(:,:) = 999999.
1094	u10m_w(:,:) = 999999.
1095	v10m_w(:,:) = 999999.
1096
1097	pblh_w(:,:) = 999999. ! Hauteur de couche limite
1098	plcl_w(:,:) = 999999. ! Niveau de condensation de la CLA
1099	capCL_w(:,:) = 999999. ! CAPE de couche limite
1100	oliqCL_w(:,:) = 999999. ! eau_liqu integree de couche limite
1101	cteiCL_w(:,:) = 999999. ! cloud top instab. crit. couche limite
1102	pblt_w(:,:) = 999999. ! T a la Hauteur de couche limite
1103	therm_w(:,:) = 999999.
1104	trmb1_w(:,:) = 999999. ! deep_cape
1105	trmb2_w(:,:) = 999999. ! inhibition
1106	trmb3_w(:,:) = 999999. ! Point Omega
1107	!!!
1108	!
1109	!!!
1110	!****************************************************************************************
1111	! 3) - Calculate pressure thickness of each layer
1112	! - Calculate the wind at first layer
1113	! - Mean calculations of albedo
1114	! - Calculate net radiance at sub-surface
1115	!****************************************************************************************
1116	DO k = 1, klev
1117	DO i = 1, klon
1118	delp(i,k) = paprs(i,k)-paprs(i,k+1)
1119	ENDDO
1120	ENDDO
1121
1122	!****************************************************************************************
1123	! Test for rugos........ from physiq.. A la fin plutot???
1124	!
1125	!****************************************************************************************
1126
1127	DO nsrf = 1, nbsrf
1128	DO i = 1, klon
1129	z0m(i,nsrf) = MAX(z0m(i,nsrf),z0min)
1130	z0h(i,nsrf) = MAX(z0h(i,nsrf),z0min)
1131	ENDDO
1132	ENDDO
1133
1134	! Mean calculations of albedo
1135	!
1136	! * alb : mean albedo for whole SW interval
1137	!
1138	! Mean albedo for grid point
1139	! * alb_m : mean albedo at whole SW interval
1140
1141	alb_dir_m(:,:) = 0.0
1142	alb_dif_m(:,:) = 0.0
1143	DO k = 1, nsw
1144	DO nsrf = 1, nbsrf
1145	DO i = 1, klon
1146	alb_dir_m(i,k) = alb_dir_m(i,k) + alb_dir(i,k,nsrf) * pctsrf(i,nsrf)
1147	alb_dif_m(i,k) = alb_dif_m(i,k) + alb_dif(i,k,nsrf) * pctsrf(i,nsrf)
1148	ENDDO
1149	ENDDO
1150	ENDDO
1151
1152	! We here suppose the fraction f1 of incoming radiance of visible radiance
1153	! as a fraction of all shortwave radiance
1154	f1 = 0.5
1155	! f1 = 1 ! put f1=1 to recreate old calculations
1156
1157	!f1 is already included with SFRWL values in each surf files
1158	alb=0.0
1159	DO k=1,nsw
1160	DO nsrf = 1, nbsrf
1161	DO i = 1, klon
1162	alb(i,nsrf) = alb(i,nsrf) + alb_dir(i,k,nsrf)*SFRWL(k)
1163	ENDDO
1164	ENDDO
1165	ENDDO
1166
1167	alb_m=0.0
1168	DO k = 1,nsw
1169	DO i = 1, klon
1170	alb_m(i) = alb_m(i) + alb_dir_m(i,k)*SFRWL(k)
1171	ENDDO
1172	ENDDO
1173	!albedo SB <<<
1174
1175
1176
1177	! Calculation of mean temperature at surface grid points
1178	ztsol(:) = 0.0
1179	DO nsrf = 1, nbsrf
1180	DO i = 1, klon
1181	ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
1182	ENDDO
1183	ENDDO
1184
1185	! Linear distrubution on sub-surface of long- and shortwave net radiance
1186	DO nsrf = 1, nbsrf
1187	DO i = 1, klon
1188	sollw(i,nsrf) = sollw_m(i) + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ts(i,nsrf))
1189
1190	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1191	! ! Martin
1192	! Apparently introduced for sisvat but not used
1193	! sollwd(i,nsrf)= sollwd_m(i)
1194	! ! Martin
1195	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1196
1197	solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
1198	ENDDO
1199	ENDDO
1200	!
1201	!<al1: second order corrections
1202	!- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)
1203	IF (iflag_order2_sollw == 1) THEN
1204	meansqT(:) = 0. ! as working buffer
1205	DO nsrf = 1, nbsrf
1206	DO i = 1, klon
1207	meansqT(i) = meansqT(i)+(ts(i,nsrf)-ztsol(i))*2 pctsrf(i,nsrf)
1208	ENDDO
1209	ENDDO
1210	DO nsrf = 1, nbsrf
1211	DO i = 1, klon
1212	sollw(i,nsrf) = sollw(i,nsrf) &
1213	+ 6.0RSIGMAztsol(i)*2 (meansqT(i)-(ztsol(i)-ts(i,nsrf))**2)
1214	ENDDO
1215	ENDDO
1216	ENDIF ! iflag_order2_sollw == 1
1217	!>al1
1218
1219	! >> PC
1220	IF (carbon_cycle_cpl .AND. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN
1221	r_co2_ppm(:) = co2_send(:)
1222	DO n=1,nbcf_out
1223	IF (cfname_out(n) == "atmco2" ) fields_out(:,n) = co2_send(:)
1224	ENDDO
1225	ENDIF
1226	IF ( .NOT. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN
1227	r_co2_ppm(:) = co2_ppm ! Constant field
1228	DO n=1,nbcf_out
1229	IF (cfname_out(n) == "atmco2" ) fields_out(:,n) = co2_ppm
1230	ENDDO
1231	ENDIF
1232	! << PC
1233
1234	!****************************************************************************************
1235	! 4) Loop over different surfaces
1236	!
1237	! Only points containing a fraction of the sub surface will be treated.
1238	!
1239	!****************************************************************************************
1240
1241	loop_nbsrf: DO nsrf = 1, nbsrf
1242	IF (prt_level >=10) print *,' Loop nsrf ',nsrf
1243
1244	! Search for index(ni) and size(knon) of domaine to treat
1245	ni(:) = 0
1246	knon = 0
1247	DO i = 1, klon
1248	IF (pctsrf(i,nsrf) > 0.) THEN
1249	knon = knon + 1
1250	ni(knon) = i
1251	ENDIF
1252	ENDDO
1253
1254	!!! jyg le 19/08/2012
1255	! IF (knon <= 0) THEN
1256	! IF (prt_level >= 10) print *,' no grid point for nsrf= ',nsrf
1257	! cycle loop_nbsrf
1258	! ENDIF
1259	!!!
1260
1261	! write index, with IOIPSL
1262	IF (debugindex .AND. mpi_size==1) THEN
1263	tabindx(:)=0.
1264	DO i=1,knon
1265	tabindx(i)=REAL(i)
1266	ENDDO
1267	debugtab(:,:) = 0.
1268	ndexbg(:) = 0
1269	CALL gath2cpl(tabindx,debugtab,knon,ni)
1270	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,nbp_lon*nbp_lat, ndexbg)
1271	ENDIF
1272
1273	!****************************************************************************************
1274	! 5) Compress variables
1275	!
1276	!****************************************************************************************
1277
1278	DO j = 1, knon
1279	i = ni(j)
1280	ypct(j) = pctsrf(i,nsrf)
1281	yts(j) = ts(i,nsrf)
1282	ysnow(j) = snow(i,nsrf)
1283	yqsurf(j) = qsurf(i,nsrf)
1284	yalb(j) = alb(i,nsrf)
1285	!albedo SB >>>
1286	yalb_vis(j) = alb_dir(i,1,nsrf)
1287	IF (nsw==6) THEN
1288	yalb_vis(j)=(alb_dir(i,1,nsrf)SFRWL(1)+alb_dir(i,2,nsrf)SFRWL(2) &
1289	+alb_dir(i,3,nsrf)*SFRWL(3))/(SFRWL(1)+SFRWL(2)+SFRWL(3))
1290	ENDIF
1291	!albedo SB <<<
1292	yrain_f(j) = rain_f(i)
1293	ysnow_f(j) = snow_f(i)
1294	yagesno(j) = agesno(i,nsrf)
1295	yfder(j) = fder(i)
1296	ylwdown(j) = lwdown_m(i)
1297	ygustiness(j) = gustiness(i)
1298	ysolsw(j) = solsw(i,nsrf)
1299	ysollw(j) = sollw(i,nsrf)
1300	yz0m(j) = z0m(i,nsrf)
1301	yz0h(j) = z0h(i,nsrf)
1302	yrugoro(j) = rugoro(i)
1303	yu1(j) = u(i,1)
1304	yv1(j) = v(i,1)
1305	ypaprs(j,klev+1) = paprs(i,klev+1)
1306	!jyg<
1307	!! ywindsp(j) = SQRT(u10m(i,nsrf)2 + v10m(i,nsrf)2 )
1308	ywindsp(j) = windsp(i,nsrf)
1309	!>jyg
1310	! Martin
1311	yzsig(j) = zsig(i)
1312	ycldt(j) = cldt(i)
1313	yrmu0(j) = rmu0(i)
1314	! Martin
1315	!!! nrlmd le 13/06/2011
1316	y_delta_tsurf(j)=delta_tsurf(i,nsrf)
1317	!!!
1318	ENDDO
1319	! >> PC
1320	!--compressing fields_out onto ORCHIDEE grid
1321	!--these fields are shared and used directly surf_land_orchidee_mod
1322	DO n = 1, nbcf_out
1323	DO j = 1, knon
1324	i = ni(j)
1325	yfields_out(j,n) = fields_out(i,n)
1326	ENDDO
1327	ENDDO
1328	! << PC
1329	DO k = 1, klev
1330	DO j = 1, knon
1331	i = ni(j)
1332	ypaprs(j,k) = paprs(i,k)
1333	ypplay(j,k) = pplay(i,k)
1334	ydelp(j,k) = delp(i,k)
1335	ENDDO
1336	ENDDO
1337	!
1338	!!! jyg le 07/02/2012 et le 10/04/2013
1339	DO k = 1, klev+1
1340	DO j = 1, knon
1341	i = ni(j)
1342	!jyg<
1343	!! ytke(j,k) = tke(i,k,nsrf)
1344	ytke(j,k) = tke_x(i,k,nsrf)
1345	ENDDO
1346	ENDDO
1347	!>jyg
1348	DO k = 1, klev
1349	DO j = 1, knon
1350	i = ni(j)
1351	y_treedrg(j,k) = treedrg(i,k,nsrf)
1352	yu(j,k) = u(i,k)
1353	yv(j,k) = v(i,k)
1354	yt(j,k) = t(i,k)
1355	yq(j,k) = q(i,k)
1356	ENDDO
1357	ENDDO
1358	!
1359	IF (iflag_split.GE.1) THEN
1360	!!! nrlmd le 02/05/2011
1361	DO k = 1, klev
1362	DO j = 1, knon
1363	i = ni(j)
1364	yu_x(j,k) = u(i,k)
1365	yv_x(j,k) = v(i,k)
1366	yt_x(j,k) = t(i,k)-wake_s(i)*wake_dlt(i,k)
1367	yq_x(j,k) = q(i,k)-wake_s(i)*wake_dlq(i,k)
1368	yu_w(j,k) = u(i,k)
1369	yv_w(j,k) = v(i,k)
1370	yt_w(j,k) = t(i,k)+(1.-wake_s(i))*wake_dlt(i,k)
1371	yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k)
1372	!!!
1373	ENDDO
1374	ENDDO
1375
1376	IF (prt_level .ge. 10) THEN
1377	print *,'pbl_surface, wake_s(1), wake_dlt(1,:) ', wake_s(1), wake_dlt(1,:)
1378	print *,'pbl_surface, wake_s(1), wake_dlq(1,:) ', wake_s(1), wake_dlq(1,:)
1379	ENDIF
1380
1381	!!! nrlmd le 02/05/2011
1382	DO k = 1, klev+1
1383	DO j = 1, knon
1384	i = ni(j)
1385	!jyg<
1386	!! ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)*wake_dltke(i,k,nsrf)
1387	!! ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))*wake_dltke(i,k,nsrf)
1388	!! ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
1389	!! ytke(j,k) = tke(i,k,nsrf)
1390	!
1391	ytke_x(j,k) = tke_x(i,k,nsrf)
1392	ytke(j,k) = tke_x(i,k,nsrf)+wake_s(i)*wake_dltke(i,k,nsrf)
1393	ytke_w(j,k) = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf)
1394	ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
1395	!>jyg
1396	ENDDO
1397	ENDDO
1398	!!!
1399	!!! jyg le 07/02/2012
1400	DO j = 1, knon
1401	i = ni(j)
1402	ywake_s(j)=wake_s(i)
1403	ywake_cstar(j)=wake_cstar(i)
1404	ywake_dens(j)=wake_dens(i)
1405	ENDDO
1406	!!!
1407	!!! nrlmd le 13/06/2011
1408	DO j=1,knon
1409	yts_x(j)=yts(j)-ywake_s(j)*y_delta_tsurf(j)
1410	yts_w(j)=yts(j)+(1.-ywake_s(j))*y_delta_tsurf(j)
1411	ENDDO
1412	!!!
1413	ENDIF ! (iflag_split .ge.1)
1414	!!!
1415	DO k = 1, nsoilmx
1416	DO j = 1, knon
1417	i = ni(j)
1418	ytsoil(j,k) = ftsoil(i,k,nsrf)
1419	ENDDO
1420	ENDDO
1421
1422	! qsol(water height in soil) only for bucket continental model
1423	IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
1424	DO j = 1, knon
1425	i = ni(j)
1426	yqsol(j) = qsol(i)
1427	ENDDO
1428	ENDIF
1429
1430	if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
1431	if (activate_ocean_skin == 2 .and. type_ocean == "couple") &
1432	ys_int(:knon) = s_int(ni(:knon))
1433	yds_ns(:knon) = ds_ns(ni(:knon))
1434	ydt_ns(:knon) = dt_ns(ni(:knon))
1435	end if
1436
1437	!****************************************************************************************
1438	! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
1439	!
1440	!****************************************************************************************
1441
1442	!!! jyg le 07/02/2012
1443	IF (iflag_split .eq.0) THEN
1444	!!!
1445	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1446	! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1447	! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1448	! yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
1449	! yts, yqsurf, yrugos, &
1450	! ycdragm, ycdragh )
1451	! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1452	DO i = 1, knon
1453	! print*,'PBL ',i,RD
1454	! print*,'PBL ',yt(i,1),ypaprs(i,1),ypplay(i,1)
1455	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1456	* (ypaprs(i,1)-ypplay(i,1))
1457	speed(i) = SQRT(yu(i,1)2+yv(i,1)2)
1458	ENDDO
1459	CALL cdrag(knon, nsrf, &
1460	speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
1461	yts, yqsurf, yz0m, yz0h, &
1462	ycdragm, ycdragh, zri1, pref )
1463
1464	! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
1465	IF (ok_prescr_ust) THEN
1466	DO i = 1, knon
1467	print *,'ycdragm avant=',ycdragm(i)
1468	vent= sqrt(yu(i,1)yu(i,1)+yv(i,1)yv(i,1))
1469	! ycdragm(i) = ustust/(1.+(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1)))
1470	! ycdragm(i) = ustust/((1.+sqrt(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1))) &
1471	! sqrt(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1)))
1472	ycdragm(i) = ust*ust/(1.+vent)/vent
1473	! print *,'ycdragm ust yu yv apres=',ycdragm(i),ust,yu(i,1),yv(i,1)
1474	ENDDO
1475	ENDIF
1476	IF (prt_level >=10) print *,'clcdrag -> ycdragh ', ycdragh
1477	ELSE !(iflag_split .eq.0)
1478
1479	! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1480	! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1481	! yu_x(:,1), yv_x(:,1), yt_x(:,1), yq_x(:,1), &
1482	! yts_x, yqsurf, yrugos, &
1483	! ycdragm_x, ycdragh_x )
1484	! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1485	DO i = 1, knon
1486	zgeo1_x(i) = RD * yt_x(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1487	* (ypaprs(i,1)-ypplay(i,1))
1488	speed_x(i) = SQRT(yu_x(i,1)2+yv_x(i,1)2)
1489	ENDDO
1490	CALL cdrag(knon, nsrf, &
1491	speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
1492	yts_x, yqsurf, yz0m, yz0h, &
1493	ycdragm_x, ycdragh_x, zri1_x, pref_x )
1494
1495	! --- special Dice. JYG+MPL 25112013
1496	IF (ok_prescr_ust) THEN
1497	DO i = 1, knon
1498	! print *,'ycdragm_x avant=',ycdragm_x(i)
1499	vent= sqrt(yu_x(i,1)yu_x(i,1)+yv_x(i,1)yv_x(i,1))
1500	ycdragm_x(i) = ust*ust/(1.+vent)/vent
1501	! print *,'ycdragm_x ust yu yv apres=',ycdragm_x(i),ust,yu_x(i,1),yv_x(i,1)
1502	ENDDO
1503	ENDIF
1504	IF (prt_level >=10) print *,'clcdrag -> ycdragh_x ', ycdragh_x
1505	!
1506	! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1507	! CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1508	! yu_w(:,1), yv_w(:,1), yt_w(:,1), yq_w(:,1), &
1509	! yts_w, yqsurf, yz0m, &
1510	! ycdragm_w, ycdragh_w )
1511	! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1512	DO i = 1, knon
1513	zgeo1_w(i) = RD * yt_w(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1514	* (ypaprs(i,1)-ypplay(i,1))
1515	speed_w(i) = SQRT(yu_w(i,1)2+yv_w(i,1)2)
1516	ENDDO
1517	CALL cdrag(knon, nsrf, &
1518	speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
1519	yts_w, yqsurf, yz0m, yz0h, &
1520	ycdragm_w, ycdragh_w, zri1_w, pref_w )
1521	!
1522	zgeo1(:) = wake_s(:)zgeo1_w(:) + (1.-wake_s(:))zgeo1_x(:)
1523
1524	! --- special Dice. JYG+MPL 25112013 puis BOMEX
1525	IF (ok_prescr_ust) THEN
1526	DO i = 1, knon
1527	! print *,'ycdragm_w avant=',ycdragm_w(i)
1528	vent= sqrt(yu_w(i,1)yu_w(i,1)+yv_w(i,1)yv_w(i,1))
1529	ycdragm_w(i) = ust*ust/(1.+vent)/vent
1530	! print *,'ycdragm_w ust yu yv apres=',ycdragm_w(i),ust,yu_w(i,1),yv_w(i,1)
1531	ENDDO
1532	ENDIF
1533	IF (prt_level >=10) print *,'clcdrag -> ycdragh_w ', ycdragh_w
1534	!!!
1535	ENDIF ! (iflag_split .eq.0)
1536	!!!
1537
1538
1539	!****************************************************************************************
1540	! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm.
1541	!
1542	!****************************************************************************************
1543
1544	!!! jyg le 07/02/2012
1545	IF (iflag_split .eq.0) THEN
1546	!!!
1547	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1548	IF (prt_level >=10) THEN
1549	print *,' args coef_diff_turb: yu ', yu
1550	print *,' args coef_diff_turb: yv ', yv
1551	print *,' args coef_diff_turb: yq ', yq
1552	print *,' args coef_diff_turb: yt ', yt
1553	print *,' args coef_diff_turb: yts ', yts
1554	print *,' args coef_diff_turb: yz0m ', yz0m
1555	print *,' args coef_diff_turb: yqsurf ', yqsurf
1556	print *,' args coef_diff_turb: ycdragm ', ycdragm
1557	print *,' args coef_diff_turb: ycdragh ', ycdragh
1558	print *,' args coef_diff_turb: ytke ', ytke
1559	ENDIF
1560	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
1561	ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
1562	ycoefm, ycoefh, ytke, y_treedrg)
1563	! ycoefm, ycoefh, ytke)
1564	!FC y_treedrg ajoute
1565	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1566	! In this case, coef_diff_turb is called for the Cd only
1567	DO k = 2, klev
1568	DO j = 1, knon
1569	i = ni(j)
1570	ycoefh(j,k) = zcoefh(i,k,nsrf)
1571	ycoefm(j,k) = zcoefm(i,k,nsrf)
1572	ENDDO
1573	ENDDO
1574	ENDIF
1575	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh
1576	!
1577	ELSE !(iflag_split .eq.0)
1578	IF (prt_level >=10) THEN
1579	print *,' args coef_diff_turb: yu_x ', yu_x
1580	print *,' args coef_diff_turb: yv_x ', yv_x
1581	print *,' args coef_diff_turb: yq_x ', yq_x
1582	print *,' args coef_diff_turb: yt_x ', yt_x
1583	print *,' args coef_diff_turb: yts_x ', yts_x
1584	print *,' args coef_diff_turb: yqsurf ', yqsurf
1585	print *,' args coef_diff_turb: ycdragm_x ', ycdragm_x
1586	print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x
1587	print *,' args coef_diff_turb: ytke_x ', ytke_x
1588	ENDIF
1589	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
1590	ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf, ycdragm_x, &
1591	ycoefm_x, ycoefh_x, ytke_x,y_treedrg)
1592	! ycoefm_x, ycoefh_x, ytke_x)
1593	!FC doit on le mettre ( on ne l utilise pas si il y a du spliting)
1594	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1595	! In this case, coef_diff_turb is called for the Cd only
1596	DO k = 2, klev
1597	DO j = 1, knon
1598	i = ni(j)
1599	ycoefh_x(j,k) = zcoefh(i,k,nsrf)
1600	ycoefm_x(j,k) = zcoefm(i,k,nsrf)
1601	ENDDO
1602	ENDDO
1603	ENDIF
1604	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x
1605	!
1606	IF (prt_level >=10) THEN
1607	print *,' args coef_diff_turb: yu_w ', yu_w
1608	print *,' args coef_diff_turb: yv_w ', yv_w
1609	print *,' args coef_diff_turb: yq_w ', yq_w
1610	print *,' args coef_diff_turb: yt_w ', yt_w
1611	print *,' args coef_diff_turb: yts_w ', yts_w
1612	print *,' args coef_diff_turb: yqsurf ', yqsurf
1613	print *,' args coef_diff_turb: ycdragm_w ', ycdragm_w
1614	print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w
1615	print *,' args coef_diff_turb: ytke_w ', ytke_w
1616	ENDIF
1617	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
1618	ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf, ycdragm_w, &
1619	ycoefm_w, ycoefh_w, ytke_w,y_treedrg)
1620	! ycoefm_w, ycoefh_w, ytke_w)
1621	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1622	! In this case, coef_diff_turb is called for the Cd only
1623	DO k = 2, klev
1624	DO j = 1, knon
1625	i = ni(j)
1626	ycoefh_w(j,k) = zcoefh(i,k,nsrf)
1627	ycoefm_w(j,k) = zcoefm(i,k,nsrf)
1628	ENDDO
1629	ENDDO
1630	ENDIF
1631	IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w
1632	!
1633	!!!jyg le 10/04/2013
1634	!! En attendant de traiter le transport des traceurs dans les poches froides, formule
1635	!! arbitraire pour ycoefh et ycoefm
1636	DO k = 2,klev
1637	DO j = 1,knon
1638	ycoefh(j,k) = ycoefh_x(j,k) + ywake_s(j)*(ycoefh_w(j,k) - ycoefh_x(j,k))
1639	ycoefm(j,k) = ycoefm_x(j,k) + ywake_s(j)*(ycoefm_w(j,k) - ycoefm_x(j,k))
1640	ENDDO
1641	ENDDO
1642	!!!
1643	ENDIF ! (iflag_split .eq.0)
1644	!!!
1645
1646	!****************************************************************************************
1647	!
1648	! 8) "La descente" - "The downhill"
1649	!
1650	! climb_hq_down and climb_wind_down calculate the coefficients
1651	! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
1652	! Only the coefficients at surface for H and Q are returned.
1653	!
1654	!****************************************************************************************
1655
1656	! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
1657	!!! jyg le 07/02/2012
1658	IF (iflag_split .eq.0) THEN
1659	!!!
1660	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1661	CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
1662	ydelp, yt, yq, dtime, &
1663	!!! jyg le 09/05/2011
1664	CcoefH, CcoefQ, DcoefH, DcoefQ, &
1665	Kcoef_hq, gama_q, gama_h, &
1666	!!!
1667	AcoefH, AcoefQ, BcoefH, BcoefQ)
1668	ELSE !(iflag_split .eq.0)
1669	CALL climb_hq_down(knon, ycoefh_x, ypaprs, ypplay, &
1670	ydelp, yt_x, yq_x, dtime, &
1671	!!! nrlmd le 02/05/2011
1672	CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
1673	Kcoef_hq_x, gama_q_x, gama_h_x, &
1674	!!!
1675	AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x)
1676	!!!
1677	IF (prt_level >=10) THEN
1678	PRINT *,'pbl_surface (climb_hq_down.x->) AcoefH_x ',AcoefH_x
1679	PRINT *,'pbl_surface (climb_hq_down.x->) AcoefQ_x ',AcoefQ_x
1680	PRINT *,'pbl_surface (climb_hq_down.x->) BcoefH_x ',BcoefH_x
1681	PRINT *,'pbl_surface (climb_hq_down.x->) BcoefQ_x ',BcoefQ_x
1682	ENDIF
1683	!
1684	CALL climb_hq_down(knon, ycoefh_w, ypaprs, ypplay, &
1685	ydelp, yt_w, yq_w, dtime, &
1686	!!! nrlmd le 02/05/2011
1687	CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
1688	Kcoef_hq_w, gama_q_w, gama_h_w, &
1689	!!!
1690	AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w)
1691	!!!
1692	IF (prt_level >=10) THEN
1693	PRINT *,'pbl_surface (climb_hq_down.w->) AcoefH_w ',AcoefH_w
1694	PRINT *,'pbl_surface (climb_hq_down.w->) AcoefQ_w ',AcoefQ_w
1695	PRINT *,'pbl_surface (climb_hq_down.w->) BcoefH_w ',BcoefH_w
1696	PRINT *,'pbl_surface (climb_hq_down.w->) BcoefQ_w ',BcoefQ_w
1697	ENDIF
1698	!!!
1699	ENDIF ! (iflag_split .eq.0)
1700	!!!
1701
1702	! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
1703	!!! jyg le 07/02/2012
1704	IF (iflag_split .eq.0) THEN
1705	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1706	CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
1707	!!! jyg le 09/05/2011
1708	CcoefU, CcoefV, DcoefU, DcoefV, &
1709	Kcoef_m, alf_1, alf_2, &
1710	!!!
1711	AcoefU, AcoefV, BcoefU, BcoefV)
1712	ELSE ! (iflag_split .eq.0)
1713	CALL climb_wind_down(knon, dtime, ycoefm_x, ypplay, ypaprs, yt_x, ydelp, yu_x, yv_x, &
1714	!!! nrlmd le 02/05/2011
1715	CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
1716	Kcoef_m_x, alf_1_x, alf_2_x, &
1717	!!!
1718	AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x)
1719	!
1720	CALL climb_wind_down(knon, dtime, ycoefm_w, ypplay, ypaprs, yt_w, ydelp, yu_w, yv_w, &
1721	!!! nrlmd le 02/05/2011
1722	CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
1723	Kcoef_m_w, alf_1_w, alf_2_w, &
1724	!!!
1725	AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w)
1726	!!!
1727	ENDIF ! (iflag_split .eq.0)
1728	!!!
1729
1730	!****************************************************************************************
1731	! 9) Small calculations
1732	!
1733	!****************************************************************************************
1734
1735	! - Reference pressure is given the values at surface level
1736	ypsref(:) = ypaprs(:,1)
1737
1738	! - CO2 field on 2D grid to be sent to ORCHIDEE
1739	! Transform to compressed field
1740	IF (carbon_cycle_cpl) THEN
1741	DO i=1,knon
1742	r_co2_ppm(i) = co2_send(ni(i))
1743	ENDDO
1744	ELSE
1745	r_co2_ppm(:) = co2_ppm ! Constant field
1746	ENDIF
1747
1748	!!! nrlmd le 02/05/2011 -----------------------On raccorde les 2 colonnes dans la couche 1
1749	!----------------------------------------------------------------------------------------
1750	!!! jyg le 07/02/2012
1751	!!! jyg le 01/02/2017
1752	IF (iflag_split .eq. 0) THEN
1753	yt1(:) = yt(:,1)
1754	yq1(:) = yq(:,1)
1755	!! ELSE IF (iflag_split .eq. 1) THEN
1756	!!!
1757	!jyg<
1758	!! CALL wx_pbl_fuse_no_dts(knon, dtime, ypplay, ywake_s, &
1759	!! yt_x, yt_w, yq_x, yq_w, &
1760	!! yu_x, yu_w, yv_x, yv_w, &
1761	!! ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
1762	!! AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
1763	!! AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
1764	!! BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
1765	!! BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
1766	!! AcoefH, AcoefQ, AcoefU, AcoefV, &
1767	!! BcoefH, BcoefQ, BcoefU, BcoefV, &
1768	!! ycdragh, ycdragm, &
1769	!! yt1, yq1, yu1, yv1 &
1770	!! )
1771	ELSE IF (iflag_split .ge. 1) THEN
1772	CALL wx_pbl0_fuse(knon, dtime, ypplay, ywake_s, &
1773	yt_x, yt_w, yq_x, yq_w, &
1774	yu_x, yu_w, yv_x, yv_w, &
1775	ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
1776	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
1777	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
1778	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
1779	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
1780	AcoefH, AcoefQ, AcoefU, AcoefV, &
1781	BcoefH, BcoefQ, BcoefU, BcoefV, &
1782	ycdragh, ycdragm, &
1783	yt1, yq1, yu1, yv1 &
1784	)
1785	!! ELSE IF (iflag_split .ge.2) THEN
1786	!!! Provisoire
1787	!! ah(:) = 0.
1788	!! bh(:) = 0.
1789	!! IF (nsrf == is_oce) THEN
1790	!! ybeta(:) = 1.
1791	!! ELSE
1792	!! ybeta(:) = beta_land
1793	!! ENDIF
1794	!! ycdragh(:) = ywake_s(:)ycdragh_w(:) + (1.-ywake_s(:))ycdragh_x(:)
1795	!! CALL wx_dts(knon, nsrf, ywake_cstar, ywake_s, ywake_dens, &
1796	!! yts, ypplay(:,1), ybeta, ycdragh , ypaprs(:,1), &
1797	!! yq(:,1), yt(:,1), yu(:,1), yv(:,1), ygustiness, &
1798	!! ah, bh &
1799	!! )
1800	!!!
1801	!! CALL wx_pbl_fuse(knon, dtime, ypplay, ywake_s, &
1802	!! yt_x, yt_w, yq_x, yq_w, &
1803	!! yu_x, yu_w, yv_x, yv_w, &
1804	!! ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
1805	!! AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
1806	!! AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
1807	!! BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
1808	!! BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
1809	!! ah, bh, &
1810	!! AcoefH, AcoefQ, AcoefU, AcoefV, &
1811	!! BcoefH, BcoefQ, BcoefU, BcoefV, &
1812	!! ycdragh, ycdragm, &
1813	!! yt1, yq1, yu1, yv1 &
1814	!! )
1815	!>jyg
1816	!!!
1817	ENDIF ! (iflag_split .eq.0)
1818	!!!
1819	IF (prt_level >=10) THEN
1820	PRINT *,'pbl_surface (fuse->): yt(1,:) ',yt(1,:)
1821	PRINT *,'pbl_surface (fuse->): yq(1,:) ',yq(1,:)
1822	PRINT *,'pbl_surface (fuse->): yu(1,:) ',yu(1,:)
1823	PRINT *,'pbl_surface (fuse->): yv(1,:) ',yv(1,:)
1824	PRINT *,'pbl_surface (fuse->): AcoefH(1) ',AcoefH(1)
1825	PRINT *,'pbl_surface (fuse->): BcoefH(1) ',BcoefH(1)
1826	ENDIF
1827
1828	!****************************************************************************************
1829	!
1830	! Calulate t2m and q2m for the case of calculation at land grid points
1831	! t2m and q2m are needed as input to ORCHIDEE
1832	!
1833	!****************************************************************************************
1834	IF (nsrf == is_ter) THEN
1835
1836	DO i = 1, knon
1837	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1838	* (ypaprs(i,1)-ypplay(i,1))
1839	ENDDO
1840
1841	! Calculate the temperature et relative humidity at 2m and the wind at 10m
1842	CALL stdlevvar(klon, knon, is_ter, zxli, &
1843	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
1844	yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
1845	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
1846
1847	ENDIF
1848
1849	!****************************************************************************************
1850	!
1851	! 10) Switch according to current surface
1852	! It is necessary to start with the continental surfaces because the ocean
1853	! needs their run-off.
1854	!
1855	!****************************************************************************************
1856	SELECT CASE(nsrf)
1857
1858	CASE(is_ter)
1859	! print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)
1860	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
1861	rlon, rlat, yrmu0, &
1862	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
1863	!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1864	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
1865	AcoefH, AcoefQ, BcoefH, BcoefQ, &
1866	AcoefU, AcoefV, BcoefU, BcoefV, &
1867	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1868	ylwdown, yq2m, yt2m, &
1869	ysnow, yqsol, yagesno, ytsoil, &
1870	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
1871	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
1872	y_flux_u1, y_flux_v1, &
1873	yveget,ylai,yheight )
1874
1875	!FC quid qd yveget ylai yheight ne sont pas definit
1876	!FC yveget,ylai,yheight, &
1877	IF (ifl_pbltree .ge. 1) THEN
1878	CALL freinage(knon, yu, yv, yt, &
1879	! yveget,ylai, yheight,ypaprs,ypplay,y_d_u_frein,y_d_v_frein)
1880	yveget,ylai, yheight,ypaprs,ypplay,y_treedrg, y_d_u_frein,y_d_v_frein)
1881	ENDIF
1882
1883
1884	! Special DICE MPL 05082013 puis BOMEX
1885	IF (ok_prescr_ust) THEN
1886	DO j=1,knon
1887	! ysnow(:)=0.
1888	! yqsol(:)=0.
1889	! yagesno(:)=50.
1890	! ytsoil(:,:)=300.
1891	! yz0_new(:)=0.001
1892	! yevap(:)=flat/RLVTT
1893	! yfluxlat(:)=-flat
1894	! yfluxsens(:)=-fsens
1895	! yqsurf(:)=0.
1896	! ytsurf_new(:)=tg
1897	! y_dflux_t(:)=0.
1898	! y_dflux_q(:)=0.
1899	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)
1900	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)
1901	ENDDO
1902	ENDIF
1903
1904	CASE(is_lic)
1905	! Martin
1906	CALL surf_landice(itap, dtime, knon, ni, &
1907	rlon, rlat, debut, lafin, &
1908	yrmu0, ylwdown, yalb, ypphi(:,1), &
1909	ysolsw, ysollw, yts, ypplay(:,1), &
1910	!!jyg ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1911	ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
1912	AcoefH, AcoefQ, BcoefH, BcoefQ, &
1913	AcoefU, AcoefV, BcoefU, BcoefV, &
1914	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1915	ysnow, yqsurf, yqsol, yagesno, &
1916	ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
1917	ytsurf_new, y_dflux_t, y_dflux_q, &
1918	yzsig, ycldt, &
1919	ysnowhgt, yqsnow, ytoice, ysissnow, &
1920	yalb3_new, yrunoff, &
1921	y_flux_u1, y_flux_v1)
1922
1923	!jyg<
1924	!! alb3_lic(:)=0.
1925	!>jyg
1926	DO j = 1, knon
1927	i = ni(j)
1928	alb3_lic(i) = yalb3_new(j)
1929	snowhgt(i) = ysnowhgt(j)
1930	qsnow(i) = yqsnow(j)
1931	to_ice(i) = ytoice(j)
1932	sissnow(i) = ysissnow(j)
1933	runoff(i) = yrunoff(j)
1934	ENDDO
1935	! Martin
1936	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
1937	IF (ok_prescr_ust) THEN
1938	DO j=1,knon
1939	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)
1940	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)
1941	ENDDO
1942	ENDIF
1943
1944	CASE(is_oce)
1945	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
1946	ywindsp, rmu0, yfder, yts, &
1947	itap, dtime, jour, knon, ni, &
1948	!!jyg ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1949	ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& ! ym missing init
1950	AcoefH, AcoefQ, BcoefH, BcoefQ, &
1951	AcoefU, AcoefV, BcoefU, BcoefV, &
1952	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1953	ysnow, yqsurf, yagesno, &
1954	yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
1955	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
1956	y_flux_u1, y_flux_v1, ysst_nff(:knon), ys_int(:knon), &
1957	yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
1958	ytkt(:knon), ytks(:knon), ytaur(:knon), ysss)
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	if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
2563	s_int = missing_val
2564	ds_ns = missing_val
2565	dt_ns = missing_val
2566	sst_nff = missing_val
2567	dter = missing_val
2568	dser = missing_val
2569	tkt = missing_val
2570	tks = missing_val
2571	taur = missing_val
2572	sss = missing_val
2573
2574	s_int(ni(:knon)) = ys_int(:knon)
2575	ds_ns(ni(:knon)) = yds_ns(:knon)
2576	dt_ns(ni(:knon)) = ydt_ns(:knon)
2577	sst_nff(ni(:knon)) = ysst_nff(:knon)
2578	dter(ni(:knon)) = ydter(:knon)
2579	dser(ni(:knon)) = ydser(:knon)
2580	tkt(ni(:knon)) = ytkt(:knon)
2581	tks(ni(:knon)) = ytks(:knon)
2582	taur(ni(:knon)) = ytaur(:knon)
2583	sss(ni(:knon)) = ysss(:knon)
2584	end if
2585
2586	!****************************************************************************************
2587	! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m
2588	! Call HBTM
2589	!
2590	!****************************************************************************************
2591	!!!
2592	!
2593	#undef T2m
2594	#define T2m
2595	#ifdef T2m
2596	! Calculations of diagnostic t,q at 2m and u, v at 10m
2597
2598	! print*,'Dans pbl OK41'
2599	! print*,'tair1,yt(:,1),y_d_t(:,1)'
2600	! print*, tair1,yt(:,1),y_d_t(:,1)
2601	!!! jyg le 07/02/2012
2602	IF (iflag_split .eq.0) THEN
2603	DO j=1, knon
2604	uzon(j) = yu(j,1) + y_d_u(j,1)
2605	vmer(j) = yv(j,1) + y_d_v(j,1)
2606	tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
2607	qair1(j) = yq(j,1) + y_d_q(j,1)
2608	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2609	* (ypaprs(j,1)-ypplay(j,1))
2610	tairsol(j) = yts(j) + y_d_ts(j)
2611	qairsol(j) = yqsurf(j)
2612	ENDDO
2613	ELSE ! (iflag_split .eq.0)
2614	DO j=1, knon
2615	uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
2616	vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
2617	tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
2618	qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
2619	zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2620	* (ypaprs(j,1)-ypplay(j,1))
2621	tairsol(j) = yts(j) + y_d_ts(j)
2622	tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
2623	qairsol(j) = yqsurf(j)
2624	ENDDO
2625	DO j=1, knon
2626	uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
2627	vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
2628	tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
2629	qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
2630	zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2631	* (ypaprs(j,1)-ypplay(j,1))
2632	tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
2633	qairsol(j) = yqsurf(j)
2634	ENDDO
2635	!!!
2636	ENDIF ! (iflag_split .eq.0)
2637	!!!
2638	DO j=1, knon
2639	! i = ni(j)
2640	! yz0h_oupas(j) = yz0m(j)
2641	! IF(nsrf.EQ.is_oce) THEN
2642	! yz0h_oupas(j) = z0m(i,nsrf)
2643	! ENDIF
2644	psfce(j)=ypaprs(j,1)
2645	patm(j)=ypplay(j,1)
2646	ENDDO
2647
2648	IF (iflag_pbl_surface_t2m_bug==1) THEN
2649	yz0h_oupas(1:knon)=yz0m(1:knon)
2650	ELSE
2651	yz0h_oupas(1:knon)=yz0h(1:knon)
2652	ENDIF
2653
2654	! print*,'Dans pbl OK42A'
2655	! print*,'tair1,yt(:,1),y_d_t(:,1)'
2656	! print*, tair1,yt(:,1),y_d_t(:,1)
2657
2658	! Calculate the temperature and relative humidity at 2m and the wind at 10m
2659	!!! jyg le 07/02/2012
2660	IF (iflag_split .eq.0) THEN
2661	CALL stdlevvar(klon, knon, nsrf, zxli, &
2662	uzon, vmer, tair1, qair1, zgeo1, &
2663	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
2664	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
2665	ELSE !(iflag_split .eq.0)
2666	CALL stdlevvar(klon, knon, nsrf, zxli, &
2667	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
2668	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
2669	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
2670	CALL stdlevvar(klon, knon, nsrf, zxli, &
2671	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
2672	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
2673	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
2674	!!!
2675	ENDIF ! (iflag_split .eq.0)
2676	!!!
2677	!!! jyg le 07/02/2012
2678	IF (iflag_split .eq.0) THEN
2679	DO j=1, knon
2680	i = ni(j)
2681	t2m(i,nsrf)=yt2m(j)
2682	q2m(i,nsrf)=yq2m(j)
2683	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2684	ustar(i,nsrf)=yustar(j)
2685	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
2686	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
2687	ENDDO
2688	ELSE !(iflag_split .eq.0)
2689	DO j=1, knon
2690	i = ni(j)
2691	t2m_x(i,nsrf)=yt2m_x(j)
2692	q2m_x(i,nsrf)=yq2m_x(j)
2693	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2694	ustar_x(i,nsrf)=yustar_x(j)
2695	u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
2696	v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
2697	ENDDO
2698	DO j=1, knon
2699	i = ni(j)
2700	t2m_w(i,nsrf)=yt2m_w(j)
2701	q2m_w(i,nsrf)=yq2m_w(j)
2702	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2703	ustar_w(i,nsrf)=yustar_w(j)
2704	u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
2705	v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
2706	!
2707	ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
2708	u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
2709	v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
2710	ENDDO
2711	!!!
2712	ENDIF ! (iflag_split .eq.0)
2713	!!!
2714
2715	! print*,'Dans pbl OK43'
2716	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
2717	!IM Ajoute dependance type surface
2718	IF (thermcep) THEN
2719	!!! jyg le 07/02/2012
2720	IF (iflag_split .eq.0) THEN
2721	DO j = 1, knon
2722	i=ni(j)
2723	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
2724	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
2725	zx_qs1 = MIN(0.5,zx_qs1)
2726	zcor1 = 1./(1.-RETV*zx_qs1)
2727	zx_qs1 = zx_qs1*zcor1
2728
2729	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
2730	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
2731	ENDDO
2732	ELSE ! (iflag_split .eq.0)
2733	DO j = 1, knon
2734	i=ni(j)
2735	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
2736	zx_qs1 = r2es * FOEEW(yt2m_x(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_x(i) = rh2m_x(i) + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
2742	qsat2m_x(i) = qsat2m_x(i) + zx_qs1 * pctsrf(i,nsrf)
2743	ENDDO
2744	DO j = 1, knon
2745	i=ni(j)
2746	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
2747	zx_qs1 = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
2748	zx_qs1 = MIN(0.5,zx_qs1)
2749	zcor1 = 1./(1.-RETV*zx_qs1)
2750	zx_qs1 = zx_qs1*zcor1
2751
2752	rh2m_w(i) = rh2m_w(i) + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
2753	qsat2m_w(i) = qsat2m_w(i) + zx_qs1 * pctsrf(i,nsrf)
2754	ENDDO
2755	!!!
2756	ENDIF ! (iflag_split .eq.0)
2757	!!!
2758	ENDIF
2759	!
2760	IF (prt_level >=10) THEN
2761	print *, 'T2m, q2m, RH2m ', &
2762	t2m, q2m, rh2m
2763	ENDIF
2764
2765	! print*,'OK pbl 5'
2766	!
2767	!!! jyg le 07/02/2012
2768	IF (iflag_split .eq.0) THEN
2769	CALL hbtm(knon, ypaprs, ypplay, &
2770	yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
2771	y_flux_t,y_flux_q,yu,yv,yt,yq, &
2772	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
2773	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
2774	IF (prt_level >=10) THEN
2775	print *,' Arg. de HBTM: yt2m ',yt2m
2776	print *,' Arg. de HBTM: yt10m ',yt10m
2777	print *,' Arg. de HBTM: yq2m ',yq2m
2778	print *,' Arg. de HBTM: yq10m ',yq10m
2779	print *,' Arg. de HBTM: yustar ',yustar
2780	print *,' Arg. de HBTM: y_flux_t ',y_flux_t
2781	print *,' Arg. de HBTM: y_flux_q ',y_flux_q
2782	print *,' Arg. de HBTM: yu ',yu
2783	print *,' Arg. de HBTM: yv ',yv
2784	print *,' Arg. de HBTM: yt ',yt
2785	print *,' Arg. de HBTM: yq ',yq
2786	ENDIF
2787	ELSE ! (iflag_split .eq.0)
2788	CALL HBTM(knon, ypaprs, ypplay, &
2789	yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
2790	y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
2791	ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
2792	ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
2793	IF (prt_level >=10) THEN
2794	print *,' Arg. de HBTM: yt2m_x ',yt2m_x
2795	print *,' Arg. de HBTM: yt10m_x ',yt10m_x
2796	print *,' Arg. de HBTM: yq2m_x ',yq2m_x
2797	print *,' Arg. de HBTM: yq10m_x ',yq10m_x
2798	print *,' Arg. de HBTM: yustar_x ',yustar_x
2799	print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x
2800	print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x
2801	print *,' Arg. de HBTM: yu_x ',yu_x
2802	print *,' Arg. de HBTM: yv_x ',yv_x
2803	print *,' Arg. de HBTM: yt_x ',yt_x
2804	print *,' Arg. de HBTM: yq_x ',yq_x
2805	ENDIF
2806	CALL HBTM(knon, ypaprs, ypplay, &
2807	yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
2808	y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
2809	ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
2810	ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
2811	!!!
2812	ENDIF ! (iflag_split .eq.0)
2813	!!!
2814
2815	!!! jyg le 07/02/2012
2816	IF (iflag_split .eq.0) THEN
2817	!!!
2818	DO j=1, knon
2819	i = ni(j)
2820	pblh(i,nsrf) = ypblh(j)
2821	wstar(i,nsrf) = ywstar(j)
2822	plcl(i,nsrf) = ylcl(j)
2823	capCL(i,nsrf) = ycapCL(j)
2824	oliqCL(i,nsrf) = yoliqCL(j)
2825	cteiCL(i,nsrf) = ycteiCL(j)
2826	pblT(i,nsrf) = ypblT(j)
2827	therm(i,nsrf) = ytherm(j)
2828	trmb1(i,nsrf) = ytrmb1(j)
2829	trmb2(i,nsrf) = ytrmb2(j)
2830	trmb3(i,nsrf) = ytrmb3(j)
2831	ENDDO
2832	IF (prt_level >=10) THEN
2833	print *, 'After HBTM: pblh ', pblh
2834	print *, 'After HBTM: plcl ', plcl
2835	print *, 'After HBTM: cteiCL ', cteiCL
2836	ENDIF
2837	ELSE !(iflag_split .eq.0)
2838	DO j=1, knon
2839	i = ni(j)
2840	pblh_x(i,nsrf) = ypblh_x(j)
2841	wstar_x(i,nsrf) = ywstar_x(j)
2842	plcl_x(i,nsrf) = ylcl_x(j)
2843	capCL_x(i,nsrf) = ycapCL_x(j)
2844	oliqCL_x(i,nsrf) = yoliqCL_x(j)
2845	cteiCL_x(i,nsrf) = ycteiCL_x(j)
2846	pblT_x(i,nsrf) = ypblT_x(j)
2847	therm_x(i,nsrf) = ytherm_x(j)
2848	trmb1_x(i,nsrf) = ytrmb1_x(j)
2849	trmb2_x(i,nsrf) = ytrmb2_x(j)
2850	trmb3_x(i,nsrf) = ytrmb3_x(j)
2851	ENDDO
2852	IF (prt_level >=10) THEN
2853	print *, 'After HBTM: pblh_x ', pblh_x
2854	print *, 'After HBTM: plcl_x ', plcl_x
2855	print *, 'After HBTM: cteiCL_x ', cteiCL_x
2856	ENDIF
2857	DO j=1, knon
2858	i = ni(j)
2859	pblh_w(i,nsrf) = ypblh_w(j)
2860	wstar_w(i,nsrf) = ywstar_w(j)
2861	plcl_w(i,nsrf) = ylcl_w(j)
2862	capCL_w(i,nsrf) = ycapCL_w(j)
2863	oliqCL_w(i,nsrf) = yoliqCL_w(j)
2864	cteiCL_w(i,nsrf) = ycteiCL_w(j)
2865	pblT_w(i,nsrf) = ypblT_w(j)
2866	therm_w(i,nsrf) = ytherm_w(j)
2867	trmb1_w(i,nsrf) = ytrmb1_w(j)
2868	trmb2_w(i,nsrf) = ytrmb2_w(j)
2869	trmb3_w(i,nsrf) = ytrmb3_w(j)
2870	ENDDO
2871	IF (prt_level >=10) THEN
2872	print *, 'After HBTM: pblh_w ', pblh_w
2873	print *, 'After HBTM: plcl_w ', plcl_w
2874	print *, 'After HBTM: cteiCL_w ', cteiCL_w
2875	ENDIF
2876	!!!
2877	ENDIF ! (iflag_split .eq.0)
2878	!!!
2879
2880	! print*,'OK pbl 6'
2881	#else
2882	! T2m not defined
2883	! No calculation
2884	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
2885	#endif
2886
2887	!****************************************************************************************
2888	! 15) End of loop over different surfaces
2889	!
2890	!****************************************************************************************
2891	ENDDO loop_nbsrf
2892
2893	!****************************************************************************************
2894	! 16) Calculate the mean value over all sub-surfaces for some variables
2895	!
2896	!****************************************************************************************
2897
2898	z0m(:,nbsrf+1) = 0.0
2899	z0h(:,nbsrf+1) = 0.0
2900	DO nsrf = 1, nbsrf
2901	DO i = 1, klon
2902	z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)
2903	z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)
2904	ENDDO
2905	ENDDO
2906
2907	! print*,'OK pbl 7'
2908	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
2909	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
2910	zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
2911	zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
2912	zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
2913	zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
2914
2915	!!! jyg le 07/02/2012
2916	IF (iflag_split .ge.1) THEN
2917	!!!
2918	!!! nrlmd & jyg les 02/05/2011, 05/02/2012
2919
2920	DO nsrf = 1, nbsrf
2921	DO k = 1, klev
2922	DO i = 1, klon
2923	zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
2924	zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
2925	zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
2926	zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
2927	!
2928	zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
2929	zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
2930	zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
2931	zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
2932	ENDDO
2933	ENDDO
2934	ENDDO
2935
2936	DO i = 1, klon
2937	zxsens_x(i) = - zxfluxt_x(i,1)
2938	zxsens_w(i) = - zxfluxt_w(i,1)
2939	ENDDO
2940	!!!
2941	ENDIF ! (iflag_split .ge.1)
2942	!!!
2943
2944	DO nsrf = 1, nbsrf
2945	DO k = 1, klev
2946	DO i = 1, klon
2947	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
2948	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
2949	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
2950	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
2951	ENDDO
2952	ENDDO
2953	ENDDO
2954
2955	DO i = 1, klon
2956	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
2957	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
2958	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
2959	ENDDO
2960	!!!
2961
2962	!
2963	! Incrementer la temperature du sol
2964	!
2965	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
2966	zt2m(:) = 0.0 ; zq2m(:) = 0.0
2967	zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
2968	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
2969	!!! jyg le 07/02/2012
2970	s_pblh_x(:) = 0.0 ; s_plcl_x(:) = 0.0
2971	s_pblh_w(:) = 0.0 ; s_plcl_w(:) = 0.0
2972	!!!
2973	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
2974	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
2975	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
2976	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
2977	wstar(:,is_ave)=0.
2978
2979	! print*,'OK pbl 9'
2980
2981	!!! nrlmd le 02/05/2011
2982	zxfluxlat_x(:) = 0.0 ; zxfluxlat_w(:) = 0.0
2983	!!!
2984
2985	DO nsrf = 1, nbsrf
2986	DO i = 1, klon
2987	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
2988
2989	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
2990	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
2991	wfbilo(i,nsrf) = (evap(i,nsrf)-(rain_f(i)+snow_f(i)))*pctsrf(i,nsrf)
2992	wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)
2993	wfrain(i,nsrf) = rain_f(i)*pctsrf(i,nsrf)
2994	wfsnow(i,nsrf) = snow_f(i)*pctsrf(i,nsrf)
2995
2996	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
2997	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
2998	ENDDO
2999	ENDDO
3000	!
3001	!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
3002	IF (iflag_order2_sollw == 1) THEN
3003	meansqT(:) = 0. ! as working buffer
3004	DO nsrf = 1, nbsrf
3005	DO i = 1, klon
3006	meansqT(i) = meansqT(i)+(ts(i,nsrf)-zxtsol(i))*2 pctsrf(i,nsrf)
3007	ENDDO
3008	ENDDO
3009	zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)
3010	ENDIF ! iflag_order2_sollw == 1
3011	!>al1
3012
3013	!!! jyg le 07/02/2012
3014	IF (iflag_split .eq.0) THEN
3015	DO nsrf = 1, nbsrf
3016	DO i = 1, klon
3017	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
3018	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
3019	zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
3020	wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
3021	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
3022	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
3023
3024	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
3025	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
3026	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
3027	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
3028	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
3029	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
3030	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
3031	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
3032	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
3033	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
3034	ENDDO
3035	ENDDO
3036	ELSE !(iflag_split .eq.0)
3037	DO nsrf = 1, nbsrf
3038	DO i = 1, klon
3039	!!! nrlmd le 02/05/2011
3040	zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
3041	zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
3042	!!!
3043	!!! jyg le 08/02/2012
3044	!! Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;
3045	!! pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
3046	!! pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
3047	!! pour les autres variables, on sort les valeurs de la region (x).
3048	zt2m(i) = zt2m(i) + (t2m_x(i,nsrf)+wake_s(i)(t2m_w(i,nsrf)-t2m_x(i,nsrf))) pctsrf(i,nsrf)
3049	zq2m(i) = zq2m(i) + q2m_x(i,nsrf) * pctsrf(i,nsrf)
3050	zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
3051	wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
3052	zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
3053	zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
3054	!
3055	s_pblh(i) = s_pblh(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
3056	s_pblh_x(i) = s_pblh_x(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
3057	s_pblh_w(i) = s_pblh_w(i) + pblh_w(i,nsrf) * pctsrf(i,nsrf)
3058	!
3059	s_plcl(i) = s_plcl(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
3060	s_plcl_x(i) = s_plcl_x(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
3061	s_plcl_w(i) = s_plcl_w(i) + plcl_w(i,nsrf) * pctsrf(i,nsrf)
3062	!
3063	s_capCL(i) = s_capCL(i) + capCL_x(i,nsrf) * pctsrf(i,nsrf)
3064	s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
3065	s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
3066	s_pblT(i) = s_pblT(i) + pblT_x(i,nsrf) * pctsrf(i,nsrf)
3067	s_therm(i) = s_therm(i) + therm_x(i,nsrf) * pctsrf(i,nsrf)
3068	s_trmb1(i) = s_trmb1(i) + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
3069	s_trmb2(i) = s_trmb2(i) + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
3070	s_trmb3(i) = s_trmb3(i) + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
3071	ENDDO
3072	ENDDO
3073	DO i = 1, klon
3074	qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
3075	ENDDO
3076	!!!
3077	ENDIF ! (iflag_split .eq.0)
3078	!!!
3079
3080	IF (check) THEN
3081	amn=MIN(ts(1,is_ter),1000.)
3082	amx=MAX(ts(1,is_ter),-1000.)
3083	DO i=2, klon
3084	amn=MIN(ts(i,is_ter),amn)
3085	amx=MAX(ts(i,is_ter),amx)
3086	ENDDO
3087	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
3088	ENDIF
3089
3090	!jg ?
3091	!!$!
3092	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
3093	!!$! sub-surfaces is distributed.
3094	!!$!
3095	!!$ DO nsrf = 1, nbsrf
3096	!!$ DO i = 1, klon
3097	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
3098	!!$ ts(i,nsrf) = zxtsol(i)
3099	!!$ t2m(i,nsrf) = zt2m(i)
3100	!!$ q2m(i,nsrf) = zq2m(i)
3101	!!$ u10m(i,nsrf) = zu10m(i)
3102	!!$ v10m(i,nsrf) = zv10m(i)
3103	!!$
3104	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
3105	!!$ pblh(i,nsrf) = s_pblh(i)
3106	!!$ plcl(i,nsrf) = s_plcl(i)
3107	!!$ capCL(i,nsrf) = s_capCL(i)
3108	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
3109	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
3110	!!$ pblT(i,nsrf) = s_pblT(i)
3111	!!$ therm(i,nsrf) = s_therm(i)
3112	!!$ trmb1(i,nsrf) = s_trmb1(i)
3113	!!$ trmb2(i,nsrf) = s_trmb2(i)
3114	!!$ trmb3(i,nsrf) = s_trmb3(i)
3115	!!$ ENDIF
3116	!!$ ENDDO
3117	!!$ ENDDO
3118
3119
3120	DO i = 1, klon
3121	fder(i) = - 4.0RSIGMAzxtsol(i)**3
3122	ENDDO
3123
3124	zxqsurf(:) = 0.0
3125	zxsnow(:) = 0.0
3126	DO nsrf = 1, nbsrf
3127	DO i = 1, klon
3128	zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
3129	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
3130	ENDDO
3131	ENDDO
3132
3133	! Premier niveau de vent sortie dans physiq.F
3134	zu1(:) = u(:,1)
3135	zv1(:) = v(:,1)
3136
3137	END SUBROUTINE pbl_surface
3138	!
3139	!****************************************************************************************
3140	!
3141	SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
3142
3143	USE indice_sol_mod
3144
3145	INCLUDE "dimsoil.h"
3146
3147	! Ouput variables
3148	!****************************************************************************************
3149	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
3150	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
3151	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
3152	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
3153
3154
3155	!****************************************************************************************
3156	! Return module variables for writing to restart file
3157	!
3158	!****************************************************************************************
3159	fder_rst(:) = fder(:)
3160	snow_rst(:,:) = snow(:,:)
3161	qsurf_rst(:,:) = qsurf(:,:)
3162	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
3163
3164	!****************************************************************************************
3165	! Deallocate module variables
3166	!
3167	!****************************************************************************************
3168	! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
3169	IF (ALLOCATED(fder)) DEALLOCATE(fder)
3170	IF (ALLOCATED(snow)) DEALLOCATE(snow)
3171	IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
3172	IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
3173
3174	!jyg<
3175	!****************************************************************************************
3176	! Deallocate variables for pbl splitting
3177	!
3178	!****************************************************************************************
3179
3180	CALL wx_pbl_final
3181	!>jyg
3182
3183	END SUBROUTINE pbl_surface_final
3184	!
3185	!****************************************************************************************
3186	!
3187
3188	!albedo SB >>>
3189	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
3190	evap, z0m, z0h, agesno, &
3191	tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)
3192	!albedo SB <<<
3193	! Give default values where new fraction has appread
3194
3195	USE indice_sol_mod
3196	use phys_state_var_mod, only: s_int, ds_ns, dt_ns
3197	use config_ocean_skin_m, only: activate_ocean_skin
3198
3199	INCLUDE "dimsoil.h"
3200	INCLUDE "clesphys.h"
3201	INCLUDE "compbl.h"
3202
3203	! Input variables
3204	!****************************************************************************************
3205	INTEGER, INTENT(IN) :: itime
3206	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
3207
3208	! InOutput variables
3209	!****************************************************************************************
3210	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
3211	!albedo SB >>>
3212	REAL, DIMENSION(klon,nsw,nbsrf), INTENT(INOUT) :: alb_dir, alb_dif
3213	INTEGER :: k
3214	!albedo SB <<<
3215	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m
3216	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: evap, agesno
3217	REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT) :: z0m,z0h
3218	REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
3219
3220	! Local variables
3221	!****************************************************************************************
3222	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
3223	CHARACTER(len=80) :: abort_message
3224	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
3225	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
3226	!
3227	! All at once !!
3228	!****************************************************************************************
3229
3230	DO nsrf = 1, nbsrf
3231	! First decide complement sub-surfaces
3232	SELECT CASE (nsrf)
3233	CASE(is_oce)
3234	nsrf_comp1=is_sic
3235	nsrf_comp2=is_ter
3236	nsrf_comp3=is_lic
3237	CASE(is_sic)
3238	nsrf_comp1=is_oce
3239	nsrf_comp2=is_ter
3240	nsrf_comp3=is_lic
3241	CASE(is_ter)
3242	nsrf_comp1=is_lic
3243	nsrf_comp2=is_oce
3244	nsrf_comp3=is_sic
3245	CASE(is_lic)
3246	nsrf_comp1=is_ter
3247	nsrf_comp2=is_oce
3248	nsrf_comp3=is_sic
3249	END SELECT
3250
3251	! Initialize all new fractions
3252	DO i=1, klon
3253	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
3254
3255	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
3256	! Use the complement sub-surface, keeping the continents unchanged
3257	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
3258	evap(i,nsrf) = evap(i,nsrf_comp1)
3259	z0m(i,nsrf) = z0m(i,nsrf_comp1)
3260	z0h(i,nsrf) = z0h(i,nsrf_comp1)
3261	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
3262	!albedo SB >>>
3263	DO k=1,nsw
3264	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp1)
3265	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp1)
3266	ENDDO
3267	!albedo SB <<<
3268	ustar(i,nsrf) = ustar(i,nsrf_comp1)
3269	u10m(i,nsrf) = u10m(i,nsrf_comp1)
3270	v10m(i,nsrf) = v10m(i,nsrf_comp1)
3271	IF (iflag_pbl > 1) THEN
3272	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
3273	ENDIF
3274	mfois(nsrf) = mfois(nsrf) + 1
3275	! F. Codron sensible default values for ocean and sea ice
3276	IF (nsrf.EQ.is_oce) THEN
3277	tsurf(i,nsrf) = 271.35
3278	! (temperature of sea water under sea ice, so that
3279	! is also the temperature of appearing sea water)
3280	DO k=1,nsw
3281	alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
3282	alb_dif(i,k,nsrf) = 0.06
3283	ENDDO
3284	if (activate_ocean_skin >= 1) then
3285	if (activate_ocean_skin == 2 &
3286	.and. type_ocean == "couple") s_int(i) = 35.
3287	ds_ns(i) = 0.
3288	dt_ns(i) = 0.
3289	end if
3290	ELSE IF (nsrf.EQ.is_sic) THEN
3291	tsurf(i,nsrf) = 271.35
3292	! (Temperature at base of sea ice. Surface
3293	! temperature could be higher, up to 0 Celsius
3294	! degrees. We set it to -1.8 Celsius degrees for
3295	! consistency with the ocean slab model.)
3296	DO k=1,nsw
3297	alb_dir(i,k,nsrf) = 0.3 ! thin ice
3298	alb_dif(i,k,nsrf) = 0.3
3299	ENDDO
3300	ENDIF
3301	ELSE
3302	! The continents have changed. The new fraction receives the mean sum of the existent fractions
3303	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3304	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3305	z0m(i,nsrf) = z0m(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0m(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3306	z0h(i,nsrf) = z0h(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0h(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3307	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3308	!albedo SB >>>
3309	DO k=1,nsw
3310	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3311	alb_dir(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3312	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3313	alb_dif(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3314	ENDDO
3315	!albedo SB <<<
3316	ustar(i,nsrf) = ustar(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3317	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3318	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3319	IF (iflag_pbl > 1) THEN
3320	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3321	ENDIF
3322
3323	! Security abort. This option has never been tested. To test, comment the following line.
3324	! abort_message='The fraction of the continents have changed!'
3325	! CALL abort_physic(modname,abort_message,1)
3326	nfois(nsrf) = nfois(nsrf) + 1
3327	ENDIF
3328	snow(i,nsrf) = 0.
3329	agesno(i,nsrf) = 0.
3330	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
3331	ELSE
3332	pfois(nsrf) = pfois(nsrf)+ 1
3333	ENDIF
3334	ENDDO
3335
3336	ENDDO
3337
3338	END SUBROUTINE pbl_surface_newfrac
3339	!
3340	!****************************************************************************************
3341	!
3342	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format