Context Navigation

pbl_surface_mod.F90 @ 1267

Last change on this file since 1267 was 1235, checked in by musat, 15 years ago
Petit bug tke: evolution tke uniquement si iflag_pbl>1 IM
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 58.1 KB

Line
1	!
2	MODULE pbl_surface_mod
3	!
4	! Planetary Boundary Layer and Surface module
5	!
6	! This module manage the calculation of turbulent diffusion in the boundary layer
7	! and all interactions towards the differents sub-surfaces.
8	!
9	!
10	USE dimphy
11	USE mod_phys_lmdz_para, ONLY : mpi_size
12	USE ioipsl
13	USE surface_data, ONLY : type_ocean, ok_veget
14	USE surf_land_mod, ONLY : surf_land
15	USE surf_landice_mod, ONLY : surf_landice
16	USE surf_ocean_mod, ONLY : surf_ocean
17	USE surf_seaice_mod, ONLY : surf_seaice
18	USE cpl_mod, ONLY : gath2cpl
19	USE climb_hq_mod, ONLY : climb_hq_down, climb_hq_up
20	USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up
21	USE coef_diff_turb_mod, ONLY : coef_diff_turb
22
23	IMPLICIT NONE
24
25	! Declaration of variables saved in restart file
26	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: qsol ! water height in the soil (mm)
27	!$OMP THREADPRIVATE(qsol)
28	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: fder ! flux drift
29	!$OMP THREADPRIVATE(fder)
30	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: snow ! snow at surface
31	!$OMP THREADPRIVATE(snow)
32	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qsurf ! humidity at surface
33	!$OMP THREADPRIVATE(qsurf)
34	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: evap ! evaporation at surface
35	!$OMP THREADPRIVATE(evap)
36	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: rugos ! rugosity at surface (m)
37	!$OMP THREADPRIVATE(rugos)
38	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: agesno ! age of snow at surface
39	!$OMP THREADPRIVATE(agesno)
40	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: ftsoil ! soil temperature
41	!$OMP THREADPRIVATE(ftsoil)
42
43	CONTAINS
44	!
45	!****************************************************************************************
46	!
47	SUBROUTINE pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst,&
48	evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
49
50	! This routine should be called after the restart file has been read.
51	! This routine initialize the restart variables and does some validation tests
52	! for the index of the different surfaces and tests the choice of type of ocean.
53
54	INCLUDE "indicesol.h"
55	INCLUDE "dimsoil.h"
56	INCLUDE "iniprint.h"
57
58	! Input variables
59	!****************************************************************************************
60	REAL, DIMENSION(klon), INTENT(IN) :: qsol_rst
61	REAL, DIMENSION(klon), INTENT(IN) :: fder_rst
62	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: snow_rst
63	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: qsurf_rst
64	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: evap_rst
65	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: rugos_rst
66	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: agesno_rst
67	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
68
69
70	! Local variables
71	!****************************************************************************************
72	INTEGER :: ierr
73	CHARACTER(len=80) :: abort_message
74	CHARACTER(len = 20) :: modname = 'pbl_surface_init'
75
76
77	!****************************************************************************************
78	! Allocate and initialize module variables with fields read from restart file.
79	!
80	!****************************************************************************************
81	ALLOCATE(qsol(klon), stat=ierr)
82	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
83
84	ALLOCATE(fder(klon), stat=ierr)
85	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
86
87	ALLOCATE(snow(klon,nbsrf), stat=ierr)
88	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
89
90	ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
91	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
92
93	ALLOCATE(evap(klon,nbsrf), stat=ierr)
94	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
95
96	ALLOCATE(rugos(klon,nbsrf), stat=ierr)
97	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
98
99	ALLOCATE(agesno(klon,nbsrf), stat=ierr)
100	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
101
102	ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
103	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
104
105
106	qsol(:) = qsol_rst(:)
107	fder(:) = fder_rst(:)
108	snow(:,:) = snow_rst(:,:)
109	qsurf(:,:) = qsurf_rst(:,:)
110	evap(:,:) = evap_rst(:,:)
111	rugos(:,:) = rugos_rst(:,:)
112	agesno(:,:) = agesno_rst(:,:)
113	ftsoil(:,:,:) = ftsoil_rst(:,:,:)
114
115
116	!****************************************************************************************
117	! Test for sub-surface indices
118	!
119	!****************************************************************************************
120	IF (is_ter /= 1) THEN
121	WRITE(lunout,)" Warning *"
122	WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
123	WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
124	abort_message="voir ci-dessus"
125	CALL abort_gcm(modname,abort_message,1)
126	ENDIF
127
128	IF ( is_oce > is_sic ) THEN
129	WRITE(lunout,)' Warning *'
130	WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
131	WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
132	WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
133	abort_message='voir ci-dessus'
134	CALL abort_gcm(modname,abort_message,1)
135	ENDIF
136
137	IF ( is_lic > is_sic ) THEN
138	WRITE(lunout,)' Warning *'
139	WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
140	WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
141	WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
142	abort_message='voir ci-dessus'
143	CALL abort_gcm(modname,abort_message,1)
144	ENDIF
145
146	!****************************************************************************************
147	! Validation of ocean mode
148	!
149	!****************************************************************************************
150
151	IF (type_ocean /= 'slab ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN
152	WRITE(lunout,)' Warning *'
153	WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean
154	abort_message='option pour l''ocean non valable'
155	CALL abort_gcm(modname,abort_message,1)
156	ENDIF
157
158	END SUBROUTINE pbl_surface_init
159	!
160	!****************************************************************************************
161	!
162
163	SUBROUTINE pbl_surface( &
164	dtime, date0, itap, jour, &
165	debut, lafin, &
166	rlon, rlat, rugoro, rmu0, &
167	rain_f, snow_f, solsw_m, sollw_m, &
168	t, q, u, v, &
169	pplay, paprs, pctsrf, &
170	ts, alb1, alb2, u10m, v10m, &
171	lwdown_m, cdragh, cdragm, zu1, zv1, &
172	alb1_m, alb2_m, zxsens, zxevap, &
173	zxtsol, zxfluxlat, zt2m, qsat2m, &
174	d_t, d_q, d_u, d_v, &
175	zcoefh, slab_wfbils, &
176	qsol_d, zq2m, s_pblh, s_plcl, &
177	s_capCL, s_oliqCL, s_cteiCL, s_pblT, &
178	s_therm, s_trmb1, s_trmb2, s_trmb3, &
179	zxrugs, zu10m, zv10m, fder_print, &
180	zxqsurf, rh2m, zxfluxu, zxfluxv, &
181	rugos_d, agesno_d, sollw, solsw, &
182	d_ts, evap_d, fluxlat, t2m, &
183	wfbils, wfbilo, flux_t, flux_u, flux_v,&
184	dflux_t, dflux_q, zxsnow, &
185	zxfluxt, zxfluxq, q2m, flux_q, tke )
186	!****************************************************************************************
187	! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
188	! Objet: interface de "couche limite" (diffusion verticale)
189	!
190	!AA REM:
191	!AA-----
192	!AA Tout ce qui a trait au traceurs est dans phytrac maintenant
193	!AA pour l'instant le calcul de la couche limite pour les traceurs
194	!AA se fait avec cltrac et ne tient pas compte de la differentiation
195	!AA des sous-fraction de sol.
196	!AA REM bis :
197	!AA----------
198	!AA Pour pouvoir extraire les coefficient d'echanges et le vent
199	!AA dans la premiere couche, 3 champs supplementaires ont ete crees
200	!AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
201	!AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
202	!AA si les informations des subsurfaces doivent etre prises en compte
203	!AA il faudra sortir ces memes champs en leur ajoutant une dimension,
204	!AA c'est a dire nbsrf (nbre de subsurface).
205	!
206	! Arguments:
207	!
208	! dtime----input-R- interval du temps (secondes)
209	! itap-----input-I- numero du pas de temps
210	! date0----input-R- jour initial
211	! t--------input-R- temperature (K)
212	! q--------input-R- vapeur d'eau (kg/kg)
213	! u--------input-R- vitesse u
214	! v--------input-R- vitesse v
215	! ts-------input-R- temperature du sol (en Kelvin)
216	! paprs----input-R- pression a intercouche (Pa)
217	! pplay----input-R- pression au milieu de couche (Pa)
218	! rlat-----input-R- latitude en degree
219	! rugos----input-R- longeur de rugosite (en m)
220	!
221	! d_t------output-R- le changement pour "t"
222	! d_q------output-R- le changement pour "q"
223	! d_u------output-R- le changement pour "u"
224	! d_v------output-R- le changement pour "v"
225	! d_ts-----output-R- le changement pour "ts"
226	! flux_t---output-R- flux de chaleur sensible (CpT) J/m2/s (W/m2)
227	! (orientation positive vers le bas)
228	! tke---input/output-R- tke (kg/m**2/s)
229	! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
230	! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
231	! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
232	! dflux_t--output-R- derive du flux sensible
233	! dflux_q--output-R- derive du flux latent
234	! zu1------output-R- le vent dans la premiere couche
235	! zv1------output-R- le vent dans la premiere couche
236	! trmb1----output-R- deep_cape
237	! trmb2----output-R- inhibition
238	! trmb3----output-R- Point Omega
239	! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL
240	! plcl-----output-R- Niveau de condensation
241	! pblh-----output-R- HCL
242	! pblT-----output-R- T au nveau HCL
243	!
244	USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
245	IMPLICIT NONE
246
247	INCLUDE "indicesol.h"
248	INCLUDE "dimsoil.h"
249	INCLUDE "YOMCST.h"
250	INCLUDE "iniprint.h"
251	INCLUDE "FCTTRE.h"
252	INCLUDE "clesphys.h"
253	INCLUDE "compbl.h"
254	INCLUDE "dimensions.h"
255	INCLUDE "YOETHF.h"
256	INCLUDE "temps.h"
257	INCLUDE "control.h"
258	! Input variables
259	!****************************************************************************************
260	REAL, INTENT(IN) :: dtime ! time interval (s)
261	REAL, INTENT(IN) :: date0 ! initial day
262	INTEGER, INTENT(IN) :: itap ! time step
263	INTEGER, INTENT(IN) :: jour ! current day of the year
264	LOGICAL, INTENT(IN) :: debut ! true if first run step
265	LOGICAL, INTENT(IN) :: lafin ! true if last run step
266	REAL, DIMENSION(klon), INTENT(IN) :: rlon ! longitudes in degrees
267	REAL, DIMENSION(klon), INTENT(IN) :: rlat ! latitudes in degrees
268	REAL, DIMENSION(klon), INTENT(IN) :: rugoro ! rugosity length
269	REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cosine of solar zenith angle
270	REAL, DIMENSION(klon), INTENT(IN) :: rain_f ! rain fall
271	REAL, DIMENSION(klon), INTENT(IN) :: snow_f ! snow fall
272	REAL, DIMENSION(klon), INTENT(IN) :: solsw_m ! net shortwave radiation at mean surface
273	REAL, DIMENSION(klon), INTENT(IN) :: sollw_m ! net longwave radiation at mean surface
274	REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K)
275	REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! water vapour (kg/kg)
276	REAL, DIMENSION(klon,klev), INTENT(IN) :: u ! u speed
277	REAL, DIMENSION(klon,klev), INTENT(IN) :: v ! v speed
278	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pression (Pa)
279	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression between layers (Pa)
280	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! sub-surface fraction
281
282	! Input/Output variables
283	!****************************************************************************************
284	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ts ! temperature at surface (K)
285	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: alb1 ! albedo in visible SW interval
286	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: alb2 ! albedo in near infra-red SW interval
287	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: u10m ! u speed at 10m
288	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: v10m ! v speed at 10m
289	REAL, DIMENSION(klon, klev+1, nbsrf), INTENT(INOUT) :: tke
290
291	! Output variables
292	!****************************************************************************************
293	REAL, DIMENSION(klon), INTENT(OUT) :: lwdown_m ! Downcoming longwave radiation
294	REAL, DIMENSION(klon), INTENT(OUT) :: cdragh ! drag coefficient for T and Q
295	REAL, DIMENSION(klon), INTENT(OUT) :: cdragm ! drag coefficient for wind
296	REAL, DIMENSION(klon), INTENT(OUT) :: zu1 ! u wind speed in first layer
297	REAL, DIMENSION(klon), INTENT(OUT) :: zv1 ! v wind speed in first layer
298	REAL, DIMENSION(klon), INTENT(OUT) :: alb1_m ! mean albedo in visible SW interval
299	REAL, DIMENSION(klon), INTENT(OUT) :: alb2_m ! mean albedo in near IR SW interval
300	REAL, DIMENSION(klon), INTENT(OUT) :: zxsens ! sensible heat flux at surface with inversed sign
301	! (=> positive sign upwards)
302	REAL, DIMENSION(klon), INTENT(OUT) :: zxevap ! water vapour flux at surface, positiv upwards
303	REAL, DIMENSION(klon), INTENT(OUT) :: zxtsol ! temperature at surface, mean for each grid point
304	REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat ! latent flux, mean for each grid point
305	REAL, DIMENSION(klon), INTENT(OUT) :: zt2m ! temperature at 2m, mean for each grid point
306	REAL, DIMENSION(klon), INTENT(OUT) :: qsat2m
307	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! change in temperature
308	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_q ! change in water vapour
309	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u ! change in u speed
310	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v ! change in v speed
311	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zcoefh ! coef for turbulent diffusion of T and Q, mean for each grid point
312
313	! Output only for diagnostics
314	REAL, DIMENSION(klon), INTENT(OUT) :: slab_wfbils! heat balance at surface only for slab at ocean points
315	REAL, DIMENSION(klon), INTENT(OUT) :: qsol_d ! water height in the soil (mm)
316	REAL, DIMENSION(klon), INTENT(OUT) :: zq2m ! water vapour at 2m, mean for each grid point
317	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh ! height of the planetary boundary layer(HPBL)
318	REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl ! condensation level
319	REAL, DIMENSION(klon), INTENT(OUT) :: s_capCL ! CAPE of PBL
320	REAL, DIMENSION(klon), INTENT(OUT) :: s_oliqCL ! liquid water intergral of PBL
321	REAL, DIMENSION(klon), INTENT(OUT) :: s_cteiCL ! cloud top instab. crit. of PBL
322	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblT ! temperature at PBLH
323	REAL, DIMENSION(klon), INTENT(OUT) :: s_therm ! thermal virtual temperature excess
324	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb1 ! deep cape, mean for each grid point
325	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb2 ! inhibition, mean for each grid point
326	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb3 ! point Omega, mean for each grid point
327	REAL, DIMENSION(klon), INTENT(OUT) :: zxrugs ! rugosity at surface (m), mean for each grid point
328	REAL, DIMENSION(klon), INTENT(OUT) :: zu10m ! u speed at 10m, mean for each grid point
329	REAL, DIMENSION(klon), INTENT(OUT) :: zv10m ! v speed at 10m, mean for each grid point
330	REAL, DIMENSION(klon), INTENT(OUT) :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
331	REAL, DIMENSION(klon), INTENT(OUT) :: zxqsurf ! humidity at surface, mean for each grid point
332	REAL, DIMENSION(klon), INTENT(OUT) :: rh2m ! relative humidity at 2m
333	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxu ! u wind tension, mean for each grid point
334	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxv ! v wind tension, mean for each grid point
335	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_d ! rugosity length (m)
336	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_d ! age of snow at surface
337	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: solsw ! net shortwave radiation at surface
338	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: sollw ! net longwave radiation at surface
339	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: d_ts ! change in temperature at surface
340	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_d ! evaporation at surface
341	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: fluxlat ! latent flux
342	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: t2m ! temperature at 2 meter height
343	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbils ! heat balance at surface
344	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbilo ! water balance at surface
345	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t ! sensible heat flux (CpT) J/m2/s (W/m2)
346	! positve orientation downwards
347	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u ! u wind tension (kg m/s)/(m**2 s) or Pascal
348	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v ! v wind tension (kg m/s)/(m**2 s) or Pascal
349
350	! Output not needed
351	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_t ! change of sensible heat flux
352	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_q ! change of water vapour flux
353	REAL, DIMENSION(klon), INTENT(OUT) :: zxsnow ! snow at surface, mean for each grid point
354	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxt ! sensible heat flux, mean for each grid point
355	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxq ! water vapour flux, mean for each grid point
356	REAL, DIMENSION(klon, nbsrf),INTENT(OUT) :: q2m ! water vapour at 2 meter height
357	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q ! water vapour flux(latent flux) (kg/m**2/s)
358
359
360	! Local variables with attribute SAVE
361	!****************************************************************************************
362	INTEGER, SAVE :: nhoridbg, nidbg ! variables for IOIPSL
363	!$OMP THREADPRIVATE(nhoridbg, nidbg)
364	LOGICAL, SAVE :: debugindex=.FALSE.
365	!$OMP THREADPRIVATE(debugindex)
366	LOGICAL, SAVE :: first_call=.TRUE.
367	!$OMP THREADPRIVATE(first_call)
368	CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
369	!$OMP THREADPRIVATE(cl_surf)
370
371	! Other local variables
372	!****************************************************************************************
373	INTEGER :: i, k, nsrf
374	INTEGER :: knon, j
375	INTEGER :: idayref
376	INTEGER , DIMENSION(klon) :: ni
377	REAL :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
378	REAL :: amn, amx
379	REAL :: f1 ! fraction de longeurs visibles parmi tout SW intervalle
380	REAL, DIMENSION(klon) :: r_co2_ppm ! taux CO2 atmosphere
381	REAL, DIMENSION(klon) :: yts, yrugos, ypct, yz0_new
382	REAL, DIMENSION(klon) :: yalb, yalb1, yalb2
383	REAL, DIMENSION(klon) :: yu1, yv1
384	REAL, DIMENSION(klon) :: ysnow, yqsurf, yagesno, yqsol
385	REAL, DIMENSION(klon) :: yrain_f, ysnow_f
386	REAL, DIMENSION(klon) :: ysolsw, ysollw
387	REAL, DIMENSION(klon) :: yfder
388	REAL, DIMENSION(klon) :: yrugoro
389	REAL, DIMENSION(klon) :: yfluxlat
390	REAL, DIMENSION(klon) :: y_d_ts
391	REAL, DIMENSION(klon) :: y_flux_t1, y_flux_q1
392	REAL, DIMENSION(klon) :: y_dflux_t, y_dflux_q
393	REAL, DIMENSION(klon) :: y_flux_u1, y_flux_v1
394	REAL, DIMENSION(klon) :: yt2m, yq2m, yu10m
395	REAL, DIMENSION(klon) :: yustar
396	REAL, DIMENSION(klon) :: ywindsp
397	REAL, DIMENSION(klon) :: yt10m, yq10m
398	REAL, DIMENSION(klon) :: ypblh
399	REAL, DIMENSION(klon) :: ylcl
400	REAL, DIMENSION(klon) :: ycapCL
401	REAL, DIMENSION(klon) :: yoliqCL
402	REAL, DIMENSION(klon) :: ycteiCL
403	REAL, DIMENSION(klon) :: ypblT
404	REAL, DIMENSION(klon) :: ytherm
405	REAL, DIMENSION(klon) :: ytrmb1
406	REAL, DIMENSION(klon) :: ytrmb2
407	REAL, DIMENSION(klon) :: ytrmb3
408	REAL, DIMENSION(klon) :: uzon, vmer
409	REAL, DIMENSION(klon) :: tair1, qair1, tairsol
410	REAL, DIMENSION(klon) :: psfce, patm
411	REAL, DIMENSION(klon) :: qairsol, zgeo1
412	REAL, DIMENSION(klon) :: rugo1
413	REAL, DIMENSION(klon) :: yfluxsens
414	REAL, DIMENSION(klon) :: AcoefH, AcoefQ, BcoefH, BcoefQ
415	REAL, DIMENSION(klon) :: AcoefU, AcoefV, BcoefU, BcoefV
416	REAL, DIMENSION(klon) :: ypsref
417	REAL, DIMENSION(klon) :: yevap, ytsurf_new, yalb1_new, yalb2_new
418	REAL, DIMENSION(klon) :: ztsol
419	REAL, DIMENSION(klon) :: alb_m ! mean albedo for whole SW interval
420	REAL, DIMENSION(klon,klev) :: y_d_t, y_d_q
421	REAL, DIMENSION(klon,klev) :: y_d_u, y_d_v
422	REAL, DIMENSION(klon,klev) :: y_flux_t, y_flux_q
423	REAL, DIMENSION(klon,klev) :: y_flux_u, y_flux_v
424	REAL, DIMENSION(klon,klev) :: ycoefh, ycoefm
425	REAL, DIMENSION(klon) :: ycdragh, ycdragm
426	REAL, DIMENSION(klon,klev) :: yu, yv
427	REAL, DIMENSION(klon,klev) :: yt, yq
428	REAL, DIMENSION(klon,klev) :: ypplay, ydelp
429	REAL, DIMENSION(klon,klev) :: delp
430	REAL, DIMENSION(klon,klev+1) :: ypaprs
431	REAL, DIMENSION(klon,klev+1) :: ytke
432	REAL, DIMENSION(klon,nsoilmx) :: ytsoil
433	CHARACTER(len=80) :: abort_message
434	CHARACTER(len=20) :: modname = 'pbl_surface'
435	LOGICAL, PARAMETER :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
436	LOGICAL, PARAMETER :: check=.FALSE.
437
438	! For debugging with IOIPSL
439	INTEGER, DIMENSION(iim*(jjm+1)) :: ndexbg
440	REAL :: zjulian
441	REAL, DIMENSION(klon) :: tabindx
442	REAL, DIMENSION(iim,jjm+1) :: zx_lon, zx_lat
443	REAL, DIMENSION(iim,jjm+1) :: debugtab
444
445
446	REAL, DIMENSION(klon,nbsrf) :: pblh ! height of the planetary boundary layer
447	REAL, DIMENSION(klon,nbsrf) :: plcl ! condensation level
448	REAL, DIMENSION(klon,nbsrf) :: capCL
449	REAL, DIMENSION(klon,nbsrf) :: oliqCL
450	REAL, DIMENSION(klon,nbsrf) :: cteiCL
451	REAL, DIMENSION(klon,nbsrf) :: pblT
452	REAL, DIMENSION(klon,nbsrf) :: therm
453	REAL, DIMENSION(klon,nbsrf) :: trmb1 ! deep cape
454	REAL, DIMENSION(klon,nbsrf) :: trmb2 ! inhibition
455	REAL, DIMENSION(klon,nbsrf) :: trmb3 ! point Omega
456	REAL, DIMENSION(klon,nbsrf) :: zx_rh2m, zx_qsat2m
457	REAL, DIMENSION(klon,nbsrf) :: zx_t1
458	REAL, DIMENSION(klon, nbsrf) :: alb ! mean albedo for whole SW interval
459	REAL, DIMENSION(klon) :: ylwdown ! jg : temporary (ysollwdown)
460
461	REAL :: zx_qs1, zcor1, zdelta1
462
463	!****************************************************************************************
464	! Declarations specifiques pour le 1D. A reprendre
465	REAL,SAVE :: fsens,flat
466	LOGICAL,SAVE :: ok_flux_surf=.FALSE.
467	!$OMP THREADPRIVATE(fsens,flat,ok_flux_surf)
468
469	!****************************************************************************************
470	! End of declarations
471	!****************************************************************************************
472
473
474	!****************************************************************************************
475	! 1) Initialisation and validation tests
476	! Only done first time entering this subroutine
477	!
478	!****************************************************************************************
479
480	IF (first_call) THEN
481	first_call=.FALSE.
482
483	! Initilize debug IO
484	IF (debugindex .AND. mpi_size==1) THEN
485	! initialize IOIPSL output
486	idayref = day_ini
487	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
488	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
489	DO i = 1, iim
490	zx_lon(i,1) = rlon(i+1)
491	zx_lon(i,jjm+1) = rlon(i+1)
492	ENDDO
493	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
494	CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
495	1,iim,1,jjm+1, &
496	itau_phy,zjulian,dtime,nhoridbg,nidbg)
497	! no vertical axis
498	cl_surf(1)='ter'
499	cl_surf(2)='lic'
500	cl_surf(3)='oce'
501	cl_surf(4)='sic'
502	DO nsrf=1,nbsrf
503	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
504	jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
505	END DO
506
507	CALL histend(nidbg)
508	CALL histsync(nidbg)
509
510	END IF
511
512	ENDIF
513
514	!****************************************************************************************
515	! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
516	! instead of ORCHIDEE)
517	IF (qsol0>0.) THEN
518	PRINT*,'WARNING : On impose qsol=',qsol0
519	qsol(:)=qsol0
520	ENDIF
521	!****************************************************************************************
522
523	!****************************************************************************************
524	! 2) Initialization to zero
525	! Done for all local variables that will be compressed later
526	! and argument with INTENT(OUT)
527	!****************************************************************************************
528	cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0
529	ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0
530	zv1 = 0.0 ; yqsurf = 0.0 ; yalb1 = 0.0 ; yalb2 = 0.0
531	yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0
532	ysollw = 0.0 ; yrugos = 0.0 ; yu1 = 0.0
533	yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0
534	ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0
535	yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
536	yrugoro = 0.0 ; ywindsp = 0.0
537	d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0
538	flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0
539	d_u = 0.0 ; d_v = 0.0 ; yqsol = 0.0
540	ytherm = 0.0 ; ytke=0.
541
542	zcoefh(:,:) = 0.0
543	zcoefh(:,1) = 999999. ! zcoefh(:,k=1) should never be used
544	ytsoil = 999999.
545
546	rh2m(:) = 0.
547	qsat2m(:) = 0.
548	!****************************************************************************************
549	! 3) - Calculate pressure thickness of each layer
550	! - Calculate the wind at first layer
551	! - Mean calculations of albedo
552	! - Calculate net radiance at sub-surface
553	!****************************************************************************************
554	DO k = 1, klev
555	DO i = 1, klon
556	delp(i,k) = paprs(i,k)-paprs(i,k+1)
557	ENDDO
558	ENDDO
559
560	!****************************************************************************************
561	! Test for rugos........ from physiq.. A la fin plutot???
562	!
563	!****************************************************************************************
564
565	zxrugs(:) = 0.0
566	DO nsrf = 1, nbsrf
567	DO i = 1, klon
568	rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015)
569	zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf)
570	ENDDO
571	ENDDO
572
573	! Mean calculations of albedo
574	!
575	! Albedo at sub-surface
576	! * alb1 : albedo in visible SW interval
577	! * alb2 : albedo in near infrared SW interval
578	! * alb : mean albedo for whole SW interval
579	!
580	! Mean albedo for grid point
581	! * alb1_m : albedo in visible SW interval
582	! * alb2_m : albedo in near infrared SW interval
583	! * alb_m : mean albedo at whole SW interval
584
585	alb1_m(:) = 0.0
586	alb2_m(:) = 0.0
587	DO nsrf = 1, nbsrf
588	DO i = 1, klon
589	alb1_m(i) = alb1_m(i) + alb1(i,nsrf) * pctsrf(i,nsrf)
590	alb2_m(i) = alb2_m(i) + alb2(i,nsrf) * pctsrf(i,nsrf)
591	ENDDO
592	ENDDO
593
594	! We here suppose the fraction f1 of incoming radiance of visible radiance
595	! as a fraction of all shortwave radiance
596	f1 = 0.5
597	! f1 = 1 ! put f1=1 to recreate old calculations
598
599	DO nsrf = 1, nbsrf
600	DO i = 1, klon
601	alb(i,nsrf) = f1alb1(i,nsrf) + (1-f1)alb2(i,nsrf)
602	ENDDO
603	ENDDO
604
605	DO i = 1, klon
606	alb_m(i) = f1alb1_m(i) + (1-f1)alb2_m(i)
607	END DO
608
609	! Calculation of mean temperature at surface grid points
610	ztsol(:) = 0.0
611	DO nsrf = 1, nbsrf
612	DO i = 1, klon
613	ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
614	ENDDO
615	ENDDO
616
617	! Linear distrubution on sub-surface of long- and shortwave net radiance
618	DO nsrf = 1, nbsrf
619	DO i = 1, klon
620	sollw(i,nsrf) = sollw_m(i) + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ts(i,nsrf))
621	solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
622	ENDDO
623	ENDDO
624
625
626	! Downwelling longwave radiation at mean surface
627	lwdown_m(:) = 0.0
628	DO i = 1, klon
629	lwdown_m(i) = sollw_m(i) + RSIGMAztsol(i)*4
630	ENDDO
631
632	!****************************************************************************************
633	! 4) Loop over different surfaces
634	!
635	! Only points containing a fraction of the sub surface will be threated.
636	!
637	!****************************************************************************************
638
639	loop_nbsrf: DO nsrf = 1, nbsrf
640
641	! Search for index(ni) and size(knon) of domaine to treat
642	ni(:) = 0
643	knon = 0
644	DO i = 1, klon
645	IF (pctsrf(i,nsrf) > 0.) THEN
646	knon = knon + 1
647	ni(knon) = i
648	ENDIF
649	ENDDO
650
651	! write index, with IOIPSL
652	IF (debugindex .AND. mpi_size==1) THEN
653	tabindx(:)=0.
654	DO i=1,knon
655	tabindx(i)=FLOAT(i)
656	END DO
657	debugtab(:,:) = 0.
658	ndexbg(:) = 0
659	CALL gath2cpl(tabindx,debugtab,knon,ni)
660	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg)
661	ENDIF
662
663	!****************************************************************************************
664	! 5) Compress variables
665	!
666	!****************************************************************************************
667
668	DO j = 1, knon
669	i = ni(j)
670	ypct(j) = pctsrf(i,nsrf)
671	yts(j) = ts(i,nsrf)
672	ysnow(j) = snow(i,nsrf)
673	yqsurf(j) = qsurf(i,nsrf)
674	yalb(j) = alb(i,nsrf)
675	yalb1(j) = alb1(i,nsrf)
676	yalb2(j) = alb2(i,nsrf)
677	yrain_f(j) = rain_f(i)
678	ysnow_f(j) = snow_f(i)
679	yagesno(j) = agesno(i,nsrf)
680	yfder(j) = fder(i)
681	ysolsw(j) = solsw(i,nsrf)
682	ysollw(j) = sollw(i,nsrf)
683	yrugos(j) = rugos(i,nsrf)
684	yrugoro(j) = rugoro(i)
685	yu1(j) = u(i,1)
686	yv1(j) = v(i,1)
687	ypaprs(j,klev+1) = paprs(i,klev+1)
688	ywindsp(j) = SQRT(u10m(i,nsrf)2 + v10m(i,nsrf)2 )
689	END DO
690
691	DO k = 1, klev
692	DO j = 1, knon
693	i = ni(j)
694	ypaprs(j,k) = paprs(i,k)
695	ypplay(j,k) = pplay(i,k)
696	ydelp(j,k) = delp(i,k)
697	ytke(j,k) = tke(i,k,nsrf)
698	yu(j,k) = u(i,k)
699	yv(j,k) = v(i,k)
700	yt(j,k) = t(i,k)
701	yq(j,k) = q(i,k)
702	ENDDO
703	ENDDO
704
705	DO k = 1, nsoilmx
706	DO j = 1, knon
707	i = ni(j)
708	ytsoil(j,k) = ftsoil(i,k,nsrf)
709	END DO
710	END DO
711
712	! qsol(water height in soil) only for bucket continental model
713	IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
714	DO j = 1, knon
715	i = ni(j)
716	yqsol(j) = qsol(i)
717	END DO
718	ENDIF
719
720	!****************************************************************************************
721	! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
722	!
723	!****************************************************************************************
724
725	CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
726	yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
727	yts, yqsurf, yrugos, &
728	ycdragm, ycdragh )
729
730	!****************************************************************************************
731	! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefm et ycoefm.
732	!
733	!****************************************************************************************
734
735	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
736	ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
737	ycoefm, ycoefh, ytke)
738
739	!****************************************************************************************
740	!
741	! 8) "La descente" - "The downhill"
742	!
743	! climb_hq_down and climb_wind_down calculate the coefficients
744	! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
745	! Only the coefficients at surface for H and Q are returned.
746	!
747	!****************************************************************************************
748
749	! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
750	CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
751	ydelp, yt, yq, dtime, &
752	AcoefH, AcoefQ, BcoefH, BcoefQ)
753
754	! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
755	CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
756	AcoefU, AcoefV, BcoefU, BcoefV)
757
758
759	!****************************************************************************************
760	! 9) Small calculations
761	!
762	!****************************************************************************************
763
764	! - Reference pressure is given the values at surface level
765	ypsref(:) = ypaprs(:,1)
766
767	! - CO2 field on 2D grid to be sent to ORCHIDEE
768	! Transform to compressed field
769	IF (carbon_cycle_cpl) THEN
770	DO i=1,knon
771	r_co2_ppm(i) = co2_send(ni(i))
772	END DO
773	ELSE
774	r_co2_ppm(:) = co2_ppm ! Constant field
775	END IF
776
777	!****************************************************************************************
778	!
779	! Calulate t2m and q2m for the case of calculation at land grid points
780	! t2m and q2m are needed as input to ORCHIDEE
781	!
782	!****************************************************************************************
783	IF (nsrf == is_ter) THEN
784
785	DO i = 1, knon
786	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
787	* (ypaprs(i,1)-ypplay(i,1))
788	END DO
789
790	! Calculate the temperature et relative humidity at 2m and the wind at 10m
791	CALL stdlevvar(klon, knon, is_ter, zxli, &
792	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
793	yts, yqsurf, yrugos, ypaprs(:,1), ypplay(:,1), &
794	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
795
796	END IF
797
798	!****************************************************************************************
799	!
800	! 10) Switch selon current surface
801	! It is necessary to start with the continental surfaces because the ocean
802	! needs their run-off.
803	!
804	!****************************************************************************************
805	SELECT CASE(nsrf)
806
807	CASE(is_ter)
808	! ylwdown : to be removed, calculation is now done at land surface in surf_land
809	ylwdown(:)=0.0
810	DO i=1,knon
811	ylwdown(i)=lwdown_m(ni(i))
812	END DO
813	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
814	rlon, rlat, &
815	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
816	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
817	AcoefH, AcoefQ, BcoefH, BcoefQ, &
818	AcoefU, AcoefV, BcoefU, BcoefV, &
819	ypsref, yu1, yv1, yrugoro, pctsrf, &
820	ylwdown, yq2m, yt2m, &
821	ysnow, yqsol, yagesno, ytsoil, &
822	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
823	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
824	y_flux_u1, y_flux_v1 )
825
826
827	CASE(is_lic)
828	CALL surf_landice(itap, dtime, knon, ni, &
829	ysolsw, ysollw, yts, ypplay(:,1), &
830	ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
831	AcoefH, AcoefQ, BcoefH, BcoefQ, &
832	AcoefU, AcoefV, BcoefU, BcoefV, &
833	ypsref, yu1, yv1, yrugoro, pctsrf, &
834	ysnow, yqsurf, yqsol, yagesno, &
835	ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
836	ytsurf_new, y_dflux_t, y_dflux_q, &
837	y_flux_u1, y_flux_v1)
838
839	CASE(is_oce)
840	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb1, &
841	yrugos, ywindsp, rmu0, yfder, yts, &
842	itap, dtime, jour, knon, ni, &
843	ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
844	AcoefH, AcoefQ, BcoefH, BcoefQ, &
845	AcoefU, AcoefV, BcoefU, BcoefV, &
846	ypsref, yu1, yv1, yrugoro, pctsrf, &
847	ysnow, yqsurf, yagesno, &
848	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
849	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
850	y_flux_u1, y_flux_v1)
851
852	CASE(is_sic)
853	CALL surf_seaice( &
854	rlon, rlat, ysolsw, ysollw, yalb1, yfder, &
855	itap, dtime, jour, knon, ni, &
856	lafin, &
857	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
858	AcoefH, AcoefQ, BcoefH, BcoefQ, &
859	AcoefU, AcoefV, BcoefU, BcoefV, &
860	ypsref, yu1, yv1, yrugoro, pctsrf, &
861	ysnow, yqsurf, yqsol, yagesno, ytsoil, &
862	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
863	ytsurf_new, y_dflux_t, y_dflux_q, &
864	y_flux_u1, y_flux_v1)
865
866
867	CASE DEFAULT
868	WRITE(lunout,*) 'Surface index = ', nsrf
869	abort_message = 'Surface index not valid'
870	CALL abort_gcm(modname,abort_message,1)
871	END SELECT
872
873
874	!****************************************************************************************
875	! 11) - Calcul the increment of surface temperature
876	!
877	!****************************************************************************************
878	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
879
880	!****************************************************************************************
881	!
882	! 12) "La remontee" - "The uphill"
883	!
884	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
885	! for X=H, Q, U and V, for all vertical levels.
886	!
887	!****************************************************************************************
888	! H and Q
889	IF (ok_flux_surf) THEN
890	PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
891	y_flux_t1(:) = fsens
892	y_flux_q1(:) = flat/RLVTT
893	yfluxlat(:) = flat
894	ELSE
895	y_flux_t1(:) = yfluxsens(:)
896	y_flux_q1(:) = -yevap(:)
897	ENDIF
898
899	CALL climb_hq_up(knon, dtime, yt, yq, &
900	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
901	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
902
903
904	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
905	y_flux_u, y_flux_v, y_d_u, y_d_v)
906
907
908	DO j = 1, knon
909	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
910	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
911	ENDDO
912
913	!****************************************************************************************
914	! 13) Transform variables for output format :
915	! - Decompress
916	! - Multiply with pourcentage of current surface
917	! - Cumulate in global variable
918	!
919	!****************************************************************************************
920
921	tke(:,:,nsrf) = 0.
922	DO k = 1, klev
923	DO j = 1, knon
924	i = ni(j)
925	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
926	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
927	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
928	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
929
930	flux_t(i,k,nsrf) = y_flux_t(j,k)
931	flux_q(i,k,nsrf) = y_flux_q(j,k)
932	flux_u(i,k,nsrf) = y_flux_u(j,k)
933	flux_v(i,k,nsrf) = y_flux_v(j,k)
934
935	tke(i,k,nsrf) = ytke(j,k)
936
937	ENDDO
938	ENDDO
939
940	evap(:,nsrf) = - flux_q(:,1,nsrf)
941
942	alb1(:, nsrf) = 0.
943	alb2(:, nsrf) = 0.
944	snow(:, nsrf) = 0.
945	qsurf(:, nsrf) = 0.
946	rugos(:, nsrf) = 0.
947	fluxlat(:,nsrf) = 0.
948	DO j = 1, knon
949	i = ni(j)
950	d_ts(i,nsrf) = y_d_ts(j)
951	alb1(i,nsrf) = yalb1_new(j)
952	alb2(i,nsrf) = yalb2_new(j)
953	snow(i,nsrf) = ysnow(j)
954	qsurf(i,nsrf) = yqsurf(j)
955	rugos(i,nsrf) = yz0_new(j)
956	fluxlat(i,nsrf) = yfluxlat(j)
957	agesno(i,nsrf) = yagesno(j)
958	cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
959	cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
960	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
961	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
962	END DO
963
964	DO k = 2, klev
965	DO j = 1, knon
966	i = ni(j)
967	zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)*ypct(j)
968	END DO
969	END DO
970
971	IF ( nsrf .EQ. is_ter ) THEN
972	DO j = 1, knon
973	i = ni(j)
974	qsol(i) = yqsol(j)
975	END DO
976	END IF
977
978	ftsoil(:,:,nsrf) = 0.
979	DO k = 1, nsoilmx
980	DO j = 1, knon
981	i = ni(j)
982	ftsoil(i, k, nsrf) = ytsoil(j,k)
983	END DO
984	END DO
985
986
987	DO k = 1, klev
988	DO j = 1, knon
989	i = ni(j)
990	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
991	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
992	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
993	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
994	END DO
995	END DO
996
997	!****************************************************************************************
998	! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m
999	! Call HBTM
1000	!
1001	!****************************************************************************************
1002	t2m(:,nsrf) = 0.
1003	q2m(:,nsrf) = 0.
1004	u10m(:,nsrf) = 0.
1005	v10m(:,nsrf) = 0.
1006
1007	pblh(:,nsrf) = 0. ! Hauteur de couche limite
1008	plcl(:,nsrf) = 0. ! Niveau de condensation de la CLA
1009	capCL(:,nsrf) = 0. ! CAPE de couche limite
1010	oliqCL(:,nsrf) = 0. ! eau_liqu integree de couche limite
1011	cteiCL(:,nsrf) = 0. ! cloud top instab. crit. couche limite
1012	pblt(:,nsrf) = 0. ! T a la Hauteur de couche limite
1013	therm(:,nsrf) = 0.
1014	trmb1(:,nsrf) = 0. ! deep_cape
1015	trmb2(:,nsrf) = 0. ! inhibition
1016	trmb3(:,nsrf) = 0. ! Point Omega
1017
1018	#undef T2m
1019	#define T2m
1020	#ifdef T2m
1021	! Calculations of diagnostic t,q at 2m and u, v at 10m
1022
1023	DO j=1, knon
1024	i = ni(j)
1025	uzon(j) = yu(j,1) + y_d_u(j,1)
1026	vmer(j) = yv(j,1) + y_d_v(j,1)
1027	tair1(j) = yt(j,1) + y_d_t(j,1)
1028	qair1(j) = yq(j,1) + y_d_q(j,1)
1029	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
1030	* (ypaprs(j,1)-ypplay(j,1))
1031	tairsol(j) = yts(j) + y_d_ts(j)
1032	rugo1(j) = yrugos(j)
1033	IF(nsrf.EQ.is_oce) THEN
1034	rugo1(j) = rugos(i,nsrf)
1035	ENDIF
1036	psfce(j)=ypaprs(j,1)
1037	patm(j)=ypplay(j,1)
1038	qairsol(j) = yqsurf(j)
1039	END DO
1040
1041
1042	! Calculate the temperature et relative humidity at 2m and the wind at 10m
1043	CALL stdlevvar(klon, knon, nsrf, zxli, &
1044	uzon, vmer, tair1, qair1, zgeo1, &
1045	tairsol, qairsol, rugo1, psfce, patm, &
1046	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
1047
1048	DO j=1, knon
1049	i = ni(j)
1050	t2m(i,nsrf)=yt2m(j)
1051	q2m(i,nsrf)=yq2m(j)
1052
1053	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
1054	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
1055	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
1056	END DO
1057
1058	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
1059	!IM Ajoute dependance type surface
1060	IF (thermcep) THEN
1061	DO j = 1, knon
1062	i=ni(j)
1063	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
1064	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
1065	zx_qs1 = MIN(0.5,zx_qs1)
1066	zcor1 = 1./(1.-RETV*zx_qs1)
1067	zx_qs1 = zx_qs1*zcor1
1068
1069	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
1070	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
1071	END DO
1072	END IF
1073
1074	CALL HBTM(knon, ypaprs, ypplay, &
1075	yt2m,yt10m,yq2m,yq10m,yustar, &
1076	y_flux_t,y_flux_q,yu,yv,yt,yq, &
1077	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
1078	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
1079
1080	DO j=1, knon
1081	i = ni(j)
1082	pblh(i,nsrf) = ypblh(j)
1083	plcl(i,nsrf) = ylcl(j)
1084	capCL(i,nsrf) = ycapCL(j)
1085	oliqCL(i,nsrf) = yoliqCL(j)
1086	cteiCL(i,nsrf) = ycteiCL(j)
1087	pblT(i,nsrf) = ypblT(j)
1088	therm(i,nsrf) = ytherm(j)
1089	trmb1(i,nsrf) = ytrmb1(j)
1090	trmb2(i,nsrf) = ytrmb2(j)
1091	trmb3(i,nsrf) = ytrmb3(j)
1092	END DO
1093
1094	#else
1095	! T2m not defined
1096	! No calculation
1097	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
1098	#endif
1099
1100	!****************************************************************************************
1101	! 15) End of loop over different surfaces
1102	!
1103	!****************************************************************************************
1104	END DO loop_nbsrf
1105
1106	!****************************************************************************************
1107	! 16) Calculate the mean value over all sub-surfaces for som variables
1108	!
1109	!****************************************************************************************
1110
1111	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
1112	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
1113	DO nsrf = 1, nbsrf
1114	DO k = 1, klev
1115	DO i = 1, klon
1116	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
1117	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
1118	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
1119	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
1120	END DO
1121	END DO
1122	END DO
1123
1124	DO i = 1, klon
1125	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1126	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1127	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
1128	ENDDO
1129
1130	!
1131	! Incrementer la temperature du sol
1132	!
1133	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
1134	zt2m(:) = 0.0 ; zq2m(:) = 0.0
1135	zu10m(:) = 0.0 ; zv10m(:) = 0.0
1136	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
1137	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
1138	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
1139	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
1140	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
1141
1142
1143	DO nsrf = 1, nbsrf
1144	DO i = 1, klon
1145	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
1146
1147	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
1148	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1149	wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
1150	pctsrf(i,nsrf)
1151
1152	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
1153	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
1154
1155	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
1156	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
1157	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
1158	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
1159
1160	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
1161	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
1162	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
1163	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
1164	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
1165	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
1166	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
1167	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
1168	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
1169	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
1170	END DO
1171	END DO
1172
1173	IF (check) THEN
1174	amn=MIN(ts(1,is_ter),1000.)
1175	amx=MAX(ts(1,is_ter),-1000.)
1176	DO i=2, klon
1177	amn=MIN(ts(i,is_ter),amn)
1178	amx=MAX(ts(i,is_ter),amx)
1179	ENDDO
1180	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
1181	ENDIF
1182
1183	!jg ?
1184	!!$!
1185	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
1186	!!$! sub-surfaces is distributed.
1187	!!$!
1188	!!$ DO nsrf = 1, nbsrf
1189	!!$ DO i = 1, klon
1190	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
1191	!!$ ts(i,nsrf) = zxtsol(i)
1192	!!$ t2m(i,nsrf) = zt2m(i)
1193	!!$ q2m(i,nsrf) = zq2m(i)
1194	!!$ u10m(i,nsrf) = zu10m(i)
1195	!!$ v10m(i,nsrf) = zv10m(i)
1196	!!$
1197	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
1198	!!$ pblh(i,nsrf) = s_pblh(i)
1199	!!$ plcl(i,nsrf) = s_plcl(i)
1200	!!$ capCL(i,nsrf) = s_capCL(i)
1201	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
1202	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
1203	!!$ pblT(i,nsrf) = s_pblT(i)
1204	!!$ therm(i,nsrf) = s_therm(i)
1205	!!$ trmb1(i,nsrf) = s_trmb1(i)
1206	!!$ trmb2(i,nsrf) = s_trmb2(i)
1207	!!$ trmb3(i,nsrf) = s_trmb3(i)
1208	!!$ ENDIF
1209	!!$ ENDDO
1210	!!$ ENDDO
1211
1212
1213	DO i = 1, klon
1214	fder(i) = - 4.0RSIGMAzxtsol(i)**3
1215	ENDDO
1216
1217	zxqsurf(:) = 0.0
1218	zxsnow(:) = 0.0
1219	DO nsrf = 1, nbsrf
1220	DO i = 1, klon
1221	zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
1222	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
1223	END DO
1224	END DO
1225
1226	! Premier niveau de vent sortie dans physiq.F
1227	zu1(:) = u(:,1)
1228	zv1(:) = v(:,1)
1229
1230	! Some of the module declared variables are returned for printing in physiq.F
1231	qsol_d(:) = qsol(:)
1232	evap_d(:,:) = evap(:,:)
1233	rugos_d(:,:) = rugos(:,:)
1234	agesno_d(:,:) = agesno(:,:)
1235
1236
1237	END SUBROUTINE pbl_surface
1238	!
1239	!****************************************************************************************
1240	!
1241	SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
1242	evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
1243
1244	INCLUDE "indicesol.h"
1245	INCLUDE "dimsoil.h"
1246
1247	! Ouput variables
1248	!****************************************************************************************
1249	REAL, DIMENSION(klon), INTENT(OUT) :: qsol_rst
1250	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
1251	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
1252	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
1253	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_rst
1254	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_rst
1255	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_rst
1256	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
1257
1258
1259	!****************************************************************************************
1260	! Return module variables for writing to restart file
1261	!
1262	!****************************************************************************************
1263	qsol_rst(:) = qsol(:)
1264	fder_rst(:) = fder(:)
1265	snow_rst(:,:) = snow(:,:)
1266	qsurf_rst(:,:) = qsurf(:,:)
1267	evap_rst(:,:) = evap(:,:)
1268	rugos_rst(:,:) = rugos(:,:)
1269	agesno_rst(:,:) = agesno(:,:)
1270	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
1271
1272	!****************************************************************************************
1273	! Deallocate module variables
1274	!
1275	!****************************************************************************************
1276	DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
1277
1278	END SUBROUTINE pbl_surface_final
1279	!
1280	!****************************************************************************************
1281	!
1282	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, tsurf, alb1, alb2, u10m, v10m, tke)
1283
1284	! Give default values where new fraction has appread
1285
1286	INCLUDE "indicesol.h"
1287	INCLUDE "dimsoil.h"
1288	INCLUDE "clesphys.h"
1289	INCLUDE "compbl.h"
1290
1291	! Input variables
1292	!****************************************************************************************
1293	INTEGER, INTENT(IN) :: itime
1294	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
1295
1296	! InOutput variables
1297	!****************************************************************************************
1298	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
1299	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: alb1, alb2
1300	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: u10m, v10m
1301	REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke
1302
1303	! Local variables
1304	!****************************************************************************************
1305	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
1306	CHARACTER(len=80) :: abort_message
1307	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
1308	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
1309	!
1310	! All at once !!
1311	!****************************************************************************************
1312
1313	DO nsrf = 1, nbsrf
1314	! First decide complement sub-surfaces
1315	SELECT CASE (nsrf)
1316	CASE(is_oce)
1317	nsrf_comp1=is_sic
1318	nsrf_comp2=is_ter
1319	nsrf_comp3=is_lic
1320	CASE(is_sic)
1321	nsrf_comp1=is_oce
1322	nsrf_comp2=is_ter
1323	nsrf_comp3=is_lic
1324	CASE(is_ter)
1325	nsrf_comp1=is_lic
1326	nsrf_comp2=is_oce
1327	nsrf_comp3=is_sic
1328	CASE(is_lic)
1329	nsrf_comp1=is_ter
1330	nsrf_comp2=is_oce
1331	nsrf_comp3=is_sic
1332	END SELECT
1333
1334	! Initialize all new fractions
1335	DO i=1, klon
1336	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
1337
1338	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
1339	! Use the complement sub-surface, keeping the continents unchanged
1340	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
1341	evap(i,nsrf) = evap(i,nsrf_comp1)
1342	rugos(i,nsrf) = rugos(i,nsrf_comp1)
1343	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
1344	alb1(i,nsrf) = alb1(i,nsrf_comp1)
1345	alb2(i,nsrf) = alb2(i,nsrf_comp1)
1346	u10m(i,nsrf) = u10m(i,nsrf_comp1)
1347	v10m(i,nsrf) = v10m(i,nsrf_comp1)
1348	if (iflag_pbl > 1) then
1349	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
1350	endif
1351	mfois(nsrf) = mfois(nsrf) + 1
1352	ELSE
1353	! The continents have changed. The new fraction receives the mean sum of the existent fractions
1354	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1355	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1356	rugos(i,nsrf) = rugos(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + rugos(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1357	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1358	alb1(i,nsrf) = alb1(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + alb1(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1359	alb2(i,nsrf) = alb2(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + alb2(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1360	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1361	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1362	if (iflag_pbl > 1) then
1363	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1364	endif
1365
1366	! Security abort. This option has never been tested. To test, comment the following line.
1367	! abort_message='The fraction of the continents have changed!'
1368	! CALL abort_gcm(modname,abort_message,1)
1369	nfois(nsrf) = nfois(nsrf) + 1
1370	END IF
1371	snow(i,nsrf) = 0.
1372	agesno(i,nsrf) = 0.
1373	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
1374	ELSE
1375	pfois(nsrf) = pfois(nsrf)+ 1
1376	END IF
1377	END DO
1378
1379	END DO
1380
1381	END SUBROUTINE pbl_surface_newfrac
1382
1383	!
1384	!****************************************************************************************
1385	!
1386
1387	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format