Context Navigation

pbl_surface_mod.F90 @ 3626

Last change on this file since 3626 was 3626, checked in by lguez, 4 years ago
Bug fix: use `missing_val` from XIOS for output with XIOS, else the operations are not done propoerly in XIOS.
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: 133.9 KB

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

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

Download in other formats:

Original Format