Context Navigation

pbl_surface_mod.F90 @ 5441

Last change on this file since 5441 was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF

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.0 KB

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

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

Download in other formats:

Original Format