Context Navigation

pbl_surface_mod.F90 @ 183

Last change on this file since 183 was 183, checked in by lfita, 10 years ago
Removing check pritings Removing nullification of MixingRatioValues?(:,:,2) = 0 and MixingTendRatioValues?(:,:,2) NO plul related to te absence of oliq at the initial conditions. If one introduce fake ones it works!
Property svn:executable set to ``*
File size: 62.3 KB

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

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

Download in other formats:

Original Format