
C**********************************************************************C
C                                                                      C
C     *** SUBROUTINE Cabon balance ***                                 C
C     Simulate GEP AND NEE for selected HUC                            C
C     WRITE GEP AND NEE TO MONTHCARBON.TXT, ANNUALCARBON.TXT, HUCCARBON.TXT        C
C     SIMULATE BIODIVERSITY FOR SELECTED HUC                           C
C     WRITE BIODIVERSITY TO HUCBIO.TXT                                 C
C                                                                      C
C**********************************************************************C
C        I=HUC; J= YEAR, M =MONTH, MNDAY= NO OF DAYS IN A MONTH

      SUBROUTINE CARBONBAL
      
      COMMON/BASIC/NGRID,NYEAR,NLC,BYEAR,IYSTART,IYEND,POP_FLAG,REGION,
     &CLIMATE_FLAG

      COMMON/HUC/ HUCAREA(2200)
      COMMON/CELLINFO/LADUSE(4000,10),HUCNO(4000),
     >                LATUDE(4000),LONGI(4000),TAREA(4000)
    
    
      COMMON/CARBON/ GEPM(2200, 200, 12),RECOM(2200, 200, 12), 
     >               NEEM(2200,200,12)
      
              
C---BIODIVERSITY CALCULATION USING ANNUAL AET, PET      
      
      COMMON/HUCPETAET/HUCAET(2200,200), HUCPET(2200,200)  
          
C --------------------------------------------------------------
      INTEGER I,J, M
      
      
C ---------------------------------------------------------------      
                     
      INTEGER HUCNO(4000), BYEAR
      
      INTEGER IYSTART, IYEND, A, R
                   
      REAL HUCAREA(2200),TAREA(4000)

      REAL ANGEP(200), ANRECO(200), ANNEE(200),  HUCGEP, HUCNEE,HUCRE
      
      REAL HUCAET(2200,200), HUCPET(2200,200)
      
      REAL HUCTRS,TRS(200), HUCMAMMALS,MAMMALS(200),HUCBIRD,BIRD(200),
     &HUCAMPHIB,AMPHIB(200),HUCREPTILES, REPTILES(200), HUCVERTEB, 
     &VERTEB(200)
     
     
      REAL GEPM(2200, 200, 12),RECOM(2200,200,12), 
     >               NEEM(2200,200,12)
      
      CHARACTER*5 CLIMATE_FLAG
      
      CHARACTER*2 REGION
      
      INTEGER BIO_FLAG
      
      REAL MOAVGGEP(2200,12),MOAVGRES(2200,12),MOAVGNEE(2200,12)

C ----------------------------------------------------------------      


        WRITE (400, 500) 
500    FORMAT ('CELL,YEAR,MONTH,GEP(gC_M2_MO),RES,NEE')


        WRITE (500, 600) 
600    FORMAT ('CELL,YEAR,GEP(gC_M2_YR),RES,NEE,', 
     >        'AET,PET')


        WRITE (600, 650) 
650    FORMAT ('CELL,YRS,GEP(gC_M2_YR),RES,NEE')

        WRITE (806, 651) 
651    FORMAT ('CELL,MONTH,GEP(gC_M2_MO),RES,NEE')

       
      BIO_FLAG = 0
       
      IF(BIO_FLAG.EQ.1) THEN

         WRITE(*,70)
70       FORMAT('  CALCULATING BIODIVERSITY.....'/)   

      ENDIF  
           
           
      DO 300 I=1, NGRID  
      
        DO 301 A=1,12
      
        MOAVGGEP(I,A)=0.
        MOAVGRES(I,A)=0.
        MOAVGNEE(I,A)=0.
                
301   CONTINUE       
                
                         
         DO 200 J=(IYSTART-BYEAR+1),(IYEND-BYEAR+1)
   
            IDY = J + BYEAR - 1 

            ANGEP(J) = 0.
            ANRECO(J) = 0.
            ANNEE(J) = 0.
                    
C----ONLY WRITE OUTPUT FOR YEARS AFTER THE MODEL WARMUP PERIOD

      IF(J.GT.(IYSTART-BYEAR)) THEN

            DO 100 M=1, 12

