b4564c818ab4d9b12a45f03b1b2510105df7df94
[MNH-git_open_source-lfs.git] / src / SURFEX / read_pgd_isban.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 READ_PGD_ISBA_n(HPROGRAM,OLAND_USE)
7 !     #########################################
8 !
9 !!****  *READ_PGD_ISBA_n* - routine to initialise ISBA physiographic variables 
10 !!
11 !!    PURPOSE
12 !!    -------
13 !!
14 !!**  METHOD
15 !!    ------
16 !!
17 !!    EXTERNAL
18 !!    --------
19 !!
20 !!
21 !!    IMPLICIT ARGUMENTS
22 !!    ------------------
23 !!
24 !!    REFERENCE
25 !!    ---------
26 !!
27 !!
28 !!    AUTHOR
29 !!    ------
30 !!      V. Masson   *Meteo France*      
31 !!
32 !!    MODIFICATIONS
33 !!    -------------
34 !!      Original    01/2003 
35 !!      P. Le Moigne  12/2004 : add type of photosynthesis
36 !!      B. Decharme      2008 : add XWDRAIN
37 !!      B. Decharme   06/2009 : add topographic index statistics
38 !!      A.L. Gibelin 04/2009 : dimension NBIOMASS for ISBA-A-gs
39 !!      B. Decharme  07/2012  : files of data for permafrost area and for SOC top and sub soil
40 !!      M. Moge      02/2015 READ_SURF
41 !-------------------------------------------------------------------------------
42 !
43 !*       0.    DECLARATIONS
44 !              ------------
45 !
46 USE MODD_TYPE_DATE_SURF
47 !
48 USE MODD_DATA_COVER_PAR, ONLY : JPCOVER
49 !
50 USE MODD_SURF_PAR,   ONLY : XUNDEF
51 USE MODD_SURF_ATM_n, ONLY : CNATURE, NSIZE_FULL
52 USE MODD_ISBA_n, ONLY : NPATCH, TTIME, XCOVER, XZS, CISBA, CPEDOTF,  &
53                           CPHOTO, LTR_ML, CRUNOFF, XCLAY, XSAND,     &
54                           XSOC, LSOCP, LNOF, XRM_PATCH,              &                          
55                           NGROUND_LAYER, NNBIOMASS,                  &
56                           XAOSIP, XAOSIM, XAOSJP, XAOSJM,            &
57                           XHO2IP, XHO2IM, XHO2JP, XHO2JM,            &
58                           XSSO_SLOPE, XSSO_STDEV, XRUNOFFB,          &
59                           XZ0EFFJPDIR, LCOVER, LECOCLIMAP, LCTI,     &
60                           XWDRAIN, XTI_MIN, XTI_MAX, XTI_MEAN,       &
61                           XTI_STD, XTI_SKEW, XSOILGRID, XPH, XFERT,  &
62                           LPERM, XPERM  
63 USE MODD_ISBA_GRID_n, ONLY : XLAT, XLON, XMESH_SIZE, CGRID, XGRID_PAR, NDIM
64 USE MODD_ISBA_PAR,    ONLY : XOPTIMGRID
65 USE MODD_GR_BIOG_n,   ONLY : XISOPOT, XMONOPOT
66 USE MODD_CH_ISBA_n,   ONLY : LCH_BIO_FLUX, LCH_NO_FLUX
67 !
68 USE MODI_INIT_IO_SURF_n
69 USE MODI_END_IO_SURF_n
70 !
71 USE MODI_READ_SURF
72 USE MODI_READ_GRID
73 USE MODI_READ_LCOVER
74 USE MODI_READ_PGD_ISBA_PAR_n
75 USE MODI_READ_PGD_TSZ0_PAR_n
76 !
77 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
78 USE PARKIND1  ,ONLY : JPRB
79 !
80 USE MODI_GET_TYPE_DIM_n
81 USE MODI_READ_LECOCLIMAP
82 !
83 USE MODI_ABOR1_SFX
84 !
85 USE MODI_GET_LUOUT
86 USE MODI_PACK_SAME_RANK
87 USE MODI_GET_SURF_MASK_n
88 !
89 IMPLICIT NONE
90 !
91 !*       0.1   Declarations of arguments
92 !              -------------------------
93 !
94  CHARACTER(LEN=6),  INTENT(IN)  :: HPROGRAM ! calling program
95 LOGICAL,           INTENT(IN)  :: OLAND_USE ! 
96 !
97 !*       0.2   Declarations of local variables
98 !              -------------------------------
99 !
100 INTEGER, DIMENSION(:), POINTER :: IMASK  ! mask for packing from complete field to nature field
101 !
102 REAL, DIMENSION(:,:), ALLOCATABLE :: ZWORK
103 !
104 INTEGER           :: IRESP          ! Error code after redding
105 !
106  CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
107 !
108 INTEGER :: ILU    ! expected physical size of full surface array
109 INTEGER :: ILUOUT ! output listing logical unit
110 INTEGER :: JLAYER  ! loop counter on layers
111 INTEGER :: IVERSION, IBUGFIX   ! surface version
112 REAL(KIND=JPRB) :: ZHOOK_HANDLE
113 !
114 !-------------------------------------------------------------------------------
115 !
116 !* 1D physical dimension
117 !
118 IF (LHOOK) CALL DR_HOOK('READ_PGD_ISBA_N',0,ZHOOK_HANDLE)
119 YRECFM='SIZE_NATURE'
120  CALL GET_TYPE_DIM_n('NATURE',NDIM)
121 !
122 YRECFM='VERSION'
123  CALL READ_SURF(HPROGRAM,YRECFM,IVERSION,IRESP)
124 !
125 YRECFM='BUG'
126  CALL READ_SURF(HPROGRAM,YRECFM,IBUGFIX,IRESP)
127 !
128 !*       2.     Dimension initializations:
129 !               -------------------------
130 !
131 !* soil scheme
132 !
133 YRECFM='ISBA'
134  CALL READ_SURF(HPROGRAM,YRECFM,CISBA,IRESP)
135 !
136 IF (IVERSION>=7) THEN
137   !
138   !* Pedo-transfert function
139   !
140   YRECFM='PEDOTF'
141   CALL READ_SURF(HPROGRAM,YRECFM,CPEDOTF,IRESP)
142   !
143 ELSE
144   CPEDOTF = 'CH78'
145 ENDIF
146 !
147 !* type of photosynthesis
148 !
149 YRECFM='PHOTO'
150  CALL READ_SURF(HPROGRAM,YRECFM,CPHOTO,IRESP)
151 !
152 !* new radiative transfert
153 !
154 IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=2) THEN
155   !
156   YRECFM='TR_ML'
157   CALL READ_SURF(HPROGRAM,YRECFM,LTR_ML,IRESP)
158   !
159 ELSE 
160   LTR_ML = .FALSE.
161 ENDIF
162 !
163 !* threshold to remove little fractions of patches
164 !
165 IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) THEN
166   !
167   YRECFM='RM_PATCH'
168   CALL READ_SURF(HPROGRAM,YRECFM,XRM_PATCH,IRESP)
169   !
170 ELSE 
171   XRM_PATCH = 0.0
172 ENDIF
173 !
174 !* number of soil layers
175 !
176 YRECFM='GROUND_LAYER'
177  CALL READ_SURF(HPROGRAM,YRECFM,NGROUND_LAYER,IRESP)
178 !
179 !* Reference grid for DIF
180 !
181 IF(CISBA=='DIF') THEN
182   ALLOCATE(XSOILGRID(NGROUND_LAYER))
183   XSOILGRID=XUNDEF
184   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=2) THEN
185     YRECFM='SOILGRID'
186     CALL READ_SURF(HPROGRAM,YRECFM,XSOILGRID,IRESP,HDIR='-')
187   ELSE
188     XSOILGRID(1:NGROUND_LAYER)=XOPTIMGRID(1:NGROUND_LAYER)
189   ENDIF
190 ELSE
191   ALLOCATE(XSOILGRID(0))
192 ENDIF
193 !
194 !* number of biomass pools
195 !
196 IF (IVERSION>=6) THEN
197   YRECFM='NBIOMASS'
198   CALL READ_SURF(HPROGRAM,YRECFM,NNBIOMASS,IRESP)
199 ELSE
200   SELECT CASE (CPHOTO)
201     CASE ('AGS','LAI','AST','LST')
202       NNBIOMASS = 1
203     CASE ('NIT')
204       NNBIOMASS = 3
205     CASE ('NCB')
206       NNBIOMASS = 6
207   END SELECT
208 ENDIF
209 !
210 !* number of tiles
211 !
212 YRECFM='PATCH_NUMBER'
213  CALL READ_SURF(HPROGRAM,YRECFM,NPATCH,IRESP)
214 !
215 !
216 !*       3.     Physiographic data fields:
217 !               -------------------------
218 !
219 !
220 !*       3.1    Cover classes :
221 !               -------------
222 !
223 ALLOCATE(LCOVER(JPCOVER))
224  CALL READ_LCOVER(HPROGRAM,LCOVER)
225 !
226 ALLOCATE(XCOVER(NDIM,JPCOVER))
227  CALL READ_SURF(HPROGRAM,'COVER',XCOVER(:,:),LCOVER,IRESP,HDIR='H')
228 !
229 !*       3.2    Orography :
230 !               ---------
231 !
232 !
233 ALLOCATE(XZS(NDIM))
234 YRECFM='ZS'
235  CALL READ_SURF(HPROGRAM,YRECFM,XZS(:),IRESP)
236 !
237 !
238 !* latitude, longitude, mesh size, and heading of JP axis (deg from N clockwise)
239 !
240 ALLOCATE(XLAT       (NDIM))
241 ALLOCATE(XLON       (NDIM))
242 ALLOCATE(XMESH_SIZE (NDIM))
243 ALLOCATE(XZ0EFFJPDIR(NDIM))
244  CALL READ_GRID(HPROGRAM,CGRID,XGRID_PAR,XLAT,XLON,XMESH_SIZE,IRESP,XZ0EFFJPDIR)
245 !
246 !* clay fraction : attention, seul un niveau est present dans le fichier
247 !* on rempli tout les niveaux de  XCLAY avec les valeurs du fichiers
248 !
249 ALLOCATE(XCLAY(NDIM,NGROUND_LAYER))
250 YRECFM='CLAY'
251  CALL READ_SURF(HPROGRAM,YRECFM,XCLAY(:,1),IRESP)
252 DO JLAYER=2,NGROUND_LAYER
253   XCLAY(:,JLAYER)=XCLAY(:,1)
254 END DO
255 !
256 !* sand fraction
257 !
258 ALLOCATE(XSAND(NDIM,NGROUND_LAYER))
259 YRECFM='SAND'
260  CALL READ_SURF(HPROGRAM,YRECFM,XSAND(:,1),IRESP)
261 DO JLAYER=2,NGROUND_LAYER
262   XSAND(:,JLAYER)=XSAND(:,1)
263 END DO
264 !
265 !* Soil organic carbon profile
266 !
267 IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) THEN
268    YRECFM='SOCP'
269    CALL READ_SURF(HPROGRAM,YRECFM,LSOCP,IRESP)
270 ELSE
271    LSOCP=.FALSE.
272 ENDIF
273 !
274 IF(LSOCP)THEN
275 !  
276   ALLOCATE(XSOC (NDIM,NGROUND_LAYER))
277 !
278   YRECFM='SOC_TOP'
279   CALL READ_SURF(HPROGRAM,YRECFM,XSOC(:,1),IRESP)
280   YRECFM='SOC_SUB'
281   CALL READ_SURF(HPROGRAM,YRECFM,XSOC(:,2),IRESP)
282 !
283   DO JLAYER=2,NGROUND_LAYER
284     XSOC (:,JLAYER)=XSOC (:,2)
285   END DO
286 !
287 ELSE
288 !  
289   ALLOCATE(XSOC (0,1))
290 !
291 ENDIF
292 !
293 !* permafrost distribution
294 !
295 IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) THEN
296    YRECFM='PERMAFROST'
297    CALL READ_SURF(HPROGRAM,YRECFM,LPERM,IRESP)
298 ELSE
299    LPERM=.FALSE.
300 ENDIF
301 !
302 IF(LPERM)THEN
303 !  
304   ALLOCATE(XPERM (NDIM))
305 !
306   YRECFM='PERM'
307   CALL READ_SURF(HPROGRAM,YRECFM,XPERM(:),IRESP)
308 !
309 ELSE
310 !  
311   ALLOCATE(XPERM (0))
312 !
313 ENDIF
314 !
315 IF (IVERSION>7 .OR. (IVERSION==7 .AND. IBUGFIX>=3)) THEN
316    YRECFM='NO'
317    CALL READ_SURF(HPROGRAM,YRECFM,LNOF,IRESP)
318 ELSE
319    LNOF = .FALSE.
320 ENDIF
321 !
322 !SOILNOX
323 !
324 IF (LCH_NO_FLUX) THEN
325   !
326   IF (LNOF) THEN
327     !
328     ALLOCATE(XPH(NDIM))
329     YRECFM='PH'
330     CALL READ_SURF(HPROGRAM,YRECFM,XPH(:),IRESP)
331     !
332     ALLOCATE(XFERT(NDIM))
333     YRECFM='FERT'
334     CALL READ_SURF(HPROGRAM,YRECFM,XFERT(:),IRESP)
335     !
336   ELSE
337     CALL ABOR1_SFX("READ_PGD_ISBAn: WITH LCH_NO_FLUX=T, PH AND FERT FIELDS ARE GIVEN AT PGD STEP")
338   ENDIF
339   !
340 ELSE
341   ALLOCATE(XPH (0))
342   ALLOCATE(XFERT(0))
343 END IF
344 !
345 !* subgrid-scale orography parameters to compute dynamical roughness length
346 !
347 ALLOCATE(XAOSIP(NDIM))
348 YRECFM='AOSIP'
349  CALL READ_SURF(HPROGRAM,YRECFM,XAOSIP,IRESP)
350 !
351 ALLOCATE(XAOSIM(NDIM))
352 YRECFM='AOSIM'
353  CALL READ_SURF(HPROGRAM,YRECFM,XAOSIM,IRESP)
354
355 ALLOCATE(XAOSJP(NDIM))
356 YRECFM='AOSJP'
357  CALL READ_SURF(HPROGRAM,YRECFM,XAOSJP,IRESP)
358 !
359 ALLOCATE(XAOSJM(NDIM))
360 YRECFM='AOSJM'
361  CALL READ_SURF(HPROGRAM,YRECFM,XAOSJM,IRESP)
362 !
363 ALLOCATE(XHO2IP(NDIM))
364 YRECFM='HO2IP'
365  CALL READ_SURF(HPROGRAM,YRECFM,XHO2IP,IRESP)
366 !
367 ALLOCATE(XHO2IM(NDIM))
368 YRECFM='HO2IM'
369  CALL READ_SURF(HPROGRAM,YRECFM,XHO2IM,IRESP)
370 !
371 ALLOCATE(XHO2JP(NDIM))
372 YRECFM='HO2JP'
373  CALL READ_SURF(HPROGRAM,YRECFM,XHO2JP,IRESP)
374 !
375 ALLOCATE(XHO2JM(NDIM))
376 YRECFM='HO2JM'
377  CALL READ_SURF(HPROGRAM,YRECFM,XHO2JM,IRESP)
378 !
379 !* orographic parameter to compute effective surface of energy exchanges
380 !
381 ALLOCATE(XSSO_SLOPE(NDIM))
382 YRECFM='SSO_SLOPE'
383  CALL READ_SURF(HPROGRAM,YRECFM,XSSO_SLOPE,IRESP)
384 !
385 !* orographic standard deviation for subgrid-scale orographic drag
386 !
387 ALLOCATE(XSSO_STDEV(NDIM))
388 YRECFM='SSO_STDEV'
389  CALL READ_SURF(HPROGRAM,YRECFM,XSSO_STDEV(:),IRESP)
390 !
391 !* orographic runoff coefficient
392 !
393 ALLOCATE(XRUNOFFB(NDIM))
394 YRECFM='RUNOFFB'
395  CALL READ_SURF(HPROGRAM,YRECFM,XRUNOFFB,IRESP)
396 !
397 !* subgrid drainage coefficient
398 !
399 ALLOCATE(XWDRAIN(NDIM))
400 IF (IVERSION<=3) THEN
401   XWDRAIN = 0.
402 ELSE
403   YRECFM='WDRAIN'
404   CALL READ_SURF(HPROGRAM,YRECFM,XWDRAIN,IRESP)
405 ENDIF
406 !
407 !* topographic index statistics
408 !
409 IF(CRUNOFF=='SGH ' .AND. IVERSION>=5) THEN 
410 !
411   YRECFM='CTI'
412   CALL READ_SURF(HPROGRAM,YRECFM,LCTI,IRESP)        
413 !
414   IF (.NOT.LCTI) CALL ABOR1_SFX("READ_PGD_ISBA_n:WITH CRUNOFF=SGH, CTI MAPS MUST BE GIVEN TO PGD")
415   !
416   ALLOCATE(XTI_MIN(NDIM))
417   ALLOCATE(XTI_MAX(NDIM))
418   ALLOCATE(XTI_MEAN(NDIM))
419   ALLOCATE(XTI_STD(NDIM))
420   ALLOCATE(XTI_SKEW(NDIM))
421 !
422   YRECFM='TI_MIN'
423   CALL READ_SURF(HPROGRAM,YRECFM,XTI_MIN,IRESP)
424 !
425   YRECFM='TI_MAX'
426   CALL READ_SURF(HPROGRAM,YRECFM,XTI_MAX,IRESP)
427 !
428   YRECFM='TI_MEAN'
429   CALL READ_SURF(HPROGRAM,YRECFM,XTI_MEAN,IRESP)
430 !
431   YRECFM='TI_STD'
432   CALL READ_SURF(HPROGRAM,YRECFM,XTI_STD,IRESP)
433 !
434   YRECFM='TI_SKEW'
435   CALL READ_SURF(HPROGRAM,YRECFM,XTI_SKEW,IRESP)
436 !
437 ELSE
438 !
439   ALLOCATE(XTI_MIN(0))
440   ALLOCATE(XTI_MAX(0))
441   ALLOCATE(XTI_MEAN(0))
442   ALLOCATE(XTI_STD(0))
443   ALLOCATE(XTI_SKEW(0))
444 !
445 ENDIF
446 !
447 !-------------------------------------------------------------------------------
448 !
449 !* biogenic chemical emissions
450 !
451 IF (LCH_BIO_FLUX) THEN
452   ALLOCATE(ZWORK(NSIZE_FULL,1))
453   !
454   CALL END_IO_SURF_n(HPROGRAM)
455   CALL INIT_IO_SURF_n(HPROGRAM,'FULL  ','SURF  ','READ ')
456   !
457   CALL GET_LUOUT(HPROGRAM,ILUOUT)
458   ALLOCATE(IMASK(NDIM))
459   ILU=0
460   CALL GET_SURF_MASK_n('NATURE',NDIM,IMASK,ILU,ILUOUT)
461   ALLOCATE(XISOPOT(NDIM))
462   ALLOCATE(XMONOPOT(NDIM))
463   !
464   ZWORK(:,:) = 0.  
465   YRECFM='E_ISOPOT'
466   CALL READ_SURF(HPROGRAM,YRECFM,ZWORK,IRESP)
467   CALL PACK_SAME_RANK(IMASK,ZWORK(:,1),XISOPOT(:))
468   !
469   ZWORK(:,:) = 0.  
470   YRECFM='E_MONOPOT'
471   CALL READ_SURF(HPROGRAM,YRECFM,ZWORK,IRESP)
472   CALL PACK_SAME_RANK(IMASK,ZWORK(:,1),XMONOPOT(:))
473   !
474   CALL END_IO_SURF_n(HPROGRAM)
475   CALL INIT_IO_SURF_n(HPROGRAM,'NATURE','ISBA  ','READ ')
476   !
477   DEALLOCATE(ZWORK)
478 ELSE
479   ALLOCATE(XISOPOT (0))
480   ALLOCATE(XMONOPOT(0))
481 END IF
482 !
483 !-------------------------------------------------------------------------------
484 !
485 !*       4.     Physiographic data fields not to be computed by ecoclimap
486 !               ---------------------------------------------------------
487 !
488  CALL READ_LECOCLIMAP(HPROGRAM,LECOCLIMAP)
489 !
490  CALL READ_PGD_ISBA_PAR_n(HPROGRAM,NDIM,OLAND_USE)
491 IF (CNATURE == 'TSZ0') CALL READ_PGD_TSZ0_PAR_n(HPROGRAM)
492 !
493 IF (LHOOK) CALL DR_HOOK('READ_PGD_ISBA_N',1,ZHOOK_HANDLE)
494 !-------------------------------------------------------------------------------
495 !
496 END SUBROUTINE READ_PGD_ISBA_n