Context Navigation

pbl_surface_mod.F90 @ 888

Last change on this file since 888 was 888, checked in by Laurent Fairhead, 17 years ago
Modifications sur l'albedo JG LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 54.1 KB

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

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

Download in other formats:

Original Format