Context Navigation

pbl_surface_mod.F90 @ 3761

Last change on this file since 3761 was 3756, checked in by oboucher, 4 years ago
The fraction of SW down radiation that is diffuse is extracted from radlwsw.F90 and passed on to pbl_surface for inclusion in fields_out for transfer to ORCHIDEE if it requested in input file coupling_fields.def. Hence there is nothing automatic here. The fraction of SW down radiation is also added as a diagnostic. To satisfy the case when this new quantity is used in ORCHIDEE, then it is added in the restart (and read in case it is there).
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: 132.4 KB

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

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

Download in other formats:

Original Format