Context Navigation

pbl_surface_mod.F90 @ 2232

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

1) Fusion des procédures clcdrag.F90 et coefcdrag.F90

dans une procédure unique cdrag.F90.
Les deux anciennes sont obsolètes mais maintenues dans le
code quelques mois pour d'éventuels tests.

2) modification du nom du module arth.F90 en arth_m.F90 pour

avoir le même nom de fichier et de module (règle adoptée
pour LMDZ et utilisée par makelmdz).

1) merging clcdrag.F90 and coefcdrag.F90 in one procedure cdrag.F90
2) renaming arth.F90 in arth_m.F90

Wang/FH

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 127.9 KB

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

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

Download in other formats:

Original Format