Context Navigation

pbl_surface_mod.F90 @ 1635

Last change on this file since 1635 was 1555, checked in by jghattas, 13 years ago
Ajout du calcul de la temperature sol pour le 1D si flux forces(Nicolas Rochetin).
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 59.1 KB

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

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

Download in other formats:

Original Format