diff --git a/src/core_atmosphere/physics/Makefile b/src/core_atmosphere/physics/Makefile index ac5ff5d6f3..55ad80e785 100644 --- a/src/core_atmosphere/physics/Makefile +++ b/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,11 +151,14 @@ 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 @@ -162,6 +166,7 @@ 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 diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F index dfa327fea8..95ab6f212d 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F +++ b/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. diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F index 3a14275826..61f8a628af 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F +++ b/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 !================================================================================================================= diff --git a/src/core_atmosphere/physics/mpas_atmphys_radiation_utils.F b/src/core_atmosphere/physics/mpas_atmphys_radiation_utils.F new file mode 100644 index 0000000000..aed506bf9b --- /dev/null +++ b/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 +!================================================================================================================= diff --git a/src/core_atmosphere/physics/mpas_atmphys_vars.F b/src/core_atmosphere/physics/mpas_atmphys_vars.F index 134009f537..30d0b101dd 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_vars.F +++ b/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: !================================================================================================================= diff --git a/src/core_atmosphere/physics/physics_wrf/module_ra_rrtmg_lw.F b/src/core_atmosphere/physics/physics_wrf/module_ra_rrtmg_lw.F index 89bc6b00b1..12e3a33100 100644 --- a/src/core_atmosphere/physics/physics_wrf/module_ra_rrtmg_lw.F +++ b/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 / diff --git a/src/core_atmosphere/physics/physics_wrf/module_ra_rrtmg_sw.F b/src/core_atmosphere/physics/physics_wrf/module_ra_rrtmg_sw.F index faa5761c9a..a0bef9afed 100644 --- a/src/core_atmosphere/physics/physics_wrf/module_ra_rrtmg_sw.F +++ b/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 /