Context Navigation

pbl_surface_mod.F90 @ 2056

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

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

Download in other formats:

Original Format