Context Navigation

pbl_surface_mod.F90 @ 3379

Last change on this file since 3379 was 3198, checked in by fhourdin, 6 years ago

Retour vers l'insensibilite au decoupage en sous domaine.
Les routines gwd_rando incluait le calcul de niveaux de reference
sur la base d'un profile pris au milieu du domaine (en klon/2).
Rempace par un test en presnivs.

Une autre intercation entre routines concernant la tke a fait apparaitre
que la tke n'était pas passee correctement au niveau klev+1 au moment
du regroupement des mailles sous les sous surface.

Ces changements garantissent la convergence numerique si
addtkeoro=0
iflag_pbl<12
et
ok_gwd_rando=n
La convergence n'est pas garantie pour les dernieres versions des physiq.def
mais les differences devraient etre mineures.

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

File size: 130.5 KB

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

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

Download in other formats:

Original Format