Add a time-dependent function for CO2 (Update PR #1173)

src/core_atmosphere/physics/Makefile

@@ -46,6 +46,7 @@ OBJS = \
 	mpas_atmphys_manager.o             \
 	mpas_atmphys_o3climatology.o       \
 	mpas_atmphys_packages.o            \
+	mpas_atmphys_radiation_utils.o     \
 	mpas_atmphys_rrtmg_lwinit.o        \
 	mpas_atmphys_rrtmg_swinit.o        \
 	mpas_atmphys_sfc_diagnostics.o     \
@@ -150,18 +151,22 @@ mpas_atmphys_driver_pbl.o: \
 	mpas_atmphys_constants.o \
 	mpas_atmphys_vars.o
 
+mpas_atmphys_radiation_utils.o: \
+	mpas_atmphys_constants.o
+
 mpas_atmphys_driver_radiation_lw.o: \
 	mpas_atmphys_camrad_init.o \
 	mpas_atmphys_constants.o \
-	mpas_atmphys_driver_radiation_sw.o \
 	mpas_atmphys_manager.o \
+	mpas_atmphys_radiation_utils.o     \
 	mpas_atmphys_rrtmg_lwinit.o \
 	mpas_atmphys_vars.o
 
 mpas_atmphys_driver_radiation_sw.o: \
 	mpas_atmphys_camrad_init.o \
 	mpas_atmphys_constants.o \
 	mpas_atmphys_manager.o \
+	mpas_atmphys_radiation_utils.o     \
 	mpas_atmphys_rrtmg_swinit.o \
 	mpas_atmphys_vars.o
 

src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F

@@ -11,12 +11,12 @@ module mpas_atmphys_driver_radiation_lw
  use mpas_pool_routines
  use mpas_timer,only: mpas_timer_start,mpas_timer_stop
 
- use mpas_atmphys_driver_radiation_sw, only: radconst
  use mpas_atmphys_constants
  use mpas_atmphys_manager, only: gmt,curr_julday,julday,year
  use mpas_atmphys_camrad_init
  use mpas_atmphys_rrtmg_lwinit
  use mpas_atmphys_vars
+ use mpas_atmphys_radiation_utils, only: radconst
 
 !wrf physics:
  use module_ra_cam            
@@ -841,6 +841,13 @@ subroutine driver_radiation_lw(xtime_s,configs,mesh,state,time_lev,diag_physics,
 !copy MPAS arrays to local arrays:
  call radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physics,atm_input,sfc_input,its,ite)
 
+! This should be OMP MASTER with barrier afterwards or OMP SINGLE, since declin and solcon
+! are global variables in mpas_atmphys_vars.F and race conditions may occur otherwise!
+!$OMP SINGLE
+!... calculates solar declination and co2 concentration:
+ call radconst(declin,solcon,curr_julday,degrad,dpd,co2_val,year)
+!$OMP END SINGLE
+
 !call to longwave radiation scheme:
  radiation_lw_select: select case (trim(radt_lw_scheme))
     case ("rrtmg_lw")
@@ -878,7 +885,7 @@ subroutine driver_radiation_lw(xtime_s,configs,mesh,state,time_lev,diag_physics,
             rre_ice  = rreice_p    , rre_snow   = rresnow_p     , lwupt     = lwupt_p    , &
             lwuptc   = lwuptc_p    , lwdnt      = lwdnt_p       , lwdntc    = lwdntc_p   , &
             lwupb    = lwupb_p     , lwupbc     = lwupbc_p      , lwdnb     = lwdnb_p    , &
-            lwdnbc   = lwdnbc_p    ,                                                       &
+            lwdnbc   = lwdnbc_p    , co2        = co2_val       ,                          &
             ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde ,        &
             ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme ,        &
             its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte          &
@@ -888,13 +895,6 @@ subroutine driver_radiation_lw(xtime_s,configs,mesh,state,time_lev,diag_physics,
     case ("cam_lw")
        xtime_m = xtime_s/60.
 
-! This should be OMP MASTER with barrier afterwards or OMP SINGLE, since declin and solcon
-! are global variables in mpas_atmphys_vars.F and race conditions may occur otherwise!
-!$OMP SINGLE
-       !... calculates solar declination:
-       call radconst(declin,solcon,curr_julday,degrad,dpd)
-!$OMP END SINGLE
-
 !... convert the radiation time_step to minutes:
        radt = dt_radtlw/60.
 

src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F

@@ -16,6 +16,7 @@ module mpas_atmphys_driver_radiation_sw
  use mpas_atmphys_camrad_init
  use mpas_atmphys_rrtmg_swinit
  use mpas_atmphys_vars
+ use mpas_atmphys_radiation_utils, only: radconst
 
 !wrf physics:
  use module_mp_thompson_aerosols
@@ -28,8 +29,7 @@ module mpas_atmphys_driver_radiation_sw
  public:: allocate_radiation_sw,   &
           deallocate_radiation_sw, &
           driver_radiation_sw,     &
-          init_radiation_sw,       &
-          radconst
+          init_radiation_sw
 
 !MPAS driver for parameterization of shortwave radiation codes.
 !Laura D. Fowler (send comments to laura@ucar.edu).
@@ -43,7 +43,6 @@ module mpas_atmphys_driver_radiation_sw
 ! driver_radiation_sw    : main driver (called from subroutine physics_driver).
 ! radiation_sw_from_MPAS : initialize local arrays.
 ! radiation_sw_to_MPAS   : copy local arrays to MPAS arrays.
-! radconst               : calculate solar declination,...
 !
 ! WRF physics called from driver_radiation_sw:
 ! --------------------------------------------
@@ -879,9 +878,10 @@ subroutine driver_radiation_sw(itimestep,configs,mesh,state,time_lev,diag_physic
 ! This should be OMP MASTER with barrier afterwards or OMP SINGLE, since declin and solcon
 ! are global variables in mpas_atmphys_vars.F and race conditions may occur otherwise!
 !$OMP SINGLE
-!... calculates solar declination:
+!... calculates solar declination and co2 concentration:
 !call radconst(declin,solcon,julday,degrad,dpd)
- call radconst(declin,solcon,curr_julday,degrad,dpd)
+ call radconst(declin,solcon,curr_julday,degrad,dpd,co2_val,year)
+!call mpas_log_write(' --- driver_radiation_sw after radconst')
 !call mpas_log_write('     ITIMESTEP   = $i', intArgs=(/itimestep/))
 !call mpas_log_write('     YEAR        = $i', intArgs=(/year/))
 !call mpas_log_write('     JULDAY      = $i', intArgs=(/julday/))
@@ -890,6 +890,7 @@ subroutine driver_radiation_sw(itimestep,configs,mesh,state,time_lev,diag_physic
 !call mpas_log_write('     CURR_JULDAY = $r', realArgs=(/curr_julday/))
 !call mpas_log_write('     SOLCON      = $r', realArgs=(/solcon/))
 !call mpas_log_write('     DECLIN      = $r', realArgs=(/declin/))
+!call mpas_log_write('     CO2_VAL     = $r', realArgs=(/co2_val/))
 !$OMP END SINGLE
 
 !... convert the radiation time_step to minutes:
@@ -919,8 +920,6 @@ subroutine driver_radiation_sw(itimestep,configs,mesh,state,time_lev,diag_physic
        call rrtmg_swrad( &
               p3d        = pres_hyd_p    , p8w      = pres2_hyd_p , pi3d     = pi_p           , &
               t3d        = t_p           , t8w      = t2_p        , dz8w     = dz_p           , &
-!             qv3d       = qv_p          , qc3d     = qc_p        , qi3d     = qi_p           , &
-!             qs3d       = qs_p          , cldfra3d = cldfrac_p   , tsk      = tsk_p          , &
               qv3d       = qvrad_p       , qc3d     = qcrad_p     , qi3d     = qirad_p        , &
               qs3d       = qsrad_p       , cldfra3d = cldfrac_p   , tsk      = tsk_p          , &
               albedo     = sfc_albedo_p  , xland    = xland_p     , xice     = xice_p         , &
@@ -938,6 +937,7 @@ subroutine driver_radiation_sw(itimestep,configs,mesh,state,time_lev,diag_physic
               swdnt      = swdnt_p       , swdntc   = swdntc_p    , swupb      = swupb_p      , &
               swupbc     = swupbc_p      , swdnb    = swdnb_p     , swdnbc     = swdnbc_p     , &
               swddir     = swddir_p      , swddni   = swddni_p    , swddif     = swddif_p     , &
+              co2        = co2_val       ,                                                      &
               ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde ,           &
               ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme ,           &
               its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte             &
@@ -1005,50 +1005,6 @@ subroutine driver_radiation_sw(itimestep,configs,mesh,state,time_lev,diag_physic
 
  end subroutine driver_radiation_sw
 
-!=================================================================================================================
- subroutine radconst(declin,solcon,julian,degrad,dpd)
-!=================================================================================================================
-
-!input arguments:
-!integer,intent(in):: julian
- real(kind=RKIND),intent(in):: julian
- real(kind=RKIND),intent(in):: degrad,dpd
-
-!output arguments:
- real(kind=RKIND),intent(out):: declin,solcon
-
-!local variables:
- real(kind=RKIND):: obecl,sinob,sxlong,arg,decdeg,djul,rjul,eccfac
-
-!-----------------------------------------------------------------------------------------------------------------
-
-   declin=0.
-   solcon=0.
-
-!obecl : obliquity = 23.5 degree.
-
-   obecl=23.5*degrad
-   sinob=sin(obecl)
-
-!calculate longitude of the sun from vernal equinox:        
-
-   if(julian.ge.80.)sxlong=dpd*(julian-80.)
-   if(julian.lt.80.)sxlong=dpd*(julian+285.)
-   sxlong=sxlong*degrad
-   arg=sinob*sin(sxlong)
-   declin=asin(arg)
-   decdeg=declin/degrad
-
-!solar constant eccentricity factor (paltridge and platt 1976)
-
-   djul=julian*360./365.
-   rjul=djul*degrad
-   eccfac=1.000110+0.034221*cos(rjul)+0.001280*sin(rjul)+0.000719*  &
-          cos(2*rjul)+0.000077*sin(2*rjul)
-   solcon=solcon_0*eccfac
-
- end subroutine radconst
-
 !=================================================================================================================
  end module mpas_atmphys_driver_radiation_sw
 !=================================================================================================================

src/core_atmosphere/physics/mpas_atmphys_radiation_utils.F

@@ -0,0 +1,107 @@
+! Copyright (c) 2013,  Los Alamos National Security, LLC (LANS)
+! and the University Corporation for Atmospheric Research (UCAR).
+!
+! Unless noted otherwise source code is licensed under the BSD license.
+! Additional copyright and license information can be found in the LICENSE file
+! distributed with this code, or at http://mpas-dev.github.com/license.html
+!
+!=================================================================================================================
+ module mpas_atmphys_radiation_utils
+
+ use mpas_kind_types, only : RKIND
+ use mpas_atmphys_constants, only : solcon_0
+ use mpas_log, only : mpas_log_write
+
+ implicit none
+ private
+ public:: radconst,     &
+          co2_estimate
+
+ contains
+
+!=================================================================================================================
+ subroutine radconst(declin,solcon,julian,degrad,dpd,co2,year)
+!=================================================================================================================
+!
+! computes solar declination
+! computes co2 concentration from WRF v4.2 (Jimy Dudhia)
+!
+!input arguments:
+ real(kind=RKIND),intent(in):: julian
+ real(kind=RKIND),intent(in):: degrad,dpd
+ integer, intent(in):: year
+
+!output arguments:
+ real(kind=RKIND),intent(out):: declin,solcon
+ real(kind=RKIND),intent(out):: co2
+
+!local variables:
+ real(kind=RKIND):: obecl,sinob,sxlong,arg,decdeg,djul,rjul,eccfac
+
+!-----------------------------------------------------------------------------------------------------------------
+
+ call mpas_log_write('--- enter subroutine radconst:')
+
+ declin=0.
+ solcon=0.
+
+!obecl : obliquity = 23.5 degree.
+
+ obecl=23.5*degrad
+ sinob=sin(obecl)
+
+!calculate longitude of the sun from vernal equinox:
+
+ if(julian.ge.80.)sxlong=dpd*(julian-80.)
+ if(julian.lt.80.)sxlong=dpd*(julian+285.)
+
+ sxlong=sxlong*degrad
+ arg=sinob*sin(sxlong)
+ declin=asin(arg)
+ decdeg=declin/degrad
+
+!solar constant eccentricity factor (paltridge and platt 1976)
+
+ djul=julian*360./365.
+ rjul=djul*degrad
+ eccfac=1.000110+0.034221*cos(rjul)+0.001280*sin(rjul)+0.000719*  &
+        cos(2*rjul)+0.000077*sin(2*rjul)
+ solcon=solcon_0*eccfac
+
+!estimate co2 concentration from an equation from WRF4.2
+
+ co2 = co2_estimate(year,julian)
+
+ end subroutine radconst
+
+!=================================================================================================================
+ real(kind=RKIND) function co2_estimate(year,julian)
+!=================================================================================================================
+
+! This function calculates CO2 from an equation in WRF by J. Dudhia
+
+ implicit none
+
+ integer, intent(in) :: year
+ real(kind=RKIND),intent(in):: julian
+
+ real(kind=RKIND) :: rdays,ryear
+ integer :: ny1, ny2, ny3
+
+!-----------------------------------------------------------------------------------------------------------------
+
+ rdays = 365._RKIND
+ ny1 = mod(year,4)
+ ny2 = mod(year,100)
+ ny3 = mod(year,400)
+ if( (ny1.eq.0 .and. ny2.ne.0) .or. ny3.eq.0 ) rdays = 366._RKIND
+ ryear = float(year)+julian/rdays
+
+! Annual function for co2 from WRF v4.2
+ co2_estimate = (280._RKIND + 90._RKIND*exp(0.02_RKIND*(ryear-2000._RKIND)))*1.e-6_RKIND
+
+ end function co2_estimate
+
+!=================================================================================================================
+ end module mpas_atmphys_radiation_utils
+!=================================================================================================================

src/core_atmosphere/physics/mpas_atmphys_vars.F

@@ -687,6 +687,13 @@ module mpas_atmphys_vars
  real(kind=RKIND),dimension(:,:,:,:),allocatable:: ssaaer_p   !aerosol single scatterin albedo in RRTMG SW     [-]
  real(kind=RKIND),dimension(:,:,:,:),allocatable:: asyaer_p   !aerosol asymmetry factor in RRTMG SW            [-]
 
+!=================================================================================================================
+!... variables and arrays related to parameterizations of long- and short-wave radiation, especially RRTMG:
+!=================================================================================================================
+
+ real(kind=RKIND):: &
+    co2_val            !co2 concentration                                                       [-]
+
 !=================================================================================================================
 !... variables and arrays related to parameterization of short-wave radiation:
 !=================================================================================================================

src/core_atmosphere/physics/physics_wrf/module_ra_rrtmg_lw.F

@@ -11582,6 +11582,7 @@ subroutine rrtmg_lwrad( &
                        lwupt,lwuptc,lwdnt,lwdntc,                 &
                        lwupb,lwupbc,lwdnb,lwdnbc,                 &
                        lwupflx, lwupflxc, lwdnflx, lwdnflxc,      &
+                       co2,                                       &
                        ids,ide, jds,jde, kds,kde,                 & 
                        ims,ime, jms,jme, kms,kme,                 &
                        its,ite, jts,jte, kts,kte                  &
@@ -11671,7 +11672,7 @@ subroutine rrtmg_lwrad( &
 !--- set trace gas volume mixing ratios, 2005 values, IPCC (2007):
 !carbon dioxide (379 ppmv)
  real :: co2
- data co2 / 379.e-6 / 
+!data co2 / 379.e-6 / 
 !methane (1774 ppbv)
  real :: ch4
  data ch4 / 1774.e-9 / 

src/core_atmosphere/physics/physics_wrf/module_ra_rrtmg_sw.F

@@ -10037,6 +10037,7 @@ subroutine rrtmg_swrad( &
                        swupb,swupbc,swdnb,swdnbc,                      &
                        swupflx, swupflxc, swdnflx, swdnflxc,           &
                        swddir,swddni,swddif,                           &
+                       co2,                                            &
                        ids,ide, jds,jde, kds,kde,                      & 
                        ims,ime, jms,jme, kms,kme,                      &
                        its,ite, jts,jte, kts,kte                       &
@@ -10142,7 +10143,7 @@ subroutine rrtmg_swrad( &
 !--- set trace gas volume mixing ratios, 2005 values, IPCC (2007):
 !carbon dioxide (379 ppmv)
  real :: co2
- data co2 / 379.e-6 / 
+!data co2 / 379.e-6 / 
 !methane (1774 ppbv)
  real :: ch4
  data ch4 / 1774.e-9 /