Philippe 07/03/2019: IO bugfix: io_set_mnhversion must be called by all the processes
[MNH-git_open_source-lfs.git] / src / SURFEX / read_isban.F90
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 READ_ISBA_n (DTCO, IO, S, NP, NPE, PCLAY, U, HPROGRAM)
7 !     ##################################
8 !
9 !!****  *READ_ISBA_n* - routine to initialise ISBA variables
10 !!                         
11 !!
12 !!    PURPOSE
13 !!    -------
14 !!
15 !!**  METHOD
16 !!    ------
17 !!
18 !!    EXTERNAL
19 !!    --------
20 !!
21 !!
22 !!    IMPLICIT ARGUMENTS
23 !!    ------------------
24 !!
25 !!    REFERENCE
26 !!    ---------
27 !!
28 !!
29 !!    AUTHOR
30 !!    ------
31 !!      V. Masson   *Meteo France*
32 !!
33 !!    MODIFICATIONS
34 !!    -------------
35 !!      Original    01/2003
36 !!
37 !!      READ_SURF for general reading : 08/2003 (S.Malardel)
38 !!      B. Decharme  2008    : Floodplains
39 !!      B. Decharme  01/2009 : Optional Arpege deep soil temperature read
40 !!      A.L. Gibelin   03/09 : modifications for CENTURY model 
41 !!      A.L. Gibelin    04/2009 : BIOMASS and RESP_BIOMASS arrays 
42 !!      A.L. Gibelin    06/2009 : Soil carbon variables for CNT option
43 !!      B. Decharme  09/2012 : suppress NWG_LAYER (parallelization problems)
44 !!      T. Aspelien  08/2013 : Read diagnostics for assimilation
45 !!      P. Samuelsson   10/2014 : MEB
46 !!
47 !-------------------------------------------------------------------------------
48 !
49 !*       0.    DECLARATIONS
50 !              ------------
51 !
52 USE MODD_DATA_COVER_n, ONLY : DATA_COVER_t
53 USE MODD_ISBA_OPTIONS_n, ONLY : ISBA_OPTIONS_t
54 USE MODD_ISBA_n, ONLY : ISBA_S_t, ISBA_NP_t, ISBA_NPE_t, ISBA_P_t, ISBA_PE_t
55 USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
56 !
57 USE MODD_CO2V_PAR,       ONLY : XANFMINIT, XCONDCTMIN
58 !                          
59 USE MODD_ASSIM,          ONLY : LASSIM,CASSIM_ISBA,XAT2M_ISBA,XAHU2M_ISBA,&
60                                 XAZON10M_ISBA,XAMER10M_ISBA,NIFIC,NVAR, &
61                                 COBS,NOBSTYPE,CVAR,LPRT,XTPRT,NIVAR,CBIO, &
62                                 XADDINFL,NENS,XSIGMA,NIE
63 !                                
64 USE MODD_SURF_PAR,       ONLY : XUNDEF, NUNDEF
65 USE MODD_SNOW_PAR,       ONLY : XZ0SN
66 USE MODD_ISBA_PAR,       ONLY : XWGMIN
67 !
68 USE MODE_READ_SURF_LAYERS
69 !
70 USE MODI_READ_SURF
71 USE MODI_MAKE_CHOICE_ARRAY
72 USE MODI_PACK_SAME_RANK
73 !
74 USE MODI_READ_GR_SNOW
75 USE MODI_ABOR1_SFX
76 USE MODI_IO_BUFF
77 !
78 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
79 USE PARKIND1  ,ONLY : JPRB
80 !
81 USE MODI_IO_BUFF_CLEAN
82 USE MODI_GET_TYPE_DIM_n
83 USE MODE_RANDOM
84 USE MODE_EKF
85 !
86 IMPLICIT NONE
87 !
88 !*       0.1   Declarations of arguments
89 !              -------------------------
90 !
91 !
92 TYPE(DATA_COVER_t), INTENT(INOUT) :: DTCO
93 TYPE(ISBA_OPTIONS_t), INTENT(INOUT) :: IO
94 TYPE(ISBA_S_t), INTENT(INOUT) :: S
95 TYPE(ISBA_NP_t), INTENT(INOUT) :: NP
96 TYPE(ISBA_NPE_t), INTENT(INOUT) :: NPE
97 REAL, DIMENSION(:,:), INTENT(IN) :: PCLAY
98 TYPE(SURF_ATM_t), INTENT(INOUT) :: U
99 !
100  CHARACTER(LEN=6),  INTENT(IN)  :: HPROGRAM ! calling program
101 !
102 !*       0.2   Declarations of local variables
103 !              -------------------------------
104 !
105 TYPE(ISBA_P_t), POINTER :: PK
106 TYPE(ISBA_PE_t), POINTER :: PEK
107 !
108 INTEGER           :: ILU          ! 1D physical dimension
109 !
110 INTEGER           :: IRESP          ! Error code after redding
111 !
112 CHARACTER(LEN=LEN_HREC) :: YRECFM         ! Name of the article to be read
113 CHARACTER(LEN=LEN_HREC) :: YCBIO          ! Name of biomass variable
114 !
115 CHARACTER(LEN=4)  :: YLVL
116 !
117 REAL, DIMENSION(:,:,:),ALLOCATABLE :: ZLAI
118 REAL, DIMENSION(:,:,:),ALLOCATABLE  :: ZWORK3D    ! 2D array to write data in file
119 REAL, DIMENSION(:,:),ALLOCATABLE  :: ZWORK      ! 2D array to write data in file
120 REAL, DIMENSION(:), ALLOCATABLE :: ZCOFSWI
121 !
122 REAL,DIMENSION(IO%NPATCH) :: ZVLAIMIN
123 REAL :: ZCOEF
124 !
125 INTEGER :: IWORK   ! Work integer
126 !
127 INTEGER :: JP, JL, JNB, JNLITTER, JNSOILCARB, JNLITTLEVS  ! loop counter on layers
128 INTEGER :: JVAR, JI
129 !
130 INTEGER           :: IVERSION       ! surface version
131 INTEGER           :: IBUGFIX
132 INTEGER           :: IIVAR
133 INTEGER           :: IOBS
134 INTEGER           :: IBSUP
135 INTEGER           :: ISIZE_LMEB_PATCH, IMASK
136 !
137 LOGICAL :: GKNOWN, GDIM
138 !
139 REAL(KIND=JPRB) :: ZHOOK_HANDLE
140 !
141 !-------------------------------------------------------------------------------
142 !
143 !
144 !* 1D physical dimension
145 !
146 IF (LHOOK) CALL DR_HOOK('READ_ISBA_N',0,ZHOOK_HANDLE)
147 YRECFM='SIZE_NATURE'
148  CALL GET_TYPE_DIM_n(DTCO, U, 'NATURE',ILU)
149 !
150 YRECFM='VERSION'
151  CALL READ_SURF(HPROGRAM,YRECFM,IVERSION,IRESP)
152 !
153 YRECFM='BUG'
154  CALL READ_SURF(HPROGRAM,YRECFM,IBUGFIX,IRESP)
155 !
156 GDIM = (IVERSION>8 .OR. IVERSION==8 .AND. IBUGFIX>0)
157 IF (GDIM) CALL READ_SURF(HPROGRAM,'SPLIT_PATCH',GDIM,IRESP)
158 !
159 YRECFM='BUG'
160  CALL READ_SURF(HPROGRAM,YRECFM,IBUGFIX,IRESP)
161
162 !*       2.     Prognostic fields:
163 !               -----------------
164 !
165 ALLOCATE(ZWORK(ILU,IO%NPATCH))
166 !* soil temperatures
167 !
168 IF(IO%LTEMP_ARP)THEN
169   IWORK=IO%NTEMPLAYER_ARP
170 ELSEIF(IO%CISBA=='DIF')THEN
171   IWORK=IO%NGROUND_LAYER
172 ELSE
173   IWORK=2 !Only 2 temperature layer in ISBA-FR
174 ENDIF
175 !
176 IF ( TRIM(CASSIM_ISBA)=="ENKF") THEN
177   DO JP = 1,IO%NPATCH
178     PK => NP%AL(JP)
179     ALLOCATE(PK%XRED_NOISE(PK%NSIZE_P,NVAR))
180     PK%XRED_NOISE(:,:) = 0.
181   ENDDO
182 ELSE
183   DO JP = 1,IO%NPATCH
184     ALLOCATE(NP%AL(JP)%XRED_NOISE(0,0))
185   ENDDO  
186 ENDIF
187 !
188 IF  ( TRIM(CASSIM_ISBA)=="ENKF" .OR. (TRIM(CASSIM_ISBA)=="EKF" .AND. LPRT) ) THEN
189   ALLOCATE(ZCOFSWI(ILU))
190   CALL COFSWI(PCLAY(:,1),ZCOFSWI)
191 ELSE
192   ALLOCATE(ZCOFSWI(0))
193 ENDIF
194 !
195 DO JP = 1,IO%NPATCH
196   PEK => NPE%AL(JP)
197   PK => NP%AL(JP)
198   !
199   ALLOCATE(PEK%XTG(PK%NSIZE_P,IWORK))
200   !
201   ALLOCATE(PEK%XWG (PK%NSIZE_P,IO%NGROUND_LAYER))
202   ALLOCATE(PEK%XWGI(PK%NSIZE_P,IO%NGROUND_LAYER))
203   !
204   PEK%XTG (:,:) = XUNDEF
205   PEK%XWG (:,:) = XUNDEF
206   PEK%XWGI(:,:) = XUNDEF
207   !
208   ALLOCATE(PEK%XWR(PK%NSIZE_P))
209   !
210 ENDDO
211 !
212 ALLOCATE(ZWORK3D(ILU,IWORK,IO%NPATCH))
213 CALL READ_SURF_LAYERS(HPROGRAM,'TG',GDIM,ZWORK3D,IRESP)
214 DO JL=1,IWORK
215   DO JP = 1,IO%NPATCH
216     CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK3D(:,JL,JP),NPE%AL(JP)%XTG(:,JL))
217   ENDDO
218 ENDDO
219 DEALLOCATE(ZWORK3D)
220 !
221 ! Perturb value if requested
222 IF ( TRIM(CASSIM_ISBA)=="EKF" .AND. LPRT ) THEN
223   !
224   DO JL=1,IWORK
225   ! read in control variable
226     IF ( (TRIM(CVAR(NIVAR))=="TG1" .AND. JL==1) .OR. (TRIM(CVAR(NIVAR))=="TG2" .AND. JL==2) ) THEN
227       DO JP = 1,IO%NPATCH
228         PEK => NPE%AL(JP)
229         WHERE ( PEK%XTG(:,JL)/=XUNDEF )
230           PEK%XTG(:,JL) = PEK%XTG(:,JL) + XTPRT(NIVAR)*PEK%XTG(:,JL)
231         ENDWHERE
232       ENDDO
233     ENDIF
234   END DO
235   !
236 ELSEIF ( TRIM(CASSIM_ISBA)=="ENKF" .AND. NIE<NENS+1 ) THEN
237   !
238   CALL MAKE_ENS_ENKF(IWORK,ILU,"TG ",ZCOFSWI,NP)
239   !
240 ENDIF
241 !
242 !
243 !* soil liquid and ice water contents
244 !
245 ALLOCATE(ZWORK3D(ILU,IO%NGROUND_LAYER,IO%NPATCH))
246 CALL READ_SURF_LAYERS(HPROGRAM,'WG',GDIM,ZWORK3D,IRESP)
247 DO JL=1,IO%NGROUND_LAYER
248   DO JP = 1,IO%NPATCH
249     CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK3D(:,JL,JP),NPE%AL(JP)%XWG(:,JL))
250   ENDDO
251 ENDDO
252 DEALLOCATE(ZWORK3D)
253 !
254 ! Perturb value if requested
255 IF ( TRIM(CASSIM_ISBA)=="EKF" .AND. LPRT ) THEN
256    !
257    DO JL=1,IO%NGROUND_LAYER
258      ! read in control variable
259      IF ( (TRIM(CVAR(NIVAR))=="WG1" .AND. JL==1) .OR. & 
260           (TRIM(CVAR(NIVAR))=="WG2" .AND. JL==2) .OR. &
261           (TRIM(CVAR(NIVAR))=="WG3" .AND. JL==3) .OR. &
262           (TRIM(CVAR(NIVAR))=="WG4" .AND. JL==4) .OR. &
263           (TRIM(CVAR(NIVAR))=="WG5" .AND. JL==5) .OR. &
264           (TRIM(CVAR(NIVAR))=="WG6" .AND. JL==6) .OR. &
265           (TRIM(CVAR(NIVAR))=="WG7" .AND. JL==7) .OR. &
266           (TRIM(CVAR(NIVAR))=="WG8" .AND. JL==8) ) THEN
267             
268        DO JP = 1,IO%NPATCH
269          PEK => NPE%AL(JP)
270          PK => NP%AL(JP)
271          DO JI = 1,PK%NSIZE_P
272            IMASK = PK%NR_P(JI)
273            IF (PEK%XWG(JI,JL)/=XUNDEF ) THEN
274              PEK%XWG(JI,JL) = PEK%XWG(JI,JL) + XTPRT(NIVAR) * ZCOFSWI(IMASK) 
275            ENDIF
276          ENDDO
277        END DO
278        !
279      ENDIF
280      !
281    END DO
282    !
283 ELSEIF ( TRIM(CASSIM_ISBA)=="ENKF" .AND. NIE<NENS+1 ) THEN
284   !
285   CALL MAKE_ENS_ENKF(IWORK,ILU,"WG ",ZCOFSWI,NP)
286   !
287 ENDIF
288 !
289 IF(IO%CISBA=='DIF')THEN
290   IWORK=IO%NGROUND_LAYER
291 ELSE
292   IWORK=2 !Only 2 soil ice layer in ISBA-FR
293 ENDIF
294 !
295 ALLOCATE(ZWORK3D(ILU,IWORK,IO%NPATCH))
296 CALL READ_SURF_LAYERS(HPROGRAM,'WGI',GDIM,ZWORK3D,IRESP)
297 DO JL=1,IWORK
298   DO JP = 1,IO%NPATCH
299     CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK3D(:,JL,JP),NPE%AL(JP)%XWGI(:,JL))
300   ENDDO
301 ENDDO
302 DEALLOCATE(ZWORK3D)
303 !
304 !* water intercepted on leaves
305 !
306 !
307 YRECFM = 'WR'
308  CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
309 DO JP = 1,IO%NPATCH
310   CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XWR(:))
311 ENDDO
312 !
313 !* Leaf Area Index
314 !
315 IF (IO%CPHOTO=='NIT' .OR. IO%CPHOTO=='NCB') THEN
316
317   YRECFM = 'LAI'
318   CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)  
319   DO JP = 1,IO%NPATCH
320     CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XLAI(:))
321   ENDDO
322
323   IF ( TRIM(CASSIM_ISBA)=="EKF" .AND. LPRT ) THEN
324     !
325     ! read in control variable
326     IF ( TRIM(CVAR(NIVAR))=="LAI" ) THEN
327       DO JP = 1,IO%NPATCH
328         PEK => NPE%AL(JP)
329         WHERE ( PEK%XLAI(:)/=XUNDEF ) 
330            PEK%XLAI(:) =  PEK%XLAI(:) + XTPRT(NIVAR)* PEK%XLAI(:)
331         ENDWHERE
332       ENDDO
333     ENDIF
334     !
335   ELSEIF ( TRIM(CASSIM_ISBA)=="ENKF" .AND. NIE<NENS+1 ) THEN
336     !
337     IF (IO%NPATCH==12) THEN
338       ZVLAIMIN = (/0.3,0.3,0.3,0.3,1.0,1.0,0.3,0.3,0.3,0.3,0.3,0.3/)
339     ELSE
340       ZVLAIMIN = (/0.3/)
341     ENDIF
342     !
343     ALLOCATE(ZLAI(ILU,1,IO%NPATCH))
344     ZLAI(:,:,:) = 0.
345     DO JP = 1,IO%NPATCH
346       ZLAI(1:NP%AL(JP)%NSIZE_P,1,JP) = NPE%AL(JP)%XLAI(:)
347     ENDDO
348     CALL MAKE_ENS_ENKF(1,ILU,"LAI",ZCOFSWI,NP,ZLAI)
349     DO JP = 1,IO%NPATCH
350       NPE%AL(JP)%XLAI(:) = MAX(ZVLAIMIN(JP),ZLAI(1:NP%AL(JP)%NSIZE_P,1,JP))
351     ENDDO
352     DEALLOCATE(ZLAI)
353     !    
354   ENDIF  
355 END IF
356 !
357 !* snow mantel
358 !
359 DO JP = 1,IO%NPATCH
360   IF (JP>1) THEN
361     NPE%AL(JP)%TSNOW%SCHEME = NPE%AL(1)%TSNOW%SCHEME
362     NPE%AL(JP)%TSNOW%NLAYER = NPE%AL(1)%TSNOW%NLAYER
363   ENDIF
364   CALL READ_GR_SNOW(HPROGRAM,'VEG','     ',ILU,NP%AL(JP)%NSIZE_P,NP%AL(JP)%NR_P,JP,&
365                     NPE%AL(JP)%TSNOW,KNPATCH=IO%NPATCH  )
366   CALL IO_BUFF_CLEAN
367 ENDDO
368 !
369 IF(IO%LGLACIER)THEN
370   DO JP = 1,IO%NPATCH
371     ALLOCATE(NPE%AL(JP)%XICE_STO(NP%AL(JP)%NSIZE_P))
372   ENDDO
373   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=2) THEN
374     YRECFM = 'ICE_STO'
375     CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
376     DO JP = 1,IO%NPATCH
377       CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XICE_STO(:))
378     ENDDO
379   ELSE
380     DO JP = 1,IO%NPATCH
381       NPE%AL(JP)%XICE_STO(:) = 0.0
382     ENDDO
383   ENDIF
384 ELSE  
385   DO JP = 1,IO%NPATCH
386     ALLOCATE(NPE%AL(JP)%XICE_STO(0))
387   ENDDO
388 ENDIF
389 !
390 !-------------------------------------------------------------------------------
391 !
392 !*       3.  MEB Prognostic or Semi-prognostic variables
393 !            -------------------------------------------
394 !
395 ISIZE_LMEB_PATCH=COUNT(IO%LMEB_PATCH(:))
396 !
397 IF (ISIZE_LMEB_PATCH>0) THEN
398 !
399   DO JP = 1,IO%NPATCH
400     PK => NP%AL(JP)
401     PEK => NPE%AL(JP)
402     ALLOCATE(PEK%XWRL (PK%NSIZE_P))
403     ALLOCATE(PEK%XWRLI(PK%NSIZE_P))
404     ALLOCATE(PEK%XWRVN(PK%NSIZE_P))
405     ALLOCATE(PEK%XTV  (PK%NSIZE_P))
406     ALLOCATE(PEK%XTL  (PK%NSIZE_P))
407     ALLOCATE(PEK%XTC  (PK%NSIZE_P))
408     ALLOCATE(PEK%XQC  (PK%NSIZE_P))
409   ENDDO
410
411 !* water intercepted on litter
412  
413  YRECFM = 'WRL'
414  CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
415 DO JP = 1,IO%NPATCH
416   CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XWRL(:))
417 ENDDO
418
419  YRECFM = 'WRLI'
420  CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
421 DO JP = 1,IO%NPATCH
422   CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XWRLI(:))
423 ENDDO 
424 !
425 !* snow intercepted on vegetation canopy leaves
426 !
427   YRECFM = 'WRVN'
428  CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
429 DO JP = 1,IO%NPATCH
430   CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XWRVN(:))
431 ENDDO   
432 !
433 !* vegetation canopy temperature
434 !
435   YRECFM = 'TV'
436  CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
437 DO JP = 1,IO%NPATCH
438   CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XTV(:))
439 ENDDO    
440 !
441 !* litter temperature
442 !
443   YRECFM = 'TL'
444  CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
445 DO JP = 1,IO%NPATCH
446   CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XTL(:))
447 ENDDO  
448 !
449 !* vegetation canopy air temperature
450 !
451   YRECFM = 'TC'
452  CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
453 DO JP = 1,IO%NPATCH
454   CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XTC(:))
455 ENDDO    
456 !
457 !* vegetation canopy air specific humidity
458 !
459   YRECFM = 'QC'
460  CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
461 DO JP = 1,IO%NPATCH
462   CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XQC(:))
463 ENDDO     
464 !
465 ENDIF
466 !
467 !-------------------------------------------------------------------------------
468 !
469 !*       4.  Semi-prognostic variables
470 !            -------------------------
471 !
472 DO JP = 1,IO%NPATCH
473   PK => NP%AL(JP)
474   PEK => NPE%AL(JP)
475
476   ALLOCATE(PEK%XRESA(PK%NSIZE_P))
477   ALLOCATE(PEK%XLE  (PK%NSIZE_P))  
478   PEK%XRESA(:) = 100.
479   PEK%XLE(:) = 0.
480   
481   IF (IO%CPHOTO/='NON') THEN
482     ALLOCATE(PEK%XANFM  (PK%NSIZE_P))
483     ALLOCATE(PEK%XAN    (PK%NSIZE_P))
484     ALLOCATE(PEK%XANDAY (PK%NSIZE_P))
485     PEK%XANFM (:) = XANFMINIT  
486     PEK%XAN   (:) = 0.
487     PEK%XANDAY(:) = 0.      
488     !
489     ALLOCATE(PEK%XBIOMASS         (PK%NSIZE_P,IO%NNBIOMASS))
490     ALLOCATE(PEK%XRESP_BIOMASS    (PK%NSIZE_P,IO%NNBIOMASS))
491     PEK%XBIOMASS(:,:) = 0.
492     PEK%XRESP_BIOMASS(:,:) = 0.    
493   END IF
494
495 ENDDO
496 !
497 !
498 !* aerodynamical resistance
499 !
500 YRECFM = 'RESA'
501  CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
502 DO JP = 1,IO%NPATCH
503   CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XRESA(:))
504 ENDDO
505 !
506 !* patch averaged radiative temperature (K)
507 !
508 ALLOCATE(S%XTSRAD_NAT(ILU))
509 IF (IVERSION<6) THEN
510   S%XTSRAD_NAT(:)=0.
511   DO JP=1,IO%NPATCH
512     DO JI = 1,NP%AL(JP)%NSIZE_P
513       IMASK = NP%AL(JP)%NR_P(JI)
514       S%XTSRAD_NAT(IMASK) = S%XTSRAD_NAT(IMASK)+NPE%AL(JP)%XTG(JI,1)
515     ENDDO
516   ENDDO
517   S%XTSRAD_NAT(:)=S%XTSRAD_NAT(:)/IO%NPATCH
518 ELSE
519   YRECFM='TSRAD_NAT'
520   CALL READ_SURF(HPROGRAM,YRECFM,S%XTSRAD_NAT(:),IRESP)
521 ENDIF
522 !
523 DO JP = 1,IO%NPATCH
524   NPE%AL(JP)%XLE(:) = XUNDEF
525 ENDDO
526 !
527 !*       5. ISBA-AGS variables
528 !
529 IF (IO%CPHOTO/='NON') THEN
530   YRECFM = 'AN'
531   CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
532   DO JP = 1,IO%NPATCH
533     CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XAN(:))
534   ENDDO
535   !
536   YRECFM = 'ANDAY'
537   CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
538   DO JP = 1,IO%NPATCH
539     CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XANDAY(:))
540   ENDDO  
541   !
542   YRECFM = 'ANFM'
543   CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
544   DO JP = 1,IO%NPATCH
545     CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XANFM(:))
546   ENDDO  
547   !
548   YRECFM = 'LE_AGS'
549   CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
550   DO JP = 1,IO%NPATCH
551     CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XLE(:))
552   ENDDO  
553 END IF
554 !
555 IF (IO%CPHOTO=='NIT'.OR.IO%CPHOTO=='NCB') THEN
556   !
557   ALLOCATE(ZWORK3D(ILU,IO%NNBIOMASS,IO%NPATCH))    
558   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) THEN
559     YRECFM='BIOMA'
560   ELSE
561     YRECFM='BIOMASS'
562   ENDIF
563   CALL READ_SURF_LAYERS(HPROGRAM,YRECFM,GDIM,ZWORK3D,IRESP)
564   DO JNB=1,IO%NNBIOMASS
565     WRITE(YLVL,'(I1)') JNB
566     IF ( TRIM(CASSIM_ISBA)=="EKF" .AND. LPRT ) THEN
567       YCBIO = YRECFM(:LEN_TRIM(YRECFM))//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
568       ! read in control variable
569       IF ( TRIM(CVAR(NIVAR)) == "LAI" .AND. TRIM(CBIO)==TRIM(YCBIO) ) THEN
570         WHERE ( ZWORK3D(:,JNB,:)/=XUNDEF ) 
571           ZWORK3D(:,JNB,:) = ZWORK3D(:,JNB,:) * ( 1. + XTPRT(NIVAR) )
572         ENDWHERE
573       ENDIF
574     ELSEIF ( TRIM(CASSIM_ISBA)=="ENKF" .AND. NIE<NENS+1 .AND. .NOT.LASSIM ) THEN
575       !
576       IF ( TRIM(CBIO)==TRIM(YRECFM) ) THEN
577         DO JVAR = 1,NVAR
578           IF (TRIM(CVAR(JVAR)) == "LAI") THEN
579             DO JI = 1,ILU
580               DO JP = 1,IO%NPATCH
581                 ZWORK3D(JI,JNB,JP) = ZWORK3D(JI,JNB,JP) + XADDINFL(JVAR)*RANDOM_NORMAL()
582               ENDDO
583             ENDDO
584             EXIT
585           ENDIF
586         ENDDO
587       ENDIF
588       !      
589     ENDIF     
590     DO JP = 1,IO%NPATCH
591       CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK3D(:,JNB,JP),NPE%AL(JP)%XBIOMASS(:,JNB))
592     ENDDO 
593   END DO
594   DEALLOCATE(ZWORK3D)
595 !
596   IWORK=0
597   IF(IO%CPHOTO=='NCB'.OR.IVERSION<8)IWORK=2
598 !
599   DO JNB=2,IO%NNBIOMASS-IWORK
600     WRITE(YLVL,'(I1)') JNB
601     IF (IVERSION>7 .OR. (IVERSION==7 .AND. IBUGFIX>=3)) THEN
602       YRECFM='RESPI'//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
603     ELSE
604       YRECFM='RESP_BIOM'//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
605     ENDIF    
606     CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
607     IF ( TRIM(CASSIM_ISBA)=="EKF" .AND. LPRT ) THEN
608       ! read in control variable
609       IF ( TRIM(CVAR(NIVAR)) == "LAI" .AND. TRIM(CBIO)==TRIM(YRECFM) ) THEN
610         WHERE ( ZWORK(:,:)/=XUNDEF ) 
611           ZWORK(:,:) = ZWORK(:,:) + XTPRT(NIVAR)*ZWORK(:,:)
612         ENDWHERE
613     ELSEIF ( TRIM(CASSIM_ISBA)=="ENKF" .AND. NIE<NENS+1 .AND. .NOT.LASSIM ) THEN
614       !
615       IF ( TRIM(CBIO)==TRIM(YRECFM) ) THEN
616         DO JVAR = 1,NVAR
617           IF (TRIM(CVAR(JVAR)) == "LAI") THEN
618             DO JI = 1,ILU
619               DO JP = 1,IO%NPATCH
620                 ZWORK(JI,JP) = ZWORK(JI,JP) + XADDINFL(JVAR)*RANDOM_NORMAL()
621               ENDDO
622             ENDDO
623             EXIT
624           ENDIF
625         ENDDO
626       ENDIF
627       !  
628       ENDIF
629     ENDIF  
630     DO JP = 1,IO%NPATCH
631       CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XRESP_BIOMASS(:,JNB))
632     ENDDO
633   END DO
634   !
635 ENDIF
636 !
637 DEALLOCATE(ZCOFSWI)
638 !
639 !*       6. Soil carbon
640 !
641 !
642 IF (IO%CRESPSL=='CNT') THEN
643   !
644   DO JP = 1,IO%NPATCH
645     PK => NP%AL(JP)
646     PEK => NPE%AL(JP)
647     ALLOCATE(PEK%XLITTER      (PK%NSIZE_P,IO%NNLITTER,IO%NNLITTLEVS))
648     ALLOCATE(PEK%XSOILCARB    (PK%NSIZE_P,IO%NNSOILCARB))
649     ALLOCATE(PEK%XLIGNIN_STRUC(PK%NSIZE_P,IO%NNLITTLEVS))
650     !
651     PEK%XLITTER(:,:,:) = 0.
652     PEK%XSOILCARB(:,:) = 0. 
653     PEK%XLIGNIN_STRUC(:,:) = 0.  
654     !
655   ENDDO
656   !
657   DO JNLITTER=1,IO%NNLITTER
658     DO JNLITTLEVS=1,IO%NNLITTLEVS
659       WRITE(YLVL,'(I1,A1,I1)') JNLITTER,'_',JNLITTLEVS
660       YRECFM='LITTER'//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
661       CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
662       DO JP = 1,IO%NPATCH
663         CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XLITTER(:,JNLITTER,JNLITTLEVS))
664       ENDDO       
665     END DO
666   END DO
667
668   DO JNSOILCARB=1,IO%NNSOILCARB
669     WRITE(YLVL,'(I4)') JNSOILCARB
670     YRECFM='SOILCARB'//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
671     CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
672     DO JP = 1,IO%NPATCH
673       CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XSOILCARB(:,JNSOILCARB))
674     ENDDO      
675   END DO
676 !
677   DO JNLITTLEVS=1,IO%NNLITTLEVS
678     WRITE(YLVL,'(I4)') JNLITTLEVS
679     YRECFM='LIGN_STR'//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
680     CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
681     DO JP = 1,IO%NPATCH
682       CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XSOILCARB(:,JNLITTLEVS))
683     ENDDO     
684   END DO
685 !
686 ENDIF
687
688 IF ( LASSIM ) THEN
689   IF ( TRIM(CASSIM_ISBA) == "OI" ) THEN
690     IF ( IO%NPATCH /= 1 ) CALL ABOR1_SFX ('Reading of diagnostical values for'&
691                        & //'assimilation at the moment only works for one patch for OI')          
692     ! Diagnostic fields for assimilation
693     IF ( .NOT. ALLOCATED(XAT2M_ISBA)) ALLOCATE(XAT2M_ISBA(ILU,1))
694     XAT2M_ISBA=XUNDEF
695     YRECFM='T2M'
696     CALL IO_BUFF(YRECFM,'R',GKNOWN)
697     CALL READ_SURF(HPROGRAM,YRECFM,XAT2M_ISBA(:,1),IRESP)
698
699     IF ( .NOT. ALLOCATED(XAHU2M_ISBA)) ALLOCATE(XAHU2M_ISBA(ILU,1))
700     XAHU2M_ISBA=XUNDEF
701     YRECFM='HU2M'
702     CALL IO_BUFF(YRECFM,'R',GKNOWN)
703     CALL READ_SURF(HPROGRAM,YRECFM,XAHU2M_ISBA(:,1),IRESP)
704
705     IF ( .NOT. ALLOCATED(XAZON10M_ISBA)) ALLOCATE(XAZON10M_ISBA(ILU,1))
706     XAZON10M_ISBA=XUNDEF
707     YRECFM='ZON10M'
708     CALL IO_BUFF(YRECFM,'R',GKNOWN)
709     CALL READ_SURF(HPROGRAM,YRECFM,XAZON10M_ISBA(:,1),IRESP)
710
711     IF ( .NOT. ALLOCATED(XAMER10M_ISBA)) ALLOCATE(XAMER10M_ISBA(ILU,1))
712     XAMER10M_ISBA=XUNDEF
713     YRECFM='MER10M'
714     CALL IO_BUFF(YRECFM,'R',GKNOWN)
715     CALL READ_SURF(HPROGRAM,YRECFM,XAMER10M_ISBA(:,1),IRESP)
716   ELSEIF ( NIFIC/=NVAR+2 ) THEN
717     ! Diagnostic fields for EKF assimilation ("observations")
718     DO IOBS = 1,NOBSTYPE
719      SELECT CASE (TRIM(COBS(IOBS)))
720        CASE("T2M")
721          IF ( .NOT. ALLOCATED(XAT2M_ISBA)) ALLOCATE(XAT2M_ISBA(ILU,1))
722          XAT2M_ISBA=XUNDEF
723          YRECFM='T2M'
724          CALL IO_BUFF(YRECFM,'R',GKNOWN)
725          CALL READ_SURF(HPROGRAM,YRECFM,XAT2M_ISBA(:,1),IRESP)
726        CASE("HU2M")
727          IF ( .NOT. ALLOCATED(XAHU2M_ISBA)) ALLOCATE(XAHU2M_ISBA(ILU,1))
728          XAHU2M_ISBA=XUNDEF
729          YRECFM='HU2M'
730          CALL IO_BUFF(YRECFM,'R',GKNOWN)
731          CALL READ_SURF(HPROGRAM,YRECFM,XAHU2M_ISBA(:,1),IRESP)
732        CASE("WG1")
733          ! This is already read above
734        CASE("WG2")
735          ! This is already read above
736        CASE("LAI")
737          ! This is already read above   
738        CASE("SWE")
739          ! This is handled independently 
740        CASE DEFAULT
741          CALL ABOR1_SFX("Mapping of "//TRIM(COBS(IOBS))//" is not defined in READ_ISBA_n!")
742      END SELECT
743     ENDDO
744   ENDIF
745 ENDIF
746 !
747 DEALLOCATE(ZWORK)
748 !
749 DO JP = 1,IO%NPATCH
750   PK => NP%AL(JP)
751   PEK => NPE%AL(JP)
752   ALLOCATE(PEK%XSNOWFREE_ALB     (PK%NSIZE_P))
753   ALLOCATE(PEK%XSNOWFREE_ALB_VEG (PK%NSIZE_P))
754   ALLOCATE(PEK%XSNOWFREE_ALB_SOIL(PK%NSIZE_P))
755 ENDDO
756 !
757 IF (LHOOK) CALL DR_HOOK('READ_ISBA_N',1,ZHOOK_HANDLE)
758 !
759 CONTAINS
760 !
761 SUBROUTINE MAKE_ENS_ENKF(KWORK,KLU,HREC,PCOFSWI,NP,PVAR)
762 !
763 USE MODD_ASSIM, ONLY : LENS_GEN, XADDTIMECORR, XADDINFL, XASSIM_WINH
764 !
765 USE MODI_ADD_NOISE
766 USE MODE_RANDOM
767 !
768 IMPLICIT NONE
769 !
770 INTEGER, INTENT(IN) :: KWORK
771 INTEGER, INTENT(IN) :: KLU
772  CHARACTER(LEN=3), INTENT(IN) :: HREC
773 REAL, DIMENSION(:), INTENT(IN) :: PCOFSWI
774 TYPE(ISBA_NP_t), INTENT(INOUT) :: NP
775 REAL, DIMENSION(:,:,:), INTENT(INOUT), OPTIONAL :: PVAR
776 !
777  CHARACTER(LEN=LEN_HREC) :: YRECFM         ! Name of the article to be read
778  CHARACTER(LEN=4) :: YLVL
779  CHARACTER(LEN=3) :: YVAR
780 REAL, DIMENSION(KLU) :: ZVAR
781 REAL :: ZWHITE_NOISE, ZVAR0
782 INTEGER :: JL, JI, JP, IVAR
783 LOGICAL :: GPASS
784 !
785 REAL(KIND=JPRB) :: ZHOOK_HANDLE
786 !
787 IF (LHOOK) CALL DR_HOOK('READ_ISBA_N:MAKE_ENS_ENKF',0,ZHOOK_HANDLE)
788 !
789 !
790 DO JL=1,KWORK
791   !
792   IF (KWORK>1) THEN
793     WRITE(YLVL,'(I4)') JL
794     YRECFM = TRIM(HREC)//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
795   ELSE
796     YRECFM = TRIM(HREC)
797   ENDIF
798   !
799   IVAR = 0
800   DO JVAR = 1,NVAR
801     GPASS = ( TRIM(CVAR(JVAR))==TRIM(YRECFM) )
802     IF (GPASS) THEN
803       IVAR = JVAR
804       EXIT
805     ENDIF
806   ENDDO
807   !
808   IF ( GPASS ) THEN
809     !
810     IF (XADDINFL(IVAR)>0.) THEN
811       !
812       IF (LASSIM .OR. (.NOT.LENS_GEN .AND. XADDTIMECORR(IVAR)>0.)) THEN
813         !
814         WRITE(YVAR,'(I3)') IVAR
815         YRECFM='RD_NS'//ADJUSTL(YVAR(:LEN_TRIM(YVAR)))
816         CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
817         DO JP = 1,IO%NPATCH
818           CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NP%AL(JP)%XRED_NOISE(:,IVAR))
819         ENDDO           
820         !
821         IF (.NOT.LASSIM) THEN
822           !
823           DO JP = 1,IO%NPATCH
824             PK => NP%AL(JP)
825             DO JI = 1,NP%AL(JP)%NSIZE_P
826               IMASK = PK%NR_P(JI)
827               ZWHITE_NOISE = XADDINFL(IVAR)*PCOFSWI(IMASK)*RANDOM_NORMAL()
828               CALL ADD_NOISE(XADDTIMECORR(IVAR),XASSIM_WINH,ZWHITE_NOISE,PK%XRED_NOISE(JI,IVAR))
829             ENDDO
830           ENDDO
831           !
832           ZCOEF = XASSIM_WINH/24.
833           !
834         ENDIF
835         !
836       ELSE
837         !
838         DO JP = 1,IO%NPATCH
839           PK => NP%AL(JP)
840           DO JI = 1,NP%AL(JP)%NSIZE_P
841             IMASK = PK%NR_P(JI)
842             NP%AL(JP)%XRED_NOISE(JI,IVAR) = XADDINFL(IVAR)*PCOFSWI(IMASK)*RANDOM_NORMAL()
843           ENDDO
844         ENDDO
845         !
846         ZCOEF = 1. 
847         !
848       ENDIF
849       !
850       IF (.NOT.LASSIM) THEN
851         !
852         DO JP = 1,IO%NPATCH
853           !
854           ZVAR(:) = 0.
855           IF (TRIM(HREC)=='TG') THEN
856             ZVAR(1:NP%AL(JP)%NSIZE_P) = NPE%AL(JP)%XTG(:,JL)
857           ELSEIF (TRIM(HREC)=='WG') THEN
858             ZVAR(1:NP%AL(JP)%NSIZE_P) = NPE%AL(JP)%XWG(:,JL)
859           ELSEIF (TRIM(HREC)=='LAI' .AND. PRESENT(PVAR)) THEN
860             ZVAR(1:NP%AL(JP)%NSIZE_P) = PVAR(1:NP%AL(JP)%NSIZE_P,JL,JP)
861           ELSE
862             CALL ABOR1_SFX("READ_ISBAn: HREC "//HREC//" not permitted")
863           ENDIF
864           !
865           DO JI = 1,NP%AL(JP)%NSIZE_P
866             IF ( ZVAR(JI)/=XUNDEF ) THEN
867               !
868               ZVAR0 = ZVAR(JI)
869               !
870               ZVAR(JI) = ZVAR(JI) + ZCOEF * NP%AL(JP)%XRED_NOISE(JI,IVAR)
871               !
872               IF (ZVAR(JI) < 0.) THEN
873                 IF (LENS_GEN) THEN
874                   ZVAR(JI) = ABS(ZVAR(JI))
875                 ELSE
876                   ZVAR(JI) = ZVAR0
877                 ENDIF
878               ENDIF
879             ENDIF
880           ENDDO
881           !
882           IF (TRIM(HREC)=='TG') THEN
883             NPE%AL(JP)%XTG(:,JL) = ZVAR(1:NP%AL(JP)%NSIZE_P)
884           ELSEIF (TRIM(HREC)=='WG') THEN
885             NPE%AL(JP)%XWG(:,JL) = ZVAR(1:NP%AL(JP)%NSIZE_P)
886           ELSEIF (TRIM(HREC)=='LAI') THEN
887             PVAR(1:NP%AL(JP)%NSIZE_P,JL,JP) = ZVAR(1:NP%AL(JP)%NSIZE_P)
888           ENDIF
889           !
890         ENDDO
891         !
892       ENDIF
893       !
894     ENDIF
895     !
896   ENDIF
897   !
898 ENDDO
899 !
900 IF (LHOOK) CALL DR_HOOK('READ_ISBA_N:MAKE_ENS_ENKF',1,ZHOOK_HANDLE)
901 !
902 END SUBROUTINE MAKE_ENS_ENKF
903 !
904 !-------------------------------------------------------------------------------
905 !
906 END SUBROUTINE READ_ISBA_n