Context Navigation

pbl_surface_mod.F90 @ 3458

Last change on this file since 3458 was 3458, checked in by lguez, 5 years ago

Introduce variable activate_ocean_skin in module config_ocean_skin_m.

Bug fix in phys_state_var_end: we need to deallocate variables for
lmdz1d (although it is useless for a 3D run).

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

File size: 133.8 KB

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

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

Download in other formats:

Original Format