Context Navigation

pbl_surface_mod.F90 @ 1139

Last change on this file since 1139 was 1069, checked in by Laurent Fairhead, 17 years ago

Modification pour utiliser albedo des deux intervalles visible et proche
infra rouge pour la repartition sur sous-surface des flux solaire.

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 56.8 KB

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

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

Download in other formats:

Original Format