Juan 8/12/2016: add management of LEN_HREC in MNH & SURFEX
[MNH-git_open_source-lfs.git] / src / SURFEX / prep_teb_extern.F90
1 !SURFEX_LIC Copyright 1994-2014 Meteo-France 
2 !SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C  licence
3 !SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SURFEX_LIC for details. version 1.
5 !!
6 !!    MODIFICATIONS
7 !!    -------------
8 !!      M.Moge    01/2016  using READ_SURF_FIELD2D/3D for 2D/3D surfex fields reads
9 !     #########
10 SUBROUTINE PREP_TEB_EXTERN(HPROGRAM,HSURF,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,KLUOUT,PFIELD)
11 !     #################################################################################
12 !
13 USE MODD_TYPE_DATE_SURF
14 !
15 USE MODI_PREP_GRID_EXTERN
16 USE MODI_READ_SURF
17 USE MODI_GET_TEB_DEPTHS
18 USE MODI_INTERP_GRID
19 USE MODI_OPEN_AUX_IO_SURF
20 USE MODI_CLOSE_AUX_IO_SURF
21 USE MODI_TOWN_PRESENCE
22 USE MODI_READ_TEB_PATCH
23 USE MODI_GET_CURRENT_TEB_PATCH
24 USE MODI_READ_SURF_FIELD2D
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  CHARACTER(LEN=6),   INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
41  CHARACTER(LEN=7),   INTENT(IN)  :: HSURF     ! type of field
42  CHARACTER(LEN=28),  INTENT(IN)  :: HFILE     ! name of file
43  CHARACTER(LEN=6),   INTENT(IN)  :: HFILETYPE ! type of input file
44  CHARACTER(LEN=28),  INTENT(IN)  :: HFILEPGD     ! name of file
45  CHARACTER(LEN=6),   INTENT(IN)  :: HFILEPGDTYPE ! type of input file
46 INTEGER,            INTENT(IN)  :: KLUOUT    ! logical unit of output listing
47 REAL,DIMENSION(:,:), POINTER    :: PFIELD    ! field to interpolate horizontally
48 !
49 !*      0.2    declarations of local variables
50 !
51 REAL, DIMENSION(:,:), ALLOCATABLE :: ZFIELD         ! field read
52 REAL, DIMENSION(:,:), ALLOCATABLE :: ZDEPTH         ! depth of each layer
53 REAL, DIMENSION(:),   ALLOCATABLE :: ZDEPTH_TOT     ! total depth of surface
54 !
55 REAL, DIMENSION(:,:),   ALLOCATABLE :: ZD  ! intermediate array
56 !
57  CHARACTER(LEN=LEN_HREC) :: YRECFM         ! Name of the article to be read
58 INTEGER           :: IRESP          ! reading return code
59 INTEGER           :: ILAYER         ! number of layers
60 INTEGER           :: JLAYER         ! loop counter
61 INTEGER           :: IVERSION       ! SURFEX version
62 INTEGER           :: IBUGFIX        ! SURFEX bug version
63 LOGICAL           :: GOLD_NAME      ! old name flag for temperatures
64  CHARACTER(LEN=4)  :: YWALL_OPT      ! option of walls
65  CHARACTER(LEN=6)  :: YSURF          ! Surface type
66  CHARACTER(LEN=3)  :: YBEM ! key of the building energy model DEF for DEFault (Masson et al. 2002) ,
67                           ! BEM for Building Energy Model (Bueno et al. 2012)
68 !
69 INTEGER           :: INI            ! total 1D dimension
70 !
71 LOGICAL                              :: GTEB      ! flag if TEB fields are present
72 INTEGER                              :: IPATCH    ! number of soil temperature patches
73 INTEGER                              :: ITEB_PATCH! number of TEB patches in file
74 INTEGER                              :: ICURRENT_PATCH! current TEB patch to be initialized
75  CHARACTER(LEN=3)                     :: YPATCH    ! indentificator for TEB patch
76 REAL(KIND=JPRB) :: ZHOOK_HANDLE
77 !-------------------------------------------------------------------------------------
78 !
79 !*      1.     Preparation of IO for reading in the file
80 !              -----------------------------------------
81 !
82 !* Note that all points are read, even those without physical meaning.
83 !  These points will not be used during the horizontal interpolation step.
84 !  Their value must be defined as XUNDEF.
85 !
86 IF (LHOOK) CALL DR_HOOK('PREP_TEB_EXTERN',0,ZHOOK_HANDLE)
87 !
88  CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'TOWN  ')
89 !
90 !* reading of version of the file being read
91  CALL READ_SURF(HFILEPGDTYPE,'VERSION',IVERSION,IRESP)
92  CALL READ_SURF(HFILEPGDTYPE,'BUG',IBUGFIX,IRESP)
93 GOLD_NAME=(IVERSION<7 .OR. (IVERSION==7 .AND. IBUGFIX<3))
94 !
95 IF (.NOT.GOLD_NAME) THEN
96    YRECFM='BEM'
97    CALL READ_SURF(HFILEPGDTYPE,YRECFM,YBEM,IRESP)
98 ELSE
99    YBEM='DEF'
100 ENDIF
101 !-------------------------------------------------------------------------------------
102 !
103 !*      2.     Reading of grid
104 !              ---------------
105 !
106 !* reads the grid
107  CALL PREP_GRID_EXTERN(HFILEPGDTYPE,KLUOUT,CINGRID_TYPE,CINTERP_TYPE,INI)
108 !
109 !
110 !* reads if TEB fields exist in the input file
111  CALL TOWN_PRESENCE(HFILEPGDTYPE,GTEB)
112 !
113 !---------------------------------------------------------------------------------------
114 !
115 !*     3.      Orography
116 !              ---------
117 !
118 IF (HSURF=='ZS     ') THEN
119   !
120   ALLOCATE(PFIELD(INI,1))
121   YRECFM='ZS'
122   CALL READ_SURF(HFILEPGDTYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='A')
123   CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
124   !
125   !---------------------------------------------------------------------------------------
126 ELSE
127 !---------------------------------------------------------------------------------------
128 !
129 !*     4.     TEB fields are read
130 !             -------------------
131 !
132   IF (GTEB) THEN
133 !
134     CALL READ_TEB_PATCH(HFILEPGDTYPE,ITEB_PATCH)
135     CALL GET_CURRENT_TEB_PATCH(ICURRENT_PATCH)
136     YPATCH='   '
137     IF (ITEB_PATCH>1) THEN
138       WRITE(YPATCH,FMT='(A,I1,A)') 'T',MIN(ICURRENT_PATCH,ITEB_PATCH),'_'
139     END IF
140 !---------------------------------------------------------------------------------------
141     SELECT CASE(HSURF)
142 !---------------------------------------------------------------------------------------
143 !
144 !*     4.1    Profile of temperatures in roads, roofs or walls
145 !             ------------------------------------------------
146 !
147     CASE('T_ROAD','T_ROOF','T_WALLA','T_WALLB','T_FLOOR','T_MASS')
148       YSURF=HSURF(1:6)
149       !* reading of number of layers
150       IF (YSURF=='T_ROAD') YRECFM='ROAD_LAYER'
151       IF (YSURF=='T_ROOF') YRECFM='ROOF_LAYER'
152       IF (YSURF=='T_WALL') YRECFM='WALL_LAYER'
153       IF (YSURF=='T_FLOO' .OR. YSURF=='T_MASS') THEN 
154         IF (YBEM=='DEF') THEN
155           YRECFM='ROAD_LAYER'
156         ELSE
157           YRECFM='FLOOR_LAYER'
158         END IF
159       END IF
160       CALL READ_SURF(HFILEPGDTYPE,YRECFM,ILAYER,IRESP)
161       !
162       ALLOCATE(ZD(INI,ILAYER))
163       IF (YSURF=='T_ROAD') CALL GET_TEB_DEPTHS(HFILEPGDTYPE,PD_ROAD=ZD)
164       IF (YSURF=='T_ROOF') CALL GET_TEB_DEPTHS(HFILEPGDTYPE,PD_ROOF=ZD)
165       IF (YSURF=='T_WALL') CALL GET_TEB_DEPTHS(HFILEPGDTYPE,PD_WALL=ZD)
166       IF (YSURF=='T_MASS') CALL GET_TEB_DEPTHS(HFILEPGDTYPE,PD_FLOOR=ZD)
167       IF (YSURF=='T_FLOO') CALL GET_TEB_DEPTHS(HFILEPGDTYPE,PD_FLOOR=ZD)
168       !
169       CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
170       CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'TOWN  ')
171       !
172       !* reading option for road orientation
173       YWALL_OPT = 'UNIF'
174       IF (YSURF =='T_WALL' .AND. .NOT. GOLD_NAME) THEN
175         CALL READ_SURF(HFILETYPE,'WALL_OPT',YWALL_OPT,IRESP)
176       END IF
177       !
178       !* reading of the profile
179       ALLOCATE(ZFIELD(INI,ILAYER))
180       DO JLAYER=1,ILAYER
181         IF (GOLD_NAME) THEN
182           WRITE(YRECFM,'(A6,I1.1)') HSURF(1:6),JLAYER
183         ELSE
184           WRITE(YRECFM,'(A1,A4,I1.1)') HSURF(1:1),HSURF(3:6),JLAYER
185           IF (YSURF =='T_WALL' .AND. YWALL_OPT/='UNIF') &
186             WRITE(YRECFM,'(A1,A5,I1.1)') HSURF(1:1),HSURF(3:7),JLAYER
187           IF ((HSURF=='T_FLOOR' .OR. HSURF=='T_MASS') .AND. YBEM=='DEF') THEN
188             IF (HSURF=='T_FLOOR' .AND. JLAYER>1) THEN 
189               WRITE(YRECFM,'(A5,I1.1)') 'TROAD',JLAYER
190             ELSE
191               WRITE(YRECFM,'(A6)') 'TI_BLD'
192             ENDIF
193           END IF
194         END IF
195         YRECFM=YPATCH//YRECFM
196         YRECFM=ADJUSTL(YRECFM)
197         CALL READ_SURF(HFILETYPE,YRECFM,ZFIELD(:,JLAYER),IRESP,HDIR='A')
198       END DO
199       CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
200       !
201       !* recovers middle layer depth (from the surface)
202       ALLOCATE(ZDEPTH    (INI,ILAYER))
203       ALLOCATE(ZDEPTH_TOT(INI))
204       ZDEPTH    (:,1)=ZD(:,1)/2.
205       ZDEPTH_TOT(:)  =ZD(:,1)
206       DO JLAYER=2,ILAYER
207         ZDEPTH    (:,JLAYER) = ZDEPTH_TOT(:) + ZD(:,JLAYER)/2.
208         ZDEPTH_TOT(:) = ZDEPTH_TOT(:) + ZD(:,JLAYER)
209       END DO
210       !
211       !* in case of wall or roof, normalizes by total wall or roof thickness
212       IF (YSURF=='T_ROOF' .OR. YSURF=='T_WALL' .OR. HSURF == 'T_FLOOR' .OR. HSURF == 'T_MASS') THEN
213         DO JLAYER=1,ILAYER
214           ZDEPTH(:,JLAYER) = ZDEPTH(:,JLAYER) / ZDEPTH_TOT(:)
215         END DO
216       END IF
217       !
218       !* interpolation on the fine vertical grid
219       IF (YSURF=='T_ROAD') THEN
220         ALLOCATE(PFIELD(SIZE(ZFIELD,1),SIZE(XGRID_ROAD)))
221         CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_ROAD,PFIELD)
222       ELSEIF (YSURF=='T_ROOF') THEN
223         ALLOCATE(PFIELD(SIZE(ZFIELD,1),SIZE(XGRID_ROOF)))
224         CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_ROOF,PFIELD)
225       ELSEIF (YSURF=='T_WALL') THEN
226         ALLOCATE(PFIELD(SIZE(ZFIELD,1),SIZE(XGRID_WALL)))
227         CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_WALL,PFIELD)
228       ELSEIF (YSURF=='T_FLOO' .OR. YSURF=='T_MASS') THEN
229         ALLOCATE(PFIELD(SIZE(ZFIELD,1),SIZE(XGRID_FLOOR)))
230         CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_FLOOR,PFIELD)
231       END IF
232       !
233       !* end
234       DEALLOCATE(ZD)
235       DEALLOCATE(ZFIELD)
236       DEALLOCATE(ZDEPTH)
237       DEALLOCATE(ZDEPTH_TOT)
238 !---------------------------------------------------------------------------------------
239 !
240 !*      4.2    Internal moisture
241 !              ---------------
242 !
243     CASE('QI_BLD ')
244       ALLOCATE(PFIELD(INI,1))
245       IF (YBEM=='BEM') THEN
246         YRECFM='QI_BLD'
247         YRECFM=YPATCH//YRECFM
248         YRECFM=ADJUSTL(YRECFM)
249         CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
250         CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'TOWN  ')
251         CALL READ_SURF(HFILETYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='A')
252         CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
253       ELSE
254         PFIELD(:,1) = XUNDEF
255       ENDIF
256 !
257 !---------------------------------------------------------------------------------------
258 !
259 !*      4.2    Other variables
260 !              ---------------
261 !
262     CASE DEFAULT
263       ALLOCATE(PFIELD(INI,1))
264       YRECFM=HSURF
265       IF (HSURF=='T_CAN  ') THEN
266         YRECFM='TCANYON'
267         IF (GOLD_NAME) YRECFM='T_CANYON'
268       ELSEIF (HSURF=='Q_CAN  ') THEN
269         YRECFM='QCANYON'
270         IF (GOLD_NAME) YRECFM='Q_CANYON'
271       ELSEIF (HSURF=='T_WIN2 ' .OR. HSURF=='T_WIN1') THEN
272         IF (YBEM=='BEM') THEN
273           YRECFM=HSURF
274         ELSE
275           YRECFM='TI_BLD'
276         ENDIF
277       ENDIF
278       YRECFM=YPATCH//YRECFM
279       YRECFM=ADJUSTL(YRECFM)
280       CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
281       CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'TOWN  ')
282       CALL READ_SURF(HFILETYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='A')
283       CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
284 !
285 !---------------------------------------------------------------------------------------
286     END SELECT
287 !---------------------------------------------------------------------------------------
288 !
289 !*     5.     Subtitutes if TEB fields do not exist
290 !             -------------------------------------
291 !
292   ELSE
293
294     CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
295
296     SELECT CASE(HSURF)
297
298     !* temperature profiles
299     CASE('T_ROAD','T_ROOF','T_WALL','T_WIN1','T_FLOOR','T_CAN','TI_ROAD')
300       YSURF=HSURF(1:6)
301       !* reading of the soil surface temperature
302       CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'NATURE')
303       CALL READ_SURF(HFILEPGDTYPE,'PATCH_NUMBER',IPATCH,IRESP)
304       CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
305       ALLOCATE(ZFIELD(INI,IPATCH))
306       CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'NATURE')
307       IF (YSURF=='T_FLOO' .OR. YSURF=='T_CAN ' .OR. YSURF=='TI_ROA') THEN
308         YRECFM='TG2'
309         CALL READ_SURF_FIELD2D(HFILETYPE,ZFIELD(:,:),YRECFM,HDIR='A')
310       ELSE
311         YRECFM='TG1'
312         CALL READ_SURF_FIELD2D(HFILETYPE,ZFIELD(:,:),YRECFM,HDIR='A')
313       ENDIF
314       CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
315       !* fills the whole temperature profile by this soil temperature
316       IF (YSURF=='T_ROAD') ILAYER=SIZE(XGRID_ROAD)
317       IF (YSURF=='T_ROOF') ILAYER=SIZE(XGRID_ROOF)
318       IF (YSURF=='T_WALL') ILAYER=SIZE(XGRID_WALL)
319       IF (YSURF=='T_FLOO') ILAYER=SIZE(XGRID_FLOOR)
320       IF (YSURF=='T_WIN1' .OR. YSURF=='T_CAN ' .OR. YSURF=='TI_ROA') ILAYER=1
321       ALLOCATE(PFIELD(INI,ILAYER))
322       IF (YSURF=='T_FLOO') THEN
323         !* sets the temperature equal to this deep soil temperature
324         PFIELD(:,1) = XTI_BLD_DEF
325       ELSE
326         PFIELD(:,1) = ZFIELD(:,1)
327       ENDIF
328       DO JLAYER=2,ILAYER
329         PFIELD(:,JLAYER) = ZFIELD(:,1)
330       END DO
331       DEALLOCATE(ZFIELD)
332
333     CASE('T_MASS','TI_BLD','T_WIN2')
334       YSURF=HSURF(1:6)
335       IF (YSURF=='T_MASS') ILAYER = SIZE(XGRID_FLOOR)
336       IF (YSURF=='TI_BLD'.OR.YSURF=='T_WIN2') ILAYER=1
337       ALLOCATE(PFIELD(INI, ILAYER))
338       PFIELD(:,:) = XTI_BLD_DEF
339  
340     !* building moisture
341     CASE('QI_BLD ')
342       ALLOCATE(PFIELD(INI,1))
343       PFIELD(:,1) = XUNDEF
344
345     !* water reservoirs
346     CASE('WS_ROOF','WS_ROAD')
347       ALLOCATE(PFIELD(INI,1))
348       IF (HSURF=='WS_ROOF') PFIELD = XWS_ROOF_DEF
349       IF (HSURF=='WS_ROAD') PFIELD = XWS_ROAD_DEF
350
351    !* other fields
352     CASE DEFAULT
353       ALLOCATE(PFIELD(INI,1))
354       PFIELD = 0.
355
356     END SELECT
357
358   END IF
359 !-------------------------------------------------------------------------------------
360 END IF
361 !-------------------------------------------------------------------------------------
362 !
363 !*      6.     End of IO
364 !              ---------
365 !
366 IF (LHOOK) CALL DR_HOOK('PREP_TEB_EXTERN',1,ZHOOK_HANDLE)
367 !
368 !---------------------------------------------------------------------------------------
369 !
370 END SUBROUTINE PREP_TEB_EXTERN