Context Navigation

pbl_surface_mod.F90 @ 5170

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

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

Download in other formats:

Original Format