Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

pbl_surface_mod.F90 @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

File size: 65.2 KB

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

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

Download in other formats:

Original Format