Context Navigation

pbl_surface_mod.F90 @ 3405

Last change on this file since 3405 was 3356, checked in by Laurent Fairhead, 6 years ago
First attempt at merging with trunk
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: 130.4 KB

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

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

Download in other formats:

Original Format