Context Navigation

pbl_surface_mod.F90 @ 3811

Last change on this file since 3811 was 3798, checked in by lguez, 4 years ago
Sync latest trunk changes to Ocean_skin
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: 135.1 KB

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

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

Download in other formats:

Original Format