Context Navigation

pbl_surface_mod.F90 @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 58.6 KB

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

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

Download in other formats:

Original Format