+
+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..14b7e89db
--- /dev/null
+++ b/src/atmos_spectral_shallow/shallow_dynamics.F90
@@ -0,0 +1,730 @@
+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
+
+use interpolator_mod, only: interpolate_type,interpolator_init,CONSTANT,interpolator
+
+!======================================================================================
+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.
+
+logical :: initial_condition_from_input_file=.false.
+character(len=64) :: init_cond_file = 'init_cond_h_vor_div'
+character(len=64) :: input_file_div_name = 'div'
+character(len=64) :: input_file_height_name = 'height'
+character(len=64) :: input_file_vor_name = 'vor'
+
+real, dimension(2) :: valid_range_v = (/-1.e3,1.e3/)
+
+type(interpolate_type),save :: init_cond_interp
+
+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, &
+ initial_condition_from_input_file, &
+ init_cond_file, &
+ input_file_div_name, &
+ input_file_height_name, &
+ input_file_vor_name
+
+
+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, allocatable, dimension(:,:) :: rad_lonb_2d, rad_latb_2d
+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.)
+allocate (rad_lonb_2d(is:ie+1, js:je+1))
+allocate (rad_latb_2d(is:ie+1, js:je+1))
+
+do i = is,ie+1
+ rad_lonb_2d(i,:) = glon_bnd(i)
+enddo
+
+do j = js,je+1
+ rad_latb_2d(:,j) = glat_bnd(j)
+enddo
+
+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
+
+if( initial_condition_from_input_file ) then
+ call interpolator_init( init_cond_interp, trim(init_cond_file)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) )
+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
+
+ if (initial_condition_from_input_file) then
+
+ call interpolator( init_cond_interp, Dyn%Grid%div(:,:,1), input_file_div_name )
+ call interpolator( init_cond_interp, Dyn%Grid%h(:,:,1), input_file_height_name )
+ call interpolator( init_cond_interp, Dyn%Grid%vor(:,:,1), input_file_vor_name )
+
+ Dyn%Grid%h(:,:,1) = Dyn%Grid%h(:,:,1)+h_0 !want to make sure that we keep h_0 consistent, so make sure mean of h input is zero, and add h_0 on afterwards...
+
+ else
+ 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 !add_initial_vortex_pair
+ enddo
+
+ endif ! initial_condition_from_input_file
+
+ 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..02dbdac8a
--- /dev/null
+++ b/src/extra/model/shallow/path_names
@@ -0,0 +1,148 @@
+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
+atmos_shared/interpolator/interpolator.F90
+shared/horiz_interp/horiz_interp_bicubic.F90
+shared/horiz_interp/horiz_interp_bilinear.F90
+shared/horiz_interp/horiz_interp_conserve.F90
+shared/horiz_interp/horiz_interp.F90
+shared/horiz_interp/horiz_interp_spherical.F90
+shared/horiz_interp/horiz_interp_type.F90
+shared/time_interp/time_interp_external.F90
+shared/time_interp/time_interp.F90
+shared/axis_utils/axis_utils.F90
diff --git a/src/extra/python/isca/__init__.py b/src/extra/python/isca/__init__.py
index 83600dfc5..830b10774 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, ColumnCodeBase, SocColumnCodeBase #, ShallowCodeBase
+from isca.codebase import IscaCodeBase, SocratesCodeBase, DryCodeBase, GreyCodeBase, ShallowCodeBase, BarotropicCodeBase, ColumnCodeBase, SocColumnCodeBase
diff --git a/src/extra/python/isca/codebase.py b/src/extra/python/isca/codebase.py
index ada901b27..c411ff928 100644
--- a/src/extra/python/isca/codebase.py
+++ b/src/extra/python/isca/codebase.py
@@ -464,8 +464,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/calendar_calc.py b/src/extra/python/scripts/calendar_calc.py
index 2c40f6673..ac3aac37e 100644
--- a/src/extra/python/scripts/calendar_calc.py
+++ b/src/extra/python/scripts/calendar_calc.py
@@ -1,4 +1,4 @@
-from cftime import utime
+import cftime
from datetime import datetime
from cmip_time import FakeDT
import numpy as np
@@ -8,9 +8,11 @@
def day_number_to_datetime_array(time_in, calendar_type, units_in):
- cdftime = utime(units_in, calendar = calendar_type)
+ # cdftime = utime(units_in, calendar = calendar_type)
- date_out = cdftime.num2date(time_in)
+ # date_out = cdftime.num2date(time_in)
+
+ date_out = cftime.num2date(time_in, units_in, calendar=calendar_type)
return date_out
diff --git a/src/extra/python/scripts/create_timeseries.py b/src/extra/python/scripts/create_timeseries.py
index 94a4a3998..b18d241d5 100644
--- a/src/extra/python/scripts/create_timeseries.py
+++ b/src/extra/python/scripts/create_timeseries.py
@@ -107,7 +107,11 @@ def create_time_arr(num_years,is_climatology,time_spacing):
return time_arr,day_number,ntime,time_units,time_bounds
-def output_to_file(data,lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,file_name,variable_name,number_dict, time_bounds=None):
+def output_multiple_variables_to_file(data_dict,lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,file_name,number_dict, time_bounds=None):
+ """Default is to accept multiple data variables to put in the file.
+ Input format in data_dict = {variable_name:variable_array}.
+ Currently only accepts all 2D fields or all 3D fields.
+ Could be updated in future to accept a mix."""
output_file = Dataset(file_name, 'w', format='NETCDF3_CLASSIC')
@@ -116,6 +120,10 @@ def output_to_file(data,lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,
else:
is_thd=True
+ if time_arr is None:
+ add_time = False
+ else:
+ add_time = True
lat = output_file.createDimension('lat', number_dict['nlat'])
lon = output_file.createDimension('lon', number_dict['nlon'])
@@ -137,7 +145,8 @@ def output_to_file(data,lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,
pfull = output_file.createDimension('pfull', number_dict['npfull'])
phalf = output_file.createDimension('phalf', number_dict['nphalf'])
- time = output_file.createDimension('time', 0) #s Key point is to have the length of the time axis 0, or 'unlimited'. This seems necessary to get the code to run properly.
+ if add_time:
+ time = output_file.createDimension('time', 0) #s Key point is to have the length of the time axis 0, or 'unlimited'. This seems necessary to get the code to run properly.
latitudes = output_file.createVariable('lat','d',('lat',))
longitudes = output_file.createVariable('lon','d',('lon',))
@@ -146,7 +155,8 @@ def output_to_file(data,lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,
pfulls = output_file.createVariable('pfull','d',('pfull',))
phalfs = output_file.createVariable('phalf','d',('phalf',))
- times = output_file.createVariable('time','d',('time',))
+ if add_time:
+ times = output_file.createVariable('time','d',('time',))
latitudes.units = 'degrees_N'.encode('utf-8')
latitudes.cartesian_axis = 'Y'
@@ -179,11 +189,11 @@ def output_to_file(data,lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,
phalfs.positive = 'down'
phalfs.long_name = 'half pressure level'
-
- times.units = time_units
- times.calendar = 'THIRTY_DAY_MONTHS'
- times.calendar_type = 'THIRTY_DAY_MONTHS'
- times.cartesian_axis = 'T'
+ if add_time:
+ times.units = time_units
+ times.calendar = 'THIRTY_DAY_MONTHS'
+ times.calendar_type = 'THIRTY_DAY_MONTHS'
+ times.cartesian_axis = 'T'
if time_bounds is not None:
@@ -198,12 +208,6 @@ def output_to_file(data,lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,
times.bounds = 'time_bounds'
-
- if is_thd:
- output_array_netcdf = output_file.createVariable(variable_name,'f4',('time','pfull', 'lat','lon',))
- else:
- output_array_netcdf = output_file.createVariable(variable_name,'f4',('time','lat','lon',))
-
latitudes[:] = lats
longitudes[:] = lons
@@ -216,14 +220,34 @@ def output_to_file(data,lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,
pfulls[:] = p_full
phalfs[:] = p_half
- if type(time_arr[0])!=np.float64 and type(time_arr[0])!=np.int64 :
- times[:] = date2num(time_arr,units='days since 0001-01-01 00:00:00.0',calendar='360_day')
- else:
- times[:] = time_arr
-
- output_array_netcdf[:] = data
+ if add_time:
+ if type(time_arr[0])!=np.float64 and type(time_arr[0])!=np.int64 :
+ times[:] = date2num(time_arr,units='days since 0001-01-01 00:00:00.0',calendar='360_day')
+ else:
+ times[:] = time_arr
+
+ for variable_name in data_dict.keys():
+ if is_thd:
+ if add_time:
+ three_d_dims = ('time','pfull', 'lat','lon',)
+ else:
+ three_d_dims = ('pfull', 'lat','lon',)
+
+ output_array_netcdf = output_file.createVariable(variable_name,'f4',three_d_dims)
+ else:
+ if add_time:
+ two_d_dims = ('time','lat','lon',)
+ else:
+ two_d_dims = ('lat','lon',)
+ output_array_netcdf = output_file.createVariable(variable_name,'f4',two_d_dims)
+
+ output_array_netcdf[:] = data_dict[variable_name]
output_file.close()
+def output_to_file(data,lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,file_name,variable_name,number_dict, time_bounds=None):
+ """Special interface for script wanting to output 1 variable only."""
+ data_dict_to_send = {variable_name:data}
+ output_multiple_variables_to_file(data_dict_to_send,lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,file_name,number_dict,time_bounds=time_bounds)
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
diff --git a/src/extra/python/scripts/shallow_water_init_conds.py b/src/extra/python/scripts/shallow_water_init_conds.py
new file mode 100644
index 000000000..2f8bca068
--- /dev/null
+++ b/src/extra/python/scripts/shallow_water_init_conds.py
@@ -0,0 +1,184 @@
+# -*- coding: utf-8 -*-s
+from typing import NoReturn
+import numpy as np
+import pdb
+import create_timeseries as cts
+import xarray as xar
+import gauss_grid as gg
+import matplotlib.pyplot as plt
+import windspharm as wsp
+import pdb
+
+def convert_to_vor_div(u_in, v_in, lat_arr, planet_radius):
+ """convert spherical polar velocities to vor and div"""
+
+ uwnd, uwnd_info = wsp.tools.prep_data(u_in, 'yx')
+ vwnd, vwnd_info = wsp.tools.prep_data(v_in, 'yx')
+
+ # It is also required that the latitude dimension is north-to-south. Again the
+ # bundled tools make this easy.
+ lat_1d_ordered, uwnd, vwnd = wsp.tools.order_latdim(lat_arr[:,0], uwnd, vwnd)
+
+ # Create a VectorWind instance to handle the computation of streamfunction and
+ # velocity potential.
+ w = wsp.standard.VectorWind(uwnd, vwnd, rsphere=planet_radius, gridtype='gaussian')
+
+ # Compute the streamfunction and velocity potential. Also use the bundled
+ # tools to re-shape the outputs to the 4D shape of the wind components as they
+ # were read off files.
+ vor = w.vorticity()
+ div = w.divergence()
+ # sf, vp = w.sfvp()
+ vor = wsp.tools.recover_data(vor, uwnd_info)
+ div = wsp.tools.recover_data(div, uwnd_info)
+
+ return vor[::-1,:], div[::-1,:] #need to reverse latitude reordering
+
+def set_u_v_height_field(lon_in, lat_in, lonb_in, latb_in, epsilon, alpha, beta, m, r_0, planet_radius, northern_hemisphere=True):
+ """Configure an initial condition for u, v and h given some
+ balance condition. Use parameters and gradient-wind balance for Saturn
+ from 10.1016/j.icarus.2017.06.006"""
+
+ deformation_scale = 3200e3 #p62 of Rostami et al 2017
+ f_0 = 3.2e-4
+ timescale = (f_0)**-1
+ velocity_scale = deformation_scale/timescale
+
+ lat_rad_2d = np.deg2rad(lat_in)
+ lon_rad_2d = np.deg2rad(lon_in)
+
+ if northern_hemisphere:
+ r_array = (planet_radius * (np.pi/2. - lat_rad_2d))/deformation_scale #non-dim
+ else:
+ r_array = (planet_radius * (np.pi/2. + lat_rad_2d))/deformation_scale #non-dim
+
+ v = np.zeros_like(lat_in)
+ u = epsilon * ((r_array - r_0)**alpha)* np.exp(-m*((r_array-r_0)**beta))
+
+ v_si_units = v * velocity_scale
+ u_si_units = u * velocity_scale
+
+ if northern_hemisphere:
+ grad_geopot = ((u_si_units**2)/(r_array* deformation_scale)) + (f_0*np.sin(lat_rad_2d)*u_si_units)
+ else:
+ #I've changed the sign of the coriolis term here. Clearly this isn't really happening, but in this funny radial coordinate system the sign of u would change for the opposite hemisphere, thus necessitating the sign change.
+ grad_geopot = ((u_si_units**2)/(r_array* deformation_scale)) - (f_0*np.sin(lat_rad_2d)*u_si_units)
+
+ geopotential = np.zeros_like(grad_geopot)
+
+ if northern_hemisphere:
+ for lat_idx in range(1, len(lat_rad_2d[:,0])):
+ geopotential[lat_idx,:] = geopotential[lat_idx-1,:] + 0.5*(grad_geopot[lat_idx-1,:]+grad_geopot[lat_idx,:])*(r_array[lat_idx]-r_array[lat_idx-1])
+ else:
+ r_array_opposite = r_array[::-1,:]
+ grad_geopot_opposite = grad_geopot[::-1,:]
+ for lat_idx in range(1, len(lat_rad_2d[:,0])):
+ geopotential[lat_idx,:] = geopotential[lat_idx-1,:] + 0.5*(grad_geopot_opposite[lat_idx-1,:]+grad_geopot_opposite[lat_idx,:])*(r_array_opposite[lat_idx]-r_array_opposite[lat_idx-1])
+
+ geopotential = geopotential[::-1,:]
+
+ #we want to pass a geopotential field that has an area-mean of zero. This is because we want to preserve the mean geopotential that the model sets as its h_0 parameter.
+
+ delta_lat_arr = np.diff(latb_in, axis=0)[:,0:-1]
+
+ area_array = np.cos(np.deg2rad(lat_in))*np.deg2rad(delta_lat_arr)
+
+ area_av_geopot = np.sum(geopotential*area_array)/np.sum(area_array)
+
+ geopotential_av_removed = geopotential-area_av_geopot
+
+ area_av_final = np.sum(geopotential_av_removed*area_array)/np.sum(area_array)
+
+ print(f'old mean = {area_av_geopot}, final area_av geopot = {area_av_final}')
+
+ geopotential_si_units = geopotential_av_removed * deformation_scale
+
+ h_0 = (deformation_scale*f_0)**2.
+
+ return u_si_units, v_si_units, geopotential_si_units, h_0, grad_geopot
+
+nlat=128
+nlon=256
+
+latitudes, latitude_bounds_2 = gg.gaussian_latitudes(int(nlat/2))
+latitude_bounds = [latitude_bound[0] for latitude_bound in latitude_bounds_2] + [latitude_bounds_2[-1][1]]
+
+longitudes = np.linspace(0., 360., nlon, endpoint=False)
+delta_lon = longitudes[1]-longitudes[0]
+longitude_bounds = [lon_val-(0.5*delta_lon) for lon_val in longitudes] + [np.max(longitudes)+(0.5*delta_lon)]
+time_arr_adj=None
+
+lon_array_2d, lat_array_2d = np.meshgrid(longitudes, latitudes)
+lonb_array_2d, latb_array_2d = np.meshgrid(longitude_bounds, latitude_bounds)
+
+#Note that in the following we're making the initial condition symmetric about the equator. This is because if you only set the initial conditions in the northern hemisphere then you end up needing a very large set of latitudinal functions to get that level of asymmetry, and the code gets very upset when translating that to a finite spectral representation. Making it symmetric gets rid of this problem, at least to some extent.
+
+epsilon = 0.15*2.
+alpha = 0.42
+beta = 1.3
+r_0 = 0.
+m_param = 1.
+planet_radius = 55000e3
+
+u_array_vortex, v_array_vortex, height_array_vortex, h_0, grad_geopot_vortex = set_u_v_height_field(lon_array_2d, lat_array_2d,lonb_array_2d, latb_array_2d, epsilon, alpha, beta, m_param, r_0, planet_radius)
+
+u_array_vortex_sp, v_array_vortex_sp, height_array_vortex_sp, h_0_sp, grad_geopot_vortex_sp = set_u_v_height_field(lon_array_2d, lat_array_2d,lonb_array_2d, latb_array_2d, epsilon, alpha, beta, m_param, r_0, planet_radius, northern_hemisphere=False)
+
+epsilon = 0.08
+alpha = 0.
+beta = 2.
+r_0 = 3.37
+m_param = 3.
+planet_radius = 55000e3
+
+u_array_jet, v_array_jet, height_array_jet, h_0, grad_geopot_jet = set_u_v_height_field(lon_array_2d, lat_array_2d,lonb_array_2d, latb_array_2d, epsilon, alpha, beta, m_param, r_0, planet_radius)
+
+u_array_jet_sp, v_array_jet_sp, height_array_jet_sp, h_0_sp, grad_geopot_jet_sp = set_u_v_height_field(lon_array_2d, lat_array_2d,lonb_array_2d, latb_array_2d, epsilon, alpha, beta, m_param, r_0, planet_radius, northern_hemisphere=False)
+
+
+u_array_total = u_array_vortex+u_array_vortex_sp + u_array_jet+u_array_jet_sp
+v_array_total = v_array_vortex+v_array_vortex_sp + v_array_jet+v_array_jet_sp
+height_array_total = height_array_vortex+height_array_vortex_sp + height_array_jet+height_array_jet_sp
+grad_geopot_total = grad_geopot_vortex + grad_geopot_vortex_sp + grad_geopot_jet + grad_geopot_jet_sp
+
+vor_array, div_array = convert_to_vor_div(u_array_total, v_array_total, height_array_total, planet_radius)
+
+
+p_full=None
+p_half=None
+
+npfull=None
+nphalf=None
+
+#Output it to a netcdf file.
+
+file_name='rostami_t85_jet_and_vortex_mk7_gg.nc'
+
+
+number_dict={}
+number_dict['nlat']=nlat
+number_dict['nlon']=nlon
+number_dict['nlatb']=nlat+1
+number_dict['nlonb']=nlon+1
+number_dict['npfull']=npfull
+number_dict['nphalf']=nphalf
+number_dict['ntime']=None
+
+data_dict = {
+ 'vor': vor_array,
+ 'height': height_array_total,
+ 'div': div_array,
+ 'ucomp': u_array_total,
+ 'vcomp': v_array_total,
+ 'grad_geopot': grad_geopot_total
+}
+
+
+time_units=None
+
+cts.output_multiple_variables_to_file(data_dict,latitudes,longitudes,latitude_bounds,longitude_bounds,p_full,p_half,time_arr_adj,time_units,file_name,number_dict)
+
+print(f'Must set h_0 parameter in code to be {h_0}')
+
+
+