Juan 8/12/2016: add management of LEN_HREC in MNH & SURFEX
[MNH-git_open_source-lfs.git] / src / SURFEX / prep_snow_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 SUBROUTINE PREP_SNOW_EXTERN(HPROGRAM,HSURF,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,&
7                             KLUOUT,PFIELD,OSNOW_IDEAL,KLAYER)
8
9 !!     ##########################################################################
10 !
11 !
12 !!****  *PREP_SNOW_EXTERN*  
13 !!
14 !!    PURPOSE
15 !!    -------
16 !       Read and prepare initial snow fields from external files
17 !     
18 !!**  METHOD
19 !!    ------
20 !
21 !!    EXTERNAL
22 !!    --------
23 !!
24 !!    IMPLICIT ARGUMENTS
25 !!    ------------------ 
26 !!
27 !!      
28 !!    REFERENCE
29 !!    ---------
30 !!
31 !!    
32 !!    AUTHOR
33 !!    ------
34 !!         * Meteo-France *
35 !!
36 !!    MODIFICATIONS
37 !!    -------------
38 !!      Original    ?
39 !!       02/2014 E. Martin : cor. for passing from from multilayer to a single layer
40 !!       2014    M.Faivre
41 !-------------------------------------------------------------------------------
42 !
43 !
44 USE MODD_TYPE_SNOW
45 USE MODD_PREP,           ONLY : CINGRID_TYPE, CINTERP_TYPE
46 USE MODD_PREP_SNOW,      ONLY : XGRID_SNOW, NGRID_LEVEL
47 USE MODD_DATA_COVER_PAR, ONLY : NVEGTYPE
48 USE MODD_SURF_PAR,       ONLY : XUNDEF
49 USE MODD_CSTS,           ONLY : XTT
50 !
51 USE MODE_SNOW3L
52 !
53 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
54 USE PARKIND1  ,ONLY : JPRB
55 !
56 USE MODI_TOWN_PRESENCE
57 USE MODI_PUT_ON_ALL_VEGTYPES
58 USE MODI_ABOR1_SFX
59 USE MODI_PREP_GRID_EXTERN
60 USE MODI_OPEN_AUX_IO_SURF
61 USE MODI_CLOSE_AUX_IO_SURF
62 USE MODI_ALLOCATE_GR_SNOW
63 USE MODI_INTERP_GRID
64 USE MODI_READ_GR_SNOW
65 USE MODI_READ_SURF
66 USE MODI_SNOW_T_WLIQ_TO_HEAT
67 USE MODI_GET_CURRENT_TEB_PATCH
68 USE MODI_READ_TEB_PATCH
69 !
70 IMPLICIT NONE
71 !
72 !*      0.1    declarations of arguments
73 !
74  CHARACTER(LEN=6),   INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
75  CHARACTER(LEN=10),  INTENT(IN)  :: HSURF     ! type of field
76  CHARACTER(LEN=28),  INTENT(IN)  :: HFILE     ! name of file
77  CHARACTER(LEN=6),   INTENT(IN)  :: HFILETYPE ! type of file
78  CHARACTER(LEN=28),  INTENT(IN)  :: HFILEPGD     ! name of file
79  CHARACTER(LEN=6),   INTENT(IN)  :: HFILEPGDTYPE ! type of file
80 INTEGER,            INTENT(IN)  :: KLUOUT    ! logical unit of output listing
81 REAL,DIMENSION(:,:,:), POINTER  :: PFIELD    ! field to interpolate horizontally
82 LOGICAL,            INTENT(IN)  :: OSNOW_IDEAL
83 INTEGER,            INTENT(IN)  :: KLAYER
84 !
85 !*      0.2    declarations of local variables
86 !
87 TYPE(SURF_SNOW)                    :: TZSNOW ! snow characteristics
88
89 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZFIELD       ! work field on input snow grid
90 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZFIELD_FINE  ! work field on fine snow grid
91 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTEMP        ! snow temperature
92 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZWLIQ        ! liquid water snow pack content
93 REAL, DIMENSION(:,:),   ALLOCATABLE :: ZD           ! total snow depth
94 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDEPTH       ! thickness of each layer (m)
95 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZGRID        ! normalized input grid
96 !
97 LOGICAL                           :: GTOWN          ! town variables written in the file
98  CHARACTER(LEN=LEN_HREC)                 :: YRECFM         ! record name
99 INTEGER                           :: IRESP          ! error return code
100 INTEGER                           :: IVERSION       ! SURFEX version
101 LOGICAL                           :: GOLD_NAME      ! old name flag 
102 INTEGER                           :: IBUGFIX        ! SURFEX bug version
103 INTEGER                           :: IVEGTYPE       ! actual number of vegtypes
104 INTEGER                           :: JLAYER         ! loop on snow vertical grids
105 INTEGER                           :: INI
106  CHARACTER(LEN=8)                  :: YAREA          ! area treated ('ROOF','ROAD','VEG ')
107  CHARACTER(LEN=3)                  :: YPREFIX        ! prefix to identify patch
108 INTEGER                           :: IPATCH         ! number of input patch
109 INTEGER                           :: ITEB_PATCH     ! number of input patch for TEB
110 INTEGER                           :: ICURRENT_TEB_PATCH ! current patch for TEB
111 INTEGER                           :: JPATCH         ! loop on patch
112  CHARACTER(LEN=6)                  :: YMASK          ! type of tile mask
113 REAL(KIND=JPRB) :: ZHOOK_HANDLE
114 !
115 !-------------------------------------------------------------------------------------
116 !
117 !*      3.     Area being treated
118 !              ------------------
119 !
120 IF (LHOOK) CALL DR_HOOK('PREP_SNOW_EXTERN',0,ZHOOK_HANDLE)
121 !
122 !-------------------------------------------------------------------------------------
123 !
124 YAREA='        '
125 YAREA(1:4) = HSURF(7:10)
126 !
127 IF (YAREA(1:4)=='VEG ') THEN
128   IVEGTYPE = NVEGTYPE
129   YMASK = 'NATURE'
130   YPREFIX = '   '  
131 ELSE
132   IVEGTYPE = 1
133   YMASK    = 'TOWN  '
134   IPATCH   = 1
135   YPREFIX = '   '  
136 END IF
137 !
138 !*      1.     Preparation of IO for reading in the file
139 !              -----------------------------------------
140 !
141 !* Note that all points are read, even those without physical meaning.
142 !  These points will not be used during the horizontal interpolation step.
143 !  Their value must be defined as XUNDEF.
144 !
145  CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,YMASK)
146 !
147 !* reading of version of the file being read
148  CALL READ_SURF(HFILEPGDTYPE,'VERSION',IVERSION,IRESP)
149  CALL READ_SURF(HFILEPGDTYPE,'BUG',IBUGFIX,IRESP)
150 GOLD_NAME=(IVERSION<7 .OR. (IVERSION==7 .AND. IBUGFIX<3))
151 !
152 IF (YAREA(1:4)=='VEG ') THEN
153   YRECFM = 'PATCH_NUMBER'
154   CALL READ_SURF(HFILEPGDTYPE,YRECFM,IPATCH,IRESP)
155 ELSE
156   IF (.NOT.GOLD_NAME) THEN
157      IF (YAREA(1:4)=='ROOF') YAREA(1:4) = 'RF  '
158      IF (YAREA(1:4)=='ROAD') YAREA(1:4) = 'RD  '
159    ENDIF
160   CALL READ_TEB_PATCH(HFILEPGDTYPE,ITEB_PATCH)
161   IF (ITEB_PATCH>1) THEN
162     CALL GET_CURRENT_TEB_PATCH(ICURRENT_TEB_PATCH)
163     WRITE(YPREFIX,FMT='(A,I1,A)') 'T',MIN(ICURRENT_TEB_PATCH,ITEB_PATCH),'_'
164   END IF  
165 END IF
166 !
167 !-------------------------------------------------------------------------------------
168 !
169 !*      2.     Reading of grid
170 !              ---------------
171 !
172  CALL PREP_GRID_EXTERN(HFILEPGDTYPE,KLUOUT,CINGRID_TYPE,CINTERP_TYPE,INI)
173 !
174 !-------------------------------------------------------------------------------------
175 !
176 !*      4.     Reading of snow data
177 !              ---------------------
178 !
179 IF (YAREA(1:2)=='RO' .OR. YAREA(1:2)=='GA' .OR. YAREA(1:2)=='RF' .OR. YAREA(1:2)=='RD') THEN
180   CALL TOWN_PRESENCE(HFILEPGDTYPE,GTOWN)
181   CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
182   IF (.NOT. GTOWN) THEN
183     TZSNOW%SCHEME='1-L'
184     CALL ALLOCATE_GR_SNOW(TZSNOW,INI,IPATCH)
185   ELSE
186     CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,YMASK)
187     CALL READ_GR_SNOW(HFILETYPE,TRIM(YAREA),YPREFIX,INI,IPATCH,TZSNOW,HDIR='A')
188     CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
189   ENDIF
190 ELSE
191   CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
192   CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,YMASK)
193   CALL READ_GR_SNOW(HFILETYPE,TRIM(YAREA),YPREFIX,INI,IPATCH,TZSNOW,HDIR='A')
194   CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
195 ENDIF
196 !
197 IF (TZSNOW%NLAYER.GT.KLAYER) THEN
198   TZSNOW%NLAYER=KLAYER
199 ELSEIF (TZSNOW%NLAYER.LT.KLAYER) THEN
200   CALL ABOR1_SFX("PREP_SNOW_EXTERN: SNOW NLAYER IN EXTERN FILE MUST BE GROWER THAN CURRENT NLAYER")
201 ENDIF
202 !
203 !
204 !-------------------------------------------------------------------------------------
205 !
206 !*      5.     Total snow content
207 !              ------------------
208 !
209 SELECT CASE (HSURF(1:3))
210   CASE ('WWW')
211     IF (OSNOW_IDEAL) THEN
212       ALLOCATE(ZFIELD(INI,KLAYER,IPATCH))
213       ZFIELD(:,:,:) = TZSNOW%WSNOW(:,1:KLAYER,:)
214       ALLOCATE(PFIELD(INI,KLAYER,IVEGTYPE))
215       CALL PUT_ON_ALL_VEGTYPES(INI,KLAYER,IPATCH,IVEGTYPE,ZFIELD,PFIELD)
216     ELSE
217       ALLOCATE(ZFIELD(INI,1,IPATCH))
218       ZFIELD = 0.
219       DO JLAYER=1,TZSNOW%NLAYER
220         ZFIELD(:,1,:) = ZFIELD(:,1,:) + TZSNOW%WSNOW(:,JLAYER,:)
221       END DO
222       ALLOCATE(PFIELD(INI,1,IVEGTYPE))
223       CALL PUT_ON_ALL_VEGTYPES(INI,1,IPATCH,IVEGTYPE,ZFIELD,PFIELD)
224     ENDIF
225     DEALLOCATE(ZFIELD)
226 !
227 !-------------------------------------------------------------------------------------
228 !
229 !*      6.     Albedo
230 !              ------
231 !
232   CASE ('ALB')
233     ALLOCATE(ZFIELD(INI,1,IPATCH))
234     ZFIELD(:,1,:) = TZSNOW%ALB(:,:)
235     !
236     ALLOCATE(PFIELD(INI,1,IVEGTYPE))
237     CALL PUT_ON_ALL_VEGTYPES(INI,1,IPATCH,IVEGTYPE,ZFIELD,PFIELD)
238     DEALLOCATE(ZFIELD)
239 !
240 !-------------------------------------------------------------------------------------
241 !
242 !*      7.     Total depth
243 !              -----------
244 !
245   CASE ('DEP')
246     IF (OSNOW_IDEAL) THEN
247       ALLOCATE(ZDEPTH(INI,KLAYER,IPATCH))
248       ZDEPTH(:,:,:) = TZSNOW%WSNOW(:,1:KLAYER,:)/TZSNOW%RHO(:,1:KLAYER,:)
249       WHERE(TZSNOW%WSNOW(:,1:KLAYER,:)==XUNDEF) ZDEPTH(:,:,:)=XUNDEF
250       ALLOCATE(PFIELD(INI,KLAYER,IVEGTYPE))
251       CALL PUT_ON_ALL_VEGTYPES(INI,KLAYER,IPATCH,IVEGTYPE,ZDEPTH,PFIELD)
252     ELSE
253       ALLOCATE(ZDEPTH(INI,1,IPATCH))
254       ZDEPTH = 0.
255       DO JPATCH=1,IPATCH
256         DO JLAYER=1,TZSNOW%NLAYER
257           ZDEPTH(:,1,JPATCH) = ZDEPTH(:,1,JPATCH) + &
258                                TZSNOW%WSNOW(:,JLAYER,JPATCH)/TZSNOW%RHO(:,JLAYER,JPATCH)
259         END DO
260         WHERE(TZSNOW%WSNOW(:,1,JPATCH)==XUNDEF) ZDEPTH(:,1,JPATCH)=XUNDEF
261       END DO
262       ALLOCATE(PFIELD(INI,1,IVEGTYPE))
263       CALL PUT_ON_ALL_VEGTYPES(INI,1,IPATCH,IVEGTYPE,ZDEPTH,PFIELD)
264     ENDIF
265     DEALLOCATE(ZDEPTH)
266 !
267 !-------------------------------------------------------------------------------------
268 !
269 !*      8.     Density or heat profile
270 !              -----------------------
271 !
272   CASE ('RHO','HEA')
273     ALLOCATE(ZFIELD     (INI,TZSNOW%NLAYER,IPATCH))
274
275     SELECT CASE (TZSNOW%SCHEME)
276       CASE ('D95','1-L','EBA')
277         ALLOCATE(ZFIELD_FINE(INI,NGRID_LEVEL,IPATCH))      
278         !* computes output physical variable
279         IF (HSURF(1:1)=='R') ZFIELD(:,1,:) = TZSNOW%RHO(:,1,:)
280         IF (HSURF(1:1)=='H') THEN
281           ALLOCATE(ZTEMP(INI,TZSNOW%NLAYER,IPATCH))
282           ALLOCATE(ZWLIQ(INI,TZSNOW%NLAYER,IPATCH))
283           IF (TZSNOW%SCHEME=='D95'.OR.TZSNOW%SCHEME=='EBA') ZTEMP(:,1,:) = XTT
284           IF (TZSNOW%SCHEME=='1-L') ZTEMP(:,1,:) = TZSNOW%T(:,1,:)
285           ZWLIQ(:,:,:) = 0.
286           CALL SNOW_T_WLIQ_TO_HEAT(ZFIELD,TZSNOW%RHO,ZTEMP,ZWLIQ)
287           DEALLOCATE(ZTEMP)
288           DEALLOCATE(ZWLIQ)
289         END IF
290         !* put profile on fine snow grid
291         DO JLAYER=1,NGRID_LEVEL
292           ZFIELD_FINE(:,JLAYER,:) = ZFIELD(:,1,:)
293         END DO
294         ALLOCATE(PFIELD(INI,NGRID_LEVEL,IVEGTYPE))
295         CALL PUT_ON_ALL_VEGTYPES(INI,NGRID_LEVEL,IPATCH,IVEGTYPE,ZFIELD_FINE,PFIELD)
296
297       CASE ('3-L','CRO')
298         !* input physical variable
299         IF (HSURF(1:1)=='R') ZFIELD(:,:,:) = TZSNOW%RHO (:,1:TZSNOW%NLAYER,:)
300         IF (HSURF(1:1)=='H') ZFIELD(:,:,:) = TZSNOW%HEAT(:,1:TZSNOW%NLAYER,:)
301         !
302         IF (OSNOW_IDEAL) THEN
303           ALLOCATE(ZFIELD_FINE(INI,KLAYER,IPATCH))
304           IF (HSURF(1:1)=='R') ZFIELD_FINE(:,:,:) = ZFIELD(:,:,:)
305           IF (HSURF(1:1)=='H') ZFIELD_FINE(:,:,:) = ZFIELD(:,:,:)
306           ALLOCATE(PFIELD(INI,KLAYER,IVEGTYPE))
307           CALL PUT_ON_ALL_VEGTYPES(INI,KLAYER,IPATCH,IVEGTYPE,ZFIELD_FINE,PFIELD)
308         ELSE
309           ALLOCATE(ZFIELD_FINE(INI,NGRID_LEVEL,IPATCH))      
310           !* total depth
311           ALLOCATE(ZD(INI,IPATCH))
312           ZD = 0.
313           DO JLAYER=1,TZSNOW%NLAYER
314             ZD(:,:) = ZD(:,:) + TZSNOW%WSNOW(:,JLAYER,:)/TZSNOW%RHO(:,JLAYER,:)
315           END DO
316           !* input snow layer thickness
317           ALLOCATE(ZDEPTH(INI,TZSNOW%NLAYER,IPATCH))
318           DO JPATCH=1,IPATCH
319             IF (KLAYER>=3) THEN
320               CALL SNOW3LGRID(ZDEPTH(:,:,JPATCH),ZD(:,JPATCH))
321             ELSEIF (KLAYER==1) THEN
322               ZDEPTH(:,1,JPATCH) = ZD(:,JPATCH)
323             ENDIF
324           END DO
325           !* input normalized grid
326           ALLOCATE(ZGRID(INI,TZSNOW%NLAYER,IPATCH))
327           WHERE(ZD(:,:)>0.)
328             ZGRID(:,1,:) = 0.5 * ZDEPTH(:,1,:) / ZD(:,:)
329           ELSEWHERE
330             ZGRID(:,1,:) = 0.5 / FLOAT(KLAYER)
331           END WHERE
332           DO JLAYER = 2,TZSNOW%NLAYER
333             WHERE (ZD(:,:)>0.)
334               ZGRID(:,JLAYER,:) = ZGRID(:,JLAYER-1,:)                     &
335                                  + 0.5 * ZDEPTH(:,JLAYER-1,:) / ZD(:,:)   &
336                                  + 0.5 * ZDEPTH(:,JLAYER  ,:) / ZD(:,:)
337             ELSEWHERE
338               ZGRID(:,JLAYER,:) = (FLOAT(JLAYER)-0.5) / FLOAT(KLAYER)
339             END WHERE
340           END DO
341           DEALLOCATE(ZD)
342           DEALLOCATE(ZDEPTH)
343           !* interpolation of profile onto fine normalized snow grid
344           DO JPATCH=1,IPATCH
345             CALL INTERP_GRID(ZGRID(:,:,JPATCH),ZFIELD(:,:,JPATCH),    &
346                              XGRID_SNOW(:),    ZFIELD_FINE(:,:,JPATCH))
347           END DO
348           DEALLOCATE(ZGRID)
349           ALLOCATE(PFIELD(INI,NGRID_LEVEL,IVEGTYPE))
350           CALL PUT_ON_ALL_VEGTYPES(INI,NGRID_LEVEL,IPATCH,IVEGTYPE,ZFIELD_FINE,PFIELD)
351         ENDIF
352       END SELECT
353     !
354     !* put field form patch to all vegtypes
355     DEALLOCATE(ZFIELD)
356     DEALLOCATE(ZFIELD_FINE)
357 !
358 !*      9.     Crocus specific fields
359 !              -----------------------
360 !
361   CASE ('SG1','SG2','HIS','AGE')
362     ALLOCATE(ZFIELD(INI,KLAYER,IPATCH))
363     IF (HSURF(1:3)=='SG1') ZFIELD(:,:,:) = TZSNOW%GRAN1(:,1:KLAYER,:)
364     IF (HSURF(1:3)=='SG2') ZFIELD(:,:,:) = TZSNOW%GRAN2(:,1:KLAYER,:)
365     IF (HSURF(1:3)=='HIS') ZFIELD(:,:,:) = TZSNOW%HIST(:,1:KLAYER,:)
366     IF (HSURF(1:3)=='AGE') ZFIELD(:,:,:) = TZSNOW%AGE(:,1:KLAYER,:)
367     ALLOCATE(PFIELD(INI,KLAYER,IVEGTYPE))
368     CALL PUT_ON_ALL_VEGTYPES(INI,KLAYER,IPATCH,IVEGTYPE,ZFIELD,PFIELD)
369     DEALLOCATE(ZFIELD)
370   !    
371 END SELECT
372 !
373 !-------------------------------------------------------------------------------------
374 !
375 !*      9.     End of IO
376 !              ---------
377 !
378 IF (LHOOK) CALL DR_HOOK('PREP_SNOW_EXTERN',1,ZHOOK_HANDLE)
379 !
380 !-------------------------------------------------------------------------------------
381 END SUBROUTINE PREP_SNOW_EXTERN