C -- WRITE MONTHLY CARBON DATA TO MONTHCARBON.TXT
       	 
               WRITE (400, 2000) HUCNO(I), IDY, M, GEPM(I,J, M), 
     >          RECOM(I,J,M), NEEM(I,J,M)         
     
          
2000           FORMAT (I12,',',I12,',',I12,',',F10.2,',',F10.2,',',
     &F10.2)
        
C--- ACCUMULATE ANNUAL GEP FROM MONTHLY VALUES (g Co2/m2/mon)

            ANGEP(J) = ANGEP(J) + GEPM(I,J,M)  
            ANRECO(J) = ANRECO(J) + RECOM(I,J,M) 
            ANNEE(J) = ANNEE(J) + NEEM(I,J,M)     
            
            MOAVGGEP(I,M)=MOAVGGEP(I,M)+GEPM(I,J,M)
            MOAVGRES(I,M)=MOAVGRES(I,M)+RECOM(I,J,M)
            MOAVGNEE(I,M)=MOAVGNEE(I,M)+NEEM(I,J,M)
                             
           
           
100         CONTINUE 

          ENDIF
           
c --- convert annual GEP to NEE. ref. Law et al (2002) USING THE AVG OF FORESTS

c            ANNEE(J) = ((285 - 0.44 * ANGEP(J)) + 
C     &(618 - 0.67*ANGEP(J)))/2.

      IF(BIO_FLAG.EQ.1) THEN

C --- CALCULATING Tree Species Richness (TSR) USING CURRIE (1987) NATURE
           
            TRS(J) = 185.8/(1.0 + EXP (3.09-0.00432*HUCAET(I,J)))
           
C --- CALCULATING  (TSR) USING CURRIE (1991) for different
C--The American Society of Naturalists Energy and Large-Scale Patterns of Animal- and Plant-Species Richness 
C --Author(s): David J. Currie Source: The American Naturalist, Vol. 137, No. 1 (Jan., 1991), pp. 27-49 Published by: The University of Chicago Press for The American Society 

            MAMMALS(J) = 1.12*(1.0 -EXP(-0.00348*HUCPET(I,J) )) + 0.653 
          
            MAMMALS(J) = 10**(MAMMALS(J)) - 1.
           
            IF (HUCPET(I,J) .LT. 525) THEN 
            
               BIRD(J) = 1.4 + 0.00159 * HUCPET(I,J)
               
            ELSE 
            
               BIRD(J) = 2.26 - 0.0000256 * HUCPET(I,J)
               
            ENDIF

            BIRD(J) = 10**(BIRD(J)) - 1.
           
            IF (HUCPET(I,J) .LE. 200) THEN 
            
               AMPHIB(J) = 0.
               
            ELSE 
            
               AMPHIB(J) = 3.07*(1.0-Exp(-0.00315*HUCPET(I,J) ))
            
            ENDIF
           
            AMPHIB(J) = exp(AMPHIB(J)) -1.

            IF (HUCPET(I,J) .LT. 400) THEN 
            
               REPTILES(J) = 0.
               
            ELSE 
            
               REPTILES(J) = 5.21*(1.0-Exp(-0.00249*HUCPET(I,J) ))
     & - 3.347
           
            ENDIF
                                
           REPTILES(J) = 10**(REPTILES(J)) -1. 

           
           VERTEB(J) = 1.49*(1.0-Exp(-0.00186*HUCPET(I,J) )) + 0.746
           
           VERTEB(J) = 10**(VERTEB(J)) -1.
           
C----ONLY WRITE OUTPUT FOR YEARS AFTER THE MODEL WARMUP PERIOD

      IF(J.GT.(IYSTART-BYEAR)) THEN
      
C ---- write annual biodiversity results TO ANNUALBIO.TXT
     
              WRITE (700, 3100) HUCNO(I), IDY, TRS(J),MAMMALS(J), 
     &BIRD(J), AMPHIB(J), REPTILES(J), VERTEB(J), 
     &HUCAET(I,J), HUCPET(I,J)
                                              
3100          FORMAT (I12, ',', I12, ',', 
     >           F10.2, ',', F10.2, ',', F10.2, ',', F10.2, ',', F10.2,
     >           ',', F10.2, ',', F8.1, ',', F8.1)    
     
      ENDIF  
      
      ENDIF
                 
           

