Context Navigation

pbl_surface_mod.F90 @ 1040

Last change on this file since 1040 was 1009, checked in by lsce, 16 years ago

Changed variable name OCEAN to type_ocean in physiq.def
Removed test of coherence between CPP_VEGET and ok_veget (later to be readded in surf_land_orchdiee)

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

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

Download in other formats:

Original Format