Context Navigation

pbl_surface_mod.F90 @ 1785

Last change on this file since 1785 was 1785, checked in by Ehouarn Millour, 11 years ago

Transformation de l'include indicesol.h en un module indice_sol_mod et modification des appels dans tous les fichiers concernés.
Aucun changement des résultats ni des sorties du modèle vs 1784.
UG

...................................................

Replacement of the indicesol.h include by a module named indice_sol_mod. Modification of the calls in every affected files.
Results and outputs of simulations are unchanged in comparison with rev 1784.
UG

Property svn:executable set to *

File size: 61.5 KB

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

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

Download in other formats:

Original Format