Context Navigation

pbl_surface_mod.F90 @ 4531

Last change on this file since 4531 was 4531, checked in by evignon, 14 months ago
petite reecriture pour neige soufflee dans pbl_surface
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: 152.1 KB

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

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

Download in other formats:

Original Format