Context Navigation

pbl_surface_mod.F90 @ 3744

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

Store delta_sst instead of sst_nff

Store as a state variable the difference between ocean-air interface
temperature and bulk SST instead of sst_nff, which can be either
interface temperature or bulk SST. This is clearer. Also, it is
analoguous to what we will do with salinity.

So replace the two dummy arguments tsurf_in and sst_nff of
procedure cpl_send_ocean_fields by a single dummy argument
delta_sst. Replace dummy argument sst_nff of procedures
ocean_cpl_noice and surf_ocean by dummy argument
delta_sst. Replace variable sst_nff of module phys_state_var_mod
by variable delta_sst. Rename local variable ysst_nff of procedure
pbl_surface to ydelta_sst. Set variable delta_sst of module
phys_state_var_mod to 0 for an appearing ocean fraction and a
missing startup field. Replace variable o_sst_nff of module
phys_output_ctrlout_mod by variable o_delta_sst.

Rename variables cpl_delta_temp and cpl_delta_temp_2D of module
cpl_mod to cpl_delta_sst and cpl_delta_sst_2D, clearer. Rename
variable ids_delta_temp of module oasis to ids_delta_sst. Change
infosend(ids_delta_temp)%name to "CODELSST".

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

File size: 134.4 KB

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

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

Download in other formats:

Original Format