Climate Models Implementation

Parent Previous Next

Find below the codification in Fortram which implement the calculation models of concentration, radiative forcing, temperature and sea level.




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!! Evaluate CO2 model 1 !!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!>

subroutine evaluateModelCO2MA()

implicit none


! Variable declarations

real*8 :: a  

logical :: b

real*8 :: C0

real*8 :: CC

real*8 :: d  

real*8 :: fert_b  

real*8 :: fnpp  

real*8 :: GCO2  

real*8 :: Ginf  

real*8 :: k1at  

real*8 :: k2at  

real*8 :: kcb1  

integer :: npp  

real*8 :: pCO2ml  

real*8 :: r  

real*8 :: sigDIC  

real*8 :: y  

integer ::i


!

! calc with australian model

! ocean

now.td.cs= sum (now.st.csv(1:9)) !(f.1,5...)

sigDIC= now.td.cs/321.21 !(f.5)

select case (trim(par.co2.pco2choice)) !

case ( 'enting') !

y= sigDIC/par.co2.dicpre !(f.7)

pCO2ml= par.co2.pCO2pre*(1+y*revelle(y)) !(f.6)

case ( {'joos', 'co2syspade', 'co2syspower'} ) !! 'co2syspade', 'co2syspower' not implemented

pCO2ml= par.co2.pCO2pre + (1000*z0(now.td.temp)*sigDIC)/(1-1000*z1(now.td.temp)*sigDIC) !(f.9)

case default

disp('Opps&, Unexpected par.use_co2 option at evaluateModelCo2')

end select

now.td.fluxco2as= par.co2.kg*(now.st.co2-pCO2ml)/par.co2.rCO2 !(f.4)

now.td.fluxco2bs= par.co2.kbs*now.st.cb1


! terrestrial carbon cycle

select case (trim( par.co2.fertchoice)) !

case ( 'log') !

npp= par.co2.npppre*(1+par.co2.beta*log(now.st.co2/par.co2.co2pre)) !(f.13)

case ( {'hyp'}) !

Ginf= 2.4

CC= 80

fnpp= 0.81

C0= par.co2.co2pre ! verificar se esta correto

d= (Ginf-1)*C0 - Ginf*CC

GCO2= Ginf*(now.st.co2 - CC)/(now.st.co2 + d) ! (f.15)

npp= par.co2.npppre*(1+fnpp*(GCO2 - 1)) !(f.14)

case ( {'hypm'}) !

CC= 31

a= 680 - CC

b= 340 - CC

r= (1+par.co2.beta*log(680/par.co2.co2pre))/(1+par.co2.beta*log(340/par.co2.co2pre)) !(f.18)

fert_b = (a - r*b)/((r-1)*a*b)

npp= par.co2.npppre*(fert_b + 1/(par.co2.npppre - CC))/(fert_b + 1/(now.st.co2 - CC)) !(f.16)

case default

disp('Opps&, Unexpected par.co2.fertchoice option at evaluateModelCo2')

end select


! atmospheric co2

! npp and fluxco2as variation with temperature

now.td.npp= npp*par.co2.q10npp**((now.td.temp - par.temp.temppre)/10) !(f.19)

k1at= par.co2.k1a*par.co2.q10npp**((now.td.temp - par.temp.temppre)/10) !(f.19);

k2at= par.co2.k2a*par.co2.q10npp**((now.td.temp - par.temp.temppre)/10) !(f.19);

now.td.fluxco2ab= now.td.npp - k1at*now.st.cb1 - k2at*now.st.cb2 !(f.21)


! stateDer

! ocean

now.stDer.csv(1:9)= par.co2.AML*(now.td.fluxco2as + now.td.fluxco2bs) - &

par.co2.ALFAML.*now.st.csv(1:9) !(f.1)

now.stDer.csv(10)= dot_product(par.co2.ALFAML,now.st.csv(1:9)) !(f.3)

! terrestrial carbon cycle

kcb1= par.co2.k1a + par.co2.k12 + par.co2.kbs + par.co2.km

now.stDer.cb1= par.co2.gamma*now.td.npp - kcb1*now.st.cb1

now.stDer.cb2= (1-par.co2.gamma)*now.td.npp + par.co2.k12*now.st.cb1 &

- par.co2.k2a*now.st.cb2

! atmospheric

now.stDer.q= now.td.ECO2_foss + now.td.ECO2_luc

select case (trim( par.co2.lucchoice)) !

