Context Navigation

pbl_surface_mod.F90 @ 2227

Last change on this file since 2227 was 2227, checked in by Laurent Fairhead, 9 years ago

New ocean albedo.

To activate the new scheme, put iflag_albedo=1 in physiq.def

To activate chlorophyll concentration effect on albedo,
put ok_chlorophyll=y in def file

and download file named chlorophyll.nc
chlorophyll.nc has the same dimension as the model grid with 12 months data,
(i=lon, j=lat, L=1:12) and can be degraded from the original file of dimension
i=1:4320 , j=1:2160 , L=1:12
ada:/workgpfs/rech/psl/rpsl949/clima/chlor_seasonal_clim_seawifs.nc

For 96X96 resolution, chlorophyll.nc file is in
ada:/workgpfs/rech/psl/rpsl949/clima/chlorophyll.nc

Baek

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: 126.6 KB

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

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

Download in other formats:

Original Format