Context Navigation

pbl_surface_mod.F90 @ 1707

Last change on this file since 1707 was 1707, checked in by Laurent Fairhead, 11 years ago

Version testing basée sur la r1706

Testing release based on r1706

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

File size: 59.6 KB

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

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

Download in other formats:

Original Format