Context Navigation

pbl_surface_mod.F90 @ 1142

Last change on this file since 1142 was 1142, checked in by yann meurdesoif, 16 years ago

Portage vers Vargas

Y.M + E.M

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,SAVE :: fsens,flat
463	LOGICAL,SAVE :: ok_flux_surf=.FALSE.
464	!$OMP THREADPRIVATE(fsens,flat,ok_flux_surf)
465
466	!****************************************************************************************
467	! End of declarations
468	!****************************************************************************************
469
470
471	!****************************************************************************************
472	! 1) Initialisation and validation tests
473	! Only done first time entering this subroutine
474	!
475	!****************************************************************************************
476
477	IF (first_call) THEN
478	first_call=.FALSE.
479
480	! Initilize debug IO
481	IF (debugindex .AND. mpi_size==1) THEN
482	! initialize IOIPSL output
483	idayref = day_ini
484	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
485	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
486	DO i = 1, iim
487	zx_lon(i,1) = rlon(i+1)
488	zx_lon(i,jjm+1) = rlon(i+1)
489	ENDDO
490	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
491	CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
492	1,iim,1,jjm+1, &
493	itau_phy,zjulian,dtime,nhoridbg,nidbg)
494	! no vertical axis
495	cl_surf(1)='ter'
496	cl_surf(2)='lic'
497	cl_surf(3)='oce'
498	cl_surf(4)='sic'
499	DO nsrf=1,nbsrf
500	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
501	jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
502	END DO
503
504	CALL histend(nidbg)
505	CALL histsync(nidbg)
506
507	END IF
508
509	ENDIF
510
511	!****************************************************************************************
512	! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
513	! instead of ORCHIDEE)
514	IF (qsol0>0.) THEN
515	PRINT*,'WARNING : On impose qsol=',qsol0
516	qsol(:)=qsol0
517	ENDIF
518	!****************************************************************************************
519
520	!****************************************************************************************
521	! 2) Initialization to zero
522	! Done for all local variables that will be compressed later
523	! and argument with INTENT(OUT)
524	!****************************************************************************************
525	cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0
526	ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0
527	zv1 = 0.0 ; yqsurf = 0.0 ; yalb1 = 0.0 ; yalb2 = 0.0
528	yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0
529	ysollw = 0.0 ; yrugos = 0.0 ; yu1 = 0.0
530	yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0
531	ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0
532	yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
533	yrugoro = 0.0 ; ywindsp = 0.0
534	d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0
535	flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0
536	d_u = 0.0 ; d_v = 0.0 ; yqsol = 0.0
537	ytherm = 0.0 ; ytke=0.
538
539	zcoefh(:,:) = 0.0
540	zcoefh(:,1) = 999999. ! zcoefh(:,k=1) should never be used
541	ytsoil = 999999.
542
543	rh2m(:) = 0.
544	qsat2m(:) = 0.
545	!****************************************************************************************
546	! 3) - Calculate pressure thickness of each layer
547	! - Calculate the wind at first layer
548	! - Mean calculations of albedo
549	! - Calculate net radiance at sub-surface
550	!****************************************************************************************
551	DO k = 1, klev
552	DO i = 1, klon
553	delp(i,k) = paprs(i,k)-paprs(i,k+1)
554	ENDDO
555	ENDDO
556
557	!****************************************************************************************
558	! Test for rugos........ from physiq.. A la fin plutot???
559	!
560	!****************************************************************************************
561
562	zxrugs(:) = 0.0
563	DO nsrf = 1, nbsrf
564	DO i = 1, klon
565	rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015)
566	zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf)
567	ENDDO
568	ENDDO
569
570	! Mean calculations of albedo
571	!
572	! Albedo at sub-surface
573	! * alb1 : albedo in visible SW interval
574	! * alb2 : albedo in near infrared SW interval
575	! * alb : mean albedo for whole SW interval
576	!
577	! Mean albedo for grid point
578	! * alb1_m : albedo in visible SW interval
579	! * alb2_m : albedo in near infrared SW interval
580	! * alb_m : mean albedo at whole SW interval
581
582	alb1_m(:) = 0.0
583	alb2_m(:) = 0.0
584	DO nsrf = 1, nbsrf
585	DO i = 1, klon
586	alb1_m(i) = alb1_m(i) + alb1(i,nsrf) * pctsrf(i,nsrf)
587	alb2_m(i) = alb2_m(i) + alb2(i,nsrf) * pctsrf(i,nsrf)
588	ENDDO
589	ENDDO
590
591	! We here suppose the fraction f1 of incoming radiance of visible radiance
592	! as a fraction of all shortwave radiance
593	f1 = 0.5
594	! f1 = 1 ! put f1=1 to recreate old calculations
595
596	DO nsrf = 1, nbsrf
597	DO i = 1, klon
598	alb(i,nsrf) = f1alb1(i,nsrf) + (1-f1)alb2(i,nsrf)
599	ENDDO
600	ENDDO
601
602	DO i = 1, klon
603	alb_m(i) = f1alb1_m(i) + (1-f1)alb2_m(i)
604	END DO
605
606	! Calculation of mean temperature at surface grid points
607	ztsol(:) = 0.0
608	DO nsrf = 1, nbsrf
609	DO i = 1, klon
610	ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
611	ENDDO
612	ENDDO
613
614	! Linear distrubution on sub-surface of long- and shortwave net radiance
615	DO nsrf = 1, nbsrf
616	DO i = 1, klon
617	sollw(i,nsrf) = sollw_m(i) + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ts(i,nsrf))
618	solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
619	ENDDO
620	ENDDO
621
622
623	! Downwelling longwave radiation at mean surface
624	lwdown_m(:) = 0.0
625	DO i = 1, klon
626	lwdown_m(i) = sollw_m(i) + RSIGMAztsol(i)*4
627	ENDDO
628
629	!****************************************************************************************
630	! 4) Loop over different surfaces
631	!
632	! Only points containing a fraction of the sub surface will be threated.
633	!
634	!****************************************************************************************
635
636	loop_nbsrf: DO nsrf = 1, nbsrf
637
638	! Search for index(ni) and size(knon) of domaine to treat
639	ni(:) = 0
640	knon = 0
641	DO i = 1, klon
642	IF (pctsrf(i,nsrf) > 0.) THEN
643	knon = knon + 1
644	ni(knon) = i
645	ENDIF
646	ENDDO
647
648	! write index, with IOIPSL
649	IF (debugindex .AND. mpi_size==1) THEN
650	tabindx(:)=0.
651	DO i=1,knon
652	tabindx(i)=FLOAT(i)
653	END DO
654	debugtab(:,:) = 0.
655	ndexbg(:) = 0
656	CALL gath2cpl(tabindx,debugtab,knon,ni)
657	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg)
658	ENDIF
659
660	!****************************************************************************************
661	! 5) Compress variables
662	!
663	!****************************************************************************************
664
665	DO j = 1, knon
666	i = ni(j)
667	ypct(j) = pctsrf(i,nsrf)
668	yts(j) = ts(i,nsrf)
669	ysnow(j) = snow(i,nsrf)
670	yqsurf(j) = qsurf(i,nsrf)
671	yalb(j) = alb(i,nsrf)
672	yalb1(j) = alb1(i,nsrf)
673	yalb2(j) = alb2(i,nsrf)
674	yrain_f(j) = rain_f(i)
675	ysnow_f(j) = snow_f(i)
676	yagesno(j) = agesno(i,nsrf)
677	yfder(j) = fder(i)
678	ysolsw(j) = solsw(i,nsrf)
679	ysollw(j) = sollw(i,nsrf)
680	yrugos(j) = rugos(i,nsrf)
681	yrugoro(j) = rugoro(i)
682	yu1(j) = u(i,1)
683	yv1(j) = v(i,1)
684	ypaprs(j,klev+1) = paprs(i,klev+1)
685	ywindsp(j) = SQRT(u10m(i,nsrf)2 + v10m(i,nsrf)2 )
686	END DO
687
688	DO k = 1, klev
689	DO j = 1, knon
690	i = ni(j)
691	ypaprs(j,k) = paprs(i,k)
692	ypplay(j,k) = pplay(i,k)
693	ydelp(j,k) = delp(i,k)
694	ytke(j,k) = tke(i,k,nsrf)
695	yu(j,k) = u(i,k)
696	yv(j,k) = v(i,k)
697	yt(j,k) = t(i,k)
698	yq(j,k) = q(i,k)
699	ENDDO
700	ENDDO
701
702	DO k = 1, nsoilmx
703	DO j = 1, knon
704	i = ni(j)
705	ytsoil(j,k) = ftsoil(i,k,nsrf)
706	END DO
707	END DO
708
709	! qsol(water height in soil) only for bucket continental model
710	IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
711	DO j = 1, knon
712	i = ni(j)
713	yqsol(j) = qsol(i)
714	END DO
715	ENDIF
716
717	!****************************************************************************************
718	! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
719	!
720	!****************************************************************************************
721
722	CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
723	yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
724	yts, yqsurf, yrugos, &
725	ycdragm, ycdragh )
726
727	!****************************************************************************************
728	! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefm et ycoefm.
729	!
730	!****************************************************************************************
731
732	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
733	ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
734	ycoefm, ycoefh, ytke)
735
736	!****************************************************************************************
737	!
738	! 8) "La descente" - "The downhill"
739	!
740	! climb_hq_down and climb_wind_down calculate the coefficients
741	! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
742	! Only the coefficients at surface for H and Q are returned.
743	!
744	!****************************************************************************************
745
746	! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
747	CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
748	ydelp, yt, yq, dtime, &
749	AcoefH, AcoefQ, BcoefH, BcoefQ)
750
751	! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
752	CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
753	AcoefU, AcoefV, BcoefU, BcoefV)
754
755
756	!****************************************************************************************
757	! 9) Small calculations
758	!
759	!****************************************************************************************
760
761	! - Reference pressure is given the values at surface level
762	ypsref(:) = ypaprs(:,1)
763
764	! - Constant CO2 is copied to global grid
765	r_co2_ppm(:) = co2_ppm
766
767	!****************************************************************************************
768	!
769	! 10) Switch selon current surface
770	! It is necessary to start with the continental surfaces because the ocean
771	! needs their run-off.
772	!
773	!****************************************************************************************
774	SELECT CASE(nsrf)
775
776	CASE(is_ter)
777	! ylwdown : to be removed, calculation is now done at land surface in surf_land
778	ylwdown(:)=0.0
779	DO i=1,knon
780	ylwdown(i)=lwdown_m(ni(i))
781	END DO
782	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
783	rlon, rlat, &
784	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
785	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
786	AcoefH, AcoefQ, BcoefH, BcoefQ, &
787	AcoefU, AcoefV, BcoefU, BcoefV, &
788	ypsref, yu1, yv1, yrugoro, pctsrf, &
789	ysnow, yqsol, yagesno, ytsoil, &
790	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
791	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
792	y_flux_u1, y_flux_v1, &
793	ylwdown)
794
795	CASE(is_lic)
796	CALL surf_landice(itap, dtime, knon, ni, &
797	ysolsw, ysollw, yts, ypplay(:,1), &
798	ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
799	AcoefH, AcoefQ, BcoefH, BcoefQ, &
800	AcoefU, AcoefV, BcoefU, BcoefV, &
801	ypsref, yu1, yv1, yrugoro, pctsrf, &
802	ysnow, yqsurf, yqsol, yagesno, &
803	ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
804	ytsurf_new, y_dflux_t, y_dflux_q, &
805	y_flux_u1, y_flux_v1)
806
807	CASE(is_oce)
808	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb1, &
809	yrugos, ywindsp, rmu0, yfder, yts, &
810	itap, dtime, jour, knon, ni, &
811	ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
812	AcoefH, AcoefQ, BcoefH, BcoefQ, &
813	AcoefU, AcoefV, BcoefU, BcoefV, &
814	ypsref, yu1, yv1, yrugoro, pctsrf, &
815	ysnow, yqsurf, yagesno, &
816	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
817	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
818	y_flux_u1, y_flux_v1)
819
820	CASE(is_sic)
821	CALL surf_seaice( &
822	rlon, rlat, ysolsw, ysollw, yalb1, yfder, &
823	itap, dtime, jour, knon, ni, &
824	lafin, &
825	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
826	AcoefH, AcoefQ, BcoefH, BcoefQ, &
827	AcoefU, AcoefV, BcoefU, BcoefV, &
828	ypsref, yu1, yv1, yrugoro, pctsrf, &
829	ysnow, yqsurf, yqsol, yagesno, ytsoil, &
830	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
831	ytsurf_new, y_dflux_t, y_dflux_q, &
832	y_flux_u1, y_flux_v1)
833
834
835	CASE DEFAULT
836	WRITE(lunout,*) 'Surface index = ', nsrf
837	abort_message = 'Surface index not valid'
838	CALL abort_gcm(modname,abort_message,1)
839	END SELECT
840
841
842	!****************************************************************************************
843	! 11) - Calcul the increment of surface temperature
844	!
845	!****************************************************************************************
846	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
847
848	!****************************************************************************************
849	!
850	! 12) "La remontee" - "The uphill"
851	!
852	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
853	! for X=H, Q, U and V, for all vertical levels.
854	!
855	!****************************************************************************************
856	! H and Q
857	IF (ok_flux_surf) THEN
858	PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
859	y_flux_t1(:) = fsens
860	y_flux_q1(:) = flat/RLVTT
861	yfluxlat(:) = flat
862	ELSE
863	y_flux_t1(:) = yfluxsens(:)
864	y_flux_q1(:) = -yevap(:)
865	ENDIF
866
867	CALL climb_hq_up(knon, dtime, yt, yq, &
868	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
869	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
870
871
872	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
873	y_flux_u, y_flux_v, y_d_u, y_d_v)
874
875
876	DO j = 1, knon
877	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
878	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
879	ENDDO
880
881	!****************************************************************************************
882	! 13) Transform variables for output format :
883	! - Decompress
884	! - Multiply with pourcentage of current surface
885	! - Cumulate in global variable
886	!
887	!****************************************************************************************
888
889	tke(:,:,nsrf) = 0.
890	DO k = 1, klev
891	DO j = 1, knon
892	i = ni(j)
893	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
894	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
895	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
896	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
897
898	flux_t(i,k,nsrf) = y_flux_t(j,k)
899	flux_q(i,k,nsrf) = y_flux_q(j,k)
900	flux_u(i,k,nsrf) = y_flux_u(j,k)
901	flux_v(i,k,nsrf) = y_flux_v(j,k)
902
903	tke(i,k,nsrf) = ytke(j,k)
904
905	ENDDO
906	ENDDO
907
908	evap(:,nsrf) = - flux_q(:,1,nsrf)
909
910	alb1(:, nsrf) = 0.
911	alb2(:, nsrf) = 0.
912	snow(:, nsrf) = 0.
913	qsurf(:, nsrf) = 0.
914	rugos(:, nsrf) = 0.
915	fluxlat(:,nsrf) = 0.
916	DO j = 1, knon
917	i = ni(j)
918	d_ts(i,nsrf) = y_d_ts(j)
919	alb1(i,nsrf) = yalb1_new(j)
920	alb2(i,nsrf) = yalb2_new(j)
921	snow(i,nsrf) = ysnow(j)
922	qsurf(i,nsrf) = yqsurf(j)
923	rugos(i,nsrf) = yz0_new(j)
924	fluxlat(i,nsrf) = yfluxlat(j)
925	agesno(i,nsrf) = yagesno(j)
926	cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
927	cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
928	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
929	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
930	END DO
931
932	DO k = 2, klev
933	DO j = 1, knon
934	i = ni(j)
935	zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)*ypct(j)
936	END DO
937	END DO
938
939	IF ( nsrf .EQ. is_ter ) THEN
940	DO j = 1, knon
941	i = ni(j)
942	qsol(i) = yqsol(j)
943	END DO
944	END IF
945
946	ftsoil(:,:,nsrf) = 0.
947	DO k = 1, nsoilmx
948	DO j = 1, knon
949	i = ni(j)
950	ftsoil(i, k, nsrf) = ytsoil(j,k)
951	END DO
952	END DO
953
954
955	DO k = 1, klev
956	DO j = 1, knon
957	i = ni(j)
958	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
959	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
960	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
961	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
962	END DO
963	END DO
964
965	!****************************************************************************************
966	! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m
967	! Call HBTM
968	!
969	!****************************************************************************************
970	t2m(:,nsrf) = 0.
971	q2m(:,nsrf) = 0.
972	u10m(:,nsrf) = 0.
973	v10m(:,nsrf) = 0.
974
975	pblh(:,nsrf) = 0. ! Hauteur de couche limite
976	plcl(:,nsrf) = 0. ! Niveau de condensation de la CLA
977	capCL(:,nsrf) = 0. ! CAPE de couche limite
978	oliqCL(:,nsrf) = 0. ! eau_liqu integree de couche limite
979	cteiCL(:,nsrf) = 0. ! cloud top instab. crit. couche limite
980	pblt(:,nsrf) = 0. ! T a la Hauteur de couche limite
981	therm(:,nsrf) = 0.
982	trmb1(:,nsrf) = 0. ! deep_cape
983	trmb2(:,nsrf) = 0. ! inhibition
984	trmb3(:,nsrf) = 0. ! Point Omega
985
986	#undef T2m
987	#define T2m
988	#ifdef T2m
989	! Calculations of diagnostic t,q at 2m and u, v at 10m
990
991	DO j=1, knon
992	i = ni(j)
993	uzon(j) = yu(j,1) + y_d_u(j,1)
994	vmer(j) = yv(j,1) + y_d_v(j,1)
995	tair1(j) = yt(j,1) + y_d_t(j,1)
996	qair1(j) = yq(j,1) + y_d_q(j,1)
997	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
998	* (ypaprs(j,1)-ypplay(j,1))
999	tairsol(j) = yts(j) + y_d_ts(j)
1000	rugo1(j) = yrugos(j)
1001	IF(nsrf.EQ.is_oce) THEN
1002	rugo1(j) = rugos(i,nsrf)
1003	ENDIF
1004	psfce(j)=ypaprs(j,1)
1005	patm(j)=ypplay(j,1)
1006	qairsol(j) = yqsurf(j)
1007	END DO
1008
1009
1010	! Calculate the temperature et relative humidity at 2m and the wind at 10m
1011	CALL stdlevvar(klon, knon, nsrf, zxli, &
1012	uzon, vmer, tair1, qair1, zgeo1, &
1013	tairsol, qairsol, rugo1, psfce, patm, &
1014	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
1015
1016	DO j=1, knon
1017	i = ni(j)
1018	t2m(i,nsrf)=yt2m(j)
1019	q2m(i,nsrf)=yq2m(j)
1020
1021	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
1022	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
1023	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
1024	END DO
1025
1026	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
1027	!IM Ajoute dependance type surface
1028	IF (thermcep) THEN
1029	DO j = 1, knon
1030	i=ni(j)
1031	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
1032	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
1033	zx_qs1 = MIN(0.5,zx_qs1)
1034	zcor1 = 1./(1.-RETV*zx_qs1)
1035	zx_qs1 = zx_qs1*zcor1
1036
1037	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
1038	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
1039	END DO
1040	END IF
1041
1042	CALL HBTM(knon, ypaprs, ypplay, &
1043	yt2m,yt10m,yq2m,yq10m,yustar, &
1044	y_flux_t,y_flux_q,yu,yv,yt,yq, &
1045	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
1046	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
1047
1048	DO j=1, knon
1049	i = ni(j)
1050	pblh(i,nsrf) = ypblh(j)
1051	plcl(i,nsrf) = ylcl(j)
1052	capCL(i,nsrf) = ycapCL(j)
1053	oliqCL(i,nsrf) = yoliqCL(j)
1054	cteiCL(i,nsrf) = ycteiCL(j)
1055	pblT(i,nsrf) = ypblT(j)
1056	therm(i,nsrf) = ytherm(j)
1057	trmb1(i,nsrf) = ytrmb1(j)
1058	trmb2(i,nsrf) = ytrmb2(j)
1059	trmb3(i,nsrf) = ytrmb3(j)
1060	END DO
1061
1062	#else
1063	! T2m not defined
1064	! No calculation
1065	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
1066	#endif
1067
1068	!****************************************************************************************
1069	! 15) End of loop over different surfaces
1070	!
1071	!****************************************************************************************
1072	END DO loop_nbsrf
1073
1074	!****************************************************************************************
1075	! 16) Calculate the mean value over all sub-surfaces for som variables
1076	!
1077	!****************************************************************************************
1078
1079	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
1080	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
1081	DO nsrf = 1, nbsrf
1082	DO k = 1, klev
1083	DO i = 1, klon
1084	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
1085	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
1086	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
1087	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
1088	END DO
1089	END DO
1090	END DO
1091
1092	DO i = 1, klon
1093	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1094	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1095	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
1096	ENDDO
1097
1098	!
1099	! Incrementer la temperature du sol
1100	!
1101	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
1102	zt2m(:) = 0.0 ; zq2m(:) = 0.0
1103	zu10m(:) = 0.0 ; zv10m(:) = 0.0
1104	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
1105	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
1106	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
1107	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
1108	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
1109
1110
1111	DO nsrf = 1, nbsrf
1112	DO i = 1, klon
1113	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
1114
1115	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
1116	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1117	wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
1118	pctsrf(i,nsrf)
1119
1120	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
1121	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
1122
1123	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
1124	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
1125	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
1126	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
1127
1128	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
1129	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
1130	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
1131	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
1132	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
1133	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
1134	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
1135	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
1136	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
1137	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
1138	END DO
1139	END DO
1140
1141	IF (check) THEN
1142	amn=MIN(ts(1,is_ter),1000.)
1143	amx=MAX(ts(1,is_ter),-1000.)
1144	DO i=2, klon
1145	amn=MIN(ts(i,is_ter),amn)
1146	amx=MAX(ts(i,is_ter),amx)
1147	ENDDO
1148	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
1149	ENDIF
1150
1151	!jg ?
1152	!!$!
1153	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
1154	!!$! sub-surfaces is distributed.
1155	!!$!
1156	!!$ DO nsrf = 1, nbsrf
1157	!!$ DO i = 1, klon
1158	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
1159	!!$ ts(i,nsrf) = zxtsol(i)
1160	!!$ t2m(i,nsrf) = zt2m(i)
1161	!!$ q2m(i,nsrf) = zq2m(i)
1162	!!$ u10m(i,nsrf) = zu10m(i)
1163	!!$ v10m(i,nsrf) = zv10m(i)
1164	!!$
1165	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
1166	!!$ pblh(i,nsrf) = s_pblh(i)
1167	!!$ plcl(i,nsrf) = s_plcl(i)
1168	!!$ capCL(i,nsrf) = s_capCL(i)
1169	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
1170	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
1171	!!$ pblT(i,nsrf) = s_pblT(i)
1172	!!$ therm(i,nsrf) = s_therm(i)
1173	!!$ trmb1(i,nsrf) = s_trmb1(i)
1174	!!$ trmb2(i,nsrf) = s_trmb2(i)
1175	!!$ trmb3(i,nsrf) = s_trmb3(i)
1176	!!$ ENDIF
1177	!!$ ENDDO
1178	!!$ ENDDO
1179
1180
1181	DO i = 1, klon
1182	fder(i) = - 4.0RSIGMAzxtsol(i)**3
1183	ENDDO
1184
1185	zxqsurf(:) = 0.0
1186	zxsnow(:) = 0.0
1187	DO nsrf = 1, nbsrf
1188	DO i = 1, klon
1189	zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
1190	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
1191	END DO
1192	END DO
1193
1194	! Premier niveau de vent sortie dans physiq.F
1195	zu1(:) = u(:,1)
1196	zv1(:) = v(:,1)
1197
1198	! Some of the module declared variables are returned for printing in physiq.F
1199	qsol_d(:) = qsol(:)
1200	evap_d(:,:) = evap(:,:)
1201	rugos_d(:,:) = rugos(:,:)
1202	agesno_d(:,:) = agesno(:,:)
1203
1204
1205	END SUBROUTINE pbl_surface
1206	!
1207	!****************************************************************************************
1208	!
1209	SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
1210	evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
1211
1212	INCLUDE "indicesol.h"
1213	INCLUDE "dimsoil.h"
1214
1215	! Ouput variables
1216	!****************************************************************************************
1217	REAL, DIMENSION(klon), INTENT(OUT) :: qsol_rst
1218	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
1219	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
1220	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
1221	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_rst
1222	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_rst
1223	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_rst
1224	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
1225
1226
1227	!****************************************************************************************
1228	! Return module variables for writing to restart file
1229	!
1230	!****************************************************************************************
1231	qsol_rst(:) = qsol(:)
1232	fder_rst(:) = fder(:)
1233	snow_rst(:,:) = snow(:,:)
1234	qsurf_rst(:,:) = qsurf(:,:)
1235	evap_rst(:,:) = evap(:,:)
1236	rugos_rst(:,:) = rugos(:,:)
1237	agesno_rst(:,:) = agesno(:,:)
1238	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
1239
1240	!****************************************************************************************
1241	! Deallocate module variables
1242	!
1243	!****************************************************************************************
1244	DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
1245
1246	END SUBROUTINE pbl_surface_final
1247	!
1248	!****************************************************************************************
1249	!
1250	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, tsurf, alb1, alb2, u10m, v10m, tke)
1251
1252	! Give default values where new fraction has appread
1253
1254	INCLUDE "indicesol.h"
1255	INCLUDE "dimsoil.h"
1256	INCLUDE "clesphys.h"
1257
1258	! Input variables
1259	!****************************************************************************************
1260	INTEGER, INTENT(IN) :: itime
1261	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
1262
1263	! InOutput variables
1264	!****************************************************************************************
1265	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
1266	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: alb1, alb2
1267	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: u10m, v10m
1268	REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke
1269
1270	! Local variables
1271	!****************************************************************************************
1272	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
1273	CHARACTER(len=80) :: abort_message
1274	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
1275	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
1276	!
1277	! All at once !!
1278	!****************************************************************************************
1279
1280	DO nsrf = 1, nbsrf
1281	! First decide complement sub-surfaces
1282	SELECT CASE (nsrf)
1283	CASE(is_oce)
1284	nsrf_comp1=is_sic
1285	nsrf_comp2=is_ter
1286	nsrf_comp3=is_lic
1287	CASE(is_sic)
1288	nsrf_comp1=is_oce
1289	nsrf_comp2=is_ter
1290	nsrf_comp3=is_lic
1291	CASE(is_ter)
1292	nsrf_comp1=is_lic
1293	nsrf_comp2=is_oce
1294	nsrf_comp3=is_sic
1295	CASE(is_lic)
1296	nsrf_comp1=is_ter
1297	nsrf_comp2=is_oce
1298	nsrf_comp3=is_sic
1299	END SELECT
1300
1301	! Initialize all new fractions
1302	DO i=1, klon
1303	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
1304
1305	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
1306	! Use the complement sub-surface, keeping the continents unchanged
1307	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
1308	evap(i,nsrf) = evap(i,nsrf_comp1)
1309	rugos(i,nsrf) = rugos(i,nsrf_comp1)
1310	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
1311	alb1(i,nsrf) = alb1(i,nsrf_comp1)
1312	alb2(i,nsrf) = alb2(i,nsrf_comp1)
1313	u10m(i,nsrf) = u10m(i,nsrf_comp1)
1314	v10m(i,nsrf) = v10m(i,nsrf_comp1)
1315	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
1316	mfois(nsrf) = mfois(nsrf) + 1
1317	ELSE
1318	! The continents have changed. The new fraction receives the mean sum of the existent fractions
1319	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1320	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1321	rugos(i,nsrf) = rugos(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + rugos(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1322	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1323	alb1(i,nsrf) = alb1(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + alb1(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1324	alb2(i,nsrf) = alb2(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + alb2(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1325	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1326	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
1327	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
1328
1329	! Security abort. This option has never been tested. To test, comment the following line.
1330	! abort_message='The fraction of the continents have changed!'
1331	! CALL abort_gcm(modname,abort_message,1)
1332	nfois(nsrf) = nfois(nsrf) + 1
1333	END IF
1334	snow(i,nsrf) = 0.
1335	agesno(i,nsrf) = 0.
1336	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
1337	ELSE
1338	pfois(nsrf) = pfois(nsrf)+ 1
1339	END IF
1340	END DO
1341
1342	END DO
1343
1344	END SUBROUTINE pbl_surface_newfrac
1345
1346	!
1347	!****************************************************************************************
1348	!
1349
1350	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format