Context Navigation

pbl_surface_mod.F90 @ 1398

Last change on this file since 1398 was 1282, checked in by Ehouarn Millour, 14 years ago

Bug fix: variables cannot be both 'save' and 'common'

(this error is tolerated by most compilers;

but not xlf)

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 58.2 KB

Line
1	!
2	! $Id: pbl_surface_mod.F90 1282 2009-12-16 11:23:05Z musat $
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 ! initialized during first_call below
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	! Initialize ok_flux_surf (for 1D model)
485	ok_flux_surf=.FALSE.
486
487	! Initilize debug IO
488	IF (debugindex .AND. mpi_size==1) THEN
489	! initialize IOIPSL output
490	idayref = day_ini
491	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
492	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
493	DO i = 1, iim
494	zx_lon(i,1) = rlon(i+1)
495	zx_lon(i,jjm+1) = rlon(i+1)
496	ENDDO
497	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
498	CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
499	1,iim,1,jjm+1, &
500	itau_phy,zjulian,dtime,nhoridbg,nidbg)
501	! no vertical axis
502	cl_surf(1)='ter'
503	cl_surf(2)='lic'
504	cl_surf(3)='oce'
505	cl_surf(4)='sic'
506	DO nsrf=1,nbsrf
507	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
508	jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
509	END DO
510
511	CALL histend(nidbg)
512	CALL histsync(nidbg)
513
514	END IF
515
516	ENDIF
517
518	!****************************************************************************************
519	! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
520	! instead of ORCHIDEE)
521	IF (qsol0>0.) THEN
522	PRINT*,'WARNING : On impose qsol=',qsol0
523	qsol(:)=qsol0
524	ENDIF
525	!****************************************************************************************
526
527	!****************************************************************************************
528	! 2) Initialization to zero
529	! Done for all local variables that will be compressed later
530	! and argument with INTENT(OUT)
531	!****************************************************************************************
532	cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0
533	ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0
534	zv1 = 0.0 ; yqsurf = 0.0 ; yalb1 = 0.0 ; yalb2 = 0.0
535	yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0
536	ysollw = 0.0 ; yrugos = 0.0 ; yu1 = 0.0
537	yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0
538	ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0
539	yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
540	yrugoro = 0.0 ; ywindsp = 0.0
541	d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0
542	flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0
543	d_u = 0.0 ; d_v = 0.0 ; yqsol = 0.0
544	ytherm = 0.0 ; ytke=0.
545
546	zcoefh(:,:) = 0.0
547	zcoefh(:,1) = 999999. ! zcoefh(:,k=1) should never be used
548	ytsoil = 999999.
549
550	rh2m(:) = 0.
551	qsat2m(:) = 0.
552	!****************************************************************************************
553	! 3) - Calculate pressure thickness of each layer
554	! - Calculate the wind at first layer
555	! - Mean calculations of albedo
556	! - Calculate net radiance at sub-surface
557	!****************************************************************************************
558	DO k = 1, klev
559	DO i = 1, klon
560	delp(i,k) = paprs(i,k)-paprs(i,k+1)
561	ENDDO
562	ENDDO
563
564	!****************************************************************************************
565	! Test for rugos........ from physiq.. A la fin plutot???
566	!
567	!****************************************************************************************
568
569	zxrugs(:) = 0.0
570	DO nsrf = 1, nbsrf
571	DO i = 1, klon
572	rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015)
573	zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf)
574	ENDDO
575	ENDDO
576
577	! Mean calculations of albedo
578	!
579	! Albedo at sub-surface
580	! * alb1 : albedo in visible SW interval
581	! * alb2 : albedo in near infrared SW interval
582	! * alb : mean albedo for whole SW interval
583	!
584	! Mean albedo for grid point
585	! * alb1_m : albedo in visible SW interval
586	! * alb2_m : albedo in near infrared SW interval
587	! * alb_m : mean albedo at whole SW interval
588
589	alb1_m(:) = 0.0
590	alb2_m(:) = 0.0
591	DO nsrf = 1, nbsrf
592	DO i = 1, klon
593	alb1_m(i) = alb1_m(i) + alb1(i,nsrf) * pctsrf(i,nsrf)
594	alb2_m(i) = alb2_m(i) + alb2(i,nsrf) * pctsrf(i,nsrf)
595	ENDDO
596	ENDDO
597
598	! We here suppose the fraction f1 of incoming radiance of visible radiance
599	! as a fraction of all shortwave radiance
600	f1 = 0.5
601	! f1 = 1 ! put f1=1 to recreate old calculations
602
603	DO nsrf = 1, nbsrf
604	DO i = 1, klon
605	alb(i,nsrf) = f1alb1(i,nsrf) + (1-f1)alb2(i,nsrf)
606	ENDDO
607	ENDDO
608
609	DO i = 1, klon
610	alb_m(i) = f1alb1_m(i) + (1-f1)alb2_m(i)
611	END DO
612
613	! Calculation of mean temperature at surface grid points
614	ztsol(:) = 0.0
615	DO nsrf = 1, nbsrf
616	DO i = 1, klon
617	ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
618	ENDDO
619	ENDDO
620
621	! Linear distrubution on sub-surface of long- and shortwave net radiance
622	DO nsrf = 1, nbsrf
623	DO i = 1, klon
624	sollw(i,nsrf) = sollw_m(i) + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ts(i,nsrf))
625	solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
626	ENDDO
627	ENDDO
628
629
630	! Downwelling longwave radiation at mean surface
631	lwdown_m(:) = 0.0
632	DO i = 1, klon
633	lwdown_m(i) = sollw_m(i) + RSIGMAztsol(i)*4
634	ENDDO
635
636	!****************************************************************************************
637	! 4) Loop over different surfaces
638	!
639	! Only points containing a fraction of the sub surface will be threated.
640	!
641	!****************************************************************************************
642
643	loop_nbsrf: DO nsrf = 1, nbsrf
644
645	! Search for index(ni) and size(knon) of domaine to treat
646	ni(:) = 0
647	knon = 0
648	DO i = 1, klon
649	IF (pctsrf(i,nsrf) > 0.) THEN
650	knon = knon + 1
651	ni(knon) = i
652	ENDIF
653	ENDDO
654
655	! write index, with IOIPSL
656	IF (debugindex .AND. mpi_size==1) THEN
657	tabindx(:)=0.
658	DO i=1,knon
659	tabindx(i)=FLOAT(i)
660	END DO
661	debugtab(:,:) = 0.
662	ndexbg(:) = 0
663	CALL gath2cpl(tabindx,debugtab,knon,ni)
664	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg)
665	ENDIF
666
667	!****************************************************************************************
668	! 5) Compress variables
669	!
670	!****************************************************************************************
671
672	DO j = 1, knon
673	i = ni(j)
674	ypct(j) = pctsrf(i,nsrf)
675	yts(j) = ts(i,nsrf)
676	ysnow(j) = snow(i,nsrf)
677	yqsurf(j) = qsurf(i,nsrf)
678	yalb(j) = alb(i,nsrf)
679	yalb1(j) = alb1(i,nsrf)
680	yalb2(j) = alb2(i,nsrf)
681	yrain_f(j) = rain_f(i)
682	ysnow_f(j) = snow_f(i)
683	yagesno(j) = agesno(i,nsrf)
684	yfder(j) = fder(i)
685	ysolsw(j) = solsw(i,nsrf)
686	ysollw(j) = sollw(i,nsrf)
687	yrugos(j) = rugos(i,nsrf)
688	yrugoro(j) = rugoro(i)
689	yu1(j) = u(i,1)
690	yv1(j) = v(i,1)
691	ypaprs(j,klev+1) = paprs(i,klev+1)
692	ywindsp(j) = SQRT(u10m(i,nsrf)2 + v10m(i,nsrf)2 )
693	END DO
694
695	DO k = 1, klev
696	DO j = 1, knon
697	i = ni(j)
698	ypaprs(j,k) = paprs(i,k)
699	ypplay(j,k) = pplay(i,k)
700	ydelp(j,k) = delp(i,k)
701	ytke(j,k) = tke(i,k,nsrf)
702	yu(j,k) = u(i,k)
703	yv(j,k) = v(i,k)
704	yt(j,k) = t(i,k)
705	yq(j,k) = q(i,k)
706	ENDDO
707	ENDDO
708
709	DO k = 1, nsoilmx
710	DO j = 1, knon
711	i = ni(j)
712	ytsoil(j,k) = ftsoil(i,k,nsrf)
713	END DO
714	END DO
715
716	! qsol(water height in soil) only for bucket continental model
717	IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
718	DO j = 1, knon
719	i = ni(j)
720	yqsol(j) = qsol(i)
721	END DO
722	ENDIF
723
724	!****************************************************************************************
725	! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
726	!
727	!****************************************************************************************
728
729	CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
730	yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
731	yts, yqsurf, yrugos, &
732	ycdragm, ycdragh )
733
734	!****************************************************************************************
735	! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefm et ycoefm.
736	!
737	!****************************************************************************************
738
739	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
740	ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
741	ycoefm, ycoefh, ytke)
742
743	!****************************************************************************************
744	!
745	! 8) "La descente" - "The downhill"
746	!
747	! climb_hq_down and climb_wind_down calculate the coefficients
748	! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
749	! Only the coefficients at surface for H and Q are returned.
750	!
751	!****************************************************************************************
752
753	! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
754	CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
755	ydelp, yt, yq, dtime, &
756	AcoefH, AcoefQ, BcoefH, BcoefQ)
757
758	! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
759	CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
760	AcoefU, AcoefV, BcoefU, BcoefV)
761
762
763	!****************************************************************************************
764	! 9) Small calculations
765	!
766	!****************************************************************************************
767
768	! - Reference pressure is given the values at surface level
769	ypsref(:) = ypaprs(:,1)
770
771	! - CO2 field on 2D grid to be sent to ORCHIDEE
772	! Transform to compressed field
773	IF (carbon_cycle_cpl) THEN
774	DO i=1,knon
775	r_co2_ppm(i) = co2_send(ni(i))
776	END DO
777	ELSE
778	r_co2_ppm(:) = co2_ppm ! Constant field
779	END IF
780
781	!****************************************************************************************
782	!
783	! Calulate t2m and q2m for the case of calculation at land grid points
784	! t2m and q2m are needed as input to ORCHIDEE
785	!
786	!****************************************************************************************
787	IF (nsrf == is_ter) THEN
788
789	DO i = 1, knon
790	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
791	* (ypaprs(i,1)-ypplay(i,1))
792	END DO
793
794	! Calculate the temperature et relative humidity at 2m and the wind at 10m
795	CALL stdlevvar(klon, knon, is_ter, zxli, &
796	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
797	yts, yqsurf, yrugos, ypaprs(:,1), ypplay(:,1), &
798	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
799
800	END IF
801
802	!****************************************************************************************
803	!
804	! 10) Switch selon current surface
805	! It is necessary to start with the continental surfaces because the ocean
806	! needs their run-off.
807	!
808	!****************************************************************************************
809	SELECT CASE(nsrf)
810
811	CASE(is_ter)
812	! ylwdown : to be removed, calculation is now done at land surface in surf_land
813	ylwdown(:)=0.0
814	DO i=1,knon
815	ylwdown(i)=lwdown_m(ni(i))
816	END DO
817	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
818	rlon, rlat, &
819	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
820	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
821	AcoefH, AcoefQ, BcoefH, BcoefQ, &
822	AcoefU, AcoefV, BcoefU, BcoefV, &
823	ypsref, yu1, yv1, yrugoro, pctsrf, &
824	ylwdown, yq2m, yt2m, &
825	ysnow, yqsol, yagesno, ytsoil, &
826	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
827	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
828	y_flux_u1, y_flux_v1 )
829
830
831	CASE(is_lic)
832	CALL surf_landice(itap, dtime, knon, ni, &
833	ysolsw, ysollw, yts, ypplay(:,1), &
834	ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
835	AcoefH, AcoefQ, BcoefH, BcoefQ, &
836	AcoefU, AcoefV, BcoefU, BcoefV, &
837	ypsref, yu1, yv1, yrugoro, pctsrf, &
838	ysnow, yqsurf, yqsol, yagesno, &
839	ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
840	ytsurf_new, y_dflux_t, y_dflux_q, &
841	y_flux_u1, y_flux_v1)
842
843	CASE(is_oce)
844	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb1, &
845	yrugos, ywindsp, rmu0, yfder, yts, &
846	itap, dtime, jour, knon, ni, &
847	ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
848	AcoefH, AcoefQ, BcoefH, BcoefQ, &
849	AcoefU, AcoefV, BcoefU, BcoefV, &
850	ypsref, yu1, yv1, yrugoro, pctsrf, &
851	ysnow, yqsurf, yagesno, &
852	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
853	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
854	y_flux_u1, y_flux_v1)
855
856	CASE(is_sic)
857	CALL surf_seaice( &
858	rlon, rlat, ysolsw, ysollw, yalb1, yfder, &
859	itap, dtime, jour, knon, ni, &
860	lafin, &
861	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
862	AcoefH, AcoefQ, BcoefH, BcoefQ, &
863	AcoefU, AcoefV, BcoefU, BcoefV, &
864	ypsref, yu1, yv1, yrugoro, pctsrf, &
865	ysnow, yqsurf, yqsol, yagesno, ytsoil, &
866	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
867	ytsurf_new, y_dflux_t, y_dflux_q, &
868	y_flux_u1, y_flux_v1)
869
870
871	CASE DEFAULT
872	WRITE(lunout,*) 'Surface index = ', nsrf
873	abort_message = 'Surface index not valid'
874	CALL abort_gcm(modname,abort_message,1)
875	END SELECT
876
877
878	!****************************************************************************************
879	! 11) - Calcul the increment of surface temperature
880	!
881	!****************************************************************************************
882	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
883
884	!****************************************************************************************
885	!
886	! 12) "La remontee" - "The uphill"
887	!
888	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
889	! for X=H, Q, U and V, for all vertical levels.
890	!
891	!****************************************************************************************
892	! H and Q
893	IF (ok_flux_surf) THEN
894	PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
895	y_flux_t1(:) = fsens
896	y_flux_q1(:) = flat/RLVTT
897	yfluxlat(:) = flat
898	ELSE
899	y_flux_t1(:) = yfluxsens(:)
900	y_flux_q1(:) = -yevap(:)
901	ENDIF
902
903	CALL climb_hq_up(knon, dtime, yt, yq, &
904	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
905	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
906
907
908	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
909	y_flux_u, y_flux_v, y_d_u, y_d_v)
910
911
912	DO j = 1, knon
913	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
914	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
915	ENDDO
916
917	!****************************************************************************************
918	! 13) Transform variables for output format :
919	! - Decompress
920	! - Multiply with pourcentage of current surface
921	! - Cumulate in global variable
922	!
923	!****************************************************************************************
924
925	tke(:,:,nsrf) = 0.
926	DO k = 1, klev
927	DO j = 1, knon
928	i = ni(j)
929	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
930	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
931	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
932	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
933
934	flux_t(i,k,nsrf) = y_flux_t(j,k)
935	flux_q(i,k,nsrf) = y_flux_q(j,k)
936	flux_u(i,k,nsrf) = y_flux_u(j,k)
937	flux_v(i,k,nsrf) = y_flux_v(j,k)
938
939	tke(i,k,nsrf) = ytke(j,k)
940
941	ENDDO
942	ENDDO
943
944	evap(:,nsrf) = - flux_q(:,1,nsrf)
945
946	alb1(:, nsrf) = 0.
947	alb2(:, nsrf) = 0.
948	snow(:, nsrf) = 0.
949	qsurf(:, nsrf) = 0.
950	rugos(:, nsrf) = 0.
951	fluxlat(:,nsrf) = 0.
952	DO j = 1, knon
953	i = ni(j)
954	d_ts(i,nsrf) = y_d_ts(j)
955	alb1(i,nsrf) = yalb1_new(j)
956	alb2(i,nsrf) = yalb2_new(j)
957	snow(i,nsrf) = ysnow(j)
958	qsurf(i,nsrf) = yqsurf(j)
959	rugos(i,nsrf) = yz0_new(j)
960	fluxlat(i,nsrf) = yfluxlat(j)
961	agesno(i,nsrf) = yagesno(j)
962	cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
963	cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
964	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
965	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
966	END DO
967
968	DO k = 2, klev
969	DO j = 1, knon
970	i = ni(j)
971	zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)*ypct(j)
972	END DO
973	END DO
974
975	IF ( nsrf .EQ. is_ter ) THEN
976	DO j = 1, knon
977	i = ni(j)
978	qsol(i) = yqsol(j)
979	END DO
980	END IF
981
982	ftsoil(:,:,nsrf) = 0.
983	DO k = 1, nsoilmx
984	DO j = 1, knon
985	i = ni(j)
986	ftsoil(i, k, nsrf) = ytsoil(j,k)
987	END DO
988	END DO
989
990
991	DO k = 1, klev
992	DO j = 1, knon
993	i = ni(j)
994	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
995	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
996	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
997	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
998	END DO
999	END DO
1000
1001	!****************************************************************************************
1002	! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m
1003	! Call HBTM
1004	!
1005	!****************************************************************************************
1006	t2m(:,nsrf) = 0.
1007	q2m(:,nsrf) = 0.
1008	u10m(:,nsrf) = 0.
1009	v10m(:,nsrf) = 0.
1010
1011	pblh(:,nsrf) = 0. ! Hauteur de couche limite
1012	plcl(:,nsrf) = 0. ! Niveau de condensation de la CLA
1013	capCL(:,nsrf) = 0. ! CAPE de couche limite
1014	oliqCL(:,nsrf) = 0. ! eau_liqu integree de couche limite
1015	cteiCL(:,nsrf) = 0. ! cloud top instab. crit. couche limite
1016	pblt(:,nsrf) = 0. ! T a la Hauteur de couche limite
1017	therm(:,nsrf) = 0.
1018	trmb1(:,nsrf) = 0. ! deep_cape
1019	trmb2(:,nsrf) = 0. ! inhibition
1020	trmb3(:,nsrf) = 0. ! Point Omega
1021
1022	#undef T2m
1023	#define T2m
1024	#ifdef T2m
1025	! Calculations of diagnostic t,q at 2m and u, v at 10m
1026
1027	DO j=1, knon
1028	i = ni(j)
1029	uzon(j) = yu(j,1) + y_d_u(j,1)
1030	vmer(j) = yv(j,1) + y_d_v(j,1)
1031	tair1(j) = yt(j,1) + y_d_t(j,1)
1032	qair1(j) = yq(j,1) + y_d_q(j,1)
1033	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
1034	* (ypaprs(j,1)-ypplay(j,1))
1035	tairsol(j) = yts(j) + y_d_ts(j)
1036	rugo1(j) = yrugos(j)
1037	IF(nsrf.EQ.is_oce) THEN
1038	rugo1(j) = rugos(i,nsrf)
1039	ENDIF
1040	psfce(j)=ypaprs(j,1)
1041	patm(j)=ypplay(j,1)
1042	qairsol(j) = yqsurf(j)
1043	END DO
1044
1045
1046	! Calculate the temperature et relative humidity at 2m and the wind at 10m
1047	CALL stdlevvar(klon, knon, nsrf, zxli, &
1048	uzon, vmer, tair1, qair1, zgeo1, &
1049	tairsol, qairsol, rugo1, psfce, patm, &
1050	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
1051
1052	DO j=1, knon
1053	i = ni(j)
1054	t2m(i,nsrf)=yt2m(j)
1055	q2m(i,nsrf)=yq2m(j)
1056
1057	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
1058	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
1059	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
1060	END DO
1061
1062	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
1063	!IM Ajoute dependance type surface
1064	IF (thermcep) THEN
1065	DO j = 1, knon
1066	i=ni(j)
1067	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
1068	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
1069	zx_qs1 = MIN(0.5,zx_qs1)
1070	zcor1 = 1./(1.-RETV*zx_qs1)
1071	zx_qs1 = zx_qs1*zcor1
1072
1073	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
1074	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
1075	END DO
1076	END IF
1077
1078	CALL HBTM(knon, ypaprs, ypplay, &
1079	yt2m,yt10m,yq2m,yq10m,yustar, &
1080	y_flux_t,y_flux_q,yu,yv,yt,yq, &
1081	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
1082	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
1083
1084	DO j=1, knon
1085	i = ni(j)
1086	pblh(i,nsrf) = ypblh(j)
1087	plcl(i,nsrf) = ylcl(j)
1088	capCL(i,nsrf) = ycapCL(j)
1089	oliqCL(i,nsrf) = yoliqCL(j)
1090	cteiCL(i,nsrf) = ycteiCL(j)
1091	pblT(i,nsrf) = ypblT(j)
1092	therm(i,nsrf) = ytherm(j)
1093	trmb1(i,nsrf) = ytrmb1(j)
1094	trmb2(i,nsrf) = ytrmb2(j)
1095	trmb3(i,nsrf) = ytrmb3(j)
1096	END DO
1097
1098	#else
1099	! T2m not defined
1100	! No calculation
1101	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
1102	#endif
1103
1104	!****************************************************************************************
1105	! 15) End of loop over different surfaces
1106	!
1107	!****************************************************************************************
1108	END DO loop_nbsrf
1109
1110	!****************************************************************************************
1111	! 16) Calculate the mean value over all sub-surfaces for som variables
1112	!
1113	!****************************************************************************************
1114
1115	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
1116	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
1117	DO nsrf = 1, nbsrf
1118	DO k = 1, klev
1119	DO i = 1, klon
1120	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
1121	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
1122	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
1123	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
1124	END DO
1125	END DO
1126	END DO
1127
1128	DO i = 1, klon
1129	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1130	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1131	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
1132	ENDDO
1133
1134	!
1135	! Incrementer la temperature du sol
1136	!
1137	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
1138	zt2m(:) = 0.0 ; zq2m(:) = 0.0
1139	zu10m(:) = 0.0 ; zv10m(:) = 0.0
1140	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
1141	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
1142	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
1143	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
1144	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
1145
1146
1147	DO nsrf = 1, nbsrf
1148	DO i = 1, klon
1149	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
1150
1151	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
1152	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1153	wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
1154	pctsrf(i,nsrf)
1155
1156	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
1157	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
1158
1159	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
1160	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
1161	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
1162	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
1163
1164	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
1165	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
1166	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
1167	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
1168	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
1169	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
1170	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
1171	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
1172	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
1173	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
1174	END DO
1175	END DO
1176
1177	IF (check) THEN
1178	amn=MIN(ts(1,is_ter),1000.)
1179	amx=MAX(ts(1,is_ter),-1000.)
1180	DO i=2, klon
1181	amn=MIN(ts(i,is_ter),amn)
1182	amx=MAX(ts(i,is_ter),amx)
1183	ENDDO
1184	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
1185	ENDIF
1186
1187	!jg ?
1188	!!$!
1189	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
1190	!!$! sub-surfaces is distributed.
1191	!!$!
1192	!!$ DO nsrf = 1, nbsrf
1193	!!$ DO i = 1, klon
1194	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
1195	!!$ ts(i,nsrf) = zxtsol(i)
1196	!!$ t2m(i,nsrf) = zt2m(i)
1197	!!$ q2m(i,nsrf) = zq2m(i)
1198	!!$ u10m(i,nsrf) = zu10m(i)
1199	!!$ v10m(i,nsrf) = zv10m(i)
1200	!!$
1201	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
1202	!!$ pblh(i,nsrf) = s_pblh(i)
1203	!!$ plcl(i,nsrf) = s_plcl(i)
1204	!!$ capCL(i,nsrf) = s_capCL(i)
1205	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
1206	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
1207	!!$ pblT(i,nsrf) = s_pblT(i)
1208	!!$ therm(i,nsrf) = s_therm(i)
1209	!!$ trmb1(i,nsrf) = s_trmb1(i)
1210	!!$ trmb2(i,nsrf) = s_trmb2(i)
1211	!!$ trmb3(i,nsrf) = s_trmb3(i)
1212	!!$ ENDIF
1213	!!$ ENDDO
1214	!!$ ENDDO
1215
1216
1217	DO i = 1, klon
1218	fder(i) = - 4.0RSIGMAzxtsol(i)**3
1219	ENDDO
1220
1221	zxqsurf(:) = 0.0
1222	zxsnow(:) = 0.0
1223	DO nsrf = 1, nbsrf
1224	DO i = 1, klon
1225	zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
1226	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
1227	END DO
1228	END DO
1229
1230	! Premier niveau de vent sortie dans physiq.F
1231	zu1(:) = u(:,1)
1232	zv1(:) = v(:,1)
1233
1234	! Some of the module declared variables are returned for printing in physiq.F
1235	qsol_d(:) = qsol(:)
1236	evap_d(:,:) = evap(:,:)
1237	rugos_d(:,:) = rugos(:,:)
1238	agesno_d(:,:) = agesno(:,:)
1239
1240
1241	END SUBROUTINE pbl_surface
1242	!
1243	!****************************************************************************************
1244	!
1245	SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
1246	evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
1247
1248	INCLUDE "indicesol.h"
1249	INCLUDE "dimsoil.h"
1250
1251	! Ouput variables
1252	!****************************************************************************************
1253	REAL, DIMENSION(klon), INTENT(OUT) :: qsol_rst
1254	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
1255	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
1256	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
1257	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_rst
1258	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_rst
1259	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_rst
1260	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
1261
1262
1263	!****************************************************************************************
1264	! Return module variables for writing to restart file
1265	!
1266	!****************************************************************************************
1267	qsol_rst(:) = qsol(:)
1268	fder_rst(:) = fder(:)
1269	snow_rst(:,:) = snow(:,:)
1270	qsurf_rst(:,:) = qsurf(:,:)
1271	evap_rst(:,:) = evap(:,:)
1272	rugos_rst(:,:) = rugos(:,:)
1273	agesno_rst(:,:) = agesno(:,:)
1274	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
1275
1276	!****************************************************************************************
1277	! Deallocate module variables
1278	!
1279	!****************************************************************************************
1280	DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
1281
1282	END SUBROUTINE pbl_surface_final
1283	!
1284	!****************************************************************************************
1285	!
1286	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, tsurf, alb1, alb2, u10m, v10m, tke)
1287
1288	! Give default values where new fraction has appread
1289
1290	INCLUDE "indicesol.h"
1291	INCLUDE "dimsoil.h"
1292	INCLUDE "clesphys.h"
1293	INCLUDE "compbl.h"
1294
1295	! Input variables
1296	!****************************************************************************************
1297	INTEGER, INTENT(IN) :: itime
1298	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
1299
1300	! InOutput variables
1301	!****************************************************************************************
1302	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
1303	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: alb1, alb2
1304	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: u10m, v10m
1305	REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke
1306
1307	! Local variables
1308	!****************************************************************************************
1309	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
1310	CHARACTER(len=80) :: abort_message
1311	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
1312	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
1313	!
1314	! All at once !!
1315	!****************************************************************************************
1316
1317	DO nsrf = 1, nbsrf
1318	! First decide complement sub-surfaces
1319	SELECT CASE (nsrf)
1320	CASE(is_oce)
1321	nsrf_comp1=is_sic
1322	nsrf_comp2=is_ter
1323	nsrf_comp3=is_lic
1324	CASE(is_sic)
1325	nsrf_comp1=is_oce
1326	nsrf_comp2=is_ter
1327	nsrf_comp3=is_lic
1328	CASE(is_ter)
1329	nsrf_comp1=is_lic
1330	nsrf_comp2=is_oce
1331	nsrf_comp3=is_sic
1332	CASE(is_lic)
1333	nsrf_comp1=is_ter
1334	nsrf_comp2=is_oce
1335	nsrf_comp3=is_sic
1336	END SELECT
1337
1338	! Initialize all new fractions
1339	DO i=1, klon
1340	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
1341
1342	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
1343	! Use the complement sub-surface, keeping the continents unchanged
1344	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
1345	evap(i,nsrf) = evap(i,nsrf_comp1)
1346	rugos(i,nsrf) = rugos(i,nsrf_comp1)
1347	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
1348	alb1(i,nsrf) = alb1(i,nsrf_comp1)
1349	alb2(i,nsrf) = alb2(i,nsrf_comp1)
1350	u10m(i,nsrf) = u10m(i,nsrf_comp1)
1351	v10m(i,nsrf) = v10m(i,nsrf_comp1)
1352	if (iflag_pbl > 1) then
1353	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
1354	endif
1355	mfois(nsrf) = mfois(nsrf) + 1
1356	ELSE
1357	! The continents have changed. The new fraction receives the mean sum of the existent fractions
1358	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1359	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1360	rugos(i,nsrf) = rugos(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + rugos(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1361	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1362	alb1(i,nsrf) = alb1(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + alb1(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1363	alb2(i,nsrf) = alb2(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + alb2(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1364	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1365	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1366	if (iflag_pbl > 1) then
1367	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1368	endif
1369
1370	! Security abort. This option has never been tested. To test, comment the following line.
1371	! abort_message='The fraction of the continents have changed!'
1372	! CALL abort_gcm(modname,abort_message,1)
1373	nfois(nsrf) = nfois(nsrf) + 1
1374	END IF
1375	snow(i,nsrf) = 0.
1376	agesno(i,nsrf) = 0.
1377	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
1378	ELSE
1379	pfois(nsrf) = pfois(nsrf)+ 1
1380	END IF
1381	END DO
1382
1383	END DO
1384
1385	END SUBROUTINE pbl_surface_newfrac
1386
1387	!
1388	!****************************************************************************************
1389	!
1390
1391	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format