Context Navigation

pbl_surface_mod.F90 @ 3817

Last change on this file since 3817 was 3817, checked in by musat, 3 years ago
Nouvelle formulation des calculs a 2m
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: 137.3 KB

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

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

Download in other formats:

Original Format