Context Navigation

pbl_surface_mod.F90 @ 1146

Last change on this file since 1146 was 1146, checked in by Laurent Fairhead, 15 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 57.6 KB

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

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

Download in other formats:

Original Format