Context Navigation

pbl_surface_mod.F90 @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 64.6 KB

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

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

Download in other formats:

Original Format