case ( 'none') !

now.td.ECO2_luc= 0

case ( {'added'}) !

!

case ( {'cb2'}) !

now.stDer.cb2= now.stDer.cb2 - now.td.ECO2_luc

case default ! gross NOT IMPLEMENTED

disp('Opps&, Unexpected par.co2.lucchoice option at evaluateModelCo2')

end select

now.stDer.co2= now.td.ECO2_foss + now.td.ECO2_luc &

+ now.td.ch42co2emiss + now.td.ECO2_volc &

- now.td.fluxco2as - now.td.fluxco2ab

now.stDer.co2= par.co2.rCO2*now.stDer.co2 ! conversion to ppm



end subroutine evaluateModelCO2MA




subroutine co2bern_model()

implicit none

!        local variables

! [i ] i,n        : index for looping

! [i ] icheck_ocean: variable to check whether only one ocean model is set

! [d ] check_max: used to check array length of internal vectors


         integer n,i,icheck_ocean

       double precision check_max


!        -----------------------------------------------------

!        set time step used to solve convolution integral (dt)

!        SHOULD BE REDUCED IF NUMERICAL INSTABILITIES OCCUR

!       -----------------------------------------------------

         co2var.dt      = 0.5d0

!

!        -------

!        CHECKS:

!        -------

!


!        1) check max array length for internal vectors:


       check_max = co2var.n_data*co2var.dt_data/co2var.dt

       

   if (check_max .ge. max_array_int) then

               stop

       endif


!        2) check if dt_data is a multiple of dt


       if ( (co2var.dt_data/co2var.dt)-floor(co2var.dt_data/co2var.dt) .gt. 1.d-10) then

               stop

       endif                


!        3) check if only one ocean model is set


       icheck_ocean =0

   

   if(control.boxdiff==1) then

           icheck_ocean = icheck_ocean + 1

   end if

   

   if(control.hildaeff==1) then

           icheck_ocean = icheck_ocean + 1

   end if


   if(control.dim2==1) then

           icheck_ocean = icheck_ocean + 1

   end if


   if(control.dim3==1) then

           icheck_ocean = icheck_ocean + 1

   end if

   

       if (icheck_ocean .ge. 2) then

               stop

       endif                



!--------------------------------------

!        INITIALIZATION

!--------------------------------------

       call co2bern_init

!---------------------------------------

!        INTEGRATION: solve convolution integral:

!        and OUTPUT to PCO2A_VEC

!---------------------------------------

             do n=1,co2var.nmax

!   --- set time, position in data vector and flag for output:


!                year at begining of integration interval

             co2var.yr = co2var.yrbegin+(n-1)*co2var.dt


               if (co2var.yr .ge. par.tvet(co2var.ipos_data+1)) then

                       co2var.ipos_data=co2var.ipos_data+1

                       co2var.flagout = .true.

       endif


       !net ecosystem production (nep) due to CO2 fertilization ----


       co2var.nep = 0.0d0


       !pulse response function of the Bern 4-box biopshere.                    

       if (control.boxdiff==1) then

                     co2var.npp_vec(n) = co2par.npp0*co2par.beta*log(co2var.pco2a_now/co2par.pco20bio)

                     co2var.hetero_resp = 0.0d0

                     do i=1,n

                   co2var.hetero_resp = co2var.hetero_resp+co2var.npp_vec(i)* co2var.respbio_vec(n+1-i)*co2var.dt

           enddo

                     co2var.nep = co2var.npp_vec(n)-co2var.hetero_resp

       end if


!     --- net air-to-sea carbon-flux (time= yr = begin of integration interval) ---


!                net air-sea flux in [mol/yr]

       co2var.fnet_air_sea_vec(n) =  co2par.a_kg*(co2var.pco2a_now - co2var.pco2s)/co2par.ppm_per_mol


!     --- determine emission (time= yr = begin of integration interval)----

               i=co2var.ipos_data

       

                 co2var.emiss_now = interpol_bern(co2var.yr,par.tvet(i),par.tvet(i+1),dat.eco2_tot(i),dat.eco2_tot(i+1))


!     ---- output at time = yr = begin of integration interval:

!                                for pco2a,pco2s,conc_oc, emission, and fluxes


               if (co2var.flagout) then

                       co2var.flagout = .false.

                       mout.st_co2(co2var.ipos_data) = co2var.pco2a_now

       endif


