Context Navigation

pbl_surface_mod.F90 @ 952

Last change on this file since 952 was 889, checked in by Laurent Fairhead, 16 years ago
Modifs pour que l'aquaplanete marche FH LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 54.5 KB

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

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

Download in other formats:

Original Format