Context Navigation

pbl_surface_mod.F90 @ 2241

Last change on this file since 2241 was 2241, checked in by fhourdin, 9 years ago

Nettoyage des anciens albedo. Elimination de alb1 et alb2
dans pbl_surface (il s'agissait de commentaires) et dans
le etats de démarrage.

Some cleaning of old albedo specification (alb1/alb2)

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

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

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

Download in other formats:

Original Format