Context Navigation

pbl_surface_mod.F90 @ 2179

Last change on this file since 2179 was 2177, checked in by Laurent Fairhead, 9 years ago
Modification pour retrouver l'exécution en OpenMP
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:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 120.6 KB

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

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

Download in other formats:

Original Format