Context Navigation

pbl_surface_mod.F90 @ 4359

Last change on this file since 4359 was 3784, checked in by evignon, 4 years ago
correction v3783 pour convergence et compilation ancienne physique, Etienne aide par Ehouarn
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: 132.7 KB

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

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

Download in other formats:

Original Format