!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS 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 FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module grey_radiation_mod ! ================================================================================== ! ================================================================================== ! use utilities_mod, only: error_mesg, open_file, file_exist, & ! check_nml_error, FATAL, get_my_pe, & ! close_file, set_domain, get_num_pes, & ! read_data, write_data, & ! check_system_clock, NOTE, & ! get_domain_decomp, check_system_clock use fms_mod, only: open_namelist_file, check_nml_error, & mpp_pe, mpp_root_pe, close_file, & write_version_number, stdlog use constants_mod, only: stefan, cp_air, grav use astronomy_mod, only: astronomy_init, daily_mean_solar, diurnal_solar use diag_manager_mod, only: register_diag_field, send_data use time_manager_mod, only: time_type, set_date, set_time, & get_time, operator(+), & operator(-), operator(/=), get_date !================================================================================== implicit none private !================================================================================== ! version information character(len=128), parameter :: version = & '$Id: grey_radiation.F90,v 17.0.4.2 2009/10/29 22:30:01 nnz Exp $' character(len=128), parameter :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' logical :: module_is_initialized = .false. !================================================================================== ! public interfaces public :: grey_radiation_init, grey_radiation, grey_radiation_end !================================================================================== ! module variables character (len=*), parameter :: module='grey_radiation_mod' logical :: initialized =.false. real, parameter :: p00 = 1000.e2 real :: solar_constant = 1360.0 real :: del_sol = 0.0 ! modif omp: winter/summer hemisphere real :: del_sw = 0.0 real :: ir_tau_eq = 4.0 real :: ir_tau_pole = 4.0 real :: atm_abs = 0.2 real :: sw_diff = 0.0 real :: long_pert = 180. real :: del_long = 30. real :: size_pert = 0. real :: linear_tau = 0.1 real :: lat_pert = 0.0 real :: lon_pert = 180.0 real :: del_lat = 30.0 real :: del_lon = 90.0 real :: fcng_pert = 0.0 real :: wave_amp = 0.0 real :: wave_lon = 180.0 real :: wave_lat = 0.0 real :: wave_del_lon = 30.0 real :: wave_del_lat = 20.0 real :: wave_period = 20.0 real :: wave_env = 80.0 logical :: wave_source = .FALSE. integer :: n_tau = 4 ! call to astronomy to include the Seasonal Cycle logical :: do_season = .false. real, save :: pi, deg_to_rad , rad_to_deg namelist/grey_radiation_nml/ solar_constant, del_sol, & ir_tau_eq, ir_tau_pole, atm_abs, sw_diff, long_pert, del_long, & size_pert, linear_tau, del_sw, & lat_pert, lon_pert, del_lat, del_lon, fcng_pert, & wave_amp, wave_lon, wave_lat, wave_del_lon, & wave_del_lat, wave_period, wave_env, wave_source, do_season !================================================================================== !-------------------- diagnostics fields ------------------------------- integer :: id_olr, id_swdn_sfc, id_swdn_toa, id_lwdn_sfc, id_lwup_sfc, & id_tdt_rad, id_flux_rad, id_flux_lw, id_flux_sw, id_entrop_rad, & id_tsurfgrey, id_tempgrey character(len=14), parameter :: mod_name = 'grey_radiation' real :: missing_value = -999. contains ! ================================================================================== ! ================================================================================== subroutine grey_radiation_init(axes, Time) !------------------------------------------------------------------------------------- integer, intent(in), dimension(4) :: axes type(time_type), intent(in) :: Time !------------------------------------------------------------------------------------- integer, dimension(3) :: half = (/1,2,4/) integer :: ierr, io, unit, logunit !----------------------------------------------------------------------------------------- ! read namelist and copy to logfile unit = open_namelist_file ( ) ierr=1 do while (ierr /= 0) read (unit, nml=grey_radiation_nml, iostat=io, end=10) ierr = check_nml_error (io, 'grey_radiation_nml') enddo 10 call close_file (unit) call write_version_number ( version, tagname ) if ( mpp_pe() == mpp_root_pe() ) then logunit = stdlog() write (logunit, nml=grey_radiation_nml) endif call close_file (unit) pi = 4.0*atan(1.) deg_to_rad = 2.*pi/360. rad_to_deg = 360.0/2./pi call astronomy_init initialized = .true. !----------------------------------------------------------------------- !------- initialize quantities for integral package ------- ! call diag_integral_field_init ('olr', 'f8.3') ! call diag_integral_field_init ('abs_sw', 'f8.3') !----------------------------------------------------------------------- !------------ initialize diagnostic fields --------------- id_olr = & register_diag_field ( mod_name, 'olr', axes(1:2), Time, & 'outgoing longwave radiation', & 'watts/m2', missing_value=missing_value ) id_swdn_sfc = & register_diag_field ( mod_name, 'swdn_sfc', axes(1:2), Time, & 'SW flux down at surface', & 'watts/m2', missing_value=missing_value ) id_swdn_toa = & register_diag_field ( mod_name, 'swdn_toa', axes(1:2), Time, & 'SW flux down at TOA', & 'watts/m2', missing_value=missing_value ) id_lwup_sfc = & register_diag_field ( mod_name, 'lwup_sfc', axes(1:2), Time, & 'LW flux up at surface', & 'watts/m2', missing_value=missing_value ) id_lwdn_sfc = & register_diag_field ( mod_name, 'lwdn_sfc', axes(1:2), Time, & 'LW flux down at surface', & 'watts/m2', missing_value=missing_value ) id_tdt_rad = & register_diag_field ( mod_name, 'tdt_rad', axes(1:3), Time, & 'Temperature tendency due to radiation', & 'K/s', missing_value=missing_value ) id_flux_rad = & register_diag_field ( mod_name, 'flux_rad', axes(half), Time, & 'Total radiative flux (positive up)', & 'W/m^2', missing_value=missing_value ) id_flux_lw = & register_diag_field ( mod_name, 'flux_lw', axes(half), Time, & 'Net longwave radiative flux (positive up)', & 'W/m^2', missing_value=missing_value ) id_flux_sw = & register_diag_field ( mod_name, 'flux_sw', axes(half), Time, & 'Net shortwave radiative flux (positive up)', & 'W/m^2', missing_value=missing_value ) id_entrop_rad = & register_diag_field ( mod_name, 'entrop_rad', axes(1:3), Time, & 'Entropy production by radiation', & '1/s', missing_value=missing_value ) ! rif:(09/09/09) New Diagnostics id_tsurfgrey = & register_diag_field ( mod_name, 't_surf_rad', axes(1:2), Time, & 'Temp surf from grey radiation', & 'K', missing_value=missing_value ) id_tempgrey = & register_diag_field ( mod_name, 'temp_rad', axes(1:3), Time, & 'Temperature from grey radiation', & 'K', missing_value=missing_value ) module_is_initialized = .true. return end subroutine grey_radiation_init ! ================================================================================== subroutine grey_radiation (is, js, Time, Time_diag, lat, lon, phalfgrey, albedo, t_surf, t, tdt, net_surf_sw_down, surf_lw_down) integer, intent(in) :: is, js type(time_type), intent(in) :: Time, Time_diag real, intent(in) , dimension(:,:) :: lat, lon, albedo real, intent(in) , dimension(:,:) :: t_surf real, intent(in) , dimension(:,:,:) :: t real, intent(in) , dimension(:,:,:) :: phalfgrey real, intent(inout), dimension(:,:,:) :: tdt real, intent(out), dimension(:,:) :: net_surf_sw_down, surf_lw_down real, dimension(size(t,2)) :: ss, ss2, ss4, ss6, ss8, solar, tau_0, solar_tau_0, p2 real, dimension(size(t,1), size(t,2)) :: b_surf real, dimension(size(t,1), size(t,2), size(t,3)) :: b, tdt_rad, entrop_rad real, dimension(size(t,1), size(t,2), size(t,3)+1) :: up, down, net, solar_down, flux_rad, flux_sw real, dimension(size(t,2), size(t,3) ) :: dtrans real, dimension(size(t,2), size(t,3)+1) :: tau, solar_tau real, dimension(size(t,1), size(t,2)) :: long_forcing, olr, swin real, dimension(size(t,1), size(t,2)) :: walker_forcing, wave_forc integer :: i, j, k, n real :: surf_albedo, cosz1, fracday1, rrsun1 real :: dist logical :: used n = size(t,3) do i=1,size(t,1) do j=1,size(t,2) long_forcing(i,j) = size_pert*exp(-(rad_to_deg*lon(i,j)-long_pert)**2./ & del_long**2.) end do end do do i = 1,size(t,1) do j = 1,size(t,2) dist = SQRT(((rad_to_deg*lat(i,j) - lat_pert)/del_lat)**2 & + ((rad_to_deg*lon(i,j) - lon_pert)/del_lon)**2) if (dist .lt. 1) then walker_forcing(i,j) = fcng_pert*cos(pi * dist*0.5)**2; else walker_forcing(i,j) = 0.0 end if end do end do ss = sin(lat(1,:)) ss2 = ss*ss p2 = (1. - 3.*ss*ss)/4. solar = 0.25*solar_constant*(1.0 + del_sol*p2 + del_sw * ss) tau_0 = ir_tau_eq +(ir_tau_pole - ir_tau_eq)*ss*ss solar_tau_0 = (1.0 - sw_diff*ss*ss)*atm_abs b = stefan*t*t*t*t b_surf = stefan*t_surf*t_surf*t_surf*t_surf do k = 1, n+1 ! modif df 1-23-04: changing profile of IR absorber ! modif rif 9-05-08: changed p_half to phalf for Bgrid model tau(:,k) = tau_0(:) * (linear_tau * phalfgrey(1,1,k)/p00 + (1.0 - linear_tau) & * (phalfgrey(1,1,k)/p00)**4) solar_tau(:,k) = solar_tau_0(:)*(phalfgrey(1,1,k)/p00)**4 end do do k = 1, n dtrans(:,k) = exp(-(tau(:,k+1)-tau(:,k))) end do up(:,:,n+1) = b_surf do k = n,1,-1 do j = 1, size(t,2) up(:,j,k) = up(:,j,k+1)*dtrans(j,k) + b(:,j,k)*(1.0 - dtrans(j,k)) end do end do down(:,:,1) = 0.0 do k = 1,n do j =1, size(t,2) down(:,j,k+1) = down(:,j,k)*dtrans(j,k) + b(:,j,k)*(1.0 - dtrans(j,k)) end do end do if (do_season) then !Seasonal Cycle do j=1,size(t,2) call daily_mean_solar(lat(1,j),time,cosz1,fracday1,rrsun1) do i=1,size(t,1) do k=1,n+1 solar_down(i,j,k) = cosz1*fracday1*solar_constant*rrsun1 end do end do end do else do i = 1, size(t,1) do j = 1, size(t,2) do k = 1,n+1 solar_down(i,j,k) = (walker_forcing(i,j)+long_forcing(i,j) & + solar(j))*exp(-solar_tau(j,k)) end do end do end do end if !if (do_season) do k = 1,n+1 net(:,:,k) = up(:,:,k)-down(:,:,k) flux_sw(:,:,k) = albedo(:,:)*solar_down(:,:,n+1) - solar_down(:,:,k) flux_rad(:,:,k) = net(:,:,k) + flux_sw(:,:,k) end do do k = 1,n tdt_rad(:,:,k) = (net(:,:,k+1) - net(:,:,k) - solar_down(:,:,k+1) + solar_down(:,:,k)) & *grav/(cp_air*(phalfgrey(1,1,k+1)-phalfgrey(1,1,k))) tdt(:,:,k) = tdt(:,:,k) + tdt_rad(:,:,k) end do surf_lw_down = down(:,:,n+1) net_surf_sw_down = solar_down(:,:,n+1)*(1. - albedo(:,:)) olr = up(:,:,1) swin = solar_down(:,:,1) !------- t_surf grey (t_surf_greyrad) ------- if ( id_tsurfgrey > 0 ) then used = send_data ( id_tsurfgrey, t_surf, Time_diag, is, js ) endif !------- temp grey (temp_greyrad) ------- if ( id_tempgrey > 0 ) then used = send_data ( id_tempgrey, t, Time_diag, is, js, 1 ) endif !------- outgoing lw flux toa (olr) ------- if ( id_olr > 0 ) then used = send_data ( id_olr, olr, Time_diag, is, js ) endif !------- downward sw flux surface ------- if ( id_swdn_sfc > 0 ) then used = send_data ( id_swdn_sfc, net_surf_sw_down, Time_diag, is, js ) endif !------- incoming sw flux toa ------- if ( id_swdn_toa > 0 ) then used = send_data ( id_swdn_toa, swin, Time_diag, is, js ) endif !------- upward lw flux surface ------- if ( id_lwup_sfc > 0 ) then used = send_data ( id_lwup_sfc, b_surf, Time_diag, is, js ) endif !------- downward lw flux surface ------- if ( id_lwdn_sfc > 0 ) then used = send_data ( id_lwdn_sfc, surf_lw_down, Time_diag, is, js ) endif !------- temperature tendency due to radiation ------------ if ( id_tdt_rad > 0 ) then used = send_data ( id_tdt_rad, tdt_rad, Time_diag, is, js, 1 ) endif !------- total radiative flux (at half levels) ----------- if ( id_flux_rad > 0 ) then used = send_data ( id_flux_rad, flux_rad, Time_diag, is, js, 1 ) endif !------- longwave radiative flux (at half levels) -------- if ( id_flux_lw > 0 ) then used = send_data ( id_flux_lw, net, Time_diag, is, js, 1 ) endif if ( id_flux_sw > 0 ) then used = send_data ( id_flux_sw, flux_sw, Time_diag, is, js, 1 ) endif if ( id_entrop_rad > 0 ) then do k=1,n entrop_rad(:,:,k) =tdt_rad(:,:,k)/t(:,:,k)*phalfgrey(1,1,n+1)/1.e5 end do used = send_data ( id_entrop_rad, entrop_rad, Time_diag, is, js, 1 ) endif return end subroutine grey_radiation ! ================================================================================== subroutine grey_radiation_end() module_is_initialized = .false. end subroutine grey_radiation_end ! ================================================================================== end module grey_radiation_mod