Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

pbl_surface_mod.F90

Last change on this file was 4743, checked in by Laurent Fairhead, 8 months ago
Merge of ACC branch with 4740 revision from trunk
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: 153.0 KB

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

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

Download in other formats:

Original Format