!============================================================================== ! This source code is part of the ! Australian Community Atmosphere Biosphere Land Exchange (CABLE) model. ! This work is licensed under the CSIRO Open Source Software License ! Agreement (variation of the BSD / MIT License). ! ! You may not use this file except in compliance with this License. ! A copy of the License (CSIRO_BSD_MIT_License_v2.0_CABLE.txt) is located ! in each directory containing CABLE code. ! ! ============================================================================== ! Purpose: module for land-use change which interacts with POP demography ! via secondary forest age-distribution, and updates casa stocks according to land-use transitions ! ! Called from: cable_serial or cable_mpimaster ! ! SUBROUTINES ! ZeroPOPLUC(POPLUC) ! execute_luc_event(from_state,to_state,frac_change_grid,g,POPLUC) ! CALCULATE_WEIGHTS(POPLUC, g) ! INCREMENT_AGE(POPLUC,g) ! POPLUCStep(POPLUC,year) ! POPLUC_weights_transfer(POPLUC,POP,LUC_EXPT) ! POP_LUC_CASA_transfer(POPLUC,POP,LUC_EXPT,casapool,casabal,casaflux,ktauday) ! POPLUC_Init(POPLUC,LUC_EXPT, casapool, casaflux, casabiome, veg, POP, np) ! POPLUC_set_patchfrac(POPLUC,LUC_EXPT) ! POPLUC_set_params(POPLUC,LUC_EXPT) ! alloc_POPLUC(POPLUC, arraysize) ! WRITE_LUC_OUTPUT_NC ( POPLUC, ctime, FINAL ) ! WRITE_LUC_RESTART_NC ( POPLUC, ctime ) ! READ_LUC_RESTART_NC (POPLUC) ! WRITE_LUC_OUTPUT_GRID_NC ( POPLUC, ctime, FINAL ) ! History: Vanessa Haverd July 2016 ! ============================================================================== !=============================================================================== MODULE POPLUC_CONSTANTS USE TYPEdef, ONLY: dp, i4b INTEGER(i4b),PARAMETER :: LENGTH_SECDF_HISTORY = 4000 INTEGER(i4b),PARAMETER :: AGE_MAX = 1000 INTEGER(i4b),PARAMETER :: disturbance_interval = 100 ! N.B. needs to be the same as veg%disturbance_interval LOGICAL, PARAMETER :: IFHARVEST=.FALSE. INTEGER(i4b), PARAMETER :: ROTATION=70 INTEGER(i4b), PARAMETER :: nLU=3 ! number of land-use tiles (pf, sf, grass) INTEGER(i4b), PARAMETER :: nTrans=4 ! number of possible gross transition types (ptog, ptos, stog, gtos) END MODULE POPLUC_CONSTANTS !=============================================================================== MODULE POPLUC_Types USE TYPEdef, ONLY: dp, i4b USE POPLUC_Constants, ONLY: LENGTH_SECDF_HISTORY, AGE_MAX TYPE POPLUC_TYPE INTEGER(i4b),POINTER :: it INTEGER(i4b),POINTER :: np INTEGER(i4b),POINTER :: firstyear INTEGER(i4b),POINTER :: thisyear INTEGER(i4b), DIMENSION(:),POINTER :: n_event ! number of secondary forest transitions REAL(dp), DIMENSION(:),POINTER :: latitude, longitude REAL(dp), DIMENSION(:),POINTER :: primf, secdf, grass, & ! land cover types ptos, ptog, STOP, stog, gtop, gtos, & ! transitions frac_primf, frac_forest REAL(dp), DIMENSION(:,:),POINTER :: freq_age_primary, freq_age_secondary, & biomass_age_primary, biomass_age_secondary REAL(dp), DIMENSION(:,:),POINTER :: age_history_secdf, area_history_secdf REAL(dp), DIMENSION(:,:),POINTER :: FNEP, Clitt, Csoil, Cbiomass REAL(dp), DIMENSION(:,:),POINTER :: FHarvest, FClearance, FTransferNet REAL(dp), DIMENSION(:,:),POINTER :: FTransferGross REAL(dp), DIMENSION(:),POINTER :: pharv, smharv, syharv REAL(dp), DIMENSION(:,:),POINTER :: HarvProd, ClearProd REAL(dp), DIMENSION(:,:),POINTER :: fracHarvProd, fracClearProd REAL(dp), DIMENSION(:,:),POINTER :: HarvProdLoss, ClearProdLoss REAL(dp), DIMENSION(:),POINTER :: fracHarvResid, fracHarvSecResid, fracClearResid END TYPE POPLUC_TYPE END MODULE POPLUC_Types !=============================================================================== MODULE POPLUC_Module !------------------------------------------------------------------------------- ! * This module contains all subroutines for POPLUC calcs at a single time step. !------------------------------------------------------------------------------- USE TYPEdef, ONLY: sp, i4b USE POPLUC_Types USE POPLUC_Constants USE casavariable, ONLY: casa_pool, casa_balance, casa_flux, casa_biome USE POP_Types, ONLY: POP_TYPE USE cable_common_module, ONLY: cable_user USE cable_IO_vars_module, ONLY: landpt, patch USE CABLE_LUC_EXPT, ONLY: LUC_EXPT_TYPE USE POPModule, ONLY: pop_init_single CONTAINS !============================================================================= SUBROUTINE ZeroPOPLUC(POPLUC) TYPE(POPLUC_TYPE), INTENT(INOUT) :: POPLUC INTEGER:: g,np np = popluc%np POPLUC%firstyear = 0 POPLUC%thisyear = 0 POPLUC%primf = 0 POPLUC%secdf = 0 POPLUC%grass = 0 POPLUC%ptos = 0 POPLUC%ptog = 0 POPLUC%stop = 0 POPLUC%stog = 0 POPLUC%gtop = 0 POPLUC%gtos = 0 POPLUC%frac_forest = 0 POPLUC%frac_primf = 0 POPLUC%area_history_secdf = 0 POPLUC%age_history_secdf = 0 POPLUC%n_event = 0 POPLUC%freq_age_secondary = 0 POPLUC%freq_age_primary = 0 POPLUC%biomass_age_primary = 0 POPLUC%biomass_age_secondary = 0 POPLUC%FNEP = 0 POPLUC%Clitt = 0 POPLUC%Csoil = 0 POPLUC%Cbiomass = 0 POPLUC%FHarvest = 0 POPLUC%FClearance = 0 POPLUC%FTransferNet = 0 POPLUC%FTransferGross = 0 POPLUC%pharv = 0 POPLUC%smharv = 0 POPLUC%syharv = 0 POPLUC%HarvProd = 0 POPLUC%HarvProdLoss = 0 POPLUC%ClearProd = 0 POPLUC%ClearProdLoss = 0 END SUBROUTINE ZeroPOPLUC !============================================================================ SUBROUTINE execute_luc_event(from_state,to_state,frac_change_grid,g,POPLUC) ! Execute a transition between land use types (states) ! frac_change_grid = fractional change in unit fraction of grid cell this year IMPLICIT NONE CHARACTER(5), INTENT(IN) :: from_state, to_state REAL(dp), INTENT(INOUT):: frac_change_grid TYPE(POPLUC_TYPE),INTENT(INOUT) :: POPLUC INTEGER(i4b), INTENT(IN) :: g ! grid cell index REAL :: frac_open_grid, remaining INTEGER(i4b) :: n ! position of new element in POPLUC%SecFor array INTEGER(i4b) :: i frac_open_grid = 1.0 -POPLUC%frac_forest(g) n = POPLUC%n_event(g) IF (from_state=='PRIMF') THEN IF(frac_change_grid.GT.POPLUC%primf(g)) THEN !PRINT*, "Warning: requested reduction in primary forest area & ! exceeds primary forest area" IF (to_state=='SECDF') POPLUC%ptos(g) = POPLUC%primf(g) IF (to_state=='C3ANN') POPLUC%ptog(g) = POPLUC%primf(g) frac_change_grid = POPLUC%primf(g) POPLUC%primf(g) = 0.0 ELSE POPLUC%primf(g) = POPLUC%primf(g) & -frac_change_grid ENDIF IF (to_state=='SECDF') THEN ! Transition from primary -> secondary forest(ptos) POPLUC%n_event(g) = POPLUC%n_event(g)+1 n = POPLUC%n_event(g) POPLUC%area_history_secdf(g,n) = frac_change_grid POPLUC%age_history_secdf(g,n) = 0 POPLUC%freq_age_secondary(g,1) = POPLUC%freq_age_secondary(g,1) + frac_change_grid ENDIF ELSEIF (from_state=='SECDF') THEN IF (to_state=='PRIMF') THEN PRINT*, "Error: cannot create primary forest from secondary forest" STOP ELSE ! Transition from secondary -> non forest (stog) ! Assumption: youngest stands cleared first remaining = frac_change_grid i = 1 DO WHILE (remaining > 0.0 .AND. i <= age_max ) IF (POPLUC%freq_age_secondary(g,i).GE.remaining) THEN POPLUC%freq_age_secondary(g,i) =POPLUC%freq_age_secondary(g,i) & - remaining remaining = 0.0 ELSE remaining = remaining - POPLUC%freq_age_secondary(g,i) POPLUC%freq_age_secondary(g,i) = 0.0 i = i+1 ENDIF ENDDO !if (remaining.gt.frac_change_grid) POPLUC%stog(g) = POPLUC%stog(g)-remaining IF (remaining.GT.0.0) POPLUC%stog(g) = POPLUC%stog(g)-remaining ENDIF ELSEIF (to_state=='SECDF') THEN POPLUC%n_event(g) = POPLUC%n_event(g)+1 n = POPLUC%n_event(g) ! Transition from non-forest to secondary forest (gtos) IF (frac_change_grid.LE.frac_open_grid) THEN POPLUC%area_history_secdf(g,n) = frac_change_grid POPLUC%age_history_secdf(g,n) = 0 POPLUC%freq_age_secondary(g,1) = POPLUC%freq_age_secondary(g,1) + frac_change_grid ELSE POPLUC%area_history_secdf(g,n) = frac_open_grid POPLUC%age_history_secdf(g,n) = 0 POPLUC%freq_age_secondary(g,1) = POPLUC%freq_age_secondary(g,1) + frac_open_grid ! gtos to frac_change_grid here! POPLUC%gtos(g) = frac_open_grid ENDIF ELSEIF (to_state=='PRIMF') THEN PRINT*, "Error: cannot create primary forest from non-forest" STOP ENDIF ENDSUBROUTINE execute_luc_event !============================================================================ SUBROUTINE CALCULATE_WEIGHTS(POPLUC, g) ! Calculates weights (fraction of total forest area on grid cell) !for primary and secondary forest stands up to specified maximum stand age IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(INOUT) :: POPLUC INTEGER(i4b), INTENT(IN) :: g REAL(dp), PARAMETER :: EPS= 1.e-9 INTEGER(i4b) :: age , i, iage REAL(dp) :: fac, disturbance_freq ! First get relative weights for primary forest disturbance_freq=1.0/REAL(disturbance_interval,dp) !fac = POPLUC%frac_primf/POPLUC%frac_forest fac = 1.0 DO iage = 1, age_max POPLUC%freq_age_primary(g,iage) = REALExponential(disturbance_freq,REAL(iage-1,dp)) POPLUC%freq_age_secondary(g,iage) = 0.0 END DO POPLUC%freq_age_primary(g,:) = POPLUC%freq_age_primary(g,:) & / SUM(POPLUC%freq_age_primary(g,:))*fac ! Loop through secondary forest stands to transfer weights fac = 1.0 DO i = 1, POPLUC%n_event(g) age = POPLUC%age_history_secdf(g,i) POPLUC%freq_age_secondary(g,age+1) = POPLUC%area_history_secdf(g,i)/fac ENDDO END SUBROUTINE CALCULATE_WEIGHTS !============================================================================ SUBROUTINE INCREMENT_AGE(POPLUC,g) IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(INOUT) :: POPLUC INTEGER(i4b) :: n_event, i, j INTEGER(i4b), INTENT(IN) :: g REAL(dp):: area, remaining n_event = POPLUC%n_event(g) POPLUC%freq_age_secondary(g,2:age_max)=POPLUC%freq_age_secondary(g,1:age_max-1) POPLUC%freq_age_secondary(g,1) = 0.0 ! adjust secondary age distribution for secondary forest harvest area IF (POPLUC%smharv(g)+POPLUC%syharv(g).GT.0) THEN remaining = POPLUC%smharv(g)+POPLUC%syharv(g) i = age_max DO WHILE (remaining.GT.0.0) IF (POPLUC%freq_age_secondary(g,i) .GT. remaining) THEN POPLUC%freq_age_secondary(g,i) = POPLUC%freq_age_secondary(g,i) - remaining POPLUC%freq_age_secondary(g,1) = POPLUC%freq_age_secondary(g,1) + remaining remaining = 0.0 ELSE POPLUC%freq_age_secondary(g,1) = POPLUC%freq_age_secondary(g,1) + & POPLUC%freq_age_secondary(g,i) remaining = remaining - POPLUC%freq_age_secondary(g,i) POPLUC%freq_age_secondary(g,i) = 0.0 i = i-1; END IF IF (i.LT.2) THEN POPLUC%smharv(g) = POPLUC%smharv(g)+POPLUC%syharv(g) - remaining POPLUC%syharv(g) = 0.0 remaining = 0.0 ENDIF ENDDO ENDIF IF (IFHARVEST .AND. POPLUC%freq_age_secondary(g,ROTATION+1).GT.0.0) THEN POPLUC%freq_age_secondary(g,1) = POPLUC%freq_age_secondary(g,1) + & POPLUC%freq_age_secondary(g,ROTATION+1) POPLUC%freq_age_secondary(g,ROTATION+1) = 0.0 ENDIF ! adjust secondary age distribution for natural disturbance i = age_max DO i = age_max, 2 , -1 POPLUC%freq_age_secondary(g,1) = POPLUC%freq_age_secondary(g,1) + & POPLUC%freq_age_secondary(g,i)/disturbance_interval POPLUC%freq_age_secondary(g,i) = POPLUC%freq_age_secondary(g,i)* & (1. - 1./disturbance_interval) ENDDO END SUBROUTINE INCREMENT_AGE !============================================================================ SUBROUTINE POPLUCStep(POPLUC,year) IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(INOUT) :: POPLUC INTEGER(i4b), INTENT(IN) :: year INTEGER(i4b) :: g,j POPLUC%it = POPLUC%it + 1 IF (year.LT.POPLUC%thisyear) THEN DO g = 1,POPLUC%np ! CALL calculate_weights(POPLUC, g) CALL increment_age(POPLUC,g) ENDDO ELSE DO g = 1,POPLUC%np IF (POPLUC%ptos(g) .GT. 0.0) & CALL execute_luc_event('PRIMF','SECDF',POPLUC%ptos(g),g,POPLUC) IF (POPLUC%ptog(g) .GT. 0.0) & CALL execute_luc_event('PRIMF','C3ANN',POPLUC%ptog(g),g,POPLUC) IF (POPLUC%stop(g) .GT.0.0) & CALL execute_luc_event('SECDF','PRIMF',POPLUC%stop(g),g,POPLUC) IF (POPLUC%stog(g) .GT.0.0) & CALL execute_luc_event('SECDF','C3ANN',POPLUC%stog(g),g,POPLUC) IF (POPLUC%gtop(g) .GT.0.0) & CALL execute_luc_event('C3ANN','PRIMF',POPLUC%gtop(g),g,POPLUC) IF (POPLUC%gtos(g) .GT.0.0) & CALL execute_luc_event('C3ANN','SECDF',POPLUC%gtos(g),g,POPLUC) POPLUC%frac_forest(g) = POPLUC%primf(g)+ SUM(POPLUC%freq_age_secondary(g,:)) CALL increment_age(POPLUC,g) ENDDO ENDIF END SUBROUTINE POPLUCStep !============================================================================= SUBROUTINE POPLUC_weights_transfer(POPLUC,POP,LUC_EXPT) !------------------------------------------------------------------------------- ! This subroutine transfers LUC-determined age distributions to POP !------------------------------------------------------------------------------- IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(IN) :: POPLUC TYPE(POP_TYPE), INTENT(INOUT) :: POP TYPE (LUC_EXPT_TYPE), INTENT(IN) :: LUC_EXPT INTEGER:: g, k, j, l REAL(dp), DIMENSION(:), ALLOCATABLE:: freq_age DO l = 1, POP%np pop%pop_grid(l)%freq_age = 0.0 ENDDO ALLOCATE (freq_age(age_max)) DO g = 1,POPLUC%np ! loop over POPLUC gricells (same as CABLE grid-cells) IF (.NOT.LUC_EXPT%prim_only(g) .AND.SUM(POPLUC%freq_age_secondary(g,:)).GT.1e-12) THEN ! check if tile is secondary forest j = landpt(g)%cstart + 1 ! cable index for secondary forest tile freq_age = POPLUC%freq_age_secondary(g,:)/SUM(POPLUC%freq_age_secondary(g,:)) DO l = 1, POP%np ! for each wooded tile IF (POP%Iwood(l).EQ.j) THEN ! check if cable index in pop structure (Iwood) is the target cable index pop%pop_grid(l)%freq_age = freq_age ENDIF ENDDO ENDIF ENDDO DEALLOCATE (freq_age) END SUBROUTINE POPLUC_weights_transfer !============================================================================= SUBROUTINE POP_LUC_CASA_transfer(POPLUC,POP,LUC_EXPT,casapool,casabal,casaflux,ktauday) !------------------------------------------------------------------------------- ! This subroutine redestributes carbon (nitrogen and phosphorous) amongst ! pools according to land-use transtions !------------------------------------------------------------------------------- IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(INOUT) :: POPLUC TYPE(POP_TYPE), INTENT(IN) :: POP TYPE (LUC_EXPT_TYPE), INTENT(IN) :: LUC_EXPT TYPE (casa_pool), INTENT(INOUT) :: casapool TYPE (casa_balance), INTENT(INOUT) :: casabal TYPE (casa_flux), INTENT(INOUT) :: casaflux INTEGER, INTENT(IN) :: ktauday ! number of cable time-steps in a day (for needed for LUC flux output) INTEGER:: g, k, j, l, idp, irp, idlu, irlu, ilu INTEGER, PARAMETER :: & ptos = 1, & ptog = 2, & stog = 3, & gtos = 4 INTEGER, PARAMETER :: & p = 1, & s = 2, & gr = 3 REAL(dp) :: dcplant(nLU,nLU,3), dclitter(nLU,nLU,3), dcsoil(nLU,nLU,3) REAL(dp) :: dcplant_r(nLU,3), dclitter_r(nLU,3), dcsoil_r(nLU,3) REAL(dp) :: dcplant_d(nLU,3), dclitter_d(nLU,3), dcsoil_d(nLU,3) REAL(dp) :: dnplant(nLU,nLU,3), dnlitter(nLU,nLU,3), dnsoil(nLU,nLU,3) REAL(dp) :: dnplant_r(nLU,3), dnlitter_r(nLU,3), dnsoil_r(nLU,3) REAL(dp) :: dnplant_d(nLU,3), dnlitter_d(nLU,3), dnsoil_d(nLU,3) REAL(dp) :: dpplant(nLU,nLU,3), dplitter(nLU,nLU,3), dpsoil(nLU,nLU,3) REAL(dp) :: dpplant_r(nLU,3), dplitter_r(nLU,3), dpsoil_r(nLU,3) REAL(dp) :: dpplant_d(nLU,3), dplitter_d(nLU,3), dpsoil_d(nLU,3) REAL(dp) :: dnsoilmin_r(nLU), dclabile_r(nLU) REAL(dp) :: dnsoilmin_d(nLU), dclabile_d(nLU) REAL(dp) :: dclabile(nLU, nLU), dnsoilmin(nLU,nLU) REAL(dp) :: dA_r(nLU), dA_d(nLU), dA(nLU), deltaA, dwood_transfer REAL(dp), ALLOCATABLE :: dcHarvCLear(:), dcHarv(:), dcClear(:), dcExpand(:), FHarvClear(:) REAL(dp) :: kHarvProd(3), kClearProd(3) ! turnover rates for harvest and clearance products (y-1) kHarvProd(1) = 1.0/1.0 kHarvProd(2) = 1.0/10.0 kHarvProd(3) = 1.0/100.0 kClearProd = kHarvProd ! zero POPLUC fluxes popluc%FtransferNet = 0.0 popluc%FtransferGross = 0.0 popluc%FHarvest = 0.0 popluc%FClearance = 0.0 ! local variable for storing sum of biomass change due to !secondary harvest, clearance and expansion, and secondary forest ! harvest and clearance fluxes ALLOCATE(FHarvClear(POPLUC%np)) ALLOCATE(dcHarvClear(POPLUC%np)) ALLOCATE(dcHarv(POPLUC%np)) ALLOCATE(dcClear(POPLUC%np)) ALLOCATE(dcExpand(POPLUC%np)) FHarvClear = 0.0 dcHarvClear = 0.0 dcHarv = 0.0 dcClear = 0.0 dcExpand = 0.0 POPLUC%FHarvest = 0.0 POPLUC%FClearance = 0.0 casaflux%FluxCtohwp = 0.0 casaflux%FluxCtoclear = 0.0 casaflux%CtransferLUC = 0.0 ! Transfer POP age-dependent biomass to POLUC variables (diagnostic only) and ! catastrophic mortality in secondary forest tiles to changes in biomass associated ! with secondary forest harvest and clearance and expansion. DO g = 1,POPLUC%np ! loop over CABLE grid-cells j = landpt(g)%cstart ! index of primary veg tile in each grid-cell DO l = 1, POP%np ! loop over all POP tiles (== wooded tiles in CABLE) IF (.NOT.LUC_EXPT%prim_only(g)) THEN ! land-use change may occur IF (POP%Iwood(l).EQ.j+1 .AND. & (patch(j+1)%frac+ POPLUC%gtos(g)+POPLUC%ptos(g)-POPLUC%stog(g)) .GT. 0.0 ) THEN ! if secondary forest and new secondary forest area > 0 ! set POLUC diagnostic 2o forest age distribution to POP age distribution POPLUC%biomass_age_secondary(g,:)=POP%pop_grid(l)%biomass_age ! change in biomass area density due to secondary forest expansion [g C m-2] dcExpand(g) = -(POPLUC%gtos(g)+POPLUC%ptos(g))*casapool%cplant(j+1,2)/ & (patch(j+1)%frac + POPLUC%gtos(g)+POPLUC%ptos(g)-POPLUC%stog(g)) IF (POP%pop_grid(l)%cmass_sum_old .GT. 0.001) THEN ! change in biomass area density due to secondary harvest and clearance ! (== POP 'cataastrophic' mortality minus natural disturbance minus ! density reduction due to expansion dcHarvClear(g) = -MAX(POP%pop_grid(l)%cat_mortality/ & (POP%pop_grid(l)%cmass_sum_old + POP%pop_grid(l)%growth) & - 1.0/ disturbance_interval ,0.0) & *casapool%cplant(j+1,2) - dcExpand(g) IF (( POPLUC%smharv(g) + POPLUC%syharv(g) + POPLUC%stog(g)).GT.1e-10) THEN ! partition dcHarvClear into harvest and clearance components dcHarv(g) = dcHarvClear(g) * (POPLUC%smharv(g)+POPLUC%syharv(g)) & /( POPLUC%smharv(g) +POPLUC%syharv(g) + POPLUC%stog(g)) dcClear(g) = dcHarvClear(g) * POPLUC%stog(g) & /( POPLUC%smharv(g) + POPLUC%syharv(g) + POPLUC%stog(g)) IF ((casapool%cplant(j+1,2) + dcExpand(g)+dcClear(g)+dcHarv(g)).GT.0.0) THEN ! flux + A0C0 = (A + dA) * (C + dC) ! flux = A0 * dC + dA( C0 + dC) ! harvest+ clearance flux (not yet corrected for carbon remaining & ! in landscape as litter) ! [g C m-2] (grid-cell basis) FHarvClear(g) = -(patch(j+1)%frac * (dcHarvClear(g) + dcExpand(g)) + & (POPLUC%ptos(g)+POPLUC%gtos(g)-POPLUC%stog(g))* & ( casapool%cplant(j+1,2) + dcHarvClear(g) + dcExpand(g) ) ) ! secondary forest harvest flux (goes to harvest wood products) & ! [g C m-2] (grid-cell basis) ! corrected for carbon remaining in landscape as litter POPLUC%FHarvest(g,2) = (1.-POPLUC%fracHarvSecResid(g))*FHarvClear(g) & * (POPLUC%smharv(g) + POPLUC%syharv(g)) & /( POPLUC%smharv(g) + POPLUC%syharv(g) + POPLUC%stog(g)) ! secondary forest clearance flux (goes to clearance pool) & ! [g C m-2] (grid-cell basis) ! corrected for carbon remaining in landscape as litter POPLUC%FClearance(g,2) = (1.-POPLUC%fracClearResid(g))* FHarvClear(g) * & POPLUC%stog(g) & /( POPLUC%smharv(g) + POPLUC%syharv(g) + POPLUC%stog(g)) ELSE dcHarv(g) = 0.0 dcClear(g) = 0.0 POPLUC%FHarvest(g,2) = 0.0 POPLUC%FClearance(g,2) = 0.0 ENDIF ENDIF ENDIF ENDIF ENDIF IF (POP%Iwood(l).EQ.j ) THEN ! set diagnostic POPLUC variables to POP equivalents & !(patch biomass and age distributions) POPLUC%biomass_age_primary(g,:)=POP%pop_grid(l)%biomass_age POPLUC%freq_age_primary(g,:)=POP%pop_grid(l)%freq_age ENDIF ENDDO ENDDO ! Calculate Carbon Pool Transfers DO g = 1,POPLUC%np ! loop over POPLUC gridcells (== CABLE gridcells) j = landpt(g)%cstart ! start index of CABLE grid-cell tiles l = landpt(g)%cend ! end index of CABLE grid-cell tiles ! initialise local variables ! _r refers to receiver tile ! _d refers to donor tile ! dA refers to change in tile area ! dC refers to change in carbon [g m-2] (grid-cell basis) ! dn refers to change in nitrogen [g m-2] (grid-cell basis) ! dp refers to change in phosphorous [g m-2] (grid-cell basis) dA_r = 0.0 dA_d = 0.0 dA = 0.0 dcsoil = 0.0 dcplant = 0.0 dclitter = 0.0 dclabile = 0.0 dcsoil_r = 0.0 dcplant_r = 0.0 dclitter_r = 0.0 dcsoil_r = 0.0 dcplant_r = 0.0 dclitter_r = 0.0 dclabile_r = 0.0 dcsoil_d = 0.0 dcplant_d = 0.0 dclitter_d = 0.0 dcsoil_d = 0.0 dcplant_d = 0.0 dclitter_d = 0.0 dclabile_d = 0.0 dnsoil = 0.0 dnplant = 0.0 dnlitter = 0.0 dnsoil_r = 0.0 dnplant_r = 0.0 dnlitter_r = 0.0 dnsoil_r = 0.0 dnplant_r = 0.0 dnlitter_r = 0.0 dnsoilmin_r = 0.0 dnlitter = 0.0 dnsoil_d = 0.0 dnplant_d = 0.0 dnlitter_d = 0.0 dnsoil_d = 0.0 dnplant_d = 0.0 dnlitter_d = 0.0 dnsoilmin_d = 0.0 dpsoil = 0.0 dpplant = 0.0 dplitter = 0.0 dpsoil_r = 0.0 dpplant_r = 0.0 dplitter_r = 0.0 dpsoil_r = 0.0 dpplant_r = 0.0 dplitter_r = 0.0 dwood_transfer = 0.0 dpsoil_d = 0.0 dpplant_d = 0.0 dplitter_d = 0.0 dpsoil_d = 0.0 dpplant_d = 0.0 dplitter_d = 0.0 IF (.NOT.LUC_EXPT%prim_only(g)) THEN ! only worry about transitions where land-use change is possible DO k = 1,nTrans ! loop over all possible transition types IF (k==1) THEN deltaA = POPLUC%ptos(g) ! transition area: primary to seondary transition idp = j ! donor patch index irp = j+1 ! receiver patch index idlu = p ! donor land use index (p = primary) irlu = s ! receiver land use index (s = secondary) ELSEIF (k==2) THEN deltaA = POPLUC%ptog(g) ! transition area: primary to grass (open) transition idp = j ! donor patch index irp = j+2 ! receiver patch index idlu = p ! donor land use index (p = primary) irlu = gr ! receiver land use index (gr = grass or open) ELSEIF (k==3) THEN deltaA = POPLUC%stog(g) ! transition area: secondary to grass (open) transition idp = j+1 ! donor patch index irp = j+2 ! receiver patch index idlu = s ! donor land use index irlu = gr ! receiver land use index ELSEIF(k==4) THEN deltaA = POPLUC%gtos(g) ! transition area: grass (open) to secondary transition idp = j+2 ! donor patch index irp = j+1 ! receiver patch index idlu = gr ! donor land use index irlu = s ! receiver land use index (s = secondary) ENDIF dcsoil = 0 dclitter = 0 dcplant = 0 dclabile = 0 dnsoil = 0 dnlitter = 0 dnplant = 0 dpsoil = 0 dplitter = 0 dpplant = 0 dwood_transfer = 0.0 ! transfer fluxes : only consider cases where gross transition area is greater than zero IF (deltaA.GT.0.0) THEN ! all soil pools ! change in carbon associated with gross transition dcsoil(irlu,idlu,:) = deltaA*casapool%csoil(idp,:) ! change in receiver carbon pool (accumulated over all gross transitions) dcsoil_r(irlu,:) = dcsoil_r(irlu,:) + deltaA*casapool%csoil(idp,:) ! change in donor carbon pool (accumulated over all gross transitions) dcsoil_d(idlu,:) = dcsoil_d(idlu,:) - deltaA*casapool%csoil(idp,:) dnsoil(irlu,idlu,:) = deltaA*casapool%nsoil(idp,:) dnsoil_r(irlu,:) = dnsoil_r(irlu,:) + deltaA*casapool%nsoil(idp,:) dnsoil_d(idlu,:) = dnsoil_d(idlu,:) - deltaA*casapool%nsoil(idp,:) dnsoilmin_r(irlu) = dnsoilmin_r(irlu) + deltaA*casapool%nsoilmin(idp) dnsoilmin_d(irlu) = dnsoilmin_d(idlu) - deltaA*casapool%nsoilmin(idp) dpsoil(irlu,idlu,:) = deltaA*casapool%psoil(idp,:) dpsoil_r(irlu,:) = dpsoil_r(irlu,:) + deltaA*casapool%psoil(idp,:) dpsoil_d(idlu,:) = dpsoil_d(idlu,:) - deltaA*casapool%psoil(idp,:) ! need to inlcude other P pools here too: occluded P and labile P) ! microbial litter dclitter(irlu,idlu,1) = deltaA*casapool%clitter(idp,1) dclitter_r(irlu,1) = dclitter_r(irlu,1) + dclitter(irlu,idlu,1) dclitter_d(idlu,1) = dclitter_d(idlu,1) - dclitter(irlu,idlu,1) dnlitter(irlu,idlu,1) = deltaA*casapool%nlitter(idp,1) dnlitter_r(irlu,1) = dnlitter_r(irlu,1) + dnlitter(irlu,idlu,1) dnlitter_d(idlu,1) = dnlitter_d(idlu,1) - dnlitter(irlu,idlu,1) dplitter(irlu,idlu,1) = deltaA*casapool%plitter(idp,1) dplitter_r(irlu,1) = dplitter_r(irlu,1) + dplitter(irlu,idlu,1) dplitter_d(idlu,1) = dplitter_d(idlu,1) - dplitter(irlu,idlu,1) ! CWD : donor pool inhreits CWD and residues from harvest/clearance IF (idlu == s .AND. casapool%cplant(idp,2).GT.1e-5 ) THEN ! seocndary forest clearance: use secondary forest clearance flux from above dclitter(irlu,idlu,3) = deltaA*casapool%clitter(idp,3) + & POPLUC%fracClearResid(g)*POPLUC%Fclearance(g,2)/ & (1.-POPLUC%fracClearResid(g)) dnlitter(irlu,idlu,3) = deltaA*casapool%nlitter(idp,3) + & POPLUC%fracClearResid(g)*POPLUC%Fclearance(g,2)/ & (1.-POPLUC%fracClearResid(g)) & *casapool%nplant(idp,2)/ casapool%cplant(idp,2) dplitter(irlu,idlu,3) = deltaA*casapool%plitter(idp,3) + & POPLUC%fracClearResid(g)*POPLUC%Fclearance(g,2)/ & (1.-POPLUC%fracClearResid(g)) & *casapool%pplant(idp,2)/ casapool%cplant(idp,2) ELSEIF (idlu == p .AND. irlu == s) THEN ! primary forest harvest: assume harvest is uniform across age-classes dclitter(irlu,idlu,3) = deltaA*(casapool%clitter(idp,3) + & POPLUC%fracHarvResid(g)*casapool%cplant(idp,2) ) dnlitter(irlu,idlu,3) = deltaA*(casapool%nlitter(idp,3) + & POPLUC%fracHarvResid(g)*casapool%nplant(idp,2) ) dplitter(irlu,idlu,3) = deltaA*(casapool%plitter(idp,3) + & POPLUC%fracHarvResid(g) *casapool%pplant(idp,2) ) ELSEIF (idlu == p .AND. irlu == gr) THEN ! primary forest clearance: assume clearance is uniform across age-classes dclitter(irlu,idlu,3) = deltaA*(casapool%clitter(idp,3) + & POPLUC%fracClearResid(g)*casapool%cplant(idp,2) ) dnlitter(irlu,idlu,3) = deltaA*(casapool%nlitter(idp,3) + & POPLUC%fracClearResid(g)*casapool%nplant(idp,2) ) dplitter(irlu,idlu,3) = deltaA*(casapool%plitter(idp,3) + & POPLUC%fracClearResid(g) *casapool%pplant(idp,2) ) ELSE dclitter(irlu,idlu,3) = deltaA*casapool%clitter(idp,3) dnlitter(irlu,idlu,3) = deltaA*casapool%nlitter(idp,3) dplitter(irlu,idlu,3) = deltaA*casapool%plitter(idp,3) ENDIF dclitter_r(irlu,3) = dclitter_r(irlu,3) + dclitter(irlu,idlu,3) dclitter_d(idlu,3) = dclitter_d(idlu,3) - dclitter(irlu,idlu,3) dnlitter_r(irlu,3) = dnlitter_r(irlu,3) + dnlitter(irlu,idlu,3) dnlitter_d(idlu,3) = dnlitter_d(idlu,3) - dnlitter(irlu,idlu,3) dplitter_r(irlu,3) = dplitter_r(irlu,3) + dplitter(irlu,idlu,3) dplitter_d(idlu,3) = dplitter_d(idlu,3) - dplitter(irlu,idlu,3) ! fine structural litter: donor pool inherits leaves and fine roots dclitter(irlu,idlu,2) = deltaA*(casapool%clitter(idp,2) + casapool%cplant(idp,1) + & casapool%cplant(idp,3)) dclitter_r(irlu,2) = dclitter_r(irlu,2) + dclitter(irlu,idlu,2) dclitter_d(idlu,2) = dclitter_d(idlu,2) - dclitter(irlu,idlu,2) dnlitter(irlu,idlu,2) = deltaA*(casapool%nlitter(idp,2) + casapool%cplant(idp,1) + & casapool%cplant(idp,3)) dnlitter_r(irlu,2) = dnlitter_r(irlu,2) + dnlitter(irlu,idlu,2) dnlitter_d(idlu,2) = dnlitter_d(idlu,2) - dnlitter(irlu,idlu,2) dplitter(irlu,idlu,2) = deltaA*(casapool%plitter(idp,2) + casapool%cplant(idp,1) + & casapool%cplant(idp,3)) dplitter_r(irlu,2) = dplitter_r(irlu,2) + dplitter(irlu,idlu,2) dplitter_d(idlu,2) = dplitter_d(idlu,2) - dplitter(irlu,idlu,2) ! labile carbon dclabile(irlu,idlu) = deltaA*casapool%clabile(idp) dclabile_r(irlu) = dclabile_r(irlu) + deltaA*casapool%clabile(idp) dclabile_d(irlu) = dclabile_d(idlu) - deltaA*casapool%clabile(idp) ! biomass: no biomass inherited dcplant(irlu,idlu,:) = 0.0 dcplant_r(irlu,:) = dcplant_r(irlu,:) + dcplant(irlu,idlu,:) dcplant_d(idlu,:) = dcplant_d(idlu,:) - dcplant(irlu,idlu,:) dnplant(irlu,idlu,:) = 0.0 dnplant_r(irlu,:) = dnplant_r(irlu,:) + dnplant(irlu,idlu,:) dnplant_d(idlu,:) = dnplant_d(idlu,:) - dnplant(irlu,idlu,:) dpplant(irlu,idlu,:) = 0.0 dpplant_r(irlu,:) = dpplant_r(irlu,:) + dpplant(irlu,idlu,:) dpplant_d(idlu,:) = dpplant_d(idlu,:) - dpplant(irlu,idlu,:) ! Gross Transfer Flux (total C-tranfer associated with k-th tranistion) popluc%FtransferGross(g,k) = & SUM(dcsoil(irlu,idlu,:) + dclitter(irlu,idlu,:) + dcplant(irlu,idlu,:)) ! Harvest Flux ! primary to secondary forest ! (note secondary harvest and clearance fluxes are already evaluated & !at the top of this subroutine) IF (idlu==p .AND. irlu ==s) THEN popluc%FHarvest(g,idlu) = (1.0 -POPLUC%fracHarvResid(g)) & *casapool%cplant(idp,2)*deltaA ENDIF ! Clearance Flux ! primary forest to grass IF ((idlu==p) .AND. irlu ==gr) THEN popluc%FClearance(g,idlu) = (1.0 -POPLUC%fracClearResid(g)) & * casapool%cplant(idp,2)*deltaA ENDIF ! transition area dA_r(irlu) = dA_r(irlu) + deltaA dA_d(idlu) = dA_d(idlu) - deltaA ENDIF ENDDO ! ntrans DO ilu=1,nLU ! update pools irp = ilu + j -1 dwood_transfer = 0.0 dA(ilu) = dA_r(ilu) + dA_d(ilu) ! Net Transfer Flux POPLUC%FTransferNet(g,ilu) = SUM(dcsoil_r(ilu,:) + dcsoil_d(ilu,:)) + & SUM(dclitter_r(ilu,:) + dclitter_d(ilu,:)) + & SUM(dcplant_r(ilu,:) + dcplant_d(ilu,:)) + & dclabile_r(ilu) + dclabile_d(ilu) IF (ilu ==2) THEN ! augment CWD pools in secondary veg tiles by harvest residues dclitter_r(ilu,3) = dclitter_r(irlu,3) + & POPLUC%fracHarvSecResid(g)*POPLUC%FHarvest(g,2)/(1.0 -POPLUC%fracHarvSecResid(g)) IF (casapool%cplant(irp,2) .GT. 1.e-5) THEN dnlitter_r(ilu,3) = dnlitter_r(irlu,3) + & POPLUC%fracHarvSecResid(g)*POPLUC%FHarvest(g,2)/(1.0 -POPLUC%fracHarvSecResid(g)) & *casapool%nplant(irp,2)/ casapool%cplant(irp,2) dplitter_r(ilu,3) = dplitter_r(irlu,3) + & POPLUC%fracHarvSecResid(g)*POPLUC%FHarvest(g,2)/(1.0 -POPLUC%fracHarvSecResid(g)) & *casapool%pplant(irp,2)/ casapool%cplant(irp,2) ENDIF ENDIF IF ((patch(irp)%frac+dA(ilu)).GT.1.e-6 ) THEN ! avoid fpe's by ensuring finite & ! new tile area casapool%nsoilmin(irp) = casapool%nsoilmin(irp) + & (dnsoilmin_r(ilu) - casapool%nsoilmin(irp)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%clabile(irp) = casapool%clabile(irp) + & (dclabile_r(ilu) - casapool%clabile(irp)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casaflux%CtransferLUC(irp) = casaflux%CtransferLUC(irp)+ & (dclabile_r(ilu) - casapool%clabile(irp)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%csoil(irp,:) = casapool%csoil(irp,:) + & (dcsoil_r(ilu,:)-casapool%csoil(irp,:)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casaflux%CtransferLUC(irp) = casaflux%CtransferLUC(irp)+ & SUM((dcsoil_r(ilu,:)-casapool%csoil(irp,:)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu))) casapool%nsoil(irp,:) = casapool%nsoil(irp,:) + & (dnsoil_r(ilu,:)-casapool%nsoil(irp,:)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%psoil(irp,:) = casapool%psoil(irp,:) + & (dpsoil_r(ilu,:)-casapool%psoil(irp,:)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%clitter(irp,:) = casapool%clitter(irp,:) + & (dclitter_r(ilu,:)-casapool%clitter(irp,:)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casaflux%CtransferLUC(irp) = casaflux%CtransferLUC(irp)+ & SUM( (dclitter_r(ilu,:)-casapool%clitter(irp,:)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu))) casapool%nlitter(irp,:) = casapool%nlitter(irp,:) + & (dnlitter_r(ilu,:)-casapool%nlitter(irp,:)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%plitter(irp,:) = casapool%plitter(irp,:) + & (dplitter_r(ilu,:)-casapool%plitter(irp,:)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%cplant(irp,1) = casapool%cplant(irp,1) + & (dcplant_r(ilu,1)-casapool%cplant(irp,1)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casaflux%CtransferLUC(irp) = casaflux%CtransferLUC(irp)+ & (dcplant_r(ilu,1)-casapool%cplant(irp,1)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%nplant(irp,1) = casapool%nplant(irp,1) + & (dnplant_r(ilu,1)-casapool%nplant(irp,1)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%pplant(irp,1) = casapool%pplant(irp,1) + & (dpplant_r(ilu,1)-casapool%pplant(irp,1)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%cplant(irp,3) = casapool%cplant(irp,3) + & (dcplant_r(ilu,3)-casapool%cplant(irp,3)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casaflux%CtransferLUC(irp) = casaflux%CtransferLUC(irp)+ & (dcplant_r(ilu,3)-casapool%cplant(irp,3)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%nplant(irp,3) = casapool%nplant(irp,3) + & (dnplant_r(ilu,3)-casapool%nplant(irp,3)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%pplant(irp,3) = casapool%pplant(irp,3) + & (dpplant_r(ilu,3)-casapool%pplant(irp,3)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%cplant(irp,2) = casapool%cplant(irp,2) + & (dcplant_r(ilu,2)-casapool%cplant(irp,2)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casaflux%CtransferLUC(irp) = casaflux%CtransferLUC(irp)+ & (dcplant_r(ilu,2)-casapool%cplant(irp,2)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%nplant(irp,2) = casapool%nplant(irp,2) + & (dnplant_r(ilu,2)-casapool%nplant(irp,2)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) casapool%pplant(irp,2) = casapool%pplant(irp,2) + & (dpplant_r(ilu,2)-casapool%pplant(irp,2)*(dA(ilu) - dA_d(ilu))) & /(patch(irp)%frac+dA(ilu)) ! account here for change in secondary forest biomass density due to: ! harvest and clearing, as well as increase in (below ground) CWD ! where secondary forest harvest occurs IF (ilu .EQ.s .AND. (casapool%cplant(irp,2) +(dcHarv(g)+dCClear(g)) ).GT.0.0 & .AND. casapool%cplant(irp,2).GT.1.e-10 ) THEN casapool%nplant(irp,2) = casapool%nplant(irp,2) + (dcHarv(g)+dcClear(g)) & * casapool%nplant(irp,2)/casapool%cplant(irp,2) casapool%pplant(irp,2) = casapool%pplant(irp,2) + (dcHarv(g)+dcClear(g)) & * casapool%pplant(irp,2)/casapool%cplant(irp,2) casapool%cplant(irp,2) = casapool%cplant(irp,2) + (dcHarv(g)+dcClear(g)) ELSEIF (ilu .EQ.s .AND. (casapool%cplant(irp,2) +(dcHarv(g)+dCClear(g)) ).LE.0.0 ) THEN POPLUC%FHarvest(g,2) = 0.0 POPLUC%FClearance(g,2) = 0.0 ENDIF ENDIF IF (patch(irp)%frac .GT. 1e-8) THEN casaflux%FluxCtohwp(irp) = POPLUC%FHarvest(g,ilu)/ patch(irp)%frac * ktauday casaflux%FluxCtoclear(irp) = POPLUC%FClearance(g,ilu)/ patch(irp)%frac * ktauday ENDIF ENDDO ! POPLUC diagnostics ! pools in gC per m2 of gridcell ! NEP in g C y-1 per m2 of gridcell POPLUC%FNEP(g,:) = casabal%Fcneeyear(j:l)*patch(j:l)%frac ! note NEE = NEP here !(+ve into surface) DO ilu=1,nLU ! update area weights irp = ilu + j -1 ! Update tile area IF (ilu == p) THEN POPLUC%primf(g) = MAX(patch(irp)%frac + dA_r(ilu) + dA_d(ilu), 0.0) ELSEIF (ilu == s) THEN POPLUC%secdf(g) = MAX(patch(irp)%frac + dA_r(ilu) + dA_d(ilu), 0.0) ELSEIF (ilu == gr) THEN POPLUC%grass(g) = MAX(patch(irp)%frac + dA_r(ilu) + dA_d(ilu), 0.0) ENDIF POPLUC%csoil(g,ilu) = SUM(casapool%csoil(irp,:))* & MAX(patch(irp)%frac + dA_r(ilu) + dA_d(ilu), 0.0) POPLUC%clitt(g,ilu) = SUM(casapool%clitter(irp,:))* & MAX(patch(irp)%frac + dA_r(ilu) + dA_d(ilu), 0.0) POPLUC%cbiomass(g,ilu) = SUM(casapool%cplant(irp,:))* & MAX(patch(irp)%frac + dA_r(ilu) + dA_d(ilu), 0.0) ENDDO ELSE POPLUC%csoil(g,:) = 0.0 POPLUC%clitt(g,:) = 0.0 POPLUC%cbiomass(g,:) = 0.0 POPLUC%FNEP(g,:) = 0.0 POPLUC%csoil(g,1) = SUM(casapool%csoil(j,:))*patch(j)%frac POPLUC%clitt(g,1) = SUM(casapool%clitter(j,:))*patch(j)%frac POPLUC%cbiomass(g,1) = SUM(casapool%cplant(j,:))*patch(j)%frac POPLUC%FNEP(g,1) = casabal%Fcneeyear(j)*patch(j)%frac POPLUC%primf(g) = patch(j)%frac POPLUC%secdf(g) = 0.0 IF (POPLUC%grass(g).GT.0.0 .AND. l.EQ.j+1) THEN POPLUC%csoil(g,3) = SUM(casapool%csoil(l,:))*patch(l)%frac POPLUC%clitt(g,3) = SUM(casapool%clitter(l,:))*patch(l)%frac POPLUC%cbiomass(g,3) = SUM(casapool%cplant(l,:))*patch(l)%frac POPLUC%FNEP(g,3) = casabal%Fcneeyear(l)*patch(l)%frac POPLUC%grass(g) = patch(l)%frac ENDIF ENDIF POPLUC%HarvProdLoss(g,:) = kHarvProd * POPLUC%HarvProd(g,:) POPLUC%ClearProdLoss(g,:) = kClearProd * POPLUC%ClearProd(g,:) DO j=1,3 POPLUC%HarvProd(g,j) = POPLUC%HarvProd(g,j) + & POPLUC%fracHarvProd(g,j)*SUM(POPLUC%FHarvest(g,:)) - POPLUC%HarvProdLoss(g,j) POPLUC%ClearProd(g,j) = POPLUC%ClearProd(g,j) + & POPLUC%fracClearProd(g,j)*SUM(POPLUC%FClearance(g,:)) - POPLUC%ClearProdLoss(g,j) ENDDO ENDDO 991 FORMAT(1166(e14.7,2x)) ! update total carbon pools and "last" pool values for use in carbon balance checks. casapool%ctot = SUM(casapool%cplant,2)+SUM(casapool%clitter,2)+ & SUM(casapool%csoil,2)+casapool%clabile casabal%cplantlast = casapool%cplant casabal%clabilelast = casapool%clabile casabal%clitterlast = casapool%clitter casabal%csoillast = casapool%csoil END SUBROUTINE POP_LUC_CASA_transfer !============================================================================ SUBROUTINE POPLUC_Init(POPLUC,LUC_EXPT, casapool, casaflux, casabiome, veg, POP, np) USE cable_def_types_mod, ONLY : veg_parameter_type USE casaparm, ONLY: LEAF, WOOD, FROOT IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(INOUT) :: POPLUC TYPE (LUC_EXPT_TYPE), INTENT(IN) :: LUC_EXPT TYPE (casa_pool), INTENT(INOUT) :: casapool TYPE (casa_flux), INTENT(INOUT) :: casaflux TYPE (casa_biome), INTENT(IN) :: casabiome TYPE (veg_parameter_type), INTENT(IN) :: veg ! vegetation parameters TYPE (POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: np INTEGER(i4b) :: g, j, k, l CALL alloc_POPLUC(POPLUC,np) POPLUC%it = 0 POPLUC%np = np CALL ZeroPOPLUC(POPLUC) CALL POPLUC_set_params(POPLUC, LUC_EXPT) IF (cable_user%POPLUC_RunType .EQ. 'init') THEN POPLUC%frac_primf = LUC_EXPT%primaryf POPLUC%primf = LUC_EXPT%primaryf POPLUC%grass = LUC_EXPT%grass WHERE ((POPLUC%primf + POPLUC%grass) > 1.0) & POPLUC%grass = 1.0 - POPLUC%primf POPLUC%frac_forest = 1- POPLUC%grass POPLUC%freq_age_secondary(:,1) = MAX(POPLUC%frac_forest - POPLUC%primf, 0.0) POPLUC%latitude = patch(landpt(:)%cstart)%latitude POPLUC%longitude = patch(landpt(:)%cstart)%longitude ! zero biomass in secondary forest tiles (both CASA and POP variables) WHERE (veg%iLU ==2) casapool%cplant(:,leaf) = 0.01 casapool%nplant(:,leaf)= casabiome%ratioNCplantmin(veg%iveg,leaf)* casapool%cplant(:,leaf) casapool%pplant(:,leaf)= casabiome%ratioPCplantmin(veg%iveg,leaf)* casapool%cplant(:,leaf) casapool%cplant(:,froot) = 0.01 casapool%nplant(:,froot)= casabiome%ratioNCplantmin(veg%iveg,froot)* & casapool%cplant(:,froot) casapool%pplant(:,froot)= casabiome%ratioPCplantmin(veg%iveg,froot)* & casapool%cplant(:,froot) casapool%cplant(:,wood) = 0.01 casapool%nplant(:,wood)= casabiome%ratioNCplantmin(veg%iveg,wood)* casapool%cplant(:,wood) casapool%pplant(:,wood)= casabiome%ratioPCplantmin(veg%iveg,wood)* casapool%cplant(:,wood) casaflux%frac_sapwood = 1.0 endwhere DO k=1,np IF (.NOT.LUC_EXPT%prim_only(k)) THEN j = landpt(k)%cstart+1 DO l=1,SIZE(POP%Iwood) IF( POP%Iwood(l) == j) THEN CALL POP_init_single(POP,veg%disturbance_interval,l) EXIT ENDIF ENDDO ENDIF ENDDO ELSEIF (cable_user%POPLUC_RunType .EQ. 'restart') THEN CALL READ_LUC_RESTART_NC (POPLUC) POPLUC%frac_primf = POPLUC%primf POPLUC%frac_forest = 1- POPLUC%grass POPLUC%latitude = patch(landpt(:)%cstart)%latitude POPLUC%longitude = patch(landpt(:)%cstart)%longitude ENDIF ! set landuse index for secondary forest POP landscapes ! (not for 'static') run type because secondary forest dynamics are only ! simulated with dynamic land-use forcing IF (TRIM(cable_user%POPLUC_RunType) .NE. 'static') THEN DO k=1,POP%np IF (veg%iLU(POP%Iwood(k)).EQ.2) THEN POP%pop_grid(k)%LU = 2 ENDIF ENDDO ENDIF END SUBROUTINE POPLUC_Init !============================================================================ SUBROUTINE POPLUC_set_patchfrac(POPLUC,LUC_EXPT) IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(IN) :: POPLUC TYPE (LUC_EXPT_TYPE), INTENT(IN) :: LUC_EXPT INTEGER(i4b) :: j, k, l DO k=1,POPLUC%np j = landpt(k)%cstart l = landpt(k)%cend IF (.NOT.LUC_EXPT%prim_only(k)) THEN patch(j)%frac = POPLUC%primf(k) patch(l)%frac = POPLUC%grass(k) patch(j+1)%frac = 1.0 - patch(j)%frac - patch(l)%frac ENDIF ENDDO END SUBROUTINE POPLUC_SET_PATCHFRAC !============================================================================ SUBROUTINE POPLUC_set_params(POPLUC,LUC_EXPT) IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(INOUT) :: POPLUC TYPE (LUC_EXPT_TYPE), INTENT(IN) :: LUC_EXPT INTEGER(i4b) :: g, np np = POPLUC%np DO g = 1,np POPLUC%fracharvProd(g, 1) = 0.9 POPLUC%fracharvProd(g, 2) = 0.04 POPLUC%fracharvProd(g, 3) = 0.06 POPLUC%fracClearProd(g, 1) = 0.597 POPLUC%fracClearProd(g, 2) = 0.403 POPLUC%fracClearProd(g, 3) = 0.00 POPLUC%fracClearResid(g) = 0.33 POPLUC%fracHarvResid(g) = 0.79 POPLUC%fracHarvSecResid(g) = 0.81 IF (LUC_EXPT%biome(g)==1 .OR. LUC_EXPT%biome(g)==2) THEN ! Tropical Evergreen and Tropical Deciduous POPLUC%fracharvProd(g, 1) = 0.9 POPLUC%fracharvProd(g, 2) = 0.04 POPLUC%fracharvProd(g, 3) = 0.06 POPLUC%fracClearProd(g, 1) = 0.597 POPLUC%fracClearProd(g, 2) = 0.403 POPLUC%fracClearProd(g, 3) = 0.00 IF (LUC_EXPT%biome(g)==1) POPLUC%fracHarvResid(g) = 0.79 IF (LUC_EXPT%biome(g)==2) POPLUC%fracHarvResid(g) = 0.86 IF (LUC_EXPT%biome(g)==1) POPLUC%fracHarvSecResid(g) = 0.71 IF (LUC_EXPT%biome(g)==2) POPLUC%fracHarvSecResid(g) = 0.81 POPLUC%fracClearResid(g) = 0.33 ELSEIF (LUC_EXPT%biome(g).GE.4 .OR. LUC_EXPT%biome(g).LE.10) THEN ! Other Forest POPLUC%fracharvProd(g, 1) = 0.4 POPLUC%fracharvProd(g, 2) = 0.24 POPLUC%fracharvProd(g, 3) = 0.36 POPLUC%fracClearProd(g, 1) = 0.597 POPLUC%fracClearProd(g, 2) = 0.2985 POPLUC%fracClearProd(g, 3) = 0.1045 IF (LUC_EXPT%ivegp(g)==2) POPLUC%fracHarvResid(g) = 0.83 IF (LUC_EXPT%ivegp(g)==1 .OR. LUC_EXPT%ivegp(g)==3 ) POPLUC%fracHarvResid(g) = 0.87 IF (LUC_EXPT%ivegp(g)==4) POPLUC%fracHarvResid(g) = 0.78 IF (LUC_EXPT%ivegp(g)==2) POPLUC%fracHarvSecResid(g) = 0.75 IF (LUC_EXPT%ivegp(g)==1 .OR. LUC_EXPT%ivegp(g)==3 ) POPLUC%fracHarvSecResid(g) = 0.82 IF (LUC_EXPT%ivegp(g)==4) POPLUC%fracHarvSecResid(g) = 0.70 POPLUC%fracClearResid(g) = 0.33 ELSEIF (LUC_EXPT%biome(g).EQ.3 .OR. LUC_EXPT%biome(g).GE.11) THEN ! savanna and shrub POPLUC%fracharvProd(g, 1) = 1.0 POPLUC%fracharvProd(g, 2) = 0.00 POPLUC%fracharvProd(g, 3) = 0.00 POPLUC%fracClearProd(g, 1) = 0.8 POPLUC%fracClearProd(g, 2) = 0.2 POPLUC%fracClearProd(g, 3) = 0.0 IF (LUC_EXPT%biome(g).EQ.13 .OR.(LUC_EXPT%biome(g).EQ.14) ) THEN POPLUC%fracHarvResid(g) = 0.78 POPLUC%fracHarvSecResid(g) = 0.70 ELSE POPLUC%fracHarvResid(g) = 0.86 POPLUC%fracHarvSecResid(g) = 0.81 ENDIF POPLUC%fracClearResid(g) = 0.5 ENDIF ENDDO END SUBROUTINE POPLUC_SET_PARAMS !============================================================================= SUBROUTINE alloc_POPLUC(POPLUC, arraysize) IMPLICIT NONE TYPE(POPLUC_TYPE),INTENT(INOUT) :: POPLUC INTEGER, INTENT(IN) :: arraysize ALLOCATE(POPLUC%it) ALLOCATE(POPLUC%np) ALLOCATE(POPLUC%firstyear) ALLOCATE(POPLUC%thisyear) ALLOCATE(POPLUC%latitude(arraysize)) ALLOCATE(POPLUC%longitude(arraysize)) ALLOCATE(POPLUC%n_event(arraysize)) ALLOCATE(POPLUC%primf(arraysize)) ALLOCATE(POPLUC%secdf(arraysize)) ALLOCATE(POPLUC%grass(arraysize)) ALLOCATE(POPLUC%ptos(arraysize)) ALLOCATE(POPLUC%ptog(arraysize)) ALLOCATE(POPLUC%stop(arraysize)) ALLOCATE(POPLUC%stog(arraysize)) ALLOCATE(POPLUC%gtop(arraysize)) ALLOCATE(POPLUC%gtos(arraysize)) ALLOCATE(POPLUC%ptog(arraysize)) ALLOCATE(POPLUC%pharv(arraysize)) ALLOCATE(POPLUC%smharv(arraysize)) ALLOCATE(POPLUC%syharv(arraysize)) ALLOCATE(POPLUC%frac_primf(arraysize)) ALLOCATE(POPLUC%frac_forest(arraysize)) ALLOCATE(POPLUC%freq_age_primary(arraysize,AGE_MAX)) ALLOCATE(POPLUC%freq_age_secondary(arraysize,AGE_MAX)) ALLOCATE(POPLUC%biomass_age_primary(arraysize,AGE_MAX)) ALLOCATE(POPLUC%biomass_age_secondary(arraysize,AGE_MAX)) ALLOCATE(POPLUC%age_history_secdf(arraysize,LENGTH_SECDF_HISTORY)) ALLOCATE(POPLUC%area_history_secdf(arraysize,LENGTH_SECDF_HISTORY)) ALLOCATE(POPLUC%FNEP(arraysize,nLU)) ALLOCATE(POPLUC%Clitt(arraysize,nLU)) ALLOCATE(POPLUC%Csoil(arraysize,nLU)) ALLOCATE(POPLUC%Cbiomass(arraysize,nLU)) ALLOCATE(POPLUC%FHarvest(arraysize,nLU)) ALLOCATE(POPLUC%FClearance(arraysize,nLU)) ALLOCATE(POPLUC%FTransferNet(arraysize,nLU)) ALLOCATE(POPLUC%FTransferGross(arraysize,nTrans)) ALLOCATE(POPLUC%HarvProd(arraysize,3)) ALLOCATE(POPLUC%ClearProd(arraysize,3)) ALLOCATE(POPLUC%HarvProdLoss(arraysize,3)) ALLOCATE(POPLUC%ClearProdLoss(arraysize,3)) ALLOCATE(POPLUC%fracHarvProd(arraysize,3)) ALLOCATE(POPLUC%fracClearProd(arraysize,3)) ALLOCATE(POPLUC%fracHarvResid(arraysize)) ALLOCATE(POPLUC%fracHarvSecResid(arraysize)) ALLOCATE(POPLUC%fracClearResid(arraysize)) END SUBROUTINE alloc_POPLUC !============================================================================= ! Exponential distribution ! Returns probability of a given time-between-events (x) ! Given a Poisson process with expected frequency (events per unit time) lambda ! Reference: http://en.wikipedia.org/wiki/Exponential_distribution ! Use to determine average age (x, years) of patches with a given random disturbance ! frequency lambda (disturbances per year) REAL(dp) FUNCTION REALExponential(lambda, x) IMPLICIT NONE REAL(dp), INTENT(IN) :: x REAL(dp), INTENT(IN) :: lambda IF (x.LT.0) THEN ! Shouldn't happen but ... REALExponential=0.0 ELSE REALExponential=lambda*EXP(-lambda*x) ENDIF END FUNCTION REALExponential !============================================================================= SUBROUTINE WRITE_LUC_OUTPUT_NC ( POPLUC, ctime, FINAL ) USE CABLE_COMMON_MODULE, ONLY: filename, cable_user USE casa_ncdf_module, ONLY: HANDLE_ERR USE netcdf IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(IN) :: POPLUC LOGICAL, INTENT(IN) :: FINAL INTEGER, INTENT(IN) :: ctime INTEGER :: STATUS INTEGER :: land_ID, age_ID, hist_ID, t_ID, nLU_ID, nTrans_ID INTEGER :: i, mp, nprod, nprod_ID LOGICAL :: CASAONLY CHARACTER :: CYEAR*4, FNAME*99,dum*50 LOGICAL, SAVE :: CALL1 = .TRUE. ! 1 dim arrays (mp ) CHARACTER(len=20),DIMENSION(2) :: A0 ! 2 dim real arrays (mp,t) CHARACTER(len=20),DIMENSION(13):: A1 ! 2 dim integer arrays (mp,t) CHARACTER(len=20),DIMENSION(1):: AI1 ! 3 dim real arrays (mp,age_max,t) CHARACTER(len=25),DIMENSION(4) :: A2 ! 3 dim real arrays (mp,LENGTH_SECDF_HISTORY,t) CHARACTER(len=20),DIMENSION(1) :: A3 ! 3 dim integer arrays (mp,LENGTH_SECDF_HISTORY,t) CHARACTER(len=20),DIMENSION(1) :: AI3 ! 3 dim real arrays (mp,nLU,t) CHARACTER(len=20),DIMENSION(7) :: A4 ! 3 dim real arrays (mp,nTrans,t) CHARACTER(len=20),DIMENSION(1) :: A5 ! 3 dim real arrays (mp,nprod,t) CHARACTER(len=20),DIMENSION(4) :: A6 INTEGER, SAVE :: VIDtime, VID0(SIZE(A0)),VID1(SIZE(A1)),VIDI1(SIZE(AI1)) INTEGER, SAVE :: VID2(SIZE(A2)),VID3(SIZE(A3)),VIDI3(SIZE(AI3)) INTEGER, SAVE :: VID4(SIZE(A4)),VID5(SIZE(A5)), VID6(SIZE(A6)) INTEGER, SAVE :: FILE_ID, CNT = 0 CHARACTER(len=50) :: RecordDimName REAL, ALLOCATABLE :: freq_age_secondary(:,:) INTEGER :: g LOGICAL :: put_age_vars mp = POPLUC%np nprod = 3 put_age_vars=.FALSE. ALLOCATE(freq_age_secondary(mp,age_max)) A0(1) = 'latitude' A0(2) = 'longitude' A1(1) = 'primf' A1(2) = 'secdf' A1(3) = 'grass' A1(4) = 'ptos' A1(5) = 'ptog' A1(6) = 'stog' A1(7) = 'gtop' A1(8) = 'gtos' A1(9) = 'frac_primf' A1(10) = 'frac_forest' A1(11) = 'pharv' A1(12) = 'smharv' A1(13) = 'syharv' AI1(1) = 'n_event' A2(1) = 'freq_age_primary' A2(2) = 'freq_age_secondary' A2(3) = 'biomass_age_primary' A2(4) = 'biomass_age_secondary' A3(1) = 'area_history_secdf' AI3(1) = 'age_history_secdf' A4(1) = 'FHarvest' A4(2) = 'FClearance' A4(3) = 'FNEP' A4(4) = 'CLitt' A4(5) = 'CSoil' A4(6) = 'CBiomass' A4(7) = 'FTransferNet' A5(1) = 'FTransferGross' A6(1) = 'HarvProd' A6(2) = 'ClearProd' A6(3) = 'HarvProdLoss' A6(4) = 'ClearProdLoss' DO g=1,mp IF (SUM(POPLUC%freq_age_secondary(g,:)).GT.1e-12 ) THEN freq_age_secondary(g,:) = POPLUC%freq_age_secondary(g,:)/SUM(POPLUC%freq_age_secondary(g,:)) ELSE freq_age_secondary(g,:) = 0.0 ENDIF ENDDO CNT = CNT + 1 IF ( CALL1 ) THEN ! Get File-Name WRITE( dum, FMT="(I4,'_',I4)")CABLE_USER%YEARSTART,CABLE_USER%YEAREND IF (CABLE_USER%YEARSTART.LT.1000.AND.CABLE_USER%YEAREND.LT.1000) THEN WRITE( dum, FMT="(I3,'_',I3)")CABLE_USER%YEARSTART,CABLE_USER%YEAREND ELSEIF (CABLE_USER%YEARSTART.LT.1000) THEN WRITE( dum, FMT="(I3,'_',I4)")CABLE_USER%YEARSTART,CABLE_USER%YEAREND ENDIF fname = TRIM(filename%path)//'/'//TRIM(cable_user%RunIden)//'_'//& TRIM(dum)//'_LUC_out.nc' ! Create NetCDF file: STATUS = NF90_create(fname, NF90_CLOBBER, FILE_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Put the file in define mode: STATUS = NF90_redef(FILE_ID) STATUS = NF90_PUT_ATT( FILE_ID, NF90_GLOBAL, "StartYear", CABLE_USER%YEARSTART ) STATUS = NF90_PUT_ATT( FILE_ID, NF90_GLOBAL, "EndYear" , CABLE_USER%YEAREND ) STATUS = NF90_PUT_ATT( FILE_ID, NF90_GLOBAL, "RunIden" , CABLE_USER%RunIden ) ! Define dimensions: ! Land (number of points) STATUS = NF90_def_dim(FILE_ID, 'land' , mp , land_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'mage' , age_max , age_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'mhist',LENGTH_SECDF_HISTORY , hist_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'mLU', nLU , nLU_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'mTrans',nTrans , nTrans_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'nprod',nprod , nprod_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'time' , NF90_UNLIMITED, t_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Define variables STATUS = NF90_def_var(FILE_ID,'time' ,NF90_INT,(/t_ID/),VIDtime ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) DO i = 1, SIZE(A0) STATUS = NF90_def_var(FILE_ID,TRIM(A0(i)) ,NF90_FLOAT,(/land_ID/),VID0(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO DO i = 1, SIZE(A1) STATUS = NF90_def_var(FILE_ID,TRIM(A1(i)) ,NF90_FLOAT,(/land_ID,t_ID/),VID1(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO DO i = 1, SIZE(AI1) STATUS = NF90_def_var(FILE_ID,TRIM(AI1(i)) ,NF90_INT,(/land_ID,t_ID/),VIDI1(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO IF(put_age_vars) THEN DO i = 1, SIZE(A2) STATUS = NF90_def_var(FILE_ID,TRIM(A2(i)) ,NF90_FLOAT,(/land_ID,age_ID,t_ID/),VID2(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO ENDIF !$ !$ DO i = 1, SIZE(A3) !$ STATUS = NF90_def_var(FILE_ID,TRIM(A3(i)) ,NF90_FLOAT,(/land_ID,hist_ID,t_ID/),VID3(i)) !$ IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) !$ END DO !$ !$ DO i = 1, SIZE(A3) !$ STATUS = NF90_def_var(FILE_ID,TRIM(AI3(i)) ,NF90_INT,(/land_ID,hist_ID,t_ID/),VIDI3(i)) !$ IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) !$ END DO DO i = 1, SIZE(A4) STATUS = NF90_def_var(FILE_ID,TRIM(A4(i)) ,NF90_FLOAT,(/land_ID,nLU_ID,t_ID/),VID4(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO DO i = 1, SIZE(A5) STATUS = NF90_def_var(FILE_ID,TRIM(A5(i)) ,NF90_FLOAT,(/land_ID,ntrans_ID,t_ID/),VID5(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO DO i = 1, SIZE(A6) STATUS = NF90_def_var(FILE_ID,TRIM(A6(i)) ,NF90_FLOAT,(/land_ID,nprod_ID,t_ID/),VID6(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO ! End define mode: STATUS = NF90_enddef(FILE_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! PUT LAT / LON ( mp ) STATUS = NF90_PUT_VAR(FILE_ID, VID0(1), POPLUC%latitude ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID0(2), POPLUC%longitude ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) CALL1 = .FALSE. ENDIF ! CALL1 ! TIME ( t ) STATUS = NF90_PUT_VAR(FILE_ID, VIDtime, ctime, start=(/ CNT /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ! PUT 2D VARS ( mp, t ) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 1), POPLUC%primf, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 2), POPLUC%secdf, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 3), POPLUC%grass, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 4), POPLUC%ptos, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 5), POPLUC%ptog, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 6), POPLUC%stog, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 7), POPLUC%gtop, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 8), POPLUC%gtos, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 9), POPLUC%frac_primf, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1(10), POPLUC%frac_forest, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 11), POPLUC%pharv, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 12), POPLUC%smharv, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 13), POPLUC%syharv, start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VIDI1(1), POPLUC%n_event, & start=(/ 1, CNT /), count=(/ mp, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) IF (put_age_vars) THEN ! PUT 3D VARS ( mp, mage, t ) STATUS = NF90_PUT_VAR(FILE_ID, VID2(1), POPLUC%freq_age_primary, & start=(/ 1,1,CNT /), count=(/ mp,age_max,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID2(2),freq_age_secondary, & start=(/ 1,1,CNT /), count=(/ mp,age_max,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID2(3), POPLUC%biomass_age_primary, & start=(/ 1,1,CNT /), count=(/ mp,age_max,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID2(4), POPLUC%biomass_age_secondary, & start=(/ 1,1,CNT /), count=(/ mp,age_max,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ENDIF ! PUT 3D VARS ( mp, nLU, t ) STATUS = NF90_PUT_VAR(FILE_ID, VID4(1), POPLUC%FHarvest, & start=(/ 1,1,CNT /), count=(/ mp,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID4(2), POPLUC%FClearance, & start=(/ 1,1,CNT /), count=(/ mp,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID4(3), POPLUC%FNEP, & start=(/ 1,1,CNT /), count=(/ mp,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID4(4), POPLUC%Clitt, & start=(/ 1,1,CNT /), count=(/ mp,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID4(5), POPLUC%CSoil, & start=(/ 1,1,CNT /), count=(/ mp,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID4(6), POPLUC%CBiomass, & start=(/ 1,1,CNT /), count=(/ mp,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID4(7), POPLUC%FTransferNet, & start=(/ 1,1,CNT /), count=(/ mp,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ! PUT 3D VARS ( mp, nTrans, t ) STATUS = NF90_PUT_VAR(FILE_ID, VID5(1), POPLUC%FTransferGross, & start=(/ 1,1,CNT /), count=(/ mp,nTrans,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ! PUT 3D VARS ( mp, nprod, t ) STATUS = NF90_PUT_VAR(FILE_ID, VID6(1), POPLUC%HarvProd, & start=(/ 1,1,CNT /), count=(/ mp,nprod,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID6(2), POPLUC%ClearProd, & start=(/ 1,1,CNT /), count=(/ mp,nprod,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID6(3), POPLUC%HarvProdLoss, & start=(/ 1,1,CNT /), count=(/ mp,nprod,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID6(4), POPLUC%ClearProdLoss, & start=(/ 1,1,CNT /), count=(/ mp,nprod,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) IF ( FINAL ) THEN ! Close NetCDF file: STATUS = NF90_close(FILE_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) WRITE(*,*) " POPLUC Output written to ",fname ENDIF END SUBROUTINE WRITE_LUC_OUTPUT_NC !============================================================================= SUBROUTINE WRITE_LUC_RESTART_NC ( POPLUC, ctime ) USE CABLE_COMMON_MODULE, ONLY: filename, cable_user USE casa_ncdf_module, ONLY: HANDLE_ERR USE netcdf IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(IN) :: POPLUC INTEGER, INTENT(IN) :: ctime INTEGER :: STATUS INTEGER :: land_ID, age_ID, nLU_ID, nTrans_ID, i, mp, nprod_ID LOGICAL :: CASAONLY CHARACTER :: CYEAR*4, FNAME*99,dum*50 LOGICAL, SAVE :: CALL1 = .TRUE. ! 1 dim arrays (mp ) CHARACTER(len=20),DIMENSION(2) :: A0 ! 2 dim real arrays (mp) CHARACTER(len=20),DIMENSION(3):: A1 ! 2 dim real arrays (mp,age_max) CHARACTER(len=25),DIMENSION(2) :: A2 ! 2 dim real arrays (mp,nprod) CHARACTER(len=25),DIMENSION(2) :: A3 INTEGER, SAVE :: VIDtime, VID0(SIZE(A0)),VID1(SIZE(A1)) INTEGER, SAVE :: VID2(SIZE(A2)), VID3(SIZE(A3)) INTEGER, SAVE :: FILE_ID, CNT = 0 CHARACTER(len=50) :: RecordDimName INTEGER :: g, nprod mp = POPLUC%np nprod = 3 ! number of product pools A0(1) = 'latitude' A0(2) = 'longitude' A1(1) = 'primf' A1(2) = 'secdf' A1(3) = 'grass' A2(1) = 'freq_age_primary' A2(2) = 'freq_age_secondary' A3(1) = 'HarvProd' A3(2) = 'ClearProd' ! Get File-Name WRITE( dum, FMT="(I4,'_',I4)")CABLE_USER%YEARSTART,CABLE_USER%YEAREND IF (CABLE_USER%YEARSTART.LT.1000.AND.CABLE_USER%YEAREND.LT.1000) THEN WRITE( dum, FMT="(I3,'_',I3)")CABLE_USER%YEARSTART,CABLE_USER%YEAREND ELSEIF (CABLE_USER%YEARSTART.LT.1000) THEN WRITE( dum, FMT="(I3,'_',I4)")CABLE_USER%YEARSTART,CABLE_USER%YEAREND ENDIF fname = TRIM(filename%path)//'/'//TRIM(cable_user%RunIden)//'_'//'LUC_rst.nc' ! Create NetCDF file: STATUS = NF90_create(fname, NF90_CLOBBER, FILE_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Put the file in define mode: STATUS = NF90_redef(FILE_ID) STATUS = NF90_PUT_ATT( FILE_ID, NF90_GLOBAL, "StartYear", CABLE_USER%YEARSTART ) STATUS = NF90_PUT_ATT( FILE_ID, NF90_GLOBAL, "EndYear" , CABLE_USER%YEAREND ) STATUS = NF90_PUT_ATT( FILE_ID, NF90_GLOBAL, "RunIden" , CABLE_USER%RunIden ) ! Define dimensions: ! Land (number of points) STATUS = NF90_def_dim(FILE_ID, 'land' , mp , land_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'mage' , age_max , age_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'nprod' , nprod , nprod_ID) ! Define variables DO i = 1, SIZE(A0) STATUS = NF90_def_var(FILE_ID,TRIM(A0(i)) ,NF90_FLOAT,(/land_ID/),VID0(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO DO i = 1, SIZE(A1) STATUS = NF90_def_var(FILE_ID,TRIM(A1(i)) ,NF90_FLOAT,(/land_ID/),VID1(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO DO i = 1, SIZE(A2) STATUS = NF90_def_var(FILE_ID,TRIM(A2(i)) ,NF90_FLOAT,(/land_ID,age_ID/),VID2(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO DO i = 1, SIZE(A3) STATUS = NF90_def_var(FILE_ID,TRIM(A3(i)) ,NF90_FLOAT,(/land_ID,nprod_ID/),VID3(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO ! End define mode: STATUS = NF90_enddef(FILE_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! PUT LAT / LON ( mp ) STATUS = NF90_PUT_VAR(FILE_ID, VID0(1), POPLUC%latitude ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID0(2), POPLUC%longitude ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ! PUT 2D VARS ( mp, t ) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 1), POPLUC%primf) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 2), POPLUC%secdf) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 3), POPLUC%grass) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ! PUT 3D VARS ( mp, mage, t ) STATUS = NF90_PUT_VAR(FILE_ID, VID2(1), POPLUC%freq_age_primary) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID2(2),POPLUC%freq_age_secondary) ! PUT 3D VARS ( mp, nprod, t ) STATUS = NF90_PUT_VAR(FILE_ID, VID3(1), POPLUC%HarvProd) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID3(2),POPLUC%ClearProd) ! Close NetCDF file: STATUS = NF90_close(FILE_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) WRITE(*,*) " POPLUC Restart written to ",fname END SUBROUTINE WRITE_LUC_RESTART_NC !============================================================================= SUBROUTINE READ_LUC_RESTART_NC (POPLUC) USE CABLE_COMMON_MODULE, ONLY: filename, cable_user USE casa_ncdf_module, ONLY: HANDLE_ERR USE netcdf IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(INOUT) :: POPLUC INTEGER :: STATUS, land_dim, mage_dim, nprod_dim INTEGER :: land_ID, age_ID, nLU_ID, nTrans_ID, i, mp, FILE_ID, dID, nprod, nprod_ID LOGICAL :: CASAONLY CHARACTER :: CYEAR*4, FNAME*99,dum*50 LOGICAL, SAVE :: CALL1 = .TRUE. REAL , ALLOCATABLE :: TMP(:), TMP2(:,:), TMP3(:,:) ! 1 dim arrays (mp ) CHARACTER(len=20),DIMENSION(2) :: A0 ! 2 dim real arrays (mp) CHARACTER(len=20),DIMENSION(3):: A1 ! 2 dim real arrays (mp,age_max) CHARACTER(len=25),DIMENSION(2) :: A2 ! 2 dim real arrays (mp,nprod) CHARACTER(len=25),DIMENSION(2) :: A3 mp = POPLUC%np nprod = 3 ALLOCATE(tmp(mp)) ALLOCATE(tmp2(mp,age_max)) ALLOCATE(tmp3(mp,nprod)) A0(1) = 'latitude' A0(2) = 'longitude' A1(1) = 'primf' A1(2) = 'secdf' A1(3) = 'grass' A2(1) = 'freq_age_primary' A2(2) = 'freq_age_secondary' A3(1) = 'HarvProd' A3(2) = 'ClearProd' !fname = TRIM(filename%path)//'/'//TRIM( cable_user%RunIden )//& ! '_LUC_rst.nc' fname = TRIM(filename%path)//'/'//TRIM(cable_user%RunIden)//'_'//'LUC_rst.nc' STATUS = NF90_OPEN( TRIM(fname), NF90_NOWRITE, FILE_ID ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! DIMS STATUS = NF90_INQ_DIMID( FILE_ID, 'land', dID ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_INQUIRE_DIMENSION( FILE_ID, dID, LEN=land_dim ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_INQ_DIMID( FILE_ID, 'mage', dID ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_INQUIRE_DIMENSION( FILE_ID, dID, LEN=mage_dim ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_INQ_DIMID( FILE_ID, 'nprod', dID ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_INQUIRE_DIMENSION( FILE_ID, dID, LEN=nprod_dim ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! READ 1-dimensional fields DO i = 1, SIZE(A1) STATUS = NF90_INQ_VARID( FILE_ID, A1(i), dID ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_GET_VAR( FILE_ID, dID, TMP ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) SELECT CASE ( TRIM(A1(i))) CASE ('primf' ) ; POPLUC%primf = TMP CASE ('secdf' ) ; POPLUC%secdf = TMP CASE ('grass' ) ; POPLUC%grass = TMP END SELECT END DO 991 FORMAT(1000(e12.4,2x)) ! READ 2-dimensional fields (mage) DO i = 1, SIZE(A2) STATUS = NF90_INQ_VARID( FILE_ID, A2(i), dID ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_GET_VAR( FILE_ID, dID, TMP2) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) SELECT CASE ( TRIM(A2(i))) CASE ('freq_age_primary' ) ; POPLUC%freq_age_primary = TMP2 CASE ('freq_age_secondary' ) ; POPLUC%freq_age_secondary = TMP2 END SELECT END DO ! READ 3-dimensional fields (nprod) DO i = 1, SIZE(A3) STATUS = NF90_INQ_VARID( FILE_ID, A3(i), dID ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_GET_VAR( FILE_ID, dID, TMP3) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) SELECT CASE ( TRIM(A3(i))) CASE ('HarvProd' ) ; POPLUC%HarvProd = TMP3 CASE ('ClearProd' ) ; POPLUC%ClearProd = TMP3 END SELECT END DO STATUS = NF90_CLOSE( FILE_ID ) END SUBROUTINE READ_LUC_RESTART_NC !============================================================================ SUBROUTINE WRITE_LUC_OUTPUT_GRID_NC ( POPLUC, ctime, FINAL ) USE cable_IO_vars_module, ONLY: mask, xdimsize, ydimsize , lat_all, lon_all USE CABLE_COMMON_MODULE, ONLY: filename, cable_user USE casa_ncdf_module, ONLY: HANDLE_ERR USE netcdf IMPLICIT NONE TYPE(POPLUC_TYPE), INTENT(IN) :: POPLUC LOGICAL, INTENT(IN) :: FINAL INTEGER, INTENT(IN) :: ctime INTEGER :: STATUS INTEGER :: land_ID, age_ID, hist_ID, t_ID, nLU_ID, nTrans_ID, xID, yID INTEGER :: xvID, yvID INTEGER :: i, mp, nprod, nprod_ID LOGICAL :: CASAONLY CHARACTER :: CYEAR*4, FNAME*99,dum*50 LOGICAL, SAVE :: CALL1 = .TRUE. ! 1 dim arrays (mp ) CHARACTER(len=20),DIMENSION(2) :: A0 ! 2 dim real arrays (mp,t) CHARACTER(len=20),DIMENSION(13):: A1 ! 2 dim integer arrays (mp,t) CHARACTER(len=20),DIMENSION(1):: AI1 ! 3 dim real arrays (mp,age_max,t) CHARACTER(len=25),DIMENSION(4) :: A2 ! 3 dim real arrays (mp,LENGTH_SECDF_HISTORY,t) CHARACTER(len=20),DIMENSION(1) :: A3 ! 3 dim integer arrays (mp,LENGTH_SECDF_HISTORY,t) CHARACTER(len=20),DIMENSION(1) :: AI3 ! 3 dim real arrays (mp,nLU,t) CHARACTER(len=20),DIMENSION(7) :: A4 ! 3 dim real arrays (mp,nTrans,t) CHARACTER(len=20),DIMENSION(1) :: A5 ! 3 dim real arrays (mp,nprod,t) CHARACTER(len=20),DIMENSION(4) :: A6 INTEGER, SAVE :: VIDtime, VID0(SIZE(A0)),VID1(SIZE(A1)),VIDI1(SIZE(AI1)) INTEGER, SAVE :: VID2(SIZE(A2)),VID3(SIZE(A3)),VIDI3(SIZE(AI3)) INTEGER, SAVE :: VID4(SIZE(A4)),VID5(SIZE(A5)), VID6(SIZE(A6)) INTEGER, SAVE :: FILE_ID, CNT = 0, latID, lonID CHARACTER(len=50) :: RecordDimName REAL, ALLOCATABLE :: freq_age_secondary(:,:) INTEGER :: g, k LOGICAL :: put_age_vars REAL(dp) :: ncmissingr = -1.0e+33 INTEGER :: ncmissingi = -9999999 ! REAL :: missing_value = -999999.0 ! for netcdf output INTEGER :: nx, ny LOGICAL,DIMENSION(:,:),ALLOCATABLE :: LandMask ! Logical landmask, true for land, false for non-land REAL(dp), DIMENSION(:,:), ALLOCATABLE :: fieldr INTEGER, DIMENSION(:,:), ALLOCATABLE :: fieldi REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: tmparr1 REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: tmparr2 REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: tmparr3 mp = POPLUC%np nx = xdimsize ny = ydimsize nprod = 3 put_age_vars=.FALSE. ALLOCATE(freq_age_secondary(mp,age_max)) ALLOCATE( landmask ( xdimsize, ydimsize) ) ! Local use in this routine (integer) ALLOCATE( fieldr ( xdimsize, ydimsize) ) ! field for UNPACK command (reals) ALLOCATE( fieldi ( xdimsize, ydimsize) ) ! field for UNPACK command (integers) ALLOCATE( tmparr1( xdimsize, ydimsize, nLU) ) ALLOCATE( tmparr2( xdimsize, ydimsize, nTrans) ) ALLOCATE( tmparr3( xdimsize, ydimsize, nprod) ) ! Convert the integer 'mask' into the logical 'landmask' WHERE (mask .EQ. 1 ) landmask = .TRUE. ELSEWHERE landmask = .FALSE. END WHERE fieldr = ncmissingr fieldi = ncmissingi A0(1) = 'local_latitude' A0(2) = 'local_longitude' A1(1) = 'primf' A1(2) = 'secdf' A1(3) = 'grass' A1(4) = 'ptos' A1(5) = 'ptog' A1(6) = 'stog' A1(7) = 'gtop' A1(8) = 'gtos' A1(9) = 'frac_primf' A1(10) = 'frac_forest' A1(11) = 'pharv' A1(12) = 'smharv' A1(13) = 'syharv' AI1(1) = 'n_event' A2(1) = 'freq_age_primary' A2(2) = 'freq_age_secondary' A2(3) = 'biomass_age_primary' A2(4) = 'biomass_age_secondary' A3(1) = 'area_history_secdf' AI3(1) = 'age_history_secdf' A4(1) = 'FHarvest' A4(2) = 'FClearance' A4(3) = 'FNEP' A4(4) = 'CLitt' A4(5) = 'CSoil' A4(6) = 'CBiomass' A4(7) = 'FTransferNet' A5(1) = 'FTransferGross' A6(1) = 'HarvProd' A6(2) = 'ClearProd' A6(3) = 'HarvProdLoss' A6(4) = 'ClearProdLoss' DO g=1,mp IF (SUM(POPLUC%freq_age_secondary(g,:)).GT.1e-12 ) THEN freq_age_secondary(g,:) = POPLUC%freq_age_secondary(g,:)/SUM(POPLUC%freq_age_secondary(g,:)) ELSE freq_age_secondary(g,:) = 0.0 ENDIF ENDDO CNT = CNT + 1 IF ( CALL1 ) THEN ! Get File-Name WRITE( dum, FMT="(I4,'_',I4)")CABLE_USER%YEARSTART,CABLE_USER%YEAREND IF (CABLE_USER%YEARSTART.LT.1000.AND.CABLE_USER%YEAREND.LT.1000) THEN WRITE( dum, FMT="(I3,'_',I3)")CABLE_USER%YEARSTART,CABLE_USER%YEAREND ELSEIF (CABLE_USER%YEARSTART.LT.1000) THEN WRITE( dum, FMT="(I3,'_',I4)")CABLE_USER%YEARSTART,CABLE_USER%YEAREND ENDIF fname = TRIM(filename%path)//'/'//TRIM(cable_user%RunIden)//'_'//& TRIM(dum)//'_LUC_out.nc' ! Create NetCDF file: STATUS = NF90_create(fname, NF90_CLOBBER, FILE_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Put the file in define mode: STATUS = NF90_redef(FILE_ID) STATUS = NF90_PUT_ATT( FILE_ID, NF90_GLOBAL, "StartYear", CABLE_USER%YEARSTART ) STATUS = NF90_PUT_ATT( FILE_ID, NF90_GLOBAL, "EndYear" , CABLE_USER%YEAREND ) STATUS = NF90_PUT_ATT( FILE_ID, NF90_GLOBAL, "RunIden" , CABLE_USER%RunIden ) ! Define dimensions: STATUS = NF90_DEF_DIM(FILE_ID, 'x', xdimsize, xID) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining x dimension in output file. '// & '(SUBROUTINE open_output_file)') STATUS = NF90_DEF_DIM(FILE_ID, 'y', ydimsize, yID) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining y dimension in output file. '// & '(SUBROUTINE open_output_file)') STATUS = NF90_def_dim(FILE_ID, 'mage' , age_max , age_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'mhist',LENGTH_SECDF_HISTORY , hist_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'mLU', nLU , nLU_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'mTrans',nTrans , nTrans_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'nprod',nprod , nprod_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_def_dim(FILE_ID, 'time' , NF90_UNLIMITED, t_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Define variables STATUS = NF90_def_var(FILE_ID,'time' ,NF90_INT,(/t_ID/),VIDtime ) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) STATUS = NF90_DEF_VAR(FILE_ID, 'latitude', NF90_FLOAT, (/xID, yID/), latID) STATUS = NF90_DEF_VAR(FILE_ID, 'longitude', NF90_FLOAT, (/xID, yID/), lonID) STATUS = NF90_PUT_ATT(FILE_ID, latID, 'units', 'degrees_north') IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining latitude variable attributes in output file. '// & '(SUBROUTINE open_output_file)') STATUS = NF90_PUT_ATT(FILE_ID, lonID, 'units', 'degrees_east') IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining longitude variable attributes in output file. '// & '(SUBROUTINE open_output_file)') ! Write "cordinate variables" to enable reading by GrADS: STATUS = NF90_DEF_VAR(FILE_ID, 'x', NF90_FLOAT, (/xID/), xvID) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining "x" variable (for GrADS) in output file. '// & '(SUBROUTINE open_output_file)') STATUS = NF90_PUT_ATT(FILE_ID, xvID, 'units', 'degrees_east') IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error writing x coordinate variable (GrADS) units in output '// & 'file. (SUBROUTINE open_output_file)') STATUS = NF90_PUT_ATT(FILE_ID, xvID, 'comment', & 'x coordinate variable for GrADS compatibility') IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error writing x variables comment in output file. '// & '(SUBROUTINE open_output_file)') STATUS = NF90_DEF_VAR(FILE_ID, 'y', NF90_FLOAT, (/yID/), yvID) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining "y" variable (for GrADS) in output file. '// & '(SUBROUTINE open_output_file)') STATUS = NF90_PUT_ATT(FILE_ID, yvID, 'units', 'degrees_north') IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error writing y coordinate variable (GrADS) units in output '// & 'file. (SUBROUTINE open_output_file)') STATUS = NF90_PUT_ATT(FILE_ID, yvID, 'comment', & 'y coordinate variable for GrADS compatibility') IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error writing y variables comment in output file. '// & '(SUBROUTINE open_output_file)') DO i = 1, SIZE(A0) STATUS = NF90_def_var(FILE_ID,TRIM(A0(i)) ,NF90_FLOAT,(/xID, yID/),VID0(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VID0(i), '_FillValue', REAL(ncmissingr,4)) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//A0(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VID0(i), 'missing_value', REAL(ncmissingr, 4)) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//A0(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') END DO DO i = 1, SIZE(A1) STATUS = NF90_def_var(FILE_ID,TRIM(A1(i)) ,NF90_FLOAT,(/xID, yID,t_ID/),VID1(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VID1(i), '_FillValue', REAL(ncmissingr, 4)) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//A1(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VID1(i), 'missing_value', REAL(ncmissingr, 4)) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//A1(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') END DO DO i = 1, SIZE(AI1) STATUS = NF90_def_var(FILE_ID,TRIM(AI1(i)) ,NF90_INT,(/xID, yID,t_ID/),VIDI1(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VIDI1(i), '_FillValue',ncmissingi) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//AI1(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VIDI1(i), 'missing_value', ncmissingi) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//AI1(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') END DO IF(put_age_vars) THEN DO i = 1, SIZE(A2) STATUS = NF90_def_var(FILE_ID,TRIM(A2(i)) ,NF90_FLOAT,(/xID, yID,age_ID,t_ID/),VID2(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) END DO ENDIF DO i = 1, SIZE(A4) STATUS = NF90_def_var(FILE_ID,TRIM(A4(i)) ,NF90_FLOAT,(/xID, yID,nLU_ID,t_ID/),VID4(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VID4(i), '_FillValue', REAL(ncmissingr, 4)) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//A4(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VID4(i), 'missing_value', REAL(ncmissingr, 4)) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//A4(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') END DO DO i = 1, SIZE(A5) STATUS = NF90_def_var(FILE_ID,TRIM(A5(i)) ,NF90_FLOAT,(/xID, yID,ntrans_ID,t_ID/),VID5(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VID5(i), '_FillValue', REAL(ncmissingr, 4)) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//A5(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VID5(i), 'missing_value', REAL(ncmissingr, 4)) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//A5(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') END DO DO i = 1, SIZE(A6) STATUS = NF90_def_var(FILE_ID,TRIM(A6(i)) ,NF90_FLOAT,(/xID, yID ,nprod_ID,t_ID/),VID6(i)) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VID6(i), '_FillValue', REAL(ncmissingr, 4)) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//A6(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') ! Define missing/fill values: STATUS = NF90_PUT_ATT(FILE_ID, VID6(i), 'missing_value', REAL(ncmissingr, 4)) IF (STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error defining '//A6(i)//' variable attributes in output file. '// & '(INTERFACE define_ovar)') END DO ! End define mode: STATUS = NF90_enddef(FILE_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) ! Write GrADS coordinate variables STATUS = NF90_PUT_VAR(FILE_ID, xvID, REAL(lon_all(:, 1), 4)) IF(STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error writing GrADS x coordinate variable to ' & //TRIM(filename%out)// ' (SUBROUTINE open_output_file)') STATUS = NF90_PUT_VAR(FILE_ID, yvID, REAL(lat_all(1, :), 4)) IF(STATUS /= NF90_NOERR) CALL handle_err & (STATUS, 'Error writing GrADS y coordinate variable to ' & //TRIM(filename%out)// ' (SUBROUTINE open_output_file)') ! PUT LAT / LON ( mp ) STATUS = NF90_PUT_VAR(FILE_ID, VID0(1), UNPACK(POPLUC%latitude,landmask, fieldr )) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID0(2), UNPACK(POPLUC%longitude, landmask, fieldr) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ! PUT LAT / LON ( mp ) STATUS = NF90_PUT_VAR(FILE_ID, latID, lat_all) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, lonID, lon_all ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) CALL1 = .FALSE. ENDIF ! CALL1 ! TIME ( t ) STATUS = NF90_PUT_VAR(FILE_ID, VIDtime, ctime, start=(/ CNT /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ! PUT 2D VARS ( mp, t ) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 1), UNPACK(POPLUC%primf,landmask, fieldr), & start=(/ 1, 1, CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 2), UNPACK(POPLUC%secdf,landmask, fieldr), & start=(/ 1, 1, CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 3), UNPACK(POPLUC%grass,landmask, fieldr), & start=(/ 1, 1, CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 4), UNPACK(POPLUC%ptos,landmask, fieldr), & start=(/ 1, 1, CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 5), UNPACK(POPLUC%ptog,landmask, fieldr), & start=(/ 1,1, CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 6), UNPACK(POPLUC%stog, landmask, fieldr), & start=(/ 1, 1,CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 7), UNPACK(POPLUC%gtop,landmask, fieldr), & start=(/ 1, 1, CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 8), UNPACK(POPLUC%gtos,landmask, fieldr), & start=(/ 1, 1,CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 9), UNPACK(POPLUC%frac_primf,landmask, fieldr), & start=(/ 1, 1,CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1(10), UNPACK(POPLUC%frac_forest,landmask, fieldr), & start=(/ 1, 1, CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 11), UNPACK(POPLUC%pharv, landmask, fieldr), & start=(/ 1, 1, CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 12), UNPACK(POPLUC%smharv,landmask, fieldr), & start=(/ 1, 1, CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VID1( 13), UNPACK(POPLUC%syharv,landmask, fieldr), & start=(/ 1, 1,CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) STATUS = NF90_PUT_VAR(FILE_ID, VIDI1(1), UNPACK(POPLUC%n_event, landmask, fieldi), & start=(/ 1,1, CNT /), count=(/ nx,ny, 1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ! PUT 3D VARS ( mp, nLU, t ) DO k=1,nLU tmparr1(:,:,k) = UNPACK(POPLUC%FHarvest(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID4(1), tmparr1, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) DO k=1,nLU tmparr1(:,:,k) = UNPACK(POPLUC%FClearance(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID4(2), tmparr1, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) DO k=1,nLU tmparr1(:,:,k) = UNPACK(POPLUC%FNEP(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID4(3), tmparr1, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) DO k=1,nLU tmparr1(:,:,k) = UNPACK(POPLUC%Clitt(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID4(4), tmparr1, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) DO k=1,nLU tmparr1(:,:,k) = UNPACK(POPLUC%Csoil(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID4(5), tmparr1, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) DO k=1,nLU tmparr1(:,:,k) = UNPACK(POPLUC%CBiomass(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID4(6), tmparr1, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) DO k=1,nLU tmparr1(:,:,k) = UNPACK(POPLUC%FTransferNet(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID4(7), tmparr1, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nLU,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ! PUT 3D VARS ( nx,ny, nTrans, t ) DO k=1,nTrans tmparr2(:,:,k) = UNPACK(POPLUC%FTransferGross(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID5(1), tmparr2, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nTrans,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) ! PUT 3D VARS ( nx,ny, nprod, t ) DO k=1,nprod tmparr3(:,:,k) = UNPACK(POPLUC%HarvProd(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID6(1), tmparr3, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nprod,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) DO k=1,nprod tmparr3(:,:,k) = UNPACK(POPLUC%ClearProd(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID6(2), tmparr3, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nprod,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) DO k=1,nprod tmparr3(:,:,k) = UNPACK(POPLUC%HarvProdLoss(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID6(3), tmparr3, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nprod,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) DO k=1,nprod tmparr3(:,:,k) = UNPACK(POPLUC%ClearProdLoss(:,k),landmask, fieldr) ENDDO STATUS = NF90_PUT_VAR(FILE_ID, VID6(4), tmparr3, & start=(/ 1,1,1,CNT /), count=(/ nx,ny,nprod,1 /) ) IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS) IF ( FINAL ) THEN ! Close NetCDF file: STATUS = NF90_close(FILE_ID) IF (STATUS /= NF90_noerr) CALL handle_err(STATUS) WRITE(*,*) " POPLUC Output written to ",fname ENDIF END SUBROUTINE WRITE_LUC_OUTPUT_GRID_NC !============================================================================= END MODULE POPLUC_MODULE