Context Navigation

pbl_surface_mod.F90 @ 2182

Last change on this file since 2182 was 2182, checked in by Laurent Fairhead, 9 years ago
Suite aux commissions de Jean-Yves, la mise en commentaire n'est plus nécéssaire
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: 123.4 KB

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

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

Download in other formats:

Original Format