Context Navigation

pbl_surface_mod.F90 @ 1763

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

New version of Mellor et Yamada pronostic TKE

... based on energy transfer from the mean state.

The new version is yamada_c.
It must be called after vertical diffusion rather than just before
since the source terms u_z w'u' + v_z w'theta' and g/theta w'theta'
are diagnosed from the vertical diffusion (as energy loss from the mean
state) rather than computed as K (u_z^2+v_z2) or g/\theta K theta_z
(where _z means vertical derivative).
The call to this version is controled by iflag_pbl.

iflag_pbl = 20 : with TKE computed at interfaces between layers
iflag_pbl = 25 : with TKE computed within the layer
In both cases, the dissipation of turbulence is translated into heat, and
passed to the physics as dtdiss (temperature tendency due to dissipation

of TKE).

The diffusion coefficient being computed after dissipation, it must be
kept for diffusion at the next time step, and thus be stored in the
restartphy file.
This coefficient must be computed and stored for each sub-surface.

A new way of managing subsurface variables is introduced.
For arrays of the form X(:,nbsrf) are extented to X(:,nbsrf+1), where
is_ave=nbsrf+1, is an additional sub-surface which contains the averaged values.

coef_diff_turb_mod.F90 : change of flags.
ener_conserv.F90 : energy conservation must not be applied in those

cases

indicesol.h : definition of is_ave
pbl_surface_mod.F90 : call to yamada_c and changes in the management of

coefh/coefm

physiq.F : Change in the initialisation of pmflxr/s and

modified calls to pbl_surface_mod (introduction
of dtdiss, initialisation
of pbl_tke and coefh in 1D).

phys_local_var_mod.F90 : declaration of d_t_diss
phys_output_mod.F90 : new outputs (bils_tke and bils_diss) and

coefh->coefh(:,:,is_ave)

phys_state_var_mod.F90 : modified declaration for coefh and coefm

(nbsrf -> nbsrf+1)

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

File size: 61.4 KB

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

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

Download in other formats:

Original Format