!   ----- ocean surface layer concentration and pco2 at time=end of integration interval


       co2var.fs = 0.0d0

       

       do i=1,n

           co2var.fs = co2var.fs+co2var.fnet_air_sea_vec(i)*co2var.respoc_vec(n+1-i)

       enddo


!                dissolved inorganic carbon perturbation in [umol/kg]

       co2var.conc_oc =co2var.fs*co2var.dt*1.d6/(co2par.aoc*co2par.hs*co2par.density)


       co2var.pco2s=co2var.conc_oc*zeta_factor(co2par.temper,co2var.conc_oc)+co2par.pco2a0


!     --- atmospheri! pco2 at time=end of integration interval:


           co2var.pco2a_now = co2var.pco2a_now +( co2var.emiss_now - co2var.nep - co2var.fnet_air_sea_vec(n) ) * co2par.ppm_per_mol*co2var.dt

           

!        call callBack(n,pco2a_now)

           

             enddo


       return

end subroutine co2bern_model

   



subroutine evaluateModelCH4 ()

implicit none

real*8 :: Fch4  

real*8 :: km  

!

! ch4

if (trim(par.use_ch4)== 'calc') then

   if (trim(par.ch4_ch4preindchoice)== 'power') then

       now.td_ch4time= par.ch4_ch4preind*(par.ch4_ch4pre/now.st_ch4)**-0.12 !(f.27)

   end if

   

   select case ( trim(par.ch4_ch4forchoice)) !

       case ( 'total') !

           Fch4= par.ch4_rCH4*now.td_ECH4

       case ( 'anthro') !

           km=1 !

           Fch4= par.ch4_rCH4*now.td_ECH4 + km*now.st_cb1 !(f.29)

       case ( 'perturbation') !

           Fch4= par.ch4_rCH4*now.td_ECH4 + par.ch4_ch4pre/now.td_ch4time !(f.30)

!        case default

!        disp('Opps&, Unexpected par.ch4_ch4preindchoice option at evaluateModelCh4')

   end select

   now.stDer_ch4= Fch4 - now.st_ch4/now.td_ch4time ! (f.28)

   ! Calc ! oxidation of CH4 to the atmospheric CO2

   now.td_fluxch4co2= par.ch4_sfch4co2*par.ch4_rCH4CO2*now.st_ch4/now.td_ch4time !(f.31)

end if

end subroutine evaluateModelCH4




subroutine evaluateModelN2O ()

implicit none

! Variable declarations

real*8 :: Fn2o    


! n2o

if (trim(par.use_n2o)== 'calc') then

   select case ( trim(par.n2o_n2oforchoice)) !

     case ( 'total') !

           Fn2o= par.n2o_rN20*now.td_EN2O

      case ( 'anthro', 'perturbation') !

           Fn2o= par.n2o_rN20*now.td_EN2O + par.n2o_n2opre/par.n2o_n2otime !(f.33, f.34)

      !case default

      !      disp('Opps&, Unexpected par.n2o_n2oforchoice option at evaluateModelCh4')

   end select

   

   now.stDer_n2o= Fn2o - now.st_n2o/par.n2o_n2otime ! (f.34)

end if

end subroutine evaluateModelN2O


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!>

subroutine evaluateModelOGAS ()

implicit none

integer ::i

! (f.35)

do i=1,par.Nogas

   if (par.use_ogas(i)== 'calc') then

       now.stDer_ogas(i)= par.ogas_conv(i)*now.td_EOGAS(i) - &

       (now.st_ogas(i)+par.ogas_pre(i))/par.ogas_time(i)

   end if

end do

end subroutine evaluateModelOGAS



subroutine evaluateModelRF()

implicit none

real*8  :: Fco2

integer :: i


! all RF calculations

select case ( trim(par.use_rfco2) ) !! RF CO2 (f.36)

   case ( 'calc', 'ccuser') !

   now.td_rfco2= par.rf_fco2*log(now.st_co2/par.co2_co2pre) !(f.36)

end select


select case ( trim(par.use_rfch4) ) !! RF CH4 (f.37)

   case ( 'calc', 'ccuser') !

   now.td_rfch4= par.rf_ch4radeff*(sqrt(now.st_ch4)-sqrt(par.ch4_ch4pre)) &

   - foverlap(now.st_ch4, par.n2o_n2opre) &

   + foverlap(par.ch4_ch4pre, par.n2o_n2opre) !(f.37)

end select


