Context Navigation

pbl_surface_mod.F90 @ 994

Last change on this file since 994 was 987, checked in by Laurent Fairhead, 17 years ago
Du nettoyage sur le parallelisme, inclusion de nouvelles interfaces pour OPA9 et ORCHIDEE YM LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 54.5 KB

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

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

Download in other formats:

Original Format