C----ONLY WRITE OUTPUT FOR YEARS AFTER THE MODEL WARMUP PERIOD

      IF(J.GT.(IYSTART-BYEAR)) THEN


C ---  write annual GEP and NEE TO ANNUALCARBON.TXT  
    

              WRITE (500, 3000) HUCNO(I), IDY, ANGEP(J),ANRECO(J),
     >                    ANNEE(J), HUCAET(I,J), HUCPET(I,J)
                                              
3000          FORMAT (I12, ',', I12, ',',F16.2, ',', 
     >           F16.2, ',', F16.2, ',', F16.2,',', F16.2)

      ENDIF

200      CONTINUE

C -- WRITE AVERAGE ANNUAL CARBON BETWEEN IYSTART AND IYEND TO AVGANNCARBON.TXT
C -- WRITE AVERAGE ANNUAL BIODIVERSITY BETWEEN IYSTART AND IYEND TO AVGANNBIO.TXT

         HUCGEP=0.
         HUCRE=0.
         HUCNEE=0.
         HUCTRS=0.
         HUCMAMMALS=0.
         HUCBIRD=0.
         HUCAMPHIB=0.
         HUCREPTILES=0.
         HUCVERTEB=0.
         
         IY=0
         
C----ONLY WRITE OUTPUT FOR YEARS AFTER THE MODEL WARMUP PERIOD

      IF(J.GT.(IYSTART-BYEAR)) THEN

         DO 75 J=IYSTART,IYEND

              HUCGEP = HUCGEP + ANGEP(J-(BYEAR-1))
              
              HUCRE = HUCRE + ANRECO(J-(BYEAR-1))
              
              HUCNEE = HUCNEE + ANNEE(J-(BYEAR-1))
            
              IF(BIO_FLAG.EQ.1) THEN
              
                 HUCTRS = HUCTRS + TRS(J-(BYEAR-1))
                 HUCMAMMALS = HUCMAMMALS + MAMMALS(J-(BYEAR-1))
                 HUCBIRD = HUCBIRD + BIRD(J-(BYEAR-1))
                 HUCAMPHIB = HUCAMPHIB + AMPHIB(J-(BYEAR-1))
                 HUCREPTILES = HUCREPTILES + REPTILES(J-(BYEAR-1))
                 HUCVERTEB = HUCVERTEB + VERTEB(J-(BYEAR-1))
              
              ENDIF

         
         IY = IY+1
         
75       CONTINUE

         DO 101 R=1,12
      
            MOAVGGEP(I,R)=MOAVGGEP(I,R)/IY
            MOAVGRES(I,R)=MOAVGRES(I,R)/IY
            MOAVGNEE(I,R)=MOAVGNEE(I,R)/IY
     
101      CONTINUE 
         
         HUCGEP = HUCGEP/IY
         HUCRE = HUCRE/IY
         HUCNEE = HUCNEE/IY
          

         WRITE (600, 4000) HUCNO(I),IY, HUCGEP, HUCRE, HUCNEE
                                              
4000     FORMAT (I12, ',', I12, ',',F14.2, ',', F14.2, ',', F14.2)       
     
         IF(BIO_FLAG.EQ.1) THEN
         
            HUCTRS = HUCTRS/IY
            HUCMAMMALS = HUCMAMMALS/IY
            HUCBIRD = HUCBIRD/IY
            HUCAMPHIB = HUCAMPHIB/IY
            HUCREPTILES = HUCREPTILES/IY
            HUCVERTEB = HUCVERTEB/IY        
         


            WRITE (800, 4200) HUCNO(I),IY, HUCTRS, HUCMAMMALS,
     >        HUCBIRD, HUCAMPHIB, HUCREPTILES, HUCVERTEB 
                                              
4200        FORMAT (I12, ',', I12, ',', 
     >           F10.2, ',', F10.2, ',', F10.2, ',', F10.2, ',', F10.2,
     >           ',', F10.2)
     
         ENDIF

       DO 252 R=1,12
      
      WRITE(806,253) HUCNO(I),R,MOAVGGEP(I,R),MOAVGRES(I,R),
     &MOAVGNEE(I,R)
     
253   FORMAT(I10,',', I10,',',F14.2, ',', F14.2, ',', F14.2)
     
252   CONTINUE

      ENDIF

300   CONTINUE

      RETURN
      END