select case ( trim(par.use_rfn2o) ) !! RF N2O (f.38)

   case ( 'calc', 'ccuser' ) !! Fórmula do artigo está errada

   now.td_rfn2o= par.rf_n2oradeff*(sqrt(now.st_n2o)-sqrt(par.n2o_n2opre)) &

   - foverlap(par.ch4_ch4pre, now.st_n2o) &

   + foverlap(par.ch4_ch4pre, par.n2o_n2opre) ! Fórmula do artigo está errada !(f.38)

end select


! OGAS

do i=1,par.Nogas

   if (trim(par.use_ogas(i))== 'calc' .or. trim(par.use_ogas(i))== 'ccuser') then

       now.td_rfogas(i)= 1e-3*par.ogas_radeff(i)*(now.st_ogas(i)+par.ogas_pre(i)) !(f.40)

   end if

end do


! AERO - Table 3 and f.42

if ( trim(par.use_rfaero)== 'calc') then

   select case (trim(par.rf_aeromodelchoice)) !

       case ( 'none') !

           par.rf_aeroval= 0

           par.rf_aeroprop= 0

           Fco2= 0

       case ( 'const') !

           par.rf_aeroprop= 0

           Fco2= 0

       case ( 'prop', 'cleanup') !

           par.rf_aeroval= 0

           Fco2= now.td_ECO2_foss

       case ( 'all') !

           Fco2= now.td_ECO2_foss

       case ( 'allluc') !

           Fco2= now.td_ECO2_foss + now.td_ECO2_luc

   end select

       

   now.td_rfaero= par.rf_aeroval + par.rf_aeroprop*Fco2* &

   min(1.0D0, exp(-(now.time-2000)/par.rf_aerotau)) ! (f.42)

   

end if


! total RF - (f.41)

if (trim(par.use_rftotal)== 'sum') then

   now.td_rftotal= par.rf_sfco2*now.td_rfco2 + par.rf_sfch4*now.td_rfch4 + &

   par.rf_sfn2o*now.td_rfn2o + par.rf_sfogas*sum(now.td_rfogas) + &

   par.rf_sfvolc*now.td_rfvolc + par.rf_sfaero*now.td_rfaero + par.rf_sfother*now.td_rfother !(f.41)

end if


end subroutine evaluateModelRF





subroutine evaluateModelTEMP ()

implicit none

! Use somemodule

! Arguments declarations

if (trim(par.use_temp)== 'calc') then

   now.stDer_temp1= par.temp_CsensT*par.temp_TP(2)/par.temp_TP(3)*now.td_rftotal - &

   now.st_temp1/par.temp_TP(3) ! (f.43)

   now.stDer_temp2= par.temp_CsensT*par.temp_TP(4)/par.temp_TP(5)*now.td_rftotal - &

   now.st_temp2/par.temp_TP(5) ! (f.43)

   now.stDer_temp3= par.temp_CsensT*par.temp_TP(6)/par.temp_TP(7)*now.td_rftotal - &

   now.st_temp3/par.temp_TP(7) ! (f.43)

   now.td_temp= now.st_temp1 + now.st_temp2 + now.st_temp3 + &

   par.temp_temppre ! (f.47)

end if

end subroutine evaluateModelTEMP





!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!>

subroutine evaluateModelSLR()

implicit none

! Use somemodule

! Arguments declarations

! Variable declarations

real*8 :: dT_dt  

real*8 :: fact  

!

! SLR

if (trim(par.use_slr)=='calc') then

   ! accc

   now.td_slra= now.st_slr1 + now.st_slr2

   fact= par.slr_slreq/par.slr_rfeq

   now.stDer_slr1= now.td_rftotal*par.slr_a1*fact - now.st_slr1/par.slr_t1

   now.stDer_slr2= now.td_rftotal*par.slr_a2*fact - now.st_slr2/par.slr_t2


! ‘vermeer’

if (trim(par.use_temp)=='calc') then

   dT_dt= now.stDer_temp1 + now.stDer_temp2 + now.stDer_temp3

else

   dT_dt= now.dTempdt

end if

now.stDer_slrv= par.slr_av*(now.td_temp - par.temp_temppre) + par.slr_bv*dT_dt

! ’jevrejeva’

now.stDer_slrj= (par.slr_aj*now.td_rftotal + par.slr_bj - now.st_slrj)/par.slr_tj

end if

end subroutine evaluateModelSLR !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!