Context Navigation

pbl_surface_mod.F90 @ 2461

Last change on this file since 2461 was 2160, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 65.8 KB

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

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

Download in other formats:

Original Format