Skip to content
Snippets Groups Projects
utils_par.f 8.88 KiB
!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!  
!*                                                               *!
!*     SUBROUTINES                                               *!
!*     - assign_CO2par                                           *!
!*     FUNCTIONS                                                 *!
!*     - CO2_annual                                              *!
!*     - CO2_hist                                                *!
!*                                                               *!
!*                  Copyright (C) 1996-2018                      *!
!*     Potsdam Institute for Climate Impact Reserach (PIK)       *!
!*          Authors and contributors see AUTHOR file             *!
!*  This file is part of 4C and is licensed under BSD-2-Clause   *!
!*                   See LICENSE file or under:                  *!
!*     http://www.https://opensource.org/licenses/BSD-2-Clause   *!
!*                           Contact:                            *!
!*       https://gitlab.pik-potsdam.de/foresee/4C                *!
!*                                                               *!
!*****************************************************************!
    
SUBROUTINE assign_CO2par
! Tables of parameters for calculation of CO2 scenarios

USE data_climate
USE data_par
USE data_simul

IMPLICIT NONE
DOUBLE PRECISION co2_annual

if (flag_co2 .ge. 250) then
    co2 = flag_co2/1000000.

else if (flag_co2 .eq. 0) then
    co2 = co2_st    

else
       ! historical CO2 increase, function fitted by Kohlmaier et al.
       p1_co2h = 0.000295
       p2_co2h = 0.027
       p3_co2h = 1860
       p4_co2h = 0.00616

    select case (flag_co2)
    case (101, 201)
       p1_co2 = 0.000295
       p2_co2 = 0.027
       p3_co2 = 1860
       p4_co2 = 0.00616

    case (102, 202)
    !IF (flag_co2==102.OR.flag_co2==202) THEN
       ! scenario function used for LTEEF II and SILVISTRAT
       p1_co2 = 0.000350
       p2_co2 = 0.0063
       p3_co2 = 1990

    case (103, 203)
       ! Mauna Loa 
       continue 

    case(110, 210)
    !ELSEIF (flag_co2==110.OR.flag_co2==210) THEN
       ! IPCC IS92a after Bern CC model, reference
       p0_co2 = -513178.349650788
       p1_co2 = 774.68578088642
       p2_co2 = -0.39065850816
       p3_co2 = 0.00006585082

    case(111,211)
    !ELSEIF (flag_co2==111.OR.flag_co2==211) THEN
       ! IPCC A1FI after Bern CC model, reference
       p0_co2 = 1818104.0398310008
       p1_co2 = -2584.5240137942828
       p2_co2 = 1.2208975180122
       p3_co2 = -0.0001915345027

    case(112, 212)
    !ELSEIF (flag_co2==112.OR.flag_co2==212) THEN
       ! IPCC A2 after Bern CC model, reference
       p0_co2 = -1045454.7878788
       p1_co2 = 1587.54094794094
       p2_co2 = -0.804265734265715
       p3_co2 = 0.00013597513597513

    case(113, 213)
    !ELSEIF (flag_co2==113.OR.flag_co2==213) THEN
       ! IPCC B1 after Bern CC model, reference
        p0_co2 = 1596094.36363588
        p1_co2 = -2362.17634032563
        p2_co2 = 1.16444055944021
        p3_co2 = -0.000191142191142135

    case(114, 214)
    !ELSEIF (flag_co2==114.OR.flag_co2==214) THEN
       ! IPCC B2 after Bern CC model, reference
       p0_co2 = 152455.527149544
       p1_co2 = - 213.773160908033
       p2_co2 = 0.0988590820945227
       p3_co2 = -0.000014997257644339

    case(115, 215)
    !ELSEIF (flag_co2==115.OR.flag_co2==215) THEN
       ! IPCC A1B after Bern CC model, reference
       p0_co2 = 1955425.02331
       p1_co2 = - 2858.095994844593
       p2_co2 = 1.39094405594
       p3_co2 = -0.00022533023

    case(116, 216)
    !ELSEIF (flag_co2==116.OR. flag_co2 == 216) THEN
       ! IPCC A1p after Bern CC model, reference
       p0_co2 = 1872081.750583
       p1_co2 = -2742.46196581203
       p2_co2 = 1.33764568765
       p3_co2 = -0.00021717172

    case(117, 217)
       ! RCP8.5, FINAL RELEASE, 26 Nov. 2009 
       ! DOCUMENTATION: M. Meinshausen, S. Smith et al., 2011: "The RCP GHG concentrations and their extension from 1765 to 2500", 
       ! Climatic Change, 109(1-2), S. 213-241.
       
       p0_co2 = 179973.892732277
       p1_co2 = -180.746725115325
       p2_co2 = 0.0454730127294021
       p3_co2 = 0.

    case(118, 218)
       ! RCP2.6, FINAL RELEASE, 26 Nov. 2009 
       ! DOCUMENTATION: M. Meinshausen, S. Smith et al., 2011: "The RCP GHG concentrations and their extension from 1765 to 2500", 
       ! Climatic Change, 109(1-2), S. 213-241.
       
       p0_co2 = 166355.928340573
       p1_co2 = -285.793245173113
       p2_co2 = 0.16016037734385
       p3_co2 = -0.00002937992775

    case(119, 219)
       ! RCP4.5, FINAL RELEASE, 26 Nov. 2009 
       ! DOCUMENTATION: M. Meinshausen, S. Smith et al., 2011: "The RCP GHG concentrations and their extension from 1765 to 2500", 
       ! Climatic Change, 109(1-2), S. 213-241.
       
       p0_co2 = 173604.151969275
       p1_co2 = -286.531541491431
       p2_co2 = 0.15501922547777
       p3_co2 = -0.00002753265703

    case(120, 220)
       ! RCP6.0, FINAL RELEASE, 26 Nov. 2009 
       ! DOCUMENTATION: M. Meinshausen, S. Smith et al., 2011: "The RCP GHG concentrations and their extension from 1765 to 2500", 
       ! Climatic Change, 109(1-2), S. 213-241.
       
       p0_co2 = 70777.
       p1_co2 = -71.604
       p2_co2 = 0.0182
       p3_co2 = 0.0


    end select
    
    co2 = co2_annual(time_cur)
