+ Contact: Isaac Held + Reviewers: Peter Phillipps ++ + + +
+ A spectral transform model for two-dimensional, non-divergent flow on the + surface of the sphere. ++ + + +
+ Integrates the barotropic vorticity equation for nondivergent flow on the + sphere using the spectral transform technique. Also allows for the + inclusion of a passive tracer advected by the same spectral advection + algorithm as the vorticity, and a gridpoint tracer advected with a finite + volume algorithm on the transform grid. The default initial condition + provided as an example is a zonal flow resembling that in the Northern + winter, plus a sinusoidal disurbance localized in midlatitudes. + + For a full description of the model and algorithms used, see barotropic.ps + + The interfaces in this module are the generic intefaces required by the + main program that can be used to drive various idealized atmospheric + models within FMS. Model resolution and related paramters are set in + namelists within the modules barotropic_xxx. ++ + + +
+ fms_mod + constants_mod + transforms_mod + time_manager_mod + diag_manager_mod + barotropic_dynamics_mod + barotropic_physics_mod + barotopic_diagnostics_mod ++ + + +
+ + use atmosphere_mod [,only: atmosphere_init, + atmosphere, + atmosphere_end] + ++ + + +
+ + There are no pubic data types + ++ + + +
+ +subroutine atmosphere_init. Initializes the model. +subroutine atmosphere. Integrates forward one time step +subroutine atmosphere_end. Terminates model, cleaning up memory and finalizing diagnostics. + + ++ + ++ + subroutine atmosphere_init(Time_init, Time, Time_step) + + input: + + type(time_type) :: Time_init -- Initial model time + + type(time_type) :: Time -- Model time + + type(time_type) :: Time_step -- Time step + + When Time=Time_init, the first time step is a forward + step rather than leap frog because a cold start is assumed. + + The FMS main program that runs the solo atmospheric models + obtains Time_init from the diag_table and Time from its namelist. + ++
+ + + subroutine atmosphere(Time) + + input: + + type(time_type) :: Time -- Model time + + Integrates forward one time step + + +
+ + subroutine atmosphere_end(Atmos) + + No calling arguments. + + Terminates model, cleaning up memory and finalizing diagnostics + + +
+
+&atmosphere_nml + + print_interval, integer : time interval in seconds + between prints of global mean energy and enstrophy to standard output ++ + + +
+ + Fatal error message if any public routine is called prior to atmosphere_init + ++ + + +
+ Contact: Isaac Held + Reviewers: Peter Phillipps + ++ + + +
+ + The diagnostics module for the model that solves the barotropic vorticity + equation on the sphere + ++ + + +
+ + Using the diagnostics manager, creates output files for the barotropic model. + Variables currently available for output are + zonal wind + meridional wind + relative vorticity + absolute vorticity + streamfunction + spectral tracer in grid domain + grid tracer + + Whether or not these fields are actually output, the location of the + output, the frequency of output, whether or not the output is averaged + in time or an instantaneous snapshot, is controlled by a + diag_table file utilized by the diagnostics manager module + + One can add other diagnostics by following the (somewhat convoluted) + pattern within the program + + ++ + + +
+ + diag_manaager_mod + transforms_mod + time_manager_mod + barotropic_dynamics_mod + barotropic_physics_mod + ++ + + +
+ + use barotropic_diagnostics_mod [,only: barotropic_diagnostics_init, + barotropic_diagnostics] + ++ + + + +
+ +subroutine barotropic_diagnostics_init +subroutine barotropic_diagnostics + + ++ + subroutine barotropic_diagnostics_init(Time, num_lon, num_lat) + + type(time_type) , intent(in) :: Time + current time + integer, intent(in) :: num_lon, num_lat + num_lon = number of longitudes in global domain + num_lat = number of latitudes in global domain + + + Initializes module + ++ +
+ + + + subroutine barotropic_diagnostics (Time, Grid, Phys, time_index) + + type(time_type), intent(in) :: Time + type(phys_type), intent(in) :: Phys + type(grid_type), intent(in) :: Grid + integer, intent(in) :: time_index + + phys_type is defined in barotropic_physics_mod; + Phys is currently empty as there is no information generated in + barotropic_physics_mod to be output; + + grid_type is defined in barotropic_dynamics_mod: + Grid contains all of the fields to be output + + many of the grid fields in grid_type are dimensioned (lon, lat, time_index) + where time_index = 1 or 2 -- the two time levels needed to update the + state of the model using a leapfrog step are toggled between (:,:,1) + and (:,:,2). The input time_index (which must equal either 1 or 2) + determines which of these two fields is output) + + (this is confusing -- the calling program needs to know what has + been placed in which slot -- it would be better to store this + information within the data type) + + + + +
+ + + diff --git a/src/atmos_spectral_barotropic/barotropic_dynamics.F90 b/src/atmos_spectral_barotropic/barotropic_dynamics.F90 new file mode 100644 index 000000000..ac18b6127 --- /dev/null +++ b/src/atmos_spectral_barotropic/barotropic_dynamics.F90 @@ -0,0 +1,574 @@ +module barotropic_dynamics_mod + +!----------------------------------------------------------------------- +! GNU General Public License +! +! This program 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; either version 2 of +! the License, or (at your option) any later version. +! +! This program 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. +! +! For the full text of the GNU General Public License, +! write to: Free Software Foundation, Inc., +! 675 Mass Ave, Cambridge, MA 02139, USA. +! or see: http://www.gnu.org/licenses/gpl.html +!----------------------------------------------------------------------- + +use fms_mod, only: open_namelist_file, & + open_restart_file, & + file_exist, & + check_nml_error, & + error_mesg, & + mpp_error, & + FATAL, WARNING, & + write_version_number, & + mpp_pe, & + mpp_root_pe, & + read_data, & + write_data, & + set_domain, & + close_file, & + stdlog + +use time_manager_mod, only : time_type, & + get_time, & + operator(==), & + operator(-) + +use constants_mod, only: radius, omega + +use transforms_mod, only: transforms_init, transforms_end, & + get_grid_boundaries, & + trans_spherical_to_grid, trans_grid_to_spherical, & + compute_laplacian, & + get_sin_lat, get_cos_lat, & + get_deg_lon, get_deg_lat, & + get_grid_domain, get_spec_domain, & + spectral_domain, grid_domain, & + vor_div_from_uv_grid, uv_grid_from_vor_div, & + horizontal_advection + +use spectral_damping_mod, only: spectral_damping_init, & + compute_spectral_damping + +use leapfrog_mod, only: leapfrog + +use fv_advection_mod, only: fv_advection_init, & + a_grid_horiz_advection + +use stirring_mod, only: stirring, stirring_end + +!=============================================================================================== +implicit none +private +!=============================================================================================== + +public :: barotropic_dynamics_init, & + barotropic_dynamics, & + barotropic_dynamics_end, & + dynamics_type, & + grid_type, & + spectral_type, & + tendency_type + + +! version information +!=================================================================== +character(len=128) :: version = '$Id: barotropic_dynamics.F90,v 17.0 2009/07/21 03:00:21 fms Exp $' +character(len=128) :: tagname = '$Name: siena_201207 $' +!=================================================================== + +type grid_type + real, pointer, dimension(:,:,:) :: u=>NULL(), v=>NULL(), vor=>NULL(), trs=>NULL(), tr=>NULL() + real, pointer, dimension(:,:) :: pv=>NULL(), stream=>NULL() + real, pointer, dimension(:) :: zonal_u_init=>NULL() +end type +type spectral_type + complex, pointer, dimension(:,:,:) :: vor=>NULL(), trs=>NULL() +end type +type tendency_type + real, pointer, dimension(:,:) :: u=>NULL(), v=>NULL(), trs=>NULL(), tr=>NULL() +end type +type dynamics_type + type(grid_type) :: grid + type(spectral_type) :: spec + type(tendency_type) :: tend + integer :: num_lon, num_lat + logical :: grid_tracer, spec_tracer +end type + +integer, parameter :: num_time_levels = 2 + +integer :: is, ie, js, je, ms, me, ns, ne + +logical :: module_is_initialized = .false. + +real, allocatable, dimension(:) :: sin_lat, cos_lat, rad_lat, rad_lon, & + deg_lat, deg_lon, & + coriolis, glon_bnd, glat_bnd + +integer :: pe, npes + +! namelist parameters with default values + +logical :: check_fourier_imag = .false. +logical :: south_to_north = .true. +logical :: triang_trunc = .true. + +real :: robert_coeff = 0.04 +real :: longitude_origin = 0.0 +real :: raw_filter_coeff = 1.0 + +character(len=64) :: damping_option = 'resolution_dependent' +integer :: damping_order = 4 +real :: damping_coeff = 1.e-04 +real :: damping_coeff_r = 0.0 + +real :: zeta_0 = 8.e-05 +integer :: m_0 = 4 +real :: eddy_width = 15.0 +real :: eddy_lat = 45.0 + +logical :: spec_tracer = .true. +logical :: grid_tracer = .true. + +integer :: num_lat = 128 +integer :: num_lon = 256 +integer :: num_fourier = 85 +integer :: num_spherical = 86 +integer :: fourier_inc = 1 +integer :: cutoff_wn = 30 + +real, dimension(2) :: valid_range_v = (/-1.e3,1.e3/) +character(len=64) :: initial_zonal_wind = 'two_jets' + +namelist /barotropic_dynamics_nml/ check_fourier_imag, south_to_north, & + triang_trunc, & + num_lon, num_lat, num_fourier, & + num_spherical, fourier_inc, & + longitude_origin, damping_option, & + damping_order, damping_coeff, & + damping_coeff_r, robert_coeff, & + spec_tracer, grid_tracer, & + eddy_lat, eddy_width, zeta_0, m_0, & + valid_range_v, initial_zonal_wind, & + cutoff_wn + +contains + +!=============================================================================================== + +subroutine barotropic_dynamics_init (Dyn, Time, Time_init, dt_real) + +type(dynamics_type), intent(inout) :: Dyn +type(time_type) , intent(in) :: Time, Time_init +real, intent(in) :: dt_real + +integer :: i, j + +real, allocatable, dimension(:) :: glon_bnd, glat_bnd +complex, allocatable, dimension(:,:) :: div +real :: xx, yy, dd + +integer :: ierr, io, unit, pe +logical :: root + +! < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > + +call write_version_number (version, tagname) + +pe = mpp_pe() +root = (pe == mpp_root_pe()) + +if (file_exist('input.nml')) then + unit = open_namelist_file () + ierr=1 + do while (ierr /= 0) + read (unit, nml=barotropic_dynamics_nml, iostat=io, end=10) + ierr = check_nml_error (io, 'barotropic_dynamics_nml') + enddo + 10 call close_file (unit) +endif + + +if (root) write (stdlog(), nml=barotropic_dynamics_nml) + +call transforms_init(radius, num_lat, num_lon, num_fourier, fourier_inc, num_spherical, & + south_to_north=south_to_north, & + triang_trunc=triang_trunc, & + longitude_origin=longitude_origin ) + +call get_grid_domain(is,ie,js,je) +call get_spec_domain(ms,me,ns,ne) + +Dyn%num_lon = num_lon +Dyn%num_lat = num_lat +Dyn%spec_tracer = spec_tracer +Dyn%grid_tracer = grid_tracer + +allocate (sin_lat (js:je)) +allocate (cos_lat (js:je)) +allocate (deg_lat (js:je)) +allocate (deg_lon (is:ie)) +allocate (rad_lat (js:je)) +allocate (rad_lon (is:ie)) +allocate (coriolis (js:je)) + +allocate (glon_bnd (num_lon + 1)) +allocate (glat_bnd (num_lat + 1)) + +call get_deg_lon (deg_lon) +call get_deg_lat (deg_lat) +call get_sin_lat (sin_lat) +call get_cos_lat (cos_lat) +call get_grid_boundaries (glon_bnd, glat_bnd, global=.true.) + +coriolis = 2*omega*sin_lat + +rad_lat = deg_lat*atan(1.0)/45.0 +rad_lon = deg_lon*atan(1.0)/45.0 + +call spectral_damping_init(damping_coeff, damping_order, damping_option, cutoff_wn, num_fourier, num_spherical, 1, 0., 0., 0., & + damping_coeff_r=damping_coeff_r) + +allocate (Dyn%spec%vor (ms:me, ns:ne, num_time_levels)) +allocate (Dyn%grid%u (is:ie, js:je, num_time_levels)) +allocate (Dyn%grid%v (is:ie, js:je, num_time_levels)) +allocate (Dyn%grid%vor (is:ie, js:je, num_time_levels)) + +allocate (Dyn%tend%u (is:ie, js:je)) +allocate (Dyn%tend%v (is:ie, js:je)) +allocate (Dyn%grid%stream (is:ie, js:je)) +allocate (Dyn%grid%pv (is:ie, js:je)) +allocate (Dyn%grid%zonal_u_init(js:je)) + +allocate (div (ms:me, ns:ne)) + +call fv_advection_init(num_lon, num_lat, glat_bnd, 360./float(fourier_inc)) +if(Dyn%grid_tracer) then + allocate(Dyn%Grid%tr (is:ie, js:je, num_time_levels)) + allocate(Dyn%Tend%tr (is:ie, js:je)) +endif + +if(Dyn%spec_tracer) then + allocate(Dyn%Grid%trs (is:ie, js:je, num_time_levels)) + allocate(Dyn%Tend%trs (is:ie, js:je)) + allocate(Dyn%Spec%trs (ms:me, ns:ne, num_time_levels)) +endif + +if(trim(initial_zonal_wind) == 'zero') then + Dyn%grid%zonal_u_init = 0.0 +else if(trim(initial_zonal_wind) == 'two_jets') then + do j = js, je + Dyn%grid%zonal_u_init(j) = 25.0*cos_lat(j) & + - 30.0*(cos_lat(j)**3) & + + 300.0*(sin_lat(j)**2)*(cos_lat(j)**6) + enddo +else + call error_mesg('barotropic_dynamics_init',trim(initial_zonal_wind)// & + ' is not a valid value of initial_zonal_wind ', FATAL) +endif + +if(Time == Time_init) then + + do j = js, je + Dyn%Grid%u(:,j,1) = Dyn%grid%zonal_u_init(j) + Dyn%Grid%v(:,j,1) = 0.0 + end do + + call vor_div_from_uv_grid(Dyn%Grid%u (:,:,1), Dyn%Grid%v (:,:,1), & + Dyn%Spec%vor(:,:,1), div) + + call trans_spherical_to_grid(Dyn%Spec%vor(:,:,1), Dyn%Grid%vor(:,:,1)) + + do j = js, je + do i = is, ie + yy = (deg_lat(j)- eddy_lat)/eddy_width + Dyn%Grid%vor(i,j,1) = Dyn%Grid%vor(i,j,1) + & + 0.5*zeta_0*cos_lat(j)*exp(-yy*yy)*cos(m_0*rad_lon(i)) + end do + end do + + call trans_grid_to_spherical(Dyn%Grid%vor(:,:,1), Dyn%Spec%vor(:,:,1)) + + div = (0.,0.) + call uv_grid_from_vor_div (Dyn%Spec%vor(:,:,1), div, & + Dyn%Grid%u (:,:,1), Dyn%Grid%v (:,:,1)) + + if(Dyn%grid_tracer) then + Dyn%Grid%tr = 0.0 + do j = js, je + if(deg_lat(j) > 10.0 .and. deg_lat(j) < 20.0) Dyn%Grid%tr(:,j,1) = 1.0 + if(deg_lat(j) > 70.0 ) Dyn%Grid%tr(:,j,1) = -1.0 + end do + endif + + if(Dyn%spec_tracer) then + Dyn%Grid%trs = 0.0 + do j = js, je + if(deg_lat(j) > 10.0 .and. deg_lat(j) < 20.0) Dyn%Grid%trs(:,j,1) = 1.0 + if(deg_lat(j) > 70.0 ) Dyn%Grid%trs(:,j,1) = -1.0 + end do + call trans_grid_to_spherical(Dyn%Grid%trs(:,:,1), Dyn%Spec%trs(:,:,1)) + endif + +else + + call read_restart(Dyn) + +endif + +module_is_initialized = .true. + +return +end subroutine barotropic_dynamics_init + +!=============================================================================================== + +subroutine barotropic_dynamics(Time, Time_init, Dyn, previous, current, future, delta_t) + +type(time_type) , intent(in) :: Time, Time_init +type(dynamics_type), intent(inout) :: Dyn +integer, intent(in ) :: previous, current, future +real, intent(in ) :: delta_t + +! < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > + +complex, dimension(ms:me, ns:ne) :: dt_vors, dt_divs, stream, zeros, spec_diss +real, dimension(is:ie, js:je) :: dt_vorg +integer :: j, seconds, days + +! < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > + +if(.not.module_is_initialized) then + call error_mesg('barotropic_dynamics','dynamics has not been initialized ', FATAL) +endif + +zeros = (0.,0.) + +do j = js, je + Dyn%grid%pv(:,j) = Dyn%grid%vor(:,j,current) + coriolis(j) +end do + +Dyn%Tend%u = Dyn%Tend%u + Dyn%grid%pv*Dyn%Grid%v(:,:,current) +Dyn%Tend%v = Dyn%Tend%v - Dyn%grid%pv*Dyn%Grid%u(:,:,current) + +call vor_div_from_uv_grid (Dyn%Tend%u, Dyn%Tend%v, dt_vors, dt_divs) + +call compute_spectral_damping(Dyn%Spec%vor(:,:,previous), dt_vors, delta_t) + +call stirring(Time, dt_vors) + +call leapfrog(Dyn%Spec%vor , dt_vors , previous, current, future, delta_t, robert_coeff, raw_filter_coeff) + +call trans_spherical_to_grid(Dyn%Spec%vor(:,:,future), Dyn%Grid%vor(:,:,future)) + +call uv_grid_from_vor_div (Dyn%Spec%vor (:,:,future), zeros, & + Dyn%Grid%u (:,:,future), Dyn%Grid%v (:,:,future)) + +if(minval(Dyn%Grid%v) < valid_range_v(1) .or. maxval(Dyn%Grid%v) > valid_range_v(2)) then + call get_time (Time, seconds, days) + call mpp_error(FATAL,'barotropic_dynamics: Meridional wind out of valid range. Model time=',days,' days ',seconds,' seconds') +endif + +if(Dyn%spec_tracer) call update_spec_tracer(Dyn%Spec%trs, Dyn%Grid%trs, Dyn%Tend%trs, & + Dyn%Grid%u, Dyn%Grid%v, previous, current, future, delta_t) + +if(Dyn%grid_tracer) call update_grid_tracer(Dyn%Grid%tr, Dyn%Tend%tr, & + Dyn%Grid%u, Dyn%Grid%v, previous, current, future, delta_t) + +stream = compute_laplacian(Dyn%Spec%vor(:,:,current), -1) +call trans_spherical_to_grid(stream, Dyn%grid%stream) + +return +end subroutine barotropic_dynamics + +!=================================================================================== + +subroutine update_spec_tracer(tr_spec, tr_grid, dt_tr, ug, vg, & + previous, current, future, delta_t) + +complex, intent(inout), dimension(ms:me, ns:ne, num_time_levels) :: tr_spec +real , intent(inout), dimension(is:ie, js:je, num_time_levels) :: tr_grid +real , intent(inout), dimension(is:ie, js:je ) :: dt_tr +real , intent(in ), dimension(is:ie, js:je, num_time_levels) :: ug, vg +real , intent(in ) :: delta_t +integer, intent(in ) :: previous, current, future + +complex, dimension(ms:me, ns:ne) :: dt_trs + +call horizontal_advection (tr_spec(:,:,current), ug(:,:,current), vg(:,:,current), dt_tr) +call trans_grid_to_spherical (dt_tr, dt_trs) +call compute_spectral_damping (tr_spec(:,:,previous), dt_trs, delta_t) +call leapfrog (tr_spec, dt_trs, previous, current, future, delta_t, robert_coeff, raw_filter_coeff) +call trans_spherical_to_grid (tr_spec(:,:,future), tr_grid(:,:,future)) + +return +end subroutine update_spec_tracer +!========================================================================== + +subroutine update_grid_tracer(tr_grid, dt_tr_grid, ug, vg, & + previous, current, future, delta_t) + +real , intent(inout), dimension(is:ie, js:je, num_time_levels) :: tr_grid +real , intent(inout), dimension(is:ie, js:je ) :: dt_tr_grid +real , intent(in ), dimension(is:ie, js:je, num_time_levels) :: ug, vg + +real , intent(in ) :: delta_t +integer, intent(in ) :: previous, current, future + +real, dimension(size(tr_grid,1),size(tr_grid,2)) :: tr_current, tr_future + +tr_future = tr_grid(:,:,previous) + delta_t*dt_tr_grid +dt_tr_grid = 0.0 +call a_grid_horiz_advection (ug(:,:,current), vg(:,:,current), tr_future, delta_t, dt_tr_grid) +tr_future = tr_future + delta_t*dt_tr_grid +tr_current = tr_grid(:,:,current) + & + robert_coeff*(tr_grid(:,:,previous) + tr_future - 2.0*tr_grid(:,:,current)) +tr_grid(:,:,current) = tr_current +tr_grid(:,:,future) = tr_future + +return +end subroutine update_grid_tracer + +!========================================================================== + +subroutine read_restart(Dyn) + +type(dynamics_type), intent(inout) :: Dyn + +integer :: unit, m, n, nt +real, dimension(ms:me, ns:ne) :: real_part, imag_part + +if(file_exist('INPUT/barotropic_dynamics.res.nc')) then + do nt = 1, 2 + call read_data('INPUT/barotropic_dynamics.res.nc', 'vors_real', real_part, spectral_domain, timelevel=nt) + call read_data('INPUT/barotropic_dynamics.res.nc', 'vors_imag', imag_part, spectral_domain, timelevel=nt) + do n=ns,ne + do m=ms,me + Dyn%Spec%vor(m,n,nt) = cmplx(real_part(m,n),imag_part(m,n)) + end do + end do + if(Dyn%spec_tracer) then + call read_data('INPUT/barotropic_dynamics.res.nc', 'trs_real', real_part, spectral_domain, timelevel=nt) + call read_data('INPUT/barotropic_dynamics.res.nc', 'trs_imag', imag_part, spectral_domain, timelevel=nt) + do n=ns,ne + do m=ms,me + Dyn%Spec%trs(m,n,nt) = cmplx(real_part(m,n),imag_part(m,n)) + end do + end do + endif + call read_data('INPUT/barotropic_dynamics.res.nc', 'u', Dyn%Grid%u (:,:,nt), grid_domain, timelevel=nt) + call read_data('INPUT/barotropic_dynamics.res.nc', 'v', Dyn%Grid%v (:,:,nt), grid_domain, timelevel=nt) + call read_data('INPUT/barotropic_dynamics.res.nc', 'vor', Dyn%Grid%vor(:,:,nt), grid_domain, timelevel=nt) + if(Dyn%spec_tracer) then + call read_data('INPUT/barotropic_dynamics.res.nc', 'trs', Dyn%Grid%trs(:,:,nt), grid_domain, timelevel=nt) + endif + if(Dyn%grid_tracer) then + call read_data('INPUT/barotropic_dynamics.res.nc', 'tr', Dyn%Grid%tr(:,:,nt), grid_domain, timelevel=nt) + endif + end do +else if(file_exist('INPUT/barotropic_dynamics.res')) then + unit = open_restart_file(file='INPUT/barotropic_dynamics.res',action='read') + + do nt = 1, 2 + call set_domain(spectral_domain) + call read_data(unit,Dyn%Spec%vor(:,:, nt)) + if(Dyn%spec_tracer) call read_data(unit,Dyn%Spec%trs(:,:, nt)) + + call set_domain(grid_domain) + call read_data(unit,Dyn%Grid%u (:,:, nt)) + call read_data(unit,Dyn%Grid%v (:,:, nt)) + call read_data(unit,Dyn%Grid%vor (:,:, nt)) + if(Dyn%spec_tracer) call read_data(unit,Dyn%Grid%trs(:,:, nt)) + if(Dyn%grid_tracer) call read_data(unit,Dyn%Grid%tr (:,:, nt)) + + end do + call close_file(unit) + +else + call error_mesg('read_restart', 'restart does not exist', FATAL) +endif + +return +end subroutine read_restart + +!==================================================================== + +subroutine write_restart(Dyn, previous, current) + +type(dynamics_type), intent(in) :: Dyn +integer, intent(in) :: previous, current + +integer :: unit, nt, nn + +do nt = 1, 2 + if(nt == 1) nn = previous + if(nt == 2) nn = current + call write_data('RESTART/barotropic_dynamics.res.nc', 'vors_real', real(Dyn%Spec%vor(:,:,nn)), spectral_domain) + call write_data('RESTART/barotropic_dynamics.res.nc', 'vors_imag', aimag(Dyn%Spec%vor(:,:,nn)), spectral_domain) + if(Dyn%spec_tracer) then + call write_data('RESTART/barotropic_dynamics.res.nc', 'trs_real', real(Dyn%Spec%trs(:,:,nn)), spectral_domain) + call write_data('RESTART/barotropic_dynamics.res.nc', 'trs_imag', aimag(Dyn%Spec%trs(:,:,nn)), spectral_domain) + endif + call write_data('RESTART/barotropic_dynamics.res.nc', 'u', Dyn%Grid%u (:,:,nn), grid_domain) + call write_data('RESTART/barotropic_dynamics.res.nc', 'v', Dyn%Grid%v (:,:,nn), grid_domain) + call write_data('RESTART/barotropic_dynamics.res.nc', 'vor', Dyn%Grid%vor(:,:,nn), grid_domain) + if(Dyn%spec_tracer) then + call write_data('RESTART/barotropic_dynamics.res.nc', 'trs', Dyn%Grid%trs(:,:,nn), grid_domain) + endif + if(Dyn%grid_tracer) then + call write_data('RESTART/barotropic_dynamics.res.nc', 'tr', Dyn%Grid%tr(:,:,nn), grid_domain) + endif +enddo + +!unit = open_restart_file(file='RESTART/barotropic_dynamics.res', action='write') + +!do nt = 1, 2 +! if(nt == 1) nn = previous +! if(nt == 2) nn = current + +! call set_domain(spectral_domain) +! call write_data(unit,Dyn%Spec%vor(:,:, nn)) +! if(Dyn%spec_tracer) call write_data(unit,Dyn%Spec%trs(:,:, nn)) + +! call set_domain(grid_domain) +! call write_data(unit,Dyn%Grid%u (:,:, nn)) +! call write_data(unit,Dyn%Grid%v (:,:, nn)) +! call write_data(unit,Dyn%Grid%vor (:,:, nn)) +! if(Dyn%spec_tracer) call write_data(unit,Dyn%Grid%trs(:,:, nn)) +! if(Dyn%grid_tracer) call write_data(unit,Dyn%Grid%tr (:,:, nn)) +!end do + +!call close_file(unit) + +end subroutine write_restart + +!==================================================================== + +subroutine barotropic_dynamics_end (Dyn, previous, current) + +type(dynamics_type), intent(inout) :: Dyn +integer, intent(in) :: previous, current + +if(.not.module_is_initialized) then + call error_mesg('barotropic_dynamics','dynamics has not been initialized ', FATAL) +endif + +call transforms_end() +call stirring_end() + +call write_restart (Dyn, previous, current) + +module_is_initialized = .false. + +return +end subroutine barotropic_dynamics_end +!=================================================================================== + +end module barotropic_dynamics_mod diff --git a/src/atmos_spectral_barotropic/barotropic_dynamics.html b/src/atmos_spectral_barotropic/barotropic_dynamics.html new file mode 100644 index 000000000..62c97898a --- /dev/null +++ b/src/atmos_spectral_barotropic/barotropic_dynamics.html @@ -0,0 +1,371 @@ + +module barotropic_dynamics_mod + + +
+ + +module barotropic_dynamics_mod
+ + ++ Contact: Isaac Held + Reviewers: Peter Phillipps + ++ + + +
+OVERVIEW
+ ++ + The dynamical core of the spectral transform model for + two-dimensional, non-divergent flow on the surface of the sphere. + ++ + + +
+DESCRIPTION
+ ++ + Integrates the barotropic vorticity equation for nondivergent flow on the + sphere using the spectral transform technique. Also allows for the + inclusion of a passive tracer advected by the same spectral advection + algorithm as the vorticity, and a gridpoint tracer advected with a finite + volume algorithm on the transform grid. The default initial condition + provided as an example is a zonal flow resembling that in the Northern + winter, plus a sinusoidal disurbance localized in midlatitudes. + + For a full description of the model and algorithms used, see + barotropic.ps + + For higher level routines for running this barotropic spectral model, + see atmosphere_mod + + ++ + + +
+OTHER MODULES USED
+ ++ + fms_mod + constants_mod + time_manager_mod + transforms_mod + spectral_damping_mod + leapfrog_mod + fv_advection_mod + ++ + + +
+PUBLIC INTERFACE
+ ++ + use barotropic_dynamics_mod [,only: barotropic_dynamics_init, + barotropic_dynamics, + barotropic_dynamics_end, + dynamics_type, + grid_type, + spectral_type, + tendency_type] + ++ + + +
+PUBLIC DATA
+ ++ ++ + + +
+ + +type grid_type + real, pointer, dimension(:,:,:) :: u, v, vor, trs, tr, pv + real, pointer, dimension(:,:) :: stream +end type + + allocated space for grid fields + + (:,:,:) => (lon, lat, time_level) + (:,:) => (lon, lat) + (lon, lat) on local computational domain + time_level stores the two time levels needed for the + leapfrog step + + u -- eastward velocity (m/s) + v -- northward velocity (m/s) + vor -- vorticity (1/s) + trs -- tracer advected spectrally + tr -- tracer advected on grid + pv -- absolute vorticity, f + vor, where f = 2*omega*sin(lat) (1/s) + stream -- streamfunction (m^2/s) at current time + + + +
+ + +type spectral_type + complex, pointer, dimension(:,:,:) :: vor, trs +end type + + allocated space for spectral fields + + (:,:,:) => (zonal, meridional, time_level) + + vor -- spectral vorticity + trs -- spectral tracer + +
+ + +type tendency_type + real, pointer, dimension(:,:) :: u, v, trs, tr +end type + + allocated space for accumulating tendencies, d/dt, in grid space, + for prognostic variables + + (:,:,:) => (lon, lat) + +
+ + +type dynamics_type + type(grid_type) :: grid + type(spectral_type) :: spec + type(tendency_type) :: tend + integer :: num_lon, num_lat ! size of global domain + logical :: grid_tracer, spec_tracer +end type + + grid_tracer = .true. => tracer with gridpoint advection is beign integrated + similarly for spec_tracer + +
+ +
+PUBLIC ROUTINES
+ ++ +subroutine barotropic_dynamics_init +subroutine barotropic _dynamics +subroutine barotropic_dynamics_end +type (grid_type) +type (spectral_type) +type (tendency_type) +type (dynamics_type) + + + ++ + + ++ + subroutine barotropic_dynamics_init(Dyn, Time, Time_init) + + type(dynamics_type), intent(inout) :: Dyn + type containing all dynamical fields and related information + (see type (dynamics_type)) + + type(time_type) , intent(in) :: Time, Time_init + current time and time at which integeration began + time_type defined by time_manager_mod + + + Initializes the module; + Reads restart from 'INPUT/barotropic_dynamics.res' if Time = Time_init; + otherwise uses default initial conditions + ++ +
+ + + + subroutine barotropic_dynamics & + (Time, Time_init, Dyn, previous, current, future, delta_t) + + type(time_type) , intent(inout) :: Time, Time_init + type(dynamics_type), intent(inout) :: Dyn + integer , intent(in ) :: previous, current, future + real , intent(in ) :: delta_t + + previous, current and future = 1 or 2 + these integers refer to the third dimension of the + three-dimensional fields in Dyn + the fields at time t - delta_t are assumed to be in (:,:,previous) + the fields at time t are assumed to be in (:,:,current) + the fields at time t + delta_t are placed in (:,:,future) + overwriting whatever is already there + + delta_t = time step in seconds + + updates dynamical fields by one time step + + + +
+ + + subroutine barotropic_dynamics_end(Dyn, previous, current) + + type(dynamics_type), intent(inout) :: Dyn + integer, intent(in) :: previous, current + + + Terminates module; + writes restart file to 'RESTART/barotropic_dynamics.res' + + + + +
+
+NAMELIST
+ ++ +&barotropic_dynamics_nml + + integer :: num_lat = 128 + number of latitudes in global grid + + integer :: num_lon = 256 + number of longitudes in global grid + should equal 2*num_lat for Triangular truncation + + integer :: num_fourier = 85 + the retained fourier wavenumber are n*fourier_inc, where + n ranges from 0 to num_fourier + + integer :: num_spherical = 86 + the maximum number of meridional modes for any zonal wavenumber + for triangular truncation, set num_spherical = num_fourier +1 + + integer :: fourier_inc = 1 + creates a "sector" model if fourier_inc > 1; integration domain is + (360 degrees longitude)/fourier_inc + + (the default values listed above define a standard T85 model) + + logical :: check_fourier_imag = .false. + if true, checks to see if fields to be transformed to grid + domain have zero imaginary part to their zonally symmetric + modes; useful for debugging + + logical :: south_to_north = .true. + true => grid runs from south to north + false => grid runs from north to south + + logical :: triangular_trunc = .true. + true => shape of truncation is triangular + false => shape of truncation is rhomboidal + + real :: robert_coeff = 0.04 + x(current) => (1-2r)*x(current) + r*(x(future)+x(previous)) + where r = robert_coeff (non-dimensional) + + real :: longitude_origin = 0.0 + longitude of first longitude, in degrees + (if you want the westgern boundary of first grid boc to be at + 0.0, set longitude_origin = 0.5*360./float(num_lon)) + + integer :: damping_option = 'resolution_dependent' + integer :: damping_order = 4 + real :: damping_coeff = 1.e-04 + + damping = nu*(del^2)^n where n = damping order + damping_option = 'resolution_dependent' or 'resolution_independent' + = 'resolution_dependent' => nu is set so that the damping rate for the + mode (m=0,n=num_spherical-1) equals damping_coeff (in 1/s) + For triangular truncation, damping_coeff is then the + rate of damping of the highest retained mode + + = 'resolution_independent' => nu = damping_coef + + + real :: zeta_0 = 8.e-05 + integer :: m_0 = 4 + real :: eddy_width = 15.0 + real :: eddy_lat = 45.0 + + eddy component of the initial condition is sinusoidal with + wavenumber m_0 and with a gaussian distribution of + vorticity in latitude, centered at eddy_lat with half-width + eddy_width + + zeta_0 ( 1/s) + eddy_width and eddy_lat (degrees) + + logical :: spec_tracer = .true. + logical :: grid_tracer = .true. + spec_tracer = true => a passive tracer is carried that is advected + spectrally, with the same algorithm as the vorticity + grid_tracer = ture => a passive tracer is carried that is advected + on the spectral transform grid by a finite-volume algorithm + (see barotropic.ps ) + Both tracers can be carried simultaeneously + +The vorticity and the tracers are initialized within subroutine + barotropic_dynamics_init + + real, dimension(2) :: valid_range_v = -1000., 1000. + A valid range for meridional wind. Model terminates if meridional wind + goes outside the valid range. Allows model to terminate gracefully when, + for example, the model becomes numerically unstable. + ++ + + +
+ERROR MESSAGES
+ ++ + "Dynamics has not been initialized" + -- barotropic_dynamics_init must be called before any other + routines in the module are called + + "restart does not exist" + -- Time is not equal to Time_init at initalization, but the file + 'INPUT/barotropic_dynamics.res' does not exit + + ++ + + +
+ + diff --git a/src/atmos_spectral_barotropic/barotropic_physics.F90 b/src/atmos_spectral_barotropic/barotropic_physics.F90 new file mode 100644 index 000000000..8590bacce --- /dev/null +++ b/src/atmos_spectral_barotropic/barotropic_physics.F90 @@ -0,0 +1,173 @@ +module barotropic_physics_mod + +!----------------------------------------------------------------------- +! GNU General Public License +! +! This program 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; either version 2 of +! the License, or (at your option) any later version. +! +! This program 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. +! +! For the full text of the GNU General Public License, +! write to: Free Software Foundation, Inc., +! 675 Mass Ave, Cambridge, MA 02139, USA. +! or see: http://www.gnu.org/licenses/gpl.html +!----------------------------------------------------------------------- + +use fms_mod, only: open_namelist_file, & + open_restart_file, & + file_exist, & + check_nml_error, & + error_mesg, & + FATAL, WARNING, & + write_version_number, & + mpp_pe, & + mpp_root_pe, & + fms_init, fms_end, & + read_data, & + write_data, & + set_domain, & + close_file, & + stdlog + +use transforms_mod, only: get_sin_lat, get_cos_lat, & + get_deg_lon, get_deg_lat, & + get_wts_lat, & + get_grid_domain, get_spec_domain, & + grid_domain + +use time_manager_mod, only: time_type + +!======================================================================== +implicit none +private +!======================================================================== + +public :: barotropic_physics_init, & + barotropic_physics, & + barotropic_physics_end, & + phys_type + +! version information +!======================================================================== +character(len=128) :: version = '$Id: barotropic_physics.F90,v 10.0 2003/10/24 22:00:58 fms Exp $' +character(len=128) :: tagname = '$Name: siena_201207 $' +!======================================================================== + +type phys_type + real, pointer, dimension(:,:) :: empty=>NULL() +end type + +logical :: module_is_initialized = .false. + +integer :: is, ie, js, je + + +integer :: pe +logical :: root + +real, allocatable, dimension(:) :: rad_lat, & + deg_lat, & + sin_lat, & + cos_lat, & + wts_lat + +! namelist +!======================================================================== + +logical :: empty + +namelist /barotropic_physics_nml/ empty +!======================================================================== + +contains + +!======================================================================== + +subroutine barotropic_physics_init(Phys) + +type(phys_type), intent(inout) :: Phys + +integer :: j, unit, ierr, io + +call write_version_number (version, tagname) + +pe = mpp_pe() +root = (pe == mpp_root_pe()) + +! read the namelist + +if (file_exist('input.nml')) then + unit = open_namelist_file () + ierr=1 + do while (ierr /= 0) + read (unit, nml=barotropic_physics_nml, iostat=io, end=10) + ierr = check_nml_error (io, 'barotropic_physics_nml') + enddo + 10 call close_file (unit) +endif + +call get_grid_domain(is,ie,js,je) + +allocate ( rad_lat (js:je) ) +allocate ( deg_lat (js:je) ) +allocate ( sin_lat (js:je) ) +allocate ( cos_lat (js:je) ) +allocate ( wts_lat (js:je) ) + +call get_wts_lat(wts_lat) +call get_deg_lat(deg_lat) +rad_lat = deg_lat*atan(1.)/45. +sin_lat = sin(rad_lat) +cos_lat = cos(rad_lat) + +module_is_initialized = .true. + +return +end subroutine barotropic_physics_init + +!======================================================================= + +subroutine barotropic_physics(Time, dt_ug, dt_vg, ug, vg, & + delta_t, previous, current, Phys) + +real, intent(inout), dimension(is:ie, js:je) :: dt_ug, dt_vg +real, intent(in) , dimension(is:ie, js:je, 2) :: ug, vg + +real , intent(in) :: delta_t +integer, intent(in) :: previous, current + +type(time_type), intent(in) :: Time +type(phys_type), intent(inout) :: Phys + +if(.not.module_is_initialized) call error_mesg('barotropic_physics', & + 'barotropic_physics is not initialized', FATAL) + +! dt_ug = dt_ug +f(ug,vg) +! dt_vg = dt_vg +f(ug,vg) +! Phys%empty = + +return +end subroutine barotropic_physics + +!====================================================================== + +subroutine barotropic_physics_end(Phys) + +type(phys_type), intent(in) :: Phys + +if(.not.module_is_initialized) call error_mesg('barotropic_physics_end', & + 'barotropic_physics is not initialized', FATAL) + +module_is_initialized = .false. +return +end subroutine barotropic_physics_end + +!====================================================================== + +end module barotropic_physics_mod diff --git a/src/atmos_spectral_barotropic/barotropic_physics.html b/src/atmos_spectral_barotropic/barotropic_physics.html new file mode 100644 index 000000000..e391ca50f --- /dev/null +++ b/src/atmos_spectral_barotropic/barotropic_physics.html @@ -0,0 +1,165 @@ + +module barotropic_physics_mod + + +
+ + +module barotropic_physics_mod
+ + ++ Contact: Isaac Held + Reviewers: Peter Phillipps + ++ + + +
+OVERVIEW
+ ++ + A module that allows one to add processes that act in the grid domain + to the dynamics of the barotropic model on the sphere + ++ + + +
+DESCRIPTION
+ ++ + A module that allows one to add processes that act in the grid domain + to the dynamics of the barotropic model on the sphere. Currently, + does nothing! + ++ + + +
+OTHER MODULES USED
+ ++ + fms_mod + transforms_mod + time_manager_mod + ++ + + +
+PUBLIC INTERFACE
+ ++ + use barotropic_physics_mod [,only: barotropic_physics_init, + barotropic_physics, + barotropic_physics_end, + phys_type] + ++ + + +
+PUBLIC DATA
+ ++ ++ + +
+ + +type phys_type + real, pointer, dimension(:,:) :: empty +end type + + fields from physics module made available for diagnostics + +
+ + +
+PUBLIC ROUTINES
+ ++ +subroutine barotropic_physics_init +subroutine barotropic_physics +subroutine barotropic_physics_end +type(phys_type) + + + ++ + subroutine barotropic_physics_init(Phys) + + type(phys_type) , intent(inout) :: Phys + + + Initializes module + ++ +
+ + + + subroutine barotropic_physics (Time, dt_ug, dt_vg, ug, vg, & + delta_t, previous, current, Phys) + + real, intent(inout), dimension(:,:) :: dt_ug, dt_vg + + the u and v tendencies onto which tendencies due to + the grid-point physics are added (m/(s^2)) + + real, intent(in) , dimension(:,:, 2) :: ug, vg + the grid zonal and meridional velocities (m/s) + the third index is the time-index used in the leapfrog step + + real , intent(in) :: delta_t + time step (s) + + integer, intent(in) :: previous, current + = 1 or 2 + ug(:,:,previous) is the velocity at t-delta_t + ug(:,:,current ) is the velocity at t + + type(time_type), intent(in) :: Time + type(phys_type), intent(inout) :: Phys + + + + +
+ + + + subroutine barotropic_physics_end (Phys) + + type(phys_type), intent(inout) :: Phys + + + + +
+ + + diff --git a/src/atmos_spectral_barotropic/stirring.F90 b/src/atmos_spectral_barotropic/stirring.F90 new file mode 100644 index 000000000..7de10c3f1 --- /dev/null +++ b/src/atmos_spectral_barotropic/stirring.F90 @@ -0,0 +1,252 @@ +module stirring_mod + +!----------------------------------------------------------------------- +! GNU General Public License +! +! This program 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; either version 2 of +! the License, or (at your option) any later version. +! +! This program 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. +! +! For the full text of the GNU General Public License, +! write to: Free Software Foundation, Inc., +! 675 Mass Ave, Cambridge, MA 02139, USA. +! or see: http://www.gnu.org/licenses/gpl.html +!----------------------------------------------------------------------- + +! Stirring is computed as described in the following paper: + +! Vallis, Gerber, Kushner, Cash, 2003: A Mechanism and Simple Dynamical Model of the North Atlantic Oscillation and Annular Modes. +! J. Atmos. Sci., 61, 264-280. + +! Stirring is not part of barotropic_physics because barotropic_physics appears to be intended for +! operations that are done completely in grid space. Stirring is computed partly in spectral space. + +use constants_mod, only: pi + +use time_manager_mod, only: time_type + +use fms_mod, only: open_namelist_file, check_nml_error, close_file, write_version_number, & + stdlog, mpp_pe, mpp_root_pe, file_exist, read_data, write_data, error_mesg, FATAL + +use transforms_mod, only: get_spec_domain, get_grid_domain, trans_spherical_to_grid, trans_grid_to_spherical, & + grid_domain, get_lon_max, get_lat_max, get_deg_lon, get_deg_lat, get_grid_boundaries, & + get_num_fourier, get_num_spherical, spectral_domain + +use diag_manager_mod, only: diag_axis_init, register_static_field, register_diag_field, send_data + +implicit none +private + +integer :: ms,me,ns,ne,is,ie,js,je +integer :: id_str_amp, id_g_stir_sqr, id_stir +logical :: used +logical, allocatable, dimension(:,:) :: wave_mask ! wave_mask(m,n) = .true. if spherical wave (m,n) is to be excited +complex, allocatable, dimension(:,:) :: s_stir ! stirring. Saved from one time step to the next +real, allocatable, dimension(:,:) :: localize ! localizes the stirring +real, allocatable, dimension(:,:) :: g_stir_sqr ! time mean of g_stir**2 over entire integration +integer, allocatable, dimension(:) :: seed ! random number seed +real :: astir, bstir +integer :: num_steps, num_fourier, num_spherical, nseed + +logical :: module_is_initialized = .false. + +character(len=128) :: version = '$Id: stirring.F90,v 17.0 2009/07/21 03:00:25 fms Exp $' +character(len=128) :: tagname = '$Name: siena_201207 $' + +public :: stirring_init, stirring, stirring_end + +real :: decay_time=2*86400, amplitude=0.0, lat0=45., widthy=12. +logical :: do_localize=.true.!Default true to allow forcing to be localized in physical space. Set to false to have forcing everywhere. + +! Set B to a non-zero value for stirring that has zonal structure. +! The strength of the stirring at latitude=lat0 is: amplitude*(1.0 + B*exp(-.5*((lon-lon0)/widthx)**2)) +real :: lon0=180., B=0.0, widthx=45., C=1.0 ! widthx +integer :: n_total_forcing_max = 15 !total wavenumbers LESS THAN this number will be forced +integer :: n_total_forcing_min = 9 !total wavenumbers GREATER THAN this number will be forced +integer :: zonal_forcing_min = 3 !Zonal wavenumbers GREATER THAN this number will be forced, subject to total wavenumber constraints + +namelist / stirring_nml / decay_time, amplitude, lat0, lon0, widthy, widthx, B, do_localize, n_total_forcing_max, n_total_forcing_min, zonal_forcing_min + +contains + +!================================================================================================================================ +subroutine stirring_init(dt, Time, id_lon, id_lat, id_lonb, id_latb) +real, intent(in) :: dt +type(time_type), intent(in) :: Time +integer, intent(in) :: id_lon, id_lat, id_lonb, id_latb +real :: xx, kk, rad_to_deg +integer :: i,j,m,n,ierr,io,unit,lon_max,lat_max +real, allocatable, dimension(:) :: ampx, ampy, lon, lat, lonb, latb +real, allocatable, dimension(:,:) :: real_part, imag_part + +if(module_is_initialized) return + +call write_version_number (version, tagname) + +if (file_exist('input.nml')) then + unit = open_namelist_file () + ierr=1 + do while (ierr /= 0) + read (unit, nml=stirring_nml, iostat=io, end=10) + ierr = check_nml_error (io, 'stirring_nml') + enddo + 10 call close_file (unit) +endif +if(mpp_pe() == mpp_root_pe()) write(stdlog(), nml=stirring_nml) + +call get_lon_max(lon_max) +call get_lat_max(lat_max) + +allocate(lon (lon_max )) ; lon = 0.0 +allocate(lat (lat_max )) ; lat = 0.0 +allocate(lonb(lon_max+1)) ; lonb = 0.0 +allocate(latb(lat_max+1)) ; latb = 0.0 + +call get_deg_lon(lon) +call get_deg_lat(lat) + +module_is_initialized = .true. +if(amplitude == 0.0) return ! stirring does nothing more unless amplitude is non-zero + +call get_spec_domain(ms,me,ns,ne) +call get_grid_domain(is,ie,js,je) +call get_num_fourier(num_fourier) +call get_num_spherical(num_spherical) + +allocate(wave_mask(ms:me,ns:ne)); wave_mask = .false. +allocate(s_stir(ms:me,ns:ne)); s_stir = cmplx(0.0,0.0) +allocate(ampx(is:ie)); ampx = 0.0 +allocate(ampy(js:je)); ampy = 0.0 +allocate(localize(is:ie,js:je)); localize = 0.0 +allocate(g_stir_sqr(is:ie,js:je)); g_stir_sqr = 0.0 + +! wave_mask is .true. when (m+n > 9) .and. (m+n < 15) .and. (m > 3) +do m=(zonal_forcing_min+1),(n_total_forcing_max-1) + if(m >= ms .and. m <= me) then + do n=(n_total_forcing_min+1)-m,(n_total_forcing_max-1)-m + if(n >= ns .and. n <= ne) then + wave_mask(m,n) = .true. + endif + enddo + endif +enddo + +astir = sqrt(1.0 - exp(-2*dt/decay_time)) +bstir = exp(-dt/decay_time) + +do i=is,ie + xx = lon(i)-lon0 + ! Make sure xx falls in the range -180. to +180. + kk = nint(xx/360.) + xx = xx - 360.*kk + ampx(i) = (1 + B*exp(-.5*(xx/widthx)**2)) +enddo +do j=js,je + ampy(j) = exp(-.5*((lat(j)-lat0)/widthy)**2) +enddo +if (do_localize) then + do j=js,je + do i=is,ie + localize(i,j) = ampx(i)*ampy(j) + enddo + enddo +else + localize = 1.0 +endif + +deallocate(ampx, ampy) + +num_steps = 0 +id_g_stir_sqr = register_static_field('stirring_mod', 'stirring_sqr', (/id_lon,id_lat/), 'stirring sqrared', '1/sec^4') +id_str_amp = register_static_field('stirring_mod', 'stirring_amp', (/id_lon,id_lat/), 'amplitude of stirring', 'none') +id_stir = register_diag_field ('stirring_mod', 'stirring', (/id_lon,id_lat/), Time, 'stirring', '1/sec^2') +used = send_data(id_str_amp, amplitude*localize) + +call random_seed(size=nseed) +allocate(seed(nseed)) + +if(file_exist('INPUT/stirring.res.nc')) then + allocate(real_part(ms:me,ns:ne), imag_part(ms:me,ns:ne)) + call read_data('INPUT/stirring.res.nc', 'stir_real', real_part, spectral_domain) + call read_data('INPUT/stirring.res.nc', 'stir_imag', imag_part, spectral_domain) + do n=ns,ne + do m=ms,me + s_stir(m,n) = cmplx(real_part(m,n),imag_part(m,n)) + end do + end do + deallocate(real_part, imag_part) + call read_data('INPUT/stirring.res.nc', 'ran_nmbr_seed', seed, no_domain=.true.) + call random_seed(put=seed) +endif + +end subroutine stirring_init +!================================================================================================================================ +subroutine stirring(Time, dt_vors) +type(time_type), intent(in) :: Time +complex, dimension(ms:me,ns:ne), intent(inout) :: dt_vors +real, dimension(is:ie,js:je) :: g_stir +complex, dimension(ms:me,ns:ne) :: new_stirring +real, dimension(0:num_fourier,0:num_spherical,2) :: ran_nmbrs +integer :: i,j,m,n +real :: x,y + +if(.not.module_is_initialized) then + call error_mesg('stirring', 'stirring_init has not been called', FATAL) +end if + +if(amplitude == 0.0) return ! stirring does nothing unless amplitude is non-zero + +call random_number(ran_nmbrs) + +do n=ns,ne +do m=ms,me + if(wave_mask(m,n)) then + new_stirring(m,n) = amplitude*astir*cmplx(2*ran_nmbrs(m,n,1)-1, 2*ran_nmbrs(m,n,2)-1) + else + new_stirring(m,n) = cmplx(0.0,0.0) + endif +enddo +enddo +call trans_spherical_to_grid(new_stirring,g_stir) +g_stir = localize*g_stir +call trans_grid_to_spherical(g_stir,new_stirring) +if(ms == 0 .and. ns == 0) then + new_stirring(0,0)=cmplx(0.0,0.0) ! A non-zero global mean is introduced by the grid space computation, but we don't want it. +endif +s_stir = bstir*s_stir + new_stirring !This is equation A.6 in Vallis et al 2004 - DOI:10.1175/1520-0469(2004)061<0264:AMASDM>2.0.CO;2 + +dt_vors = dt_vors + s_stir +call trans_spherical_to_grid(s_stir,g_stir) +g_stir_sqr = g_stir_sqr + g_stir*g_stir +num_steps = num_steps + 1 +used = send_data(id_stir, g_stir, Time) + +end subroutine stirring +!================================================================================================================================ +subroutine stirring_end + +if(.not.module_is_initialized) return + +if(amplitude == 0.0) return ! stirring does nothing unless amplitude is non-zero + +g_stir_sqr = g_stir_sqr/num_steps +used = send_data(id_g_stir_sqr, g_stir_sqr) + +call write_data('RESTART/stirring.res.nc', 'stir_real', real(s_stir), spectral_domain) +call write_data('RESTART/stirring.res.nc', 'stir_imag', aimag(s_stir), spectral_domain) +call random_seed(get=seed) +call write_data('RESTART/stirring.res.nc', 'ran_nmbr_seed', seed, no_domain=.true.) + +deallocate(wave_mask, s_stir, localize, g_stir_sqr) +module_is_initialized = .false. + +end subroutine stirring_end +!================================================================================================================================ + +end module stirring_mod diff --git a/src/atmos_spectral_shallow/atmosphere.F90 b/src/atmos_spectral_shallow/atmosphere.F90 new file mode 100644 index 000000000..7336f9b72 --- /dev/null +++ b/src/atmos_spectral_shallow/atmosphere.F90 @@ -0,0 +1,281 @@ +Module atmosphere_mod + +!----------------------------------------------------------------------- +! GNU General Public License +! +! This program 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; either version 2 of +! the License, or (at your option) any later version. +! +! This program 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. +! +! For the full text of the GNU General Public License, +! write to: Free Software Foundation, Inc., +! 675 Mass Ave, Cambridge, MA 02139, USA. +! or see: http://www.gnu.org/licenses/gpl.html +!----------------------------------------------------------------------- + +!========================================================================= + +use fms_mod, only: open_namelist_file, & + open_restart_file, & + file_exist, & + check_nml_error, & + error_mesg, & + FATAL, WARNING, & + write_version_number, & + mpp_pe, & + mpp_root_pe, & + close_file, & + stdlog + +use mpp_mod, only: mpp_max + +use constants_mod, only: radius_earth => radius, & + omega_earth => omega + +use transforms_mod, only : get_deg_lon, & + get_deg_lat, & + get_grid_boundaries, & + get_grid_domain, & + get_spec_domain, & + area_weighted_global_mean, & + atmosphere_domain + +use time_manager_mod, only : time_type, & + set_time, & + get_time, & + interval_alarm, & + operator(+), & + operator(<), & + operator(==) + +use shallow_dynamics_mod, only : shallow_dynamics_init, & + shallow_dynamics, & + shallow_dynamics_end, & + dynamics_type + +use shallow_physics_mod, only : shallow_physics_init, & + shallow_physics, & + shallow_physics_end, & + phys_type + +use shallow_diagnostics_mod, only : shallow_diagnostics_init, & + shallow_diagnostics + +use stirring_mod, only: stirring_init + +!======================================================================== +implicit none +private +!======================================================================== + +! version information +!======================================================================== +character(len=128) :: version = '$Id: atmosphere.F90,v 14.0 2007/03/15 22:13:18 fms Exp $' +character(len=128) :: tagname = '$Name: siena_201207 $' +!======================================================================== + +public :: atmosphere_init, & + atmosphere, & + atmosphere_end, & + atmosphere_domain + +!======================================================================== + +integer, parameter :: num_time_levels = 2 + +integer :: unit, seconds, days +integer :: pe, npes +integer :: previous, current, future +logical :: root + +integer :: dt_integer +real :: dt_real +type(time_type) :: dt_time_type, Time_init, Time_step + +real :: delta_t ! = 2*dt_real for leapfrog step + +integer, dimension(2) :: axes + +type(phys_type), save :: Phys +type(dynamics_type), save :: Dyn + +integer :: is, ie, js, je, ms, me, ns, ne +integer :: num_lon, num_lat +integer, dimension(4) :: axis_id ! axes identifiers + +logical :: module_is_initialized =.false. + + +integer :: print_interval +! namelist +!======================================================================== +namelist /atmosphere_nml/ print_interval +!======================================================================== + +contains +!======================================================================= + +subroutine atmosphere_init(Time_init_in, Time, Time_step_in) + +type (time_type), intent(in) :: Time_init_in, Time, Time_step_in + +integer :: i, j, n, nn, ierr, io, unit, id_lon, id_lat, id_lonb, id_latb +integer :: nlon, nlat + +pe = mpp_pe() +root = (pe == mpp_root_pe()) + +Time_step = Time_step_in +call get_time(Time_step, seconds, days) +dt_integer = 86400*days + seconds +dt_real = float(dt_integer) +dt_time_type = Time_step +Time_init = Time_init_in + +! read the namelist + +if (file_exist('input.nml')) then + unit = open_namelist_file () + ierr=1 + do while (ierr /= 0) + read (unit, nml=atmosphere_nml, iostat=io, end=10) + ierr = check_nml_error (io, 'atmosphere_nml') + enddo + 10 call close_file (unit) +endif +call write_version_number(version, tagname) +if (root) write (stdlog(), nml=atmosphere_nml) + +call shallow_dynamics_init (Dyn, Time, Time_init, dt_real) + +call get_grid_domain(is,ie,js,je) +call get_spec_domain(ms,me,ns,ne) + +num_lon = Dyn%num_lon +num_lat = Dyn%num_lat + +nlon = ie+1-is ! size of grid on each processor +nlat = je+1-js + +call shallow_physics_init(Phys) +call shallow_diagnostics_init(Time, num_lon, num_lat, id_lon, id_lat, id_lonb, id_latb) +call stirring_init(dt_real, Time, id_lon, id_lat, id_lonb, id_latb) + +if(Time == Time_init) then + previous = 1 + current = 1 +else + previous = 1 + current = 2 +endif + +module_is_initialized = .true. + +return +end subroutine atmosphere_init + +!===================================================================== + +subroutine atmosphere(Time) + +type (time_type), intent(in) :: Time +integer :: day, second, dt + + +if(.not.module_is_initialized) then + call error_mesg('atmosphere', & + 'atmosphere_init has not been called', FATAL) +end if + +call get_time(Time_step, second, day) +dt = second + 86400*day + +Dyn%Tend%u = 0.0 +Dyn%Tend%v = 0.0 +Dyn%Tend%h = 0.0 +if(Dyn%grid_tracer) Dyn%Tend%tr = 0.0 +if(Dyn%spec_tracer) Dyn%Tend%trs = 0.0 + +if(Time == Time_init) then + delta_t = dt_real + future = 2 +else + delta_t = 2.0*dt_real + future = previous +endif + +call shallow_physics(Time, & + Dyn%Tend%u, Dyn%Tend%v, Dyn%Tend%h, & + Dyn%Grid%u, Dyn%Grid%v, Dyn%Grid%h, & + delta_t, previous, current, & + Phys) + +call shallow_dynamics(Time, Time_init, & + Dyn, previous, current, future, delta_t) + +previous = current +current = future + +call shallow_diagnostics (Time+Time_step, Dyn%Grid, Phys, current) + +call get_time(Time+Time_step, second, day) +if(mod(second+86400*day, print_interval) < dt) & + call global_diag(second, day, current) + +return +end subroutine atmosphere + +!======================================================================================= + +subroutine global_diag(second, day, current) + +integer, intent(in) :: second, day, current + +real :: enstrophy, div_squared, max_Froude +real, dimension(size(Dyn%Grid%u,1), size(Dyn%Grid%u,2)) :: speed + +enstrophy = & +area_weighted_global_mean(Dyn%grid%vor(:,:,current)*Dyn%grid%vor(:,:,current)) + +div_squared = & +area_weighted_global_mean(Dyn%grid%div(:,:,current)*Dyn%grid%div(:,:,current)) + +speed = Dyn%Grid%u(:,:,current)*Dyn%Grid%u(:,:,current) +& + Dyn%Grid%v(:,:,current)*Dyn%Grid%v(:,:,current) +max_Froude = maxval(speed/Dyn%Grid%h(:,:,current)) +call mpp_max(max_Froude) + +if(root) then + write(*,1000) day, second, enstrophy, div_squared, max_Froude +end if +1000 format(1x, 'day =',i6,2x,'second =', i6, & + 2x,'enstrophy = ',e13.6,3x,'div_squared = ',e13.6, 3x, & + 'max_Froude = ', e10.3) + +return +end subroutine global_diag + +!=============================================================================== +subroutine atmosphere_end + +if(.not.module_is_initialized) then + call error_mesg('atmosphere_end', & + 'atmosphere_init has not been called.', FATAL) +end if + +call shallow_physics_end (Phys) +call shallow_dynamics_end (Dyn, previous, current) + +module_is_initialized = .false. + +return +end subroutine atmosphere_end + +!======================================================================================= +end module atmosphere_mod diff --git a/src/atmos_spectral_shallow/atmosphere.html b/src/atmos_spectral_shallow/atmosphere.html new file mode 100644 index 000000000..13e269faf --- /dev/null +++ b/src/atmos_spectral_shallow/atmosphere.html @@ -0,0 +1,164 @@ + +module atmosphere_mod + + +
+ +module atmosphere_mod
+ + ++ Contact: Isaac Held + Reviewers: Peter Phillipps ++ + + +
+OVERVIEW
+ ++ A spectral transform model for the shallow water equations on the sphere ++ + + +
+DESCRIPTION
+ ++ Integrates the shallow water equations for hydrostatic flow in a thin layer + of homogeneous fluid on the sphere, using the spectral transform technique. + Also allows for the inclusion of a passive tracer advected by the + spectral advection algorithm as the vorticity, and a gridpoint tracer + advected with a finite volume algorithm on the transform grid. + The default experiment is forced by a "monsoonal" mass source, starting + from a state of rest. + + For a full description of the model and algorithms used, see shallow.ps + + The interfaces in this module are the generic intefaces required by the + main program that can be used to drive various idealized atmospheric + models within FMS. Model resolution and related parameters are set in + namelists within the modules shallow_xxx. ++ + + +
+OTHER MODULES USED
+ ++ fms_mod + constants_mod + transforms_mod + time_manager_mod + diag_manager_mod + shallow_dynamics_mod + shallow_physics_mod + shallow_diagnostics_mod ++ + + +
+PUBLIC INTERFACE
+ ++ use atmosphere_mod [,only: atmosphere_init, + atmosphere, + atmosphere_end] ++ + + +
+PUBLIC DATA
+ ++ There are no public data types ++ + + +
+PUBLIC ROUTINES
+ ++subroutine atmosphere_init. Initializes the model. +subroutine atmosphere. Integrates forward one time step +subroutine atmosphere_end. Terminates model, cleaning up memory and finalizing diagnostics. + + ++ + ++ subroutine atmosphere_init(Time_init, Time, Time_step) + + input: + + type(time_type) :: Time_init -- Initial model time + + type(time_type) :: Time -- Model time + + type(time_type) :: Time_step -- Time step + + When Time=Time_init, the first time step is a forward + step rather than leap frog because a cold start is assumed. + + The FMS main program that runs the solo atmospheric models + obtains Time_init from the diag_table and Time from its namelist. + ++
+ + + subroutine atmosphere(Time) + + input: + + type(time_type) :: Time -- Model time + + Integrates forward one time step + + +
+ + subroutine atmosphere_end + + No calling arguments. + + Terminates model, cleaning up memory and finalizing diagnostics + +
+NAMELIST
+ ++&atmosphere_nml + + print_interval, integer : time interval in seconds + between prints of global mean energy and enstrophy to standard output ++ + + +
+ERROR MESSAGES
+ ++ Fatal error message if subroutine atmosphere or atmosphere_end + is called prior to atmosphere_init. ++ + +
+ + diff --git a/src/atmos_spectral_shallow/shallow.pdf b/src/atmos_spectral_shallow/shallow.pdf new file mode 100644 index 000000000..3ae5dbd1a Binary files /dev/null and b/src/atmos_spectral_shallow/shallow.pdf differ diff --git a/src/atmos_spectral_shallow/shallow_diagnostics.F90 b/src/atmos_spectral_shallow/shallow_diagnostics.F90 new file mode 100644 index 000000000..3ea5dc0a8 --- /dev/null +++ b/src/atmos_spectral_shallow/shallow_diagnostics.F90 @@ -0,0 +1,259 @@ + +module shallow_diagnostics_mod + +!----------------------------------------------------------------------- +! GNU General Public License +! +! This program 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; either version 2 of +! the License, or (at your option) any later version. +! +! This program 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. +! +! For the full text of the GNU General Public License, +! write to: Free Software Foundation, Inc., +! 675 Mass Ave, Cambridge, MA 02139, USA. +! or see: http://www.gnu.org/licenses/gpl.html +!----------------------------------------------------------------------- + +use fms_mod, only: write_version_number + +use transforms_mod, only: get_grid_boundaries, & + get_deg_lon, & + get_deg_lat, & + get_grid_domain, & + get_spec_domain, & + grid_domain, & + area_weighted_global_mean + +use diag_manager_mod, only: diag_axis_init, & + register_diag_field, & + register_static_field, & + send_data + +use time_manager_mod, only: time_type, & + get_time + + +use shallow_physics_mod, only: phys_type +use shallow_dynamics_mod, only: grid_type + + +implicit none +private + +public :: shallow_diagnostics_init, & + shallow_diagnostics + +character(len=84), parameter :: version = '$Id: shallow_diagnostics.F90,v 10.0 2003/10/24 22:01:02 fms Exp $' +character(len=84), parameter :: tagname = '$Name: siena_201207 $' +character(len=8) :: axiset = 'shallow' +character(len=84) :: mod_name = 'shallow_diagnostics' + +logical :: module_is_initialized = .false. + +integer :: id_vor, id_stream, id_pv, id_u, id_v, id_div, id_h, id_trs, id_tr, id_d_geopot, id_u_sqd, id_v_sqd, id_h_sqd, id_u_sqd_mean, id_v_sqd_mean, id_h_sqd_mean, id_ekin, id_ekin_density, id_eq_geopot, id_e_kin_real_units, id_e_pot_real_units, id_e_tot_real_units, id_u_rms, id_vcomp_vor, id_ucomp_vcomp + +integer :: is, ie, js, je + +contains + +!----------------------------------------------------------------------------------------------------------------- +subroutine shallow_diagnostics_init(Time, lon_max, lat_max, id_lon, id_lat, id_lonb, id_latb) + +type(time_type), intent(in) :: Time +integer, intent(in) :: lon_max, lat_max +integer, intent(out):: id_lon, id_lat, id_lonb, id_latb + +real, dimension(lon_max ) :: lon +real, dimension(lon_max+1) :: lonb +real, dimension(lat_max ) :: lat +real, dimension(lat_max+1) :: latb + +integer, dimension(2) :: axis_2d + +integer :: log_unit +integer :: namelist_unit, ierr, io +real :: rad_to_deg +logical :: used + +call write_version_number(version, tagname) + +call get_grid_domain(is, ie, js, je) + +rad_to_deg = 45./atan(1.) +call get_grid_boundaries(lonb,latb,global=.true.) +call get_deg_lon(lon) +call get_deg_lat(lat) + +id_lonb=diag_axis_init('lonb', rad_to_deg*lonb, 'degrees_E', 'x', 'longitude edges', set_name=axiset, Domain2=grid_domain) +id_latb=diag_axis_init('latb', rad_to_deg*latb, 'degrees_N', 'y', 'latitude edges', set_name=axiset, Domain2=grid_domain) +id_lon =diag_axis_init('lon', lon, 'degrees_E', 'x', 'longitude', set_name=axiset, Domain2=grid_domain, edges=id_lonb) +id_lat =diag_axis_init('lat', lat, 'degrees_N', 'y', 'latitude', set_name=axiset, Domain2=grid_domain, edges=id_latb) + +axis_2d(1) = id_lon +axis_2d(2) = id_lat + +id_u = register_diag_field(mod_name, 'ucomp' , axis_2d, Time, 'u_wind' , 'm/s' ) +id_v = register_diag_field(mod_name, 'vcomp' , axis_2d, Time, 'v_wind' , 'm/s' ) +id_vor = register_diag_field(mod_name, 'vor' , axis_2d, Time, 'relative vorticity' , '1/s' ) +id_div = register_diag_field(mod_name, 'div' , axis_2d, Time, 'divergence' , '1/s' ) +id_h = register_diag_field(mod_name, 'h' , axis_2d, Time, 'geopotential' , 'm2/s2' ) +id_pv = register_diag_field(mod_name, 'pv_corrected' , axis_2d, Time, 'potential vorticity' , 's/m2' ) +id_stream = register_diag_field(mod_name, 'stream', axis_2d, Time, 'streamfunction' , 'm^2/s' ) +id_trs = register_diag_field(mod_name, 'trs' , axis_2d, Time, 'spectral tracer' , 'none' ) +id_tr = register_diag_field(mod_name, 'tr' , axis_2d, Time, 'grid tracer' , 'none' ) +id_d_geopot = register_diag_field(mod_name, 'deep_geopot', axis_2d, Time, 'deep_geopot' , 'm2/s2') + +id_u_sqd = register_diag_field(mod_name, 'ucomp_sqd' , axis_2d, Time, 'u_wind_sqd' , 'm^2/s^2' ) +id_v_sqd = register_diag_field(mod_name, 'vcomp_sqd' , axis_2d, Time, 'v_wind_sqd' , 'm^2/s^2' ) +id_h_sqd = register_diag_field(mod_name, 'h_sqd' , axis_2d, Time, 'geopotential_sqd' , 'm4/s4' ) + +id_u_sqd_mean = register_diag_field(mod_name, 'ucomp_sqd_mean' , Time, 'u_wind_sqd_mean' , 'm^2/s^2' ) +id_v_sqd_mean = register_diag_field(mod_name, 'vcomp_sqd_mean' , Time, 'v_wind_sqd_mean' , 'm^2/s^2' ) +id_h_sqd_mean = register_diag_field(mod_name, 'h_sqd_mean' , Time, 'geopotential_sqd_mean' , 'm4/s4' ) + +id_ekin = register_diag_field(mod_name, 'e_kin' , Time, 'kinetic_energy' , 'm^2/s^2' ) +id_ekin_density = register_diag_field(mod_name, 'e_kin_density' , Time, 'kinetic_energy_density' , 'm^3/s^2' ) + +id_eq_geopot = register_diag_field(mod_name, 'eq_geopot' , Time, 'equilibrium_geopotential' , 'm^2/s^2' ) +id_e_kin_real_units = register_diag_field(mod_name, 'e_kin_real_units' , Time, 'e_kin_real_units' , 'J/kg' ) +id_e_pot_real_units = register_diag_field(mod_name, 'e_pot_real_units' , Time, 'e_pot_real_units' , 'J/kg' ) +id_e_tot_real_units = register_diag_field(mod_name, 'e_tot_real_units' , Time, 'e_tot_real_units' , 'J/kg' ) + +id_u_rms = register_diag_field(mod_name, 'u_rms' , Time, 'r_rms' , 'm/s' ) + +id_vcomp_vor = register_diag_field(mod_name, 'vcomp_vor' , axis_2d, Time, 'vcomp * relative vorticity' , 'm/s^2' ) + +id_ucomp_vcomp = register_diag_field(mod_name, 'ucomp_vcomp' , axis_2d, Time, 'ucomp * vcomp' , 'm^2/s^2' ) + +module_is_initialized = .true. + +return +end subroutine shallow_diagnostics_init + +!-------------------------------------------------------------------------------------------- + +subroutine shallow_diagnostics(Time, Grid, Phys, time_index) + +type(time_type), intent(in) :: Time +type(phys_type), intent(in) :: Phys +type(grid_type), intent(in) :: Grid +integer, intent(in) :: time_index + +real :: e_kin_real_units, e_pot_real_units, e_tot_real_units, eq_geopot + +logical :: used + +if(id_u > 0) used = send_data(id_u , Grid%u (:,:, time_index) , time) +if(id_v > 0) used = send_data(id_v , Grid%v (:,:, time_index) , time) +if(id_vor > 0) used = send_data(id_vor , Grid%vor (:,:, time_index) , time) +if(id_div > 0) used = send_data(id_div , Grid%div (:,:, time_index) , time) +if(id_h > 0) used = send_data(id_h , Grid%h (:,:, time_index) , time) +if(id_pv > 0) used = send_data(id_pv , Grid%pv (:,:) , time) +if(id_stream > 0) used = send_data(id_stream , Grid%stream (:,:) , time) +if(id_tr > 0) used = send_data(id_tr , Grid%tr (:,:, time_index) , time) +if(id_trs > 0) used = send_data(id_trs , Grid%trs (:,:, time_index) , time) +if(id_d_geopot > 0) used = send_data(id_d_geopot, Grid%deep_geopot (:,:) , time) + +if (id_u_sqd > 0) then + used = send_data(id_u_sqd , Grid%u (:,:, time_index)**2 , time) +endif + +if (id_v_sqd > 0) then + used = send_data(id_v_sqd , Grid%v (:,:, time_index)**2 , time) +endif + +if (id_h_sqd > 0) then + used = send_data(id_h_sqd , Grid%h (:,:, time_index)**2 , time) +endif + +if (id_u_sqd_mean > 0) then + used = send_data(id_u_sqd_mean , area_weighted_global_mean(Grid%u (:,:, time_index)**2) , time) +endif + +if (id_v_sqd_mean > 0) then + used = send_data(id_v_sqd_mean , area_weighted_global_mean(Grid%v (:,:, time_index)**2) , time) +endif + +if (id_h_sqd_mean > 0) then + used = send_data(id_h_sqd_mean , area_weighted_global_mean(Grid%h (:,:, time_index)**2) , time) +endif + +if (id_ekin > 0) then + used = send_data(id_ekin , 0.5*(area_weighted_global_mean(Grid%u (:,:, time_index)**2) + area_weighted_global_mean(Grid%v (:,:, time_index)**2)) , time) +endif + +if (id_ekin_density > 0) then + used = send_data(id_ekin_density , 0.5*(area_weighted_global_mean(Grid%h (:,:, time_index)*(Grid%u (:,:, time_index)**2)) + area_weighted_global_mean(Grid%h (:,:, time_index)*(Grid%v (:,:, time_index)**2))) , time) +endif + +eq_geopot = 0. +e_kin_real_units = 0. +e_pot_real_units = 0. + +if (id_eq_geopot > 0) then + + eq_geopot = area_weighted_global_mean(Grid%h (:,:, time_index)) + used = send_data(id_eq_geopot , eq_geopot, time) + +endif + + +if (id_e_kin_real_units > 0) then + + if (eq_geopot == 0.) eq_geopot = area_weighted_global_mean(Grid%h (:,:, time_index)) + + e_kin_real_units = 0.5*(area_weighted_global_mean(Grid%h (:,:, time_index)*(Grid%u (:,:, time_index)**2)) + area_weighted_global_mean(Grid%h (:,:, time_index)*(Grid%v (:,:, time_index)**2))) / eq_geopot + + used = send_data(id_e_kin_real_units , e_kin_real_units, time) + +endif + +if (id_e_pot_real_units > 0) then + + if (eq_geopot == 0.) eq_geopot = area_weighted_global_mean(Grid%h (:,:, time_index)) + + e_pot_real_units = 0.5*(area_weighted_global_mean(Grid%h (:,:, time_index)**2.)) / eq_geopot + + used = send_data(id_e_pot_real_units , e_pot_real_units, time) + +endif + +if (id_e_tot_real_units > 0) then + + if (eq_geopot == 0.) eq_geopot = area_weighted_global_mean(Grid%h (:,:, time_index)) + + if (e_kin_real_units == 0.) then + e_kin_real_units = 0.5*(area_weighted_global_mean(Grid%h (:,:, time_index)*(Grid%u (:,:, time_index)**2)) + area_weighted_global_mean(Grid%h (:,:, time_index)*(Grid%v (:,:, time_index)**2))) / eq_geopot + endif + + if (e_pot_real_units == 0.) then + e_pot_real_units = 0.5*(area_weighted_global_mean(Grid%h (:,:, time_index)**2.)) / eq_geopot + endif + + e_tot_real_units = e_kin_real_units + e_pot_real_units + + used = send_data(id_e_tot_real_units , e_tot_real_units, time) + +endif + +if (id_u_rms > 0) then + + used = send_data(id_u_rms , (area_weighted_global_mean(Grid%u (:,:, time_index)**2) + area_weighted_global_mean(Grid%v (:,:, time_index)**2))**0.5, time) + +endif + + +if(id_vcomp_vor > 0) used = send_data(id_vcomp_vor , Grid%v (:,:, time_index) * Grid%vor (:,:, time_index) , time) +if(id_ucomp_vcomp > 0) used = send_data(id_ucomp_vcomp , Grid%u (:,:, time_index) * Grid%v (:,:, time_index) , time) + +return +end subroutine shallow_diagnostics +!-------------------------------------------------------------------------------------------- + +end module shallow_diagnostics_mod diff --git a/src/atmos_spectral_shallow/shallow_diagnostics.html b/src/atmos_spectral_shallow/shallow_diagnostics.html new file mode 100644 index 000000000..65db99f52 --- /dev/null +++ b/src/atmos_spectral_shallow/shallow_diagnostics.html @@ -0,0 +1,150 @@ + +module shallow_diagnostics_mod + + +
+ + +module shallow_diagnostics_mod
+ + ++ Contact: Isaac Held + Reviewers: Peter Phillipps + ++ + + +
+OVERVIEW
+ ++ + The diagnostics module for the model that solves the shallow water + equation on the sphere + ++ + + +
+DESCRIPTION
+ ++ + Using the diagnostics manager, creates output files for the shallow model. + Variables currently available for output are + zonal wind + meridional wind + relative vorticity + absolute vorticity + streamfunction + spectral tracer in grid domain + grid tracer + + Whether or not these fields are actually output, the location of the + output, the frequency of output, whether or not the output is averaged + in time or an instantaneous snapshot, is controlled by a + diag_table file utilized by the diagnostics manager module + + One can add other diagnostics by following the (somewhat convoluted) + pattern within the program + + ++ + + +
+OTHER MODULES USED
+ ++ + diag_manaager_mod + transforms_mod + time_manager_mod + shallow_dynamics_mod + shallow_physics_mod + ++ + + +
+PUBLIC INTERFACE
+ ++ + use shallow_diagnostics_mod [,only: shallow_diagnostics_init, + shallow_diagnostics] + ++ + + + +
+PUBLIC ROUTINES
+ ++ +subroutine shallow_diagnostics_init +subroutine shallow_diagnostics + + ++ + subroutine shallow_diagnostics_init(Time, num_lon, num_lat) + + type(time_type) , intent(in) :: Time + current time + integer, intent(in) :: num_lon, num_lat + num_lon = number of longitudes in global domain + num_lat = number of latitudes in global domain + + + Initializes module + ++ +
+ + + + subroutine shallow_diagnostics (Time, Grid, Phys, time_index) + + type(time_type), intent(in) :: Time + type(phys_type), intent(in) :: Phys + type(grid_type), intent(in) :: Grid + integer, intent(in) :: time_index + + phys_type is defined in shallow_physics_mod; ; + + grid_type is defined in shallow_dynamics_mod: + Grid contains all of the fields to be output + + many of the grid fields in grid_type are dimensioned (lon, lat, time_index) + where time_index = 1 or 2 -- the two time levels needed to update the + state of the model using a leapfrog step are toggled between (:,:,1) + and (:,:,2). The input time_index (which must equal either 1 or 2) + determines which of these two fields is output) + + (this is confusing -- the calling program needs to know what has + been placed in which slot -- it would be better to store this + information within the data type) + + + + +
+ + + diff --git a/src/atmos_spectral_shallow/shallow_dynamics.F90 b/src/atmos_spectral_shallow/shallow_dynamics.F90 new file mode 100644 index 000000000..7ba1ec3f8 --- /dev/null +++ b/src/atmos_spectral_shallow/shallow_dynamics.F90 @@ -0,0 +1,690 @@ +module shallow_dynamics_mod + +!----------------------------------------------------------------------- +! GNU General Public License +! +! This program 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; either version 2 of +! the License, or (at your option) any later version. +! +! This program 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. +! +! For the full text of the GNU General Public License, +! write to: Free Software Foundation, Inc., +! 675 Mass Ave, Cambridge, MA 02139, USA. +! or see: http://www.gnu.org/licenses/gpl.html +!----------------------------------------------------------------------- + +use fms_mod, only: open_namelist_file, & + open_restart_file, & + file_exist, & + check_nml_error, & + error_mesg, & + FATAL, & + write_version_number, & + mpp_pe, & + mpp_root_pe, & + read_data, & + write_data, & + set_domain, & + close_file, & + stdlog + +use time_manager_mod, only : time_type, & + get_time, & + operator(==), & + operator(-) + +use constants_mod, only : radius, omega, DEG_TO_RAD + +use transforms_mod, only: transforms_init, transforms_end, & + get_grid_boundaries, horizontal_advection, & + trans_spherical_to_grid, trans_grid_to_spherical, & + compute_laplacian, get_eigen_laplacian, & + get_sin_lat, get_cos_lat, & + get_deg_lon, get_deg_lat, & + get_grid_domain, get_spec_domain, & + spectral_domain, grid_domain, & + vor_div_from_uv_grid, uv_grid_from_vor_div, & + area_weighted_global_mean + +use spectral_damping_mod, only: spectral_damping_init, compute_spectral_damping + +use leapfrog_mod, only: leapfrog + +use fv_advection_mod, only : fv_advection_init, a_grid_horiz_advection + +use stirring_mod, only : stirring, stirring_end + +!====================================================================================== +implicit none +private +!====================================================================================== + +public :: shallow_dynamics_init, shallow_dynamics, shallow_dynamics_end, & + dynamics_type, grid_type, spectral_type, tendency_type + + +! version information +!=================================================================== +character(len=128) :: version = '$Id: shallow_dynamics.F90,v 10.0 2003/10/24 22:01:02 fms Exp $' +character(len=128) :: tagname = '$Name: siena_201207 $' +!=================================================================== + +type grid_type + real, pointer, dimension(:,:,:) :: u=>NULL(), v=>NULL(), vor=>NULL(), div=>NULL(), h=>NULL(), trs=>NULL(), tr=>NULL() + real, pointer, dimension(:,:) :: stream=>NULL(), pv=>NULL(), deep_geopot=>NULL() +end type +type spectral_type + complex, pointer, dimension(:,:,:) :: vor=>NULL(), div=>NULL(), h=>NULL(), trs=>NULL() +end type +type tendency_type + real, pointer, dimension(:,:) :: u=>NULL(), v=>NULL(), h=>NULL(), trs=>NULL(), tr=>NULL() +end type +type dynamics_type + type(grid_type) :: grid + type(spectral_type) :: spec + type(tendency_type) :: tend + integer :: num_lon, num_lat + logical :: grid_tracer, spec_tracer +end type + + + +integer, parameter :: num_time_levels = 2 + +integer :: is, ie, js, je, ms, me, ns, ne + +logical :: module_is_initialized = .false. + +real, allocatable, dimension(:) :: sin_lat, cos_lat, rad_lat, deg_lat, deg_lon, & + coriolis + +real, allocatable, dimension(:,:) :: eigen + +integer :: pe, npes + + +! namelist parameters with default values + +integer :: num_lon = 256 +integer :: num_lat = 128 +integer :: num_fourier = 85 +integer :: num_spherical = 86 +integer :: fourier_inc = 1 +integer :: cutoff_wn = 30 +! (these define a standard T85 model) + +logical :: check_fourier_imag = .false. +logical :: south_to_north = .true. +logical :: triang_trunc = .true. + +real :: robert_coeff = 0.04 +real :: robert_coeff_tracer = 0.04 +real :: longitude_origin = 0.0 +real :: raw_filter_coeff = 1.0 + +character(len=64) :: damping_option = 'resolution_dependent' +integer :: damping_order = 4 +real :: damping_coeff = 1.e-04 +real :: h_0 = 3.e04 + +real :: u_deep_mag = 0. +real :: n_merid_deep_flow = 3. +real :: u_upper_mag_init = 0. + +logical :: spec_tracer = .true. +logical :: grid_tracer = .true. + +!Options for injecting an initial vortex pair +real :: lon_centre_init_cyc = 0. +real :: lat_centre_init_cyc = 60. +real :: lon_centre_init_acyc = 180. +real :: lat_centre_init_acyc = 60. +real :: init_vortex_radius_deg = 5. +real :: init_vortex_vor_f = 0.5 +real :: init_vortex_h_h_0 = 0.1 +logical :: add_initial_vortex_pair = .false. +logical :: add_initial_vortex_as_height = .true. + + +real, dimension(2) :: valid_range_v = (/-1.e3,1.e3/) + +namelist /shallow_dynamics_nml/ check_fourier_imag, & + south_to_north, triang_trunc, & + num_lon, num_lat, num_fourier, & + num_spherical, fourier_inc, & + longitude_origin, damping_option, & + damping_order, damping_coeff, & + robert_coeff, robert_coeff_tracer, & + h_0, spec_tracer, grid_tracer, & + valid_range_v, cutoff_wn, & + raw_filter_coeff, & + u_deep_mag, n_merid_deep_flow, & + u_upper_mag_init, & + lon_centre_init_cyc, & + lat_centre_init_cyc, & + lon_centre_init_acyc, & + lat_centre_init_acyc, & + init_vortex_radius_deg, & + init_vortex_vor_f, & + init_vortex_h_h_0, & + add_initial_vortex_pair, & + add_initial_vortex_as_height + + +contains + +!======================================================================================= + +subroutine shallow_dynamics_init (Dyn, Time, Time_init, dt_real) + +type(dynamics_type), intent(inout) :: Dyn +type(time_type) , intent(in) :: Time, Time_init +real , intent(in) :: dt_real + +integer :: i, j + +real, allocatable, dimension(:) :: glon_bnd, glat_bnd +real :: xx, yy, dd, deep_geopot_global_mean, radius_loc_cyc, radius_loc_acyc + +integer :: ierr, io, unit, id_lon, id_lat, id_lonb, id_latb +logical :: root + +! < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > + +call write_version_number (version, tagname) + +pe = mpp_pe() +root = (pe == mpp_root_pe()) + +if (file_exist('input.nml')) then + unit = open_namelist_file () + ierr=1 + do while (ierr /= 0) + read (unit, nml=shallow_dynamics_nml, iostat=io, end=10) + ierr = check_nml_error (io, 'shallow_dynamics_nml') + enddo + 10 call close_file (unit) +endif + +if (root) write (stdlog(), nml=shallow_dynamics_nml) + +call transforms_init(radius, num_lat, num_lon, num_fourier, fourier_inc, num_spherical, & + south_to_north=south_to_north, & + triang_trunc=triang_trunc, & + longitude_origin=longitude_origin ) + +call get_grid_domain(is,ie,js,je) +call get_spec_domain(ms,me,ns,ne) + +Dyn%num_lon = num_lon +Dyn%num_lat = num_lat +Dyn%spec_tracer = spec_tracer +Dyn%grid_tracer = grid_tracer + +allocate (sin_lat (js:je)) +allocate (cos_lat (js:je)) +allocate (deg_lat (js:je)) +allocate (deg_lon (is:ie)) +allocate (coriolis (js:je)) + +call get_deg_lon (deg_lon) +call get_deg_lat (deg_lat) +call get_sin_lat (sin_lat) +call get_cos_lat (cos_lat) + +allocate (glon_bnd (num_lon + 1)) +allocate (glat_bnd (num_lat + 1)) +call get_grid_boundaries (glon_bnd, glat_bnd, global=.true.) + +coriolis = 2*omega*sin_lat + +call spectral_damping_init(damping_coeff, damping_order, damping_option, cutoff_wn, num_fourier, num_spherical, 1, 0., 0., 0.) + +allocate(eigen(ms:me,ns:ne)) +call get_eigen_laplacian(eigen) + +allocate (Dyn%spec%vor (ms:me, ns:ne, num_time_levels)) +allocate (Dyn%spec%div (ms:me, ns:ne, num_time_levels)) +allocate (Dyn%spec%h (ms:me, ns:ne, num_time_levels)) + +allocate (Dyn%grid%u (is:ie, js:je, num_time_levels)) +allocate (Dyn%grid%v (is:ie, js:je, num_time_levels)) +allocate (Dyn%grid%vor (is:ie, js:je, num_time_levels)) +allocate (Dyn%grid%div (is:ie, js:je, num_time_levels)) +allocate (Dyn%grid%h (is:ie, js:je, num_time_levels)) + +allocate (Dyn%tend%u (is:ie, js:je)) +allocate (Dyn%tend%v (is:ie, js:je)) +allocate (Dyn%tend%h (is:ie, js:je)) +allocate (Dyn%grid%stream (is:ie, js:je)) +allocate (Dyn%grid%pv (is:ie, js:je)) +allocate (Dyn%grid%deep_geopot(is:ie, js:je)) + + +call fv_advection_init(num_lon, num_lat, glat_bnd, 360./float(fourier_inc)) +if(Dyn%grid_tracer) then + allocate(Dyn%Grid%tr (is:ie, js:je, num_time_levels)) + allocate(Dyn%Tend%tr (is:ie, js:je)) +endif + +if(Dyn%spec_tracer) then + allocate(Dyn%Grid%trs (is:ie, js:je, num_time_levels)) + allocate(Dyn%Tend%trs (is:ie, js:je)) + allocate(Dyn%Spec%trs (ms:me, ns:ne, num_time_levels)) +endif + +do i = is, ie + Dyn%grid%deep_geopot(i, js:je) = -2.*omega * u_deep_mag * radius * (1./(1.-n_merid_deep_flow**2.))*(-cos(n_merid_deep_flow*DEG_TO_RAD*deg_lat(js:je))*cos(DEG_TO_RAD*deg_lat(js:je)) - n_merid_deep_flow * (sin(n_merid_deep_flow*DEG_TO_RAD*deg_lat(js:je))*sin(DEG_TO_RAD*deg_lat(js:je))-sin(n_merid_deep_flow*(2.*atan(1.))))) +enddo + +deep_geopot_global_mean = area_weighted_global_mean(Dyn%grid%deep_geopot(:,:)) +Dyn%grid%deep_geopot(:,:) = Dyn%grid%deep_geopot(:,:)-deep_geopot_global_mean + +if(Time == Time_init) then + + Dyn%Grid%div(:,:,1) = 0.0 + Dyn%Grid%h (:,:,1) = h_0 - Dyn%grid%deep_geopot(:,:) + + do i = is, ie + Dyn%Grid%vor(i,js:je,1) = -((u_upper_mag_init * n_merid_deep_flow)/radius) * sin(DEG_to_RAD * deg_lat(js:je)) + + if (add_initial_vortex_pair) then + + do j=js, je + + radius_loc_cyc = ((min((deg_lon(i)-lon_centre_init_cyc)**2., (deg_lon(i)-lon_centre_init_cyc-360.)**2.)+(deg_lat(j)-lat_centre_init_cyc)**2.)**0.5)/init_vortex_radius_deg + radius_loc_acyc = ((min((deg_lon(i)-lon_centre_init_acyc)**2., (deg_lon(i)-lon_centre_init_acyc-360.)**2.)+(deg_lat(j)-lat_centre_init_acyc)**2.)**0.5)/init_vortex_radius_deg + + + if(radius_loc_cyc.le.1.0 .and. radius_loc_acyc.le.1.0) then + call error_mesg('shallow_dynamics','Cannot initialise cyclone and anticyclone in same grid box ', FATAL) + endif + + if(add_initial_vortex_as_height) then + if (radius_loc_cyc.le.2.0) then + Dyn%Grid%h(i,j,1) = Dyn%Grid%h(i,j,1) + init_vortex_h_h_0 * -h_0 * exp(-radius_loc_cyc**2.) + elseif (radius_loc_acyc.le.2.0) then + Dyn%Grid%h(i,j,1) = Dyn%Grid%h(i,j,1) + init_vortex_h_h_0 * h_0 * exp(-radius_loc_acyc**2.) + endif + else + if (radius_loc_cyc.le.1.0) then + Dyn%Grid%vor(i,j,1) = init_vortex_vor_f * 2.*omega + elseif (radius_loc_acyc.le.1.0) then + Dyn%Grid%vor(i,j,1) = init_vortex_vor_f * -2.*omega + endif + + endif + + enddo + + + endif + + enddo + + call trans_grid_to_spherical(Dyn%Grid%vor(:,:,1), Dyn%Spec%vor(:,:,1)) + call trans_grid_to_spherical(Dyn%Grid%div(:,:,1), Dyn%Spec%div(:,:,1)) + call trans_grid_to_spherical(Dyn%Grid%h (:,:,1), Dyn%Spec%h (:,:,1)) + + call uv_grid_from_vor_div (Dyn%Spec%vor(:,:,1), Dyn%Spec%div(:,:,1), & + Dyn%Grid%u (:,:,1), Dyn%Grid%v (:,:,1)) + + if(Dyn%grid_tracer) then + Dyn%Grid%tr = 0.0 + do j = js, je + if(deg_lat(j) > 10.0 .and. deg_lat(j) < 20.0) Dyn%Grid%tr(:,j,1) = 1.0 + if(deg_lat(j) > 70.0 ) Dyn%Grid%tr(:,j,1) = -1.0 + end do + endif + + if(Dyn%spec_tracer) then + Dyn%Grid%trs = 0.0 + do j = js, je + if(deg_lat(j) > 10.0 .and. deg_lat(j) < 20.0) Dyn%Grid%trs(:,j,1) = 1.0 + if(deg_lat(j) > 70.0 ) Dyn%Grid%trs(:,j,1) = -1.0 + end do + call trans_grid_to_spherical(Dyn%Grid%trs(:,:,1), Dyn%Spec%trs(:,:,1)) + endif + +else + + call read_restart(Dyn) + +endif + +module_is_initialized = .true. + +return +end subroutine shallow_dynamics_init + +!======================================================================================== + +subroutine shallow_dynamics(Time, Time_init, Dyn, previous, current, future, delta_t) + +type(time_type) , intent(in) :: Time, Time_init +type(dynamics_type), intent(inout) :: Dyn +integer, intent(in ) :: previous, current, future +real, intent(in ) :: delta_t + +! < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > + +complex, dimension(ms:me, ns:ne) :: dt_vors, dt_divs, dt_hs, stream, bs, work +real, dimension(is:ie, js:je) :: vorg, bg, h_future, h_dt, dt_vorg +integer :: j + +! < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > + +if(.not.module_is_initialized) then + call error_mesg('shallow_dynamics','dynamics has not been initialized ', FATAL) +endif + + +do j = js,je + vorg(:,j) = Dyn%Grid%vor(:,j,current) + coriolis(j) +end do +Dyn%Tend%u = Dyn%Tend%u + vorg*Dyn%Grid%v(:,:,current) +Dyn%Tend%v = Dyn%Tend%v - vorg*Dyn%Grid%u(:,:,current) + +call vor_div_from_uv_grid (Dyn%Tend%u, Dyn%Tend%v, dt_vors, dt_divs) + +call horizontal_advection (Dyn%Spec%h(:,:,current), & + Dyn%Grid%u(:,:,current), Dyn%Grid%v(:,:,current), Dyn%Tend%h) + +Dyn%Tend%h = Dyn%Tend%h - Dyn%Grid%h(:,:,current)*Dyn%Grid%div(:,:,current) + +call trans_grid_to_spherical (Dyn%Tend%h, dt_hs) + +bg = (Dyn%Grid%h(:,:,current) + Dyn%grid%deep_geopot(:,:) + & + 0.5*(Dyn%Grid%u(:,:,current)**2 + Dyn%Grid%v(:,:,current)**2)) + +call trans_grid_to_spherical(bg, bs) +dt_divs = dt_divs - compute_laplacian(bs) + +call implicit_correction (dt_divs, dt_hs, Dyn%Spec%div, Dyn%Spec%h, & + delta_t, previous, current) + +call compute_spectral_damping(Dyn%Spec%vor(:,:,previous), dt_vors, delta_t) +call compute_spectral_damping(Dyn%Spec%div(:,:,previous), dt_divs, delta_t) +call compute_spectral_damping(Dyn%Spec%h (:,:,previous), dt_hs , delta_t) + +call stirring(Time, dt_vors) + +call leapfrog(Dyn%Spec%vor , dt_vors , previous, current, future, delta_t, robert_coeff, raw_filter_coeff) +call leapfrog(Dyn%Spec%div , dt_divs , previous, current, future, delta_t, robert_coeff, raw_filter_coeff) +call leapfrog(Dyn%Spec%h , dt_hs , previous, current, future, delta_t, robert_coeff, raw_filter_coeff) + +call trans_spherical_to_grid(Dyn%Spec%vor(:,:,future), Dyn%Grid%vor(:,:,future)) +call trans_spherical_to_grid(Dyn%Spec%div(:,:,future), Dyn%Grid%div(:,:,future)) + + +call uv_grid_from_vor_div (Dyn%Spec%vor (:,:,future), Dyn%Spec%div(:,:,future), & + Dyn%Grid%u (:,:,future), Dyn%Grid%v (:,:,future)) + +call trans_spherical_to_grid (Dyn%Spec%h (:,:,future), Dyn%Grid%h(:,:,future)) + +if(minval(Dyn%Grid%v) < valid_range_v(1) .or. maxval(Dyn%Grid%v) > valid_range_v(2)) then + call error_mesg('shallow_dynamics','meridional wind out of valid range', FATAL) +endif + +if(Dyn%spec_tracer) call update_spec_tracer(Dyn%Spec%trs, Dyn%Grid%trs, Dyn%Tend%trs, & + Dyn%Grid%u, Dyn%Grid%v, previous, current, future, delta_t) + +if(Dyn%grid_tracer) call update_grid_tracer(Dyn%Grid%tr, Dyn%Tend%tr, & + Dyn%Grid%u, Dyn%Grid%v, previous, current, future, delta_t) + + +! for diagnostics + +stream = compute_laplacian(Dyn%Spec%vor(:,:,current), -1) ! for diagnostic purposes +call trans_spherical_to_grid(stream, Dyn%grid%stream) + +Dyn%Grid%pv = vorg/(Dyn%Grid%h(:,:,current)) + +return +end subroutine shallow_dynamics +!================================================================================ + +subroutine implicit_correction(dt_divs, dt_hs, divs, hs, delta_t, previous, current) + +complex, intent(inout), dimension(ms:,ns:) :: dt_divs, dt_hs +complex, intent(in), dimension(ms:,ns:,:) :: divs, hs +real , intent(in) :: delta_t +integer, intent(in) :: previous, current + +real :: xi, mu, mu2 + +xi = 0.5 ! centered implicit (for backwards implicit, set xi = 1.0) + +mu = xi*delta_t +mu2 = mu*mu + +dt_hs = dt_hs + h_0*(divs(:,:,current) - divs(:,:,previous)) +dt_divs = dt_divs - eigen*(hs(:,:,current) - hs(:,:,previous)) + +dt_divs = (dt_divs + mu*eigen*dt_hs)/(1.0 + mu2*eigen*h_0) +dt_hs = dt_hs - mu*h_0*dt_divs + +return +end subroutine implicit_correction + +!=================================================================================== + +subroutine update_spec_tracer(tr_spec, tr_grid, dt_tr, ug, vg, & + previous, current, future, delta_t) + +complex, intent(inout), dimension(ms:me, ns:ne, num_time_levels) :: tr_spec +real , intent(inout), dimension(is:ie, js:je, num_time_levels) :: tr_grid +real , intent(inout), dimension(is:ie, js:je ) :: dt_tr +real , intent(in ), dimension(is:ie, js:je, num_time_levels) :: ug, vg +real , intent(in ) :: delta_t +integer, intent(in ) :: previous, current, future + +complex, dimension(ms:me, ns:ne) :: dt_trs + +call horizontal_advection (tr_spec(:,:,current), ug(:,:,current), vg(:,:,current), dt_tr) +call trans_grid_to_spherical (dt_tr, dt_trs) +call compute_spectral_damping (tr_spec(:,:,previous), dt_trs, delta_t) +call leapfrog (tr_spec, dt_trs, previous, current, future, delta_t, robert_coeff, raw_filter_coeff) +call trans_spherical_to_grid (tr_spec(:,:,future), tr_grid(:,:,future)) + +return +end subroutine update_spec_tracer +!========================================================================== + +subroutine update_grid_tracer(tr_grid, dt_tr_grid, ug, vg, & + previous, current, future, delta_t) + +real , intent(inout), dimension(is:ie, js:je, num_time_levels) :: tr_grid +real , intent(inout), dimension(is:ie, js:je ) :: dt_tr_grid +real , intent(in ), dimension(is:ie, js:je, num_time_levels) :: ug, vg + +real , intent(in ) :: delta_t +integer, intent(in ) :: previous, current, future + +real, dimension(is:ie,js:je) :: tr_current, tr_future + +tr_future = tr_grid(:,:,previous) + delta_t*dt_tr_grid +dt_tr_grid = 0.0 +call a_grid_horiz_advection (ug(:,:,current), vg(:,:,current), tr_future, delta_t, dt_tr_grid) +tr_future = tr_future + delta_t*dt_tr_grid +tr_current = tr_grid(:,:,current) + & + robert_coeff_tracer*(tr_grid(:,:,previous) + tr_future - 2.0*tr_grid(:,:,current)) +tr_grid(:,:,current) = tr_current +tr_grid(:,:,future) = tr_future + +return +end subroutine update_grid_tracer + +!========================================================================== + +subroutine read_restart(Dyn) + +type(dynamics_type), intent(inout) :: Dyn + +integer :: unit, m, n, nt +real, dimension(ms:me, ns:ne) :: real_part, imag_part + +if(file_exist('INPUT/shallow_dynamics.res.nc')) then + do nt = 1, 2 + call read_data('INPUT/shallow_dynamics.res.nc', 'vors_real', real_part, spectral_domain, timelevel=nt) + call read_data('INPUT/shallow_dynamics.res.nc', 'vors_imag', imag_part, spectral_domain, timelevel=nt) + do n=ns,ne + do m=ms,me + Dyn%Spec%vor(m,n,nt) = cmplx(real_part(m,n),imag_part(m,n)) + end do + end do + call read_data('INPUT/shallow_dynamics.res.nc', 'divs_real', real_part, spectral_domain, timelevel=nt) + call read_data('INPUT/shallow_dynamics.res.nc', 'divs_imag', imag_part, spectral_domain, timelevel=nt) + do n=ns,ne + do m=ms,me + Dyn%Spec%div(m,n,nt) = cmplx(real_part(m,n),imag_part(m,n)) + end do + end do + call read_data('INPUT/shallow_dynamics.res.nc', 'hs_real', real_part, spectral_domain, timelevel=nt) + call read_data('INPUT/shallow_dynamics.res.nc', 'hs_imag', imag_part, spectral_domain, timelevel=nt) + do n=ns,ne + do m=ms,me + Dyn%Spec%h(m,n,nt) = cmplx(real_part(m,n),imag_part(m,n)) + end do + end do + if(Dyn%spec_tracer) then + call read_data('INPUT/shallow_dynamics.res.nc', 'trs_real', real_part, spectral_domain, timelevel=nt) + call read_data('INPUT/shallow_dynamics.res.nc', 'trs_imag', imag_part, spectral_domain, timelevel=nt) + do n=ns,ne + do m=ms,me + Dyn%Spec%trs(m,n,nt) = cmplx(real_part(m,n),imag_part(m,n)) + end do + end do + endif + call read_data('INPUT/shallow_dynamics.res.nc', 'u', Dyn%Grid%u (:,:,nt), grid_domain, timelevel=nt) + call read_data('INPUT/shallow_dynamics.res.nc', 'v', Dyn%Grid%v (:,:,nt), grid_domain, timelevel=nt) + call read_data('INPUT/shallow_dynamics.res.nc', 'vor', Dyn%Grid%vor(:,:,nt), grid_domain, timelevel=nt) + call read_data('INPUT/shallow_dynamics.res.nc', 'div', Dyn%Grid%div(:,:,nt), grid_domain, timelevel=nt) + call read_data('INPUT/shallow_dynamics.res.nc', 'h', Dyn%Grid%h (:,:,nt), grid_domain, timelevel=nt) + if(Dyn%spec_tracer) then + call read_data('INPUT/shallow_dynamics.res.nc', 'trs', Dyn%Grid%trs(:,:,nt), grid_domain, timelevel=nt) + endif + if(Dyn%grid_tracer) then + call read_data('INPUT/shallow_dynamics.res.nc', 'tr', Dyn%Grid%tr(:,:,nt), grid_domain, timelevel=nt) + endif + end do +else if(file_exist('INPUT/shallow_dynamics.res')) then + unit = open_restart_file(file='INPUT/shallow_dynamics.res',action='read') + + do nt = 1, 2 + call set_domain(spectral_domain) + call read_data(unit,Dyn%Spec%vor(:,:, nt)) + call read_data(unit,Dyn%Spec%div(:,:, nt)) + call read_data(unit,Dyn%Spec%h (:,:, nt)) + if(Dyn%spec_tracer) call read_data(unit,Dyn%Spec%trs(:,:, nt)) + + call set_domain(grid_domain) + call read_data(unit,Dyn%Grid%u (:,:, nt)) + call read_data(unit,Dyn%Grid%v (:,:, nt)) + call read_data(unit,Dyn%Grid%vor (:,:, nt)) + call read_data(unit,Dyn%Grid%div (:,:, nt)) + call read_data(unit,Dyn%Grid%h (:,:, nt)) + if(Dyn%spec_tracer) call read_data(unit,Dyn%Grid%trs(:,:, nt)) + if(Dyn%grid_tracer) call read_data(unit,Dyn%Grid%tr (:,:, nt)) + + end do + call close_file(unit) + +else + call error_mesg('read_restart', 'restart does not exist', FATAL) +endif + +return +end subroutine read_restart + +!==================================================================== + +subroutine write_restart(Dyn, previous, current) + +type(dynamics_type), intent(in) :: Dyn +integer, intent(in) :: previous, current + +integer :: unit, nt, nn + +do nt = 1, 2 + if(nt == 1) nn = previous + if(nt == 2) nn = current + call write_data('RESTART/shallow_dynamics.res.nc', 'vors_real', real(Dyn%Spec%vor(:,:,nn)), spectral_domain) + call write_data('RESTART/shallow_dynamics.res.nc', 'vors_imag', aimag(Dyn%Spec%vor(:,:,nn)), spectral_domain) + call write_data('RESTART/shallow_dynamics.res.nc', 'divs_real', real(Dyn%Spec%div(:,:,nn)), spectral_domain) + call write_data('RESTART/shallow_dynamics.res.nc', 'divs_imag', aimag(Dyn%Spec%div(:,:,nn)), spectral_domain) + call write_data('RESTART/shallow_dynamics.res.nc', 'hs_real', real(Dyn%Spec%h (:,:,nn)), spectral_domain) + call write_data('RESTART/shallow_dynamics.res.nc', 'hs_imag', aimag(Dyn%Spec%h (:,:,nn)), spectral_domain) + if(Dyn%spec_tracer) then + call write_data('RESTART/shallow_dynamics.res.nc', 'trs_real', real(Dyn%Spec%trs(:,:,nn)), spectral_domain) + call write_data('RESTART/shallow_dynamics.res.nc', 'trs_imag', aimag(Dyn%Spec%trs(:,:,nn)), spectral_domain) + endif + call write_data('RESTART/shallow_dynamics.res.nc', 'u', Dyn%Grid%u (:,:,nn), grid_domain) + call write_data('RESTART/shallow_dynamics.res.nc', 'v', Dyn%Grid%v (:,:,nn), grid_domain) + call write_data('RESTART/shallow_dynamics.res.nc', 'vor', Dyn%Grid%vor(:,:,nn), grid_domain) + call write_data('RESTART/shallow_dynamics.res.nc', 'div', Dyn%Grid%div(:,:,nn), grid_domain) + call write_data('RESTART/shallow_dynamics.res.nc', 'h', Dyn%Grid%h (:,:,nn), grid_domain) + if(Dyn%spec_tracer) then + call write_data('RESTART/shallow_dynamics.res.nc', 'trs', Dyn%Grid%trs(:,:,nn), grid_domain) + endif + if(Dyn%grid_tracer) then + call write_data('RESTART/shallow_dynamics.res.nc', 'tr', Dyn%Grid%tr(:,:,nn), grid_domain) + endif +enddo + +!unit = open_restart_file(file='RESTART/shallow_dynamics.res', action='write') + +!do n = 1, 2 +! if(n == 1) nn = previous +! if(n == 2) nn = current +! +! call set_domain(spectral_domain) +! call write_data(unit,Dyn%Spec%vor(:,:, nn)) +! call write_data(unit,Dyn%Spec%div(:,:, nn)) +! call write_data(unit,Dyn%Spec%h (:,:, nn)) +! if(Dyn%spec_tracer) call write_data(unit,Dyn%Spec%trs(:,:, nn)) +! +! call set_domain(grid_domain) +! call write_data(unit,Dyn%Grid%u (:,:, nn)) +! call write_data(unit,Dyn%Grid%v (:,:, nn)) +! call write_data(unit,Dyn%Grid%vor (:,:, nn)) +! call write_data(unit,Dyn%Grid%div (:,:, nn)) +! call write_data(unit,Dyn%Grid%h (:,:, nn)) +! if(Dyn%spec_tracer) call write_data(unit,Dyn%Grid%trs(:,:, nn)) +! if(Dyn%grid_tracer) call write_data(unit,Dyn%Grid%tr (:,:, nn)) +! +!end do + +!call close_file(unit) + +end subroutine write_restart + +!==================================================================== + +subroutine shallow_dynamics_end (Dyn, previous, current) + +type(dynamics_type), intent(inout) :: Dyn +integer, intent(in) :: previous, current + +if(.not.module_is_initialized) then + call error_mesg('shallow_dynamics_end','dynamics has not been initialized ', FATAL) +endif + +call write_restart (Dyn, previous, current) + +call transforms_end +call stirring_end + +module_is_initialized = .false. + +return +end subroutine shallow_dynamics_end +!=================================================================================== + +end module shallow_dynamics_mod diff --git a/src/atmos_spectral_shallow/shallow_dynamics.html b/src/atmos_spectral_shallow/shallow_dynamics.html new file mode 100644 index 000000000..600933a30 --- /dev/null +++ b/src/atmos_spectral_shallow/shallow_dynamics.html @@ -0,0 +1,370 @@ + +module shallow_dynamics_mod + + +
+ + +module shallow_dynamics_mod
+ + ++ Contact: Isaac Held + Reviewers: Peter Phillipps + ++ + + +
+OVERVIEW
+ ++ + The dynamical core of the spectral transform model for + the shallow water equations on the sphere. + ++ + + +
+DESCRIPTION
+ ++ + Integrates the shallow water equation for hydrostatic flow of a homgeoneous, + incompressible fluid on the + sphere using the spectral transform technique. Also allows for the + inclusion of a passive tracer advected by the the spectral advection + algorithm, and a gridpoint tracer advected with a finite + volume algorithm on the transform grid. Thinking of the model as one of + the upper tropopsheric flow, the default experiment involves relaxation of + the geopotential to an "equilibrium value" with maxima (whose amplitude + and shape are controlled from the namelist) along the equator and in the + subtropicals. + + For a full description of the model and algorithms used, see + shallow.ps + + For higher level routines for running this shallow water model, + see atmosphere_mod + + ++ + + +
+OTHER MODULES USED
+ ++ + fms_mod + constants_mod + time_manager_mod + transforms_mod + spectral_damping_mod + leapfrog_mod + fv_advection_mod + ++ + + +
+PUBLIC INTERFACE
+ ++ + use shallow_dynamics_mod [,only: shallow_dynamics_init, + shallow_dynamics, + shallow_dynamics_end, + dynamics_type, + grid_type, + spectral_type, + tendency_type] + ++ + + +
+PUBLIC DATA
+ ++ ++ + + +
+ + +type grid_type + real, pointer, dimension(:,:,:) :: u, v, vor, div, h, trs, tr + real, pointer, dimension(:,:) :: stream, pv +end type + + allocated space for grid fields + + (:,:,:) => (lon, lat, time_level) + (:,:) => (lon, lat) + (lon, lat) on local computational domain + time_level stores the two time levels needed for the + leapfrog step + + u -- eastward velocity (m/s) + v -- northward velocity (m/s) + vor -- vorticity (1/s) + div -- divergence (1/s) + h -- geopotential (m^2/s^2) + trs -- tracer advected spectrally + tr -- tracer advected on grid + pv -- (f + vor)/h, where f = 2*omega*sin(lat) (s/m^2) + stream -- streamfunction (m^2/s) at current time + + + +
+ + +type spectral_type + complex, pointer, dimension(:,:,:) :: vor, div, h, trs +end type + + allocated space for spectral fields + + (:,:,:) => (zonal, meridional, time_level) + + vor -- spectral vorticity + div -- spectral divergence + h -- spectral geopotential + trs -- spectral tracer + +
+ + +type tendency_type + real, pointer, dimension(:,:) :: u, v, h, trs, tr +end type + + allocated space for accumulating tendencies, d/dt, in grid space, + for prognostic variables + + (:,:,:) => (lon, lat) + +
+ + +type dynamics_type + type(grid_type) :: grid + type(spectral_type) :: spec + type(tendency_type) :: tend + integer :: num_lon, num_lat ! size of global domain + logical :: grid_tracer, spec_tracer +end type + + grid_tracer = .true. => tracer with gridpoint advection is beign integrated + similarly for spec_tracer + +
+ +
+PUBLIC ROUTINES
+ ++ +subroutine shallow_dynamics_init +subroutine shallow _dynamics +subroutine shallow_dynamics_end +type (grid_type) +type (spectral_type) +type (tendency_type) +type (dynamics_type) + + + ++ + + ++ + subroutine shallow_dynamics_init(Dyn, Time, Time_init) + + type(dynamics_type), intent(inout) :: Dyn + type containing all dynamical fields and related information + (see type (dynamics_type)) + + type(time_type) , intent(in) :: Time, Time_init + current time and time at which integeration began + time_type defined by time_manager_mod + + + Initializes the module; + Reads restart from 'INPUT/shallow_dynamics.res' if Time = Time_init; + otherwise uses default initial conditions + ++ +
+ + + + subroutine shallow_dynamics & + (Time, Time_init, Dyn, previous, current, future, delta_t) + + type(time_type) , intent(inout) :: Time, Time_init + type(dynamics_type), intent(inout) :: Dyn + integer , intent(in ) :: previous, current, future + real , intent(in ) :: delta_t + + previous, current and future = 1 or 2 + these integers refer to the third dimension of the + three-dimensional fields in Dyn + the fields at time t - delta_t are assumed to be in (:,:,previous) + the fields at time t are assumed to be in (:,:,current) + the fields at time t + delta_t are placed in (:,:,future) + overwriting whatever is already there + + delta_t = time step in seconds + + updates dynamical fields by one time step + + + +
+ + + subroutine shallow_dynamics_end(Dyn, previous, current) + + type(dynamics_type), intent(inout) :: Dyn + integer, intent(in) :: previous, current + + + Terminates module; + writes restart file to 'RESTART/shallow_dynamics.res' + + + + +
+
+NAMELIST
+ ++ +&shallow_dynamics_nml + + integer :: num_lat = 128 + number of latitudes in global grid + + integer :: num_lon = 256 + number of longitudes in global grid + should equal 2*num_lat for Triangular truncation + + integer :: num_fourier = 85 + the retained fourier wavenumber are n*fourier_inc, where + n ranges from 0 to num_fourier + + integer :: num_spherical = 86 + the maximum number of meridional modes for any zonal wavenumber + for triangular truncation, set num_spherical = num_fourier +1 + + integer :: fourier_inc = 1 + creates a "sector" model if fourier_inc > 1; integration domain is + (360 degrees longitude)/fourier_inc + + (the default values listed above define a standard T85 model) + + logical :: check_fourier_imag = .false. + if true, checks to see if fields to be transformed to grid + domain have zero imaginary part to their zonally symmetric + modes; useful for debugging + + logical :: south_to_north = .true. + true => grid runs from south to north + false => grid runs from north to south + + logical :: triangular_trunc = .true. + true => shape of truncation is triangular + false => shape of truncation is rhomboidal + + real :: robert_coeff = 0.04 + x(current) => (1-2r)*x(current) + r*(x(future)+x(previous)) + where r = robert_coeff (non-dimensional) + + real :: robert_coeff_tracer = 0.04 + (same as robert_coeff, but for grid tracer) + + real :: longitude_origin = 0.0 + longitude of first longitude, in degrees + (if you want the westgern boundary of first grid boc to be at + 0.0, set longitude_origin = 0.5*360./float(num_lon)) + + integer :: damping_option = 'resolution_dependent' + integer :: damping_order = 4 + real :: damping_coeff = 1.e-04 + + damping = nu*(del^2)^n where n = damping order + damping_option = 'resolution_dependent' or 'resolution_independent' + = 'resolution_dependent' => nu is set so that the damping rate for the + mode (m=0,n=num_spherical-1) equals damping_coeff (in 1/s) + For triangular truncation, damping_coeff is then the + rate of damping of the highest retained mode + + = 'resolution_independent' => nu = damping_coef + + real :: h_0 = 3.e04 + (m^2)/(s^2) + the initial condition is a state of rest with geopotential = h_0 + (h_0 is also used to determine the part of the divergence equation + that is integrated implicitly) + + logical :: spec_tracer = .true. + logical :: grid_tracer = .true. + spec_tracer = true => a passive tracer is carried that is advected + spectrally, with the same algorithm as the vorticity + grid_tracer = ture => a passive tracer is carried that is advected + on the spectral transform grid by a finite-volume algorithm + (see shallow.ps ) + Both tracers can be carried simultaeneously + + real, dimension(2) :: valid_range_v = -1000., 1000. + A valid range for meridional wind. Model terminates if meridional wind + goes outside the valid range. Allows model to terminate gracefully when, + for example, the model becomes numerically unstable. + ++ + + +
+ERROR MESSAGES
+ ++ + "Dynamics has not been initialized" + -- shallow_dynamics_init must be called before any other + routines in the module are called + + "restart does not exist" + -- Time is not equal to Time_init at initalization, but the file + 'INPUT/shallow_dynamics.res' does not exit + + ++ + + +
+ + diff --git a/src/atmos_spectral_shallow/shallow_physics.F90 b/src/atmos_spectral_shallow/shallow_physics.F90 new file mode 100644 index 000000000..8f8896d2e --- /dev/null +++ b/src/atmos_spectral_shallow/shallow_physics.F90 @@ -0,0 +1,232 @@ +module shallow_physics_mod + +!----------------------------------------------------------------------- +! GNU General Public License +! +! This program 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; either version 2 of +! the License, or (at your option) any later version. +! +! This program 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. +! +! For the full text of the GNU General Public License, +! write to: Free Software Foundation, Inc., +! 675 Mass Ave, Cambridge, MA 02139, USA. +! or see: http://www.gnu.org/licenses/gpl.html +!----------------------------------------------------------------------- + +use fms_mod, only: open_namelist_file, & + open_restart_file, & + file_exist, & + check_nml_error, & + error_mesg, & + FATAL, WARNING, & + write_version_number, & + mpp_pe, & + mpp_root_pe, & + fms_init, fms_end, & + read_data, & + write_data, & + set_domain, & + close_file, & + stdlog + +use transforms_mod, only: get_sin_lat, get_cos_lat, & + get_deg_lon, get_deg_lat, & + get_wts_lat, & + get_grid_domain, get_spec_domain, & + grid_domain + +use time_manager_mod, only: time_type + +!======================================================================== +implicit none +private +!======================================================================== + +public :: shallow_physics_init, & + shallow_physics, & + shallow_physics_end, & + phys_type + + +! version information +!======================================================================== +character(len=128) :: version = '$Id: shallow_physics.F90,v 10.0 2003/10/24 22:01:02 fms Exp $' +character(len=128) :: tagname = '$Name: siena_201207 $' +!======================================================================== + +type phys_type + real, pointer, dimension(:,:) :: empty=>NULL() +end type + +logical :: module_is_initialized = .false. + +integer :: is, ie, js, je + +integer :: pe +logical :: root + +real, allocatable, dimension(:) :: rad_lat, deg_lat, deg_lon, & + sin_lat, cos_lat, wts_lat + +real, allocatable, dimension(:,:) :: h_eq + +real :: kappa_m, kappa_t + + + +! namelist +!======================================================================== + +real :: fric_damp_time = -20.0 +real :: therm_damp_time = -10.0 +real :: del_h = 1.e04 +real :: h_0 = 3.e04 +real :: h_amp = 2.e04 +real :: h_lon = 90.0 +real :: h_lat = 25.0 +real :: h_width = 15.0 +real :: h_itcz = 1.e05 +real :: itcz_width = 4.0 + +namelist /shallow_physics_nml/ fric_damp_time, therm_damp_time, del_h, h_0, & + h_amp, h_lon, h_lat, h_width, & + itcz_width, h_itcz +!======================================================================== + +contains + +!======================================================================== + +subroutine shallow_physics_init(Phys) + +type(phys_type), intent(inout) :: Phys + +integer :: i, j, unit, ierr, io + +real :: xx, yy, dd + +call write_version_number(version, tagname) + +pe = mpp_pe() +root = (pe == mpp_root_pe()) + +! read the namelist + +if (file_exist('input.nml')) then + unit = open_namelist_file () + ierr=1 + do while (ierr /= 0) + read (unit, nml=shallow_physics_nml, iostat=io, end=10) + ierr = check_nml_error (io, 'shallow_physics_nml') + enddo + 10 call close_file (unit) +endif + +if(fric_damp_time < 0.0) fric_damp_time = - fric_damp_time*86400 +if(therm_damp_time < 0.0) therm_damp_time = - therm_damp_time*86400 + +kappa_m = 0.0 +kappa_t = 0.0 +if( fric_damp_time .ne. 0.0) kappa_m = 1./fric_damp_time +if(therm_damp_time .ne. 0.0) kappa_t = 1./therm_damp_time + +call get_grid_domain(is,ie,js,je) + +allocate ( rad_lat (js:je) ) +allocate ( deg_lat (js:je) ) +allocate ( sin_lat (js:je) ) +allocate ( cos_lat (js:je) ) +allocate ( wts_lat (js:je) ) +allocate ( deg_lon (is:ie) ) +allocate ( h_eq (is:ie,js:je) ) + +call get_wts_lat(wts_lat) +call get_deg_lat(deg_lat) +call get_deg_lon(deg_lon) +rad_lat = deg_lat*atan(1.)/45. +sin_lat = sin(rad_lat) +cos_lat = cos(rad_lat) + + +do j = js, je + do i = is, ie + xx = (deg_lon(i) - h_lon)/(h_width*2.0) + yy = (deg_lat(j) - h_lat)/h_width + dd = xx*xx + yy*yy + h_eq(i,j) = h_0 + h_amp*max(1.e-10, exp(-dd)) + end do +end do + +do j = js, je + yy = deg_lat(j)/itcz_width + dd = yy*yy + h_eq(:,j) = h_eq(:,j) + h_itcz*exp(-dd) +end do + +!if(file_exist('INPUT/shallow_physics.res')) then +! unit = open_restart_file(file='INPUT/shallow_physics.res',action='read') +! call set_domain(grid_domain) +! call close_file(unit) +!else + +!endif + +module_is_initialized = .true. + +return +end subroutine shallow_physics_init + +!======================================================================= + +subroutine shallow_physics(Time, dt_ug, dt_vg, dt_hg, ug, vg, hg, & + delta_t, previous, current, Phys) + +real, intent(inout), dimension(is:ie, js:je) :: dt_ug, dt_vg, dt_hg +real, intent(in) , dimension(is:ie, js:je, 2) :: ug, vg, hg + +real , intent(in) :: delta_t +integer, intent(in) :: previous, current + +type(time_type), intent(in) :: Time +type(phys_type), intent(inout) :: Phys + +dt_ug = dt_ug - kappa_m*ug(:,:,previous) +dt_vg = dt_vg - kappa_m*vg(:,:,previous) +dt_hg = dt_hg - kappa_t*(hg(:,:,previous) - h_eq) + + +return +end subroutine shallow_physics + +!====================================================================== + +subroutine shallow_physics_end(Phys) + +type(phys_type), intent(in) :: Phys + +integer :: unit + +if(.not.module_is_initialized) then + call error_mesg('shallow_physics_end','physics has not been initialized ', FATAL) +endif + +!unit = open_restart_file(file='RESTART/shallow_physics.res', action='write') + +!call set_domain(grid_domain) + +!call close_file(unit) + +module_is_initialized = .false. + +return +end subroutine shallow_physics_end + +!====================================================================== + +end module shallow_physics_mod diff --git a/src/atmos_spectral_shallow/shallow_physics.html b/src/atmos_spectral_shallow/shallow_physics.html new file mode 100644 index 000000000..d3ad34980 --- /dev/null +++ b/src/atmos_spectral_shallow/shallow_physics.html @@ -0,0 +1,210 @@ + +module shallow_physics_mod + + +
+ + +module shallow_physics_mod
+ + ++ Contact: Isaac Held + Reviewers: Peter Phillipps + ++ + + +
+OVERVIEW
+ ++ + A module that allows one to add processes that act in the grid domain + to the dynamics of the shallow model on the sphere + ++ + + +
+DESCRIPTION
+ ++ + A module that allows one to add processes that act in the grid domain + to the dynamics of the shallow model on the sphere. Currently adds + a relaxation to a specified "equilibrium geopotential" and relaxes + the winds to zero ++ + + +
+OTHER MODULES USED
+ ++ + fms_mod + transforms_mod + time_manager_mod + ++ + + +
+PUBLIC INTERFACE
+ ++ + use shallow_physics_mod [,only: shallow_physics_init, + shallow_physics, + shallow_physics_end, + phys_type] + ++ + + +
+PUBLIC DATA
+ ++ ++ + + +
+ + +type phys_type + real, pointer, dimension(:,:) :: empty +end type + + fields from physics module made available for diagnostics + +
+ + +
+PUBLIC ROUTINES
+ ++ +subroutine shallow_physics_init +subroutine shallow_physics +subroutine shallow_physics_end +type(phys_type) + + + ++ + + + ++ + subroutine shallow_physics_init(Phys) + + type(phys_type) , intent(inout) :: Phys + + + Initializes module + ++ +
+ + + + subroutine shallow_physics (Time, dt_ug, dt_vg, dt_hg, ug, vg, hg, & + delta_t, previous, current, Phys) + + real, intent(inout), dimension(:,:) :: dt_ug, dt_vg, dt_hg + + the u, v and geopotential tendencies onto which tendencies due to + the grid-point physics are added (m/(s^2) for dt_ug, dt_vg; + (m^2)/(s^3) for dt_hg) + + real, intent(in) , dimension(:,:, 2) :: ug, vg, hg + the grid zonal and meridional velocities (m/s) and + geopotential (m^2/s^2) + the third index is the time-index used in the leapfrog step + + real , intent(in) :: delta_t + time step (s) + + integer, intent(in) :: previous, current + = 1 or 2 + ug(:,:,previous) is the velocity at t-delta_t + ug(:,:,current ) is the velocity at t + + type(time_type), intent(in) :: Time + type(phys_type), intent(inout) :: Phys + + + + +
+ + + + subroutine shallow_physics_end (Phys) + + type(phys_type), intent(inout) :: Phys + + + +
+
+NAMELIST
+ ++ +&shallow_physics_nml + +real :: fric_damp_time = -20.0 + rate at which ua nd v are relaxed to zero (seconds) + (if negative, units are days instead -- negative sign is ignored) + +real :: therm_damp_time = -10.0 + rate at which geopotential is relaxed to h_eq + (units as above) + +real :: h_0 = 3.e04 (m^2/s^2) +real :: h_amp = 2.e04 (m^2/s^2) +real :: h_lon = 90.0 degrees +real :: h_lat = 25.0 degrees +real :: h_width = 15.0 degrees +real :: h_itcz = 1.e05 (m^2/s^2) +real :: itcz_width = 4.0 degrees + + h_eq is defined as + h_0 + h_amp*exp(-r^2) + h_itcz*exp(-d^2) + + where r^2 = xx^2 + yy^2 + xx = (lon - h_lon)/(2*h_width) + yy = (lat - h_lat)/h_width + + and d = lat/itcz_width + + ++ + + +
+ + + diff --git a/src/extra/model/barotropic/field_table b/src/extra/model/barotropic/field_table new file mode 100644 index 000000000..e69de29bb diff --git a/src/extra/model/barotropic/path_names b/src/extra/model/barotropic/path_names new file mode 100644 index 000000000..e7c619d86 --- /dev/null +++ b/src/extra/model/barotropic/path_names @@ -0,0 +1,138 @@ +atmos_solo/atmos_model.F90 +atmos_spectral_barotropic/atmosphere.F90 +atmos_spectral_barotropic/barotropic_diagnostics.F90 +atmos_spectral_barotropic/barotropic_dynamics.F90 +atmos_spectral_barotropic/barotropic_physics.F90 +atmos_spectral_barotropic/stirring.F90 +atmos_spectral/model/fv_advection.F90 +atmos_spectral/model/leapfrog.F90 +atmos_spectral/model/spectral_damping.F90 +atmos_spectral/tools/gauss_and_legendre.F90 +atmos_spectral/tools/grid_fourier.F90 +atmos_spectral/tools/spec_mpp.F90 +atmos_spectral/tools/spherical.F90 +atmos_spectral/tools/spherical_fourier.F90 +atmos_spectral/tools/transforms.F90 +shared/constants/constants.F90 +shared/diag_manager/diag_axis.F90 +shared/diag_manager/diag_data.F90 +shared/diag_manager/diag_grid.F90 +shared/diag_manager/diag_manager.F90 +shared/diag_manager/diag_output.F90 +shared/diag_manager/diag_table.F90 +shared/diag_manager/diag_util.F90 +shared/fft/fft99.F90 +shared/fft/fft.F90 +shared/field_manager/field_manager.F90 +shared/field_manager/fm_util.F90 +shared/field_manager/parse.inc +shared/fms/fms.F90 +shared/fms/fms_io.F90 +shared/fms/read_data_2d.inc +shared/fms/read_data_3d.inc +shared/fms/read_data_4d.inc +shared/fms/test_fms_io.F90 +shared/fms/write_data.inc +shared/include/fms_platform.h +shared/memutils/memuse.c +shared/memutils/memutils.F90 +shared/mosaic/constant.h +shared/mosaic/create_xgrid.c +shared/mosaic/create_xgrid.h +shared/mosaic/gradient_c2l.c +shared/mosaic/gradient_c2l.h +shared/mosaic/gradient.F90 +shared/mosaic/grid.F90 +shared/mosaic/interp.c +shared/mosaic/interp.h +shared/mosaic/mosaic.F90 +shared/mosaic/mosaic_util.c +shared/mosaic/mosaic_util.h +shared/mosaic/read_mosaic.c +shared/mosaic/read_mosaic.h +shared/mpp/affinity.c +shared/mpp/include/mpp_chksum.h +shared/mpp/include/mpp_chksum_int.h +shared/mpp/include/mpp_chksum_scalar.h +shared/mpp/include/mpp_comm.inc +shared/mpp/include/mpp_comm_mpi.inc +shared/mpp/include/mpp_comm_nocomm.inc +shared/mpp/include/mpp_comm_sma.inc +shared/mpp/include/mpp_data_mpi.inc +shared/mpp/include/mpp_data_nocomm.inc +shared/mpp/include/mpp_data_sma.inc +shared/mpp/include/mpp_define_nest_domains.inc +shared/mpp/include/mpp_do_check.h +shared/mpp/include/mpp_do_checkV.h +shared/mpp/include/mpp_do_get_boundary.h +shared/mpp/include/mpp_do_global_field.h +shared/mpp/include/mpp_domains_comm.inc +shared/mpp/include/mpp_domains_define.inc +shared/mpp/include/mpp_domains_misc.inc +shared/mpp/include/mpp_domains_reduce.inc +shared/mpp/include/mpp_domains_util.inc +shared/mpp/include/mpp_do_redistribute.h +shared/mpp/include/mpp_do_update_ad.h +shared/mpp/include/mpp_do_update.h +shared/mpp/include/mpp_do_update_nest.h +shared/mpp/include/mpp_do_update_nonblock.h +shared/mpp/include/mpp_do_updateV_ad.h +shared/mpp/include/mpp_do_updateV.h +shared/mpp/include/mpp_do_updateV_nonblock.h +shared/mpp/include/mpp_error_a_a.h +shared/mpp/include/mpp_error_a_s.h +shared/mpp/include/mpp_error_s_a.h +shared/mpp/include/mpp_error_s_s.h +shared/mpp/include/mpp_gather.h +shared/mpp/include/mpp_get_boundary.h +shared/mpp/include/mpp_global_field.h +shared/mpp/include/mpp_global_reduce.h +shared/mpp/include/mpp_global_sum_ad.h +shared/mpp/include/mpp_global_sum.h +shared/mpp/include/mpp_global_sum_tl.h +shared/mpp/include/mpp_io_connect.inc +shared/mpp/include/mpp_io_misc.inc +shared/mpp/include/mpp_io_read.inc +shared/mpp/include/mpp_io_util.inc +shared/mpp/include/mpp_io_write.inc +shared/mpp/include/mpp_read_2Ddecomp.h +shared/mpp/include/mpp_reduce_mpi.h +shared/mpp/include/mpp_reduce_nocomm.h +shared/mpp/include/mpp_reduce_sma.h +shared/mpp/include/mpp_sum.inc +shared/mpp/include/mpp_sum_mpi.h +shared/mpp/include/mpp_sum_nocomm.h +shared/mpp/include/mpp_sum_sma.h +shared/mpp/include/mpp_transmit.inc +shared/mpp/include/mpp_transmit_mpi.h +shared/mpp/include/mpp_transmit_nocomm.h +shared/mpp/include/mpp_transmit_sma.h +shared/mpp/include/mpp_update_domains2D_ad.h +shared/mpp/include/mpp_update_domains2D.h +shared/mpp/include/mpp_update_domains2D_nonblock.h +shared/mpp/include/mpp_update_nest_domains.h +shared/mpp/include/mpp_util.inc +shared/mpp/include/mpp_util_mpi.inc +shared/mpp/include/mpp_util_nocomm.inc +shared/mpp/include/mpp_util_sma.inc +shared/mpp/include/mpp_write_2Ddecomp.h +shared/mpp/include/mpp_write.h +shared/mpp/include/system_clock.h +shared/mpp/mpp_data.F90 +shared/mpp/mpp_domains.F90 +shared/mpp/mpp.F90 +shared/mpp/mpp_io.F90 +shared/mpp/mpp_memutils.F90 +shared/mpp/mpp_parameter.F90 +shared/mpp/mpp_pset.F90 +shared/mpp/mpp_utilities.F90 +shared/mpp/nsclock.c +shared/mpp/test_mpp_domains.F90 +shared/mpp/test_mpp.F90 +shared/mpp/test_mpp_io.F90 +shared/mpp/test_mpp_pset.F90 +shared/mpp/threadloc.c +shared/platform/platform.F90 +shared/time_manager/get_cal_time.F90 +shared/time_manager/time_manager.F90 +shared/tracer_manager/tracer_manager.F90 diff --git a/src/extra/model/shallow/path_names b/src/extra/model/shallow/path_names new file mode 100644 index 000000000..14f1c56f1 --- /dev/null +++ b/src/extra/model/shallow/path_names @@ -0,0 +1,138 @@ +atmos_solo/atmos_model.F90 +atmos_spectral/model/fv_advection.F90 +atmos_spectral/model/leapfrog.F90 +atmos_spectral/model/spectral_damping.F90 +atmos_spectral_shallow/atmosphere.F90 +atmos_spectral_shallow/shallow_diagnostics.F90 +atmos_spectral_shallow/shallow_dynamics.F90 +atmos_spectral_shallow/shallow_physics.F90 +atmos_spectral_barotropic/stirring.F90 +atmos_spectral/tools/gauss_and_legendre.F90 +atmos_spectral/tools/grid_fourier.F90 +atmos_spectral/tools/spec_mpp.F90 +atmos_spectral/tools/spherical.F90 +atmos_spectral/tools/spherical_fourier.F90 +atmos_spectral/tools/transforms.F90 +shared/constants/constants.F90 +shared/diag_manager/diag_axis.F90 +shared/diag_manager/diag_data.F90 +shared/diag_manager/diag_grid.F90 +shared/diag_manager/diag_manager.F90 +shared/diag_manager/diag_output.F90 +shared/diag_manager/diag_table.F90 +shared/diag_manager/diag_util.F90 +shared/fft/fft99.F90 +shared/fft/fft.F90 +shared/field_manager/field_manager.F90 +shared/field_manager/fm_util.F90 +shared/field_manager/parse.inc +shared/fms/fms.F90 +shared/fms/fms_io.F90 +shared/fms/read_data_2d.inc +shared/fms/read_data_3d.inc +shared/fms/read_data_4d.inc +shared/fms/test_fms_io.F90 +shared/fms/write_data.inc +shared/include/fms_platform.h +shared/memutils/memuse.c +shared/memutils/memutils.F90 +shared/mosaic/constant.h +shared/mosaic/create_xgrid.c +shared/mosaic/create_xgrid.h +shared/mosaic/gradient_c2l.c +shared/mosaic/gradient_c2l.h +shared/mosaic/gradient.F90 +shared/mosaic/grid.F90 +shared/mosaic/interp.c +shared/mosaic/interp.h +shared/mosaic/mosaic.F90 +shared/mosaic/mosaic_util.c +shared/mosaic/mosaic_util.h +shared/mosaic/read_mosaic.c +shared/mosaic/read_mosaic.h +shared/mpp/affinity.c +shared/mpp/include/mpp_chksum.h +shared/mpp/include/mpp_chksum_int.h +shared/mpp/include/mpp_chksum_scalar.h +shared/mpp/include/mpp_comm.inc +shared/mpp/include/mpp_comm_mpi.inc +shared/mpp/include/mpp_comm_nocomm.inc +shared/mpp/include/mpp_comm_sma.inc +shared/mpp/include/mpp_data_mpi.inc +shared/mpp/include/mpp_data_nocomm.inc +shared/mpp/include/mpp_data_sma.inc +shared/mpp/include/mpp_define_nest_domains.inc +shared/mpp/include/mpp_do_check.h +shared/mpp/include/mpp_do_checkV.h +shared/mpp/include/mpp_do_get_boundary.h +shared/mpp/include/mpp_do_global_field.h +shared/mpp/include/mpp_domains_comm.inc +shared/mpp/include/mpp_domains_define.inc +shared/mpp/include/mpp_domains_misc.inc +shared/mpp/include/mpp_domains_reduce.inc +shared/mpp/include/mpp_domains_util.inc +shared/mpp/include/mpp_do_redistribute.h +shared/mpp/include/mpp_do_update_ad.h +shared/mpp/include/mpp_do_update.h +shared/mpp/include/mpp_do_update_nest.h +shared/mpp/include/mpp_do_update_nonblock.h +shared/mpp/include/mpp_do_updateV_ad.h +shared/mpp/include/mpp_do_updateV.h +shared/mpp/include/mpp_do_updateV_nonblock.h +shared/mpp/include/mpp_error_a_a.h +shared/mpp/include/mpp_error_a_s.h +shared/mpp/include/mpp_error_s_a.h +shared/mpp/include/mpp_error_s_s.h +shared/mpp/include/mpp_gather.h +shared/mpp/include/mpp_get_boundary.h +shared/mpp/include/mpp_global_field.h +shared/mpp/include/mpp_global_reduce.h +shared/mpp/include/mpp_global_sum_ad.h +shared/mpp/include/mpp_global_sum.h +shared/mpp/include/mpp_global_sum_tl.h +shared/mpp/include/mpp_io_connect.inc +shared/mpp/include/mpp_io_misc.inc +shared/mpp/include/mpp_io_read.inc +shared/mpp/include/mpp_io_util.inc +shared/mpp/include/mpp_io_write.inc +shared/mpp/include/mpp_read_2Ddecomp.h +shared/mpp/include/mpp_reduce_mpi.h +shared/mpp/include/mpp_reduce_nocomm.h +shared/mpp/include/mpp_reduce_sma.h +shared/mpp/include/mpp_sum.inc +shared/mpp/include/mpp_sum_mpi.h +shared/mpp/include/mpp_sum_nocomm.h +shared/mpp/include/mpp_sum_sma.h +shared/mpp/include/mpp_transmit.inc +shared/mpp/include/mpp_transmit_mpi.h +shared/mpp/include/mpp_transmit_nocomm.h +shared/mpp/include/mpp_transmit_sma.h +shared/mpp/include/mpp_update_domains2D_ad.h +shared/mpp/include/mpp_update_domains2D.h +shared/mpp/include/mpp_update_domains2D_nonblock.h +shared/mpp/include/mpp_update_nest_domains.h +shared/mpp/include/mpp_util.inc +shared/mpp/include/mpp_util_mpi.inc +shared/mpp/include/mpp_util_nocomm.inc +shared/mpp/include/mpp_util_sma.inc +shared/mpp/include/mpp_write_2Ddecomp.h +shared/mpp/include/mpp_write.h +shared/mpp/include/system_clock.h +shared/mpp/mpp_data.F90 +shared/mpp/mpp_domains.F90 +shared/mpp/mpp.F90 +shared/mpp/mpp_io.F90 +shared/mpp/mpp_memutils.F90 +shared/mpp/mpp_parameter.F90 +shared/mpp/mpp_pset.F90 +shared/mpp/mpp_utilities.F90 +shared/mpp/nsclock.c +shared/mpp/test_mpp_domains.F90 +shared/mpp/test_mpp.F90 +shared/mpp/test_mpp_io.F90 +shared/mpp/test_mpp_pset.F90 +shared/mpp/threadloc.c +shared/platform/platform.F90 +shared/time_manager/get_cal_time.F90 +shared/time_manager/time_manager.F90 +shared/tracer_manager/tracer_manager.F90 diff --git a/src/extra/python/isca/__init__.py b/src/extra/python/isca/__init__.py index 5bf28c1a9..28da6f612 100644 --- a/src/extra/python/isca/__init__.py +++ b/src/extra/python/isca/__init__.py @@ -83,4 +83,4 @@ def emit(self, event, *args, **kwargs): from isca.experiment import Experiment, DiagTable, Namelist, FailedRunError -from isca.codebase import IscaCodeBase, SocratesCodeBase, DryCodeBase, GreyCodeBase #, ShallowCodeBase +from isca.codebase import IscaCodeBase, SocratesCodeBase, DryCodeBase, GreyCodeBase, ShallowCodeBase, BarotropicCodeBase diff --git a/src/extra/python/isca/codebase.py b/src/extra/python/isca/codebase.py index 09353e7ce..0a08a9307 100644 --- a/src/extra/python/isca/codebase.py +++ b/src/extra/python/isca/codebase.py @@ -389,8 +389,14 @@ class DryCodeBase(GreyCodeBase): -# class ShallowCodeBase(CodeBase): -# """The Shallow Water Equations. -# """ -# name = 'shallow' -# executable_name = 'shallow.x' +class ShallowCodeBase(CodeBase): + """The Shallow Water Equations. + """ + name = 'shallow' + executable_name = 'shallow.x' + +class BarotropicCodeBase(CodeBase): + """The Barotropic vorticity equations. + """ + name = 'barotropic' + executable_name = 'barotropic_isca.x' \ No newline at end of file diff --git a/src/extra/python/scripts/remove_certain_restart_and_data_files.py b/src/extra/python/scripts/remove_certain_restart_and_data_files.py index 98d658a95..582a05c9c 100644 --- a/src/extra/python/scripts/remove_certain_restart_and_data_files.py +++ b/src/extra/python/scripts/remove_certain_restart_and_data_files.py @@ -1,6 +1,7 @@ import sh import os import pdb +from glob import glob P = os.path.join @@ -12,11 +13,15 @@ def __init__(self, basedir, workdir, datadir, exp_name): self.expname = exp_name -def create_exp_object(exp_name): +def create_exp_object(exp_name, data_directory=None): + + if data_directory is None: + datadir = os.environ['GFDL_DATA'] + else: + datadir = data_directory - basedir = os.environ['GFDL_BASE'] workdir = os.environ['GFDL_WORK'] - datadir = os.environ['GFDL_DATA'] + basedir = os.environ['GFDL_BASE'] expname = '/'+exp_name+'/' exp_object = temporary_exp_object(basedir, workdir, datadir, exp_name) @@ -24,34 +29,42 @@ def create_exp_object(exp_name): return exp_object -def keep_only_certain_restart_files(exp_object, max_num_files, interval=12): +# def keep_only_certain_restart_files(exp_object, max_num_files, interval=12): - # sh.ls(sh.glob(P(self.workdir,'restarts','res_*.cpio'))) #TODO get max_num_files calculated in line, rather than a variable to pass. +# # sh.ls(sh.glob(P(self.workdir,'restarts','res_*.cpio'))) #TODO get max_num_files calculated in line, rather than a variable to pass. - #First defines a list of ALL the restart file numbers - files_to_remove=list(range(0,max_num_files)) +# #First defines a list of ALL the restart file numbers +# files_to_remove=list(range(0,max_num_files)) - #Then defines a list of the ones we want to KEEP - files_to_keep =list(range(0,max_num_files,interval)) +# #Then defines a list of the ones we want to KEEP +# files_to_keep =list(range(0,max_num_files,interval)) - #Then we remove the files we want to keep from the list of all files, giving a list of those we wish to remove - for x in files_to_keep: - files_to_remove.remove(x) +# #Then we remove the files we want to keep from the list of all files, giving a list of those we wish to remove +# for x in files_to_keep: +# files_to_remove.remove(x) - #Then we remove them. - for entry in files_to_remove: - try: - sh.rm(P(exp_object.workdir,exp_object.expname,'restarts','res_'+str(entry)+'.cpio')) -# print P(exp_object.workdir,exp_object.expname,'restarts','res_'+str(entry)+'.cpio') +# #Then we remove them. +# for entry in files_to_remove: +# try: +# sh.rm(P(exp_object.workdir,exp_object.expname,'restarts','res_'+str(entry)+'.cpio')) +# # print P(exp_object.workdir,exp_object.expname,'restarts','res_'+str(entry)+'.cpio') - except sh.ErrorReturnCode_1: - pass -# print 'Tried to remove some restart files, but number '+str(entry)+' does not exist' +# except sh.ErrorReturnCode_1: +# pass +# # print 'Tried to remove some restart files, but number '+str(entry)+' does not exist' -def keep_only_certain_restart_files_data_dir(exp_object, max_num_files, interval=12): +def keep_only_certain_restart_files_data_dir(exp_object, max_num_files=None, interval=12): # sh.ls(sh.glob(P(self.workdir,'restarts','res_*.cpio'))) #TODO get max_num_files calculated in line, rather than a variable to pass. + if max_num_files is None: + month_list = glob(P(exp_object.datadir,exp_object.expname, 'restarts')+'/res*.tar.gz') + if len(month_list)==0: + return + else: + final_month = month_list[-1].split('/res') + max_num_files = int(final_month[-1].split('.tar.gz')[0]) + #First defines a list of ALL the restart file numbers files_to_remove=list(range(0,max_num_files)) @@ -62,13 +75,27 @@ def keep_only_certain_restart_files_data_dir(exp_object, max_num_files, interval for x in files_to_keep: files_to_remove.remove(x) + first_to_be_removed = True + number_removed = 0 + number_not_removed = 0 #Then we remove them. - for entry in files_to_remove: + for entry in files_to_remove[1:-1]: try: - sh.rm(P(exp_object.datadir,exp_object.expname,'run%03d' % entry,'INPUT','res')) -# print 'would be removing ' + P(exp_object.datadir,exp_object.expname,'run'+str(entry),'INPUT','res') + file_to_remove = P(exp_object.datadir,exp_object.expname, 'restarts', 'res%04d.tar.gz' % entry) + if os.path.isfile(file_to_remove) and first_to_be_removed: + first_to_be_removed=False + number_not_removed+=1 + # print('would have removed '+file_to_remove+' but wanted to make sure not to delete the first restart') + else: + sh.rm(file_to_remove) + number_removed+=1 + # print('have removed ' + file_to_remove) except sh.ErrorReturnCode_1: + number_not_removed+=1 + # print('could not remove ' + file_to_remove) pass + print(P(exp_object.datadir,exp_object.expname), 'number removed '+str(number_removed), 'number not removed '+str(number_not_removed)) + # print 'Tried to remove some restart files, but number '+str(entry)+' does not exist' def keep_only_certain_daily_data_uninterp(exp_object, max_num_files, interval=None, file_name = 'atmos_daily.nc'): @@ -88,8 +115,8 @@ def keep_only_certain_daily_data_uninterp(exp_object, max_num_files, interval=No #Then we remove them. for entry in files_to_remove: try: - sh.rm(P(exp_object.datadir,exp_object.expname,'run%03d' % entry,file_name)) - print(('Removed '+P(exp_object.datadir,exp_object.expname,'run%03d' % entry,file_name))) + sh.rm(P(exp_object.datadir,exp_object.expname,'run%04d' % entry,file_name)) + print(('Removed '+P(exp_object.datadir,exp_object.expname,'run%04d' % entry,file_name))) except sh.ErrorReturnCode_1: pass # print 'Tried to remove some atmos_daily files, but number '+str(entry)+' does not exist' @@ -98,36 +125,18 @@ def keep_only_certain_daily_data_uninterp(exp_object, max_num_files, interval=No if __name__=="__main__": - max_num_files_input = 325 + max_num_files_input = None -# exp_name_list=['simple_continents_post_princeton_qflux_anoms_'+str(x) for x in range(31,32)] - -# exp_name_list=['aquaplanet_qflux_anoms_'+str(x) for x in [12,18,23,32,8]] - -# exp_name_list = ['simple_continents_post_princeton_qflux_control_1','simple_continents_post_princeton_fixed_sst_1', 'simple_continents_post_princeton_qflux_control_nod_1', 'simple_continents_post_princeton_qflux_control_scf_1'] -# -# exp_name_list = ['annual_mean_ice_princeton_qflux_control_matrix_qflux_2017_code_1', 'annual_mean_ice_post_princeton_fixed_sst_1', 'annual_mean_ice_princeton_fixed_sst_1'] -# -# exp_name_list.extend(['annual_mean_ice_post_princeton_fixed_sst_el_nino_1']) - -# exp_name_list = ['simple_continents_post_princeton_qflux_control_1'] - -# exp_name_list = ['annual_mean_ice_princeton_qflux_control_1']#, 'annual_mean_ice_post_princeton_qflux_control_1'] - -# exp_name_list = ['annual_mean_ice_post_princeton_fixed_sst_TEST_1', 'annual_mean_ice_princeton_qflux_control_matrix_qflux_1'] - -# exp_name_list.extend(['simple_continents_post_princeton_fixed_sst_1']) - -# exp_name_list = ['giant_drag_exp_chai_values_without_dc_bug_latest_1'] -# exp_name_list = ['aquaplanet_qflux_control_1'] + # exp_name_list = [''] + exp_name_list = glob('/disca/share/sit204/data_isca_from_gv5/frierson_post_soc_fix_*/') - exp_name_list = ['giant_drag_exp_chai_values_with_dc_bug_latest_start_to_finish_1', 'giant_drag_exp_chai_values_without_dc_bug_latest_start_to_finish_1'] for exp_name_input in exp_name_list: - temp_obj = create_exp_object(exp_name_input) - keep_only_certain_restart_files(temp_obj, max_num_files_input) + print('Percentage progress through list:'+str(exp_name_list.index(exp_name_input)/len(exp_name_list))) + temp_obj = create_exp_object(exp_name_input, data_directory='/disca/share/sit204/data_from_isca_cpu/') + # keep_only_certain_restart_files(temp_obj, max_num_files_input) keep_only_certain_restart_files_data_dir(temp_obj, max_num_files_input) - keep_only_certain_daily_data_uninterp(temp_obj, max_num_files_input, file_name = 'fms_moist.x') + # keep_only_certain_daily_data_uninterp(temp_obj, max_num_files_input, file_name = 'fms_moist.x') # keep_only_certain_daily_data_uninterp(temp_obj, max_num_files_input) \ No newline at end of file