Context Navigation

pbl_surface_mod.F90 @ 2209

Last change on this file since 2209 was 2209, checked in by Ehouarn Millour, 9 years ago
Update of the slab ocean by Francis Codron. There are now 3 possibilities for the "version_ocean" slab type: sicOBS = prescribed ice fraction. Water temperature nearby is set to -1.8°C and cannot become lower. sicNO = ignore sea ice. One can prescribe a fraction, but the nearby ocean evolves freely, depending on surface fluxes: temperature can go below freezing point or above... sicINT = interactive sea ice. EM
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.1 KB

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

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

Download in other formats:

Original Format