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