Context Navigation

pbl_surface_mod.F90 @ 4584

Last change on this file since 4584 was 4584, checked in by Laurent Fairhead, 12 months ago
Merged trunk revisions from r4443 to r4582 (HEAD) into branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 153.1 KB

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

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

Download in other formats:

Original Format