!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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 spectral_physics_mod use fms_mod, only: mpp_pe, mpp_root_pe, error_mesg, FATAL, write_version_number, set_domain, file_exist, & field_size, read_data, write_data, fms_init use constants_mod, only: grav, pi use time_manager_mod, only: time_type, set_time, get_time, operator(-), operator(/=), time_manager_init use press_and_geopot_mod, only: pressure_variables, compute_pressures_and_heights use transforms_mod, only: get_grid_boundaries, get_deg_lon, get_deg_lat, get_wts_lat, & get_grid_domain, get_lon_max, get_lat_max, grid_domain, area_weighted_global_mean use global_integral_mod, only: mass_weighted_global_integral use spectral_dynamics_mod, only: get_reference_sea_level_press, get_num_levels use physics_driver_mod, only: physics_driver_init, physics_driver_down, physics_driver_up, physics_driver_end, & surf_diff_type, do_moist_in_phys_up, get_diff_t, get_radturbten, zero_radturbten use moist_processes_mod, only: moist_processes_init, moist_processes, moist_processes_end use mcm_moist_processes_mod, only: mcm_moist_processes, mcm_moist_processes_init, mcm_moist_processes_end use tracer_type_mod, only: tracer_type use field_manager_mod, only: MODEL_ATMOS use tracer_manager_mod, only: get_number_tracers, get_tracer_index, NO_TRACER implicit none private public :: spectral_physics_init, spectral_physics_down, spectral_physics_up, & spectral_physics_end, surf_diff_type, spectral_physics_moist character(len=128), parameter :: version = & '$Id: spectral_physics.F90,v 13.0 2006/03/28 21:17:25 fms Exp $' character(len=128), parameter :: tagname = & '$Name: mom4p1_pubrel_dec2009_nnz $' integer, parameter :: num_time_levels=2 real, allocatable, dimension(:,: ) :: rad_lon_2d, rad_lat_2d, area_2d real, allocatable, dimension(:,:,: ) :: diff_cu_mo logical, allocatable, dimension(:,: ) :: convect real, allocatable, dimension(:,:,:,: ) :: diag_tracers integer :: num_levels, num_tracers, nhum integer :: is, ie, js, je logical :: module_is_initialized = .false. logical :: do_mcm_moist_processes contains !------------------------------------------------------------------------------------------------ subroutine spectral_physics_init(Time, axes, Surf_diff, nhum_in, p_half, do_mcm_moist_processes_in) type(time_type), intent(in) :: Time integer, intent(in), dimension(:) :: axes type(surf_diff_type), intent(inout) :: Surf_diff integer, intent(in) :: nhum_in real, intent(in), dimension(:,:,:) :: p_half logical, intent(in) :: do_mcm_moist_processes_in real, allocatable, dimension(:,:,:,:) :: grid_tracers real, allocatable, dimension(:) :: rad_lon, rad_lat, wts_lat, lon_boundaries, lat_boundaries real, dimension(2) :: radiation_ref_press_surf = (/ 101325., 81060. /) real, allocatable, dimension(:,:) :: radiation_ref_press, rconvect real, allocatable, dimension(:) :: p_half_1d, ln_p_half_1d, shalf real, allocatable, dimension(:) :: p_full_1d, ln_p_full_1d, sfull real :: reference_sea_level_press integer :: i, j, num_diag, lon_max, lat_max character(len=4) :: ch1, ch2, ch3, ch4, ch5, ch6 character(len=64) :: file integer, dimension(4) :: siz logical :: do_donner if(module_is_initialized) return call write_version_number(version, tagname) call fms_init call time_manager_init nhum = nhum_in do_mcm_moist_processes = do_mcm_moist_processes_in call get_grid_domain(is, ie, js, je) allocate(rad_lon(is:ie), rad_lat(js:je), wts_lat(js:je)) allocate(lon_boundaries(ie-is+2), lat_boundaries(je-js+2)) call get_num_levels(num_levels) allocate(radiation_ref_press(num_levels+1,2)) allocate(p_half_1d(num_levels+1), ln_p_half_1d(num_levels+1), shalf(num_levels+1)) allocate(p_full_1d(num_levels ), ln_p_full_1d(num_levels ), sfull(num_levels )) allocate (rad_lon_2d(is:ie,js:je)) allocate (rad_lat_2d(is:ie,js:je)) allocate ( area_2d(is:ie,js:je)) call pressure_variables(p_half_1d,ln_p_half_1d,radiation_ref_press(1:num_levels,1),ln_p_full_1d,radiation_ref_press_surf(1)) call pressure_variables(p_half_1d,ln_p_half_1d,radiation_ref_press(1:num_levels,2),ln_p_full_1d,radiation_ref_press_surf(2)) radiation_ref_press(num_levels+1,:) = radiation_ref_press_surf call get_reference_sea_level_press(reference_sea_level_press) call pressure_variables(p_half_1d, ln_p_half_1d, p_full_1d, ln_p_full_1d, reference_sea_level_press) shalf = p_half_1d/reference_sea_level_press sfull = p_full_1d/reference_sea_level_press call get_lon_max(lon_max) call get_lat_max(lat_max) call get_deg_lon(rad_lon) call get_deg_lat(rad_lat) call get_wts_lat(wts_lat) rad_lon = pi*rad_lon/180. rad_lat = pi*rad_lat/180. do j=js,je rad_lat_2d(:,j) = rad_lat(j) area_2d(:,j) = wts_lat(j)/(2.*lon_max) enddo do i=is,ie rad_lon_2d(i,:) = rad_lon(i) enddo call get_grid_boundaries(lon_boundaries, lat_boundaries) call get_number_tracers(MODEL_ATMOS, num_diag=num_diag) if(num_diag > 0) then call error_mesg('spectral_physics_init', & 'This version of the spectral atmospheric model not coded to handle diagnostic tracers',FATAL) endif allocate (diag_tracers(is:ie,js:je,num_levels,num_diag)) diag_tracers = 0. call get_number_tracers(MODEL_ATMOS, num_prog=num_tracers) call set_domain(grid_domain) if(do_mcm_moist_processes) then call mcm_moist_processes_init( ie-is+1, je-js+1, num_levels, axes, Time) else call moist_processes_init(ie-is+1, je-js+1, num_levels, lon_boundaries, lat_boundaries, & radiation_ref_press(:,1), axes, Time, do_donner) endif allocate(grid_tracers(is:ie, js:je, num_levels, num_tracers)) grid_tracers = 0. call physics_driver_init(Time, lon_boundaries, lat_boundaries, axes, radiation_ref_press, grid_tracers, Surf_diff, p_half) if(sum(grid_tracers) /= 0.) then call error_mesg('spectral_physics_init','This version of the spectral atmospheric model not coded to handle'// & ' initialization of tracer fields by physics_driver_init',FATAL) endif deallocate(rad_lon, rad_lat, wts_lat, lon_boundaries, lat_boundaries, grid_tracers) deallocate(radiation_ref_press, p_half_1d, ln_p_half_1d, shalf, p_full_1d, ln_p_full_1d, sfull) allocate (diff_cu_mo(is:ie,js:je,num_levels)) allocate ( convect(is:ie,js:je)) file = 'INPUT/spectral_physics.res.nc' if(file_exist(trim(file))) then call field_size(trim(file), 'diff_cu_mo', siz) if(siz(1) /= lon_max .or. siz(2) /= lat_max .or. siz(3) /= num_levels) then write(ch1,'(i4)') siz(1) write(ch2,'(i4)') siz(2) write(ch3,'(i4)') siz(3) write(ch4,'(i4)') lon_max write(ch5,'(i4)') lat_max write(ch6,'(i4)') num_levels call error_mesg('spectral_physics_init','Resolution of restart data is incorrect.'// & ' Restart data: lon_max='//ch1//', lat_max='//ch2//', num_levels='//ch3// & ' Should be: lon_max='//ch4//', lat_max='//ch5//', num_levels='//ch6, FATAL) endif call read_data(trim(file), 'diff_cu_mo', diff_cu_mo, grid_domain) allocate (rconvect(is:ie,js:je)) call read_data(trim(file), 'convect', rconvect, grid_domain) ! No interface for reading/writing netcdf logicals where(rconvect == 1.) convect = .true. elsewhere convect = .false. endwhere deallocate (rconvect) else diff_cu_mo = 0. convect = .false. endif module_is_initialized = .true. return end subroutine spectral_physics_init !------------------------------------------------------------------------------------------------ subroutine spectral_physics_down(Time_prev, Time, Time_next, previous, current, & p_half, p_full, z_half, z_full, psg, ug, vg, tg, grid_tracers, & frac_land, rough_mom, albedo, t_surf, u_star, b_star, q_star, dtau_du, dtau_dv, tau_x, tau_y, & albedo_vis_dir, albedo_nir_dir, albedo_vis_dif, albedo_nir_dif, & dt_ug, dt_vg, dt_tg, dt_tracers, flux_sw, flux_sw_dir, flux_sw_dif, flux_sw_down_vis_dir, & flux_sw_down_vis_dif, flux_sw_down_total_dir, flux_sw_down_total_dif, flux_sw_vis, & flux_sw_vis_dir, flux_sw_vis_dif, flux_lw, coszen, gust, Surf_diff) type(time_type), intent(in) :: Time_prev, Time, Time_next integer, intent(in) :: previous, current real, intent(in), dimension(:,:,: ) :: p_full, z_full real, intent(in), dimension(:,:,: ) :: p_half, z_half real, intent(in), dimension(:,:,: ) :: psg real, intent(in), dimension(:,:,:,: ) :: ug, vg, tg real, intent(inout), dimension(:,:,:,:,:) :: grid_tracers real, intent(in), dimension(:,: ) :: frac_land, rough_mom, albedo, t_surf, u_star, b_star, q_star, dtau_du, dtau_dv real, intent(inout), dimension(:,: ) :: tau_x, tau_y real, intent(in), dimension(:,: ) :: albedo_vis_dir, albedo_nir_dir, albedo_vis_dif, albedo_nir_dif real, intent(inout), dimension(:,:,: ) :: dt_ug, dt_vg, dt_tg real, intent(inout), dimension(:,:,:,: ) :: dt_tracers real, intent(out), dimension(:,: ) :: flux_sw, flux_sw_dir, flux_sw_dif, flux_sw_down_vis_dir real, intent(out), dimension(:,: ) :: flux_sw_down_vis_dif, flux_sw_down_total_dir, flux_sw_down_total_dif real, intent(out), dimension(:,: ) :: flux_sw_vis, flux_sw_vis_dir, flux_sw_vis_dif, flux_lw, coszen, gust type(surf_diff_type), intent(inout) :: Surf_diff real, dimension(size(dt_tracers,4)) :: gavg_rrv integer :: nco2 !************************************************************************************** if(.not.module_is_initialized) then call error_mesg('spectral_physics_down','spectral_physics module is not initialized', FATAL) end if nco2 = get_tracer_index(MODEL_ATMOS, 'co2') gavg_rrv = 0.0 if(nco2 /= NO_TRACER) then gavg_rrv(nco2) = mass_weighted_global_integral(grid_tracers(:,:,:,current,nco2), psg(:,:,current))/ & (area_weighted_global_mean(psg(:,:,current))/grav) endif if(do_moist_in_phys_up()) then call physics_driver_down(1, ie-is+1, 1, je-js+1, Time_prev, Time, Time_next, & rad_lat_2d, rad_lon_2d, area_2d, & p_half, p_full, z_half, z_full, & ug(:,:,:,current), vg(:,:,:,current), & tg(:,:,:,current), grid_tracers(:,:,:,current,nhum), & grid_tracers(:,:,:,current,:), ug(:,:,:,previous), vg(:,:,:,previous), & tg(:,:,:,previous), grid_tracers(:,:,:,previous,nhum), & grid_tracers(:,:,:,previous,:), & frac_land, rough_mom, albedo, & albedo_vis_dir, albedo_nir_dir, albedo_vis_dif, albedo_nir_dif, & t_surf, u_star, b_star, q_star, & dtau_du, dtau_dv, tau_x, tau_y, & dt_ug, dt_vg, dt_tg, & dt_tracers(:,:,:,nhum), dt_tracers, flux_sw(:,:), & flux_sw_dir, flux_sw_dif, flux_sw_down_vis_dir, & flux_sw_down_vis_dif, flux_sw_down_total_dir, & flux_sw_down_total_dif, flux_sw_vis, & flux_sw_vis_dir, flux_sw_vis_dif, & flux_lw, coszen, gust, Surf_diff, gavg_rrv) else call physics_driver_down(1, ie-is+1, 1, je-js+1, Time_prev, Time, Time_next, & rad_lat_2d, rad_lon_2d, area_2d, & p_half, p_full, z_half, z_full, & ug(:,:,:,current), vg(:,:,:,current), & tg(:,:,:,current), grid_tracers(:,:,:,current,nhum), & grid_tracers(:,:,:,current,:), ug(:,:,:,previous), vg(:,:,:,previous), & tg(:,:,:,previous), grid_tracers(:,:,:,previous,nhum), & grid_tracers(:,:,:,previous,:), & frac_land, rough_mom, albedo, & albedo_vis_dir, albedo_nir_dir, albedo_vis_dif, albedo_nir_dif, & t_surf, u_star, b_star, q_star, & dtau_du, dtau_dv, tau_x, tau_y, & dt_ug, dt_vg, dt_tg, & dt_tracers(:,:,:,nhum), dt_tracers, flux_sw, & flux_sw_dir, flux_sw_dif, flux_sw_down_vis_dir, & flux_sw_down_vis_dif, flux_sw_down_total_dir, & flux_sw_down_total_dif, flux_sw_vis, & flux_sw_vis_dir, flux_sw_vis_dif, & flux_lw, coszen, gust, Surf_diff, gavg_rrv, diff_cum_mom=diff_cu_mo, & moist_convect=convect) endif return end subroutine spectral_physics_down !------------------------------------------------------------------------------------------------ subroutine spectral_physics_up(Time_prev, Time, Time_next, previous, current, p_half, p_full, & z_half, z_full, wg_full, ug, vg, tg, grid_tracers, & frac_land, dt_ug, dt_vg, dt_tg, dt_tracers, Surf_diff, lprec, fprec, gust) type(time_type), intent(in) :: Time_prev, Time, Time_next integer, intent(in) :: previous, current real, intent(in), dimension(:,:,: ) :: p_full, z_full, wg_full real, intent(in), dimension(:,:,: ) :: p_half, z_half real, intent(in), dimension(:,:,:,: ) :: ug, vg, tg real, intent(in), dimension(:,:,:,:,:) :: grid_tracers real, intent(in), dimension(:,: ) :: frac_land real, intent(inout), dimension(:,:,: ) :: dt_ug, dt_vg, dt_tg real, intent(inout), dimension(:,:,:,: ) :: dt_tracers type(surf_diff_type), intent(inout) :: Surf_diff real, intent(out), dimension(:,:) :: lprec, fprec real, intent(inout), dimension(:,:) :: gust if(.not.module_is_initialized) then call error_mesg('spectral_physics_up','spectral_physics module is not initialized', FATAL) end if call physics_driver_up(1, ie-is+1, 1, je-js+1, Time_prev, Time, Time_next, & rad_lat_2d, rad_lon_2d, area_2d, & p_half, p_full, z_half, z_full, & wg_full, ug(:,:,:,current), vg(:,:,:,current), & tg(:,:,:,current), grid_tracers(:,:,:,current,nhum), & grid_tracers(:,:,:,current,:), & ug(:,:,:,previous), vg(:,:,:,previous), & tg(:,:,:,previous), grid_tracers(:,:,:,previous,nhum), & grid_tracers(:,:,:,previous,:), & frac_land, dt_ug, dt_vg, dt_tg, & dt_tracers(:,:,:,nhum), dt_tracers, Surf_diff, & lprec, fprec, gust) return end subroutine spectral_physics_up !------------------------------------------------------------------------------------------------ subroutine spectral_physics_moist(Time_next, delta_t, frac_land, p_half, p_full, z_half, z_full, wg_full, & tg, qg, ug, vg, tracers, lprec, fprec, gust) type(time_type), intent(in) :: Time_next real, intent(in) :: delta_t real, intent(in), dimension(:,:) :: frac_land real, intent(in), dimension(:,:,:) :: p_half, p_full, z_half, z_full, wg_full real, intent(inout), dimension(:,:,:) :: tg, qg, ug, vg real, intent(inout), dimension(:,:,:,:) :: tracers real, intent(inout), dimension(:,:) :: lprec, fprec, gust real, dimension(size(tg,1), size(tg,2), size(tg,3)) :: dt_tg, dt_qg, dt_ug, dt_vg real, dimension(size(tracers,1), size(tracers,2), size(tracers,3), size(tracers,4)) :: dt_tracers real, dimension(size(gust,1), size(gust,2)) :: gust_cv dt_tg=0.; dt_qg=0.; dt_ug=0.; dt_vg=0.; dt_tracers=0. if(do_moist_in_phys_up()) then if ( do_mcm_moist_processes ) then call error_mesg('spectral_physics_moist','do_mcm_moist_processes cannot be .true. when moist_processes'// & ' is called by physics_driver_up', FATAL) endif else if ( do_mcm_moist_processes ) then call mcm_moist_processes(1, ie-is+1, 1, je-js+1, Time_next, delta_t, p_half, p_full, tg, tracers(:,:,:,nhum), lprec, fprec) else call moist_processes(1, ie-is+1, 1, je-js+1, Time_next, delta_t, frac_land, p_half, p_full, z_half, z_full, wg_full, & get_diff_t(), get_radturbten(), & tg, qg, tracers, ug, vg, tg, qg, tracers, ug, vg, dt_tg, dt_qg, dt_tracers, dt_ug, dt_vg, & diff_cu_mo, convect, lprec, fprec, gust_cv, area_2d, rad_lat_2d) gust = sqrt( gust*gust + gust_cv*gust_cv) call zero_radturbten() tg = tg + dt_tg*delta_t qg = qg + dt_qg*delta_t ug = ug + dt_ug*delta_t vg = vg + dt_vg*delta_t tracers = tracers + dt_tracers*delta_t endif endif return end subroutine spectral_physics_moist !------------------------------------------------------------------------------------------------ subroutine spectral_physics_end(Time) character(len=64) :: file type(time_type), intent(in) :: Time real, dimension(is:ie,js:je) :: rconvect if(.not.module_is_initialized) return file = 'RESTART/spectral_physics.res.nc' call write_data(trim(file), 'diff_cu_mo', diff_cu_mo, grid_domain) where(convect) rconvect = 1. elsewhere rconvect = 0. endwhere call write_data(trim(file), 'convect', rconvect, grid_domain) ! pjp: No interface for reading/writing netcdf logicals deallocate(rad_lon_2d, rad_lat_2d, area_2d, diff_cu_mo, convect) call physics_driver_end(Time) if(do_mcm_moist_processes) then call mcm_moist_processes_end else call moist_processes_end endif module_is_initialized = .false. return end subroutine spectral_physics_end !------------------------------------------------------------------------------------------------ end module spectral_physics_mod