Context Navigation

pbl_surface_mod.F90 @ 5416

Last change on this file since 5416 was 1942, checked in by jghattas, 11 years ago
Passing the cosine of solar zenith angle from pbl_surface into orchidee. Matthew McGrath?(LSCE)
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 61.5 KB

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

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

Download in other formats:

Original Format