diff --git a/core/biogeochem/POP.F90 b/core/biogeochem/POP.F90
index 5c62059cb..9bec02b1c 100755
--- a/core/biogeochem/POP.F90
+++ b/core/biogeochem/POP.F90
@@ -270,6 +270,7 @@ MODULE POP_Types
      REAL(dp), DIMENSION(NYEAR_HISTORY) :: cat_mortality_history
      REAL(dp), DIMENSION(AGEMAX) :: freq_age ! age weighting (by age in y: 0:AGE_MAX-1)
      REAL(dp), DIMENSION(AGEMAX) :: biomass_age
+     REAL(dp) :: rkill
   END TYPE Landscape
 
   TYPE POP_TYPE
@@ -2844,7 +2845,7 @@ SUBROUTINE INTERPOLATE_FIREMORTALITY(pop, disturbance_interval,it,g)
                   pop%pop_grid(g)%patch(:)%age(1).LT.age_max)
           endif
 
-
+          tmp5_min = 0.0_dp
           if (i_min.ne.0.and.age_min.ne.0) then
              MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_min
              where (MASK)
@@ -2852,19 +2853,24 @@ SUBROUTINE INTERPOLATE_FIREMORTALITY(pop, disturbance_interval,it,g)
              elsewhere
                 tmp_array = 0.0_dp
              endwhere
-             tmp5_min = SUM(pop%pop_grid(g)%patch(:)%fire_mortality,MASK)/SUM(tmp_array)
+             if ( any(MASK) ) then
+               tmp5_min = SUM(pop%pop_grid(g)%patch(:)%fire_mortality,MASK)/SUM(tmp_array)
+             endif
           else
              tmp5_min = 0.0_dp
           endif
 
+          tmp5_max = 0.0_dp
           MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_max
           where (MASK)
              tmp_array = 1.0_dp
           elsewhere
              tmp_array = 0.0_dp
           endwhere
-          tmp5_max = SUM(pop%pop_grid(g)%patch(:)%fire_mortality,MASK)/SUM(tmp_array)
-
+          if ( any(MASK) ) then
+            tmp5_max = sum(pop%pop_grid(g)%patch(:)%fire_mortality,MASK)/sum(tmp_array)
+          endif
+          
           fire_mort_age(iage) = tmp5_min + (tmp5_max-tmp5_min)/real(age_max-age_min,dp)* &
                real(age(iage)-age_min,dp)
 
@@ -2955,6 +2961,17 @@ SUBROUTINE ADJUST_POP_FOR_FIRE(pop,disturbance_interval, burned_area, FLI)
        ! creates new value for  POP%pop_grid(g)%fire_mortality
        CALL INTERPOLATE_FIREMORTALITY(pop, disturbance_interval,it,g)
 
+       !CLN Kill ratio to be used within BLAZE to compute fluxes
+       POP%pop_grid(g)%rkill = 0.
+       if ( POP%pop_grid(g)%cmass_sum .gt. 0.) then
+          POP%pop_grid(g)%rkill = POP%pop_grid(g)%fire_mortality / POP%pop_grid(g)%cmass_sum
+       else
+          POP%pop_grid(g)%rkill = 0.
+       endif
+       if (POP%pop_grid(g)%rkill .GT. 0.) then
+          WRITE(*,*) "driver rkill",g,POP%pop_grid(g)%rkill
+       endif
+
        POP%pop_grid(g)%cmass_sum = POP%pop_grid(g)%cmass_sum - POP%pop_grid(g)%fire_mortality
 
 
