!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! GNU General Public License !!
!! !!
!! This file is part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS is free software; you can redistribute it and/or modify !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation. !!
!! !!
!! FMS is distributed in the hope that it will be useful, !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!
!! GNU General Public License for more details. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module polvani_2004_mod
! Polvani, L. M., R. K. Scott, and S. J. Thomas, 2004:
! Numerically Converged Solutions of the Global Primitive Equations
! for Testing the Dynamical Core of Atmospheric GCMs.
! Monthly Weather Review, 132, 2539-2552.
use mpp_mod, only: mpp_error, FATAL
use fms_mod, only: mpp_pe, mpp_root_pe, stdlog, &
write_version_number, close_file, check_nml_error, open_namelist_file
use constants_mod, only: PI, &
GRAV, & ! GRAV=9.800, not 9.806 as specified by Polvani
RDGAS, & ! RDGAS=287.04, not 287.00 as specified by Polvani
OMEGA, & ! same as specified by Polvani
RADIUS,& ! same as specified by Polvani
KAPPA ! same as specified by Polvani
use vert_coordinate_mod, only: compute_vert_coord
use transforms_mod, only: get_grid_boundaries, get_sin_lat, get_cos_lat, trans_grid_to_spherical, trans_spherical_to_grid,&
get_wts_lat, get_lon_max, get_lat_max, get_deg_lon, get_deg_lat, &
vor_div_from_uv_grid, uv_grid_from_vor_div, get_grid_domain
use press_and_geopot_mod, only: press_and_geopot_init, pressure_variables
use diag_manager_mod, only: diag_axis_init, register_static_field, send_data
implicit none
private
character(len=128), parameter :: version = &
'$Id: polvani_2004.F90,v 17.0 2009/07/21 03:00:37 fms Exp $'
character(len=128), parameter :: tagname = &
'$Name: mom4p1_pubrel_dec2009_nnz $'
public :: polvani_2004
integer, parameter :: num_standard_levels = 8
real, parameter, dimension(num_standard_levels ) :: z_standard=(/0., 11.e3, 20.e3, 32.e3, 47.e3, 51.e3, 71.e3, 80.e3/)
real, parameter, dimension(num_standard_levels-1) :: lapse_standard=(/-6.5e-3, 0.0, 1.e-3, 2.8e-3, 0.0, -2.8e-3, -2.0e-3/)
real, dimension(num_standard_levels ) :: T_standard
real :: lambda0, phi0, alpha, beta ! These determine the position of the initial temperature perturbation.
real, parameter :: P00=1.e5 ! Used to compute potential temperature. sea_level_press may be the same, but one should not assume that.
real, parameter :: twopi = 2*PI
real :: xx
integer :: kk
real :: H = 7.340e3 ! meters
real :: z0 = 22.e3 ! meters
real :: delta_z0 = 5.e3 ! meters
real :: z1 = 30.e3 ! meters
real :: u0 = 50. ! m/sec
real :: perturb_amp = 1.0 ! deg K
namelist / polvani_2004_nml / H, z0, delta_z0, z1, u0, perturb_amp
contains
!=========================================================================================================================
subroutine polvani_2004(sea_level_press, triang_trunc, vert_difference_option, pk, bk, &
vors, divs, ts, ln_ps, ug, vg, tg, psg, vorg, divg, surf_geopotential)
real, intent(in) :: sea_level_press
logical, intent(in) :: triang_trunc
character(len=*), intent(in) :: vert_difference_option
real, intent(out), dimension(:) :: pk, bk
complex, intent(out), dimension(:,:,:) :: vors, divs, ts
complex, intent(out), dimension(:,: ) :: ln_ps
real, intent(out), dimension(:,:,:) :: ug, vg, tg
real, intent(out), dimension(:,: ) :: psg
real, intent(out), dimension(:,:,:) :: vorg, divg
real, intent(out), dimension(:,: ) :: surf_geopotential
integer :: unit, ierr, io, i, j, k, ks, lon_max, lat_max, num_levels
integer :: id_lon, id_lat, id_phalf, id_pfull, id_basic_flow, id_term1_eq10, id_basic_temp, id_pot_temp, id_perturbation
integer :: is, ie, js, je
logical :: used
real, dimension(size(ug,1), size(ug,2)) :: ln_psg
real, dimension(size(ug,3) ) :: p_full, ln_p_full, T0
real, dimension(size(ug,3)+1) :: p_half, ln_p_half
real :: F, du_dz, dF_dz, ff1, ff2, dzz1_dz, dzz2_dz, dff1_dz, dff2_dz, dTdy, dTdy_prev, ln_slp
real, dimension(size(ug,3)) :: z, zz1, zz2, gmean_term1
! The allocatable arrays below must be allocated for the global domain
real, allocatable, dimension(:) :: deg_lon, lonb, deg_lat, latb, rad_lon, rad_lat, wts_lat, sin_lat, cos_lat, tan_lat, coriolis
real, allocatable, dimension(:,:) :: basic_flow, term1_eq10, basic_temp, pot_temp
real, allocatable, dimension(:) :: lon_factor, lat_factor
real, allocatable, dimension(:,:) :: perturbation
!------------------------------------------------------------------------------------------------
unit = open_namelist_file()
ierr=1
do while (ierr /= 0)
read(unit, nml=polvani_2004_nml, iostat=io, end=20)
ierr = check_nml_error (io, 'polvani_2004_nml')
enddo
20 call close_file (unit)
call write_version_number(version, tagname)
if(mpp_pe() == mpp_root_pe()) write (stdlog(), nml=polvani_2004_nml)
T_standard(1) = 288.15
do ks=2,num_standard_levels
T_standard(ks) = T_standard(ks-1) + lapse_standard(ks-1)*(z_standard(ks)-z_standard(ks-1))
enddo
call compute_vert_coord ('even_sigma', 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, pk, bk)
surf_geopotential = 0.0
call press_and_geopot_init(pk, bk, .false., vert_difference_option)
call pressure_variables(p_half, ln_p_half, p_full, ln_p_full, sea_level_press)
call get_lon_max(lon_max)
call get_lat_max(lat_max)
num_levels = size(ug,3)
if(num_levels /= 20) then
call mpp_error(FATAL,'Number of vertical levels must be 20 when using polvani_2007 initialization')
endif
allocate(rad_lon(lon_max), deg_lon(lon_max), lon_factor(lon_max))
allocate(rad_lat(lat_max), deg_lat(lat_max), lat_factor(lat_max))
allocate(wts_lat(lat_max), sin_lat(lat_max), cos_lat(lat_max), tan_lat(lat_max))
allocate(coriolis(lat_max), latb(lat_max+1), lonb(lon_max+1))
allocate(basic_flow(lat_max,num_levels), term1_eq10(lat_max,num_levels))
allocate(basic_temp(lat_max,num_levels), pot_temp(lat_max,num_levels))
allocate(perturbation(lon_max,lat_max))
call get_deg_lon(deg_lon)
rad_lon = PI*deg_lon/180.
call get_deg_lat(deg_lat)
rad_lat = PI*deg_lat/180.
call get_wts_lat(wts_lat)
call get_sin_lat(sin_lat)
call get_cos_lat(cos_lat)
tan_lat = sin_lat/cos_lat
call get_grid_boundaries(lonb, latb, global=.true.) ! units = radians
coriolis = 2*OMEGA*sin_lat
ln_slp = log(sea_level_press)
do k=1,num_levels
z(k) = H*(ln_slp - ln_p_full(k))
zz1(k) = (z(k)-z0)/delta_z0
zz2(k) = PI*z(k)/z1
enddo
do k=1,num_levels
if(z(k) > z_standard(num_standard_levels)) then
T0(k) = T_standard(num_standard_levels)
cycle
endif
do ks=2,num_standard_levels
if(z(k) <= z_standard(ks) .and. z(k) >= z_standard(ks-1)) then
T0(k)=((z_standard(ks)-z(k))*T_standard(ks-1) + (z(k)-z_standard(ks-1))*T_standard(ks))/(z_standard(ks)-z_standard(ks-1))
exit
else
cycle
endif
enddo
enddo
do k=1,num_levels
ff1 = 1 - tanh(zz1(k))**3
ff2 = sin(zz2(k))
F = .5*ff1*ff2
do j=1,lat_max
if(sin_lat(j) > 0.0) then
basic_flow(j,k) = u0*F*(sin(PI*sin_lat(j)**2))**3
else
basic_flow(j,k) = 0.0
endif
enddo
enddo
! The k loop below computes the integral of equation 10 in the paper.
! The inner loop over latitude does the actual integration.
! Prior to that the vertical derivative of F (equation 6) must be computed.
do k=1,num_levels
ff1 = 1 - tanh(zz1(k))**3
ff2 = sin(zz2(k))
dzz1_dz = 1/delta_z0
dzz2_dz = PI/z1
dff1_dz = -3 * ( tanh(zz1(k))/cosh(zz1(k)) )**2 * dzz1_dz
dff2_dz = cos(zz2(k)) * dzz2_dz
dF_dz = .5*(ff1*dff2_dz + dff1_dz*ff2)
gmean_term1(k) = 0.0
do j=1,lat_max
if(sin_lat(j) > 0.0) then
du_dz = u0 * (sin(PI*sin_lat(j)**2))**3 * dF_dz
else
du_dz = 0.0
endif
dTdy = -(H/RDGAS) * (RADIUS*coriolis(j) + 2*basic_flow(j,k)*tan_lat(j)) * du_dz
if(j == 1) then
! term1_eq10(j,k) = .5*wts_lat(j)*dTdy/cos_lat(j)
term1_eq10(j,k) = (rad_lat(j)-latb(j))*dTdy/cos_lat(j)
else
! term1_eq10(j,k) = term1_eq10(j-1,k) + .5*(wts_lat(j-1)*dTdy/cos_lat(j-1) + wts_lat(j)*dTdy/cos_lat(j))
term1_eq10(j,k) = term1_eq10(j-1,k) + (latb(j)-rad_lat(j-1))*dTdy_prev + (rad_lat(j)-latb(j))*dTdy
endif
gmean_term1(k) = gmean_term1(k) + .5*wts_lat(j)*term1_eq10(j,k)
dTdy_prev = dTdy
enddo
enddo
! The temperature perturbation is computed below.
lambda0 = 0.0
phi0 = PI/4.
alpha = 1./3.
beta = 1./6.
do i=1,lon_max
xx = rad_lon(i)-lambda0
! Make sure xx falls in the range -PI to +PI.
kk = nint(xx/twopi)
xx = xx - twopi*kk
lon_factor(i) = 1./cosh(xx/alpha)**2
enddo
do j=1,lat_max
lat_factor(j) = 1./cosh((rad_lat(j)-phi0)/beta)**2
enddo
do j=1,lat_max
do i=1,lon_max
perturbation(i,j) = perturb_amp*lon_factor(i)*lat_factor(j)
enddo
enddo
id_phalf = diag_axis_init('phalf',.01*p_half,'hPa','z','approx half pressure level',direction=-1)
id_pfull = diag_axis_init('pfull',.01*p_full,'hPa','z','approx full pressure level',direction=-1,edges=id_phalf)
id_lon = diag_axis_init('lon',deg_lon,'degrees E','x','longitude')
id_lat = diag_axis_init('lat',deg_lat,'degrees N','y','latitude')
id_basic_flow = register_static_field('polvani_2004', 'basic_flow', (/id_lat,id_pfull/),'initial zonal wind', 'm/sec')
id_term1_eq10 = register_static_field('polvani_2004', 'term1_eq10', (/id_lat,id_pfull/),'first term of equation 10', 'K')
id_basic_temp = register_static_field('polvani_2004', 'basic_temp', (/id_lat,id_pfull/),'initial basic_temp', 'K')
id_pot_temp = register_static_field('polvani_2004', 'pot_temp', (/id_lat,id_pfull/),'initial potential basic_temp', 'K')
id_perturbation = register_static_field('polvani_2004', 'perturbation',(/id_lon,id_lat/), 'initial temperature perturbation','K')
if(id_basic_flow > 0) used = send_data(id_basic_flow, basic_flow)
if(id_term1_eq10 > 0) used = send_data(id_term1_eq10, term1_eq10)
do k=1,num_levels
do j=1,lat_max
basic_temp(j,k) = term1_eq10(j,k) - gmean_term1(k) + T0(k)
pot_temp(j,k) = basic_temp(j,k)*(P00/p_full(k))**KAPPA
enddo
enddo
if(id_basic_temp > 0) used = send_data(id_basic_temp, basic_temp)
if(id_pot_temp > 0) used = send_data(id_pot_temp, pot_temp)
if(id_perturbation>0) used = send_data(id_perturbation, perturbation)
call get_grid_domain(is, ie, js, je)
do k=1,num_levels
do i=1,size(ug,1)
ug(i,:,k) = basic_flow(js:je,k)
tg(i,:,k) = basic_temp(js:je,k)
enddo
enddo
do j=1,size(ug,2)
do i=1,size(ug,1)
tg(i,j,:) = tg(i,j,:) + perturbation(is+i-1,js+j-1)
enddo
enddo
vg = 0.0
! Transform grid fields to spherical waves then transform
! back to ensure that spectral and grid representations match.
call trans_grid_to_spherical(tg, ts)
call trans_spherical_to_grid(ts, tg)
call vor_div_from_uv_grid(ug, vg, vors, divs, triang=triang_trunc)
call uv_grid_from_vor_div(vors, divs, ug, vg)
call trans_spherical_to_grid(vors, vorg)
call trans_spherical_to_grid(divs, divg)
psg = sea_level_press
ln_psg = alog(psg)
call trans_grid_to_spherical(ln_psg, ln_ps)
call trans_spherical_to_grid(ln_ps, ln_psg)
psg = exp(ln_psg)
deallocate(deg_lon, deg_lat, rad_lon, rad_lat, wts_lat, sin_lat, cos_lat, tan_lat, coriolis)
deallocate(basic_flow, term1_eq10, basic_temp, pot_temp)
deallocate(lon_factor, lat_factor)
deallocate(perturbation)
return
end subroutine polvani_2004
!================================================================================
end module polvani_2004_mod