Context Navigation

pbl_surface_mod.F90 @ 1064

Last change on this file since 1064 was 1064, checked in by Laurent Fairhead, 16 years ago
Correction: initialisation rh2m, qsat2m en dehors boucle nbsrf IM
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 56.8 KB

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

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

Download in other formats:

Original Format