SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
prep_teb_extern.F90
Go to the documentation of this file.
1 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SFX_LIC for details. version 1.
5 ! #########
6 SUBROUTINE prep_teb_extern (DTCO, &
7  hprogram,hsurf,hfile,hfiletype,hfilepgd,hfilepgdtype,kluout,kpatch,pfield)
8 ! #################################################################################
9 !
10 !
11 !
12 !
14 !
16 !
17 USE modi_prep_grid_extern
19 USE modi_get_teb_depths
21 USE modi_open_aux_io_surf
22 USE modi_close_aux_io_surf
23 USE modi_town_presence
24 USE modi_read_teb_patch
25 !
26 USE modd_prep, ONLY : cingrid_type, cinterp_type
27 USE modd_prep_teb, ONLY : xgrid_road, xgrid_wall, xgrid_roof, &
28  xgrid_floor, xws_roof, xws_road, &
29  xti_bld_def, xws_roof_def, xws_road_def, xhui_bld_def
30 USE modd_data_cover_par, ONLY : jpcover
31 USE modd_surf_par, ONLY: xundef
32 !
33 USE yomhook ,ONLY : lhook, dr_hook
34 USE parkind1 ,ONLY : jprb
35 !
36 IMPLICIT NONE
37 !
38 !* 0.1 declarations of arguments
39 !
40 !
41 TYPE(data_cover_t), INTENT(INOUT) :: dtco
42 !
43  CHARACTER(LEN=6), INTENT(IN) :: hprogram ! program calling surf. schemes
44  CHARACTER(LEN=7), INTENT(IN) :: hsurf ! type of field
45  CHARACTER(LEN=28), INTENT(IN) :: hfile ! name of file
46  CHARACTER(LEN=6), INTENT(IN) :: hfiletype ! type of input file
47  CHARACTER(LEN=28), INTENT(IN) :: hfilepgd ! name of file
48  CHARACTER(LEN=6), INTENT(IN) :: hfilepgdtype ! type of input file
49 INTEGER, INTENT(IN) :: kluout ! logical unit of output listing
50 INTEGER, INTENT(IN) :: kpatch
51 REAL,DIMENSION(:,:), POINTER :: pfield ! field to interpolate horizontally
52 !
53 !* 0.2 declarations of local variables
54 !
55 REAL, DIMENSION(:,:), ALLOCATABLE :: zfield ! field read
56 REAL, DIMENSION(:,:), ALLOCATABLE :: zdepth ! depth of each layer
57 REAL, DIMENSION(:), ALLOCATABLE :: zdepth_tot ! total depth of surface
58 !
59 REAL, DIMENSION(:,:), ALLOCATABLE :: zd ! intermediate array
60 !
61 REAL, DIMENSION(:), ALLOCATABLE :: zmask
62 !
63  CHARACTER(LEN=12) :: yrecfm ! Name of the article to be read
64 INTEGER :: iresp ! reading return code
65 INTEGER :: ilayer ! number of layers
66 INTEGER :: jlayer ! loop counter
67 INTEGER :: iversion ! SURFEX version
68 INTEGER :: ibugfix ! SURFEX bug version
69 LOGICAL :: gold_name ! old name flag for temperatures
70  CHARACTER(LEN=4) :: ywall_opt ! option of walls
71  CHARACTER(LEN=6) :: ysurf ! Surface type
72  CHARACTER(LEN=3) :: ybem ! key of the building energy model DEF for DEFault (Masson et al. 2002) ,
73  ! BEM for Building Energy Model (Bueno et al. 2012)
74 !
75 INTEGER :: ini ! total 1D dimension
76 !
77 LOGICAL :: gteb ! flag if TEB fields are present
78 INTEGER :: ipatch ! number of soil temperature patches
79 INTEGER :: iteb_patch! number of TEB patches in file
80  CHARACTER(LEN=3) :: ypatch ! indentificator for TEB patch
81 REAL(KIND=JPRB) :: zhook_handle
82 !-------------------------------------------------------------------------------------
83 !
84 !* 1. Preparation of IO for reading in the file
85 ! -----------------------------------------
86 !
87 !* Note that all points are read, even those without physical meaning.
88 ! These points will not be used during the horizontal interpolation step.
89 ! Their value must be defined as XUNDEF.
90 !
91 IF (lhook) CALL dr_hook('PREP_TEB_EXTERN',0,zhook_handle)
92 !
93 !
94 !* reading of version of the file being read
95  CALL open_aux_io_surf(&
96  hfilepgd,hfilepgdtype,'FULL ')
97  CALL read_surf(&
98  hfilepgdtype,'VERSION',iversion,iresp)
99  CALL read_surf(&
100  hfilepgdtype,'BUG',ibugfix,iresp)
101  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
102 gold_name=(iversion<7 .OR. (iversion==7 .AND. ibugfix<3))
103 !
104 !-------------------------------------------------------------------------------------
105 !
106 !* 2. Reading of grid
107 ! ---------------
108 !
109  CALL open_aux_io_surf(&
110  hfilepgd,hfilepgdtype,'FULL ')
111 !* reads the grid
112  CALL prep_grid_extern(&
113  hfilepgdtype,kluout,cingrid_type,cinterp_type,ini)
114 !
115 yrecfm='VERSION'
116  CALL read_surf(hfilepgdtype,yrecfm,iversion,iresp)
117 !
118 ALLOCATE(zmask(ini))
119 IF (iversion>=7) THEN
120  yrecfm='FRAC_TOWN'
121  CALL read_surf(hfilepgdtype,yrecfm,zmask,iresp,hdir='A')
122 ELSE
123  zmask(:) = 1.
124 ENDIF
125 !
126 !* reads if TEB fields exist in the input file
127  CALL town_presence(&
128  hfilepgdtype,gteb)
129  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
130 !
131 IF (.NOT.gold_name.AND.gteb) THEN
132  yrecfm='BEM'
133  CALL open_aux_io_surf(&
134  hfilepgd,hfilepgdtype,'TOWN ')
135  CALL read_surf(&
136  hfilepgdtype,yrecfm,ybem,iresp)
137  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
138 ELSE
139  ybem='DEF'
140 ENDIF
141 !---------------------------------------------------------------------------------------
142 !
143 !* 3. Orography
144 ! ---------
145 !
146 IF (hsurf=='ZS ') THEN
147  !
148  ALLOCATE(pfield(ini,1))
149  yrecfm='ZS'
150  CALL open_aux_io_surf(&
151  hfilepgd,hfilepgdtype,'FULL ')
152  CALL read_surf(&
153  hfilepgdtype,yrecfm,pfield(:,1),iresp,hdir='A')
154  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
155  !
156  !---------------------------------------------------------------------------------------
157 ELSE
158 !---------------------------------------------------------------------------------------
159 !
160 !* 4. TEB fields are read
161 ! -------------------
162 !
163  IF (gteb) THEN
164 !
165  CALL read_teb_patch(&
166  hfilepgd,hfilepgdtype,iteb_patch)
167  ypatch=' '
168  IF (iteb_patch>1) THEN
169  WRITE(ypatch,fmt='(A,I1,A)') 'T',min(kpatch,iteb_patch),'_'
170  END IF
171 !
172  CALL open_aux_io_surf(&
173  hfilepgd,hfilepgdtype,'TOWN ')
174 !
175 !---------------------------------------------------------------------------------------
176  SELECT CASE(hsurf)
177 !---------------------------------------------------------------------------------------
178 !
179 !* 4.1 Profile of temperatures in roads, roofs or walls
180 ! ------------------------------------------------
181 !
182  CASE('T_ROAD','T_ROOF','T_WALLA','T_WALLB','T_FLOOR','T_MASS')
183  ysurf=hsurf(1:6)
184  !* reading of number of layers
185  IF (ysurf=='T_ROAD') yrecfm='ROAD_LAYER'
186  IF (ysurf=='T_ROOF') yrecfm='ROOF_LAYER'
187  IF (ysurf=='T_WALL') yrecfm='WALL_LAYER'
188  IF (ysurf=='T_WALLA') yrecfm='WALL_LAYER'
189  IF (ysurf=='T_WALLB') yrecfm='WALL_LAYER'
190  IF (ysurf=='T_FLOO' .OR. ysurf=='T_MASS') THEN
191  IF (ybem=='DEF') THEN
192  yrecfm='ROAD_LAYER'
193  ELSE
194  yrecfm='FLOOR_LAYER'
195  END IF
196  END IF
197  CALL read_surf(&
198  hfilepgdtype,yrecfm,ilayer,iresp)
199  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
200  !
201  ALLOCATE(zd(ini,ilayer))
202  IF (ysurf=='T_ROAD') CALL get_teb_depths(&
203  dtco, &
204  hfilepgd,hfilepgdtype,pd_road=zd)
205  IF (ysurf=='T_ROOF') CALL get_teb_depths(&
206  dtco, &
207  hfilepgd,hfilepgdtype,pd_roof=zd)
208  IF (ysurf=='T_WALL') CALL get_teb_depths(&
209  dtco, &
210  hfilepgd,hfilepgdtype,pd_wall=zd)
211  IF (ysurf=='T_WALLA') CALL get_teb_depths(&
212  dtco, &
213  hfilepgd,hfilepgdtype,pd_wall=zd)
214  IF (ysurf=='T_WALLB') CALL get_teb_depths(&
215  dtco, &
216  hfilepgd,hfilepgdtype,pd_wall=zd)
217  IF (ysurf=='T_MASS') CALL get_teb_depths(&
218  dtco, &
219  hfilepgd,hfilepgdtype,pd_floor=zd)
220  IF (ysurf=='T_FLOO') CALL get_teb_depths(&
221  dtco, &
222  hfilepgd,hfilepgdtype,pd_floor=zd)
223  !
224  CALL open_aux_io_surf(&
225  hfile,hfiletype,'TOWN ')
226  !
227  !* reading option for road orientation
228  ywall_opt = 'UNIF'
229  IF (ysurf =='T_WALL' .AND. .NOT. gold_name) THEN
230  CALL read_surf(&
231  hfiletype,'WALL_OPT',ywall_opt,iresp)
232  END IF
233  IF (ysurf =='T_WALLA' .AND. .NOT. gold_name) THEN
234  CALL read_surf(&
235  hfiletype,'WALL_OPT',ywall_opt,iresp)
236  END IF
237  IF (ysurf =='T_WALLB' .AND. .NOT. gold_name) THEN
238  CALL read_surf(&
239  hfiletype,'WALL_OPT',ywall_opt,iresp)
240  END IF
241  !
242  !* reading of the profile
243  ALLOCATE(zfield(ini,ilayer))
244  DO jlayer=1,ilayer
245  IF (gold_name) THEN
246  WRITE(yrecfm,'(A6,I1.1)') hsurf(1:6),jlayer
247  ELSE
248  WRITE(yrecfm,'(A1,A4,I1.1)') hsurf(1:1),hsurf(3:6),jlayer
249  IF (ysurf =='T_WALL' .AND. ywall_opt/='UNIF') &
250  WRITE(yrecfm,'(A1,A5,I1.1)') hsurf(1:1),hsurf(3:7),jlayer
251  IF (ysurf =='T_WALLA' .AND. ywall_opt/='UNIF') &
252  WRITE(yrecfm,'(A1,A5,I1.1)') hsurf(1:1),hsurf(3:7),jlayer
253  IF (ysurf =='T_WALLB' .AND. ywall_opt/='UNIF') &
254  WRITE(yrecfm,'(A1,A5,I1.1)') hsurf(1:1),hsurf(3:7),jlayer
255  IF ((hsurf=='T_FLOOR' .OR. hsurf=='T_MASS') .AND. ybem=='DEF') THEN
256  IF (hsurf=='T_FLOOR' .AND. jlayer>1) THEN
257  WRITE(yrecfm,'(A5,I1.1)') 'TROAD',jlayer
258  ELSE
259  WRITE(yrecfm,'(A6)') 'TI_BLD'
260  ENDIF
261  END IF
262  END IF
263  yrecfm=ypatch//yrecfm
264  yrecfm=adjustl(yrecfm)
265  CALL read_surf(&
266  hfiletype,yrecfm,zfield(:,jlayer),iresp,hdir='A')
267  END DO
268  CALL close_aux_io_surf(hfile,hfiletype)
269  DO jlayer=1,SIZE(zfield,2)
270  WHERE (zmask(:)==0.) zfield(:,jlayer) = xundef
271  ENDDO
272  !
273  !* recovers middle layer depth (from the surface)
274  ALLOCATE(zdepth(ini,ilayer))
275  ALLOCATE(zdepth_tot(ini))
276  zdepth(:,1)=zd(:,1)/2.
277  zdepth_tot(:) =zd(:,1)
278  DO jlayer=2,ilayer
279  zdepth(:,jlayer) = zdepth_tot(:) + zd(:,jlayer)/2.
280  zdepth_tot(:) = zdepth_tot(:) + zd(:,jlayer)
281  END DO
282  !
283  !* in case of wall or roof, normalizes by total wall or roof thickness
284  IF (ysurf=='T_ROOF' .OR. ysurf=='T_WALL' .OR. hsurf == 'T_FLOOR' &
285  &.OR. hsurf == 'T_MASS'.OR. ysurf=='T_WALLA' .OR. hsurf == 'T_WALLB') THEN
286  DO jlayer=1,ilayer
287  zdepth(:,jlayer) = zdepth(:,jlayer) / zdepth_tot(:)
288  END DO
289  END IF
290  !
291  !* interpolation on the fine vertical grid
292  IF (ysurf=='T_ROAD') THEN
293  ALLOCATE(pfield(SIZE(zfield,1),SIZE(xgrid_road)))
294  CALL interp_grid(zdepth,zfield,xgrid_road,pfield)
295  ELSEIF (ysurf=='T_ROOF') THEN
296  ALLOCATE(pfield(SIZE(zfield,1),SIZE(xgrid_roof)))
297  CALL interp_grid(zdepth,zfield,xgrid_roof,pfield)
298  ELSEIF (ysurf=='T_WALL') THEN
299  ALLOCATE(pfield(SIZE(zfield,1),SIZE(xgrid_wall)))
300  CALL interp_grid(zdepth,zfield,xgrid_wall,pfield)
301  ELSEIF (ysurf=='T_WALLA') THEN
302  ALLOCATE(pfield(SIZE(zfield,1),SIZE(xgrid_wall)))
303  CALL interp_grid(zdepth,zfield,xgrid_wall,pfield)
304  ELSEIF (ysurf=='T_WALLB') THEN
305  ALLOCATE(pfield(SIZE(zfield,1),SIZE(xgrid_wall)))
306  CALL interp_grid(zdepth,zfield,xgrid_wall,pfield)
307  ELSEIF (ysurf=='T_FLOO' .OR. ysurf=='T_MASS') THEN
308  ALLOCATE(pfield(SIZE(zfield,1),SIZE(xgrid_floor)))
309  CALL interp_grid(zdepth,zfield,xgrid_floor,pfield)
310  END IF
311  !
312  !* end
313  DEALLOCATE(zd)
314  DEALLOCATE(zfield)
315  DEALLOCATE(zdepth)
316  DEALLOCATE(zdepth_tot)
317 !---------------------------------------------------------------------------------------
318 !
319 !* 4.2 Internal moisture
320 ! ---------------
321 !
322  CASE('QI_BLD ')
323  ALLOCATE(pfield(ini,1))
324  IF (ybem=='BEM') THEN
325  yrecfm='QI_BLD'
326  yrecfm=ypatch//yrecfm
327  yrecfm=adjustl(yrecfm)
328  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
329  CALL open_aux_io_surf(&
330  hfile,hfiletype,'TOWN ')
331  CALL read_surf(&
332  hfiletype,yrecfm,pfield(:,1),iresp,hdir='A')
333  CALL close_aux_io_surf(hfile,hfiletype)
334  WHERE (zmask(:)==0.) pfield(:,1) = xundef
335  ELSE
336  pfield(:,1) = xundef
337  ENDIF
338 !
339 !---------------------------------------------------------------------------------------
340 !
341 !* 4.2 Other variables
342 ! ---------------
343 !
344  CASE default
345  ALLOCATE(pfield(ini,1))
346  yrecfm=hsurf
347  IF (hsurf=='T_CAN ') THEN
348  yrecfm='TCANYON'
349  IF (gold_name) yrecfm='T_CANYON'
350  ELSEIF (hsurf=='Q_CAN ') THEN
351  yrecfm='QCANYON'
352  IF (gold_name) yrecfm='Q_CANYON'
353  ELSEIF (hsurf=='T_WIN2 ' .OR. hsurf=='T_WIN1') THEN
354  IF (ybem=='BEM') THEN
355  yrecfm=hsurf
356  ELSE
357  yrecfm='TI_BLD'
358  ENDIF
359  ENDIF
360  yrecfm=ypatch//yrecfm
361  yrecfm=adjustl(yrecfm)
362  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
363  CALL open_aux_io_surf(&
364  hfile,hfiletype,'TOWN ')
365  CALL read_surf(&
366  hfiletype,yrecfm,pfield(:,1),iresp,hdir='A')
367  CALL close_aux_io_surf(hfile,hfiletype)
368  WHERE (zmask(:)==0.) pfield(:,1) = xundef
369 !
370 !---------------------------------------------------------------------------------------
371  END SELECT
372 !---------------------------------------------------------------------------------------
373 !
374 !* 5. Subtitutes if TEB fields do not exist
375 ! -------------------------------------
376 !
377  ELSE
378 
379  SELECT CASE(hsurf)
380 
381  !* temperature profiles
382  CASE('T_ROAD','T_ROOF','T_WALL','T_WIN1','T_FLOOR','T_CAN','TI_ROAD','T_WALLA','T_WALLB')
383  ysurf=hsurf(1:6)
384  !* reading of the soil surface temperature
385  CALL open_aux_io_surf(&
386  hfilepgd,hfilepgdtype,'NATURE')
387  CALL read_surf(&
388  hfilepgdtype,'PATCH_NUMBER',ipatch,iresp)
389  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
390  ALLOCATE(zfield(ini,ipatch))
391  CALL open_aux_io_surf(&
392  hfile,hfiletype,'NATURE')
393  IF (ysurf=='T_FLOO' .OR. ysurf=='T_CAN ' .OR. ysurf=='TI_ROA') THEN
394  CALL read_surf(&
395  hfiletype,'TG2',zfield(:,:),iresp,hdir='A')
396  ELSE
397  CALL read_surf(&
398  hfiletype,'TG1',zfield(:,:),iresp,hdir='A')
399  ENDIF
400  CALL close_aux_io_surf(hfile,hfiletype)
401  DO jlayer=1,SIZE(zfield,2)
402  WHERE (zmask(:)==0.) zfield(:,jlayer) = xundef
403  ENDDO
404  !* fills the whole temperature profile by this soil temperature
405  IF (ysurf=='T_ROAD') ilayer=SIZE(xgrid_road)
406  IF (ysurf=='T_ROOF') ilayer=SIZE(xgrid_roof)
407  IF (ysurf=='T_WALL') ilayer=SIZE(xgrid_wall)
408  IF (ysurf=='T_WALLA') ilayer=SIZE(xgrid_wall)
409  IF (ysurf=='T_WALLB') ilayer=SIZE(xgrid_wall)
410  IF (ysurf=='T_FLOO') ilayer=SIZE(xgrid_floor)
411  IF (ysurf=='T_WIN1' .OR. ysurf=='T_CAN ' .OR. ysurf=='TI_ROA') ilayer=1
412  ALLOCATE(pfield(ini,ilayer))
413  IF (ysurf=='T_FLOO') THEN
414  !* sets the temperature equal to this deep soil temperature
415  pfield(:,1) = xti_bld_def
416  ELSE
417  pfield(:,1) = zfield(:,1)
418  ENDIF
419  DO jlayer=2,ilayer
420  pfield(:,jlayer) = zfield(:,1)
421  END DO
422  DEALLOCATE(zfield)
423 
424  CASE('T_MASS','TI_BLD','T_WIN2')
425  ysurf=hsurf(1:6)
426  IF (ysurf=='T_MASS') ilayer = SIZE(xgrid_floor)
427  IF (ysurf=='TI_BLD'.OR.ysurf=='T_WIN2') ilayer=1
428  ALLOCATE(pfield(ini, ilayer))
429  pfield(:,:) = xti_bld_def
430 
431  !* building moisture
432  CASE('QI_BLD ')
433  ALLOCATE(pfield(ini,1))
434  pfield(:,1) = xundef
435 
436  !* water reservoirs
437  CASE('WS_ROOF','WS_ROAD')
438  ALLOCATE(pfield(ini,1))
439  IF (hsurf=='WS_ROOF') pfield = xws_roof_def
440  IF (hsurf=='WS_ROAD') pfield = xws_road_def
441 
442  !* other fields
443  CASE default
444  ALLOCATE(pfield(ini,1))
445  pfield = 0.
446 
447  END SELECT
448 
449  END IF
450 !-------------------------------------------------------------------------------------
451 END IF
452 !-------------------------------------------------------------------------------------
453 !
454 DEALLOCATE(zmask)
455 !
456 !* 6. End of IO
457 ! ---------
458 !
459 IF (lhook) CALL dr_hook('PREP_TEB_EXTERN',1,zhook_handle)
460 !
461 !---------------------------------------------------------------------------------------
462 !
463 END SUBROUTINE prep_teb_extern
subroutine close_aux_io_surf(HFILE, HFILETYPE)
subroutine open_aux_io_surf(HFILE, HFILETYPE, HMASK)
subroutine get_teb_depths(DTCO, HFILEPGD, HFILEPGDTYPE, PD_ROOF, PD_ROAD, PD_WALL, PD_FLOOR)
subroutine read_teb_patch(HFILEPGD, HFILEPGDTYPE, KTEB_PATCH)
subroutine prep_teb_extern(DTCO, HPROGRAM, HSURF, HFILE, HFILETYPE, HFILEPGD, HFILEPGDTYPE, KLUOUT, KPATCH, PFIELD)
subroutine town_presence(HFILETYPE, OTEB)
subroutine prep_grid_extern(HFILETYPE, KLUOUT, HGRIDTYPE, HINTERP_TYPE, KNI)