Context Navigation

pbl_surface_mod.F90 @ 1787

Last change on this file since 1787 was 1787, checked in by idelkadi, 11 years ago
Correction pour prendre en compte le cas AMMA
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 61.4 KB

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

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

Download in other formats:

Original Format