Context Navigation

pbl_surface_mod.F90 @ 1780

Last change on this file since 1780 was 1780, checked in by Laurent Fairhead, 11 years ago

Inclusion des cas AMMA et Sandu dans la configuration 1D
MPL

AMMA and Sandu cases included in 1D configuration
MPL

Property svn:executable set to *

File size: 61.6 KB

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

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

Download in other formats:

Original Format