end if

END 

!*****************************************************************************

DOUBLE PRECISION FUNCTION CO2_annual(int_time)
! calculates annual atmospheric CO2 mixing ratio for scenarios

USE data_climate
USE data_par
USE data_simul

IMPLICIT NONE

Integer int_time
REAL x_time
DOUBLE PRECISION CO2_hist

! variable x_time foreseen for choice of year of step change
! help variable
x_time = real(int_time)

! first set of functions for continuous scenarios, flag_co2 values < 200
IF(flag_co2<200) THEN
  select case (flag_co2)
  case(101)
!  IF (flag_co2==101) THEN
! historical increase (Kohlmaier)
    CO2_annual=(p4_co2*(exp(p2_co2*(x_time-p3_co2))-1.)+1.)*p1_co2    

  case(102)
!  ELSE IF (flag_co2==102) THEN
! LTEEF, SILVISTRAT Scenarios
     CO2_annual=p1_co2*exp(p2_co2*((x_time-1)-p3_co2))     

  case(103)
!  ELSE IF (flag_co2==103) THEN
! Mauna Loa
        CO2_annual = CO2_hist(x_time)
 
  case(111, 112, 113, 114, 115, 116, 117, 118,119,120) 
!  ELSE IF (flag_co2==111.OR.flag_co2==114.OR.flag_co2==112.OR.flag_co2==113.OR.flag_co2==115) THEN
! Sceanrio A1F1, B1, B2, A2 
     IF(x_time > year_CO2) THEN
        CO2_annual = p0_co2 + p1_co2*x_time + p2_co2*x_time*x_time + p3_co2*x_time*x_time*x_time
        CO2_annual = CO2_annual/1000000.
     ELSE
        CO2_annual = CO2_hist(x_time)
     ENDIF
     
 case(131)
      
       CO2_annual = RCP_2p6(x_time)/1000000.
 case(132)
       if((time .gt. 0) .and. (x_time.le.2150)) then
         CO2_annual = RCP_6p0(x_time)/1000000.
       else
         CO2_annual = RCP_6p0(2150)/1000000.
       end if
 case (133)
        if((time .gt. 0) .and. (x_time .le. 2005)) then
           CO2_annual = RCP_2p6(x_time)/1000000.
        else
           CO2_annual= 378.81/1000000.
        end if
        
  end select   ! flag_co2
ELSE
! second set of functions for step change scenarios, flag_co2 values > 200
! step change in the middle of the simulation period
! in the first (second) half the CO2 partial pressure of the start (end) year is used
  IF (flag_co2==201) THEN
     IF(time < year/2) THEN
        CO2_annual=(p4_co2*(exp(p2_co2*(time_b-p3_co2))-1.)+1.)*p1_co2
     ELSE
        CO2_annual=(p4_co2*(exp(p2_co2*(year+time_b-p3_co2))-1.)+1.)*p1_co2
     ENDIF
  ELSE IF (flag_co2==202) THEN
     IF(time < year/2) THEN
        CO2_annual=p1_co2*exp(p2_co2*((time_b-1)-p3_co2))
     ELSE
        CO2_annual=p1_co2*exp(p2_co2*((year+time_b-1)-p3_co2))
     ENDIF
! Hyytil 1995 - 2009
  ELSE IF (flag_co2 == 199) THEN
        if(x_time .gt.1994 .and. x_time .lt.2010) then
		   CO2_annual = (1.8607 * (x_time-1) -3353)/ 1000000.
		else
		   CO2_annual = (1.8607 * 2009 -3353)/ 1000000.
        end if
  ENDIF
  
ENDIF

END function   ! CO2_annual

!*****************************************************************************

DOUBLE PRECISION FUNCTION CO2_hist(x_time1)
! calculates annual atmospheric CO2 mixing ratio for historical times

USE data_climate
USE data_out
USE data_par
USE data_simul

IMPLICIT NONE

REAL x_time1
integer i_time1

! Mauna Loa
   IF(x_time1 > year_CO2 .or. x_time1 < 1959) then
   write (unit_err,*)'Mauna Loa can only be used for 1959 through ', year_CO2
   write (unit_err,*)'calendar year: ',x_time1,' is outside this range'
   write (unit_err,*)'Application of historical increase'
        CO2_hist=(p4_co2h*(exp(p2_co2h*(x_time1-p3_co2h))-1.)+1.)*p1_co2h
   ELSE
        i_time1 = int(x_time1)
        CO2_hist=Mauna_Loa_CO2(i_time1)
   ENDIF

END function   ! CO2_hist