source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/stomate_check.f90 @ 6737

Last change on this file since 6737 was 5149, checked in by chao.yue, 6 years ago

Add mass balance check for land use change
File size: 12.0 KB
Line 
1! =================================================================================================================================
2! MODULE       : stomate_check
3!
4! CONTACT      : orchidee-help _at_ ipsl.jussieu.fr
5!
6! LICENCE      : IPSL (2006)
7! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8!
9!>\BRIEF       Subroutines to check mass conservation
10!!     
11!!\n DESCRIPTION: None
12!!
13!! RECENT CHANGE(S): None
14!!
15!! SVN :
16!! $HeadURL: $
17!! $Date: $
18!! $Revision: $
19!! \n
20!================================================================================================================================
21
22
23MODULE stomate_check
24
25  USE xios_orchidee
26  USE constantes
27  USE stomate_data
28  USE ioipsl_para
29  USE constantes_soil
30  USE ieee_arithmetic 
31
32  IMPLICIT NONE
33
34  ! private routines
35  PRIVATE
36
37  ! public routines
38  PUBLIC stomate_check_mass_values, stomate_check_cons_mass
39
40  INTERFACE stomate_check_cons_mass
41    MODULE PROCEDURE stomate_check_cons_mass_r2d, stomate_check_cons_mass_r3d
42  END INTERFACE
43
44
45CONTAINS 
46
47
48!! ================================================================================================================================
49!! SUBROUTINE   : stomate_check_mass
50!!
51!>\BRIEF       Check for biomass variable consistency 
52!!
53!! DESCRIPTION  : Check for biomass values: negatives, very big values and Nan
54!!
55!! RECENT CHANGE(S) : None
56!!
57!! MAIN OUTPUT VARIABLE(S): None
58!!
59!! REFERENCES   : None
60!!
61!! FLOWCHART    : None
62!! \n
63!_ ================================================================================================================================
64  SUBROUTINE stomate_check_mass_values(npts,biomass,identifier)
65
66  INTEGER(i_std), INTENT(in)                    :: npts         !! Domain size (unitless)
67  REAL(r_std), DIMENSION(:,:,:,:),INTENT(in)    :: biomass      !! Biomass components of the model tree 
68 
69  CHARACTER(LEN=*), INTENT(in),OPTIONAL         ::  identifier  !! string with to identify where this routine was called form
70 
71  ! locals
72  INTEGER                                       :: nelement, npart, pft, idx
73  LOGICAL                                       :: is_crash 
74  CHARACTER(LEN=:), ALLOCATABLE                 :: ident_routine
75! ===================================================================================================
76
77
78   ident_routine='No information'
79   IF (PRESENT(identifier)) THEN
80     ident_routine=identifier
81   ENDIF
82
83   ! 1. check for negative biomass pools:
84   IF (printlev >= 1) WRITE (numout,*) '=== BIOMASS CHECK ==='
85
86   is_crash=.FALSE.
87   DO nelement=1,nelements
88     DO npart=1,nparts
89        ! dont check crops and crops, it does not work yet
90        DO pft=1,nvm-2
91           DO idx=1,npts
92             IF ((biomass(idx,pft,npart,nelement) .LT. zero).OR. & ! negative
93                   (biomass(idx,pft,npart,nelement) .GT. large_value).OR. & ! infinity
94                   (isnan(biomass(idx,pft,npart,nelement))))   THEN ! NaN
95                        WRITE (numout,*) 'FATAL negative biomass detected'
96                        WRITE (numout,*) 'routine:'  , ident_routine
97                        WRITE (numout,*) 'biomass = ', biomass(idx,pft,npart,nelement)
98                        WRITE (numout,*) 'box     = ',idx
99                        WRITE (numout,*) 'PFT     = ',pft
100                        WRITE (numout,*) 'npart   = ',npart
101                        WRITE (numout,*) '1=leaf,2=sapabov,3=sapbelo,4=heartabov,5=heatbelow,6=root,7=fruit,8=carbres,9=ilabile'
102                        WRITE (numout,*) 'element = ',nelement
103                        is_crash=.TRUE.
104             ENDIF
105           ENDDO
106        ENDDO
107     ENDDO
108   ENDDO
109
110
111   IF (is_crash) THEN
112      CALL ipslerr_p ( 3, 'stomate_check_mass_values',   &
113                          'Negative biomass pool was detected:', &
114                          'Failed for subroutine ', &
115                          identifier)
116   ELSE
117      IF (printlev >= 1) WRITE (numout,*) '     passed          '
118   ENDIF
119
120  END SUBROUTINE stomate_check_mass_values
121
122!! ================================================================================================================================
123!! SUBROUTINE   : stomate_check_cons_mass_r3d
124!!
125!>\BRIEF       
126!!
127!! DESCRIPTION  : None
128!!
129!! RECENT CHANGE(S) : None
130!!
131!! MAIN OUTPUT VARIABLE(S): None
132!!
133!! REFERENCES   : None
134!!
135!! FLOWCHART    : None
136!! \n
137!_ ================================================================================================================================
138  SUBROUTINE stomate_check_cons_mass_r3d(lalo, mass_before,      &
139                          mass_after,       &
140                          mass_change,      &
141                          identifier)
142
143  REAL(r_std),DIMENSION(:,:),INTENT(in)              :: lalo          !! Geographical coordinates (latitude,longitude) for pixels (degrees)
144  REAL(r_std), DIMENSION(:,:,:),INTENT(in)           :: mass_before   !! Biomass old
145  REAL(r_std), DIMENSION(:,:,:),INTENT(in)           :: mass_after    !! Biomass new (npts,nvm,nelements)
146  REAL(r_std), DIMENSION(:,:,:),INTENT(in)           :: mass_change   !! Biomass change
147  CHARACTER(LEN=*), INTENT(in),OPTIONAL              ::  identifier   !! string with to identify where this routine was called form
148
149  ! locals
150  INTEGER                                            :: nelement, pft, idx, ierr
151  LOGICAL,ALLOCATABLE, DIMENSION(:,:,:)              :: is_crash
152  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)         :: mass_conserv  !! conservation
153  CHARACTER(LEN=:), ALLOCATABLE                      :: ident_routine
154  ! this subroutine checks for mass conservation
155  INTEGER(i_std)                                     :: npts          !! Domain size (unitless)
156  INTEGER(i_std)                                     :: nvm           !! Domain size (unitless)
157  INTEGER(i_std)                                     :: nelements     !! Domain size (unitless)
158 
159! ===================================================================================================
160
161  npts = SIZE(mass_before, DIM=1)
162  nvm = SIZE(mass_before, DIM=2)
163  nelements = SIZE(mass_before, DIM=3)
164
165  IF (PRESENT(identifier)) THEN
166      ident_routine=identifier
167  ELSE
168       ident_routine=' no information'
169  ENDIF
170
171  ALLOCATE(is_crash(npts, nvm, nelements), stat=ierr)
172  IF (ierr /= 0) CALL ipslerr_p(3, 'stomate_chekc_cons_mass_r3d', 'Memory allocation problem with ', 'is_crash variable', '')
173
174  ALLOCATE(mass_conserv(npts, nvm, nelements), stat=ierr)
175  IF (ierr /= 0) CALL ipslerr_p(3, 'stomate_check_cons_mass_r3d', 'Memory allocation problem with ', 'mass_conserv variable', '')
176
177  is_crash(:,:,:)=.FALSE.
178  ! bookkeeping:
179  mass_conserv(:,:,:)          = mass_before(:,:,:) - mass_after(:,:,:) + mass_change(:,:,:)
180
181  ! check the bookkeeping
182  DO nelement=1, nelements-1
183    ! take out crops and grass; till they are fixed
184    DO pft=1, nvm-2
185      DO idx=1, npts
186         IF (ABS(mass_conserv(idx,pft,nelement)).GT. (min_stomate*100.)) THEN
187                WRITE (numout,*) 'FATAL mass conservation failed (positive value = leak)'
188                WRITE (numout,*) 'routine:', ident_routine
189                ! the limit is given by the allocation routine, which accepts
190                ! mass leaks<ABS(10E-6)
191                WRITE (numout,*) ' limit: ',min_stomate*100.
192                WRITE (numout,*) ' Coordinates :', lalo(idx,:) 
193                WRITE (numout,*) ' for element :', nelement, ' of nelements: ',nelements
194                WRITE (numout,*) ' gridpoint: ',idx , ' of ngrids: ',npts
195                WRITE (numout,*) ' PFT: ',pft , ' of npfts: ',nvm
196                WRITE (numout,*) ' mismatch =', mass_conserv(idx,pft,nelement)
197                WRITE (numout,*) ' mass(before) =', mass_before(idx,pft,nelement)
198                WRITE (numout,*) ' mass(after) =', mass_after(idx,pft,nelement)
199                is_crash(idx,pft,nelement)=.TRUE.
200         ENDIF 
201      ENDDO
202    ENDDO
203  ENDDO
204
205  IF (ANY(is_crash(:,:,:))) THEN
206     CALL ipslerr_p ( 3, 'stomate_check_cons_mass_r3d', &
207                    'Mass conservation is not preserved', &
208          &         'in subroutine', &
209                    ident_routine)
210  ENDIF
211
212  DEALLOCATE(mass_conserv)
213  DEALLOCATE(is_crash)
214
215  END SUBROUTINE stomate_check_cons_mass_r3d
216
217!! ================================================================================================================================
218!! SUBROUTINE   : stomate_check_cons_mass_r2d
219!!
220!>\BRIEF       
221!!
222!! DESCRIPTION  : None
223!!
224!! RECENT CHANGE(S) : None
225!!
226!! MAIN OUTPUT VARIABLE(S): None
227!!
228!! REFERENCES   : None
229!!
230!! FLOWCHART    : None
231!! \n
232!_ ================================================================================================================================
233  SUBROUTINE stomate_check_cons_mass_r2d(lalo, mass_before,      &
234                          mass_after,       &
235                          mass_change,      &
236                          identifier)
237
238
239  REAL(r_std),DIMENSION(:,:),INTENT(in)              :: lalo          !! Geographical coordinates (latitude,longitude) for pixels (degrees)
240  ! biomass pools to be checked
241  REAL(r_std), DIMENSION(:,:),INTENT(in)             :: mass_before   !! Biomass old
242  REAL(r_std), DIMENSION(:,:),INTENT(in)             :: mass_after    !! Biomass new (npts,nvm)
243  REAL(r_std), DIMENSION(:,:),INTENT(in)             :: mass_change   !! Biomass change
244
245  CHARACTER(LEN=*), INTENT(in),OPTIONAL              ::  identifier   !! string with to identify where this routine was called form
246
247  ! locals
248  INTEGER(i_std)                                     :: iele, idx, ierr
249  LOGICAL,ALLOCATABLE,DIMENSION(:,:)                 :: is_crash
250  REAL(r_std),ALLOCATABLE, DIMENSION(:,:)            :: mass_conserv  !! conservation
251  CHARACTER(LEN=:), ALLOCATABLE                      :: ident_routine
252  ! this subroutine checks for mass conservation
253  INTEGER(i_std)                                     :: npts          !! Domain size (unitless)
254  INTEGER(i_std)                                     :: nelements           !! Domain size (unitless)
255! ===================================================================================================
256  npts = SIZE(mass_before, DIM=1)
257  nelements = SIZE(mass_before, DIM=2)
258
259  ALLOCATE(is_crash(npts, nelements), stat=ierr)
260  IF (ierr /= 0) CALL ipslerr_p(3, 'stomate_check_cons_mass_r2d', 'Memory allocation problem with ', 'is_crash variable', '')
261
262  ALLOCATE(mass_conserv(npts, nelements), stat=ierr)
263  IF (ierr /= 0) CALL ipslerr_p(3, 'stomate_check_cons_mass_r2d', 'Memory allocation problem with ', 'mass_conserv variable', '')
264
265  IF (PRESENT(identifier)) THEN
266      ident_routine=identifier
267  ELSE
268       ident_routine=' no information'
269  ENDIF
270
271
272  is_crash(:,:)=.FALSE.
273  ! bookkeeping:
274  mass_conserv(:,:)          = mass_before(:,:) - mass_after(:,:) + mass_change(:,:)
275
276  ! check the bookkeeping
277  ! take out crops and grass; till they are fixed
278  DO iele=1,nelements
279    DO idx=1,npts
280       IF (ABS(mass_conserv(idx,iele)).GT. (min_stomate*100.)) THEN
281              WRITE (numout,*) 'FATAL mass conservation failed (positive value = leak)'
282              WRITE (numout,*) 'routine:', ident_routine
283                ! the limit is given by the allocation routine, which accepts
284                ! mass leaks<ABS(10E-6)
285              WRITE (numout,*) ' limit: '       ,min_stomate*100.
286              WRITE (numout,*) ' Coordinates :' , lalo(idx,:) 
287              WRITE (numout,*) ' gridpoint: '   , idx , ' of ngrids: ',npts
288              WRITE (numout,*) ' element: '     , iele , ' of nelements: ',nelements
289              WRITE (numout,*) ' mismatch ='    , mass_conserv(idx,iele)
290              WRITE (numout,*) ' mass(before) =', mass_before(idx,iele)
291              WRITE (numout,*) ' mass(after) =' , mass_after(idx,iele)
292              is_crash(idx,iele)=.TRUE.
293       ENDIF 
294    ENDDO
295  ENDDO
296
297  IF (ANY(is_crash(:,:))) THEN
298     CALL ipslerr_p ( 3, 'stomate_check_cons_mass_r2d', &
299                    'Mass convervation is not preserved', &
300          &         'in subroutine', &
301                    ident_routine)
302  ENDIF
303
304  DEALLOCATE(mass_conserv)
305  DEALLOCATE(is_crash)
306
307  END SUBROUTINE stomate_check_cons_mass_r2d
308
309
310END MODULE stomate_check
Note: See TracBrowser for help on using the repository browser.