Context Navigation

pbl_surface_mod.F90 @ 3605

Last change on this file since 3605 was 3605, checked in by lguez, 4 years ago
Merge revisions 3427:3600 of trunk into branch `Ocean_skin`
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:keywords set to `Author Date Id Revision`
File size: 133.9 KB

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

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

Download in other formats:

Original Format