Context Navigation

pbl_surface_mod.F90 @ 1850

Last change on this file since 1850 was 1816, checked in by idelkadi, 11 years ago
Inclusion of the convective scale velocity w* for tracers Concerns : hbtm.F, pbl_surface_mod.F90, physiq.F, phys_output_ctrlout_mod.F90, phys_output_write_mod.F90
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 61.7 KB

Line
1	!
2	! $Id: pbl_surface_mod.F90 1816 2013-07-25 08:33:44Z idelkadi $
3	!
4	MODULE pbl_surface_mod
5	!
6	! Planetary Boundary Layer and Surface module
7	!
8	! This module manage the calculation of turbulent diffusion in the boundary layer
9	! and all interactions towards the differents sub-surfaces.
10	!
11	!
12	USE dimphy
13	USE mod_phys_lmdz_para, ONLY : mpi_size
14	USE ioipsl
15	USE surface_data, ONLY : type_ocean, ok_veget
16	USE surf_land_mod, ONLY : surf_land
17	USE surf_landice_mod, ONLY : surf_landice
18	USE surf_ocean_mod, ONLY : surf_ocean
19	USE surf_seaice_mod, ONLY : surf_seaice
20	USE cpl_mod, ONLY : gath2cpl
21	USE climb_hq_mod, ONLY : climb_hq_down, climb_hq_up
22	USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up
23	USE coef_diff_turb_mod, ONLY : coef_diff_turb
24	USE control_mod
25
26
27	IMPLICIT NONE
28
29	! Declaration of variables saved in restart file
30	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: qsol ! water height in the soil (mm)
31	!$OMP THREADPRIVATE(qsol)
32	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: fder ! flux drift
33	!$OMP THREADPRIVATE(fder)
34	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: snow ! snow at surface
35	!$OMP THREADPRIVATE(snow)
36	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qsurf ! humidity at surface
37	!$OMP THREADPRIVATE(qsurf)
38	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: evap ! evaporation at surface
39	!$OMP THREADPRIVATE(evap)
40	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: rugos ! rugosity at surface (m)
41	!$OMP THREADPRIVATE(rugos)
42	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: agesno ! age of snow at surface
43	!$OMP THREADPRIVATE(agesno)
44	! Correction pour le cas AMMA (PRIVATE)
45	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ftsoil ! soil temperature
46	!$OMP THREADPRIVATE(ftsoil)
47
48	CONTAINS
49	!
50	!****************************************************************************************
51	!
52	SUBROUTINE pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst,&
53	evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
54
55	! This routine should be called after the restart file has been read.
56	! This routine initialize the restart variables and does some validation tests
57	! for the index of the different surfaces and tests the choice of type of ocean.
58
59	USE indice_sol_mod
60
61	INCLUDE "dimsoil.h"
62	INCLUDE "iniprint.h"
63
64	! Input variables
65	!****************************************************************************************
66	REAL, DIMENSION(klon), INTENT(IN) :: qsol_rst
67	REAL, DIMENSION(klon), INTENT(IN) :: fder_rst
68	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: snow_rst
69	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: qsurf_rst
70	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: evap_rst
71	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: rugos_rst
72	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: agesno_rst
73	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
74
75
76	! Local variables
77	!****************************************************************************************
78	INTEGER :: ierr
79	CHARACTER(len=80) :: abort_message
80	CHARACTER(len = 20) :: modname = 'pbl_surface_init'
81
82
83	!****************************************************************************************
84	! Allocate and initialize module variables with fields read from restart file.
85	!
86	!****************************************************************************************
87	ALLOCATE(qsol(klon), stat=ierr)
88	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
89
90	ALLOCATE(fder(klon), stat=ierr)
91	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
92
93	ALLOCATE(snow(klon,nbsrf), stat=ierr)
94	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
95
96	ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
97	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
98
99	ALLOCATE(evap(klon,nbsrf), stat=ierr)
100	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
101
102	ALLOCATE(rugos(klon,nbsrf), stat=ierr)
103	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
104
105	ALLOCATE(agesno(klon,nbsrf), stat=ierr)
106	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
107
108	ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
109	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
110
111
112	qsol(:) = qsol_rst(:)
113	fder(:) = fder_rst(:)
114	snow(:,:) = snow_rst(:,:)
115	qsurf(:,:) = qsurf_rst(:,:)
116	evap(:,:) = evap_rst(:,:)
117	rugos(:,:) = rugos_rst(:,:)
118	agesno(:,:) = agesno_rst(:,:)
119	ftsoil(:,:,:) = ftsoil_rst(:,:,:)
120
121
122	!****************************************************************************************
123	! Test for sub-surface indices
124	!
125	!****************************************************************************************
126	IF (is_ter /= 1) THEN
127	WRITE(lunout,)" Warning *"
128	WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
129	WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
130	abort_message="voir ci-dessus"
131	CALL abort_gcm(modname,abort_message,1)
132	ENDIF
133
134	IF ( is_oce > is_sic ) THEN
135	WRITE(lunout,)' Warning *'
136	WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
137	WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
138	WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
139	abort_message='voir ci-dessus'
140	CALL abort_gcm(modname,abort_message,1)
141	ENDIF
142
143	IF ( is_lic > is_sic ) THEN
144	WRITE(lunout,)' Warning *'
145	WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
146	WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
147	WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
148	abort_message='voir ci-dessus'
149	CALL abort_gcm(modname,abort_message,1)
150	ENDIF
151
152	!****************************************************************************************
153	! Validation of ocean mode
154	!
155	!****************************************************************************************
156
157	IF (type_ocean /= 'slab ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN
158	WRITE(lunout,)' Warning *'
159	WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean
160	abort_message='option pour l''ocean non valable'
161	CALL abort_gcm(modname,abort_message,1)
162	ENDIF
163
164	END SUBROUTINE pbl_surface_init
165	!
166	!****************************************************************************************
167	!
168
169	SUBROUTINE pbl_surface( &
170	dtime, date0, itap, jour, &
171	debut, lafin, &
172	rlon, rlat, rugoro, rmu0, &
173	rain_f, snow_f, solsw_m, sollw_m, &
174	t, q, u, v, &
175	pplay, paprs, pctsrf, &
176	ts, alb1, alb2,ustar, u10m, v10m,wstar, &
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+1), INTENT(INOUT) :: wstar ! w* (m/s)
297	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: u10m ! u speed at 10m
298	REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: v10m ! v speed at 10m
299	REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke
300
301	! Output variables
302	!****************************************************************************************
303	REAL, DIMENSION(klon), INTENT(OUT) :: lwdown_m ! Downcoming longwave radiation
304	REAL, DIMENSION(klon), INTENT(OUT) :: cdragh ! drag coefficient for T and Q
305	REAL, DIMENSION(klon), INTENT(OUT) :: cdragm ! drag coefficient for wind
306	REAL, DIMENSION(klon), INTENT(OUT) :: zu1 ! u wind speed in first layer
307	REAL, DIMENSION(klon), INTENT(OUT) :: zv1 ! v wind speed in first layer
308	REAL, DIMENSION(klon), INTENT(OUT) :: alb1_m ! mean albedo in visible SW interval
309	REAL, DIMENSION(klon), INTENT(OUT) :: alb2_m ! mean albedo in near IR SW interval
310	REAL, DIMENSION(klon), INTENT(OUT) :: zxsens ! sensible heat flux at surface with inversed sign
311	! (=> positive sign upwards)
312	REAL, DIMENSION(klon), INTENT(OUT) :: zxevap ! water vapour flux at surface, positiv upwards
313	REAL, DIMENSION(klon), INTENT(OUT) :: zxtsol ! temperature at surface, mean for each grid point
314	REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat ! latent flux, mean for each grid point
315	REAL, DIMENSION(klon), INTENT(OUT) :: zt2m ! temperature at 2m, mean for each grid point
316	REAL, DIMENSION(klon), INTENT(OUT) :: qsat2m
317	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! change in temperature
318	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t_diss ! change in temperature
319	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_q ! change in water vapour
320	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u ! change in u speed
321	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v ! change in v speed
322	REAL, DIMENSION(klon, klev,nbsrf+1), INTENT(OUT) :: zcoefh ! coef for turbulent diffusion of T and Q, mean for each grid point
323	REAL, DIMENSION(klon, klev,nbsrf+1), INTENT(OUT) :: zcoefm ! coef for turbulent diffusion of U and V (?), mean for each grid point
324
325	! Output only for diagnostics
326	REAL, DIMENSION(klon), INTENT(OUT) :: slab_wfbils! heat balance at surface only for slab at ocean points
327	REAL, DIMENSION(klon), INTENT(OUT) :: qsol_d ! water height in the soil (mm)
328	REAL, DIMENSION(klon), INTENT(OUT) :: zq2m ! water vapour at 2m, mean for each grid point
329	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh ! height of the planetary boundary layer(HPBL)
330	REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl ! condensation level
331	REAL, DIMENSION(klon), INTENT(OUT) :: s_capCL ! CAPE of PBL
332	REAL, DIMENSION(klon), INTENT(OUT) :: s_oliqCL ! liquid water intergral of PBL
333	REAL, DIMENSION(klon), INTENT(OUT) :: s_cteiCL ! cloud top instab. crit. of PBL
334	REAL, DIMENSION(klon), INTENT(OUT) :: s_pblT ! temperature at PBLH
335	REAL, DIMENSION(klon), INTENT(OUT) :: s_therm ! thermal virtual temperature excess
336	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb1 ! deep cape, mean for each grid point
337	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb2 ! inhibition, mean for each grid point
338	REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb3 ! point Omega, mean for each grid point
339	REAL, DIMENSION(klon), INTENT(OUT) :: zxrugs ! rugosity at surface (m), mean for each grid point
340	REAL, DIMENSION(klon), INTENT(OUT) :: zustar ! u*
341	REAL, DIMENSION(klon), INTENT(OUT) :: zu10m ! u speed at 10m, mean for each grid point
342	REAL, DIMENSION(klon), INTENT(OUT) :: zv10m ! v speed at 10m, mean for each grid point
343	REAL, DIMENSION(klon), INTENT(OUT) :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
344	REAL, DIMENSION(klon), INTENT(OUT) :: zxqsurf ! humidity at surface, mean for each grid point
345	REAL, DIMENSION(klon), INTENT(OUT) :: rh2m ! relative humidity at 2m
346	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxu ! u wind tension, mean for each grid point
347	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxv ! v wind tension, mean for each grid point
348	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_d ! rugosity length (m)
349	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_d ! age of snow at surface
350	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: solsw ! net shortwave radiation at surface
351	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: sollw ! net longwave radiation at surface
352	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: d_ts ! change in temperature at surface
353	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_d ! evaporation at surface
354	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: fluxlat ! latent flux
355	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: t2m ! temperature at 2 meter height
356	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbils ! heat balance at surface
357	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbilo ! water balance at surface
358	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t ! sensible heat flux (CpT) J/m2/s (W/m2)
359	! positve orientation downwards
360	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u ! u wind tension (kg m/s)/(m**2 s) or Pascal
361	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v ! v wind tension (kg m/s)/(m**2 s) or Pascal
362
363	! Output not needed
364	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_t ! change of sensible heat flux
365	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_q ! change of water vapour flux
366	REAL, DIMENSION(klon), INTENT(OUT) :: zxsnow ! snow at surface, mean for each grid point
367	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxt ! sensible heat flux, mean for each grid point
368	REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxq ! water vapour flux, mean for each grid point
369	REAL, DIMENSION(klon, nbsrf),INTENT(OUT) :: q2m ! water vapour at 2 meter height
370	REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q ! water vapour flux(latent flux) (kg/m**2/s)
371
372
373	! Local variables with attribute SAVE
374	!****************************************************************************************
375	INTEGER, SAVE :: nhoridbg, nidbg ! variables for IOIPSL
376	!$OMP THREADPRIVATE(nhoridbg, nidbg)
377	LOGICAL, SAVE :: debugindex=.FALSE.
378	!$OMP THREADPRIVATE(debugindex)
379	LOGICAL, SAVE :: first_call=.TRUE.
380	!$OMP THREADPRIVATE(first_call)
381	CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
382	!$OMP THREADPRIVATE(cl_surf)
383
384	! Other local variables
385	!****************************************************************************************
386	INTEGER :: i, k, nsrf
387	INTEGER :: knon, j
388	INTEGER :: idayref
389	INTEGER , DIMENSION(klon) :: ni
390	REAL :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
391	REAL :: amn, amx
392	REAL :: f1 ! fraction de longeurs visibles parmi tout SW intervalle
393	REAL, DIMENSION(klon) :: r_co2_ppm ! taux CO2 atmosphere
394	REAL, DIMENSION(klon) :: yts, yrugos, ypct, yz0_new
395	REAL, DIMENSION(klon) :: yalb, yalb1, yalb2
396	REAL, DIMENSION(klon) :: yu1, yv1,ytoto
397	REAL, DIMENSION(klon) :: ysnow, yqsurf, yagesno, yqsol
398	REAL, DIMENSION(klon) :: yrain_f, ysnow_f
399	REAL, DIMENSION(klon) :: ysolsw, ysollw
400	REAL, DIMENSION(klon) :: yfder
401	REAL, DIMENSION(klon) :: yrugoro
402	REAL, DIMENSION(klon) :: yfluxlat
403	REAL, DIMENSION(klon) :: y_d_ts
404	REAL, DIMENSION(klon) :: y_flux_t1, y_flux_q1
405	REAL, DIMENSION(klon) :: y_dflux_t, y_dflux_q
406	REAL, DIMENSION(klon) :: y_flux_u1, y_flux_v1
407	REAL, DIMENSION(klon) :: yt2m, yq2m, yu10m
408	REAL, DIMENSION(klon) :: yustar
409	REAL, DIMENSION(klon) :: ywstar
410	REAL, DIMENSION(klon) :: ywindsp
411	REAL, DIMENSION(klon) :: yt10m, yq10m
412	REAL, DIMENSION(klon) :: ypblh
413	REAL, DIMENSION(klon) :: ylcl
414	REAL, DIMENSION(klon) :: ycapCL
415	REAL, DIMENSION(klon) :: yoliqCL
416	REAL, DIMENSION(klon) :: ycteiCL
417	REAL, DIMENSION(klon) :: ypblT
418	REAL, DIMENSION(klon) :: ytherm
419	REAL, DIMENSION(klon) :: ytrmb1
420	REAL, DIMENSION(klon) :: ytrmb2
421	REAL, DIMENSION(klon) :: ytrmb3
422	REAL, DIMENSION(klon) :: uzon, vmer
423	REAL, DIMENSION(klon) :: tair1, qair1, tairsol
424	REAL, DIMENSION(klon) :: psfce, patm
425	REAL, DIMENSION(klon) :: qairsol, zgeo1
426	REAL, DIMENSION(klon) :: rugo1
427	REAL, DIMENSION(klon) :: yfluxsens
428	REAL, DIMENSION(klon) :: AcoefH, AcoefQ, BcoefH, BcoefQ
429	REAL, DIMENSION(klon) :: AcoefU, AcoefV, BcoefU, BcoefV
430	REAL, DIMENSION(klon) :: ypsref
431	REAL, DIMENSION(klon) :: yevap, ytsurf_new, yalb1_new, yalb2_new
432	REAL, DIMENSION(klon) :: ztsol
433	REAL, DIMENSION(klon) :: alb_m ! mean albedo for whole SW interval
434	REAL, DIMENSION(klon,klev) :: y_d_t, y_d_q, y_d_t_diss
435	REAL, DIMENSION(klon,klev) :: y_d_u, y_d_v
436	REAL, DIMENSION(klon,klev) :: y_flux_t, y_flux_q
437	REAL, DIMENSION(klon,klev) :: y_flux_u, y_flux_v
438	REAL, DIMENSION(klon,klev) :: ycoefh, ycoefm,ycoefq
439	REAL, DIMENSION(klon) :: ycdragh, ycdragm
440	REAL, DIMENSION(klon,klev) :: yu, yv
441	REAL, DIMENSION(klon,klev) :: yt, yq
442	REAL, DIMENSION(klon,klev) :: ypplay, ydelp
443	REAL, DIMENSION(klon,klev) :: delp
444	REAL, DIMENSION(klon,klev+1) :: ypaprs
445	REAL, DIMENSION(klon,klev+1) :: ytke
446	REAL, DIMENSION(klon,nsoilmx) :: ytsoil
447	CHARACTER(len=80) :: abort_message
448	CHARACTER(len=20) :: modname = 'pbl_surface'
449	LOGICAL, PARAMETER :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
450	LOGICAL, PARAMETER :: check=.FALSE.
451	REAL, DIMENSION(klon) :: Kech_h ! Coefficient d'echange pour l'energie
452
453	! For debugging with IOIPSL
454	INTEGER, DIMENSION(iim*(jjm+1)) :: ndexbg
455	REAL :: zjulian
456	REAL, DIMENSION(klon) :: tabindx
457	REAL, DIMENSION(iim,jjm+1) :: zx_lon, zx_lat
458	REAL, DIMENSION(iim,jjm+1) :: debugtab
459
460
461	REAL, DIMENSION(klon,nbsrf) :: pblh ! height of the planetary boundary layer
462	REAL, DIMENSION(klon,nbsrf) :: plcl ! condensation level
463	REAL, DIMENSION(klon,nbsrf) :: capCL
464	REAL, DIMENSION(klon,nbsrf) :: oliqCL
465	REAL, DIMENSION(klon,nbsrf) :: cteiCL
466	REAL, DIMENSION(klon,nbsrf) :: pblT
467	REAL, DIMENSION(klon,nbsrf) :: therm
468	REAL, DIMENSION(klon,nbsrf) :: trmb1 ! deep cape
469	REAL, DIMENSION(klon,nbsrf) :: trmb2 ! inhibition
470	REAL, DIMENSION(klon,nbsrf) :: trmb3 ! point Omega
471	REAL, DIMENSION(klon,nbsrf) :: zx_rh2m, zx_qsat2m
472	REAL, DIMENSION(klon,nbsrf) :: zx_t1
473	REAL, DIMENSION(klon, nbsrf) :: alb ! mean albedo for whole SW interval
474	REAL, DIMENSION(klon) :: ylwdown ! jg : temporary (ysollwdown)
475
476	REAL :: zx_qs1, zcor1, zdelta1
477
478	!****************************************************************************************
479	! Declarations specifiques pour le 1D. A reprendre
480	!****************************************************************************************
481	REAL :: fsens,flat
482	LOGICAL :: ok_flux_surf ! initialized during first_call below
483	COMMON /flux_arp/fsens,flat,ok_flux_surf
484	! End of declarations
485	!****************************************************************************************
486
487
488	!****************************************************************************************
489	! 1) Initialisation and validation tests
490	! Only done first time entering this subroutine
491	!
492	!****************************************************************************************
493
494	IF (first_call) THEN
495	first_call=.FALSE.
496
497	! Initialize ok_flux_surf (for 1D model)
498	if (klon>1) ok_flux_surf=.FALSE.
499
500	! Initilize debug IO
501	IF (debugindex .AND. mpi_size==1) THEN
502	! initialize IOIPSL output
503	idayref = day_ini
504	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
505	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
506	DO i = 1, iim
507	zx_lon(i,1) = rlon(i+1)
508	zx_lon(i,jjm+1) = rlon(i+1)
509	ENDDO
510	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
511	CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
512	1,iim,1,jjm+1, &
513	itau_phy,zjulian,dtime,nhoridbg,nidbg)
514	! no vertical axis
515	cl_surf(1)='ter'
516	cl_surf(2)='lic'
517	cl_surf(3)='oce'
518	cl_surf(4)='sic'
519	DO nsrf=1,nbsrf
520	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
521	jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
522	END DO
523
524	CALL histend(nidbg)
525	CALL histsync(nidbg)
526
527	END IF
528
529	ENDIF
530
531	!****************************************************************************************
532	! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
533	! instead of ORCHIDEE)
534	IF (qsol0>0.) THEN
535	PRINT*,'WARNING : On impose qsol=',qsol0
536	qsol(:)=qsol0
537	ENDIF
538	!****************************************************************************************
539
540	!****************************************************************************************
541	! 2) Initialization to zero
542	! Done for all local variables that will be compressed later
543	! and argument with INTENT(OUT)
544	!****************************************************************************************
545	cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0
546	ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0
547	zv1 = 0.0 ; yqsurf = 0.0 ; yalb1 = 0.0 ; yalb2 = 0.0
548	yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0
549	ysollw = 0.0 ; yrugos = 0.0 ; yu1 = 0.0
550	yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0
551	ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0
552	yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
553	yrugoro = 0.0 ; ywindsp = 0.0
554	d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0
555	flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0
556	d_t_diss= 0.0 ;d_u = 0.0 ; d_v = 0.0 ; yqsol = 0.0
557	ytherm = 0.0 ; ytke=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	END DO
715
716	DO k = 1, klev
717	DO j = 1, knon
718	i = ni(j)
719	ypaprs(j,k) = paprs(i,k)
720	ypplay(j,k) = pplay(i,k)
721	ydelp(j,k) = delp(i,k)
722	ytke(j,k) = tke(i,k,nsrf)
723	yu(j,k) = u(i,k)
724	yv(j,k) = v(i,k)
725	yt(j,k) = t(i,k)
726	yq(j,k) = q(i,k)
727	ENDDO
728	ENDDO
729
730	DO k = 1, nsoilmx
731	DO j = 1, knon
732	i = ni(j)
733	ytsoil(j,k) = ftsoil(i,k,nsrf)
734	END DO
735	END DO
736
737	! qsol(water height in soil) only for bucket continental model
738	IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
739	DO j = 1, knon
740	i = ni(j)
741	yqsol(j) = qsol(i)
742	END DO
743	ENDIF
744
745	!****************************************************************************************
746	! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
747	!
748	!****************************************************************************************
749
750	CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
751	yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
752	yts, yqsurf, yrugos, &
753	ycdragm, ycdragh )
754
755	!****************************************************************************************
756	! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefm et ycoefm.
757	!
758	!****************************************************************************************
759
760	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
761	ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
762	ycoefm, ycoefh, ytke)
763
764	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
765	! In this case, coef_diff_turb is called for the Cd only
766	DO k = 2, klev
767	DO j = 1, knon
768	i = ni(j)
769	ycoefh(j,k) = zcoefh(i,k,nsrf)
770	ycoefm(j,k) = zcoefm(i,k,nsrf)
771	ENDDO
772	ENDDO
773	ENDIF
774
775	!****************************************************************************************
776	!
777	! 8) "La descente" - "The downhill"
778	!
779	! climb_hq_down and climb_wind_down calculate the coefficients
780	! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
781	! Only the coefficients at surface for H and Q are returned.
782	!
783	!****************************************************************************************
784
785	! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
786	CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
787	ydelp, yt, yq, dtime, &
788	AcoefH, AcoefQ, BcoefH, BcoefQ)
789
790	! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
791	CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
792	AcoefU, AcoefV, BcoefU, BcoefV)
793
794
795	!****************************************************************************************
796	! 9) Small calculations
797	!
798	!****************************************************************************************
799
800	! - Reference pressure is given the values at surface level
801	ypsref(:) = ypaprs(:,1)
802
803	! - CO2 field on 2D grid to be sent to ORCHIDEE
804	! Transform to compressed field
805	IF (carbon_cycle_cpl) THEN
806	DO i=1,knon
807	r_co2_ppm(i) = co2_send(ni(i))
808	END DO
809	ELSE
810	r_co2_ppm(:) = co2_ppm ! Constant field
811	END IF
812
813	!****************************************************************************************
814	!
815	! Calulate t2m and q2m for the case of calculation at land grid points
816	! t2m and q2m are needed as input to ORCHIDEE
817	!
818	!****************************************************************************************
819	IF (nsrf == is_ter) THEN
820
821	DO i = 1, knon
822	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
823	* (ypaprs(i,1)-ypplay(i,1))
824	END DO
825
826	! Calculate the temperature et relative humidity at 2m and the wind at 10m
827	CALL stdlevvar(klon, knon, is_ter, zxli, &
828	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
829	yts, yqsurf, yrugos, ypaprs(:,1), ypplay(:,1), &
830	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
831
832	END IF
833
834	!****************************************************************************************
835	!
836	! 10) Switch selon current surface
837	! It is necessary to start with the continental surfaces because the ocean
838	! needs their run-off.
839	!
840	!****************************************************************************************
841	SELECT CASE(nsrf)
842
843	CASE(is_ter)
844	! ylwdown : to be removed, calculation is now done at land surface in surf_land
845	ylwdown(:)=0.0
846	DO i=1,knon
847	ylwdown(i)=lwdown_m(ni(i))
848	END DO
849	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
850	rlon, rlat, &
851	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
852	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
853	AcoefH, AcoefQ, BcoefH, BcoefQ, &
854	AcoefU, AcoefV, BcoefU, BcoefV, &
855	ypsref, yu1, yv1, yrugoro, pctsrf, &
856	ylwdown, yq2m, yt2m, &
857	ysnow, yqsol, yagesno, ytsoil, &
858	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
859	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
860	y_flux_u1, y_flux_v1 )
861
862
863	CASE(is_lic)
864	CALL surf_landice(itap, dtime, knon, ni, &
865	ysolsw, ysollw, yts, ypplay(:,1), &
866	ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
867	AcoefH, AcoefQ, BcoefH, BcoefQ, &
868	AcoefU, AcoefV, BcoefU, BcoefV, &
869	ypsref, yu1, yv1, yrugoro, pctsrf, &
870	ysnow, yqsurf, yqsol, yagesno, &
871	ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
872	ytsurf_new, y_dflux_t, y_dflux_q, &
873	y_flux_u1, y_flux_v1)
874
875	CASE(is_oce)
876	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb1, &
877	yrugos, ywindsp, rmu0, yfder, yts, &
878	itap, dtime, jour, knon, ni, &
879	ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
880	AcoefH, AcoefQ, BcoefH, BcoefQ, &
881	AcoefU, AcoefV, BcoefU, BcoefV, &
882	ypsref, yu1, yv1, yrugoro, pctsrf, &
883	ysnow, yqsurf, yagesno, &
884	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
885	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
886	y_flux_u1, y_flux_v1)
887
888	CASE(is_sic)
889	CALL surf_seaice( &
890	rlon, rlat, ysolsw, ysollw, yalb1, yfder, &
891	itap, dtime, jour, knon, ni, &
892	lafin, &
893	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
894	AcoefH, AcoefQ, BcoefH, BcoefQ, &
895	AcoefU, AcoefV, BcoefU, BcoefV, &
896	ypsref, yu1, yv1, yrugoro, pctsrf, &
897	ysnow, yqsurf, yqsol, yagesno, ytsoil, &
898	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
899	ytsurf_new, y_dflux_t, y_dflux_q, &
900	y_flux_u1, y_flux_v1)
901
902
903	CASE DEFAULT
904	WRITE(lunout,*) 'Surface index = ', nsrf
905	abort_message = 'Surface index not valid'
906	CALL abort_gcm(modname,abort_message,1)
907	END SELECT
908
909
910	!****************************************************************************************
911	! 11) - Calcul the increment of surface temperature
912	!
913	!****************************************************************************************
914	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
915
916	!****************************************************************************************
917	!
918	! 12) "La remontee" - "The uphill"
919	!
920	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
921	! for X=H, Q, U and V, for all vertical levels.
922	!
923	!****************************************************************************************
924	! H and Q
925	IF (ok_flux_surf) THEN
926	PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
927	y_flux_t1(:) = fsens
928	y_flux_q1(:) = flat/RLVTT
929	yfluxlat(:) = flat
930
931	Kech_h(:) = ycdragh(:) * (1.0+SQRT(yu(:,1)2+yv(:,1)2)) * &
932	ypplay(:,1)/(RD*yt(:,1))
933	ytoto(:)=(1./RCPD)(AcoefH(:)+BcoefH(:)y_flux_t1(:)*dtime)
934	ytsurf_new(:)=ytoto(:)-y_flux_t1(:)/(Kech_h(:)*RCPD)
935	y_d_ts(:) = ytsurf_new(:) - yts(:)
936
937	ELSE
938	y_flux_t1(:) = yfluxsens(:)
939	y_flux_q1(:) = -yevap(:)
940	ENDIF
941
942	CALL climb_hq_up(knon, dtime, yt, yq, &
943	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
944	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
945
946
947	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
948	y_flux_u, y_flux_v, y_d_u, y_d_v)
949
950
951	y_d_t_diss(:,:)=0.
952	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
953	CALL yamada_c(knon,dtime,ypaprs,ypplay &
954	& ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
955	& ,iflag_pbl,nsrf)
956	ENDIF
957	! print*,'yamada_c OK'
958
959	DO j = 1, knon
960	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
961	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
962	ENDDO
963
964	!****************************************************************************************
965	! 13) Transform variables for output format :
966	! - Decompress
967	! - Multiply with pourcentage of current surface
968	! - Cumulate in global variable
969	!
970	!****************************************************************************************
971
972	DO k = 1, klev
973	DO j = 1, knon
974	i = ni(j)
975	y_d_t_diss(j,k) = y_d_t_diss(j,k) * ypct(j)
976	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
977	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
978	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
979	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
980
981	flux_t(i,k,nsrf) = y_flux_t(j,k)
982	flux_q(i,k,nsrf) = y_flux_q(j,k)
983	flux_u(i,k,nsrf) = y_flux_u(j,k)
984	flux_v(i,k,nsrf) = y_flux_v(j,k)
985
986
987	ENDDO
988	ENDDO
989
990	! print*,'Dans pbl OK1'
991
992	evap(:,nsrf) = - flux_q(:,1,nsrf)
993
994	alb1(:, nsrf) = 0.
995	alb2(:, nsrf) = 0.
996	snow(:, nsrf) = 0.
997	qsurf(:, nsrf) = 0.
998	rugos(:, nsrf) = 0.
999	fluxlat(:,nsrf) = 0.
1000	DO j = 1, knon
1001	i = ni(j)
1002	d_ts(i,nsrf) = y_d_ts(j)
1003	alb1(i,nsrf) = yalb1_new(j)
1004	alb2(i,nsrf) = yalb2_new(j)
1005	snow(i,nsrf) = ysnow(j)
1006	qsurf(i,nsrf) = yqsurf(j)
1007	rugos(i,nsrf) = yz0_new(j)
1008	fluxlat(i,nsrf) = yfluxlat(j)
1009	agesno(i,nsrf) = yagesno(j)
1010	cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
1011	cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
1012	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
1013	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
1014	END DO
1015
1016	! print*,'Dans pbl OK2'
1017
1018	DO k = 2, klev
1019	DO j = 1, knon
1020	i = ni(j)
1021	tke(i,k,nsrf) = ytke(j,k)
1022	zcoefh(i,k,nsrf) = ycoefh(j,k)
1023	zcoefm(i,k,nsrf) = ycoefm(j,k)
1024	tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
1025	zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
1026	zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
1027	END DO
1028	END DO
1029
1030	! print*,'Dans pbl OK3'
1031
1032	IF ( nsrf .EQ. is_ter ) THEN
1033	DO j = 1, knon
1034	i = ni(j)
1035	qsol(i) = yqsol(j)
1036	END DO
1037	END IF
1038
1039	ftsoil(:,:,nsrf) = 0.
1040	DO k = 1, nsoilmx
1041	DO j = 1, knon
1042	i = ni(j)
1043	ftsoil(i, k, nsrf) = ytsoil(j,k)
1044	END DO
1045	END DO
1046
1047
1048	DO k = 1, klev
1049	DO j = 1, knon
1050	i = ni(j)
1051	d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
1052	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
1053	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
1054	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
1055	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
1056	END DO
1057	END DO
1058
1059	! print*,'Dans pbl OK4'
1060
1061	!****************************************************************************************
1062	! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m
1063	! Call HBTM
1064	!
1065	!****************************************************************************************
1066	t2m(:,nsrf) = 0.
1067	q2m(:,nsrf) = 0.
1068	ustar(:,nsrf) = 0.
1069	wstar(:,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,ywstar, &
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	wstar(i,nsrf) = ywstar(j)
1160	plcl(i,nsrf) = ylcl(j)
1161	capCL(i,nsrf) = ycapCL(j)
1162	oliqCL(i,nsrf) = yoliqCL(j)
1163	cteiCL(i,nsrf) = ycteiCL(j)
1164	pblT(i,nsrf) = ypblT(j)
1165	therm(i,nsrf) = ytherm(j)
1166	trmb1(i,nsrf) = ytrmb1(j)
1167	trmb2(i,nsrf) = ytrmb2(j)
1168	trmb3(i,nsrf) = ytrmb3(j)
1169	END DO
1170
1171	! print*,'OK pbl 6'
1172	#else
1173	! T2m not defined
1174	! No calculation
1175	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
1176	#endif
1177
1178	!****************************************************************************************
1179	! 15) End of loop over different surfaces
1180	!
1181	!****************************************************************************************
1182	END DO loop_nbsrf
1183
1184	!****************************************************************************************
1185	! 16) Calculate the mean value over all sub-surfaces for som variables
1186	!
1187	!****************************************************************************************
1188
1189	! print*,'OK pbl 7'
1190	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
1191	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
1192	DO nsrf = 1, nbsrf
1193	DO k = 1, klev
1194	DO i = 1, klon
1195	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
1196	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
1197	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
1198	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
1199	END DO
1200	END DO
1201	END DO
1202
1203	! print*,'OK pbl 8'
1204	DO i = 1, klon
1205	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1206	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1207	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
1208	ENDDO
1209
1210	!
1211	! Incrementer la temperature du sol
1212	!
1213	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
1214	zt2m(:) = 0.0 ; zq2m(:) = 0.0
1215	zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
1216	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
1217	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
1218	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
1219	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
1220	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
1221	wstar(:,is_ave)=0.
1222
1223	! print*,'OK pbl 9'
1224
1225	DO nsrf = 1, nbsrf
1226	DO i = 1, klon
1227	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
1228
1229	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
1230	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1231	wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
1232	pctsrf(i,nsrf)
1233
1234	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
1235	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
1236
1237	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
1238	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
1239	zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
1240	wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
1241	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
1242	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
1243
1244	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
1245	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
1246	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
1247	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
1248	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
1249	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
1250	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
1251	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
1252	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
1253	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
1254	END DO
1255	END DO
1256	! print*,'OK pbl 10'
1257
1258	IF (check) THEN
1259	amn=MIN(ts(1,is_ter),1000.)
1260	amx=MAX(ts(1,is_ter),-1000.)
1261	DO i=2, klon
1262	amn=MIN(ts(i,is_ter),amn)
1263	amx=MAX(ts(i,is_ter),amx)
1264	ENDDO
1265	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
1266	ENDIF
1267
1268	!jg ?
1269	!!$!
1270	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
1271	!!$! sub-surfaces is distributed.
1272	!!$!
1273	!!$ DO nsrf = 1, nbsrf
1274	!!$ DO i = 1, klon
1275	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
1276	!!$ ts(i,nsrf) = zxtsol(i)
1277	!!$ t2m(i,nsrf) = zt2m(i)
1278	!!$ q2m(i,nsrf) = zq2m(i)
1279	!!$ u10m(i,nsrf) = zu10m(i)
1280	!!$ v10m(i,nsrf) = zv10m(i)
1281	!!$
1282	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
1283	!!$ pblh(i,nsrf) = s_pblh(i)
1284	!!$ plcl(i,nsrf) = s_plcl(i)
1285	!!$ capCL(i,nsrf) = s_capCL(i)
1286	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
1287	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
1288	!!$ pblT(i,nsrf) = s_pblT(i)
1289	!!$ therm(i,nsrf) = s_therm(i)
1290	!!$ trmb1(i,nsrf) = s_trmb1(i)
1291	!!$ trmb2(i,nsrf) = s_trmb2(i)
1292	!!$ trmb3(i,nsrf) = s_trmb3(i)
1293	!!$ ENDIF
1294	!!$ ENDDO
1295	!!$ ENDDO
1296
1297
1298	DO i = 1, klon
1299	fder(i) = - 4.0RSIGMAzxtsol(i)**3
1300	ENDDO
1301
1302	zxqsurf(:) = 0.0
1303	zxsnow(:) = 0.0
1304	DO nsrf = 1, nbsrf
1305	DO i = 1, klon
1306	zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
1307	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
1308	END DO
1309	END DO
1310
1311	! Premier niveau de vent sortie dans physiq.F
1312	zu1(:) = u(:,1)
1313	zv1(:) = v(:,1)
1314
1315	! Some of the module declared variables are returned for printing in physiq.F
1316	qsol_d(:) = qsol(:)
1317	evap_d(:,:) = evap(:,:)
1318	rugos_d(:,:) = rugos(:,:)
1319	agesno_d(:,:) = agesno(:,:)
1320
1321
1322	END SUBROUTINE pbl_surface
1323	!
1324	!****************************************************************************************
1325	!
1326	SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
1327	evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
1328
1329	USE indice_sol_mod
1330
1331	INCLUDE "dimsoil.h"
1332
1333	! Ouput variables
1334	!****************************************************************************************
1335	REAL, DIMENSION(klon), INTENT(OUT) :: qsol_rst
1336	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
1337	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
1338	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
1339	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_rst
1340	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_rst
1341	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_rst
1342	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
1343
1344
1345	!****************************************************************************************
1346	! Return module variables for writing to restart file
1347	!
1348	!****************************************************************************************
1349	qsol_rst(:) = qsol(:)
1350	fder_rst(:) = fder(:)
1351	snow_rst(:,:) = snow(:,:)
1352	qsurf_rst(:,:) = qsurf(:,:)
1353	evap_rst(:,:) = evap(:,:)
1354	rugos_rst(:,:) = rugos(:,:)
1355	agesno_rst(:,:) = agesno(:,:)
1356	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
1357
1358	!****************************************************************************************
1359	! Deallocate module variables
1360	!
1361	!****************************************************************************************
1362	! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
1363	IF (ALLOCATED(qsol)) DEALLOCATE(qsol)
1364	IF (ALLOCATED(fder)) DEALLOCATE(fder)
1365	IF (ALLOCATED(snow)) DEALLOCATE(snow)
1366	IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
1367	IF (ALLOCATED(evap)) DEALLOCATE(evap)
1368	IF (ALLOCATED(rugos)) DEALLOCATE(rugos)
1369	IF (ALLOCATED(agesno)) DEALLOCATE(agesno)
1370	IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
1371
1372	END SUBROUTINE pbl_surface_final
1373	!
1374	!****************************************************************************************
1375	!
1376	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, tsurf, alb1, alb2, ustar, u10m, v10m, tke)
1377
1378	! Give default values where new fraction has appread
1379
1380	USE indice_sol_mod
1381
1382	INCLUDE "dimsoil.h"
1383	INCLUDE "clesphys.h"
1384	INCLUDE "compbl.h"
1385
1386	! Input variables
1387	!****************************************************************************************
1388	INTEGER, INTENT(IN) :: itime
1389	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
1390
1391	! InOutput variables
1392	!****************************************************************************************
1393	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
1394	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: alb1, alb2
1395	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m
1396	REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke
1397
1398	! Local variables
1399	!****************************************************************************************
1400	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
1401	CHARACTER(len=80) :: abort_message
1402	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
1403	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
1404	!
1405	! All at once !!
1406	!****************************************************************************************
1407
1408	DO nsrf = 1, nbsrf
1409	! First decide complement sub-surfaces
1410	SELECT CASE (nsrf)
1411	CASE(is_oce)
1412	nsrf_comp1=is_sic
1413	nsrf_comp2=is_ter
1414	nsrf_comp3=is_lic
1415	CASE(is_sic)
1416	nsrf_comp1=is_oce
1417	nsrf_comp2=is_ter
1418	nsrf_comp3=is_lic
1419	CASE(is_ter)
1420	nsrf_comp1=is_lic
1421	nsrf_comp2=is_oce
1422	nsrf_comp3=is_sic
1423	CASE(is_lic)
1424	nsrf_comp1=is_ter
1425	nsrf_comp2=is_oce
1426	nsrf_comp3=is_sic
1427	END SELECT
1428
1429	! Initialize all new fractions
1430	DO i=1, klon
1431	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
1432
1433	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
1434	! Use the complement sub-surface, keeping the continents unchanged
1435	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
1436	evap(i,nsrf) = evap(i,nsrf_comp1)
1437	rugos(i,nsrf) = rugos(i,nsrf_comp1)
1438	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
1439	alb1(i,nsrf) = alb1(i,nsrf_comp1)
1440	alb2(i,nsrf) = alb2(i,nsrf_comp1)
1441	ustar(i,nsrf) = ustar(i,nsrf_comp1)
1442	u10m(i,nsrf) = u10m(i,nsrf_comp1)
1443	v10m(i,nsrf) = v10m(i,nsrf_comp1)
1444	if (iflag_pbl > 1) then
1445	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
1446	endif
1447	mfois(nsrf) = mfois(nsrf) + 1
1448	ELSE
1449	! The continents have changed. The new fraction receives the mean sum of the existent fractions
1450	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1451	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1452	rugos(i,nsrf) = rugos(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + rugos(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1453	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1454	alb1(i,nsrf) = alb1(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + alb1(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1455	alb2(i,nsrf) = alb2(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + alb2(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1456	ustar(i,nsrf) = ustar(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1457	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1458	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1459	if (iflag_pbl > 1) then
1460	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1461	endif
1462
1463	! Security abort. This option has never been tested. To test, comment the following line.
1464	! abort_message='The fraction of the continents have changed!'
1465	! CALL abort_gcm(modname,abort_message,1)
1466	nfois(nsrf) = nfois(nsrf) + 1
1467	END IF
1468	snow(i,nsrf) = 0.
1469	agesno(i,nsrf) = 0.
1470	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
1471	ELSE
1472	pfois(nsrf) = pfois(nsrf)+ 1
1473	END IF
1474	END DO
1475
1476	END DO
1477
1478	END SUBROUTINE pbl_surface_newfrac
1479
1480	!
1481	!****************************************************************************************
1482	!
1483
1484	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format