Context Navigation

pbl_surface_mod.F90 @ 2298

Last change on this file since 2298 was 2298, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes -r2237:2291 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 123.7 KB

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

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

Download in other formats:

Original Format