Context Navigation

pbl_surface_mod.F90 @ 2930

Last change on this file since 2930 was 2898, checked in by fhourdin, 7 years ago
Corrections de bugs sur les deux dernieres commissions
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:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 127.8 KB

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

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

Download in other formats:

Original Format