Context Navigation

pbl_surface_mod.F90 @ 4947

Last change on this file since 4947 was 4485, checked in by evignon, 20 months ago
premier commit pour l'ajout de la neige soufflee sur la nouvelle branche
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.0 KB

Line
1	!
2	! $Id: pbl_surface_mod.F90 4485 2023-03-30 07:07:40Z 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 diffuvisity coefficient for
1911	! suspended particles is larger than Km by a factor zeta_bs
1912	! which is equal to 3 by default
1913	ycoefqbs=ycoefm*zeta_bs
1914	CALL climb_qbs_down(knon, ycoefqbs, ypaprs, ypplay, &
1915	ydelp, yt, yqbs, dtime, &
1916	CcoefQBS, DcoefQBS, &
1917	Kcoef_qbs, gama_qbs, &
1918	AcoefQBS, BcoefQBS)
1919	ENDIF
1920
1921	!****************************************************************************************
1922	! 9) Small calculations
1923	!
1924	!****************************************************************************************
1925
1926	! - Reference pressure is given the values at surface level
1927	ypsref(:) = ypaprs(:,1)
1928
1929	! - CO2 field on 2D grid to be sent to ORCHIDEE
1930	! Transform to compressed field
1931	IF (carbon_cycle_cpl) THEN
1932	DO i=1,knon
1933	r_co2_ppm(i) = co2_send(ni(i))
1934	ENDDO
1935	ELSE
1936	r_co2_ppm(:) = co2_ppm ! Constant field
1937	ENDIF
1938
1939	!!! nrlmd le 02/05/2011 -----------------------On raccorde les 2 colonnes dans la couche 1
1940	!----------------------------------------------------------------------------------------
1941	!!! jyg le 07/02/2012
1942	!!! jyg le 01/02/2017
1943	IF (iflag_split .eq. 0) THEN
1944	yt1(:) = yt(:,1)
1945	yq1(:) = yq(:,1)
1946	ELSE IF (iflag_split .ge. 1) THEN
1947	!
1948	! Cdragq computation
1949	! ------------------
1950	!******************************************************************************
1951	! Cdragq computed from cdrag
1952	! The difference comes only from a factor (f_z0qh_oce) on z0, so that
1953	! it can be computed inside wx_pbl0_merge
1954	! More complicated appraches may require the propagation through
1955	! pbl_surface of an independant cdragq variable.
1956	!******************************************************************************
1957	!
1958	IF ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) THEN
1959	! Si on suit les formulations par exemple de Tessel, on
1960	! a z0h=0.4nu/u, z0q=0.62nu/u, d'ou f_z0qh_oce=0.62/0.4=1.55
1961	!! ycdragq_x(1:knon)=ycdragh_x(1:knon)* &
1962	!! log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
1963	!! ycdragq_w(1:knon)=ycdragh_w(1:knon)* &
1964	!! log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
1965	!
1966	DO j = 1,knon
1967	z1lay = zgeo1(j)/RG
1968	fact_cdrag = log(z1lay/yz0h(j))/log(z1lay/(f_z0qh_oce*yz0h(j)))
1969	ycdragq_x(j)=ycdragh_x(j)*fact_cdrag
1970	ycdragq_w(j)=ycdragh_w(j)*fact_cdrag
1971	!! Print *,'YYYYpbl0: fact_cdrag ', fact_cdrag
1972	ENDDO ! j = 1,knon
1973	!
1974	!! Print *,'YYYYpbl0: z1lay, yz0h, f_z0qh_oce, ycdragh_w, ycdragq_w ', &
1975	!! z1lay, yz0h(1:knon), f_z0qh_oce, ycdragh_w(1:knon), ycdragq_w(1:knon)
1976	ELSE
1977	ycdragq_x(1:knon)=ycdragh_x(1:knon)
1978	ycdragq_w(1:knon)=ycdragh_w(1:knon)
1979	ENDIF ! ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce)
1980	!
1981	CALL wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, ywake_s, &
1982	yts, y_delta_tsurf, ygustiness, &
1983	yt_x, yt_w, yq_x, yq_w, &
1984	yu_x, yu_w, yv_x, yv_w, &
1985	ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
1986	ycdragm_x, ycdragm_w, &
1987	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
1988	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
1989	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
1990	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
1991	Kech_h_x, Kech_h_w, Kech_h &
1992	)
1993	CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, &
1994	BcoefQ_x, BcoefQ_w &
1995	)
1996	CALL wx_pbl0_merge(knon, ypplay, ypaprs, &
1997	ywake_s, ydTs0, ydqs0, &
1998	yt_x, yt_w, yq_x, yq_w, &
1999	yu_x, yu_w, yv_x, yv_w, &
2000	ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
2001	ycdragm_x, ycdragm_w, &
2002	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
2003	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
2004	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
2005	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
2006	AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
2007	BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
2008	ycdragh, ycdragq, ycdragm, &
2009	yt1, yq1, yu1, yv1 &
2010	)
2011	IF (iflag_split .eq. 2 .AND. nsrf .ne. is_oce) THEN
2012	CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
2013	ywake_s, ybeta, ywake_cstar, ywake_dens, &
2014	AcoefH_x, AcoefH_w, &
2015	BcoefH_x, BcoefH_w, &
2016	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2017	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2018	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2019	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2020	yg_T, yg_Q, &
2021	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2022	ydTs_ins, ydqs_ins &
2023	)
2024	ELSE !
2025	AcoefH(:) = AcoefH_0(:)
2026	AcoefQ(:) = AcoefQ_0(:)
2027	BcoefH(:) = BcoefH_0(:)
2028	BcoefQ(:) = BcoefQ_0(:)
2029	yg_T(:) = 0.
2030	yg_Q(:) = 0.
2031	yGamma_dTs_phiT(:) = 0.
2032	yGamma_dQs_phiQ(:) = 0.
2033	ydTs_ins(:) = 0.
2034	ydqs_ins(:) = 0.
2035	ENDIF ! (iflag_split .eq. 2)
2036	ENDIF ! (iflag_split .eq.0)
2037	!!!
2038	IF (prt_level >=10) THEN
2039	PRINT *,'pbl_surface (merge->): yt(1,:) ',yt(1,:)
2040	PRINT *,'pbl_surface (merge->): yq(1,:) ',yq(1,:)
2041	PRINT *,'pbl_surface (merge->): yu(1,:) ',yu(1,:)
2042	PRINT *,'pbl_surface (merge->): yv(1,:) ',yv(1,:)
2043	PRINT *,'pbl_surface (merge->): AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1) ', &
2044	AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1)
2045	PRINT *,'pbl_surface (merge->): BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1) ', &
2046	BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1)
2047
2048	ENDIF
2049
2050	! Save initial value of z0h for use in evappot (z0h wiil be computed again in the surface models)
2051	yz0h_old(1:knon) = yz0h(1:knon)
2052	!
2053	!****************************************************************************************
2054	!
2055	! Calulate t2m and q2m for the case of calculation at land grid points
2056	! t2m and q2m are needed as input to ORCHIDEE
2057	!
2058	!****************************************************************************************
2059	IF (nsrf == is_ter) THEN
2060
2061	DO i = 1, knon
2062	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
2063	* (ypaprs(i,1)-ypplay(i,1))
2064	ENDDO
2065
2066	! Calculate the temperature et relative humidity at 2m and the wind at 10m
2067	IF (iflag_new_t2mq2m==1) THEN
2068	CALL stdlevvarn(klon, knon, is_ter, zxli, &
2069	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
2070	yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
2071	yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
2072	yn2mout(:, nsrf, :))
2073	ELSE
2074	CALL stdlevvar(klon, knon, is_ter, zxli, &
2075	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
2076	yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
2077	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
2078	ENDIF
2079
2080	ENDIF
2081
2082	!****************************************************************************************
2083	!
2084	! 10) Switch according to current surface
2085	! It is necessary to start with the continental surfaces because the ocean
2086	! needs their run-off.
2087	!
2088	!****************************************************************************************
2089	SELECT CASE(nsrf)
2090
2091	CASE(is_ter)
2092	! print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)
2093	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
2094	rlon, rlat, yrmu0, &
2095	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
2096	!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
2097	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
2098	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2099	AcoefU, AcoefV, BcoefU, BcoefV, &
2100	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
2101	ylwdown, yq2m, yt2m, &
2102	ysnow, yqsol, yagesno, ytsoil, &
2103	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,yfluxbs,&
2104	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
2105	y_flux_u1, y_flux_v1, &
2106	yveget,ylai,yheight )
2107
2108	!FC quid qd yveget ylai yheight ne sont pas definit
2109	!FC yveget,ylai,yheight, &
2110	IF (ifl_pbltree .ge. 1) THEN
2111	CALL freinage(knon, yu, yv, yt, &
2112	! yveget,ylai, yheight,ypaprs,ypplay,y_d_u_frein,y_d_v_frein)
2113	yveget,ylai, yheight,ypaprs,ypplay,y_treedrg, y_d_u_frein,y_d_v_frein)
2114	ENDIF
2115
2116
2117	! Special DICE MPL 05082013 puis BOMEX
2118	IF (ok_prescr_ust) THEN
2119	DO j=1,knon
2120	! ysnow(:)=0.
2121	! yqsol(:)=0.
2122	! yagesno(:)=50.
2123	! ytsoil(:,:)=300.
2124	! yz0_new(:)=0.001
2125	! yevap(:)=flat/RLVTT
2126	! yfluxlat(:)=-flat
2127	! yfluxsens(:)=-fsens
2128	! yqsurf(:)=0.
2129	! ytsurf_new(:)=tg
2130	! y_dflux_t(:)=0.
2131	! y_dflux_q(:)=0.
2132	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)
2133	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)
2134	ENDDO
2135	ENDIF
2136
2137	CASE(is_lic)
2138	! Martin
2139
2140	IF (landice_opt .LT. 2) THEN
2141	! Land ice is treated by LMDZ and not by ORCHIDEE
2142	CALL surf_landice(itap, dtime, knon, ni, &
2143	rlon, rlat, debut, lafin, &
2144	yrmu0, ylwdown, yalb, zgeo1, &
2145	ysolsw, ysollw, yts, ypplay(:,1), &
2146	ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
2147	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2148	AcoefU, AcoefV, BcoefU, BcoefV, &
2149	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
2150	ysnow, yqsurf, yqsol,yqbs1, yagesno, &
2151	ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
2152	yfluxbs, ytsurf_new, y_dflux_t, y_dflux_q, &
2153	yzmea, yzsig, ycldt, &
2154	ysnowhgt, yqsnow, ytoice, ysissnow, &
2155	yalb3_new, yrunoff, &
2156	y_flux_u1, y_flux_v1)
2157
2158	!jyg<
2159	!! alb3_lic(:)=0.
2160	!>jyg
2161	DO j = 1, knon
2162	i = ni(j)
2163	alb3_lic(i) = yalb3_new(j)
2164	snowhgt(i) = ysnowhgt(j)
2165	qsnow(i) = yqsnow(j)
2166	to_ice(i) = ytoice(j)
2167	sissnow(i) = ysissnow(j)
2168	runoff(i) = yrunoff(j)
2169	ENDDO
2170	! Martin
2171	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2172	IF (ok_prescr_ust) THEN
2173	DO j=1,knon
2174	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)
2175	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)
2176	ENDDO
2177	ENDIF
2178
2179	END IF
2180
2181	CASE(is_oce)
2182	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
2183	ywindsp, rmu0, yfder, yts, &
2184	itap, dtime, jour, knon, ni, &
2185	!!jyg ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
2186	ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt(:,1), yq(:,1),& ! ym missing init
2187	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2188	AcoefU, AcoefV, BcoefU, BcoefV, &
2189	ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
2190	ysnow, yqsurf, yagesno, &
2191	yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
2192	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
2193	y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
2194	yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
2195	ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss)
2196	IF (prt_level >=10) THEN
2197	print *,'arg de surf_ocean: ycdragh ',ycdragh
2198	print *,'arg de surf_ocean: ycdragm ',ycdragm
2199	print *,'arg de surf_ocean: yt ', yt
2200	print *,'arg de surf_ocean: yq ', yq
2201	print *,'arg de surf_ocean: yts ', yts
2202	print *,'arg de surf_ocean: AcoefH ',AcoefH
2203	print *,'arg de surf_ocean: AcoefQ ',AcoefQ
2204	print *,'arg de surf_ocean: BcoefH ',BcoefH
2205	print *,'arg de surf_ocean: BcoefQ ',BcoefQ
2206	print *,'arg de surf_ocean: yevap ',yevap
2207	print *,'arg de surf_ocean: yfluxsens ',yfluxsens
2208	print *,'arg de surf_ocean: yfluxlat ',yfluxlat
2209	print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new
2210	ENDIF
2211	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2212	IF (ok_prescr_ust) THEN
2213	DO j=1,knon
2214	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)
2215	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)
2216	ENDDO
2217	ENDIF
2218
2219	CASE(is_sic)
2220	CALL surf_seaice( &
2221	!albedo SB >>>
2222	rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, &
2223	!albedo SB <<<
2224	itap, dtime, jour, knon, ni, &
2225	lafin, &
2226	!!jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
2227	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
2228	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2229	AcoefU, AcoefV, BcoefU, BcoefV, &
2230	ypsref, yu1, yv1, ygustiness, pctsrf, &
2231	ysnow, yqsurf, yqsol, yagesno, ytsoil, &
2232	!albedo SB >>>
2233	yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
2234	!albedo SB <<<
2235	ytsurf_new, y_dflux_t, y_dflux_q, &
2236	y_flux_u1, y_flux_v1)
2237
2238	! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2239	IF (ok_prescr_ust) THEN
2240	DO j=1,knon
2241	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)
2242	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)
2243	ENDDO
2244	ENDIF
2245
2246	CASE DEFAULT
2247	WRITE(lunout,*) 'Surface index = ', nsrf
2248	abort_message = 'Surface index not valid'
2249	CALL abort_physic(modname,abort_message,1)
2250	END SELECT
2251
2252
2253	!****************************************************************************************
2254	! 11) - Calcul the increment of surface temperature
2255	!
2256	!****************************************************************************************
2257
2258	IF (evap0>=0.) THEN
2259	yevap(:)=evap0
2260	yevap(:)=RLVTT*evap0
2261	ENDIF
2262
2263	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
2264
2265	!****************************************************************************************
2266	!
2267	! 12) "La remontee" - "The uphill"
2268	!
2269	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
2270	! for X=H, Q, U and V, for all vertical levels.
2271	!
2272	!****************************************************************************************
2273	!!
2274	!!!
2275	!!! jyg le 10/04/2013 et EV 10/2020
2276
2277	IF (ok_forc_tsurf) THEN
2278	DO j=1,knon
2279	ytsurf_new(j)=tg
2280	y_d_ts(j) = ytsurf_new(j) - yts(j)
2281	ENDDO
2282	ENDIF ! ok_forc_tsurf
2283
2284	!!!
2285	IF (ok_flux_surf) THEN
2286	IF (prt_level >=10) THEN
2287	PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
2288	ENDIF
2289	y_flux_t1(:) = fsens
2290	y_flux_q1(:) = flat/RLVTT
2291	yfluxlat(:) = flat
2292	!
2293	!! Test sur iflag_split retire le 2/02/2018, sans vraiment comprendre la raison de ce test. (jyg)
2294	!! IF (iflag_split .eq.0) THEN
2295	DO j=1,knon
2296	Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)2+yv(j,1)2)) * &
2297	ypplay(j,1)/(RD*yt(j,1))
2298	ENDDO
2299	!! ENDIF ! (iflag_split .eq.0)
2300
2301	DO j = 1, knon
2302	yt1_new=(1./RCPD)(AcoefH(j)+BcoefH(j)y_flux_t1(j)*dtime)
2303	ytsurf_new(j)=yt1_new-y_flux_t1(j)/(Kech_h(j)*RCPD)
2304	ENDDO
2305
2306	DO j=1,knon
2307	y_d_ts(j) = ytsurf_new(j) - yts(j)
2308	ENDDO
2309
2310	ELSE ! (ok_flux_surf)
2311	DO j=1,knon
2312	y_flux_t1(j) = yfluxsens(j)
2313	y_flux_q1(j) = -yevap(j)
2314	ENDDO
2315	ENDIF ! (ok_flux_surf)
2316
2317	! flux of blowing snow at the first level
2318	IF (ok_bs) THEN
2319	DO j=1,knon
2320	y_flux_bs(j)=yfluxbs(j)
2321	ENDDO
2322	ENDIF
2323	!
2324	! ------------------------------------------------------------------------------
2325	! 12a) Splitting
2326	! ------------------------------------------------------------------------------
2327
2328	IF (iflag_split .GE. 1) THEN
2329	!
2330	IF (nsrf .ne. is_oce) THEN
2331	!
2332	! Compute potential evaporation and aridity factor (jyg, 20200328)
2333	ybeta_prev(:) = ybeta(:)
2334	DO j = 1, knon
2335	yqa(j) = AcoefQ(j) - BcoefQ(j)yevap(j)dtime
2336	ENDDO
2337	!
2338	CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot)
2339	!
2340	ybeta(1:knon) = min(yevap(1:knon)/yevap_pot(1:knon), 1.)
2341
2342	IF (prt_level >=10) THEN
2343	DO j=1,knon
2344	print*,'y_flux_t1,yfluxlat,wakes' &
2345	& , y_flux_t1(j), yfluxlat(j), ywake_s(j)
2346	print*,'beta_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j)
2347	print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
2348	ENDDO
2349	ENDIF ! (prt_level >=10)
2350	!
2351	! Second call to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account
2352	! the update of the aridity coeficient beta.
2353	!
2354	CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta, &
2355	BcoefQ_x, BcoefQ_w &
2356	)
2357	CALL wx_pbl0_merge(knon, ypplay, ypaprs, &
2358	ywake_s, ydTs0, ydqs0, &
2359	yt_x, yt_w, yq_x, yq_w, &
2360	yu_x, yu_w, yv_x, yv_w, &
2361	ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
2362	ycdragm_x, ycdragm_w, &
2363	AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
2364	AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
2365	BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
2366	BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
2367	AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
2368	BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
2369	ycdragh, ycdragq, ycdragm, &
2370	yt1, yq1, yu1, yv1 &
2371	)
2372	IF (iflag_split .eq. 2) THEN
2373	CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
2374	ywake_s, ybeta, ywake_cstar, ywake_dens, &
2375	AcoefH_x, AcoefH_w, &
2376	BcoefH_x, BcoefH_w, &
2377	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2378	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2379	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2380	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2381	yg_T, yg_Q, &
2382	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2383	ydTs_ins, ydqs_ins &
2384	)
2385	ELSE !
2386	AcoefH(:) = AcoefH_0(:)
2387	AcoefQ(:) = AcoefQ_0(:)
2388	BcoefH(:) = BcoefH_0(:)
2389	BcoefQ(:) = BcoefQ_0(:)
2390	yg_T(:) = 0.
2391	yg_Q(:) = 0.
2392	yGamma_dTs_phiT(:) = 0.
2393	yGamma_dQs_phiQ(:) = 0.
2394	ydTs_ins(:) = 0.
2395	ydqs_ins(:) = 0.
2396	ENDIF ! (iflag_split .eq. 2)
2397	!
2398	ELSE ! (nsrf .ne. is_oce)
2399	ybeta(1:knon) = 1.
2400	yevap_pot(1:knon) = yevap(1:knon)
2401	AcoefH(:) = AcoefH_0(:)
2402	AcoefQ(:) = AcoefQ_0(:)
2403	BcoefH(:) = BcoefH_0(:)
2404	BcoefQ(:) = BcoefQ_0(:)
2405	yg_T(:) = 0.
2406	yg_Q(:) = 0.
2407	yGamma_dTs_phiT(:) = 0.
2408	yGamma_dQs_phiQ(:) = 0.
2409	ydTs_ins(:) = 0.
2410	ydqs_ins(:) = 0.
2411	ENDIF ! (nsrf .ne. is_oce)
2412	!
2413	CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, &
2414	yg_T, yg_Q, &
2415	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2416	ydTs_ins, ydqs_ins, &
2417	y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
2418	!!!! HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
2419	phiQ0_b, phiT0_b, &
2420	y_flux_t1_x, y_flux_t1_w, &
2421	y_flux_q1_x, y_flux_q1_w, &
2422	y_flux_u1_x, y_flux_u1_w, &
2423	y_flux_v1_x, y_flux_v1_w, &
2424	yfluxlat_x, yfluxlat_w, &
2425	y_delta_qsats, &
2426	y_delta_tsurf_new, y_delta_qsurf &
2427	)
2428	!
2429	CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
2430	yTs, y_delta_tsurf, &
2431	yqsurf, yTsurf_new, &
2432	y_delta_tsurf_new, y_delta_qsats, &
2433	AcoefH_x, AcoefH_w, &
2434	BcoefH_x, BcoefH_w, &
2435	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2436	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2437	y_flux_t1, y_flux_q1, &
2438	y_flux_t1_x, y_flux_t1_w, &
2439	y_flux_q1_x, y_flux_q1_w)
2440	!
2441	IF (nsrf .ne. is_oce) THEN
2442	CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
2443	yTs, y_delta_tsurf, &
2444	yqsurf, yTsurf_new, &
2445	y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf, &
2446	AcoefH_x, AcoefH_w, &
2447	BcoefH_x, BcoefH_w, &
2448	AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0, &
2449	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2450	HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2451	phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2452	yg_T, yg_Q, &
2453	yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2454	ydTs_ins, ydqs_ins, &
2455	y_flux_t1, y_flux_q1, &
2456	y_flux_t1_x, y_flux_t1_w, &
2457	y_flux_q1_x, y_flux_q1_w )
2458	ENDIF ! (nsrf .ne. is_oce)
2459	!
2460	ELSE ! (iflag_split .ge. 1)
2461	ybeta(1:knon) = 1.
2462	yevap_pot(1:knon) = yevap(1:knon)
2463	ENDIF ! (iflag_split .ge. 1)
2464	!
2465	IF (prt_level >= 10) THEN
2466	print *,'pbl_surface, ybeta , yevap, yevap_pot ', &
2467	ybeta , yevap, yevap_pot
2468	ENDIF ! (prt_level >= 10)
2469	!
2470	!>jyg
2471	!
2472
2473	!!jyg!! A reprendre apres reflexion ===============================================
2474	!!jyg!!
2475	!!jyg!! DO j=1,knon
2476	!!jyg!!!!! nrlmd le 13/06/2011
2477	!!jyg!!
2478	!!jyg!!!----Diffusion dans le sol dans le cas continental seulement
2479	!!jyg!! IF (nsrf.eq.is_ter) THEN
2480	!!jyg!!!----Calcul du coefficient delta_coeff
2481	!!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)))
2482	!!jyg!!
2483	!!jyg!!! delta_coef(j)=dtime/(inertia*sqrt(tau_eq(j)))
2484	!!jyg!! delta_coef(j)=facteur*sqrt(tau_eq(j))/inertia
2485	!!jyg!!! delta_coef(j)=0.
2486	!!jyg!! ELSE
2487	!!jyg!! delta_coef(j)=0.
2488	!!jyg!! ENDIF
2489	!!jyg!!
2490	!!jyg!!!----Calcul de delta_tsurf
2491	!!jyg!! y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j)
2492	!!jyg!!
2493	!!jyg!!!----Si il n'y a pas des poches...
2494	!!jyg!! IF (wake_cstar(j).le.0.01) THEN
2495	!!jyg!! y_delta_tsurf(j)=0.
2496	!!jyg!! y_delta_flux_t1(j)=0.
2497	!!jyg!! ENDIF
2498	!!jyg!!
2499	!!jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle)
2500	!!jyg!!!!!!! jyg le 23/02/2012
2501	!!jyg!!!!!!!
2502	!!jyg!!!! ybeta(j)=y_flux_q1(j) / &
2503	!!jyg!!!! & (Kech_h(j)*(yq(j,1)-yqsatsurf(j)))
2504	!!jyg!!!!!! ybeta(j)=-1.*yevap(j) / &
2505	!!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)))
2506	!!jyg!!!!!!! fin jyg
2507	!!jyg!!!!!
2508	!!jyg!!
2509	!!jyg!! ENDDO
2510	!!jyg!!
2511	!!jyg!!!!! fin nrlmd le 13/06/2011
2512	!!jyg!!
2513	IF (iflag_split .ge. 1) THEN
2514	IF (prt_level >=10) THEN
2515	DO j = 1, knon
2516	print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
2517	print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
2518	print*,'t1x, t1w, t1, t1_ancien', &
2519	& yt_x(j,1), yt_w(j,1), yt(j,1), t(j,1)
2520	print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j)
2521	ENDDO
2522
2523	DO j=1,knon
2524	print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
2525	& , 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)
2526	print*,'beta, ytsurf_new ', ybeta(j), ytsurf_new(j)
2527	print*,'inertia, facteur, cstar', inertia, facteur,wake_cstar(j)
2528	ENDDO
2529	ENDIF ! (prt_level >=10)
2530
2531	!!! jyg le 07/02/2012
2532	ENDIF ! (iflag_split .ge.1)
2533	!!!
2534
2535	!!! jyg le 07/02/2012
2536	IF (iflag_split .eq.0) THEN
2537	!!!
2538	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2539	CALL climb_hq_up(knon, dtime, yt, yq, &
2540	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
2541	!!! jyg le 07/02/2012
2542	AcoefH, AcoefQ, BcoefH, BcoefQ, &
2543	CcoefH, CcoefQ, DcoefH, DcoefQ, &
2544	Kcoef_hq, gama_q, gama_h, &
2545	!!!
2546	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
2547	ELSE !(iflag_split .eq.0)
2548	CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
2549	y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &
2550	!!! nrlmd le 02/05/2011
2551	AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &
2552	CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
2553	Kcoef_hq_x, gama_q_x, gama_h_x, &
2554	!!!
2555	y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:))
2556	!
2557	CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
2558	y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &
2559	!!! nrlmd le 02/05/2011
2560	AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &
2561	CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
2562	Kcoef_hq_w, gama_q_w, gama_h_w, &
2563	!!!
2564	y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:))
2565	!!!
2566	ENDIF ! (iflag_split .eq.0)
2567	!!!
2568
2569	!!! jyg le 07/02/2012
2570	IF (iflag_split .eq.0) THEN
2571	!!!
2572	!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2573	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
2574	!!! jyg le 07/02/2012
2575	AcoefU, AcoefV, BcoefU, BcoefV, &
2576	CcoefU, CcoefV, DcoefU, DcoefV, &
2577	Kcoef_m, &
2578	!!!
2579	y_flux_u, y_flux_v, y_d_u, y_d_v)
2580	y_d_t_diss(:,:)=0.
2581	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2582	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2583	& ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
2584	& ,iflag_pbl)
2585	ENDIF
2586	! print*,'yamada_c OK'
2587
2588	ELSE !(iflag_split .eq.0)
2589	CALL climb_wind_up(knon, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &
2590	!!! nrlmd le 02/05/2011
2591	AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &
2592	CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
2593	Kcoef_m_x, &
2594	!!!
2595	y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)
2596	!
2597	y_d_t_diss_x(:,:)=0.
2598	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2599	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2600	& ,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 &
2601	,ycoefq_x,y_d_t_diss_x,yustar_x &
2602	& ,iflag_pbl)
2603	ENDIF
2604	! print*,'yamada_c OK'
2605
2606	CALL climb_wind_up(knon, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &
2607	!!! nrlmd le 02/05/2011
2608	AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &
2609	CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
2610	Kcoef_m_w, &
2611	!!!
2612	y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)
2613	!!!
2614	y_d_t_diss_w(:,:)=0.
2615	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2616	CALL yamada_c(knon,dtime,ypaprs,ypplay &
2617	& ,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 &
2618	,ycoefq_w,y_d_t_diss_w,yustar_w &
2619	& ,iflag_pbl)
2620	ENDIF
2621	! print*,'yamada_c OK'
2622	!
2623	IF (prt_level >=10) THEN
2624	print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &
2625	yfluxlat_x, yfluxlat_w
2626	ENDIF
2627	!
2628	ENDIF ! (iflag_split .eq.0)
2629
2630	IF (ok_bs) THEN
2631	CALL climb_qbs_up(knon, dtime, yqbs, &
2632	y_flux_bs, ypaprs, ypplay, &
2633	AcoefQBS, BcoefQBS, &
2634	CcoefQBS, DcoefQBS, &
2635	Kcoef_qbs, gama_qbs, &
2636	y_flux_qbs(:,:), y_d_qbs(:,:))
2637	ENDIF
2638
2639	!!!
2640	!!
2641	!! DO j = 1, knon
2642	!! y_dflux_t(j) = y_dflux_t(j) * ypct(j)
2643	!! y_dflux_q(j) = y_dflux_q(j) * ypct(j)
2644	!! ENDDO
2645	!!
2646	!****************************************************************************************
2647	! 13) Transform variables for output format :
2648	! - Decompress
2649	! - Multiply with pourcentage of current surface
2650	! - Cumulate in global variable
2651	!
2652	!****************************************************************************************
2653
2654
2655	!!! jyg le 07/02/2012
2656	IF (iflag_split.EQ.0) THEN
2657	!!!
2658	DO k = 1, klev
2659	DO j = 1, knon
2660	i = ni(j)
2661	y_d_t_diss(j,k) = y_d_t_diss(j,k) * ypct(j)
2662	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
2663	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
2664	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
2665	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
2666	!FC
2667	IF (nsrf .EQ. is_ter .and. ifl_pbltree .GE. 1) THEN
2668	! if (y_d_u_frein(j,k).ne.0. ) then
2669	! print, nsrf,'IS_TER ++', y_d_u_frein(j,k)ypct(j),y_d_u(j,k),j,k
2670	! ENDIF
2671	y_d_u(j,k) =y_d_u(j,k) + y_d_u_frein(j,k)*ypct(j)
2672	y_d_v(j,k) =y_d_v(j,k) + y_d_v_frein(j,k)*ypct(j)
2673	treedrg(i,k,nsrf)=y_treedrg(j,k)
2674	ELSE
2675	treedrg(i,k,nsrf)=0.
2676	ENDIF
2677	!FC
2678	flux_t(i,k,nsrf) = y_flux_t(j,k)
2679	flux_q(i,k,nsrf) = y_flux_q(j,k)
2680	flux_u(i,k,nsrf) = y_flux_u(j,k)
2681	flux_v(i,k,nsrf) = y_flux_v(j,k)
2682	ENDDO
2683	ENDDO
2684
2685	ELSE !(iflag_split .eq.0)
2686
2687	! Tendances hors poches
2688	DO k = 1, klev
2689	DO j = 1, knon
2690	i = ni(j)
2691	y_d_t_diss_x(j,k) = y_d_t_diss_x(j,k) * ypct(j)
2692	y_d_t_x(j,k) = y_d_t_x(j,k) * ypct(j)
2693	y_d_q_x(j,k) = y_d_q_x(j,k) * ypct(j)
2694	y_d_u_x(j,k) = y_d_u_x(j,k) * ypct(j)
2695	y_d_v_x(j,k) = y_d_v_x(j,k) * ypct(j)
2696
2697	flux_t_x(i,k,nsrf) = y_flux_t_x(j,k)
2698	flux_q_x(i,k,nsrf) = y_flux_q_x(j,k)
2699	flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
2700	flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)
2701	ENDDO
2702	ENDDO
2703
2704	! Tendances dans les poches
2705	DO k = 1, klev
2706	DO j = 1, knon
2707	i = ni(j)
2708	y_d_t_diss_w(j,k) = y_d_t_diss_w(j,k) * ypct(j)
2709	y_d_t_w(j,k) = y_d_t_w(j,k) * ypct(j)
2710	y_d_q_w(j,k) = y_d_q_w(j,k) * ypct(j)
2711	y_d_u_w(j,k) = y_d_u_w(j,k) * ypct(j)
2712	y_d_v_w(j,k) = y_d_v_w(j,k) * ypct(j)
2713
2714	flux_t_w(i,k,nsrf) = y_flux_t_w(j,k)
2715	flux_q_w(i,k,nsrf) = y_flux_q_w(j,k)
2716	flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
2717	flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)
2718	ENDDO
2719	ENDDO
2720
2721	! Flux, tendances et Tke moyenne dans la maille
2722	DO k = 1, klev
2723	DO j = 1, knon
2724	i = ni(j)
2725	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))
2726	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))
2727	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))
2728	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))
2729	ENDDO
2730	ENDDO
2731	DO j=1,knon
2732	yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j))
2733	ENDDO
2734	IF (prt_level >=10) THEN
2735	print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &
2736	nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf)
2737	ENDIF
2738
2739	DO k = 1, klev
2740	DO j = 1, knon
2741	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))
2742	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))
2743	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))
2744	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))
2745	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))
2746	ENDDO
2747	ENDDO
2748
2749	ENDIF ! (iflag_split .eq.0)
2750	!!!
2751
2752	! tendencies of blowing snow
2753	IF (ok_bs) THEN
2754	DO k = 1, klev
2755	DO j = 1, knon
2756	i = ni(j)
2757	y_d_qbs(j,k)=y_d_qbs(j,k) * ypct(j)
2758	flux_qbs(i,k,nsrf) = y_flux_qbs(j,k)
2759	ENDDO
2760	ENDDO
2761	ENDIF
2762
2763
2764	DO j = 1, knon
2765	i = ni(j)
2766	evap(i,nsrf) = - flux_q(i,1,nsrf) !jyg
2767	if (ok_bs) then ; snowerosion(i,nsrf)=flux_qbs(i,1,nsrf); endif
2768	beta(i,nsrf) = ybeta(j) !jyg
2769	d_ts(i,nsrf) = y_d_ts(j)
2770	!albedo SB >>>
2771	DO k=1,nsw
2772	alb_dir(i,k,nsrf) = yalb_dir_new(j,k)
2773	alb_dif(i,k,nsrf) = yalb_dif_new(j,k)
2774	ENDDO
2775	!albedo SB <<<
2776	snow(i,nsrf) = ysnow(j)
2777	qsurf(i,nsrf) = yqsurf(j)
2778	z0m(i,nsrf) = yz0m(j)
2779	z0h(i,nsrf) = yz0h(j)
2780	fluxlat(i,nsrf) = yfluxlat(j)
2781	agesno(i,nsrf) = yagesno(j)
2782	cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
2783	cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
2784	dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)
2785	dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)
2786	ENDDO
2787
2788	! print*,'Dans pbl OK2'
2789
2790	!!! jyg le 07/02/2012
2791	IF (iflag_split .ge.1) THEN
2792	!!!
2793	!!! nrlmd le 02/05/2011
2794	DO j = 1, knon
2795	i = ni(j)
2796	fluxlat_x(i,nsrf) = yfluxlat_x(j)
2797	fluxlat_w(i,nsrf) = yfluxlat_w(j)
2798	!!!
2799	!!! nrlmd le 13/06/2011
2800	!!jyg20170131 delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
2801	!!jyg20210118 delta_tsurf(i,nsrf)=y_delta_tsurf(j)
2802	delta_tsurf(i,nsrf)=y_delta_tsurf_new(j)
2803	!
2804	delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j)*ypct(j)
2805	!
2806	cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
2807	cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)
2808	cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j)
2809	cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j)
2810	kh(i) = kh(i) + Kech_h(j)*ypct(j)
2811	kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j)
2812	kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j)
2813	!!!
2814	ENDDO
2815	!!!
2816	ENDIF ! (iflag_split .ge.1)
2817	!!!
2818	!!! nrlmd le 02/05/2011
2819	!!jyg le 20/02/2011
2820	!! tke_x(:,:,nsrf)=0.
2821	!! tke_w(:,:,nsrf)=0.
2822	!!jyg le 20/02/2011
2823	!! DO k = 1, klev+1
2824	!! DO j = 1, knon
2825	!! i = ni(j)
2826	!! wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2827	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2828	!! ENDDO
2829	!! ENDDO
2830	!!jyg le 20/02/2011
2831	!! DO k = 1, klev+1
2832	!! DO j = 1, knon
2833	!! i = ni(j)
2834	!! tke(i,k,nsrf)=(1.-ywake_s(j))tke_x(i,k,nsrf)+ywake_s(j)tke_w(i,k,nsrf)
2835	!! ENDDO
2836	!! ENDDO
2837	!!!
2838	IF (iflag_split .eq.0) THEN
2839	wake_dltke(:,:,nsrf) = 0.
2840	DO k = 1, klev+1
2841	DO j = 1, knon
2842	i = ni(j)
2843	!jyg<
2844	!! tke(i,k,nsrf) = ytke(j,k)
2845	!! tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
2846	tke_x(i,k,nsrf) = ytke(j,k)
2847	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
2848
2849	!>jyg
2850	ENDDO
2851	ENDDO
2852
2853	ELSE ! (iflag_split .eq.0)
2854	DO k = 1, klev+1
2855	DO j = 1, knon
2856	i = ni(j)
2857	wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2858	!jyg<
2859	!! tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2860	!! tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
2861	tke_x(i,k,nsrf) = ytke_x(j,k)
2862	tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)
2863	wake_dltke(i,k,is_ave) = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
2864
2865
2866	!>jyg
2867	ENDDO
2868	ENDDO
2869	ENDIF ! (iflag_split .eq.0)
2870	!!!
2871	DO k = 2, klev
2872	DO j = 1, knon
2873	i = ni(j)
2874	zcoefh(i,k,nsrf) = ycoefh(j,k)
2875	zcoefm(i,k,nsrf) = ycoefm(j,k)
2876	zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
2877	zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
2878	ENDDO
2879	ENDDO
2880
2881	! print*,'Dans pbl OK3'
2882
2883	IF ( nsrf .EQ. is_ter ) THEN
2884	DO j = 1, knon
2885	i = ni(j)
2886	qsol(i) = yqsol(j)
2887	ENDDO
2888	ENDIF
2889
2890	!jyg<
2891	!! ftsoil(:,:,nsrf) = 0.
2892	!>jyg
2893	DO k = 1, nsoilmx
2894	DO j = 1, knon
2895	i = ni(j)
2896	ftsoil(i, k, nsrf) = ytsoil(j,k)
2897	ENDDO
2898	ENDDO
2899
2900	!!! jyg le 07/02/2012
2901	IF (iflag_split .ge.1) THEN
2902	!!!
2903	!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
2904	DO k = 1, klev
2905	DO j = 1, knon
2906	i = ni(j)
2907	d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)
2908	d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)
2909	d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)
2910	d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)
2911	d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)
2912	!
2913	d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)
2914	d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)
2915	d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)
2916	d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
2917	d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
2918	!
2919	!! d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
2920	!! d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)
2921	ENDDO
2922	ENDDO
2923	!!!
2924	ENDIF ! (iflag_split .ge.1)
2925	!!!
2926
2927	DO k = 1, klev
2928	DO j = 1, knon
2929	i = ni(j)
2930	d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
2931	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
2932	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
2933	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
2934	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
2935	ENDDO
2936	ENDDO
2937
2938
2939	IF (ok_bs) THEN
2940	DO k = 1, klev
2941	DO j = 1, knon
2942	i = ni(j)
2943	d_qbs(i,k) = d_qbs(i,k) + y_d_qbs(j,k)
2944	ENDDO
2945	ENDDO
2946	ENDIF
2947
2948	! print*,'Dans pbl OK4'
2949
2950	IF (prt_level >=10) THEN
2951	print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
2952	d_t_w(:,1), d_t_x(:,1), d_t(:,1)
2953	ENDIF
2954
2955	if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
2956	delta_sal = missing_val
2957	ds_ns = missing_val
2958	dt_ns = missing_val
2959	delta_sst = missing_val
2960	dter = missing_val
2961	dser = missing_val
2962	tkt = missing_val
2963	tks = missing_val
2964	taur = missing_val
2965	sss = missing_val
2966
2967	delta_sal(ni(:knon)) = ydelta_sal(:knon)
2968	ds_ns(ni(:knon)) = yds_ns(:knon)
2969	dt_ns(ni(:knon)) = ydt_ns(:knon)
2970	delta_sst(ni(:knon)) = ydelta_sst(:knon)
2971	dter(ni(:knon)) = ydter(:knon)
2972	dser(ni(:knon)) = ydser(:knon)
2973	tkt(ni(:knon)) = ytkt(:knon)
2974	tks(ni(:knon)) = ytks(:knon)
2975	taur(ni(:knon)) = ytaur(:knon)
2976	sss(ni(:knon)) = ysss(:knon)
2977
2978	if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
2979	dt_ds = missing_val
2980	dt_ds(ni(:knon)) = ydt_ds(:knon)
2981	end if
2982	end if
2983
2984
2985	!****************************************************************************************
2986	! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m
2987	! Call HBTM
2988	!
2989	!****************************************************************************************
2990	!!!
2991	!
2992	#undef T2m
2993	#define T2m
2994	#ifdef T2m
2995	! Calculations of diagnostic t,q at 2m and u, v at 10m
2996
2997	! print*,'Dans pbl OK41'
2998	! print*,'tair1,yt(:,1),y_d_t(:,1)'
2999	! print*, tair1,yt(:,1),y_d_t(:,1)
3000	!!! jyg le 07/02/2012
3001	IF (iflag_split .eq.0) THEN
3002	DO j=1, knon
3003	uzon(j) = yu(j,1) + y_d_u(j,1)
3004	vmer(j) = yv(j,1) + y_d_v(j,1)
3005	tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
3006	qair1(j) = yq(j,1) + y_d_q(j,1)
3007	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
3008	* (ypaprs(j,1)-ypplay(j,1))
3009	tairsol(j) = yts(j) + y_d_ts(j)
3010	qairsol(j) = yqsurf(j)
3011	ENDDO
3012	ELSE ! (iflag_split .eq.0)
3013	DO j=1, knon
3014	uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
3015	vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
3016	tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
3017	qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
3018	zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
3019	* (ypaprs(j,1)-ypplay(j,1))
3020	tairsol(j) = yts(j) + y_d_ts(j)
3021	!! tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
3022	tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf_new(j)
3023	qairsol(j) = yqsurf(j)
3024	ENDDO
3025	DO j=1, knon
3026	uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
3027	vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
3028	tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
3029	qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
3030	zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
3031	* (ypaprs(j,1)-ypplay(j,1))
3032	tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
3033	qairsol(j) = yqsurf(j)
3034	ENDDO
3035	!!!
3036	ENDIF ! (iflag_split .eq.0)
3037	!!!
3038	DO j=1, knon
3039	! i = ni(j)
3040	! yz0h_oupas(j) = yz0m(j)
3041	! IF(nsrf.EQ.is_oce) THEN
3042	! yz0h_oupas(j) = z0m(i,nsrf)
3043	! ENDIF
3044	psfce(j)=ypaprs(j,1)
3045	patm(j)=ypplay(j,1)
3046	ENDDO
3047
3048	IF (iflag_pbl_surface_t2m_bug==1) THEN
3049	yz0h_oupas(1:knon)=yz0m(1:knon)
3050	ELSE
3051	yz0h_oupas(1:knon)=yz0h(1:knon)
3052	ENDIF
3053
3054	! print*,'Dans pbl OK42A'
3055	! print*,'tair1,yt(:,1),y_d_t(:,1)'
3056	! print*, tair1,yt(:,1),y_d_t(:,1)
3057
3058	! Calculate the temperature and relative humidity at 2m and the wind at 10m
3059	!!! jyg le 07/02/2012
3060	IF (iflag_split .eq.0) THEN
3061	IF (iflag_new_t2mq2m==1) THEN
3062	CALL stdlevvarn(klon, knon, nsrf, zxli, &
3063	uzon, vmer, tair1, qair1, zgeo1, &
3064	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3065	yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
3066	yn2mout(:, nsrf, :))
3067	ELSE
3068	CALL stdlevvar(klon, knon, nsrf, zxli, &
3069	uzon, vmer, tair1, qair1, zgeo1, &
3070	tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3071	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
3072	ENDIF
3073	ELSE !(iflag_split .eq.0)
3074	IF (iflag_new_t2mq2m==1) THEN
3075	CALL stdlevvarn(klon, knon, nsrf, zxli, &
3076	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
3077	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3078	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, &
3079	yn2mout_x(:, nsrf, :))
3080	CALL stdlevvarn(klon, knon, nsrf, zxli, &
3081	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
3082	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3083	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, &
3084	yn2mout_w(:, nsrf, :))
3085	ELSE
3086	CALL stdlevvar(klon, knon, nsrf, zxli, &
3087	uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
3088	tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3089	yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
3090	CALL stdlevvar(klon, knon, nsrf, zxli, &
3091	uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
3092	tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3093	yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
3094	ENDIF
3095	!!!
3096	ENDIF ! (iflag_split .eq.0)
3097	!!!
3098	!!! jyg le 07/02/2012
3099	IF (iflag_split .eq.0) THEN
3100	DO j=1, knon
3101	i = ni(j)
3102	t2m(i,nsrf)=yt2m(j)
3103	q2m(i,nsrf)=yq2m(j)
3104	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3105	ustar(i,nsrf)=yustar(j)
3106	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
3107	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
3108	!
3109	DO k = 1, 6
3110	n2mout(i,nsrf,k) = yn2mout(j,nsrf,k)
3111	END DO
3112	!
3113	ENDDO
3114	ELSE !(iflag_split .eq.0)
3115	DO j=1, knon
3116	i = ni(j)
3117	t2m_x(i,nsrf)=yt2m_x(j)
3118	q2m_x(i,nsrf)=yq2m_x(j)
3119	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3120	ustar_x(i,nsrf)=yustar_x(j)
3121	u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
3122	v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)2+vmer_x(j)2)
3123	!
3124	DO k = 1, 6
3125	n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k)
3126	END DO
3127	!
3128	ENDDO
3129	DO j=1, knon
3130	i = ni(j)
3131	t2m_w(i,nsrf)=yt2m_w(j)
3132	q2m_w(i,nsrf)=yq2m_w(j)
3133	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3134	ustar_w(i,nsrf)=yustar_w(j)
3135	u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
3136	v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)2+vmer_w(j)2)
3137	!
3138	ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
3139	u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
3140	v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
3141	!
3142	DO k = 1, 6
3143	n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k)
3144	END DO
3145	!
3146	ENDDO
3147	!!!
3148	ENDIF ! (iflag_split .eq.0)
3149	!!!
3150
3151	! print*,'Dans pbl OK43'
3152	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
3153	!IM Ajoute dependance type surface
3154	IF (thermcep) THEN
3155	!!! jyg le 07/02/2012
3156	IF (iflag_split .eq.0) THEN
3157	DO j = 1, knon
3158	i=ni(j)
3159	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
3160	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
3161	zx_qs1 = MIN(0.5,zx_qs1)
3162	zcor1 = 1./(1.-RETV*zx_qs1)
3163	zx_qs1 = zx_qs1*zcor1
3164
3165	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
3166	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
3167	ENDDO
3168	ELSE ! (iflag_split .eq.0)
3169	DO j = 1, knon
3170	i=ni(j)
3171	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
3172	zx_qs1 = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)
3173	zx_qs1 = MIN(0.5,zx_qs1)
3174	zcor1 = 1./(1.-RETV*zx_qs1)
3175	zx_qs1 = zx_qs1*zcor1
3176
3177	rh2m_x(i) = rh2m_x(i) + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
3178	qsat2m_x(i) = qsat2m_x(i) + zx_qs1 * pctsrf(i,nsrf)
3179	ENDDO
3180	DO j = 1, knon
3181	i=ni(j)
3182	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
3183	zx_qs1 = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
3184	zx_qs1 = MIN(0.5,zx_qs1)
3185	zcor1 = 1./(1.-RETV*zx_qs1)
3186	zx_qs1 = zx_qs1*zcor1
3187
3188	rh2m_w(i) = rh2m_w(i) + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
3189	qsat2m_w(i) = qsat2m_w(i) + zx_qs1 * pctsrf(i,nsrf)
3190	ENDDO
3191	!!!
3192	ENDIF ! (iflag_split .eq.0)
3193	!!!
3194	ENDIF
3195	!
3196	IF (prt_level >=10) THEN
3197	print *, 'T2m, q2m, RH2m ', &
3198	t2m, q2m, rh2m
3199	ENDIF
3200
3201	! print*,'OK pbl 5'
3202	!
3203	!!! jyg le 07/02/2012
3204	IF (iflag_split .eq.0) THEN
3205	CALL hbtm(knon, ypaprs, ypplay, &
3206	yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
3207	y_flux_t,y_flux_q,yu,yv,yt,yq, &
3208	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
3209	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
3210	IF (prt_level >=10) THEN
3211	print *,' Arg. de HBTM: yt2m ',yt2m
3212	print *,' Arg. de HBTM: yt10m ',yt10m
3213	print *,' Arg. de HBTM: yq2m ',yq2m
3214	print *,' Arg. de HBTM: yq10m ',yq10m
3215	print *,' Arg. de HBTM: yustar ',yustar
3216	print *,' Arg. de HBTM: y_flux_t ',y_flux_t
3217	print *,' Arg. de HBTM: y_flux_q ',y_flux_q
3218	print *,' Arg. de HBTM: yu ',yu
3219	print *,' Arg. de HBTM: yv ',yv
3220	print *,' Arg. de HBTM: yt ',yt
3221	print *,' Arg. de HBTM: yq ',yq
3222	ENDIF
3223	ELSE ! (iflag_split .eq.0)
3224	CALL HBTM(knon, ypaprs, ypplay, &
3225	yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
3226	y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
3227	ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
3228	ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
3229	IF (prt_level >=10) THEN
3230	print *,' Arg. de HBTM: yt2m_x ',yt2m_x
3231	print *,' Arg. de HBTM: yt10m_x ',yt10m_x
3232	print *,' Arg. de HBTM: yq2m_x ',yq2m_x
3233	print *,' Arg. de HBTM: yq10m_x ',yq10m_x
3234	print *,' Arg. de HBTM: yustar_x ',yustar_x
3235	print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x
3236	print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x
3237	print *,' Arg. de HBTM: yu_x ',yu_x
3238	print *,' Arg. de HBTM: yv_x ',yv_x
3239	print *,' Arg. de HBTM: yt_x ',yt_x
3240	print *,' Arg. de HBTM: yq_x ',yq_x
3241	ENDIF
3242	CALL HBTM(knon, ypaprs, ypplay, &
3243	yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
3244	y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
3245	ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
3246	ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
3247	!!!
3248	ENDIF ! (iflag_split .eq.0)
3249	!!!
3250
3251	!!! jyg le 07/02/2012
3252	IF (iflag_split .eq.0) THEN
3253	!!!
3254	DO j=1, knon
3255	i = ni(j)
3256	pblh(i,nsrf) = ypblh(j)
3257	wstar(i,nsrf) = ywstar(j)
3258	plcl(i,nsrf) = ylcl(j)
3259	capCL(i,nsrf) = ycapCL(j)
3260	oliqCL(i,nsrf) = yoliqCL(j)
3261	cteiCL(i,nsrf) = ycteiCL(j)
3262	pblT(i,nsrf) = ypblT(j)
3263	therm(i,nsrf) = ytherm(j)
3264	trmb1(i,nsrf) = ytrmb1(j)
3265	trmb2(i,nsrf) = ytrmb2(j)
3266	trmb3(i,nsrf) = ytrmb3(j)
3267	ENDDO
3268	IF (prt_level >=10) THEN
3269	print *, 'After HBTM: pblh ', pblh
3270	print *, 'After HBTM: plcl ', plcl
3271	print *, 'After HBTM: cteiCL ', cteiCL
3272	ENDIF
3273	ELSE !(iflag_split .eq.0)
3274	DO j=1, knon
3275	i = ni(j)
3276	pblh_x(i,nsrf) = ypblh_x(j)
3277	wstar_x(i,nsrf) = ywstar_x(j)
3278	plcl_x(i,nsrf) = ylcl_x(j)
3279	capCL_x(i,nsrf) = ycapCL_x(j)
3280	oliqCL_x(i,nsrf) = yoliqCL_x(j)
3281	cteiCL_x(i,nsrf) = ycteiCL_x(j)
3282	pblT_x(i,nsrf) = ypblT_x(j)
3283	therm_x(i,nsrf) = ytherm_x(j)
3284	trmb1_x(i,nsrf) = ytrmb1_x(j)
3285	trmb2_x(i,nsrf) = ytrmb2_x(j)
3286	trmb3_x(i,nsrf) = ytrmb3_x(j)
3287	ENDDO
3288	IF (prt_level >=10) THEN
3289	print *, 'After HBTM: pblh_x ', pblh_x
3290	print *, 'After HBTM: plcl_x ', plcl_x
3291	print *, 'After HBTM: cteiCL_x ', cteiCL_x
3292	ENDIF
3293	DO j=1, knon
3294	i = ni(j)
3295	pblh_w(i,nsrf) = ypblh_w(j)
3296	wstar_w(i,nsrf) = ywstar_w(j)
3297	plcl_w(i,nsrf) = ylcl_w(j)
3298	capCL_w(i,nsrf) = ycapCL_w(j)
3299	oliqCL_w(i,nsrf) = yoliqCL_w(j)
3300	cteiCL_w(i,nsrf) = ycteiCL_w(j)
3301	pblT_w(i,nsrf) = ypblT_w(j)
3302	therm_w(i,nsrf) = ytherm_w(j)
3303	trmb1_w(i,nsrf) = ytrmb1_w(j)
3304	trmb2_w(i,nsrf) = ytrmb2_w(j)
3305	trmb3_w(i,nsrf) = ytrmb3_w(j)
3306	ENDDO
3307	IF (prt_level >=10) THEN
3308	print *, 'After HBTM: pblh_w ', pblh_w
3309	print *, 'After HBTM: plcl_w ', plcl_w
3310	print *, 'After HBTM: cteiCL_w ', cteiCL_w
3311	ENDIF
3312	!!!
3313	ENDIF ! (iflag_split .eq.0)
3314	!!!
3315
3316	! print*,'OK pbl 6'
3317	#else
3318	! T2m not defined
3319	! No calculation
3320	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
3321	#endif
3322
3323	!****************************************************************************************
3324	! 15) End of loop over different surfaces
3325	!
3326	!****************************************************************************************
3327	ENDDO loop_nbsrf
3328	!
3329	!----------------------------------------------------------------------------------------
3330	! Reset iflag_split
3331	!
3332	iflag_split=iflag_split_ref
3333
3334	!****************************************************************************************
3335	! 16) Calculate the mean value over all sub-surfaces for some variables
3336	!
3337	!****************************************************************************************
3338
3339	z0m(:,nbsrf+1) = 0.0
3340	z0h(:,nbsrf+1) = 0.0
3341	DO nsrf = 1, nbsrf
3342	DO i = 1, klon
3343	z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)
3344	z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)
3345	ENDDO
3346	ENDDO
3347
3348	! print*,'OK pbl 7'
3349	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
3350	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
3351	zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
3352	zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
3353	zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
3354	zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
3355
3356	!!! jyg le 07/02/2012
3357	IF (iflag_split .ge.1) THEN
3358	!!!
3359	!!! nrlmd & jyg les 02/05/2011, 05/02/2012
3360
3361	DO nsrf = 1, nbsrf
3362	DO k = 1, klev
3363	DO i = 1, klon
3364	zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
3365	zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
3366	zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
3367	zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
3368	!
3369	zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
3370	zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
3371	zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
3372	zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
3373	ENDDO
3374	ENDDO
3375	ENDDO
3376
3377	DO i = 1, klon
3378	zxsens_x(i) = - zxfluxt_x(i,1)
3379	zxsens_w(i) = - zxfluxt_w(i,1)
3380	ENDDO
3381	!!!
3382	ENDIF ! (iflag_split .ge.1)
3383	!!!
3384
3385	DO nsrf = 1, nbsrf
3386	DO k = 1, klev
3387	DO i = 1, klon
3388	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
3389	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
3390	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
3391	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
3392	ENDDO
3393	ENDDO
3394	ENDDO
3395
3396	DO i = 1, klon
3397	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
3398	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
3399	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
3400	ENDDO
3401
3402	! if blowing snow
3403	if (ok_bs) then
3404	DO nsrf = 1, nbsrf
3405	DO k = 1, klev
3406	DO i = 1, klon
3407	zxfluxqbs(i,k) = zxfluxqbs(i,k) + flux_qbs(i,k,nsrf) * pctsrf(i,nsrf)
3408	ENDDO
3409	ENDDO
3410	ENDDO
3411
3412	DO i = 1, klon
3413	zxsnowerosion(i) = zxfluxqbs(i,1) ! blowings snow flux at the surface
3414	END DO
3415	endif
3416
3417	!!!
3418
3419	!
3420	! Incrementer la temperature du sol
3421	!
3422	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
3423	zt2m(:) = 0.0 ; zq2m(:) = 0.0 ; zn2mout(:,:) = 0
3424	zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
3425	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
3426	!!! jyg le 07/02/2012
3427	s_pblh_x(:) = 0.0 ; s_plcl_x(:) = 0.0
3428	s_pblh_w(:) = 0.0 ; s_plcl_w(:) = 0.0
3429	!!!
3430	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
3431	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
3432	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
3433	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
3434	wstar(:,is_ave)=0.
3435
3436	! print*,'OK pbl 9'
3437
3438	!!! nrlmd le 02/05/2011
3439	zxfluxlat_x(:) = 0.0 ; zxfluxlat_w(:) = 0.0
3440	!!!
3441
3442	DO nsrf = 1, nbsrf
3443	DO i = 1, klon
3444	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
3445
3446	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
3447	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
3448	wfbilo(i,nsrf) = (evap(i,nsrf)-(rain_f(i)+snow_f(i)))*pctsrf(i,nsrf)
3449	wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)
3450	wfrain(i,nsrf) = rain_f(i)*pctsrf(i,nsrf)
3451	wfsnow(i,nsrf) = snow_f(i)*pctsrf(i,nsrf)
3452
3453	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
3454	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
3455	ENDDO
3456	ENDDO
3457	!
3458	!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
3459	IF (iflag_order2_sollw == 1) THEN
3460	meansqT(:) = 0. ! as working buffer
3461	DO nsrf = 1, nbsrf
3462	DO i = 1, klon
3463	meansqT(i) = meansqT(i)+(ts(i,nsrf)-zxtsol(i))*2 pctsrf(i,nsrf)
3464	ENDDO
3465	ENDDO
3466	zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)
3467	ENDIF ! iflag_order2_sollw == 1
3468	!>al1
3469
3470	!!! jyg le 07/02/2012
3471	IF (iflag_split .eq.0) THEN
3472	DO nsrf = 1, nbsrf
3473	DO i = 1, klon
3474	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
3475	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
3476	!
3477	DO k = 1, 6
3478	zn2mout(i,k) = zn2mout(i,k) + n2mout(i,nsrf,k) * pctsrf(i,nsrf)
3479	ENDDO
3480	!
3481	zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
3482	wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
3483	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
3484	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
3485
3486	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
3487	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
3488	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
3489	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
3490	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
3491	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
3492	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
3493	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
3494	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
3495	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
3496	ENDDO
3497	ENDDO
3498	ELSE !(iflag_split .eq.0)
3499	DO nsrf = 1, nbsrf
3500	DO i = 1, klon
3501	!!! nrlmd le 02/05/2011
3502	zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
3503	zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
3504	!!!
3505	!!! jyg le 08/02/2012
3506	!! Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;
3507	!! pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
3508	!! pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
3509	!! pour les autres variables, on sort les valeurs de la region (x).
3510	zt2m(i) = zt2m(i) + (t2m_x(i,nsrf)+wake_s(i)(t2m_w(i,nsrf)-t2m_x(i,nsrf))) pctsrf(i,nsrf)
3511	zq2m(i) = zq2m(i) + q2m_x(i,nsrf) * pctsrf(i,nsrf)
3512	!
3513	DO k = 1, 6
3514	zn2mout(i,k) = zn2mout(i,k) + n2mout_x(i,nsrf,k) * pctsrf(i,nsrf)
3515	ENDDO
3516	!
3517	zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
3518	wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
3519	zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
3520	zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
3521	!
3522	s_pblh(i) = s_pblh(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
3523	s_pblh_x(i) = s_pblh_x(i) + pblh_x(i,nsrf) * pctsrf(i,nsrf)
3524	s_pblh_w(i) = s_pblh_w(i) + pblh_w(i,nsrf) * pctsrf(i,nsrf)
3525	!
3526	s_plcl(i) = s_plcl(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
3527	s_plcl_x(i) = s_plcl_x(i) + plcl_x(i,nsrf) * pctsrf(i,nsrf)
3528	s_plcl_w(i) = s_plcl_w(i) + plcl_w(i,nsrf) * pctsrf(i,nsrf)
3529	!
3530	s_capCL(i) = s_capCL(i) + capCL_x(i,nsrf) * pctsrf(i,nsrf)
3531	s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
3532	s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
3533	s_pblT(i) = s_pblT(i) + pblT_x(i,nsrf) * pctsrf(i,nsrf)
3534	s_therm(i) = s_therm(i) + therm_x(i,nsrf) * pctsrf(i,nsrf)
3535	s_trmb1(i) = s_trmb1(i) + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
3536	s_trmb2(i) = s_trmb2(i) + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
3537	s_trmb3(i) = s_trmb3(i) + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
3538	ENDDO
3539	ENDDO
3540	DO i = 1, klon
3541	qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
3542	ENDDO
3543	!!!
3544	ENDIF ! (iflag_split .eq.0)
3545	!!!
3546
3547	IF (check) THEN
3548	amn=MIN(ts(1,is_ter),1000.)
3549	amx=MAX(ts(1,is_ter),-1000.)
3550	DO i=2, klon
3551	amn=MIN(ts(i,is_ter),amn)
3552	amx=MAX(ts(i,is_ter),amx)
3553	ENDDO
3554	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
3555	ENDIF
3556
3557	!jg ?
3558	!!$!
3559	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
3560	!!$! sub-surfaces is distributed.
3561	!!$!
3562	!!$ DO nsrf = 1, nbsrf
3563	!!$ DO i = 1, klon
3564	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
3565	!!$ ts(i,nsrf) = zxtsol(i)
3566	!!$ t2m(i,nsrf) = zt2m(i)
3567	!!$ q2m(i,nsrf) = zq2m(i)
3568	!!$ u10m(i,nsrf) = zu10m(i)
3569	!!$ v10m(i,nsrf) = zv10m(i)
3570	!!$
3571	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
3572	!!$ pblh(i,nsrf) = s_pblh(i)
3573	!!$ plcl(i,nsrf) = s_plcl(i)
3574	!!$ capCL(i,nsrf) = s_capCL(i)
3575	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
3576	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
3577	!!$ pblT(i,nsrf) = s_pblT(i)
3578	!!$ therm(i,nsrf) = s_therm(i)
3579	!!$ trmb1(i,nsrf) = s_trmb1(i)
3580	!!$ trmb2(i,nsrf) = s_trmb2(i)
3581	!!$ trmb3(i,nsrf) = s_trmb3(i)
3582	!!$ ENDIF
3583	!!$ ENDDO
3584	!!$ ENDDO
3585
3586
3587	DO i = 1, klon
3588	fder(i) = - 4.0RSIGMAzxtsol(i)**3
3589	ENDDO
3590
3591	zxqsurf(:) = 0.0
3592	zxsnow(:) = 0.0
3593	DO nsrf = 1, nbsrf
3594	DO i = 1, klon
3595	zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
3596	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
3597	ENDDO
3598	ENDDO
3599
3600	! Premier niveau de vent sortie dans physiq.F
3601	zu1(:) = u(:,1)
3602	zv1(:) = v(:,1)
3603
3604	END SUBROUTINE pbl_surface
3605	!
3606	!****************************************************************************************
3607	!
3608	SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
3609
3610	USE indice_sol_mod
3611
3612	INCLUDE "dimsoil.h"
3613
3614	! Ouput variables
3615	!****************************************************************************************
3616	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
3617	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
3618	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
3619	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
3620
3621
3622	!****************************************************************************************
3623	! Return module variables for writing to restart file
3624	!
3625	!****************************************************************************************
3626	fder_rst(:) = fder(:)
3627	snow_rst(:,:) = snow(:,:)
3628	qsurf_rst(:,:) = qsurf(:,:)
3629	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
3630
3631	!****************************************************************************************
3632	! Deallocate module variables
3633	!
3634	!****************************************************************************************
3635	! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
3636	IF (ALLOCATED(fder)) DEALLOCATE(fder)
3637	IF (ALLOCATED(snow)) DEALLOCATE(snow)
3638	IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
3639	IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
3640	IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
3641	IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
3642
3643	!jyg<
3644	!****************************************************************************************
3645	! Deallocate variables for pbl splitting
3646	!
3647	!****************************************************************************************
3648
3649	CALL wx_pbl_final
3650	!>jyg
3651
3652	END SUBROUTINE pbl_surface_final
3653	!
3654	!****************************************************************************************
3655	!
3656
3657	!albedo SB >>>
3658	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
3659	evap, z0m, z0h, agesno, &
3660	tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)
3661	!albedo SB <<<
3662	! Give default values where new fraction has appread
3663
3664	USE indice_sol_mod
3665	use phys_state_var_mod, only: delta_sal, ds_ns, dt_ns, delta_sst, dter, &
3666	dser, dt_ds
3667	use config_ocean_skin_m, only: activate_ocean_skin
3668
3669	INCLUDE "dimsoil.h"
3670	INCLUDE "clesphys.h"
3671	INCLUDE "compbl.h"
3672
3673	! Input variables
3674	!****************************************************************************************
3675	INTEGER, INTENT(IN) :: itime
3676	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
3677
3678	! InOutput variables
3679	!****************************************************************************************
3680	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
3681	!albedo SB >>>
3682	REAL, DIMENSION(klon,nsw,nbsrf), INTENT(INOUT) :: alb_dir, alb_dif
3683	INTEGER :: k
3684	!albedo SB <<<
3685	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m
3686	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: evap, agesno
3687	REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT) :: z0m,z0h
3688	REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
3689
3690	! Local variables
3691	!****************************************************************************************
3692	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
3693	CHARACTER(len=80) :: abort_message
3694	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
3695	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
3696	!
3697	! All at once !!
3698	!****************************************************************************************
3699
3700	DO nsrf = 1, nbsrf
3701	! First decide complement sub-surfaces
3702	SELECT CASE (nsrf)
3703	CASE(is_oce)
3704	nsrf_comp1=is_sic
3705	nsrf_comp2=is_ter
3706	nsrf_comp3=is_lic
3707	CASE(is_sic)
3708	nsrf_comp1=is_oce
3709	nsrf_comp2=is_ter
3710	nsrf_comp3=is_lic
3711	CASE(is_ter)
3712	nsrf_comp1=is_lic
3713	nsrf_comp2=is_oce
3714	nsrf_comp3=is_sic
3715	CASE(is_lic)
3716	nsrf_comp1=is_ter
3717	nsrf_comp2=is_oce
3718	nsrf_comp3=is_sic
3719	END SELECT
3720
3721	! Initialize all new fractions
3722	DO i=1, klon
3723	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
3724
3725	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
3726	! Use the complement sub-surface, keeping the continents unchanged
3727	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
3728	evap(i,nsrf) = evap(i,nsrf_comp1)
3729	z0m(i,nsrf) = z0m(i,nsrf_comp1)
3730	z0h(i,nsrf) = z0h(i,nsrf_comp1)
3731	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
3732	!albedo SB >>>
3733	DO k=1,nsw
3734	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp1)
3735	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp1)
3736	ENDDO
3737	!albedo SB <<<
3738	ustar(i,nsrf) = ustar(i,nsrf_comp1)
3739	u10m(i,nsrf) = u10m(i,nsrf_comp1)
3740	v10m(i,nsrf) = v10m(i,nsrf_comp1)
3741	IF (iflag_pbl > 1) THEN
3742	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
3743	ENDIF
3744	mfois(nsrf) = mfois(nsrf) + 1
3745	! F. Codron sensible default values for ocean and sea ice
3746	IF (nsrf.EQ.is_oce) THEN
3747	tsurf(i,nsrf) = 271.35
3748	! (temperature of sea water under sea ice, so that
3749	! is also the temperature of appearing sea water)
3750	DO k=1,nsw
3751	alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
3752	alb_dif(i,k,nsrf) = 0.06
3753	ENDDO
3754	if (activate_ocean_skin >= 1) then
3755	if (activate_ocean_skin == 2 &
3756	.and. type_ocean == "couple") then
3757	delta_sal(i) = 0.
3758	delta_sst(i) = 0.
3759	dter(i) = 0.
3760	dser(i) = 0.
3761	dt_ds(i) = 0.
3762	end if
3763
3764	ds_ns(i) = 0.
3765	dt_ns(i) = 0.
3766	end if
3767	ELSE IF (nsrf.EQ.is_sic) THEN
3768	tsurf(i,nsrf) = 271.35
3769	! (Temperature at base of sea ice. Surface
3770	! temperature could be higher, up to 0 Celsius
3771	! degrees. We set it to -1.8 Celsius degrees for
3772	! consistency with the ocean slab model.)
3773	DO k=1,nsw
3774	alb_dir(i,k,nsrf) = 0.3 ! thin ice
3775	alb_dif(i,k,nsrf) = 0.3
3776	ENDDO
3777	ENDIF
3778	ELSE
3779	! The continents have changed. The new fraction receives the mean sum of the existent fractions
3780	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3781	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3782	z0m(i,nsrf) = z0m(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0m(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3783	z0h(i,nsrf) = z0h(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + z0h(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3784	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3785	!albedo SB >>>
3786	DO k=1,nsw
3787	alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3788	alb_dir(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3789	alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3790	alb_dif(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3791	ENDDO
3792	!albedo SB <<<
3793	ustar(i,nsrf) = ustar(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3794	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3795	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
3796	IF (iflag_pbl > 1) THEN
3797	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
3798	ENDIF
3799
3800	! Security abort. This option has never been tested. To test, comment the following line.
3801	! abort_message='The fraction of the continents have changed!'
3802	! CALL abort_physic(modname,abort_message,1)
3803	nfois(nsrf) = nfois(nsrf) + 1
3804	ENDIF
3805	snow(i,nsrf) = 0.
3806	agesno(i,nsrf) = 0.
3807	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
3808	ELSE
3809	pfois(nsrf) = pfois(nsrf)+ 1
3810	ENDIF
3811	ENDDO
3812
3813	ENDDO
3814
3815	END SUBROUTINE pbl_surface_newfrac
3816	!
3817	!****************************************************************************************
3818	!
3819	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format