Context Navigation

pbl_surface_mod.F90 @ 1865

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

Inclusion de la bibliothèque SISVAT/MAR à LMDZ pour le traitement des surfaces
"land ice"

Menegoz

Integration of the SISVAT/MAR library to LMDZ to model the land ice surfaces

Menegoz

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

File size: 64.9 KB

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

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

Download in other formats:

Original Format