-
Notifications
You must be signed in to change notification settings - Fork 0
/
gemmes_mod.f90
1433 lines (1203 loc) · 48.2 KB
/
gemmes_mod.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
!
!> VUA and IPSL/LSCE by the iLOVECLIM / LUDUS coding group / Within the LUDUS code environement
!
! LICENSING TERMS:
!> \copyright
!! This file is part of the iLOVECLIM coupled climate model under the LUDUS framework.
!! global_constants_mod is free software: you can redistribute it and/or modify it under the terms of the GNU General Public
!! License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later
!! version.
!!
!! global_constants_mod is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
!! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License along with Foobar.
!! If not, see <http://www.gnu.org/licenses/>.
!
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
! MODULE: gemmes_mod
!
!> @author Hugo Martin, Timothée Nicolas, Aurélien Quiquet, Didier M. Roche (dmr)
!
!
!> @brief This module gemmes_mod is used to couple iLOVECLIM with the GEMMES Integrated Assessment Model
!
!> @date Creation date: December, 4th, 2019
!> @date Last modification: $LastChangedDate$
!> @author Last modified by : ham
!
!> @version This is svn version: $LastChangedRevision$
!
!> Here add the long_description of the module ...
!! Now this module contains all the gemmes code.
!
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
! iLVC = 0 : Run Gemmes (Coping2) as an independant program.
! = 1 : Run the iLOVECLIM/GEMMES(Coping2) model full fortran.
! = 2 : Run the iLOVECLIM/GEMMES(Coping1) model fortran/R.
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
#define iLVC 0
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
! dmr stdin, stdout definition ... global variables and subroutines
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
module gemmes_mod
#if ( iLVC == 0 )
implicit none
public :: init, read_namelist, initial_conditions, &
solve, output, endgemmes
private
#else
use global_constants_mod, only: dblp=>dp
implicit none
public :: gemmes_init, gemmes_step, gemmes_accum_tglob, &
gemmes_recup_emissions, gemmes_emissions, timelov, &
init, read_namelist, initial_conditions, &
solve, output, endgemmes
private
intrinsic :: get_environment_variable,execute_command_line
real(kind=dblp),save :: timelov
real(kind=dblp),save :: glob_t2m_init
real(kind=dblp),save :: glob_t2m_accum
real(kind=dblp),save :: glob_t2m_anom
real(kind=dblp),save :: t2m_to_gemmes
real(kind=dblp),save :: gemmes_emissions
character(len=128), parameter :: gemmespath = "/home/giraud/&
&Code_Hugo/iloveclim_gemmes/gemmes/"
#endif
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
! Declaration of variables and constants.
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
! module state
real(8) :: capital
real(8) :: npop
real(8) :: debt
real(8) :: wage
real(8) :: productivity
real(8) :: price
real(8) :: eland
real(8) :: sigma
real(8) :: gsigma
real(8) :: co2at
real(8) :: co2up
real(8) :: co2lo
real(8) :: temp
real(8) :: temp0
real(8) :: pbs
real(8) :: pcar
!end module state
!module other_vars
real(8) :: n_red_fac ! reduction factor (n)
real(8) :: abat ! abatement (A)
real(8) :: dam ! damage
real(8) :: dam_k ! capital damage
real(8) :: dam_y ! output damage
real(8) :: gdp ! GDP (Y)
real(8) :: gdp0 ! GDP before abatement and damage (Y0)
real(8) :: workforce ! Workforce (L)
real(8) :: omega ! wage share wL/pY
real(8) :: lambda ! employent rate L/N
real(8) :: debtratio ! debta ratio d = D/(pY)
real(8) :: eind ! industrial emissions
real(8) :: inflation ! inflation
real(8) :: rcb ! central bank interest rate
real(8) :: pi ! profits
real(8) :: smallpi ! profit ratio pi/(pY)
real(8) :: smallpi_K ! return on assets
real(8) :: leverage ! leverage
real(8) :: cr
real(8) :: forcing ! radiative forcing
real(8) :: find ! industrial forcing
real(8) :: fexo ! exogenous forcing
real(8) :: emissions ! total emissions
real(8) :: id ! demand for investment
real(8) :: debtd ! demand for debt
real(8) :: pir ! remaining profits after payment of dividends
real(8) :: investment ! investment
!end module other_vars
!module model_pars
real(8) :: alpha ! growth rate of productivity
real(8) :: deltanpop ! leading growth rate of workforce
real(8) :: npopbar ! maximum population in the logistic evolution
real(8) :: nu ! Constant capital-to-output ratio
real(8) :: pi1 ! damage function parameter
real(8) :: pi2 ! damage function parameter
real(8) :: pi3 ! damage function parameter
real(8) :: zeta3 ! damage function parameter
real(8) :: fdamk ! (paper = 1./3) fraction of environmental damage allocated to the stock of capital
real(8) :: delta ! capital depreciation rate
real(8) :: sa ! Fraction of abatement costs that are subsidized
real(8) :: apc ! carbon price parameter
real(8) :: bpc ! carbon price parameter
real(8) :: conv10to15 ! conversion factor
real(8) :: convco2toc=1./3.666 ! conversion factor
real(8) :: deltagsigma ! dynamics of emissivity
real(8) :: eta ! relaxation parameter of inflation
real(8) :: etar ! relaxation parameter of the interest rate
real(8) :: mu ! markup of prices over the average cost
real(8) :: omitted ! offset for the production cost in the inflation
real(8) :: rstar ! Long-term interest rate target of the economy
real(8) :: phitaylor ! parameter characterizing the reactivity of the monetary policy
real(8) :: istar !interest rate targeted by the monetary policy
real(8) :: srep ! Fraction of the outstanding debt repaid yearly
real(8) :: climate_sens ! Climate sensitivity
real(8) :: gammastar ! Heat exchange coefficient between temperature layers
real(8) :: f2co2 ! Change in radiative forcing resulting from doubling of CO2
real(8) :: cat_pind ! CO2 preindustrial concentration in atmosphere
real(8) :: cup_pind ! CO2 preindustrial concentration in upper layer of ocean and biosphere
real(8) :: clo_pind ! CO2 preindustrial concentration in bottom layer of the ocean
real(8) :: fexo0 ! Initial value of the exogenous radiative forcing
real(8) :: fexo1 ! value of the exogenous radiative forcing in 2100
real(8) :: phi12 ! Transfer coefficient for carbon from the atmosphere to the upper ocean
real(8) :: phi23 ! Transfer coefficient for carbon from the upper ocean/biosphere to the lower ocean
real(8) :: deltapbs ! Exogenous growth rate of the back-stop technology price
real(8) :: theta ! parameter of the abatement cost function
real(8) :: phi0 ! Constant of the short-term Philips curve
real(8) :: phi1 ! Slope of the short-term Philips curve
real(8) :: kappa0 ! Constant of the investment function
real(8) :: kappa1 ! Slope of the investment function
real(8) :: kappamin ! Minimum of the investment function
real(8) :: kappamax ! Maximum of the investment function
real(8) :: div0 ! Constant of the dividend function
real(8) :: div1 ! Slope of the dividend function
real(8) :: divmin ! Minimum of the dividend function
real(8) :: divmax ! Maximum of the dividend function
real(8) :: heat_cap_at ! Heat capacity of the atmosphere, biosphere and upper ocean (C)
real(8) :: heat_cap_lo ! Heat capacity of the deeper ocean (C0)
real(8) :: delta_eland ! Growth rate of land-use change CO2 of emissions
real(8) :: cr0 ! Constant of the leverage function
real(8) :: crlev ! slope of the leverage function
real(8) :: g0 ! growth rate of gdp0
integer :: dam_type ! Type of damage (1 = 'No', 2 = 'Q', 3 = 'Extreme')
integer :: infla_type ! Type of inflation 1 constant, 2 with omitted
real(8), dimension(3,3) :: phimat
!end module model_pars
!module ini_cond
real(8) :: co2at_ini
real(8) :: co2up_ini
real(8) :: co2lo_ini
real(8) :: d_ini
real(8) :: eind_ini
real(8) :: eland_ini
real(8) :: gsigma_ini
real(8) :: pbs_ini
real(8) :: n_ini
real(8) :: npop_ini
real(8) :: temp_ini
real(8) :: temp0_ini
real(8) :: y_ini
real(8) :: lambda_ini
real(8) :: omega_ini
real(8) :: r_ini
real(8) :: t_ini ! time
!end module ini_cond
!module num_pars
real(8) :: dt
real(8) :: tmax
!end module num_pars
!module solution
integer :: nts
real(8), allocatable, dimension(:,:) :: sol
real(8), allocatable, dimension(:) :: time
!end module solution
#if ( iLVC == 1 )
integer :: counter_sol
#endif
contains
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
! Declaration of subroutines.
!-----|--1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2----+----3-|
! subroutine init
subroutine init
implicit none
!!!! Default variables
! Workforce
deltanpop=0.0305 ! leading growth rate of workforce
npopbar=7.056 ! maximum population in the logistic evolution
phi0=-0.292 ! Constant of the short-term Philips curve
phi1=0.469 ! Slope of the short-term Philips curve
alpha=0.02 !* growth rate of productivity
! Capital
delta=0.04 ! capital depreciation rate
nu=2.7 !* Constant capital-to-output ratio
! Dividends
div0=0.051 ! Constant of the dividend function
div1=0.4729 ! Slope of the dividend function
divmin=0. ! Minimum of the dividend function
divmax=1. ! Maximum of the dividend function; last value 1.
! Investment
kappa0=0.0397 ! Constant of the investment function
kappa1=0.719 ! Slope of the investment function
kappamax=0.3 ! Maximum of the investment function
kappamin=0. ! Minimum of the investment function
! Inflation
infla_type = 1 ! 1 constant, 2 with omitted, 3 with WACC
eta=0.192 ! relaxation parameter of inflation
mu=1.875 ! markup of prices over the average cost
omitted=0. ! offset for the production cost in the inflation
istar=0.02 ! inflation rate targeted by the monetary policy
! Interest rate
phitaylor=0.5 ! param characterizing the reactivity of the monetary policy
etar=10. ! relaxation parameter of the interest rate
rstar=0.01 ! Long-term interest rate target of the economy
! This is also the rate value if constant
srep=0.1 ! Fraction of the outstanding debt repaid yearly
! Abatement and control prices
theta=2.6 ! parameter of the abatement cost function
deltapbs=-0.005 ! Exo growth rate of the back-stop technology price
! Carbon price
conv10to15=1.160723971/1000. ! conversion factor
apc=0.125 ! carbon price parameter
bpc=0.625 ! carbon price parameter
! public - private sector
sa=0. ! Fraction of abatement costs that are subsidized
! leverage
!H pas de leverage dans le code R
cr0=0. ! Constant of the leverage function
crlev=0.! slope of the leverage function
! Climate damages
pi1=0. ! damage function parameter
pi2=0.00236 ! damage function parameter
pi3=0.0000819 ! damage function parameter
zeta3=6.754 ! damage function parameter
fdamk=1./3 ! (paper = 1./3) fraction of environmental damage allocated
! to the stock of capital
dam_type = 1 ! by default, 1: no damage
! 2: medium damage
! 3: high damage
! CO2 emissions
delta_eland=-0.022 ! Growth rate of land-use change CO2 of emissions
deltagsigma=-0.001 ! dynamics of emissivity
! Climate constants for climate module
climate_sens=3.1 ! Climate sensitivity
gammastar=0.0176 ! Heat exchange coefficient between temperature layers
f2co2=3.681 ! Change in radiative forcing resulting from doubling of CO2
cat_pind=588. ! CO2 preind conc in atmosphere
cup_pind=360. ! CO2 preind conc in upper layer of ocean and biosphere
clo_pind=1720. ! CO2 preindustrial concentration in bottom layer of the ocean
fexo0=0.5 ! Initial value of the exogenous radiative forcing
fexo1=1. ! value of the exogenous radiative forcing in 2100
phi12=0.024 ! Transfer coef for carbon from the atmo to the upper ocean
phi23=0.001 ! Transfer coef for carbon from the upper ocean/biosphere
! to the lower ocean
heat_cap_at=1/.098 ! Heat capacity of the atmo biosphere and
! upper ocean
heat_cap_lo=3.52 ! Heat capacity of the deeper ocean
! Initial conditions, values with * are equal as in the code R
y_ini=59.7387 !
npop_ini=4.8 !
lambda_ini=0.675 !
omega_ini=0.578 !
d_ini=1.53 !
eland_ini=2.6 !
eind_ini=35.85 !
n_ini=0.03 !
gsigma_ini=-0.0152 !
pbs_ini=547.22 !
co2at_ini=851. !
co2up_ini=460. !
co2lo_ini=1740. !
temp_ini=0.85 !
temp0_ini=0.0068 !
r_ini=0.10627649250000004 !H a quoi sert r_ini
! initial time
t_ini = 2015.
! Numerical parameters
dt = 0.05
tmax = 2100.
end subroutine init
! subroutine read_namelist
subroutine read_namelist
implicit none
integer :: mp
integer :: i
integer :: j
real(8) :: catup,cuplo
!namelist /model_parameters/ &
namelist /gemmes_mod/ &
alpha, &
deltanpop, &
npopbar, &
nu, &
pi1, &
pi2, &
pi3, &
zeta3, &
fdamk, &
delta, &
sa, &
apc, &
bpc, &
conv10to15, &
deltagsigma, &
eta, &
etar, &
mu, &
omitted, &
rstar, &
phitaylor, &
istar, &
srep, &
climate_sens, &
gammastar, &
f2co2, &
cat_pind, &
cup_pind, &
clo_pind, &
fexo0, &
fexo1, &
phi12, &
phi23, &
deltapbs, &
theta, &
phi0, &
phi1, &
kappa0, &
kappa1, &
kappamin, &
kappamax, &
div0, &
div1, &
divmin, &
divmax, &
heat_cap_at, &
heat_cap_lo, &
delta_eland, &
cr0, &
crlev, &
dam_type, &
infla_type, &
!namelist /initial_conditions/ &
co2at_ini, &
co2up_ini, &
co2lo_ini, &
d_ini, &
eind_ini, &
eland_ini, &
gsigma_ini, &
pbs_ini, &
n_ini, &
npop_ini, &
temp_ini, &
temp0_ini, &
y_ini, &
lambda_ini, &
omega_ini, &
r_ini, &
!namelist /numerical_parameters/ &
dt, &
tmax
mp=1
open(mp, &
#if ( iLVC == 0 )
& file='gemmes.dat', &
#elif ( iLVC == 1 )
& file='../../../gemmes_cpl/sources/gemmes.dat', &
#endif
& delim= 'apostrophe', &
& form = 'formatted', &
& action = 'read', &
& status = 'old')
read(mp,gemmes_mod)
close(mp)
catup = cat_pind/cup_pind
cuplo = cup_pind/clo_pind
phimat(1,1) = -phi12
phimat(1,2) = phi12*catup
phimat(2,1) = phi12
phimat(2,2) = -phi12*catup - phi23
phimat(2,3) = phi23*cuplo
phimat(3,2) = phi23
phimat(3,3) = -phi23*cuplo
#if (iLVC == 1)
t_ini = 2015
tmax = 3015
nts = 1000
dt = 1.
#else
nts = floor((tmax-t_ini+1)/dt)+1
#endif
allocate(sol(nts,33),time(nts))
do i=1,nts
time(i) = t_ini + dt*real(i-1,8)
do j=1,33
sol(i,j) = 0.
enddo
enddo
end subroutine read_namelist
! subroutine initial_conditions
subroutine initial_conditions
implicit none
real(8) :: tc
real(8) :: transfers
real(8) :: kappapi
real(8) :: deltapik
real(8) :: philam
real(8) :: deltad
real(8) :: costprod
real(8) :: WACC
! carbon price
pcar = pbs_ini*n_ini**(theta - 1.)
! emissivity
sigma = eind_ini/((1.-n_ini)*y_ini)
! reduction emission factor and abatement
abat = 0.001*sigma*pbs_ini*n_ini**theta/theta
select case(dam_type)
case(1)
dam = 0.
case(2)
dam = 1. - 1./(1 + pi1*temp_ini + pi2*temp_ini**2)
case(3)
dam = 1.-1./(1 + pi1*temp_ini + pi2*temp_ini**2 &
& + pi3*temp_ini**zeta3)
end select
! Damage
dam_k = fdamk*dam
dam_y = 1. - (1.-dam)/(1.-dam_k)
deltad = delta + dam_k
! Total cost of climate change
tc = (1.-dam_y)*(1.-abat)
! GDP
gdp0 = y_ini/tc
capital = gdp0*nu
! workforce
workforce = lambda_ini*npop_ini
! Wage share
productivity = gdp0/workforce
! Initial price index
price = 1.
! wages
wage = omega_ini*price*y_ini/workforce
! Net transfers between public and private sectors
transfers = sa*abat*gdp0 - pcar*conv10to15*eind_ini
! Nominal debt
debt = d_ini*price*y_ini
! inflation for infla_type=1, 2 and rate
if (infla_type<3) then
if (infla_type.eq.1) then
inflation = istar
rcb = rstar
else
inflation = eta*(mu*(omega_ini + omitted) - 1.)
! central bank interest rate
call Taylor(inflation,rcb)
endif
else
rcb = rstar
endif
! profits
pi = price*y_ini - wage*workforce - rcb*debt &
& - deltad*price*capital + price*transfers
smallpi_k = pi/(price*capital)
smallpi = pi/(price*y_ini)
! Dividends
call Dividends(smallpi_k,deltapik)
! inflation for infla_type = 3
if (infla_type.eq.3) then
WACC = (rcb*debt + deltapik*price*capital)/price/capital
costprod = 1./tc*(wage/price/productivity &
& + pcar*conv10to15*sigma*(1-n_ini) &
& + nu*(WACC + deltad))
inflation = eta*(mu*costprod - 1.)
endif
! leverage
leverage = debtratio*tc/nu
call Tau(leverage,cr)
! emissions
find = f2co2*log(co2at_ini/cat_pind)/log(2.)
fexo = fexo0
forcing = find + fexo
emissions = eind_ini + eland_ini
! leverage
leverage = d_ini*tc/nu
call Tau(leverage,cr)
! growth rate of gpd0
call kappa(smallpi,kappapi)
g0 = ((1.-cr)*(1.25*kappapi*tc/nu - deltad) &
& +cr*(smallpi_k - deltapik - tc*srep*d_ini/nu))
! initial carbon concentrations
co2at = co2at_ini
co2up = co2up_ini
co2lo = co2lo_ini
! other variables that were not set yet
npop = npop_ini
eland = eland_ini
gsigma = gsigma_ini
temp = temp_ini
temp0 = temp0_ini
pbs = pbs_ini
end subroutine initial_conditions
! subroutine solve
subroutine solve
implicit none
integer :: i
real(8) :: t
real(8), dimension(16) :: u
! initialisation of solution
time(1) = t_ini
sol(1,1) = capital
sol(1,2) = npop_ini
sol(1,3) = debt
sol(1,4) = wage
sol(1,5) = productivity
sol(1,6) = price
sol(1,7) = eland_ini
sol(1,8) = sigma
sol(1,9) = gsigma_ini
sol(1,10) = co2at_ini
sol(1,11) = co2up_ini
sol(1,12) = co2lo_ini
sol(1,13) = temp_ini
sol(1,14) = temp0_ini
sol(1,15) = pbs_ini
sol(1,16) = pcar
sol(1,17) = omega_ini
sol(1,18) = lambda_ini
sol(1,19) = d_ini
sol(1,20) = gdp0
sol(1,21) = y_ini
sol(1,22) = eind_ini
sol(1,23) = inflation
sol(1,24) = abat
sol(1,25) = n_ini
sol(1,26) = smallpi
sol(1,27) = smallpi_k
sol(1,28) = dam
sol(1,29) = dam_k
sol(1,30) = dam_y
sol(1,31) = fexo
sol(1,32) = find
sol(1,33) = rcb
! state init
u(1) = capital
u(2) = npop
u(3) = debt
u(4) = wage
u(5) = productivity
u(6) = price
u(7) = eland
u(8) = sigma
u(9) = gsigma
u(10) = co2at
u(11) = co2up
u(12) = co2lo
u(13) = temp
u(14) = temp0
u(15) = pbs
u(16) = pcar
do i=2,nts
t = time(i-1)
call rk4(t,u)
time(i) = t_ini+dt*real(i-1,8)
sol(i,1:16) = u
sol(i,17) = omega
sol(i,18) = lambda
sol(i,19) = debtratio
sol(i,20) = gdp0
sol(i,21) = gdp
sol(i,22) = eind
sol(i,23) = inflation
sol(i,24) = abat
sol(i,25) = n_red_fac
sol(i,26) = smallpi
sol(i,27) = smallpi_k
sol(i,28) = dam
sol(i,29) = dam_k
sol(i,30) = dam_y
sol(i,31) = fexo
sol(i,32) = find
sol(i,33) = rcb
end do
end subroutine solve
! subroutine rk4
subroutine rk4(t,u)
implicit none
real(8), intent(in) :: t
real(8) :: t_aux
real(8), dimension(16), intent(inout) :: u
real(8), dimension(16) :: u_aux,k1,k2,k3,k4
# if ( iLVC == 1 )
real(8) :: temp_aux
temp_aux = u(13)
#endif
call system_function(t,u,k1)
u_aux = u + k1*0.5*dt
t_aux = t + 0.5*dt
# if ( iLVC == 1 )
u(13) = temp_aux
#endif
call system_function(t_aux,u_aux,k2)
u_aux = u + k2*0.5*dt
# if ( iLVC == 1 )
u(13) = temp_aux
#endif
call system_function(t_aux,u_aux,k3)
u_aux = u + k3*dt
t_aux = t + dt
# if ( iLVC == 1 )
u(13) = temp_aux
#endif
call system_function(t_aux,u_aux,k4)
u = u + (k1+2.*(k2+k3)+k4)*dt/6.
# if ( iLVC == 1 )
u(13) = temp_aux
#endif
!H carbon price if like in R code
if (u(16).gt.u(15)) then
u(16) = u(15)
end if
end subroutine rk4
! subroutine system_function
subroutine system_function(t,u,k)
implicit none
real(8) :: t
real(8), dimension(16), intent(in) :: u
real(8), dimension(16), intent(out) :: k
! Intermediate variables
real(8) :: t2015, t2100
real(8) :: tc
real(8) :: ttilde
real(8) :: transfers
real(8) :: kappapi
real(8) :: deltapik
real(8) :: philam
real(8) :: deltad
real(8) :: beta
real(8) :: rho
real(8) :: costprod
real(8) :: WACC
! dstate/dt
real(8) :: kdot
real(8) :: ndot
real(8) :: ddot
real(8) :: wdot
real(8) :: adot
real(8) :: pdot
real(8) :: elanddot
real(8) :: sigmadot
real(8) :: gsigmadot
real(8) :: co2atdot
real(8) :: co2updot
real(8) :: co2lodot
real(8) :: tdot
real(8) :: t0dot
real(8) :: pbsdot
real(8) :: pcdot
rho = f2co2/climate_sens
! expand state
capital = u(1)
npop = u(2)
debt = u(3)
wage = u(4)
productivity = u(5)
price = u(6)
eland = u(7)
sigma = u(8)
gsigma = u(9)
co2at = u(10)
co2up = u(11)
co2lo = u(12)
temp = u(13)
temp0 = u(14)
pbs = u(15)
pcar = u(16)
!H if like in R code
if (pcar.gt.pbs) then
pcar = pbs
end if
! time intermediate variables
t2015 = 2015.
t2100 = 2100.
! Say's law
call Sayslaw(npop,capital,productivity,gdp0,workforce,lambda)
! reduction emission factor and abatement
n_red_fac = min((pcar/((1.-sa)*pbs))**(1./(theta-1.)),1.)
abat = 0.001*sigma*pbs*n_red_fac**theta/theta
select case(dam_type)
case(1)
dam = 0.
case(2)
dam = 1. - 1./(1 + pi1*temp + pi2*temp**2)
case(3)
dam = 1. - 1./(1 + pi1*temp + pi2*temp**2 &
& + pi3*temp**zeta3)
end select
! Damage
dam_k = fdamk*dam
dam_y = 1. - (1.-dam)/(1.-dam_k)
deltad = delta + dam_k
! Total cost of climate change
tc = (1.-dam_y)*(1.-abat)
! GDP
gdp = gdp0*tc
! Wage share
omega = wage*workforce/(price*gdp)
! debt ratio
debtratio = debt/(price*gdp)
! Industrial emission
eind = gdp0*sigma*(1.-n_red_fac)
! inflation and rate for infla_type=1, 2
if (infla_type<3) then
if (infla_type.eq.1) then
inflation = istar
rcb = rstar
else
inflation = eta*(mu*omega - 1.)
! central bank interest rate
call Taylor(inflation,rcb)
endif
else
rcb = rstar
endif
! Net transfers between public and private sectors
transfers = sa*abat*gdp0 - pcar*conv10to15*eind
! profits
pi = price*gdp - wage*workforce - rcb*debt &
& - deltad*price*capital + price*transfers
smallpi_k = pi/(price*capital)
smallpi = pi/(price*gdp)
! Dividends
call Dividends(smallpi_k,deltapik)
! inflation for infla_type=3
if (infla_type.eq.3) then
WACC = (rcb*debt + deltapik*price*capital)/price/capital
costprod = 1./tc*(wage/price/productivity &
& + pcar*conv10to15*sigma*(1-n_red_fac) &
& + nu*(WACC + deltad))
inflation = eta*(mu*costprod - 1.)
endif
! leverage
leverage = debtratio*tc/nu
call Tau(leverage,cr)
! population growth
beta = deltanpop*(1.-npop/npopbar)
! emissions
find = f2co2*log(co2at/cat_pind)/log(2.)
ttilde = (t-t2015)/(t2100-t2015)
fexo = min(fexo1*ttilde+fexo0*(1-ttilde), fexo1)
forcing = find + fexo
emissions = eind + eland
! investment
call kappa(smallpi,kappapi)
id = kappapi*gdp
pir = pi - deltapik*price*capital
investment = cr*(pir/price + deltad*capital &
& - srep*debt/price) + (1. - cr)*id
! debt demand
debtd = price*id - pir + srep*debt - deltad*price*capital
! computation of final results
kdot = investment - deltad*capital
ndot = beta*npop
ddot = (1.-cr)*debtd - srep*debt
call Phi(lambda,philam)
wdot = philam*wage
adot = productivity*alpha
pdot = inflation*price
elanddot = delta_eland*eland
sigmadot = gsigma*sigma
gsigmadot = deltagsigma*gsigma
call co2dot(emissions,co2at,co2up, &
& co2lo,co2atdot,co2updot,co2lodot)
tdot = (forcing - rho*temp &
& - gammastar*(temp-temp0))/heat_cap_at
t0dot = gammastar*(temp-temp0)/heat_cap_lo
pbsdot = pbs*deltapbs
pcdot = pcar*(apc + bpc/t)
g0 = ((1.-cr)*(1.25*kappapi*tc/nu - deltad) &
& +cr*(smallpi_k - deltapik - tc*srep*debtratio/nu))
k(1) = kdot
k(2) = ndot
k(3) = ddot
k(4) = wdot
k(5) = adot
k(6) = pdot
k(7) = elanddot
k(8) = sigmadot
k(9) = gsigmadot
k(10) = co2atdot
k(11) = co2updot
k(12) = co2lodot
k(13) = tdot
k(14) = t0dot
k(15) = pbsdot
k(16) = pcdot
end subroutine system_function
! subroutine Taylor
subroutine Taylor(inflation,rcb)
implicit none
real(8) :: inflation
real(8) :: rcb
rcb = min(max( &
& rstar + inflation + phitaylor*(inflation - istar) &
& , 0.), 1.)
end subroutine Taylor
! subroutine Sayslaw
subroutine Sayslaw(npop,capital,productivity,gdp0, &
& workforce,lambda)
implicit none
real(8) :: npop
real(8) :: capital
real(8) :: productivity
real(8) :: gdp0
real(8) :: workforce
real(8) :: lambda
! GDP
gdp0 = capital/nu
! workforce
workforce = gdp0/productivity
! check Say's law and set employment parameter
if (workforce>npop) then
lambda = 1
workforce = npop
gdp0 = productivity*workforce
capital = nu*gdp0
else
lambda = workforce/npop
endif
end subroutine Sayslaw