Context Navigation

pbl_surface_mod.F90 @ 2261

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

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

Download in other formats:

Original Format