diff --git a/core/biogeochem/casa_inout.F90 b/core/biogeochem/casa_inout.F90
index 83665f1dd..1be7c861d 100644
--- a/core/biogeochem/casa_inout.F90
+++ b/core/biogeochem/casa_inout.F90
@@ -1608,43 +1608,45 @@ subroutine write_casa_output_nc(veg, casamet, casapool, casabal, casaflux, casao
     a1(22) = 'CallocLeaf'
     a1(23) = 'CallocWood'
     a1(24) = 'CallocFineRoot'
-    na1 = 24
+    a1(25) = 'fluxCtoCO2_plant_fire'
+    a1(26) = 'fluxCtoCO2_litter_fire'
+    na1 = 26
     ! N
-    a1(25) = 'sumnbal'
-    a1(26) = 'Nminfix'
-    a1(27) = 'Nmindep'
-    a1(28) = 'Nminloss'
-    a1(29) = 'Nminleach'
-    a1(30) = 'Nupland'
-    a1(31) = 'Nlittermin'
-    a1(32) = 'Nsmin'
-    a1(33) = 'Nsimm'
-    a1(34) = 'Nsnet'
-    a1(35) = 'fNMinloss'
-    a1(36) = 'Nsoilmin'
-    if (icycle==2) na1 = 36
+    a1(27) = 'sumnbal'
+    a1(28) = 'Nminfix'
+    a1(29) = 'Nmindep'
+    a1(30) = 'Nminloss'
+    a1(31) = 'Nminleach'
+    a1(32) = 'Nupland'
+    a1(33) = 'Nlittermin'
+    a1(34) = 'Nsmin'
+    a1(35) = 'Nsimm'
+    a1(36) = 'Nsnet'
+    a1(37) = 'fNMinloss'
+    a1(38) = 'Nsoilmin'
+    if (icycle==2) na1 = 38
     ! P
-    a1(37) = 'psoillab'
-    a1(38) = 'psoilsorb'
-    a1(39) = 'psoilocc'
-    a1(40) = 'sumpbal'
-    a1(41) = 'Plabuptake'
-    a1(42) = 'Pdep'
-    a1(43) = 'pwea'
-    a1(44) = 'Pleach'
-    a1(45) = 'Ploss'
-    a1(46) = 'Pupland'
-    a1(47) = 'Plittermin'
-    a1(48) = 'Psmin'
-    a1(49) = 'Psimm'
-    a1(50) = 'Psnet'
-    a1(51) = 'fPleach'
-    a1(52) = 'kPlab'
-    a1(53) = 'kPsorb'
-    a1(54) = 'kpocc'
-    a1(55) = 'kmlabP'
-    a1(56) = 'Psorbmax'
-    if (icycle==3) na1 = 56
+    a1(39) = 'psoillab'
+    a1(40) = 'psoilsorb'
+    a1(41) = 'psoilocc'
+    a1(42) = 'sumpbal'
+    a1(43) = 'Plabuptake'
+    a1(44) = 'Pdep'
+    a1(45) = 'pwea'
+    a1(46) = 'Pleach'
+    a1(47) = 'Ploss'
+    a1(48) = 'Pupland'
+    a1(49) = 'Plittermin'
+    a1(50) = 'Psmin'
+    a1(51) = 'Psimm'
+    a1(52) = 'Psnet'
+    a1(53) = 'fPleach'
+    a1(54) = 'kPlab'
+    a1(55) = 'kPsorb'
+    a1(56) = 'kpocc'
+    a1(57) = 'kmlabP'
+    a1(58) = 'Psorbmax'
+    if (icycle==3) na1 = 58
 
     ! C
     a2(1) = 'cplant'
@@ -1917,74 +1919,79 @@ subroutine write_casa_output_nc(veg, casamet, casapool, casabal, casaflux, casao
     if(status /= nf90_noerr) call handle_err(status)
     status = nf90_put_var(file_id, vid1(24), real(casaflux%CallocFineRoot,sp), start=(/1,cnt/), count=(/mp,1/) )
     if(status /= nf90_noerr) call handle_err(status)
+    !BLAZE outputs
+    status = nf90_put_var(file_id, vid1(25), real(casaflux%fluxCtoCO2_plant_fire,sp), start=(/1,cnt/), count=(/mp,1/) )
+    if(status /= nf90_noerr) call handle_err(status)
+    status = nf90_put_var(file_id, vid1(26), real(casaflux%fluxCtoCO2_litter_fire,sp), start=(/1,cnt/), count=(/mp,1/) )
+    if(status /= nf90_noerr) call handle_err(status)
     ! N
     if (icycle > 1) then
-       status = nf90_put_var(file_id, vid1(25), real(casabal%sumnbal,sp),      start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(27), real(casabal%sumnbal,sp),      start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(26), real(casaflux%nminfix,sp),     start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(28), real(casaflux%nminfix,sp),     start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(27), real(casaflux%nmindep,sp),     start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(29), real(casaflux%nmindep,sp),     start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(28), real(casaflux%nminloss,sp),    start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(30), real(casaflux%nminloss,sp),    start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(29), real(casaflux%nminleach,sp),   start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(31), real(casaflux%nminleach,sp),   start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(30), real(casaflux%nupland,sp),     start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(32), real(casaflux%nupland,sp),     start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(31), real(casaflux%nlittermin,sp),  start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(33), real(casaflux%nlittermin,sp),  start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(32), real(casaflux%nsmin,sp),       start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(34), real(casaflux%nsmin,sp),       start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(33), real(casaflux%nsimm,sp),       start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(35), real(casaflux%nsimm,sp),       start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(34), real(casaflux%nsnet,sp),       start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(36), real(casaflux%nsnet,sp),       start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(35), real(casaflux%fnminloss,sp),   start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(37), real(casaflux%fnminloss,sp),   start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(36), real(casapool%nsoilmin,sp),    start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(38), real(casapool%nsoilmin,sp),    start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
     endif
     ! P
     if (icycle > 2) then
-       status = nf90_put_var(file_id, vid1(37), real(casapool%psoillab,sp),    start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(39), real(casapool%psoillab,sp),    start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(38), real(casapool%psoilsorb,sp),   start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(40), real(casapool%psoilsorb,sp),   start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(39), real(casapool%psoilocc,sp),    start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(41), real(casapool%psoilocc,sp),    start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(40), real(casabal%sumpbal,sp),      start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(42), real(casabal%sumpbal,sp),      start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(41), real(casaflux%plabuptake,sp),  start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(43), real(casaflux%plabuptake,sp),  start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(42), real(casaflux%pdep,sp),        start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(44), real(casaflux%pdep,sp),        start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(43), real(casaflux%pwea,sp),        start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(45), real(casaflux%pwea,sp),        start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(44), real(casaflux%pleach,sp),      start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(46), real(casaflux%pleach,sp),      start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(45), real(casaflux%ploss,sp),       start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(47), real(casaflux%ploss,sp),       start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(46), real(casaflux%pupland,sp),     start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(48), real(casaflux%pupland,sp),     start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(47), real(casaflux%plittermin,sp),  start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(49), real(casaflux%plittermin,sp),  start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(48), real(casaflux%psmin,sp),       start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(50), real(casaflux%psmin,sp),       start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(49), real(casaflux%psimm,sp),       start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(51), real(casaflux%psimm,sp),       start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(50), real(casaflux%psnet,sp),       start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(52), real(casaflux%psnet,sp),       start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(51), real(casaflux%fpleach,sp),     start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(53), real(casaflux%fpleach,sp),     start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(52), real(casaflux%kplab,sp),       start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(54), real(casaflux%kplab,sp),       start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(53), real(casaflux%kpsorb,sp),      start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(55), real(casaflux%kpsorb,sp),      start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(54), real(casaflux%kpocc,sp),       start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(56), real(casaflux%kpocc,sp),       start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(55), real(casaflux%kmlabp,sp),      start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(57), real(casaflux%kmlabp,sp),      start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
-       status = nf90_put_var(file_id, vid1(56), real(casaflux%psorbmax,sp),    start=(/1,cnt/), count=(/mp,1/) )
+       status = nf90_put_var(file_id, vid1(58), real(casaflux%psorbmax,sp),    start=(/1,cnt/), count=(/mp,1/) )
        if(status /= nf90_noerr) call handle_err(status)
     endif
 
@@ -2211,7 +2218,9 @@ subroutine write_casa_output_grid_nc(veg, casamet, casapool, casabal, casaflux,
     a1(22) = 'CallocLeaf'
     a1(23) = 'CallocWood'
     a1(24) = 'CallocFineRoot'
-    na1 = 24
+    a1(25) = 'fluxCtoCO2_plant_fire'
+    a1(26) = 'fluxCtoCO2_litter_fire'
+    na1 = 26
 
     u1(3)  = 'm2 m-2' 
     u1(4)  = 'g C m-2 d-1'
@@ -2235,26 +2244,26 @@ subroutine write_casa_output_grid_nc(veg, casamet, casapool, casabal, casaflux,
     u1(22) = 'g C m-2 d-1'
     u1(23) = 'g C m-2 d-1'
     u1(24) = 'g C m-2 d-1'
+    u1(25) = 'g C m-2 d-1'
+    u1(26) = 'g C m-2 d-1'
 
 
     ! N
-    a1(25) = 'sumnbal'
-    a1(26) = 'Nminfix'
-    a1(27) = 'Nmindep'
-    a1(28) = 'Nminloss'
-    a1(29) = 'Nminleach'
-    a1(30) = 'Nupland'
-    a1(31) = 'Nlittermin'
-    a1(32) = 'Nsmin'
-    a1(33) = 'Nsimm'
-    a1(34) = 'Nsnet'
-    a1(35) = 'fNMinloss'
-    a1(36) = 'Nsoilmin'
-    if (icycle==2) na1 = 36
-
-    u1(25) = ''
-    u1(26) = 'g N m-2 d-1'
-    u1(27) = 'g N m-2 d-1'
+    a1(27) = 'sumnbal'
+    a1(28) = 'Nminfix'
+    a1(29) = 'Nmindep'
+    a1(30) = 'Nminloss'
+    a1(31) = 'Nminleach'
+    a1(32) = 'Nupland'
+    a1(33) = 'Nlittermin'
+    a1(34) = 'Nsmin'
+    a1(35) = 'Nsimm'
+    a1(36) = 'Nsnet'
+    a1(37) = 'fNMinloss'
+    a1(38) = 'Nsoilmin'
+    if (icycle==2) na1 = 38
+
+    u1(27) = ''
     u1(28) = 'g N m-2 d-1'
     u1(29) = 'g N m-2 d-1'
     u1(30) = 'g N m-2 d-1'
@@ -2264,37 +2273,37 @@ subroutine write_casa_output_grid_nc(veg, casamet, casapool, casabal, casaflux,
     u1(34) = 'g N m-2 d-1'
     u1(35) = 'g N m-2 d-1'
     u1(36) = 'g N m-2 d-1'
+    u1(37) = 'g N m-2 d-1'
+    u1(38) = 'g N m-2 d-1'
 
 
     ! P
-    a1(37) = 'psoillab'
-    a1(38) = 'psoilsorb'
-    a1(39) = 'psoilocc'
-    a1(40) = 'sumpbal'
-    a1(41) = 'Plabuptake'
-    a1(42) = 'Pdep'
-    a1(43) = 'pwea'
-    a1(44) = 'Pleach'
-    a1(45) = 'Ploss'
-    a1(46) = 'Pupland'
-    a1(47) = 'Plittermin'
-    a1(48) = 'Psmin'
-    a1(49) = 'Psimm'
-    a1(50) = 'Psnet'
-    a1(51) = 'fPleach'
-    a1(52) = 'kPlab'
-    a1(53) = 'kPsorb'
-    a1(54) = 'kpocc'
-    a1(55) = 'kmlabP'
-    a1(56) = 'Psorbmax'
-    if (icycle==3) na1 = 56
-
-    u1(37) = 'g P m-2 d-1'
-    u1(38) = 'g P m-2 d-1'
+    a1(39) = 'psoillab'
+    a1(40) = 'psoilsorb'
+    a1(41) = 'psoilocc'
+    a1(42) = 'sumpbal'
+    a1(43) = 'Plabuptake'
+    a1(44) = 'Pdep'
+    a1(45) = 'pwea'
+    a1(46) = 'Pleach'
+    a1(47) = 'Ploss'
+    a1(48) = 'Pupland'
+    a1(49) = 'Plittermin'
+    a1(50) = 'Psmin'
+    a1(51) = 'Psimm'
+    a1(52) = 'Psnet'
+    a1(53) = 'fPleach'
+    a1(54) = 'kPlab'
+    a1(55) = 'kPsorb'
+    a1(56) = 'kpocc'
+    a1(57) = 'kmlabP'
+    a1(58) = 'Psorbmax'
+    if (icycle==3) na1 = 58
+
     u1(39) = 'g P m-2 d-1'
-    u1(40) = 'sumpbal'
+    u1(40) = 'g P m-2 d-1'
     u1(41) = 'g P m-2 d-1'
-    u1(42) = 'g P m-2 d-1'
+    u1(42) = 'sumpbal'
     u1(43) = 'g P m-2 d-1'
     u1(44) = 'g P m-2 d-1'
     u1(45) = 'g P m-2 d-1'
@@ -2303,12 +2312,14 @@ subroutine write_casa_output_grid_nc(veg, casamet, casapool, casabal, casaflux,
     u1(48) = 'g P m-2 d-1'
     u1(49) = 'g P m-2 d-1'
     u1(50) = 'g P m-2 d-1'
-    u1(51) = ''
-    u1(52) = ''
+    u1(51) = 'g P m-2 d-1'
+    u1(52) = 'g P m-2 d-1'
     u1(53) = ''
     u1(54) = ''
     u1(55) = ''
     u1(56) = ''
+    u1(57) = ''
+    u1(58) = ''
 
     ! C
     a2(1) = 'cplant'
@@ -2625,45 +2636,48 @@ subroutine write_casa_output_grid_nc(veg, casamet, casapool, casabal, casaflux,
     call put_casa_var_grid_patch(file_id, vid1(22), real(casaflux%CallocLeaf,sp), cnt, otmp4patch)
     call put_casa_var_grid_patch(file_id, vid1(23), real(casaflux%CallocWood,sp), cnt, otmp4patch)
     call put_casa_var_grid_patch(file_id, vid1(24), real(casaflux%CallocFineRoot,sp), cnt, otmp4patch)
+    !BLAZE outputs
+    call put_casa_var_grid_patch(file_id, vid1(25), real(casaflux%fluxCtoCO2_plant_fire,sp), cnt, otmp4patch)
+    call put_casa_var_grid_patch(file_id, vid1(26), real(casaflux%fluxCtoCO2_litter_fire,sp), cnt, otmp4patch)
 
     ! N
     if (icycle > 1) then
-       call put_casa_var_grid_patch(file_id, vid1(25), real(casabal%sumnbal,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(26), real(casaflux%nminfix,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(27), real(casaflux%nmindep,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(28), real(casaflux%nminloss,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(29), real(casaflux%nminleach,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(30), real(casaflux%nupland,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(31), real(casaflux%nlittermin,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(32), real(casaflux%nsmin,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(33), real(casaflux%nsimm,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(34), real(casaflux%nsnet,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(35), real(casaflux%fnminloss,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(36), real(casapool%nsoilmin,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(27), real(casabal%sumnbal,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(28), real(casaflux%nminfix,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(29), real(casaflux%nmindep,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(30), real(casaflux%nminloss,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(31), real(casaflux%nminleach,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(32), real(casaflux%nupland,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(33), real(casaflux%nlittermin,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(34), real(casaflux%nsmin,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(35), real(casaflux%nsimm,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(36), real(casaflux%nsnet,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(37), real(casaflux%fnminloss,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(38), real(casapool%nsoilmin,sp), cnt, otmp4patch )
     endif
 
     ! P
     if (icycle > 2) then
-       call put_casa_var_grid_patch(file_id, vid1(37), real(casapool%psoillab,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(38), real(casapool%psoilsorb,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(39), real(casapool%psoilocc,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(40), real(casabal%sumpbal,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(41), real(casaflux%plabuptake,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(42), real(casaflux%pdep,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(43), real(casaflux%pwea,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(44), real(casaflux%pleach,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(45), real(casaflux%ploss,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(46), real(casaflux%pupland,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(47), real(casaflux%plittermin,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(48), real(casaflux%psmin,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(49), real(casaflux%psimm,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(50), real(casaflux%psnet,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(51), real(casaflux%fpleach,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(52), real(casaflux%kplab,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(53), real(casaflux%kpsorb,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(54), real(casaflux%kpocc,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(55), real(casaflux%kmlabp,sp), cnt, otmp4patch )
-       call put_casa_var_grid_patch(file_id, vid1(56), real(casaflux%psorbmax,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(39), real(casapool%psoillab,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(40), real(casapool%psoilsorb,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(41), real(casapool%psoilocc,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(42), real(casabal%sumpbal,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(43), real(casaflux%plabuptake,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(44), real(casaflux%pdep,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(45), real(casaflux%pwea,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(46), real(casaflux%pleach,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(47), real(casaflux%ploss,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(48), real(casaflux%pupland,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(49), real(casaflux%plittermin,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(50), real(casaflux%psmin,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(51), real(casaflux%psimm,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(52), real(casaflux%psnet,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(53), real(casaflux%fpleach,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(54), real(casaflux%kplab,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(55), real(casaflux%kpsorb,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(56), real(casaflux%kpocc,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(57), real(casaflux%kmlabp,sp), cnt, otmp4patch )
+       call put_casa_var_grid_patch(file_id, vid1(58), real(casaflux%psorbmax,sp), cnt, otmp4patch )
     endif
 
     ! 4D vars (x,y,mplant,t)
@@ -2968,7 +2982,9 @@ subroutine write_casa_output_patch_nc(veg, casamet, casapool, casabal, casaflux,
     a1(22) = 'CallocLeaf'
     a1(23) = 'CallocWood'
     a1(24) = 'CallocFineRoot'
-    na1 = 24
+    a1(25) = 'fluxCtoCO2_plant_fire'
+    a1(26) = 'fluxCtoCO2_litter_fire'
+    na1 = 26
 
     u1(3)  = 'm2 m-2' 
     u1(4)  = 'g C m-2 d-1'
@@ -2992,25 +3008,25 @@ subroutine write_casa_output_patch_nc(veg, casamet, casapool, casabal, casaflux,
     u1(22) = 'g C m-2 d-1'
     u1(23) = 'g C m-2 d-1'
     u1(24) = 'g C m-2 d-1'
+    u1(25) = 'g C m-2 d-1'
+    u1(26) = 'g C m-2 d-1'
 
     ! N
-    a1(25) = 'sumnbal'
-    a1(26) = 'Nminfix'
-    a1(27) = 'Nmindep'
-    a1(28) = 'Nminloss'
-    a1(29) = 'Nminleach'
-    a1(30) = 'Nupland'
-    a1(31) = 'Nlittermin'
-    a1(32) = 'Nsmin'
-    a1(33) = 'Nsimm'
-    a1(34) = 'Nsnet'
-    a1(35) = 'fNMinloss'
-    a1(36) = 'Nsoilmin'
-    if (icycle==2) na1 = 36
-
-    u1(25) = ''
-    u1(26) = 'g N m-2 d-1'
-    u1(27) = 'g N m-2 d-1'
+    a1(27) = 'sumnbal'
+    a1(28) = 'Nminfix'
+    a1(29) = 'Nmindep'
+    a1(30) = 'Nminloss'
+    a1(31) = 'Nminleach'
+    a1(32) = 'Nupland'
+    a1(33) = 'Nlittermin'
+    a1(34) = 'Nsmin'
+    a1(35) = 'Nsimm'
+    a1(36) = 'Nsnet'
+    a1(37) = 'fNMinloss'
+    a1(38) = 'Nsoilmin'
+    if (icycle==2) na1 = 38
+
+    u1(27) = ''
     u1(28) = 'g N m-2 d-1'
     u1(29) = 'g N m-2 d-1'
     u1(30) = 'g N m-2 d-1'
@@ -3020,37 +3036,37 @@ subroutine write_casa_output_patch_nc(veg, casamet, casapool, casabal, casaflux,
     u1(34) = 'g N m-2 d-1'
     u1(35) = 'g N m-2 d-1'
     u1(36) = 'g N m-2 d-1'
+    u1(37) = 'g N m-2 d-1'
+    u1(38) = 'g N m-2 d-1'
 
 
     ! P
-    a1(37) = 'psoillab'
-    a1(38) = 'psoilsorb'
-    a1(39) = 'psoilocc'
-    a1(40) = 'sumpbal'
-    a1(41) = 'Plabuptake'
-    a1(42) = 'Pdep'
-    a1(43) = 'pwea'
-    a1(44) = 'Pleach'
-    a1(45) = 'Ploss'
-    a1(46) = 'Pupland'
-    a1(47) = 'Plittermin'
-    a1(48) = 'Psmin'
-    a1(49) = 'Psimm'
-    a1(50) = 'Psnet'
-    a1(51) = 'fPleach'
-    a1(52) = 'kPlab'
-    a1(53) = 'kPsorb'
-    a1(54) = 'kpocc'
-    a1(55) = 'kmlabP'
-    a1(56) = 'Psorbmax'
-    if (icycle==3) na1 = 56
-
-    u1(37) = 'g P m-2 d-1'
-    u1(38) = 'g P m-2 d-1'
+    a1(39) = 'psoillab'
+    a1(40) = 'psoilsorb'
+    a1(41) = 'psoilocc'
+    a1(42) = 'sumpbal'
+    a1(43) = 'Plabuptake'
+    a1(44) = 'Pdep'
+    a1(45) = 'pwea'
+    a1(46) = 'Pleach'
+    a1(47) = 'Ploss'
+    a1(48) = 'Pupland'
+    a1(49) = 'Plittermin'
+    a1(50) = 'Psmin'
+    a1(51) = 'Psimm'
+    a1(52) = 'Psnet'
+    a1(53) = 'fPleach'
+    a1(54) = 'kPlab'
+    a1(55) = 'kPsorb'
+    a1(56) = 'kpocc'
+    a1(57) = 'kmlabP'
+    a1(58) = 'Psorbmax'
+    if (icycle==3) na1 = 58
+
     u1(39) = 'g P m-2 d-1'
-    u1(40) = 'sumpbal'
+    u1(40) = 'g P m-2 d-1'
     u1(41) = 'g P m-2 d-1'
-    u1(42) = 'g P m-2 d-1'
+    u1(42) = 'sumpbal'
     u1(43) = 'g P m-2 d-1'
     u1(44) = 'g P m-2 d-1'
     u1(45) = 'g P m-2 d-1'
@@ -3059,12 +3075,14 @@ subroutine write_casa_output_patch_nc(veg, casamet, casapool, casabal, casaflux,
     u1(48) = 'g P m-2 d-1'
     u1(49) = 'g P m-2 d-1'
     u1(50) = 'g P m-2 d-1'
-    u1(51) = ''
-    u1(52) = ''
+    u1(51) = 'g P m-2 d-1'
+    u1(52) = 'g P m-2 d-1'
     u1(53) = ''
     u1(54) = ''
     u1(55) = ''
     u1(56) = ''
+    u1(57) = ''
+    u1(58) = ''
 
     ! C
     a2(1) = 'cplant'
@@ -3399,44 +3417,46 @@ subroutine write_casa_output_patch_nc(veg, casamet, casapool, casabal, casaflux,
     call put_casa_var_patch(file_id, vid1(22), real(casaflux%CallocLeaf,sp), cnt, otmp3patch)
     call put_casa_var_patch(file_id, vid1(23), real(casaflux%CallocWood,sp), cnt, otmp3patch)
     call put_casa_var_patch(file_id, vid1(24), real(casaflux%CallocFineRoot,sp), cnt, otmp3patch)
+    call put_casa_var_patch(file_id, vid1(25), real(casaflux%fluxCtoCO2_plant_fire,sp), cnt, otmp3patch)
+    call put_casa_var_patch(file_id, vid1(26), real(casaflux%fluxCtoCO2_litter_fire,sp), cnt, otmp3patch)
 
     ! N
     if (icycle > 1) then
-       call put_casa_var_patch(file_id, vid1(25), real(casabal%sumnbal,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(26), real(casaflux%nminfix,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(27), real(casaflux%nmindep,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(28), real(casaflux%nminloss,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(29), real(casaflux%nminleach,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(30), real(casaflux%nupland,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(31), real(casaflux%nlittermin,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(32), real(casaflux%nsmin,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(33), real(casaflux%nsimm,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(34), real(casaflux%nsnet,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(35), real(casaflux%fnminloss,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(36), real(casapool%nsoilmin,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(27), real(casabal%sumnbal,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(28), real(casaflux%nminfix,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(29), real(casaflux%nmindep,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(30), real(casaflux%nminloss,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(31), real(casaflux%nminleach,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(32), real(casaflux%nupland,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(33), real(casaflux%nlittermin,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(34), real(casaflux%nsmin,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(35), real(casaflux%nsimm,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(36), real(casaflux%nsnet,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(37), real(casaflux%fnminloss,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(38), real(casapool%nsoilmin,sp), cnt, otmp3patch )
     endif
     ! P
     if (icycle > 2) then
-       call put_casa_var_patch(file_id, vid1(37), real(casapool%psoillab,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(38), real(casapool%psoilsorb,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(39), real(casapool%psoilocc,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(40), real(casabal%sumpbal,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(41), real(casaflux%plabuptake,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(42), real(casaflux%pdep,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(43), real(casaflux%pwea,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(44), real(casaflux%pleach,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(45), real(casaflux%ploss,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(46), real(casaflux%pupland,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(47), real(casaflux%plittermin,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(48), real(casaflux%psmin,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(49), real(casaflux%psimm,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(50), real(casaflux%psnet,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(51), real(casaflux%fpleach,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(52), real(casaflux%kplab,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(53), real(casaflux%kpsorb,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(54), real(casaflux%kpocc,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(55), real(casaflux%kmlabp,sp), cnt, otmp3patch )
-       call put_casa_var_patch(file_id, vid1(56), real(casaflux%psorbmax,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(39), real(casapool%psoillab,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(40), real(casapool%psoilsorb,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(41), real(casapool%psoilocc,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(42), real(casabal%sumpbal,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(43), real(casaflux%plabuptake,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(44), real(casaflux%pdep,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(45), real(casaflux%pwea,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(46), real(casaflux%pleach,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(47), real(casaflux%ploss,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(48), real(casaflux%pupland,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(49), real(casaflux%plittermin,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(50), real(casaflux%psmin,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(51), real(casaflux%psimm,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(52), real(casaflux%psnet,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(53), real(casaflux%fpleach,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(54), real(casaflux%kplab,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(55), real(casaflux%kpsorb,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(56), real(casaflux%kpocc,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(57), real(casaflux%kmlabp,sp), cnt, otmp3patch )
+       call put_casa_var_patch(file_id, vid1(58), real(casaflux%psorbmax,sp), cnt, otmp3patch )
     endif
     ! 4D vars (land,patch,mplant,t)
     ! C
diff --git a/core/blaze/blaze.F90 b/core/blaze/blaze.F90
index c80b9bd49..ce01a8878 100644
--- a/core/blaze/blaze.F90
+++ b/core/blaze/blaze.F90
@@ -20,6 +20,7 @@ MODULE BLAZE_MOD
    CHARACTER(len=8) :: BLAZE_TSTEP    = "annually"  ! Call frequency ("daily","monthly","annually")
    CHARACTER(len=6) :: SIMFIRE_REGION = "GLOBAL"    ! either GLOBAL, EUROPE, ANZ
    CHARACTER(len=7) :: BURNT_AREA_SRC = "SIMFIRE"   ! either SIMFIRE or NONE !CLN for now!
+   INTEGER                              :: IAM ! number of master/worker for diagnostic output reasons
 END TYPE TYPE_BLAZE
 
 TYPE TYPE_TURNOVER
@@ -29,6 +30,7 @@ MODULE BLAZE_MOD
 REAL, DIMENSION(:,:), ALLOCATABLE  :: BLAZEFLX ! To BLAZE!!!
 
 INTEGER, PARAMETER :: NTO     = 7 ! Number of TurnOver Parameters ,i.e. #lines below
+INTEGER, PARAMETER :: NFLUXES = 14! Number of FLUXES in BLAZE output
 INTEGER, PARAMETER :: LEAF    = 1
 INTEGER, PARAMETER :: WOOD    = 2
 INTEGER, PARAMETER :: FROOT   = 3
@@ -142,7 +144,7 @@ SUBROUTINE INI_BLAZE ( np, LAT, LON, BLAZE)
   ALLOCATE ( BLAZE%POP_TO  ( np ) )
   ALLOCATE ( BLAZE%POP_CWD ( np ) )
   ALLOCATE ( BLAZE%POP_STR ( np ) )
-  ALLOCATE ( BLAZE%FLUXES  ( np, 13 ) )
+  ALLOCATE ( BLAZE%FLUXES  ( np, NFLUXES ) )
 
   ALLOCATE( BLAZE%AGLit_g(np,NPOOLS),BLAZE%AGLit_w(np,NPOOLS) )
   ALLOCATE( BLAZE%BGLit_g(np,NPOOLS),BLAZE%BGLit_w(np,NPOOLS) )
@@ -175,6 +177,7 @@ SUBROUTINE INI_BLAZE ( np, LAT, LON, BLAZE)
   BLAZE%AB        = 0.
   BLAZE%DEADWOOD  = 0.
   BLAZE%FSTEP     = 'DAILY'
+  !NESP2pt9 has this as =TRIM(blazeTStep)
 
   BLAZE%LAT       = LAT
   BLAZE%LON       = LON
@@ -374,7 +377,7 @@ SUBROUTINE BLAZE_TURNOVER(AB, CPLANT_g, CPLANT_w, AGL_g, AGL_w, &
   REAL,               INTENT(IN)    :: AB, shootfrac
   REAL,               INTENT(INOUT) :: CPLANT_w(3) , CPLANT_g(3)
   REAL, DIMENSION(3), INTENT(INOUT) :: AGL_w, AGL_g, BGL_w, BGL_g
-  REAL,               INTENT(OUT)   :: BT(13)
+  REAL,               INTENT(OUT)   :: BT(NFLUXES)
   REAL,     OPTIONAL, INTENT(IN)    :: POP_TO
   TYPE(TYPE_TURNOVER)               :: MTO(7)
   REAL    :: fAB
@@ -396,22 +399,23 @@ SUBROUTINE BLAZE_TURNOVER(AB, CPLANT_g, CPLANT_w, AGL_g, AGL_w, &
   MTO(SLIT )%TO_ATM = fAB * TO(SLIT )%TO_ATM * AGL_w(STR  )
   MTO(CLIT )%TO_ATM = fAB * TO(CLIT )%TO_ATM * AGL_w(CWD  )
 
-  ! Diagnostics
+  ! Diagnostics (new order acc to blaze_driver.F90 ll 215ff)
   BT( 1) = MTO(LEAF )%TO_ATM
-  BT( 2) = MTO(WOOD )%TO_ATM
-  BT( 3) = MTO(LEAF )%TO_STR
-  BT( 4) = MTO(WOOD )%TO_STR
-  BT( 5) = MTO(WOOD )%TO_CWD
-  BT( 6) = MTO(MLIT )%TO_ATM
-  BT( 7) = MTO(SLIT )%TO_ATM
-  BT( 8) = MTO(CLIT )%TO_ATM
-  BT( 9) = MTO(FROOT)%TO_ATM
-  BT(10) = MTO(FROOT)%TO_STR
-  BT(11) = AB * CPLANT_g (LEAF)
-  BT(12) = AB * AGL_g(METB)
-  BT(13) = AB * AGL_g(STR )
+  BT( 3) = MTO(WOOD )%TO_ATM
+  BT( 7) = MTO(LEAF )%TO_STR
+  BT( 9) = MTO(WOOD )%TO_STR
+  BT(10) = MTO(WOOD )%TO_CWD
 
+  BT( 4) = MTO(MLIT )%TO_ATM
+  BT( 5) = MTO(SLIT )%TO_ATM
+  BT( 6) = MTO(CLIT )%TO_ATM
+  BT( 2) = MTO(FROOT)%TO_ATM
+  BT( 8) = MTO(FROOT)%TO_STR
 
+  BT(11) = AB * CPLANT_g (LEAF)
+  BT(12) = AB * CPLANT_g (FROOT)
+  BT(13) = AB * AGL_g(METB)
+  BT(14) = AB * AGL_g(STR )
 
 !CLN  PRINT*,'BT( 1) =  MTO(LEAF )%TO_ATM    ',MTO(LEAF)%TO_ATM
 !CLN  PRINT*,'BT( 2) =  MTO(WOOD )%TO_ATM    ',MTO(WOOD)%TO_ATM
@@ -802,7 +806,7 @@ SUBROUTINE RUN_BLAZE(BLAZE, SF, CPLANT_g, CPLANT_w, tstp, YYYY, doy, TO , climat
      ! CALL SIMFIRE DAILY FOR ACOUNTING OF PARAMETERS
      CALL SIMFIRE ( SF, RAINF, TMAX, TMIN, DOY,MM, YYYY, BLAZE%AB , climate)
 
-PRINT*,"BLAZE doy,AB",doy,BLAZE%AB
+!PRINT*,"BLAZE doy,AB",doy,BLAZE%AB
      DO np = 1, BLAZE%NCELLS
         IF ( AVAIL_FUEL(1, CPLANT_w(np,:), CPLANT_g(np,:),BLAZE%AGLit_w(np,:),BLAZE%AGLit_g(np,:) ) .LE. MIN_FUEL ) &
              BLAZE%AB(np) = 0.
@@ -814,24 +818,22 @@ SUBROUTINE RUN_BLAZE(BLAZE, SF, CPLANT_g, CPLANT_w, tstp, YYYY, doy, TO , climat
      STOP -1
   ENDIF
 
+  !CALL PRINT_BLAZE(BLAZE)
+
   ! Apply half of former deadwood to atm now How to distribut (str
   ! set following Fraver 2013 pinus rosinosa (hardwood/decid. wood to be added
 
   DO np = 1, BLAZE%NCELLS
 
-
-
      CALL COMBUST( BLAZE, np, CPLANT_g(np,:), CPLANT_w(np,:),TO(np,:),BLAZE%AB(np).GT.0 )
 
-
-
      BLAZE%DFLI(np) = BLAZE%FLI(np)
      BLAZE%TO(np,:) = 0.
 
      IF (BLAZE%AB(np) .GT. 0. ) THEN
         CALL BLAZE_TURNOVER( BLAZE%AB(np), CPLANT_g(np,:), CPLANT_w(np,:), &
-             AGL_g(np,:), AGL_w(np,:),BGL_g(np,:), BGL_w(np,:), &
-             BLAZE%shootfrac(np), TO(np,:), BLAZE%FLUXES(np,:), BLAZE%BURNMODE )
+            BLAZE%AGLit_g(np,:), BLAZE%AGLit_w(np,:),BLAZE%BGLit_g(np,:), BLAZE%BGLit_w(np,:), &
+            BLAZE%shootfrac(np), TO(np,:), BLAZE%FLUXES(np,:), BLAZE%BURNMODE)
      ENDIF
 
 
@@ -866,7 +868,7 @@ SUBROUTINE WRITE_BLAZE_OUTPUT_NC ( BLAZE, FINAL )
     !INTEGER, INTENT(IN)    :: ctime
 
     INTEGER   :: STATUS
-    INTEGER   :: land_ID, t_ID
+    INTEGER   :: land_ID, t_ID, f_ID
     INTEGER   :: i, mp
     CHARACTER :: FNAME*99,dum*50
     LOGICAL, SAVE :: CALL1 = .TRUE.
@@ -877,8 +879,10 @@ SUBROUTINE WRITE_BLAZE_OUTPUT_NC ( BLAZE, FINAL )
     CHARACTER(len=20),DIMENSION(25):: A1
     ! 2 dim integer arrays (mp,t)
     CHARACTER(len=20),DIMENSION(1):: AI1
+    ! 3 dim real arrays (mp,t,nfluxes)
+    CHARACTER(len=20),DIMENSION(1):: A2
 
-    INTEGER, SAVE :: VIDtime, VID0(SIZE(A0)),VID1(SIZE(A1)),VIDI1(SIZE(AI1))
+    INTEGER, SAVE :: VIDtime, VID0(SIZE(A0)),VID1(SIZE(A1)),VIDI1(SIZE(AI1)),VID2(SIZE(A2))
     INTEGER, SAVE :: FILE_ID, CNT = 0
 
     A0(1) = 'latitude'
@@ -913,6 +917,8 @@ SUBROUTINE WRITE_BLAZE_OUTPUT_NC ( BLAZE, FINAL )
 
     AI1(1) = 'DaysSinceLastRain'
 
+    A2(1)  = 'BLAZE_fluxes'
+
     CNT = CNT + 1
 
     mp = BLAZE%ncells
@@ -949,6 +955,9 @@ SUBROUTINE WRITE_BLAZE_OUTPUT_NC ( BLAZE, FINAL )
        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)
+
+       STATUS = NF90_def_dim(FILE_ID, 'fluxes' , NFLUXES, f_ID)
        IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
 
         ! Define variables
@@ -982,6 +991,12 @@ SUBROUTINE WRITE_BLAZE_OUTPUT_NC ( BLAZE, FINAL )
            write(*,*) 'def AI1'
         END DO
 
+        DO i = 1, SIZE(A2)
+            STATUS = NF90_def_var(FILE_ID,TRIM(A2(i)) ,NF90_FLOAT,(/land_ID,f_ID ,t_ID/),VID2(i))
+            IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
+            write(*,*) 'def A2'
+        END DO
+
        ! End define mode:
        STATUS = NF90_enddef(FILE_ID)
        IF (STATUS /= NF90_noerr) CALL handle_err(STATUS)
@@ -1092,6 +1107,12 @@ SUBROUTINE WRITE_BLAZE_OUTPUT_NC ( BLAZE, FINAL )
     STATUS = NF90_PUT_VAR(FILE_ID, VIDI1(1), BLAZE%DSLR, start=(/ 1, CNT /), count=(/ mp, 1 /)  )
     IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)
 
+    ! Float 3D
+    DO i=1,NFLUXES 
+      STATUS = NF90_PUT_VAR(FILE_ID, VID2(1),BLAZE%FLUXES(:,i), start=(/ 1, i, CNT /), count=(/ mp, 1, 1 /))  
+      IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)
+    END DO
+
     IF(STATUS /= NF90_NoErr) CALL handle_err(STATUS)
 
     IF ( FINAL ) THEN
diff --git a/core/blaze/blaze_driver.F90 b/core/blaze/blaze_driver.F90
index 747396a8d..d40dcbb6d 100644
--- a/core/blaze/blaze_driver.F90
+++ b/core/blaze/blaze_driver.F90
@@ -37,9 +37,9 @@ SUBROUTINE BLAZE_DRIVER(NCELLS, BLAZE, SF, casapool,  casaflux, casamet, &
     TYPE(TYPE_TURNOVER)   ,ALLOCATABLE,SAVE :: TO(:,:)
     REAL,   DIMENSION(NCELLS) :: POP_TO
 
-    INTEGER       :: MM, DD, i, np, j, patch_index, p
+    INTEGER       :: MM, DD, i, np, j, patch_index, p, pidx
     REAL          :: TSTP
-    REAL          :: ag_litter_frac
+    REAL          :: ag_litter_frac, twto, rkill
     REAL          :: CPLANT_g (ncells,3),CPLANT_w (ncells,3)
     REAL          :: CLITTER_g(ncells,3),CLITTER_w(ncells,3)
     LOGICAL       :: EOY
@@ -169,20 +169,8 @@ SUBROUTINE BLAZE_DRIVER(NCELLS, BLAZE, SF, casapool,  casaflux, casamet, &
     DO i = 1, BLAZE%NCELLS
        DO p = 1, landpt(i)%nap  ! loop over number of active patches
           patch_index = landpt(i)%cstart + p - 1 ! patch index in CABLE vector
-!!$        IF ( casamet%lnonwood(patch_index) == 0 ) THEN ! Here woody patches
           veg%disturbance_intensity(patch_index,1) = BLAZE%AB(i) ! needed for ADJUST_POP_FOR_FIRE
           veg%disturbance_intensity(patch_index,2) = BLAZE%FLI(i) ! needed for ADJUST_POP_FOR_FIRE
-
-!!$           DO nh = 1,mheights
-!!$              hgt = real(casaflux%fire_mortality_vs_height(patch_index,nh,1))
-!!$              casaflux%fire_mortality_vs_height(patch_index,nh,2) = &
-!!$                   (1._r_2 - real(p_surv_OzSavanna(hgt,BLAZE%FLI(i)),r_2)) * BLAZE%AB(i) !*
-!!$
-!!$
-!!$              ! (1._r_2 - casaflux%fire_mortality_vs_height(patch_index,nh,2))  &
-!!$
-!!$           ENDDO
-!!$        ENDIF
        ENDDO
     ENDDO
 
@@ -191,32 +179,32 @@ SUBROUTINE BLAZE_DRIVER(NCELLS, BLAZE, SF, casapool,  casaflux, casamet, &
     ! POP cohorts and then interpolating fire_mortality between POP patches
     if (.NOT.Allocated(Iw)) allocate(Iw(POP%np))
     Iw = POP%Iwood
-    ! print*, 'IW', POP%Iwood
-    ! print*, 'fire_mort',  casaflux%fire_mortality_vs_height(Iw,:,:)
-    ! print*, 'BA: ', BLAZE%AB
-    ! print*, 'dist: ', int(veg%disturbance_interval(Iw,:), i4b)
 
     POP%pop_grid(:)%fire_mortality = 0.0_dp
-!!$  CALL ADJUST_POP_FOR_FIRE(pop,int(veg%disturbance_interval(Iw,:), i4b), &
-!!$       casaflux%fire_mortality_vs_height(Iw,:,:))
     CALL ADJUST_POP_FOR_FIRE(pop,int(veg%disturbance_interval(Iw,:), i4b), &
          veg%disturbance_intensity(Iw,1), veg%disturbance_intensity(Iw,2)  )
-    print*,"CLN ADJUST_POP_FOR_FIRE" ,int(veg%disturbance_interval(Iw,:), i4b), &
-         veg%disturbance_intensity(Iw,1), veg%disturbance_intensity(Iw,2)
-    casaflux%kplant_fire(Iw,WOOD) = max(min(POP%pop_grid(:)%fire_mortality/POP%pop_grid(:)%cmass_sum, &
+   
+   !this line is not in NESP2pt9_BLAZE but it looks sensible
+   casaflux%kplant_fire(Iw,WOOD) = max(min(POP%pop_grid(:)%fire_mortality/POP%pop_grid(:)%cmass_sum, &
          0.99_dp),0.0_dp)
 
-    !write(556,*) casaflux%kplant_fire(1,WOOD),  casaflux%kplant(1,WOOD), BLAZE%AB
-    !write(557,*) BLAZE%AB, casaflux%fire_mortality_vs_height(1,:,2)
-
-    !print*,  'fire mort', POP%pop_grid(:)%fire_mortality, casaflux%kplant_fire(Iw,WOOD)
-
-    DO i = 1, BLAZE%NCELLS
+   ! Apply turn-overs to biomass killed by fire in POP
+   BLAZE%FLUXES(:,:) = 0.
+   ! pop_grid index
+   pidx = 1
+   DO i = 1, BLAZE%NCELLS
        DO p = 1, landpt(i)%nap  ! loop over number of active patches
           patch_index = landpt(i)%cstart + p - 1 ! patch index in CABLE vector
 
           IF ( casamet%lnonwood(patch_index) == 1 ) THEN ! Here non-wood
 
+            !THe following section is in NESP2pt9 but commented out
+            !CLN wrong for kplant. Take out cmass here
+            !CLNcasaflux%kplant_fire(patch_index,LEAF)  = real(BLAZE%AB(i), r_2) &
+            !CLN     * real(casapool%cplant(patch_index,LEAF )*patch(patch_index)%frac, r_2)
+            !CLNcasaflux%kplant_fire(patch_index,FROOT) = real(BLAZE%AB(i), r_2) &
+            !CLN     * real(casapool%cplant(patch_index,FROOT)*patch(patch_index)%frac, r_2)
+
              casaflux%kplant_fire(patch_index,LEAF)  = real(BLAZE%AB(i), r_2) !CLN ???
              casaflux%kplant_fire(patch_index,FROOT) =  real(BLAZE%AB(i), r_2) !CLN ???
              casaflux%kplant_fire(patch_index,WOOD)  = 0.0_r_2
@@ -237,38 +225,107 @@ SUBROUTINE BLAZE_DRIVER(NCELLS, BLAZE, SF, casapool,  casaflux, casamet, &
              casaflux%fromPtoL_fire(patch_index,CWD,FROOT) = 0.0_r_2
              casaflux%fromPtoL_fire(patch_index,CWD,WOOD)  = 0.0_r_2
 
-          ELSEIF ( casamet%lnonwood(patch_index) == 0 ) THEN ! Here woody patches
-             casaflux%kplant_fire(patch_index,LEAF) =  real(BLAZE%AB(i) * &
-                  (TO(i, LEAF )%TO_ATM + TO(i, LEAF )%TO_STR), r_2)
-             casaflux%kplant_fire(patch_index,FROOT) = real(BLAZE%AB(i) * &
-                  (TO(i, FROOT )%TO_ATM + TO(i, FROOT )%TO_STR), r_2)
-
-
-             casaflux%klitter_fire(patch_index,METB) = real(BLAZE%AB(i) * &
-                  (TO(i, MLIT )%TO_ATM) * ag_litter_frac, r_2)
-             casaflux%klitter_fire(patch_index,STR) =  real(BLAZE%AB(i) * &
-                  (TO(i, SLIT )%TO_ATM) * ag_litter_frac, r_2)
-             casaflux%klitter_fire(patch_index,CWD) =  real(BLAZE%AB(i) * &
-                  (TO(i, CLIT )%TO_ATM) * ag_litter_frac, r_2)
-
-             casaflux%fromPtoL_fire(patch_index,STR,LEAF) = real(TO(i, LEAF )%TO_STR/ &
-                  MAX((TO(i, LEAF )%TO_STR + TO(i, LEAF )%TO_ATM),1.e-5), r_2)
-
-             casaflux%fromPtoL_fire(patch_index,STR,FROOT) = real(TO(i, FROOT )%TO_STR/ &
-                  MAX((TO(i, FROOT )%TO_STR + TO(i, FROOT )%TO_ATM),1.e-5), r_2)
-
-             ! maintain look-up table values for partitioning flux from woody biomass, but assume
-             ! flux to atm only applies to above-ground wood (70% of total)
-
-             casaflux%fromPtoL_fire(patch_index,CWD,WOOD) = real(TO(i, WOOD )%TO_CWD/ &
-                  MAX((TO(i, WOOD )%TO_CWD + TO(i, WOOD )%TO_ATM*0.7 +  TO(i, WOOD )%TO_STR ),1.e-5), r_2)
-
-             casaflux%fromPtoL_fire(patch_index,STR,WOOD) = real(TO(i, WOOD )%TO_STR/ &
-                  MAX((TO(i, WOOD )%TO_CWD + TO(i, WOOD )%TO_ATM + &
-                  TO(i, WOOD )%TO_STR ),1e-5), r_2)
-
+             ! BLAZE fluxes
+            BLAZE%FLUXES(i,11) = BLAZE%FLUXES(i,11) + casaflux%kplant_fire(patch_index,LEAF ) &
+               * real(casapool%cplant(patch_index,LEAF )*patch(patch_index)%frac, r_2)
+            BLAZE%FLUXES(i,12) = BLAZE%FLUXES(i,12) + casaflux%kplant_fire(patch_index,FROOT) &
+               * real(casapool%cplant(patch_index,FROOT)*patch(patch_index)%frac, r_2)
 
+            BLAZE%FLUXES(i,13) = BLAZE%FLUXES(i,13) + casaflux%klitter_fire(patch_index,METB) &
+               * real(casapool%clitter(patch_index,METB)*patch(patch_index)%frac, r_2)
+            BLAZE%FLUXES(i,14) = BLAZE%FLUXES(i,14) + casaflux%klitter_fire(patch_index,STR ) &
+               * real(casapool%clitter(patch_index,STR )*patch(patch_index)%frac, r_2)
 
+          ELSEIF ( casamet%lnonwood(patch_index) == 0 ) THEN ! Here woody patches
+            !LArge differences between CABLE-POP_TRENDY and NESP2pt9
+            !this is copied straight from NESP
+            
+            !CLN increment iwood for pop_grid
+            rkill=POP%pop_grid(pidx)%rkill
+            pidx = pidx + 1
+            
+            ! Check if there is mortality and COMBUST has only computed non-woody TO
+            ! When POP is involved these fluxes need to sum up to 1, assuming that
+            ! all that is not going to ATM or STR will be going to CWD (DEADWOOD)
+            !IF (CALL_POP) THEN
+            TO(i, WOOD)%TO_CWD = 1. - TO(i, WOOD)%TO_ATM - TO(i, WOOD)%TO_STR
+            TO(i, FROOT)%TO_STR   = MAX( TO(i, WOOD)%TO_ATM + TO(i, WOOD)%TO_STR + TO(i, WOOD)%TO_CWD - &
+                 TO(i, FROOT)%TO_ATM , 0. )
+            ! ENDIF
+ 
+            ! Total wood turn-over
+            twto = MAX(TO(i, WOOD)%TO_ATM * 0.7 + TO(i, WOOD )%TO_CWD + TO(i, WOOD )%TO_STR,1.e-7)
+ 
+            casaflux%kplant_fire(patch_index,LEAF)  = real(BLAZE%AB(i) * TO(i, LEAF )%TO_ATM, r_2)
+            casaflux%kplant_fire(patch_index,FROOT) = real(BLAZE%AB(i) * (1.-rkill) * TO(i, FROOT)%TO_ATM, r_2)
+            casaflux%kplant_fire(patch_index,WOOD ) = real(rkill/twto  * TO(i, WOOD )%TO_ATM * 0.7, r_2)
+ 
+            casaflux%klitter_fire(patch_index,METB) = real(BLAZE%AB(i) * TO(i, MLIT )%TO_ATM &
+                 * ag_litter_frac, r_2)
+            casaflux%klitter_fire(patch_index,STR)  = real(BLAZE%AB(i) * TO(i, SLIT )%TO_ATM &
+                 * ag_litter_frac, r_2)
+            casaflux%klitter_fire(patch_index,CWD)  = real(BLAZE%AB(i) * TO(i, CLIT )%TO_ATM &
+                 * ag_litter_frac, r_2)
+ 
+            casaflux%fromPtoL_fire(patch_index,STR,LEAF) = real(BLAZE%AB(i) * TO(i, LEAF )%TO_STR, r_2)
+            casaflux%fromPtoL_fire(patch_index,STR,FROOT)= real(rkill       * TO(i, FROOT)%TO_STR, r_2)
+            casaflux%fromPtoL_fire(patch_index,STR,WOOD) = real(rkill/twto  * TO(i, WOOD )%TO_STR, r_2)
+ 
+            casaflux%fromPtoL_fire(patch_index,CWD,WOOD) = real(rkill/twto  * TO(i, WOOD )%TO_CWD, r_2)
+ 
+            ! BLAZE fluxes 
+            BLAZE%FLUXES(i, 1) = BLAZE%FLUXES(i, 1) + casaflux%kplant_fire(patch_index,LEAF ) &
+                 * real(casapool%cplant(patch_index,LEAF )*patch(patch_index)%frac, r_2)
+            BLAZE%FLUXES(i, 2) = BLAZE%FLUXES(i, 2) + casaflux%kplant_fire(patch_index,FROOT) &
+                 * real(casapool%cplant(patch_index,FROOT)*patch(patch_index)%frac, r_2)
+            BLAZE%FLUXES(i, 3) = BLAZE%FLUXES(i, 3) + casaflux%kplant_fire(patch_index,WOOD ) &
+                 * real(casapool%cplant(patch_index,WOOD )*patch(patch_index)%frac, r_2) 
+ 
+            BLAZE%FLUXES(i, 4) = BLAZE%FLUXES(i, 4) + casaflux%klitter_fire(patch_index,METB) &
+                 * real(casapool%clitter(patch_index,METB)*patch(patch_index)%frac, r_2)
+            BLAZE%FLUXES(i, 5) = BLAZE%FLUXES(i, 5) + casaflux%klitter_fire(patch_index,STR ) &
+                 * real(casapool%clitter(patch_index,STR )*patch(patch_index)%frac, r_2)
+            BLAZE%FLUXES(i, 6) = BLAZE%FLUXES(i, 6) + casaflux%klitter_fire(patch_index,CWD ) &
+                 * real(casapool%clitter(patch_index,CWD )*patch(patch_index)%frac, r_2)
+ 
+            BLAZE%FLUXES(i, 7) = BLAZE%FLUXES(i, 7) + casaflux%fromPtoL_fire(patch_index,STR,LEAF )  &
+                 * real(casapool%clitter(patch_index,LEAF )*patch(patch_index)%frac, r_2)
+            BLAZE%FLUXES(i, 8) = BLAZE%FLUXES(i, 8) + casaflux%fromPtoL_fire(patch_index,STR,FROOT)  &
+                 * real(casapool%clitter(patch_index,FROOT)*patch(patch_index)%frac, r_2)
+            BLAZE%FLUXES(i, 9) = BLAZE%FLUXES(i, 9) + casaflux%fromPtoL_fire(patch_index,STR,WOOD )  &
+                 * real(casapool%clitter(patch_index,WOOD )*patch(patch_index)%frac, r_2)
+            BLAZE%FLUXES(i,10) = BLAZE%FLUXES(i,10) + casaflux%fromPtoL_fire(patch_index,CWD,WOOD )  &
+                 * real(casapool%clitter(patch_index,WOOD)*patch(patch_index)%frac, r_2)
+
+            !old code
+!             casaflux%kplant_fire(patch_index,LEAF) =  real(BLAZE%AB(i) * &
+!                  (TO(i, LEAF )%TO_ATM + TO(i, LEAF )%TO_STR), r_2)
+!             casaflux%kplant_fire(patch_index,FROOT) = real(BLAZE%AB(i) * &
+!                  (TO(i, FROOT )%TO_ATM + TO(i, FROOT )%TO_STR), r_2)
+!
+!             casaflux%klitter_fire(patch_index,METB) = real(BLAZE%AB(i) * &
+!                  (TO(i, MLIT )%TO_ATM) * ag_litter_frac, r_2)
+!             casaflux%klitter_fire(patch_index,STR) =  real(BLAZE%AB(i) * &
+!                  (TO(i, SLIT )%TO_ATM) * ag_litter_frac, r_2)
+!             casaflux%klitter_fire(patch_index,CWD) =  real(BLAZE%AB(i) * &
+!                  (TO(i, CLIT )%TO_ATM) * ag_litter_frac, r_2)
+!
+!             casaflux%fromPtoL_fire(patch_index,STR,LEAF) = real(TO(i, LEAF )%TO_STR/ &
+!                  MAX((TO(i, LEAF )%TO_STR + TO(i, LEAF )%TO_ATM),1.e-5), r_2)
+!
+!             casaflux%fromPtoL_fire(patch_index,STR,FROOT) = real(TO(i, FROOT )%TO_STR/ &
+!                  MAX((TO(i, FROOT )%TO_STR + TO(i, FROOT )%TO_ATM),1.e-5), r_2)
+!
+!             ! maintain look-up table values for partitioning flux from woody biomass, but assume
+!             ! flux to atm only applies to above-ground wood (70% of total)
+!
+!             casaflux%fromPtoL_fire(patch_index,CWD,WOOD) = real(TO(i, WOOD )%TO_CWD/ &
+!                 MAX((TO(i, WOOD )%TO_CWD + TO(i, WOOD )%TO_ATM*0.7 +  TO(i, WOOD )%TO_STR ),1.e-5), r_2)
+!
+!             casaflux%fromPtoL_fire(patch_index,STR,WOOD) = real(TO(i, WOOD )%TO_STR/ &
+!                  MAX((TO(i, WOOD )%TO_CWD + TO(i, WOOD )%TO_ATM + &
+!                  TO(i, WOOD )%TO_STR ),1e-5), r_2)
+!
           ENDIF
 
        ENDDO ! number of active patches
diff --git a/core/blaze/simfire.F90 b/core/blaze/simfire.F90
index 6bc751bbb..ef414eca4 100644
--- a/core/blaze/simfire.F90
+++ b/core/blaze/simfire.F90
@@ -77,7 +77,7 @@ SUBROUTINE INI_SIMFIRE( NCELLS, SF, modis_igbp )
 
   ! SF%RES    = 1./12.
 
-  ! SF%RES = 0.5
+  SF%RES = 0.5       ! CLN should stay 0.5 deg becasue SIMFIRE was trained on it. 
   SF%NCELLS = NCELLS
 
   ALLOCATE( SF%IGBP        (NCELLS) )
@@ -133,8 +133,8 @@ SUBROUTINE INI_SIMFIRE( NCELLS, SF, modis_igbp )
         SF%BIOME(i) = IGBP2BIOME(SF%IGBP(i),1)
      ENDIF
 
-     WRITE(*,FMT='(A5,I2,A17,I1,(1X,F7.2))')"IGBP ",SF%IGBP(i)," -> SIMFIREBIOME ", &
-          SF%BIOME(i),SF%LAT(i), SF%LON(i)
+     !WRITE(*,FMT='(A5,I2,A17,I1,(1X,F7.2))')"IGBP ",SF%IGBP(i)," -> SIMFIREBIOME ", &
+     !     SF%BIOME(i),SF%LAT(i), SF%LON(i)
 
      IF ( TRIM(SIMFIRE_REGION) == "GLOBAL" ) THEN ! GLOBAL
         SF%REGION(i) = 1
@@ -275,12 +275,13 @@ SUBROUTINE GET_POPDENS ( SF, YEAR )
   ELSE
      ISTEP = 5
   END IF
-  SF%RES    = HYRES
+
+  !CLN WROOONGGG!!!! SF%RES    = HYRES  !see note about calibration of SIMFIRE
   IF ( CALL1 ) THEN
      RF = NINT(SF%RES/HYRES)
      ! Check for Res being an integral multiple of 5' [RES] = fract. deg
      IF ( REAL(RF) .NE. SF%RES/HYRES .OR. SF%RES .LT. HYRES ) THEN
-        WRITE(*,*) 'Spatial resolution must be integral multiple of HYDE res. '
+        WRITE(*,*) 'Spatial resolution must be integer multiple of HYDE res. '
         WRITE(*,*) "RES:",SF%RES,"/ HYDE:",HYRES," = ",SF%RES/HYRES
         STOP "get_popdens in simfire_mod.f90"
      END IF
@@ -320,7 +321,7 @@ SUBROUTINE GET_POPDENS ( SF, YEAR )
         Y(i) = INT( (SF%LAT(i) +  90.) / SF%RES ) + 1
      END DO
 
-     write(*,*) 'X,Y:', X, Y,SF%RES
+     !write(*,*) 'X,Y:', X, Y,SF%RES
      NREAD = 2
      CALL1 = .FALSE.
   ELSE IF ( YEAR .GE. FINAL_YEAR ) THEN
@@ -336,7 +337,7 @@ SUBROUTINE GET_POPDENS ( SF, YEAR )
      NREAD = 0
   END IF
 
-  SF%RES = 1./12.
+  !CLN SF%RES = 1./12.  !see note about calibration of SIMFIRE
   RF = NINT(SF%RES/HYRES)
 
 
@@ -377,7 +378,7 @@ SUBROUTINE GET_POPDENS ( SF, YEAR )
            READ(iu,*)
         END DO
         ! Read data backwards ( llcorner is -180E, 90N )
-        write(*,*) NLAT, NLON
+        !write(*,*) NLAT, NLON
         DO j = NLAT, 1, -1
            READ(iu,*)(RVAL(i,j),i=1,NLON)
           ! write(3334,"(4320e16.6)") (RVAL(i,j),i=1,NLON)
@@ -400,7 +401,7 @@ SUBROUTINE GET_POPDENS ( SF, YEAR )
                  wPOPD = wPOPD + RVAL(ix,jy) * LAND_AREA(ix,jy)
                  wTOT  = wTOT  + LAND_AREA(ix,jy)
 
-                 write(*,*) 'RVAL: ',   RVAL(ix,jy), ix, jy
+                 !write(*,*) 'RVAL: ',   RVAL(ix,jy), ix, jy
               END DO
            END DO
 
diff --git a/offline/cable_driver.F90 b/offline/cable_driver.F90
index 0af2233a4..8c6b0b11a 100644
--- a/offline/cable_driver.F90
+++ b/offline/cable_driver.F90
@@ -1070,7 +1070,7 @@ PROGRAM cable_offline_driver
 
                  if (liseod) then ! end of day
                     if (cable_user%CALL_BLAZE) then
-                       print*, "CLN CAlling BLAZE"
+                       !print*, "CLN CAlling BLAZE"
                        call BLAZE_ACCOUNTING(BLAZE, climate, ktau, dels, YYYY, idoy)
                        call blaze_driver(blaze%ncells, blaze, simfire, casapool, casaflux, &
                             casamet, climate, rshootfrac, idoy, YYYY, 1, POP, veg)
diff --git a/offline/cable_iovars.F90 b/offline/cable_iovars.F90
index 6f14161f1..8976165ee 100644
--- a/offline/cable_iovars.F90
+++ b/offline/cable_iovars.F90
@@ -313,6 +313,7 @@ MODULE cable_IO_vars_module
          !! vh_js !! additional casa variables
 
          NBP = .FALSE., &
+         NEP = .FALSE., &
          dCdt = .FALSE., &
          TotSoilCarb = .FALSE.,   &
          TotLivBiomass = .FALSE., &
diff --git a/offline/cable_output.F90 b/offline/cable_output.F90
index 336fc5285..dfca2e126 100644
--- a/offline/cable_output.F90
+++ b/offline/cable_output.F90
@@ -70,7 +70,7 @@ MODULE cable_output_module
           snm, snmsl, tran, albs, albsn, cw, lqsn, lwsnl, mrfsofr, &
           mrlqso, mrlsl, snc, snd, snw, snwc, tcs, tgs, ts, tsl, &
           tsn, tsns, &
-          NBP, TotSoilCarb, TotLivBiomass, &
+          NBP, NEP, TotSoilCarb, TotLivBiomass, &
           TotLittCarb, SoilCarbFast, SoilCarbSlow, SoilCarbPassive, &
           LittCarbMetabolic, LittCarbStructural, LittCarbCWD, &
           PlantCarbLeaf, PlantCarbFineRoot, PlantCarbWood, &
@@ -202,6 +202,7 @@ MODULE cable_output_module
 
      ! [umol/m2/s]
      REAL(KIND=r_1), POINTER, DIMENSION(:) :: NBP => null()
+     REAL(KIND=r_1), POINTER, DIMENSION(:) :: NEP => null()
      REAL(KIND=r_1), POINTER, DIMENSION(:) :: dCdt => null()
      ! [kg C /m2]
      REAL(KIND=r_1), POINTER, DIMENSION(:) :: TotSoilCarb => null()
@@ -1238,13 +1239,22 @@ SUBROUTINE open_output_file(dels, soil, veg, bgc, rough)
 
     IF(output%carbon .OR. output%NBP) THEN
        CALL define_ovar(ncid_out, ovid%NBP, 'NBP', 'umol/m^2/s', &
-            'Net Biosphere Production (uptake +ve)' &
+            'Net Biosphere Production (uptake +ve, , includes fire and harvest fluxes)' &
             , patchout%NBP, &
             'dummy', xID, yID, zID, landID, patchID, tID)
        ALLOCATE(out%NBP(mp))
        out%NBP = zero4 ! initialise
     ENDIF
 
+    IF(output%carbon .OR. output%NEP) THEN
+       CALL define_ovar(ncid_out, ovid%NEP, 'NEP', 'umol/m^2/s', &
+            'Net Ecosystem Production (uptake +ve)'           &
+            , patchout%NEP, &
+            'dummy', xID, yID, zID, landID, patchID, tID)
+       ALLOCATE(out%NEP(mp))
+       out%NEP = zero4 ! initialise
+    ENDIF
+
     IF (output%casa) THEN
 
        CALL define_ovar(ncid_out, ovid%dCdt, 'dCdt', 'umol/m^2/s', &
@@ -3283,18 +3293,25 @@ SUBROUTINE write_output(dels, ktau, met, canopy, casaflux, casapool, casamet, ss
     END IF
 
     ! NBP and turnover fluxes [umol/m^2/s]
+    ! use variable NEP if wish to access NBP without the disturbance fluxes
     IF (output%carbon .OR. output%NBP) THEN
        ! Add current timestep's value to total of temporary output variable:
        IF (cable_user%POPLUC) THEN
           out%NBP = out%NBP + (-toreal4((casaflux%Crsoil-casaflux%cnpp &
-               - casapool%dClabiledt) * gd2umols)) !- &
-          !REAL((casaflux%FluxCtohwp + casaflux%FluxCtoclear  )/86400.0_r_2 &
-          !/ 1.201E-5_r_2, 4)
+               - casapool%dClabiledt +                                 &
+               REAL(casaflux%FluxCtohwp + casaflux%FluxCtoclear))      &
+               * gd2umols))
        ELSE
           out%NBP = out%NBP + (-toreal4((casaflux%Crsoil-casaflux%cnpp &
                - casapool%dClabiledt) * gd2umols))
        ENDIF
 
+       !add on fire fluxes iff BLAZE active
+       IF (cable_user%call_blaze) THEN
+          out%NBP = out%NBP - toreal4((casaflux%fluxCtoCO2_plant_fire          &
+               + casaflux%fluxCtoCO2_litter_fire) * gd2umols)
+       ENDIF
+
        IF (writenow) THEN
           ! Divide accumulated variable by number of accumulated time steps:
           out%NBP = out%NBP * rinterval
@@ -3305,6 +3322,24 @@ SUBROUTINE write_output(dels, ktau, met, canopy, casaflux, casapool, casamet, ss
           out%NBP = zero4
        END IF
     ENDIF
+
+    IF (output%carbon .OR. output%NEP) THEN
+          ! Add current timestep's value to total of temporary output variable:
+          ! NOte NEP doesn't include fire or harvest fluxes   
+          out%NEP = out%NEP + (-toreal4((casaflux%Crsoil-casaflux%cnpp &
+               - casapool%dClabiledt) * gd2umols))
+     
+          IF (writenow) THEN
+               ! Divide accumulated variable by number of accumulated time steps:
+               out%NEP = out%NEP * rinterval
+               ! Write value to file:
+               CALL write_ovar(out_timestep, ncid_out, ovid%NEP, 'NEP', out%NEP,    &
+                    ranges%NEE, patchout%NEP, 'default', met)
+               ! Reset temporary output variable:
+               out%NEP = zero4
+          END IF
+     ENDIF
+
     IF (output%casa) THEN
        ! Add current timestep's value to total of temporary output variable:
        out%dCdt = out%dCdt + toreal4((casapool%ctot-casapool%ctot_0) * gd2umols)
diff --git a/scripts/run_cable-pop-actblaze.sh b/scripts/run_cable-pop-actblaze.sh
new file mode 100755
index 000000000..05194636c
--- /dev/null
+++ b/scripts/run_cable-pop-actblaze.sh
@@ -0,0 +1,1941 @@
+#!/usr/bin/env bash
+
+# Explor / Pearcey (launch with:  sbatch --ignore-pbs)
+# https://slurm.schedmd.com/sbatch.html
+# Name - 8 letters and digits
+#SBATCH -J x0003-s0
+#SBATCH -o %x-%j.out
+#SBATCH -e %x-%j.out
+# Explor partitions (sinfo): std (2x16, parallel), sky (2x16, parallel, AVX512), hf (2x4, serial),
+#                            P100 (2x16, GPU), GTX (2x16, GPU), ivy (2x8, parallel), k20 (2x8, GPU)
+#SBATCH -p std
+# -N Nodes / -n tasks (mpiexec, srun, ...) / -c cpus_per_task (OpenMP, make-j, ...)
+#SBATCH -N 1
+#SBATCH -n 1
+# #SBATCH --ntasks-per-node=32
+# Check memory on *nix with /usr/bin/time -v ./prog
+# time (day-hh:mm:ss) / memory (optional, units K,M,G,T)
+#SBATCH -t 03:59:59
+#SBATCH --mem=100G
+# notify: Valid type values are NONE,BEGIN,END,FAIL,REQUEUE,ALL,STAGE_OUT,TIME_LIMIT,TIME_LIMIT_90/80/50,ARRAY_TASKS
+#SBATCH --mail-type=FAIL,STAGE_OUT,TIME_LIMIT
+#SBATCH --mail-user=matthias.cuntz@inrae.fr
+
+# Gadi
+# https://opus.nci.org.au/display/Help/4.+PBS+Jobs
+#PBS -N x0001-s1
+#PBS -P x45
+# express / normal / copyq (2x24, cascadelake)
+# expressbw / normalbw (2x14, broadwell) / normalsl (2x16, skylake)- ex-Raijin nodes
+#PBS -q normal
+# Global run or Australia at 0.25: 192 GB memory and 48 cpus, ca. 12 hours walltime
+#PBS -l walltime=01:00:00
+#PBS -l mem=48GB
+#PBS -l ncpus=1
+# #PBS -l jobfs=1GB
+#PBS -l storage=gdata/x45+gdata/rp23
+#PBS -l software=netCDF:MPI:Intel:GNU
+#PBS -r y
+#PBS -l wd
+#PBS -j oe
+#PBS -S /bin/bash
+#PBS -M ian.harman@csiro.au
+#PBS -m ae
+
+# cuntz@explor, cuntz@mc16, cuntz@mcinra, moc801@gadi cuntz@gadi
+# kna016@pearcey knauer@pearcey, jk8585@gadi knauer@gadi
+# bri220@pearcey pcb599@gadi, briggs@gadi
+# vil029@percey yc3714@gadi villalobos@gadi
+# nieradzik@aurora
+system=inh599@gadi
+
+# MPI run or single processor run
+# nproc should fit with job tasks
+dompi=0   # 0: normal run: ./cable
+          # 1: MPI run: mpiexec -n ${nproc} ./cable_mpi
+nproc=1   # Number of cores for MPI runs
+          # must be same as above: SBATCH -n nproc or PBS -l ncpus=nproc
+
+# --------------------------------------------------------------------
+#
+# Full Cable run with biomass spinup, POP, land-use change, etc.
+#
+# This script uses CRU-JRA forcing.
+#
+# Global meteo and land-use change data can be used with a mask giving land points.
+# In step 0, the land mask can be extracted for a single point, an area
+# or a number of random points can be chosen.
+# Alternatively, single site met, LUH2, forcing and mask can be extracted from the global data sets.
+#
+# The run sequence is as follows:
+#   1. Create a climate restart file using Cable's default vegetation distribution.
+#   2. First phase of spinup with static land use, fixed atmospheric CO2 and
+#      N deposition from 1700, and 30 years of repeated meteorology.
+#      Zero initial biomass stocks, POP and climate.
+#   3. Bring biomass stocks into equilibrium, restricting labile P and mineral N pools.
+#      Repeat a and b several times.
+#      a) Start from restart files.
+#      b) Use dump files for the biophysics. Quasi-equilibrium of soil and litter pools
+#         using an analytic steady-state solution.
+#   4. Same as 3 but without any restriction on labile P and mineral N pools.
+#      Repeat a and b several times.
+#   5. Second phase of spinup with dynamic land use, atmospheric CO2 and N deposition.
+#      a) Dynamic land use from 1580 to 1699, using still fixed atmospheric CO2 and
+#         N deposition from 1700, and 30 years of repeated meteorology.
+#      b) Run from 1700 to 1899 with dynmic land use, varying atmospheric CO2 and N deposition,
+#         but still with 30 years of repeated meteorology.
+#   6. Final historical run, everything dynamic from 1900 to 2017.
+#   7. Future run, everything dynamic (not all met types).
+#
+# Written,  Matthias Cuntz, Aug 2019, following the run scripts and namelists provided by Vanessa Haverd
+# Modified, Jurgen Knauer, 2020      - gm_explicit, coordination, acclimation
+#                                    - bios, plume, future runs
+#           Matthias Cuntz, Mar 2021 - functions into run_cable-pop_lib.sh
+#
+# --------------------------------------------------------------------
+
+#ASKJK - changes in comparison to gm_acclim_coord
+# 1. Removed switch Test_parallel=1 -> used [[ dompi -eq 1 ]]
+# 2. ScriptsPath="${workpath}/scripts"
+#    but
+#    namelistpath="$(dirname ${workpath})/namelists"
+#    with
+#      workpath="/home/599/jk8585/CABLE_run/gm_acclim_coord/global_runs"
+#      cablehome="/home/599/jk8585/CABLE_code"
+#    -> changed to similar of ScriptsPath="$(dirname ${workpath})/scripts"
+# 3. Need plume.nml, bios.nml, gm_LUT_*.nc
+# 4. output%grid = "mask" (in gm_acclim_coord) or "land" (default before)
+# 5. What should be for Run in bios.nml after 1. Climate restart?
+# 6. Why is YearEnd different for plume (1849) compared to cru (1699) in 5a. First dynamic land use?
+# 7. Do we need chunking in 5a, 5b, 6, and 7?
+# 8. Do we need cropping output to latlon region at the end: is this not in step 0 with ${doextractsite} -eq 1?
+#ASKJK - changes in comparison to gm_acclim_coord
+
+# --------------------------------------------------------------------
+# Sequence switches
+#
+
+imeteo=0        # 0: Use global meteo, land use and mask
+                # 1: Use local mask, but global meteo and land use (doextractsite=1)
+                # 2: Use local meteo, land use and mask (doextractsite=2)
+# Step 0
+purge_restart=0 # 1/0: Do/Do not delete all restart files (completely new run, e.g. if settings changed)
+doextractsite=0 # 0: Do not extract local meteo, land use nor mask
+                # 1: Do extract only mask at specific site/region (imeteo=1)
+                # 2: Do extract meteo, land use and mask at specific site/region (imeteo=2)
+                #    Does not work with randompoints /= 0 but with latlon
+    experiment=alledits_blazeoff
+    randompoints=0   # <0: use -1*randompoints from file ${LandMaskFilePath}/${experiment}_points.csv if existing
+                     # 0:  use latlon
+                     # >0: generate and use randompoints random grid points from GlobalLandMaskFile
+    # lat,lon  or  latmin,latmax,lonmin,lonmax   # must have . in numbers otherwise indexes taken
+    latlon=42.536875,-72.172602
+    # latlon=-34.5,-33.5,149.5,156.5
+    # latlon=42.5,43.5,109.5,110.5
+    # latlon=-44.0,-10.0,110.0,155.0  # Australia
+# Step 1
+doclimate=1     # 1/0: Do/Do not create climate restart file
+# Step 2
+dofromzero=1    # 1/0: Do/Do not first spinup phase from zero biomass stocks
+# Step 3
+doequi1=1       # 1/0: Do/Do not bring biomass stocks into quasi-equilibrium with restricted P and N pools
+nequi1=1        #      number of times to repeat steps in doequi1
+# Step 4
+doequi2=0       # 1/0: Do/Do not bring biomass stocks into quasi-equilibrium with unrestricted P and N pools
+nequi2=2        #      number of times to repeat steps in doequi2
+# Step 5a
+doiniluc=0      # 1/0: Do/Do not spinup with dynamic land use (5a)
+# Step 5b
+doinidyn=0      # 1/0: Do/Do not full dynamic spinup from 1700 to 1899 (5b)
+# Step 6
+dofinal=0       # 1/0: Do/Do not final run from 1900 to 2017
+# Step 7
+dofuture=0      # 1/0: Do/Do not future runs (plume only)
+
+# --------------------------------------------------------------------
+# Other switches
+#
+
+# MetType
+mettype="bios"       # "cru", "plume", "bios"
+degrees=0.05        #  bios only: resolution of met and luc (0.05 or 0.25)
+metmodel="hadgem2"  # "hadgem2", "ipsl" (only used if mettype="plume")
+RCP="hist"          # "hist", "2.6", "4.5", "6.0", "8.5" (no future runs if RCP="hist")
+
+# Cable
+# Accounting for mesophyll conductance
+# 1: explicit (finite)
+# 0: implicit
+explicit_gm=1
+# Parameter conversion accounting for gm (only used if explicit_gm=1)
+# 1: Use lookup table for parameter conversion accounting for gm
+# 0: Do not use lookup table
+use_LUTgm=1
+# Parameters for RubisCO temperature function
+# "Bernacchi_2002" or "Walker_2013"
+Rubisco_params="Bernacchi_2002"
+# Coordination of photosynthesis
+# 1: Coordinate photosynthetic rates
+# 0: Do not coordinate
+coordinate_photosyn=1
+# Coordination type (only used if coordinate_photosyn=1)
+# T: force coordinated RubisCO- and energy-limited rates
+# F: optimise photosynthesis
+coord=F
+# Acclimation of photosynthesis
+# 1: Acclimate photosynthesis
+# 0: Do not acclimate
+acclimate_photosyn=1
+# POP
+# 1: Use POP, population dynamics model coupled to CASA
+# 0: Do not use POP
+call_pop=1
+# Soil and snow model
+# "default" or "sli"
+soil_struc="default"
+# 13C
+# 1: Calculate 13C
+# 0: no 13C
+doc13o2=0
+# 13C leaf discrimination
+# 1: use simple 13C leaf discrimination (a+(b-a)*Ci/Ca)
+# 0: calculate full 13C leaf discrimination
+c13o2_simple_disc=0
+
+# --------------------------------------------------------------------
+# Setup
+#
+
+set -e
+
+trap cleanup 1 2 3 6
+
+pid=$$
+isdir="${PWD}"
+prog=$0
+pprog=$(basename ${prog})
+pdir=$(dirname ${prog})
+tmp=${TMPDIR:-"/tmp"}
+
+system=$(echo ${system} | tr A-Z a-z)
+sys=${system#*@}
+user=${system%@*}
+
+#
+# Special things on specific computer system such as loading modules
+#
+
+export mpiexecdir=
+if [[ "${sys}" == "explor" ]] ; then
+    # prog is slurm_script
+    pdir=${isdir}
+    # # INTELMPI - load mpi module first, otherwise intel module will not prepend LD_LIBRARY_PATH
+    # module load intelmpi/2018.5.274
+    # module load intel/2018.5
+    # unset I_MPI_PMI_LIBRARY
+    # export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${HOME}/local/lib:${HOME}/local/netcdf-fortran-4.4.4-ifort2018.0/lib
+    # export mpiexecdir=/soft/env/soft/all/intel/2018.3/compilers_and_libraries_2018.5.274/linux/mpi/intel64/bin
+    # # INTEL / OpenMPI - load mpi module first, otherwise intel module will not pre-pend LD_LIBRARY_PATH
+    # module load openmpi/3.0.0/intel18
+    # module load intel/2018.5
+    # export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${HOME}/local/lib:${HOME}/local/netcdf-fortran-4.4.4-ifort2018.0/lib
+    # export mpiexecdir=/opt/soft/hf/openmpi-3.0.0-intel18/bin
+    # GNU / OpenMPI
+    module load gcc/6.3.0
+    module load openmpi/3.0.1/gcc/6.3.0
+    export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${HOME}/local/lib:${HOME}/local/netcdf-fortran-4.4.4-gfortran63/lib
+    export mpiexecdir=/opt/soft/hf/openmpi/3.0.1/gcc/6.3.0/bin
+    if [[ ${doextractsite} -ge 1 ]] ; then module load python/intel/2019/3 ; fi
+elif [[ "${sys}" == "mc16" ]] ; then
+    export mpiexecdir=/usr/local/openmpi-4.1.6-gfortran/bin
+    # export mpiexecdir=/usr/local/openmpi-4.1.1-ifort/bin
+elif [[ "${sys}" == "mcinra" ]] ; then
+    export mpiexecdir=/usr/local/openmpi-3.1.4-gfortran/bin
+    # export mpiexecdir=/usr/local/openmpi-3.1.5-ifort/bin
+elif [[ "${sys}" == "pearcey" ]] ; then
+    # prog is slurm_script
+    pdir=${isdir}
+    module del intel-cc intel-fc
+    module add intel-cc/16.0.1.150 intel-fc/16.0.1.150
+    module unload intel-mpi/5.0.1.035
+    module add netcdf/4.3.3.1 openmpi/1.8.8
+elif [[ "${sys}" == "gadi" ]] ; then
+    pdir=${isdir}
+    . /etc/bashrc
+    module purge
+    #module load intel-compiler/2019.5.281
+    #module load intel-mpi/2019.5.281
+    #module load netcdf/4.6.3
+    module load intel-compiler/2021.5.0
+    module load intel-mpi/2021.5.1
+    module load netcdf/4.8.0
+    # module load hdf5/1.10.5
+    if [[ ${doextractsite} -ge 1 ]] ; then
+        module load python3/3.7.4
+        export PYTHONPATH=${PYTHONPATH}:/g/data/x45/python/lib/python3.7/site-packages
+    fi
+    if [[ ${randompoints} -eq 0 ]] ; then module load nco/4.9.2 ; fi  # needed for cropping outputs
+    export mpiexecdir=/apps/intel-mpi/2019.5.281/intel64/bin
+elif [[ "${sys}" == "aurora" ]] ; then
+    pdir=${isdir}
+    module purge
+    module load intel/2020a
+    module load netCDF-Fortran/4.5.2
+    if [[ ${doextractsite} -ge 1 ]] ; then
+        module load GCCcore/10.3.0 Python/3.9.5
+        export PYTHONPATH=${PYTHONPATH}  #CLN not copied!#:/g/data/x45/python/lib/python3.7/site-packages
+	      echo 'Extraction not possible at the moment. Set Py-Path!'
+	      exit -666
+    fi
+    export mpiexecdir=""
+fi
+if [[ ! -z ${mpiexecdir} ]] ; then export mpiexecdir="${mpiexecdir}/" ; fi
+
+#
+# Directories of things
+#
+
+# Naming of experiments for gm_acclim_coord experiments of Jurgen
+if [[ ("${system}" == "jk8585@gadi") || ("${system}" == "knauer@gadi") ]] ; then
+    experiment=test_${mettype}
+    if [[ "${mettype}" == "plume" ]] ; then
+        # (TEST_)(exp|imp)_(no)coord(T|F)(_PSacclim)_(B2002|W2013)
+        if [[ ${explicit_gm} -eq 1 ]] ; then
+            if [[ ${coordinate_photosyn} -eq 1 ]] ; then
+                if [[ ${acclimate_photosyn} -eq 1 ]] ; then
+                    experiment=exp_coord${coord}_PSacclim
+                else
+	                  experiment=exp_coord${coord}
+                fi
+            else
+	              if [[ ${acclimate_photosyn} -eq 1 ]] ; then
+	                  experiment=exp_nocoord_PSacclim
+	              else
+	                  experiment=exp_nocoord
+	              fi
+	          fi
+        else
+	          if [[ ${coordinate_photosyn} -eq 1 ]] ; then
+	              if [[ ${acclimate_photosyn} -eq 1 ]] ; then
+		                experiment=imp_coord${coord}_PSacclim
+	              else
+		                experiment=imp_coord${coord}
+	              fi
+	          else
+	              if [[ ${acclimate_photosyn} -eq 1 ]] ; then
+		                experiment=imp_nocoord_PSacclim
+	              else
+		                experiment=imp_nocoord
+	              fi
+            fi
+        fi
+
+        if [[ "${Rubisco_params}" == "Bernacchi_2002" ]] ; then
+            experiment=${experiment}_B2002
+        elif [[ "${Rubisco_params}" == "Walker_2013" ]] ; then
+            experiment=${experiment}_W2013
+        else
+	          echo "unknown option for 'Rubisco_params', exiting..."
+	          exit 1
+        fi
+
+        if [[ ${dompi} -ne 0 ]] ; then
+	          experiment=TEST_${experiment}
+        fi
+    fi
+fi
+
+#
+# Relative directories must be relative to the directory of this script,
+#   not relative to the directory from which this script is launched (if different)
+#   nor relative to the run path.
+#
+#if [[ "${system}" == "cuntz@explor" ]] ; then
+#    # Run directory: runpath="${sitepath}/run_xxx"
+#    cablebase="/home/oqx29/zzy20/prog/cable"
+#    sitepath="${cablebase}/runs/single_sites/${experiment}"
+#    cablehome="${cablebase}/branches/CABLE-POP_TRENDY"
+#    workpath=${cablehome}
+#    # Cable executable
+#    if [[ ${dompi} -eq 1 ]] ; then
+#        exe="${cablehome}/offline/cable-mpi"
+#    else
+#        exe="${cablehome}/offline/cable"
+#    fi
+#    # CABLE-AUX directory (uses offline/gridinfo_CSIRO_1x1.nc and offline/modis_phenology_csiro.txt)
+#    aux="${cablebase}/CABLE-AUX"
+#    # Global Mask - for create_landmask.py
+#    GlobalLandMaskFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"
+#    SurfaceFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"   # note that SurfaceFile does not need subsetting
+#    # # Global Mask - for global run
+#    # GlobalLandMaskFile="/home/oqx29/zzy20/data/crujra/daily_1deg/glob_ipsl_1x1.nc"
+#    # Global CRU
+#    GlobalMetPath="/home/oqx29/zzy20/data/crujra/daily_1deg"
+#    MetVersion="CRUJRA_2018"
+#    # Global LUC
+#    GlobalTransitionFilePath="/home/oqx29/zzy20/data/LUH2_v3_1deg/"
+#elif [[ "${system}" == "cuntz@mc16" || "${system}" == "cuntz@mcinra" ]] ; then
+#    # Run directory: runpath="${sitepath}/run_xxx"
+#    cablebase="/Users/cuntz/prog/vanessa/cable"
+#    sitepath="${cablebase}/runs/single_sites/${experiment}"
+#    cablehome="${cablebase}/branches/CABLE-POP_TRENDY"
+#    workpath=${cablehome}
+#    # Cable executable
+#    if [[ ${dompi} -eq 1 ]] ; then
+#        # exe="${cablehome}/offline/cable-mpi-gfortran"
+#        exe="${cablehome}/offline/cable-mpi-ifort"
+#    else
+#        exe="${cablehome}/offline/cable"
+#    fi
+#    # CABLE-AUX directory (uses offline/gridinfo_CSIRO_1x1.nc and offline/modis_phenology_csiro.txt)
+#    aux="${cablebase}/CABLE-AUX"
+#    # Global Mask, CRU, LUC
+#    GlobalLandMaskFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"
+#    SurfaceFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"   # note that SurfaceFile does not need subsetting
+#    GlobalMetPath=
+#    MetVersion="CRUJRA_2018"
+#    GlobalTransitionFilePath=
+#elif [[ "${system}" == "moc801@gadi" || "${system}" == "cuntz@gadi" ]] ; then
+#    # Run directory: runpath="${sitepath}/run_xxx"
+#    # sitepath="/home/801/moc801/prog/cable/runs/single_sites/${experiment}"
+#    sitepath="/scratch/x45/moc801/cable/c13"
+#    cablehome="/home/801/moc801/prog/cable/branches/CABLE-POP_TRENDY"
+#    workpath=${cablehome}
+#    # Cable executable
+#    if [[ ${dompi} -eq 1 ]] ; then
+#        exe="${cablehome}/offline/cable-mpi"
+#    else
+#        exe="${cablehome}/offline/cable"
+#    fi
+#    # CABLE-AUX directory (uses offline/gridinfo_CSIRO_1x1.nc and offline/modis_phenology_csiro.txt)
+#    aux="/g/data/x45/CABLE-AUX"
+#    # BLAZE directory
+#    BlazeDataPath="/g/data/x45/Data_BLAZE"
+#    # Global MET, MASK, LUH
+#    if [[ "${mettype}" == "cru" ]] ; then
+#        # Global MASK
+#        SurfaceFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"   # note that SurfaceFile does not need subsetting
+#        # Global MET
+#        # GlobalLandMaskFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"
+#        GlobalLandMaskFile="/g/data/x45/ipbes/masks/glob_ipsl_1x1.nc"
+#        MetVersion="CRUJRA_2021"
+#        GlobalMetPath="/g/data/x45/CRUJRA2021/daily_1deg"
+#        # Global LUC
+#        GlobalTransitionFilePath="/g/data/x45/LUH2/GCB_2019/1deg/EXTRACT"
+#    elif [[ "${mettype}" == "plume" ]] ; then
+#        # Global MASK
+#        SurfaceFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"   # note that SurfaceFile does not need subsetting
+#        # Global MET
+#	      # GlobalLandMaskFile="/g/data/x45/ipbes/masks/glob_ipsl_1x1.nc"
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/gridinfo_CSIRO_1x1.nc"
+#        MetVersion=
+#	      GlobalMetPath="/g/data/x45/ipbes/${metmodel}/1deg"
+#        # only in plume.nml
+#        CO2Path="/g/data/x45/ipbes/co2"
+#	      NdepPath="/g/data/x45/ipbes/ndep"
+#        # Global LUC
+#        GlobalTransitionFilePath="/g/data/x45/LUH2/GCB_2019/1deg/EXTRACT"
+#    elif [[ "${mettype}" == "bios" ]] ; then
+#        # Global MASK
+#        SurfaceFile="${aux}/offline/gridinfo_CSIRO_CRU05x05_4tiles.nc"   # note that SurfaceFile does not need subsetting
+#        # Global MET
+#        GlobalLandMaskFile="/g/data/x45/BIOS3_forcing/acttest9/acttest9"  # no file extension
+#        MetVersion=
+#        GlobalMetPath="/g/data/x45/BIOS3_forcing/acttest9/met/"  # last slash is needed
+#	      ParamPath="/g/data/x45/BIOS3_forcing/acttest9/params/"  # only in bios.nml
+#        # Global LUC
+#        GlobalTransitionFilePath="/g/data/x45/LUH2/v3h/${degrees}deg_aust/EXTRACT"
+#    fi
+#elif [[ "${system}" == "kna016@pearcey" || "${system}" == "knauer@pearcey" ]] ; then
+#    # Run directory: runpath="${sitepath}/run_xxx"
+#    sitepath="/OSM/CBR/OA_GLOBALCABLE/work/Juergen/CABLE_run/parallel_runs/${experiment}"
+#    cablehome="/OSM/CBR/OA_GLOBALCABLE/work/Juergen/CABLE_code/CABLE-POP_TRENDY"
+#    workpath=${cablehome}
+#    # Cable executable
+#    if [[ ${dompi} -eq 1 ]] ; then
+#        exe="${cablehome}/offline/cable-mpi"
+#    else
+#        exe="${cablehome}/offline/cable"
+#    fi
+#    # CABLE-AUX directory (uses offline/gridinfo_CSIRO_1x1.nc and offline/modis_phenology_csiro.txt)
+#    aux="/OSM/CBR/OA_GLOBALCABLE/work/Vanessa/CABLE-AUX"
+#    # Global Mask
+#    GlobalLandMaskFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"
+#    SurfaceFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"   # note that SurfaceFile does not need subsetting
+#    # Global CRU
+#    GlobalMetPath="/OSM/CBR/OA_GLOBALCABLE/work/CRU-JRA55/crujra/daily_1deg"
+#    MetVersion="CRUJRA_2021"
+#    # Global LUC
+#    GlobalTransitionFilePath="/OSM/CBR/OA_GLOBALCABLE/work/LUH2/v3/1deg"
+#elif [[ "${system}" == "jk8585@gadi" || "${system}" == "knauer@gadi" ]] ; then
+#    # Run directory: runpath="${sitepath}/run"
+#    sitepath="/g/data/x45/jk8585/global_runs/gm_acclim_coord/${experiment}"
+#    cablecode="/home/599/jk8585/CABLE_code/gm_acclim_coord"
+#    cablehome="/home/599/jk8585/CABLE_run/gm_acclim_coord"
+#    workpath=${cablehome}
+#    # Cable executable
+#    if [[ ${dompi} -eq 1 ]] ; then
+#        exe="${cablecode}/offline/cable-mpi"
+#    else
+#        exe="${cablecode}/offline/cable"
+#    fi
+#    # CABLE-AUX directory (uses offline/gridinfo_CSIRO_1x1.nc and offline/modis_phenology_csiro.txt)
+#    aux="/g/data/x45/CABLE-AUX"
+#    # Global Mask
+#    SurfaceFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"   # note that SurfaceFile does not need subsetting
+#    # Global Met
+#    if [[ "${mettype}" == "cru" ]] ; then
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/glob_ipsl_1x1.nc"
+#        GlobalMetPath="/g/data/x45/CRUJRA2020/daily_1deg"
+#        MetVersion="CRUJRA_2021"
+#    elif [[ "${mettype}" == "plume" ]] ; then
+#	      # GlobalLandMaskFile="/g/data/x45/ipbes/masks/glob_ipsl_1x1.nc"
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/gridinfo_CSIRO_1x1.nc"
+#	      GlobalMetPath="/g/data/x45/ipbes/${metmodel}/1deg"
+#        MetVersion=
+#        # only in plume.nml
+#        CO2Path="/g/data/x45/ipbes/co2"
+#	      NdepPath="/g/data/x45/ipbes/ndep"
+#    elif [[ "${mettype}" == "bios" ]] ; then
+#	      GlobalLandMaskFile="/g/data/x45/BIOS3_forcing/acttest9/acttest9" # no file extension
+#        GlobalMetPath="/g/data/x45/BIOS3_forcing/acttest9/met/"  # last slash is needed
+#        MetVersion=
+#        # only in bios.nml
+#	      ParamPath="/g/data/x45/BIOS3_forcing/acttest9/params/"
+#    fi
+#    # Global LUC
+#    GlobalTransitionFilePath="/g/data/x45/LUH2/GCB_2019/1deg/EXTRACT"
+#elif [[ "${system}" == "yc3714@gadi" || "${system}" == "villalobos@gadi" ]] ; then
+#    # Run directory: runpath="${sitepath}/run"
+#    sitepath="/g/data/x45/BIOS3_output/${experiment}" # Results 
+#    cablehome="/home/563/yc3714/CSIRO/CABLE_code/9011/NESP2pt9_BLAZE" # model home
+#    workpath="/home/563/yc3714/CSIRO/CABLE_run/BIOS3_blaze" # run directory
+#    # Cable executable
+#    if [[ ${dompi} -eq 1 ]] ; then
+#        exe="${cablehome}/offline/cable-mpi"
+#    else
+#        exe="${cablehome}/offline/cable"
+#    fi
+#    # CABLE-AUX directory (uses offline/gridinfo_CSIRO_1x1.nc and offline/modis_phenology_csiro.txt)
+#    aux="/g/data/x45/CABLE-AUX"
+#    # Global Mask
+#    SurfaceFile="${aux}/offline/gridinfo_CSIRO_CRU05x05_4tiles.nc"   # note that SurfaceFile does not need subsetting
+#    # Global Met
+#    if [[ "${mettype}" == "cru" ]] ; then
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/glob_ipsl_1x1.nc"
+#        MetVersion="CRUJRA_2021"
+#        GlobalMetPath="/g/data/x45/CRUJRA2021/daily_1deg"
+#    elif [[ "${mettype}" == "plume" ]] ; then
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/gridinfo_CSIRO_1x1.nc"
+#        MetVersion=
+#	      GlobalMetPath="/g/data/x45/ipbes/${metmodel}/1deg"
+#        # only in plume.nml
+#        CO2Path="/g/data/x45/ipbes/co2"
+#	      NdepPath="/g/data/x45/ipbes/ndep"
+#    elif [[ "${mettype}" == "bios" ]] ; then
+#    	  GlobalLandMaskFile="/g/data/x45/BIOS3_forcing/acttest9/acttest9" # no file extension
+#        MetVersion=
+#        GlobalMetPath="/g/data/x45/BIOS3_forcing/acttest9/met/"          # last slash is needed
+#	      ParamPath="/g/data/x45/BIOS3_forcing/acttest9/params/"           # only in bios.nml
+#        GlobalTransitionFilePath="/g/data/x45/LUH2/v3h/${degrees}deg_aust/EXTRACT"
+#    fi
+#    # Global LUC
+#    # GlobalTransitionFilePath="/g/data/x45/LUH2/GCB_2019/1deg/EXTRACT"
+#elif [[ "${system}" == "bri220@pearcey" || "${system}" == "briggs@pearcey" ]] ; then
+#    # Run directory: runpath="${sitepath}/run_xxx"
+#    sitepath="/OSM/CBR/OA_GLOBALCABLE/work/Peter/CABLE_run/parallel_runs/${experiment}"
+#    cablehome="/OSM/CBR/OA_GLOBALCABLE/work/Peter/CABLE_code"
+#    workpath=${cablehome}
+#    # Cable executable
+#    if [[ ${dompi} -eq 1 ]] ; then
+#        exe="${cablehome}/NESP2pt9_BLAZE/offline/cable-mpi"
+#    else
+#        exe="${cablehome}/NESP2pt9_BLAZE/offline/cable"
+#    fi
+#    # CABLE-AUX directory (uses offline/gridinfo_CSIRO_1x1.nc and offline/modis_phenology_csiro.txt)
+#    aux="/OSM/CBR/OA_GLOBALCABLE/work/Vanessa/CABLE-AUX"
+#    # Global Mask
+#    GlobalLandMaskFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"
+#    # Global CRU
+#    MetVersion="CRUJRA_2021"
+#    GlobalMetPath="/OSM/CBR/OA_GLOBALCABLE/work/CRU-JRA55/crujra/daily_1deg"
+#    # Global LUC
+#    GlobalTransitionFilePath="/OSM/CBR/OA_GLOBALCABLE/work/LUH2/v3/1deg"
+#elif [[ "${system}" == "pcb599@gadi" || "${system}" == "briggs@gadi" ]] ; then
+#    # Run directory: runpath="${sitepath}/run_xxx"
+#    sitepath="/g/data/x45/BIOS3_output/${experiment}"
+#    cablehome="/home/599/pcb599/CABLE_code/NESP2pt9_BLAZE"
+#    workpath="/home/599/pcb599/CABLE_run/BIOS3/run"
+#    # Cable executable
+#    if [[ ${dompi} -eq 1 ]] ; then
+#	      exe="${cablehome}/offline/cable-mpi"
+#    else
+#	      exe="${cablehome}/offline/cable"
+#    fi
+#    # CABLE-AUX directory (uses offline/gridinfo_CSIRO_1x1.nc and offline/modis_phenology_csiro.txt)
+#    aux="/g/data/x45/CABLE-AUX"
+#    # Global Mask
+#    GlobalLandMaskFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"
+#    SurfaceFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"   # note that SurfaceFile does not need subsetting
+#    # Global LUC
+#    GlobalTransitionFilePath="/g/data/x45/LUH2/GCB_2019/1deg/EXTRACT"
+#    # Global Met
+#    if [[ "${mettype}" == 'plume' ]] ; then
+#	      GlobalMetPath="/g/data/x45/ipbes/${metmodel}/1deg"
+#        CO2Path="/g/data/x45/ipbes/co2"
+#	      NdepPath="/g/data/x45/ipbes/ndep"
+#        MetVersion=
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/glob_ipsl_1x1.nc"
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/gridinfo_CSIRO_1x1.nc"
+#    elif [[ "${mettype}" == 'bios' ]] ; then
+#	      GlobalLandMaskFile="/g/data/x45/BIOS3_forcing/acttest9/acttest9" # no file extension
+#        GlobalMetPath="/g/data/x45/BIOS3_forcing/acttest9/met/"  # last slash is needed
+#        MetVersion=
+#	      ParamPath="/g/data/x45/BIOS3_forcing/acttest9/params/"
+#        GlobalTransitionFilePath="/g/data/x45/LUH2/v3h/${degrees}deg_aust/EXTRACT"
+#    elif [[ "${mettype}" == 'cru' ]] ; then
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/glob_ipsl_1x1.nc"
+#        MetVersion="CRUJRA_2021"
+#        GlobalMetPath="/g/data/x45/CRUJRA201/daily_1deg"
+#    fi
+#elif [[ "${system}" == "nieradzik@aurora" ]] ; then
+#    # Run directory: runpath="${sitepath}/run_xxx"
+#    sitepath="/home/x_larni/STOREDIR/RUNDIR/CABLE/${experiment}"
+#    cablehome="/home/x_larni/SRC/CABLE/NESP2pt9_BLAZE"
+#    workpath=${cablehome}
+#    # Cable executable
+#    if [[ ${dompi} -eq 1 ]] ; then
+#	      exe="${cablehome}/offline/cable-mpi"
+#    else
+#	      exe="${cablehome}/offline/cable"
+#    fi
+#    # CABLE-AUX directory (uses offline/gridinfo_CSIRO_1x1.nc and offline/modis_phenology_csiro.txt)
+#    aux="/home/x_larni/STOREDIR/DATA/CABLE-AUX"
+#    # Global Mask
+#    GlobalLandMaskFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"
+#    SurfaceFile="${aux}/offline/gridinfo_CSIRO_1x1.nc"   # note that SurfaceFile does not need subsetting
+#    # Global LUC
+#    GlobalTransitionFilePath="/g/data/x45/LUH2/GCB_2019/1deg/EXTRACT"
+#    # Global Met
+#    if [[ "${mettype}" == 'plume' ]] ; then
+#        MetVersion=
+#	      GlobalMetPath="/g/data/x45/ipbes/${metmodel}/1deg"
+#        CO2Path="/g/data/x45/ipbes/co2"
+#	      NdepPath="/g/data/x45/ipbes/ndep"
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/glob_ipsl_1x1.nc"
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/gridinfo_CSIRO_1x1.nc"
+#    elif [[ "${mettype}" == 'bios' ]] ; then
+#	      GlobalLandMaskFile="/home/x_larni/STOREDIR/DATA/CABLE_INPUT/acttest9/acttest9" # no file extension
+#        MetVersion=
+#        GlobalMetPath="/home/x_larni/STOREDIR/DATA/CABLE_INPUT/acttest9/met/"  # last slash is needed
+#	      ParamPath="/home/x_larni/STOREDIR/DATA/CABLE_INPUT/acttest9/params/"
+#        GlobalTransitionFilePath="/home/x_larni/STOREDIR/DATA/CABLE_INPUT/LUH2/v3h/${degrees}deg_aust/EXTRACT"
+#    elif [[ "${mettype}" == 'cru' ]] ; then
+#	      GlobalLandMaskFile="/g/data/x45/ipbes/masks/glob_ipsl_1x1.nc"
+#        MetVersion="CRUJRA_2021"
+#        GlobalMetPath="/g/data/x45/CRUJRA2021/daily_1deg"
+#    fi
+#else
+#    echo "System not known."
+#    exit 1
+#fi
+
+#this is the section specifically for inh599 on gadi
+sitepath="/scratch/x45/inh599/BIOStests/${experiment}" # Results
+workpath="/home/599/inh599/JAC/CABLE/newBLAZE" # run directory
+cablehome="/home/599/inh599/JAC/git/CABLE-POP-BLAZE/CABLE/"
+# Cable executable
+if [[ ${dompi} -eq 1 ]] ; then
+    exe="${cablehome}/offline/cable-mpi"   #not valid
+else
+    exe="${cablehome}/bin/cable"
+fi
+# CABLE-AUX directory (uses offline/gridinfo_CSIRO_1x1.nc and offline/modis_phenology_csiro.txt)
+aux="/g/data/x45/CABLE-AUX"
+BlazeDataPath="/g/data/x45/Data_BLAZE"
+# Global Mask
+SurfaceFile="${aux}/offline/gridinfo_CSIRO_CRU05x05_4tiles.nc"
+GlobalLandMaskFile="/g/data/x45/BIOS3_forcing/acttest9/acttest9" # no file extension
+GlobalMetPath="/g/data/x45/BIOS3_forcing/acttest9/met/"          # last slash is needed - updated 29/3/2024
+ParamPath="/g/data/x45/BIOS3_forcing/acttest9/params/"           # only in bios.nml
+#GlobalTransitionFilePath="/g/data/x45/LUH2/v3h/${degrees}deg_aust/EXTRACT"
+GlobalTransitionFilePath="/g/data/rp23/experiments/2024-07-01_TRENDYv13/input/luc/"
+# end of inh599 stuff
+
+# Run directory
+runpath="${sitepath}/run"
+
+# Cable parameters - inh599 edits here
+if [[ "${mettype}" == 'bios' ]] ; then
+    namelistpath="${workpath}/namelists_bios"
+    namelistcrupath="${workpath}/namelists_bios"
+    filename_veg="${cablehome}/params/v11/def_veg_params.txt"
+    filename_soil="${cablehome}/params/v11/def_soil_params.txt"
+    casafile_cnpbiome="${cablehome}/params/v11/pftlookup.csv"
+else
+    namelistpath="${workpath}/namelists"
+    namelistcrupath="${workpath}/namelists/v12"
+    filename_veg="${workpath}/params/def_veg_params.txt"
+    filename_soil="${workpath}/params/def_soil_params.txt"
+    casafile_cnpbiome="${workpath}/params/pftlookup.csv"
+fi
+# Other scripts
+ScriptsPath="${cablehome}/scripts"
+# Mask (not used unless imeteo > 0)
+LandMaskFile="${sitepath}/mask/${experiment}_landmask.nc"
+# Met (not used unless imeteo > 0)
+if [[ "${mettype}" == 'bios' ]] ; then
+    MetPath="${sitepath}/met/bios_${degrees}deg"
+    ClimateFile="${sitepath}/mask/bios_climate_rst.nc"
+    if [[ "${sys}" == "gadi" ]]; then
+        if [[ (${doclimate} -eq 0) && (! -f ${ClimateFile}) ]] ; then
+            ClimateFile="/g/data/x45/BIOS3_output/bio_climate_acttest9/bios_climate_rst.nc"
+        fi
+    fi
+else
+    MetPath="${sitepath}/met/cru_jra_1deg"
+    ClimateFile="${sitepath}/mask/cru_climate_rst.nc"
+    if [[ "${sys}" == "gadi" ]]; then
+        if [[ (${doclimate} -eq 0) && (! -f ${ClimateFile}) ]] ; then
+            ClimateFile="/g/data/x45/ipbes/cable_climate/ipsl_climate_rst_glob_1deg.nc"
+        fi
+    fi
+fi
+# LUC
+TransitionFilePath="${sitepath}/LUH2/v3/1deg"
+# gm lookup tables
+gm_lut_bernacchi_2002="${cablehome}/params/gm_LUT_351x3601x7_1pt8245_Bernacchi2002.nc"
+gm_lut_walker_2013="${cablehome}/params/gm_LUT_351x3601x7_1pt8245_Walker2013.nc"
+# 13C
+filename_d13c_atm="${cablehome}/params/graven_et_al_gmd_2017-table_s1-delta_13c-1700-2025.txt"
+
+# --------------------------------------------------------------------
+# Start Script
+# --------------------------------------------------------------------
+
+# --------------------------------------------------------------------
+# Helper functions, most functions are in plumber_cable-pop_lib.sh
+#
+
+source ${pdir}/run_cable-pop_lib.sh
+
+# usage of script
+function usage()
+{
+    printf "${pprog} [-h]\n"
+    printf "Runs Cable on a single grid cell with spinup, POP, land-use change, etc.\n"
+    printf "Behaviour of the script is controlled by switches at the top of the script (ca. line 101ff).\n"
+    printf "\n"
+    printf "Options\n"
+    printf "    -h    Prints this help screen.\n"
+}
+
+# cleanup at end or at trap
+function cleanup()
+{
+    \rm -f ${tmp}/*.${pid}*
+    exit 1
+}
+
+# returns argument to extract lat and lon with ncks
+function nckslatlon()
+{
+    vars=$(ncvarlist ${1})
+    if [[ -z $(isin latitude ${vars}) ]] ; then ilat="lat" ; else ilat="latitude" ; fi
+    if [[ -z $(isin longitude ${vars}) ]] ; then ilon="lon" ; else ilon="longitude" ; fi
+    if [[ -z $(echo ${2} | cut -f 3 -d ",") || -z $(echo ${2} | cut -f 4 -d ",") ]] ; then
+        iilat=$(echo ${2} | cut -f 1 -d ",")
+        iilon=$(echo ${2} | cut -f 2 -d ",")
+        echo "-d ${ilat},${iilat} -d ${ilon},${iilon}"
+    else
+        iilat1=$(echo ${2} | cut -f 1 -d ",")
+        iilat2=$(echo ${2} | cut -f 2 -d ",")
+        iilon1=$(echo ${2} | cut -f 3 -d ",")
+        iilon2=$(echo ${2} | cut -f 4 -d ",")
+        echo "-d ${ilat},${iilat1},${iilat2} -d ${ilon},${iilon1},${iilon2}"
+    fi
+}
+
+# --------------------------------------------------------------------------------------------------
+# Preparation
+#
+# Get options
+while getopts "h" option ; do
+    case ${option} in
+        h) usage; exit;;
+        *) printf "Error ${pprog}: unimplemented option.\n\n" 1>&2;  usage 1>&2; exit 1;;
+    esac
+done
+shift $((${OPTIND} - 1))
+
+#
+# get directories
+mkdir -p ${sitepath}/mask
+pdir=$(abspath ${pdir})
+cd ${pdir}
+adir=$(abspath ${aux})
+exe=$(absfile ${exe})
+mkdir -p ${runpath}
+rdir=$(abspath ${runpath})
+ndir=$(abspath ${namelistpath})
+ncdir=$(abspath ${namelistcrupath})
+sdir=$(abspath ${ScriptsPath})
+
+#
+# prepare run directory
+cd ${rdir}
+mkdir -p logs
+mkdir -p outputs
+mkdir -p restart
+ln -sf ${adir}
+# ln -sf ${exe}
+cp ${exe} ./
+iexe=$(basename ${exe})
+cd ${pdir}
+
+#
+# set stacksize to unlimited if permitted, otherwise to 15 bit if possible
+set +e
+ulimit -s unlimited 2> /dev/null || ulimit -s 32768
+set -e
+
+# --------------------------------------------------------------------
+# Info
+#
+t1=$(date +%s)
+printf "Started at %s\n" "$(date)"
+
+printf "\nSetup\n"
+printf "    Serial / Parallel\n"
+printf "        dompi=${dompi}\n"
+printf "            nproc=${nproc}\n"
+printf "\n"
+printf "    Sequence\n"
+printf "        imeteo=${imeteo}\n"
+printf "        doextractsite=${doextractsite}\n"
+printf "            experiment=${experiment}\n"
+printf "            randompoints=${randompoints}\n"
+printf "            latlon=${latlon}\n"
+printf "        doclimate=${doclimate}\n"
+printf "        dofromzero=${dofromzero}\n"
+printf "        doequi1=${doequi1}\n"
+printf "            nequi1=${nequi1}\n"
+printf "        doequi2=${doequi2}\n"
+printf "            nequi2=${nequi2}\n"
+printf "        doiniluc=${doiniluc}\n"
+printf "        doinidyn=${doinidyn}\n"
+printf "        dofinal=${dofinal}\n"
+printf "        dofuture=${dofuture}\n"
+printf "\n"
+printf "    Options\n"
+printf "        mettype=${mettype}\n"
+printf "        metmodel=${metmodel}\n"
+printf "        RCP=${RCP}\n"
+printf "        explicit_gm=${explicit_gm}\n"
+printf "        use_LUTgm=${use_LUTgm}\n"
+printf "        Rubisco_params=${Rubisco_params}\n"
+printf "        coordinate_photosyn=${coordinate_photosyn}\n"
+printf "        coord=${coord}\n"
+printf "        acclimate_photosyn=${acclimate_photosyn}\n"
+printf "        call_pop=${call_pop}\n"
+printf "        soil_struc=${soil_struc}\n"
+printf "        doc13o2=${doc13o2}\n"
+printf "        c13o2_simple_disc=${c13o2_simple_disc}\n"
+printf "\n"
+printf "    Directories\n"
+printf "        sitepath=${sitepath}\n"
+printf "        cablehome=${cablehome}\n"
+printf "        exe=${exe}\n"
+printf "        aux=${aux}\n"
+printf "        GlobalLandMaskFile=${GlobalLandMaskFile}\n"
+printf "        SurfaceFile=${SurfaceFile}\n"
+printf "        GlobalMetPath=${GlobalMetPath}\n"
+printf "        GlobalTransitionFilePath=${GlobalTransitionFilePath}\n"
+printf "        runpath=${runpath}\n"
+printf "        namelistpath=${namelistpath}\n"
+printf "        filename_veg=${filename_veg}\n"
+printf "        filename_soil=${filename_soil}\n"
+printf "        casafile_cnpbiome=${casafile_cnpbiome}\n"
+printf "        LandMaskFile=${LandMaskFile}\n"
+printf "        MetPath=${MetPath}\n"
+printf "        ClimateFile=${ClimateFile}\n"
+printf "        TransitionFilePath=${TransitionFilePath}\n"
+printf "        gm_lut_bernacchi_2002=${gm_lut_bernacchi_2002}\n"
+printf "        gm_lut_walker_2013=${gm_lut_walker_2013}\n"
+printf "        filename_d13c_atm=${filename_d13c_atm}\n"
+printf "\n"
+
+# --------------------------------------------------------------------
+# Prep input
+#
+
+# 0. Extract meteo, land use and mask for one specific site from global files
+if [[ ${doextractsite} -ge 2 ]] ; then
+    set +e
+    xnco=$(which ncks 2> /dev/null)
+    set -e
+    if [[ -z ${xnco} ]] ; then module load nco ; fi
+fi
+
+if [[ ${doextractsite} -eq 1 ]] ; then
+    echo "0. Set local mask"
+    cd ${pdir}
+    # mask
+    LandMaskFilePath=$(dirname ${LandMaskFile})
+    mkdir -p ${LandMaskFilePath}
+    LandMaskFile=$(absfile ${LandMaskFile})
+    # echo $(basename ${LandMaskFile})
+    # generate random points if ${randompoints} > 0
+    if [[ ${randompoints} -ne 0 ]] ; then
+        dogeneraterandom=1
+        rpoints=${randompoints}
+        if [[ ${randompoints} -lt 0 ]] ; then
+            if [[ -f ${LandMaskFilePath}/${experiment}_points.csv ]] ; then dogeneraterandom = 0 ; fi
+            rpoints=$(( ${randompoints} * -1 ))
+        fi
+        # generate random points
+        cat > ${tmp}/sedtmp.${pid} << EOF
+            basepath      = "${sitepath}"
+            gridinfo_file = "${GlobalLandMaskFile}"
+            outname       = "${LandMaskFilePath}/${experiment}_points.csv"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${sdir}/generate_latlonlist.py ${LandMaskFilePath}/generate_latlonlist.py
+        python3 ${LandMaskFilePath}/generate_latlonlist.py ${rpoints}
+
+        # set mask to generated random points
+        cat > ${tmp}/sedtmp.${pid} << EOF
+            path          = "${LandMaskFilePath}"
+            maskfname     = "${LandMaskFile}"
+            latlonfile    = "${LandMaskFilePath}/${experiment}_points.csv"
+            gridinfo_file = "${GlobalLandMaskFile}"
+EOF
+	      if [[ "${mettype}" == "bios" ]] ; then
+            echo "res = 0.25" >> ${tmp}/sedtmp.${pid}
+	      fi
+        applysed ${tmp}/sedtmp.${pid} ${sdir}/create_landmask.py ${LandMaskFilePath}/create_landmask.py
+        sed -i -e "s|from lnutils.*|sys.path.insert(1,'${sdir}'); from lnutils import latlon2ixjy|" ${LandMaskFilePath}/create_landmask.py
+        python3 ${LandMaskFilePath}/create_landmask.py
+    else  # no random points
+        # # cdo -s -f nc4 -z zip sellonlatbox,-72.5,-72.0,42.5,43.0 ${GlobalLandMaskFile} ${LandMaskFile}
+        # ncks -O $(nckslatlon ${GlobalLandMaskFile} ${latlon}) ${GlobalLandMaskFile} ${LandMaskFile}
+        cat > ${tmp}/sedtmp.${pid} << EOF
+            path          = "${LandMaskFilePath}"
+            maskfname     = "${LandMaskFile}"
+            gridinfo_file = "${GlobalLandMaskFile}"
+EOF
+	      if [[ "${mettype}" == "bios" ]] ; then
+            echo "res = 0.25" >> ${tmp}/sedtmp.${pid}
+	      fi
+        applysed ${tmp}/sedtmp.${pid} ${sdir}/create_landmask.py ${LandMaskFilePath}/create_landmask.py
+        sed -i -e "s|from lnutils.*|sys.path.insert(1,'${sdir}'); from lnutils import latlon2ixjy|" ${LandMaskFilePath}/create_landmask.py
+        python3 ${LandMaskFilePath}/create_landmask.py ${latlon}
+    fi
+elif [[ ${doextractsite} -eq 2 ]] ; then
+    echo "0. Extract local meteo and mask"
+    cd ${pdir}
+    # meteorology
+    met_list="pre pres dlwrf dswrf spfh tmax tmin ugrd vgrd ndep"
+    mkdir -p ${MetPath}
+    MetPath=$(abspath ${MetPath})
+    for mm in ${met_list} ; do
+        echo ${GlobalMetPath}/${mm}
+        mkdir -p ${MetPath}/${mm}
+        for nc in ${GlobalMetPath}/${mm}/*.nc ; do
+            ff=$(basename ${nc})
+            echo "    ${ff}"
+            # cdo -s -f nc4 -z zip sellonlatbox,-72.5,-72.0,42.5,43.0 ${nc} ${MetPath}/${mm}/${ff}
+            ncks -O $(nckslatlon ${nc} ${latlon}) ${nc} ${MetPath}/${mm}/${ff}
+        done
+    done
+    echo ${GlobalMetPath}/co2/*.csv
+    mkdir -p ${MetPath}/co2
+    cp ${GlobalMetPath}/co2/*.csv ${MetPath}/co2/
+
+    # land use
+    mkdir -p ${TransitionFilePath}
+    TransitionFilePath=$(abspath ${TransitionFilePath})
+    echo ${GlobalTransitionFilePath}
+    for nc in ${GlobalTransitionFilePath}/*.nc ; do
+        ff=$(basename ${nc})
+        echo "    ${ff}"
+        # cdo -s -f nc4 -z zip sellonlatbox,-72.5,-72.0,42.5,43.0 ${nc} ${TransitionFilePath}/${ff}
+        ncks -O $(nckslatlon ${nc} ${latlon}) ${nc} ${TransitionFilePath}/${ff}
+    done
+
+    # mask
+    LandMaskFilePath=$(dirname ${LandMaskFile})
+    mkdir -p ${LandMaskFilePath}
+    LandMaskFile=$(absfile ${LandMaskFile})
+    echo $(basename ${LandMaskFile})
+    # cdo -s -f nc4 -z zip sellonlatbox,-72.5,-72.0,42.5,43.0 ${GlobalLandMaskFile} ${LandMaskFile}
+    ncks -O $(nckslatlon ${GlobalLandMaskFile} ${latlon}) ${GlobalLandMaskFile} ${LandMaskFile}
+fi
+
+
+# --------------------------------------------------------------------
+# Prepare sequence
+#
+
+# Choose meteo, land use and mask directories and files
+if [[ ${imeteo} -eq 0 ]] ; then
+    MetPath=$(abspath ${GlobalMetPath})
+    TransitionFilePath=$(abspath ${GlobalTransitionFilePath})
+    LandMaskFile=$(absfile ${GlobalLandMaskFile})
+elif [[ ${imeteo} -eq 1 ]] ; then
+    MetPath=$(abspath ${GlobalMetPath})
+    TransitionFilePath=$(abspath ${GlobalTransitionFilePath})
+    LandMaskFile=$(absfile ${LandMaskFile})
+elif [[ ${imeteo} -eq 2 ]] ; then
+    MetPath=$(abspath ${MetPath})
+    TransitionFilePath=$(abspath ${TransitionFilePath})
+    LandMaskFile=$(absfile ${LandMaskFile})
+else
+    printf "Error ${pprog}: imeteo option unknown: ${imeteo}.\n\n"
+    exit 1
+fi
+
+
+# absolute pathes of other parameter files
+ClimateFile=$(absfile ${ClimateFile})
+filename_veg=$(absfile ${filename_veg})
+filename_soil=$(absfile ${filename_soil})
+casafile_cnpbiome=$(absfile ${casafile_cnpbiome})
+gm_lut_bernacchi_2002=$(absfile ${gm_lut_bernacchi_2002})
+gm_lut_walker_2013=$(absfile ${gm_lut_walker_2013})
+filename_d13c_atm=$(absfile ${filename_d13c_atm})
+if [[ "${Rubisco_params}" == "Bernacchi_2002" ]] ; then
+    filename_gm_lut=${gm_lut_bernacchi_2002}
+elif [[ "${Rubisco_params}" == "Walker_2013" ]] ; then
+    filename_gm_lut=${gm_lut_walker_2013}
+else
+    filename_gm_lut=""
+fi
+
+
+# delete all restart files if required
+if [[ ${purge_restart} -eq 1 ]] ; then
+    rm -f ${rdir}/restart/*
+    # rm -f ${ClimateFile}
+fi
+
+
+# Write standard namelists with options that are common to all steps of the sequence.
+# They can, however, be overwritten in later steps.
+
+# global meteo namelist file
+if [[ "${mettype}" == "cru" ]] ; then
+    cat > ${tmp}/sedtmp.${pid} << EOF
+        BasePath     = "${MetPath}"
+        MetPath      = "${MetPath}"
+        MetVersion   = "${MetVersion}"
+        LandMaskFile = "${LandMaskFile}"
+        Run          = "S0_TRENDY"
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${ncdir}/cru.nml ${rdir}/cru_${experiment}.nml
+elif [[ "${mettype}" == "plume" ]] ; then
+    #TODO: NDEPfile, CO2file
+    cat > ${tmp}/sedtmp.${pid} << EOF
+        BasePath     = "${MetPath}"
+        Forcing      = "hadgem2-es"
+    	  LandMaskFile = "${LandMaskFile}"
+    	  Run          = "spinup"
+    	  RCP          = "hist"
+    	  CO2          = "static1850"
+    	  CO2file      = "${CO2Path}/co2_1850_2005_hist.dat"
+    	  NDEP         = "static1850"
+        NDEPfile     = "${NdepPath}/NOy_plus_NHx_dry_plus_wet_deposition_hist_1850_2015_annual_1deg.nc"
+EOF
+    if [[ "${metmodel}" == "ipsl" ]] ; then
+        sed -i -e "/Forcing/s/=.*/= 'ipsl-cm5a-lr'/" ${tmp}/sedtmp.${pid}
+    fi
+    applysed ${tmp}/sedtmp.${pid} ${ndir}/plume.nml ${rdir}/plume_${experiment}.nml
+elif [[ "${mettype}" == "bios" ]] ; then
+    #using updated meteorology
+    cat > ${tmp}/sedtmp.${pid} << EOF
+        met_path         = "${MetPath}/"
+        param_path       = "${ParamPath}"
+        landmaskflt_file = "${GlobalLandMaskFile}.flt"
+        landmaskhdr_file = "${GlobalLandMaskFile}.hdr"
+        rain_file        = "1900010120231231_rain_recal_b2405.bin"
+        swdown_file      = "1900010120231231_rad_b2405.bin"
+        tairmax_file     = "1900010120231231_tmax_noclim_b2405.bin"
+        tairmin_file     = "1900010120231231_tmin_noclim_b2405.bin"
+        wind_file        = "1900010120231231_windspeed_ms_b2405.bin"
+        vp0900_file      = "1900010120231231_vph09_b2405.bin"
+        vp1500_file      = "1900010120231231_vph15_b2405.bin"
+        co2_file         = "1700_2023_trendy_global_co2_ann.bin"
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${ndir}/bios.nml ${rdir}/bios_${experiment}.nml
+fi
+
+# global landuse change namelist
+cat > ${tmp}/sedtmp.${pid} << EOF
+    TransitionFilePath = "${TransitionFilePath}"
+    ClimateFile        = "${ClimateFile}"
+    YearStart          = 1700
+    YearEnd            = 2017
+EOF
+applysed ${tmp}/sedtmp.${pid} ${ncdir}/luc.nml ${rdir}/luc_${experiment}.nml
+
+# Blaze namelist !CLN CHECK
+cat > ${tmp}/sedtmp.${pid} << EOF
+    blazeTStep       = "annually"  ! Call frequency ("daily", "monthly", "annually")
+    BurnedAreaSource = "SIMFIRE"   ! Burnt Area ("PRESCRIBED", "SIMFIRE", "GFED4")
+    BurnedAreaFile   = "${BlazeDataPath}/BA_Aust_2001-2019.nc"  ! used for Prescribed fires !CLN not available for now!
+    SIMFIRE_REGION   = "ANZ"       ! ("ANZ", "EUROPE", "GLOBAL")
+    HydePath         = "${BlazeDataPath}/HYDE3.1"  ! Path to Hyde3.1 population density data
+    BurnedAreaClimatologyFile = "${BlazeDataPath}/simfire_monthly_ba.nc"  ! BA climatology file (needed when blazeTStep!="annually")
+EOF
+applysed ${tmp}/sedtmp.${pid} ${ndir}/blaze.nml ${rdir}/blaze_${experiment}.nml
+cp ${rdir}/blaze_${experiment}.nml ${rdir}/blaze.nml
+
+# casafile%cnpipool                  = "restart/${mettype}_casa_rst.nc"
+# casafile%cnpepool                  = "restart/${mettype}_casa_rst.nc"
+# global Cable namelist
+cat > ${tmp}/sedtmp.${pid} << EOF
+    filename%met                       = "${mettype}"
+    filename%veg                       = "${filename_veg}"
+    filename%soil                      = "${filename_soil}"
+    filename%type                      = "${SurfaceFile}"
+    filename%out                       = "outputs/${mettype}_out_cable.nc"
+    filename%restart_in                = "restart/${mettype}_cable_rst.nc"
+    filename%restart_out               = "restart/${mettype}_cable_rst.nc"
+    casafile%cnpbiome                  = "${casafile_cnpbiome}"
+    casafile%out                       = "outputs/${mettype}_out_casa.nc"
+    casafile%cnpipool                  = "restart/${mettype}_casa"
+    casafile%cnpepool                  = "restart/${mettype}_casa"
+    cable_user%CASA_OUT_FREQ           = "monthly"
+    cable_user%POP_restart_in          = "restart/pop_${mettype}_ini.nc"
+    cable_user%POP_restart_out         = "restart/pop_${mettype}_ini.nc"
+    cable_user%LUC_restart_in          = "restart/${mettype}_luc_rst.nc"
+    cable_user%LUC_restart_out         = "restart/${mettype}_luc_rst.nc"
+    cable_user%LUC_outfile             = "outputs/${mettype}_out_luc.nc"
+    cable_user%climate_restart_in      = "restart/${mettype}_climate_rst.nc"
+    cable_user%climate_restart_out     = "restart/${mettype}_climate_rst.nc"
+    cable_user%RunIden                 = "${mettype}"
+    cable_user%MetType                 = "${mettype}"
+    output%averaging                   = "monthly"
+    output%grid                        = "land"
+    leaps                              = .false.
+    cable_user%SOIL_STRUC              = "${soil_struc}"
+    cable_user%Rubisco_parameters      = "${Rubisco_params}"
+    cable_user%CALL_POP                = .false.
+    cable_user%coordinate_photosyn     = .false.
+    cable_user%acclimate_photosyn      = .false.
+    cable_user%explicit_gm             = .false.
+    cable_user%gm_LUT_file             = "${filename_gm_lut}"
+    cable_user%CALL_BLAZE              = .false.
+    cable_user%c13o2                   = .false.
+    cable_user%c13o2_simple_disc       = .true.
+    cable_user%c13o2_delta_atm_file    = "${filename_d13c_atm}"
+    cable_user%c13o2_outfile           = "outputs/${mettype}_out_casa_c13o2.nc"
+    cable_user%c13o2_restart_in_flux   = "restart/${mettype}_c13o2_flux_rst.nc"
+    cable_user%c13o2_restart_out_flux  = "restart/${mettype}_c13o2_flux_rst.nc"
+    cable_user%c13o2_restart_in_pools  = "restart/${mettype}_c13o2_pools_rst.nc"
+    cable_user%c13o2_restart_out_pools = "restart/${mettype}_c13o2_pools_rst.nc"
+    cable_user%c13o2_restart_in_luc    = "restart/${mettype}_c13o2_luc_rst.nc"
+    cable_user%c13o2_restart_out_luc   = "restart/${mettype}_c13o2_luc_rst.nc"
+EOF
+if [[ ${call_pop} -eq 1 ]] ; then
+    sed -i -e "/cable_user%CALL_POP/s/=.*/= .true./" ${tmp}/sedtmp.${pid}
+fi
+if [[ ${coordinate_photosyn} -eq 1 ]] ; then
+    sed -i -e "/cable_user%coordinate_photosyn/s/=.*/= .true./" ${tmp}/sedtmp.${pid}
+fi
+if [[ ${acclimate_photosyn} -eq 1 ]] ; then
+    sed -i -e "/cable_user%acclimate_photosyn/s/=.*/= .true./" ${tmp}/sedtmp.${pid}
+fi
+if [[ ${explicit_gm} -eq 1 ]] ; then
+    sed -i -e "/cable_user%explicit_gm/s/=.*/= .true./" ${tmp}/sedtmp.${pid}
+fi
+if [[ ${doc13o2} -eq 1 ]] ; then
+    sed -i -e "/cable_user%c13o2/s/=.*/= .true./" ${tmp}/sedtmp.${pid}
+    if [[ ${c13o2_simple_disc} -eq 1 ]] ; then
+        sed -i -e "/cable_user%c13o2_simple_disc/s/=.*/= .true./" ${tmp}/sedtmp.${pid}
+    fi
+fi
+if [[ ${call_blaze} -eq 1 ]] ; then
+    sed -i -e "cable_user%CALL_BLAZE/s/=.*/= .true./" ${tmp}/sedtmp.${pid}
+fi
+applysed ${tmp}/sedtmp.${pid} ${ncdir}/cable.nml ${rdir}/cable_${experiment}.nml
+
+
+## BLAZE (no changes at the moment)
+if [[ ${call_blaze} -eq 1 ]] ; then
+    cp ${ndir}/blaze.nml ${rdir}/blaze.nml
+fi
+
+
+# --------------------------------------------------------------------
+# Sequence
+#
+
+
+# --------------------------------------------------------------------
+# 1. Create climate restart file
+if [[ ${doclimate} -eq 1 ]] ; then
+    echo "1. Create climate restart file"
+    rid="climate_restart"
+    # Met forcing
+    if [[ "${mettype}" == "cru" ]] ; then
+        cp ${rdir}/cru_${experiment}.nml ${rdir}/cru.nml
+    elif [[ "${mettype}" == "plume" ]] ; then
+        cp ${rdir}/plume_${experiment}.nml ${rdir}/plume.nml
+    elif [[ "${mettype}" == "bios" ]] ; then
+        cat > ${tmp}/sedtmp.${pid} << EOF
+	          Run = "spinup"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/bios_${experiment}.nml ${rdir}/bios.nml
+    fi
+    # LUC
+    cp ${rdir}/luc_${experiment}.nml ${rdir}/luc.nml
+    # Cable
+    #   do not calculate 13C because there is no 13C in the climate restart file
+    #MCTEST
+    # cable_user%YearEnd             = 1861
+    # cable_user%CASA_SPIN_ENDYEAR   = 1861
+    # cable_user%YearEnd = 1889
+    # cable_user%CASA_SPIN_ENDYEAR = 1869
+    #MCTEST
+    cat > ${tmp}/sedtmp.${pid} << EOF
+        filename%restart_in            = ""
+        cable_user%CLIMATE_fromZero    = .true.
+        cable_user%YearStart           = 1860
+        cable_user%YearEnd             = 1889
+        icycle                         = 2
+        spincasa                       = .false.
+        cable_user%CASA_fromZero       = .true.
+        cable_user%CASA_DUMP_READ      = .false.
+        cable_user%CASA_DUMP_WRITE     = .true.
+        cable_user%CASA_SPIN_STARTYEAR = 1860
+        cable_user%CASA_SPIN_ENDYEAR   = 1869
+        cable_user%limit_labile        = .true.
+        casafile%cnpipool              = ""
+        cable_user%POP_fromZero        = .true.
+        cable_user%POP_out             = "ini"
+        cable_user%POP_restart_in      = ""
+        cable_user%POPLUC              = .false.
+        cable_user%POPLUC_RunType      = "static"
+        cable_user%c13o2               = .false.
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/cable_${experiment}.nml ${rdir}/cable.nml
+    # run model
+    cd ${rdir}
+    irm logs/log_cable.txt logs/log_out_cable.txt
+    if [[ ${dompi} -eq 1 ]] ; then
+        ${mpiexecdir}mpiexec -n ${nproc} ./${iexe} > logs/log_out_cable.txt
+    else
+        ./${iexe} > logs/log_out_cable.txt
+    fi
+    # save output
+    renameid ${rid} ${mettype}.nml luc.nml cable.nml
+    imv *_${rid}.nml restart/
+    cd logs
+    renameid ${rid} log_cable.txt log_out_cable.txt
+    cd ../restart
+    copyid ${rid} ${mettype}_climate_rst.nc
+    cp ${mettype}_climate_rst.nc ${ClimateFile}
+    cd ../outputs
+    renameid ${rid} ${mettype}_out_cable.nc ${mettype}_out_casa.nc
+    cd ..
+    cd ${pdir}
+fi
+
+
+# --------------------------------------------------------------------
+# 2. First spinup phase from zero biomass
+if [[ ${dofromzero} -eq 1 ]] ; then
+    echo "2. First spinup from zero biomass"
+    rid="zero_biomass"
+    # Met forcing
+    if [[ "${mettype}" == "cru" ]] ; then
+        cp ${rdir}/cru_${experiment}.nml ${rdir}/cru.nml
+    elif [[ "${mettype}" == "plume" ]] ; then
+        cp ${rdir}/plume_${experiment}.nml ${rdir}/plume.nml
+    elif [[ "${mettype}" == "bios" ]] ; then
+        cat > ${tmp}/sedtmp.${pid} << EOF
+	          Run = "spinup"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/bios_${experiment}.nml ${rdir}/bios.nml
+    fi
+    # LUC
+    cp ${rdir}/luc_${experiment}.nml ${rdir}/luc.nml
+    # Cable
+    #MCTEST
+    # cable_user%YearEnd                = 1861
+    # cable_user%CASA_SPIN_ENDYEAR      = 1861
+    # cable_user%YearEnd = 1889
+    # cable_user%CASA_SPIN_ENDYEAR = 1869
+    # remove
+    #   cable_user%CASA_OUT_FREQ
+    #   output%averaging
+    #MCTEST
+    cat > ${tmp}/sedtmp.${pid} << EOF
+        filename%restart_in               = ""
+        cable_user%CLIMATE_fromZero       = .true.
+        cable_user%YearStart              = 1860
+        cable_user%YearEnd                = 1889
+        icycle                            = 2
+        spincasa                          = .false.
+        cable_user%CASA_OUT_FREQ          = "monthly"
+        cable_user%CASA_fromZero          = .true.
+        cable_user%CASA_DUMP_READ         = .false.
+        cable_user%CASA_DUMP_WRITE        = .true.
+        output%averaging                  = "all"
+        cable_user%CASA_SPIN_STARTYEAR    = 1860
+        cable_user%CASA_SPIN_ENDYEAR      = 1869
+        cable_user%limit_labile           = .true.
+        casafile%cnpipool                 = ""
+        cable_user%POP_fromZero           = .true.
+        cable_user%POP_out                = "ini"
+        cable_user%POP_restart_in         = ""
+        cable_user%POPLUC                 = .true.
+        cable_user%POPLUC_RunType         = "static"
+        cable_user%c13o2_restart_in_flux  = ""
+        cable_user%c13o2_restart_in_pools = ""
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/cable_${experiment}.nml ${rdir}/cable.nml
+    # run model
+    cd ${rdir}
+    irm logs/log_cable.txt logs/log_out_cable.txt
+    if [[ ${dompi} -eq 1 ]] ; then
+        ${mpiexecdir}mpiexec -n ${nproc} ./${iexe} > logs/log_out_cable.txt
+    else
+        ./${iexe} > logs/log_out_cable.txt
+    fi
+    # save output
+    renameid ${rid} ${mettype}.nml luc.nml cable.nml
+    imv *_${rid}.nml restart/
+    cd logs
+    renameid ${rid} log_cable.txt log_out_cable.txt
+    cd ../restart
+    # copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_rst.nc ${mettype}_cable_rst.nc pop_${mettype}_ini.nc
+    copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_met_rst.nc ${mettype}_casa_bal_rst.nc ${mettype}_casa_biome_rst.nc ${mettype}_casa_pool_rst.nc ${mettype}_casa_flux_rst.nc ${mettype}_casa_phen_rst.nc ${mettype}_cable_rst.nc pop_${mettype}_ini.nc
+    copyid ${rid} ${mettype}_c13o2_flux_rst.nc ${mettype}_c13o2_pools_rst.nc
+    cd ../outputs
+    renameid ${rid} ${mettype}_out_cable.nc ${mettype}_out_casa.nc ${mettype}_out_casa_c13o2.nc
+    cd ..
+    cd ${pdir}
+fi
+
+
+# --------------------------------------------------------------------
+# 3. Biomass into quasi-equilibrium with unrestricted N and P pools
+if [[ ${doequi1} -eq 1 ]] ; then
+    echo "3. Bring biomass into quasi-equilibrium with unrestricted N and P pools"
+    for ((iequi1=1; iequi1<=${nequi1}; iequi1++)) ; do
+        # 3a. 30 year run starting from restart files
+        echo "   3a. 30 year spinup from accumulated biomass; iequi1=${iequi1}/${nequi1}"
+        rid="spinup_limit_labile"
+        # rid="spinup_limit_labile${iequi}"
+        # Met forcing
+        if [[ "${mettype}" == "cru" ]] ; then
+            cp ${rdir}/cru_${experiment}.nml ${rdir}/cru.nml
+        elif [[ "${mettype}" == "plume" ]] ; then
+            cp ${rdir}/plume_${experiment}.nml ${rdir}/plume.nml
+        elif [[ "${mettype}" == "bios" ]] ; then
+            cat > ${tmp}/sedtmp.${pid} << EOF
+    	          Run = "spinup"
+EOF
+            applysed ${tmp}/sedtmp.${pid} ${rdir}/bios_${experiment}.nml ${rdir}/bios.nml
+        fi
+        # LUC
+        cp ${rdir}/luc_${experiment}.nml ${rdir}/luc.nml
+        # Cable
+        #MCTEST
+        # cable_user%YearEnd             = 1841
+        # cable_user%CASA_SPIN_ENDYEAR   = 1861
+        # cable_user%YearEnd = 1859
+        # cable_user%CASA_SPIN_ENDYEAR = 1869
+        #MCTEST
+        cat > ${tmp}/sedtmp.${pid} << EOF
+            cable_user%CLIMATE_fromZero    = .false.
+            cable_user%YearStart           = 1840
+            cable_user%YearEnd             = 1859
+            icycle                         = 2
+            spincasa                       = .false.
+            cable_user%CASA_fromZero       = .false.
+            cable_user%CASA_DUMP_READ      = .false.
+            cable_user%CASA_DUMP_WRITE     = .true.
+            cable_user%CASA_SPIN_STARTYEAR = 1860
+            cable_user%CASA_SPIN_ENDYEAR   = 1869
+            cable_user%limit_labile        = .true.
+            cable_user%POP_fromZero        = .false.
+            cable_user%POP_out             = "ini"
+            cable_user%POPLUC              = .true.
+            cable_user%POPLUC_RunType      = "static"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/cable_${experiment}.nml ${rdir}/cable.nml
+        # run model
+        cd ${rdir}
+        irm logs/log_cable.txt logs/log_out_cable.txt
+        if [[ ${dompi} -eq 1 ]] ; then
+            ${mpiexecdir}mpiexec -n ${nproc} ./${iexe} > logs/log_out_cable.txt
+        else
+            ./${iexe} > logs/log_out_cable.txt
+        fi
+        # save output
+        renameid ${rid} ${mettype}.nml luc.nml cable.nml
+        mv *_${rid}.nml restart/
+        cd logs
+        renameid ${rid} log_cable.txt log_out_cable.txt
+        cd ../restart
+        # copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_rst.nc ${mettype}_cable_rst.nc pop_${mettype}_ini.nc
+        copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_met_rst.nc ${mettype}_casa_bal_rst.nc ${mettype}_casa_biome_rst.nc ${mettype}_casa_pool_rst.nc ${mettype}_casa_flux_rst.nc ${mettype}_casa_phen_rst.nc ${mettype}_cable_rst.nc pop_${mettype}_ini.nc
+        copyid ${rid} ${mettype}_c13o2_flux_rst.nc ${mettype}_c13o2_pools_rst.nc
+        cd ../outputs
+        renameid ${rid} ${mettype}_out_cable.nc ${mettype}_out_casa.nc ${mettype}_out_casa_c13o2.nc
+        cd ..
+        cd ${pdir}
+        #
+        # 3b. analytic quasi-equilibrium of biomass pools
+        echo "   3b. Analytic solution of biomass pools"
+        rid="spinup_analytic_limit_labile"
+        # rid="spin_casa_limit_labile${iequi}"
+        # Met forcing
+        if [[ "${mettype}" == "cru" ]] ; then
+            cp ${rdir}/cru_${experiment}.nml ${rdir}/cru.nml
+        elif [[ "${mettype}" == "plume" ]] ; then
+            cp ${rdir}/plume_${experiment}.nml ${rdir}/plume.nml
+        elif [[ "${mettype}" == "bios" ]] ; then
+            cat > ${tmp}/sedtmp.${pid} << EOF
+    	          Run = "spinup"
+EOF
+            applysed ${tmp}/sedtmp.${pid} ${rdir}/bios_${experiment}.nml ${rdir}/bios.nml
+        fi
+        # LUC
+        cp ${rdir}/luc_${experiment}.nml ${rdir}/luc.nml
+        # Cable
+        #MCTEST
+        # cable_user%YearEnd             = 1841
+        # cable_user%CASA_SPIN_ENDYEAR   = 1841
+        # cable_user%YearEnd = 1859
+        # cable_user%CASA_SPIN_ENDYEAR = 1859
+        #MCTEST
+        cat > ${tmp}/sedtmp.${pid} << EOF
+            cable_user%CLIMATE_fromZero    = .false.
+            cable_user%YearStart           = 1840
+            cable_user%YearEnd             = 1859
+            icycle                         = 12
+            spincasa                       = .true.
+            cable_user%CASA_fromZero       = .false.
+            cable_user%CASA_DUMP_READ      = .true.
+            cable_user%CASA_DUMP_WRITE     = .false.
+            cable_user%CASA_SPIN_STARTYEAR = 1840
+            cable_user%CASA_SPIN_ENDYEAR   = 1859
+            cable_user%limit_labile        = .true.
+            cable_user%POP_fromZero        = .false.
+            cable_user%POP_out             = "ini"
+            cable_user%POPLUC              = .true.
+            cable_user%POPLUC_RunType      = "static"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/cable_${experiment}.nml ${rdir}/cable.nml
+        # run model
+        cd ${rdir}
+        irm logs/log_cable.txt logs/log_out_cable.txt
+        if [[ ${dompi} -eq 1 ]] ; then
+            ${mpiexecdir}mpiexec -n ${nproc} ./${iexe} > logs/log_out_cable.txt
+        else
+            ./${iexe} > logs/log_out_cable.txt
+        fi
+        # save output
+        renameid ${rid} ${mettype}.nml luc.nml cable.nml
+        mv *_${rid}.nml restart/
+        cd logs
+        renameid ${rid} log_cable.txt log_out_cable.txt
+        cd ../restart
+        # copyid ${rid} ${mettype}_casa_rst.nc pop_${mettype}_ini.nc
+        copyid ${rid} ${mettype}_casa_met_rst.nc ${mettype}_casa_bal_rst.nc ${mettype}_casa_biome_rst.nc ${mettype}_casa_pool_rst.nc ${mettype}_casa_flux_rst.nc ${mettype}_casa_phen_rst.nc pop_${mettype}_ini.nc
+        copyid ${rid} ${mettype}_c13o2_flux_rst.nc ${mettype}_c13o2_pools_rst.nc
+        if [[ ${dompi} -eq 0 ]] ; then # no output only restart if MPI
+            cd ../outputs
+            renameid ${rid} ${mettype}_out_casa.nc ${mettype}_out_casa_c13o2.nc
+            cd ..
+        fi
+        cd ${pdir}
+    done
+fi
+
+
+# --------------------------------------------------------------------
+# 4. Biomass into quasi-equilibrium with restricted N and P pools
+if [[ ${doequi2} -eq 1 ]] ; then
+    echo "4. Bring biomass into quasi-equilibrium with restricted N and P pools"
+    for ((iequi2=1; iequi2<=${nequi2}; iequi2++)) ; do
+        # 4a. 30 year run starting from restart files
+        echo "   4a. 30 year spinup from accumulated biomass; iequi2=${iequi2}/${nequi2}"
+        rid="spinup"
+        # rid="spinup${iequi}"
+        # Met forcing
+        if [[ "${mettype}" == "cru" ]] ; then
+            cp ${rdir}/cru_${experiment}.nml ${rdir}/cru.nml
+        elif [[ "${mettype}" == "plume" ]] ; then
+            cp ${rdir}/plume_${experiment}.nml ${rdir}/plume.nml
+        elif [[ "${mettype}" == "bios" ]] ; then
+            cat > ${tmp}/sedtmp.${pid} << EOF
+    	          Run = "spinup"
+EOF
+            applysed ${tmp}/sedtmp.${pid} ${rdir}/bios_${experiment}.nml ${rdir}/bios.nml
+        fi
+        # LUC
+        cp ${rdir}/luc_${experiment}.nml ${rdir}/luc.nml
+        # Cable
+        #MCTEST
+        # cable_user%YearEnd             = 1841
+        # cable_user%CASA_SPIN_ENDYEAR   = 1861
+        # cable_user%YearEnd = 1859
+        # cable_user%CASA_SPIN_ENDYEAR = 1869
+        #MCTEST
+        cat > ${tmp}/sedtmp.${pid} << EOF
+            cable_user%CLIMATE_fromZero    = .false.
+            cable_user%YearStart           = 1840
+            cable_user%YearEnd             = 1859
+            icycle                         = 2
+            spincasa                       = .false.
+            cable_user%CASA_fromZero       = .false.
+            cable_user%CASA_DUMP_READ      = .false.
+            cable_user%CASA_DUMP_WRITE     = .true.
+            cable_user%CASA_SPIN_STARTYEAR = 1860
+            cable_user%CASA_SPIN_ENDYEAR   = 1869
+            cable_user%limit_labile        = .false.
+            cable_user%POP_fromZero        = .false.
+            cable_user%POP_out             = "ini"
+            cable_user%POPLUC              = .true.
+            cable_user%POPLUC_RunType      = "static"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/cable_${experiment}.nml ${rdir}/cable.nml
+        # run model
+        cd ${rdir}
+        irm logs/log_cable.txt logs/log_out_cable.txt
+        if [[ ${dompi} -eq 1 ]] ; then
+            ${mpiexecdir}mpiexec -n ${nproc} ./${iexe} > logs/log_out_cable.txt
+        else
+            ./${iexe} > logs/log_out_cable.txt
+        fi
+        # save output
+        renameid ${rid} ${mettype}.nml luc.nml cable.nml
+        mv *_${rid}.nml restart/
+        cd logs
+        renameid ${rid} log_cable.txt log_out_cable.txt
+        cd ../restart
+        # copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_rst.nc ${mettype}_cable_rst.nc pop_${mettype}_ini.nc
+        copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_met_rst.nc ${mettype}_casa_bal_rst.nc ${mettype}_casa_biome_rst.nc ${mettype}_casa_pool_rst.nc ${mettype}_casa_flux_rst.nc ${mettype}_casa_phen_rst.nc ${mettype}_cable_rst.nc pop_${mettype}_ini.nc
+        copyid ${rid} ${mettype}_c13o2_flux_rst.nc ${mettype}_c13o2_pools_rst.nc
+        cd ../outputs
+        renameid ${rid} ${mettype}_out_cable.nc ${mettype}_out_casa.nc ${mettype}_out_casa_c13o2.nc
+        cd ..
+        cd ${pdir}
+        #
+        # 4b. analytic quasi-equilibrium of biomass pools
+        echo "   4b. Analytic solution of biomass pools"
+        rid="spinup_analytic"
+        # rid="spin_casa${iequi}"
+        # Met forcing
+        if [[ "${mettype}" == "cru" ]] ; then
+            cp ${rdir}/cru_${experiment}.nml ${rdir}/cru.nml
+        elif [[ "${mettype}" == "plume" ]] ; then
+            cp ${rdir}/plume_${experiment}.nml ${rdir}/plume.nml
+        elif [[ "${mettype}" == "bios" ]] ; then
+            cat > ${tmp}/sedtmp.${pid} << EOF
+    	          Run = "spinup"
+EOF
+            applysed ${tmp}/sedtmp.${pid} ${rdir}/bios_${experiment}.nml ${rdir}/bios.nml
+        fi
+        # LUC
+        cp ${rdir}/luc_${experiment}.nml ${rdir}/luc.nml
+        # Cable
+        #MCTEST
+        # cable_user%YearEnd             = 1841
+        # cable_user%CASA_SPIN_ENDYEAR   = 1841
+        # cable_user%YearEnd = 1859
+        # cable_user%CASA_SPIN_ENDYEAR = 1859
+        #MCTEST
+        cat > ${tmp}/sedtmp.${pid} << EOF
+            cable_user%CLIMATE_fromZero    = .false.
+            cable_user%YearStart           = 1840
+            cable_user%YearEnd             = 1859
+            icycle                         = 12
+            spincasa                       = .true.
+            cable_user%CASA_fromZero       = .false.
+            cable_user%CASA_DUMP_READ      = .true.
+            cable_user%CASA_DUMP_WRITE     = .false.
+            cable_user%CASA_SPIN_STARTYEAR = 1840
+            cable_user%CASA_SPIN_ENDYEAR   = 1859
+            cable_user%limit_labile        = .false.
+            cable_user%POP_fromZero        = .false.
+            cable_user%POP_out             = "ini"
+            cable_user%POPLUC              = .true.
+            cable_user%POPLUC_RunType      = "static"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/cable_${experiment}.nml ${rdir}/cable.nml
+        # run model
+        cd ${rdir}
+        irm logs/log_cable.txt logs/log_out_cable.txt
+        if [[ ${dompi} -eq 1 ]] ; then
+            ${mpiexecdir}mpiexec -n ${nproc} ./${iexe} > logs/log_out_cable.txt
+        else
+            ./${iexe} > logs/log_out_cable.txt
+        fi
+        # save output
+        renameid ${rid} ${mettype}.nml luc.nml cable.nml
+        mv *_${rid}.nml restart/
+        cd logs
+        renameid ${rid} log_cable.txt log_out_cable.txt
+        cd ../restart
+        # copyid ${rid} ${mettype}_casa_rst.nc pop_${mettype}_ini.nc
+        copyid ${rid} ${mettype}_casa_met_rst.nc ${mettype}_casa_bal_rst.nc ${mettype}_casa_biome_rst.nc ${mettype}_casa_pool_rst.nc ${mettype}_casa_flux_rst.nc ${mettype}_casa_phen_rst.nc pop_${mettype}_ini.nc
+        copyid ${rid} ${mettype}_c13o2_flux_rst.nc ${mettype}_c13o2_pools_rst.nc
+        if [[ ${dompi} -eq 0 ]] ; then # no output only restart if MPI
+            cd ../outputs
+            renameid ${rid} ${mettype}_out_casa.nc ${mettype}_out_casa_c13o2.nc
+            cd ..
+        fi
+        cd ${pdir}
+    done
+fi
+
+# --------------------------------------------------------------------
+# 5a. First dynamic land use (initialise land use)
+if [[ ${doiniluc} -eq 1 ]] ; then
+    echo "5a. First dynamic land use (initialise land use)"
+    # Met forcing
+    if [[ "${mettype}" == "cru" ]] ; then
+        YearStart=1580
+        YearEnd=1699
+        cp ${rdir}/cru_${experiment}.nml ${rdir}/cru.nml
+    elif [[ "${mettype}" == "plume" ]] ; then
+	      YearStart=1580
+        YearEnd=1849
+        cp ${rdir}/plume_${experiment}.nml ${rdir}/plume.nml
+    elif [[ "${mettype}" == "bios" ]] ; then
+        YearStart=1580
+        YearEnd=1699
+        cat > ${tmp}/sedtmp.${pid} << EOF
+	          Run = "spinup"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/bios_${experiment}.nml ${rdir}/bios.nml
+    fi
+    rid=${YearStart}_${YearEnd}
+    # LUC
+    cat > ${tmp}/sedtmp.${pid} << EOF
+         YearStart = ${YearStart}
+         YearEnd   = ${YearEnd}
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/luc_${experiment}.nml ${rdir}/luc.nml
+    # Cable
+    #MCTEST
+    # cable_user%CASA_SPIN_ENDYEAR    = 1859
+    # cable_user%YearEnd              = ${YearEnd}
+    # cable_user%CASA_SPIN_ENDYEAR    = 1841
+    # cable_user%YearEnd              = $(( ${YearStart} + 1 ))
+    #MCTEST
+    cat > ${tmp}/sedtmp.${pid} << EOF
+        cable_user%CLIMATE_fromZero     = .false.
+        cable_user%YearStart            = ${YearStart}
+        cable_user%YearEnd              = ${YearEnd}
+        icycle                          = 12
+        spincasa                        = .false.
+        cable_user%CASA_OUT_FREQ        = "annually"
+        cable_user%CASA_fromZero        = .false.
+        cable_user%CASA_DUMP_READ       = .true.
+        cable_user%CASA_DUMP_WRITE      = .false.
+        cable_user%CASA_SPIN_STARTYEAR  = 1840
+        cable_user%CASA_SPIN_ENDYEAR    = 1859
+        cable_user%limit_labile         = .false.
+        cable_user%POP_fromZero         = .false.
+        cable_user%POP_out              = "ini"
+        cable_user%POPLUC               = .true.
+        cable_user%POPLUC_RunType       = "init"
+        cable_user%LUC_restart_in       = ""
+        cable_user%c13o2_restart_in_luc = ""
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/cable_${experiment}.nml ${rdir}/cable.nml
+    # run model
+    cd ${rdir}
+    irm logs/log_cable.txt logs/log_out_cable.txt
+    if [[ ${dompi} -eq 1 ]] ; then
+        ${mpiexecdir}mpiexec -n ${nproc} ./${iexe} > logs/log_out_cable.txt
+    else
+        ./${iexe} > logs/log_out_cable.txt
+    fi
+    # save output
+    renameid ${rid} ${mettype}.nml luc.nml cable.nml
+    mv *_${rid}.nml restart/
+    cd logs
+    renameid ${rid} log_cable.txt log_out_cable.txt
+    cd ../restart
+    # copyid ${rid} ${mettype}_casa_rst.nc ${mettype}_luc_rst.nc pop_${mettype}_ini.nc
+    copyid ${rid} ${mettype}_casa_met_rst.nc ${mettype}_casa_bal_rst.nc ${mettype}_casa_biome_rst.nc ${mettype}_casa_pool_rst.nc ${mettype}_casa_flux_rst.nc ${mettype}_casa_phen_rst.nc ${mettype}_luc_rst.nc pop_${mettype}_ini.nc
+    copyid ${rid} ${mettype}_c13o2_pools_rst.nc ${mettype}_c13o2_luc_rst.nc
+    # cd ../outputs
+    # renameid ${rid} ${mettype}_out_luc.nc
+    # cd ..
+    cd ${pdir}
+fi
+
+
+# --------------------------------------------------------------------
+# 5b. Second full dynamic spinup
+if [[ ${doinidyn} -eq 1 ]] ; then
+    echo "5b. Full dynamic spinup"
+    # Met forcing
+    if [[ "${mettype}" == "cru" ]] ; then
+        YearStart=1700
+        YearEnd=1899
+        cat > ${tmp}/sedtmp.${pid} << EOF
+            Run = "S1_TRENDY"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/cru_${experiment}.nml ${rdir}/cru.nml
+    elif [[ "${mettype}" == "plume" ]] ; then
+	      YearStart=1850
+        YearEnd=1900
+        cat > ${tmp}/sedtmp.${pid} << EOF
+    	      Run  = "${YearStart}_${YearEnd}"
+    	      CO2  = "varying"
+    	      NDEP = "varying"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/plume_${experiment}.nml ${rdir}/plume.nml
+    elif [[ "${mettype}" == "bios" ]] ; then
+        YearStart=1700
+        YearEnd=1899
+        cat > ${tmp}/sedtmp.${pid} << EOF
+	          Run = "premet"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/bios_${experiment}.nml ${rdir}/bios.nml
+    fi
+    rid=${YearStart}_${YearEnd}
+    # LUC
+    cat > ${tmp}/sedtmp.${pid} << EOF
+         YearStart = ${YearStart}
+         YearEnd   = ${YearEnd}
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/luc_${experiment}.nml ${rdir}/luc.nml
+    # Cable
+    #MCTEST
+    # cable_user%CASA_SPIN_ENDYEAR   = 1851
+    # cable_user%YearEnd             = $(( ${YearStart} + 1 ))
+    # cable_user%CASA_SPIN_ENDYEAR = 1859
+    # cable_user%YearEnd             = ${YearEnd}
+    #MCTEST
+    cat > ${tmp}/sedtmp.${pid} << EOF
+        cable_user%CLIMATE_fromZero    = .false.
+        cable_user%YearStart           = ${YearStart}
+        cable_user%YearEnd             = ${YearEnd}
+        icycle                         = 2
+        spincasa                       = .false.
+        cable_user%CASA_fromZero       = .false.
+        cable_user%CASA_DUMP_READ      = .false.
+        cable_user%CASA_DUMP_WRITE     = .false.
+        cable_user%CASA_SPIN_STARTYEAR = 1850
+        cable_user%CASA_SPIN_ENDYEAR   = 1859
+        cable_user%limit_labile        = .false.
+        cable_user%POP_fromZero        = .false.
+        cable_user%POP_out             = "ini"
+        cable_user%POPLUC              = .true.
+        cable_user%POPLUC_RunType      = "restart"
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/cable_${experiment}.nml ${rdir}/cable.nml
+    # run model
+    cd ${rdir}
+    irm logs/log_cable.txt logs/log_out_cable.txt
+    if [[ ${dompi} -eq 1 ]] ; then
+        ${mpiexecdir}mpiexec -n ${nproc} ./${iexe} > logs/log_out_cable.txt
+    else
+        ./${iexe} > logs/log_out_cable.txt
+    fi
+    # save output
+    renameid ${rid} ${mettype}.nml luc.nml cable.nml
+    mv *_${rid}.nml restart/
+    cd logs
+    renameid ${rid} log_cable.txt log_out_cable.txt
+    cd ../restart
+    # copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_rst.nc ${mettype}_cable_rst.nc
+    copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_met_rst.nc ${mettype}_casa_bal_rst.nc ${mettype}_casa_biome_rst.nc ${mettype}_casa_pool_rst.nc ${mettype}_casa_flux_rst.nc ${mettype}_casa_phen_rst.nc ${mettype}_cable_rst.nc
+    copyid ${rid} ${mettype}_luc_rst.nc pop_${mettype}_ini.nc
+    copyid ${rid} ${mettype}_c13o2_flux_rst.nc ${mettype}_c13o2_pools_rst.nc ${mettype}_c13o2_luc_rst.nc
+    cd ../outputs
+    renameid ${rid} ${mettype}_out_cable.nc ${mettype}_out_casa.nc ${mettype}_out_luc.nc
+    renameid ${rid} ${mettype}_out_casa_c13o2.nc
+    cd ..
+    cd ${pdir}
+fi
+
+
+# --------------------------------------------------------------------
+# 6. Final centennial run
+if [[ ${dofinal} -eq 1 ]] ; then
+    echo "6. Final centennial run"
+    # Met forcing
+    if [[ "${mettype}" == "cru" ]] ; then
+        YearStart=1900
+        YearEnd=2017
+        cat > ${tmp}/sedtmp.${pid} << EOF
+            Run = "S2_TRENDY"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/cru_${experiment}.nml ${rdir}/cru.nml
+    elif [[ "${mettype}" == "plume" ]] ; then
+	      YearStart=1901
+        YearEnd=2005
+        cat > ${tmp}/sedtmp.${pid} << EOF
+    	      Run  = "${YearStart}_${YearEnd}"
+    	      CO2  = "varying"
+    	      NDEP = "varying"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/plume_${experiment}.nml ${rdir}/plume.nml
+    elif [[ "${mettype}" == "bios" ]] ; then
+        YearStart=1900
+        YearEnd=2019
+        cat > ${tmp}/sedtmp.${pid} << EOF
+	          Run = "standard"
+EOF
+        applysed ${tmp}/sedtmp.${pid} ${rdir}/bios_${experiment}.nml ${rdir}/bios.nml
+    fi
+    rid=${YearStart}_${YearEnd}
+    # LUC
+    cat > ${tmp}/sedtmp.${pid} << EOF
+         YearStart = ${YearStart}
+         YearEnd   = ${YearEnd}
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/luc_${experiment}.nml ${rdir}/luc.nml
+    # Cable
+    #MCTEST
+    # cable_user%CASA_SPIN_ENDYEAR   = 1851
+    # cable_user%YearStart           = $(( ${YearStart} + 1 ))
+    # cable_user%YearEnd             = $(( ${YearStart} + 2 ))
+    # cable_user%CASA_SPIN_ENDYEAR = 1859
+    # cable_user%YearStart           = ${YearStart}
+    # cable_user%YearEnd             = ${YearEnd}
+    #MCTEST
+    cat > ${tmp}/sedtmp.${pid} << EOF
+        cable_user%CLIMATE_fromZero    = .false.
+        cable_user%YearStart           = ${YearStart}
+        cable_user%YearEnd             = ${YearEnd}
+        icycle                         = 2
+        spincasa                       = .false.
+        cable_user%CASA_fromZero       = .false.
+        cable_user%CASA_DUMP_READ      = .false.
+        cable_user%CASA_DUMP_WRITE     = .false.
+        cable_user%CASA_SPIN_STARTYEAR = 1850
+        cable_user%CASA_SPIN_ENDYEAR   = 1859
+        cable_user%limit_labile        = .false.
+        cable_user%POP_fromZero        = .false.
+        cable_user%POP_out             = "ini"
+        cable_user%POPLUC              = .true.
+        cable_user%POPLUC_RunType      = "restart"
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/cable_${experiment}.nml ${rdir}/cable.nml
+    # run model
+    cd ${rdir}
+    irm logs/log_cable.txt logs/log_out_cable.txt
+    if [[ ${dompi} -eq 1 ]] ; then
+        ${mpiexecdir}mpiexec -n ${nproc} ./${iexe} > logs/log_out_cable.txt
+    else
+        ./${iexe} > logs/log_out_cable.txt
+    fi
+    # save output
+    renameid ${rid} ${mettype}.nml luc.nml cable.nml
+    mv *_${rid}.nml restart/
+    cd logs
+    renameid ${rid} log_cable.txt log_out_cable.txt
+    cd ../restart
+    # copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_rst.nc ${mettype}_cable_rst.nc
+    copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_met_rst.nc ${mettype}_casa_bal_rst.nc ${mettype}_casa_biome_rst.nc ${mettype}_casa_pool_rst.nc ${mettype}_casa_flux_rst.nc ${mettype}_casa_phen_rst.nc ${mettype}_cable_rst.nc
+    copyid ${rid} ${mettype}_luc_rst.nc pop_${mettype}_ini.nc
+    copyid ${rid} ${mettype}_c13o2_flux_rst.nc ${mettype}_c13o2_pools_rst.nc ${mettype}_c13o2_luc_rst.nc
+    cd ../outputs
+    renameid ${rid} ${mettype}_out_cable.nc ${mettype}_out_casa.nc ${mettype}_out_luc.nc
+    renameid ${rid} ${mettype}_out_casa_c13o2.nc
+    cd ..
+    cd ${pdir}
+fi
+
+
+# --------------------------------------------------------------------
+# 7. Future run
+if [[ (${dofuture} -eq 1) && ("${mettype}" == "plume") && ("${RCP}" != "hist") ]] ; then
+    echo "7. Future run"
+    rcpnd=$(echo ${RCP} | sed "s|\.||")
+    rcpd=$(echo ${RCP} | sed "s|\.|p|")
+    # Met forcing
+	  YearStart=2006
+    YearEnd=2099
+    cat > ${tmp}/sedtmp.${pid} << EOF
+        Run          = "${YearStart}_${YearEnd}"
+        RCP          = "${RCP}"
+        CO2          = "varying"
+        CO2file      = "${CO2Path}/co2_2006_2100_rcp${rcpnd}.dat"
+        NDEP         = "varying"
+        NDEPfile     = "${NdepPath}/RCP${rcpnd}/ndep_total_2000-2109_1.0x1.0_FD.nc"
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/plume_${experiment}.nml ${rdir}/plume.nml
+    rid=${YearStart}_${YearEnd}_${rcpd}
+    # LUC
+    cat > ${tmp}/sedtmp.${pid} << EOF
+         YearStart = ${YearStart}
+         YearEnd   = ${YearEnd}
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/luc_${experiment}.nml ${rdir}/luc.nml
+    # Cable
+    #MCTEST
+    # cable_user%CASA_SPIN_ENDYEAR = 1859
+    # cable_user%YearEnd             = ${YearEnd}
+    #MCTEST
+    cat > ${tmp}/sedtmp.${pid} << EOF
+        cable_user%CLIMATE_fromZero    = .false.
+        cable_user%YearStart           = ${YearStart}
+        cable_user%YearEnd             = $(( ${YearStart} + 1 ))
+        icycle                         = 2
+        spincasa                       = .false.
+        cable_user%CASA_fromZero       = .false.
+        cable_user%CASA_DUMP_READ      = .false.
+        cable_user%CASA_DUMP_WRITE     = .false.
+        cable_user%CASA_SPIN_STARTYEAR = 1850
+        cable_user%CASA_SPIN_ENDYEAR   = 1851
+        cable_user%limit_labile        = .false.
+        cable_user%POP_fromZero        = .false.
+        cable_user%POP_out             = "ini"
+        cable_user%POPLUC              = .true.
+        cable_user%POPLUC_RunType      = "static"
+EOF
+    applysed ${tmp}/sedtmp.${pid} ${rdir}/cable_${experiment}.nml ${rdir}/cable.nml
+    # run model
+    cd ${rdir}
+    irm logs/log_cable.txt logs/log_out_cable.txt
+    if [[ ${dompi} -eq 1 ]] ; then
+        ${mpiexecdir}mpiexec -n ${nproc} ./${iexe} > logs/log_out_cable.txt
+    else
+        ./${iexe} > logs/log_out_cable.txt
+    fi
+    # save output
+    renameid ${rid} ${mettype}.nml luc.nml cable.nml
+    mv *_${rid}.nml restart/
+    cd logs
+    renameid ${rid} log_cable.txt log_out_cable.txt
+    cd ../restart
+    # copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_rst.nc ${mettype}_cable_rst.nc pop_${mettype}_ini.nc
+    copyid ${rid} ${mettype}_climate_rst.nc ${mettype}_casa_met_rst.nc ${mettype}_casa_bal_rst.nc ${mettype}_casa_biome_rst.nc ${mettype}_casa_pool_rst.nc ${mettype}_casa_flux_rst.nc ${mettype}_casa_phen_rst.nc ${mettype}_cable_rst.nc pop_${mettype}_ini.nc
+    copyid ${rid} ${mettype}_c13o2_flux_rst.nc ${mettype}_c13o2_pools_rst.nc
+    cd ../outputs
+    renameid ${rid} ${mettype}_out_cable.nc ${mettype}_out_casa.nc
+    renameid ${rid} ${mettype}_out_casa_c13o2.nc
+    cd ..
+    cd ${pdir}
+fi
+
+
+# --------------------------------------------------------------------
+# Finish
+#
+cd ${isdir}
+
+t2=$(date +%s)
+dt=$((t2-t1))
+printf "\n"
+if [[ ${dt} -lt 60 ]] ; then
+    printf "Finished at %s   in %i seconds.\n" "$(date)" ${dt}
+else
+    dm=$(echo "(${t2}-${t1})/60." | bc -l)
+    printf "Finished at %s   in %.2f minutes.\n" "$(date)" ${dm}
+fi
+
+exit