Context Navigation

pbl_surface_mod.F90 @ 2243

Last change on this file since 2243 was 2243, checked in by fhourdin, 9 years ago

Revisite de la formule des flux de surface
(en priorité sur l'océan) en tenant compte des bourrasques de
vent et de la différence entre les hauteurs de rugosités pour
la quantité de mouvement, l'enthalpie et éventuellement l'humidité.

Etape 2 :

Séparation des z0 pour la quantité de mouvement et l'enthalpie.

rugs (ou frugs, rugos, yrugos ...) disparait au profit de z0m, z0h.
Les variables qui étaient à la fois dans pbl_surface_init et

dans l'interface de pbl_surface sont suprimées de pbl_surface_init.

On travaille directement pour ces variables (evap, z0, qsol, agesno)
avec les versions de phys_state_var_mod (qui étaient
précédemment dans phys_local_var_mod

Nouveaux paramètres de contrôle :
- iflag_z0_oce (par défaut 0, et seule option active jusque là)
- z0m_seaice_omp, z0h_seaice_omp, comme leur nom l'indique (utilisées dans surf_landice
- z0min appliqué sur z0m et z0h dans pbl_surface
Introduction des fonction phyeta0_get et phyetat0_srf pour lire

les conditions de initiales dans startphy.
Du coup une seule ligne suffit pour lire et contrôler d'éventuels
problèmes.

Pour la variable fxrugs, elle est remplacée par z0m(:,nbsrf+1)

Ce choix déjà utilisé pour d'autres variables pourrait être
systématiser pour alléger l'interface de pbl_surface_mod.

Dans les sorties, les variables rugs* ont été remplacées par

des z0m* et z0h*

Nettoyage des anciens alb1/alb2 dans les lectures/écritures

des états de redémarrage (et dans pbl_surface_mod.F90).

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: 123.1 KB

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

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

Download in other formats:

Original Format