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 mpp_mod, only: input_nml_file use fms_mod, only: open_namelist_file, check_nml_error, & mpp_pe, mpp_root_pe, close_file, & write_version_number, stdlog, file_exist 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 19.0 2012/01/06 20:09:59 fms Exp $' character(len=128), parameter :: tagname = '$Name: tikal $' 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 #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=grey_radiation_nml, iostat=io) ierr = check_nml_error(io,'grey_radiation_nml') #else if ( file_exist('input.nml')) then 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) endif #endif call write_version_number(version, tagname) if ( mpp_pe() == mpp_root_pe() ) then logunit = stdlog() write (logunit, nml=grey_radiation_nml) endif 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, 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 integer :: i, j, k, n real :: 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