module polvani_2007_mod

! Lorenzo Polvani, and J. Esler:
! Transport and mixing of chemical airmasses in idealized baroclinic life cycles
! J. Geophys. Res., 112, D23102, 10.1029/2007JD008555, 2007.

#ifdef INTERNAL_FILE_NML
use mpp_mod, only: input_nml_file
#else
use fms_mod, only: open_namelist_file
#endif

use               mpp_mod, only: mpp_error, FATAL

use               fms_mod, only: mpp_pe, mpp_root_pe, stdlog, stdout, lowercase, &
                                 write_version_number, close_file, check_nml_error

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, grid_domain

use  press_and_geopot_mod, only: press_and_geopot_init, pressure_variables, compute_pressures_and_heights

use      diag_manager_mod, only: diag_axis_init, register_static_field, send_data

use       tracer_type_mod, only: tracer_type

use     field_manager_mod, only: MODEL_ATMOS

use    tracer_manager_mod, only: get_number_tracers, get_tracer_names

implicit none
private

character(len=128), parameter :: version='$Id: polvani_2007.F90,v 19.0 2012/01/06 20:35:37 fms Exp $'
character(len=128), parameter :: tagname='$Name: tikal $'

public :: polvani_2007, polvani_2007_tracer_init, get_polvani_2007_tracers

character(len=12), parameter :: mod_name='polvani_2007'

real,    parameter :: P00=1.e5        ! Used to compute potential temperature.
integer, parameter :: num_iter=10     ! Number of iterations in computation of surface pressure. (Overkill -- only 3 or 4 would do)
integer, parameter :: num_p2007_tr=11 ! Number of tracers associated with polvani_2007 initialization.
                                      ! Additional tracers may be defined elsewhere in the model. There is no conflict. 

integer :: id_LC2_Pstar, id_LC2_Tstar, id_LC2_U1star, id_LC2_U2star, id_LC2_Ustar, id_LC2_lapse00
logical :: used
integer :: i, j, k, lon_max, lat_max, is, ie, js, je, num_levels
integer :: id_LC1_u, id_LC1_pot_t, id_LC1_t, id_perturb, id_Tss, id_Uss, id_LC2_u, id_LC2_pot_t, id_LC2_t
real    :: du_dz, ln_slp, ztmp

real, allocatable, dimension(:,:) :: LC1_u, LC1_t, LC1_pot_t, LC2_u, LC2_t, LC2_pot_t
real, allocatable, dimension(:,:) :: dTdy, Uss, Tss, perturbation
real, allocatable, dimension(:)   :: sin_lat, tan_lat, coriolis, latb, rad_lat, fy1, fy2
real, allocatable, dimension(:)   :: LC2_Tstar, LC2_U1star, LC2_U2star, LC2_Ustar, LC2_Zstar
real, allocatable, dimension(:)   :: deg_lon, lonb, deg_lat, rad_lon, wts_lat, cos_lat
real, allocatable, dimension(:)   :: fz, dfz_dz, lon_factor, lat_factor, z
real, allocatable, dimension(:)   :: LC1_psurf, LC2_psurf
real, allocatable, dimension(:)   :: LC1_p_full, LC1_ln_p_full, LC1_p_half, LC1_ln_p_half
real, allocatable, dimension(:,:) :: LC2_p_full, LC2_ln_p_full, LC2_p_half, LC2_ln_p_half
real, allocatable, dimension(:)   :: trop_strat_boundary

real, allocatable, dimension(:,:,:) :: p2007_tr ! The tracers associated with polvani_2007 initialization.

character(len=15), parameter, dimension(num_p2007_tr) :: p2007_tr_name =   &
(/'polvani_2007_ov','polvani_2007_s1','polvani_2007_s2','polvani_2007_s3', &
  'polvani_2007_t1','polvani_2007_t2','polvani_2007_t3','polvani_2007_un', &
  'polvani_2007_bp','polvani_2007_b0','polvani_2007_bm'/)

character(len=31), parameter, dimension(num_p2007_tr) :: p2007_tr_longname =   &
(/'overworld tracer               ','high middleworld strat tracer  ','center middleworld strat tracer', &
  'low middleworld strat tracer   ','high middleworld trop tracer   ','center middleworld trop tracer ', &
  'low middleworld trop tracer    ','underworld tracer              ', 'high boundary layer tracer     ',&
  'canonical boundary layer tracer','low boundary layer tracer      '/)

! press_diag is the pressure levels used for diagnostic meta data.
! In the LC2 case they are not the actual model pressure levels.
real, allocatable, dimension(:) :: press_diag

character(len=8) :: type_of_init = 'LC1     ' ! 'LC1' or 'LC2'
real :: T_hat = 1.0      ! deg K ------------- Amplitude of temperature perturbation
integer :: m = 6         ! dimensionless ----- Zonal wavenumber of temperature perturbation
real :: theta_hat = 45.  ! degrees latitude -- Latitude center of temperature perturbation
real :: H  = 7.5e3       ! meters
real :: U0 = 45.         ! m/sec
real :: sigma_top = .02  ! dimensionless
real :: zt = 13.e3       ! meters
real :: lapse = -6.5e-3  ! deg K/meter
real :: T0 = 300.        ! deg K
real :: alpha = 10.      ! dimensionless
real :: Us = 45.         ! m/sec
real :: zs = 1.e4        ! meters
real :: theta_s = 35.    ! degrees latitude -- Center of the second term of the zonal wind profile for LC2 ("us" in equation (8) of the paper)
real :: delta_s = 20.    ! degrees latitude -- Width of same

namelist / polvani_2007_nml / type_of_init, T_hat, m, theta_hat, H, U0, &
                              sigma_top, zt, lapse, T0, alpha, Us, zs, theta_s, delta_s

contains

!=================================================================================================================================

subroutine polvani_2007(triang_trunc, vert_difference_option, pk, bk, &
                   vors, divs, ts, ln_ps, ug, vg, tg, psg, vorg, divg, surf_geopotential)

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

real, dimension(size(ug,1), size(ug,2)) :: ln_psg
integer :: id_lon, id_lat, id_press, ntmp, id_iter, unit, ierr, io
real :: ln_sigma_top, riter(num_iter)

#ifdef INTERNAL_FILE_NML
    read (input_nml_file, nml=polvani_2007_nml, iostat=io)
    ierr = check_nml_error(io, 'polvani_2007_nml')
#else
    unit = open_namelist_file()
    ierr=1
    do while (ierr /= 0)
      read(unit, nml=polvani_2007_nml, iostat=io, end=20)
      ierr = check_nml_error (io, 'polvani_2007_nml')
    enddo
20  call close_file (unit)
#endif
call write_version_number(version, tagname)
if(mpp_pe() == mpp_root_pe()) write (stdlog(), nml=polvani_2007_nml)

call get_lon_max(lon_max)
call get_lat_max(lat_max)
num_levels = size(ug,3)
ntmp = nint(.46875*lat_max)
if(num_levels /= ntmp) then
  call mpp_error(FATAL,'Number of vertical levels must be ',ntmp, &
                ' when using polvani_2007 initialization with ',lat_max,' latitude rows')
endif
pk = 0.0
bk(1) = sigma_top
ln_sigma_top = alog(bk(1))
do k=2,num_levels
  bk(k) = exp((1.0 - (k-1)/float(num_levels))*ln_sigma_top)
enddo
bk(num_levels+1) = 1.0

surf_geopotential = 0.0

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), trop_strat_boundary(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), z(num_levels+1))
allocate(LC1_u(lat_max,num_levels+1), dTdy(lat_max,num_levels+1), Uss(lat_max,num_levels+1))
allocate(LC1_pot_t(lat_max,num_levels+1), LC1_t(lat_max,num_levels+1), Tss(lat_max,num_levels+1))
allocate(LC2_u(lat_max,num_levels+1), LC2_pot_t(lat_max,num_levels+1), LC2_t(lat_max,num_levels+1))
allocate(perturbation(lon_max,lat_max), press_diag(num_levels+1))
allocate(fz(num_levels+1), dfz_dz(num_levels+1), fy1(lat_max), fy2(lat_max))
allocate(LC2_Tstar(lat_max), LC2_U1star(lat_max), LC2_U2star(lat_max), LC2_Ustar(lat_max), LC2_Zstar(lat_max))
allocate(LC1_psurf(lat_max), LC2_psurf(lat_max))
allocate(LC1_p_full(num_levels),   LC1_ln_p_full(num_levels))
allocate(LC1_p_half(num_levels+1), LC1_ln_p_half(num_levels+1))
allocate(LC2_p_full(lat_max,num_levels),   LC2_ln_p_full(lat_max,num_levels))
allocate(LC2_p_half(lat_max,num_levels+1), LC2_ln_p_half(lat_max,num_levels+1))
allocate(p2007_tr(lat_max,num_levels,num_p2007_tr))

call press_and_geopot_init(pk, bk, .false., vert_difference_option)
call pressure_variables(LC1_p_half, LC1_ln_p_half, LC1_p_full, LC1_ln_p_full, P00)
call get_deg_lon(deg_lon)
rad_lon = PI*deg_lon/180.
call get_deg_lat(deg_lat)
rad_lat = PI*deg_lat/180.
trop_strat_boundary = 1500*(55-deg_lat) ! units=meters (only used in the narrow latitude range of about 45 to 52 degrees north)
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(P00)

riter = (/(i,i=1,num_iter)/)
press_diag(1:num_levels) = .01*LC1_p_full
press_diag(num_levels+1) = .01*LC1_p_half(num_levels+1)
id_press = diag_axis_init('p_full',press_diag,'hPa','z','pressure level',direction=-1)
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_iter  = diag_axis_init('iter',riter,'hPa','z','iteration')
id_LC1_u     = register_static_field(mod_name,'LC1_zonal_wind', (/id_lat,id_press/),'initial LC1 zonal wind', 'm/sec')
id_LC1_t     = register_static_field(mod_name,'LC1_temperature',(/id_lat,id_press/),'initial LC1 temperature', 'K')
id_LC1_pot_t = register_static_field(mod_name,'LC1_pot_temp',   (/id_lat,id_press/),'initial LC1 potential temperature','K')
id_Tss       = register_static_field(mod_name,'Tss',            (/id_lat,id_press/),'Tss temperature', 'K')
id_Uss       = register_static_field(mod_name,'Uss',            (/id_lat,id_press/),'Uss zonal wind', 'm/sec')
id_LC2_u     = register_static_field(mod_name,'LC2_zonal_wind', (/id_lat,id_press/),'initial LC2 zonal wind', 'm/sec')
id_LC2_t     = register_static_field(mod_name,'LC2_temperature',(/id_lat,id_press/),'initial LC2 temperature', 'K')
id_LC2_pot_t = register_static_field(mod_name,'LC2_pot_temp',   (/id_lat,id_press/),'initial LC2 potential temperature','K')
id_perturb   = register_static_field(mod_name,'perturbation',   (/id_lon,id_lat/),  'initial temperature perturbation','K')
id_LC2_Pstar = register_static_field(mod_name,'LC2_Pstar',      (/id_lat/),         'initial surface pressure','hPa')
id_LC2_Tstar = register_static_field(mod_name,'LC2_Tstar',      (/id_lat/),         'initial surface temperature','K')
id_LC2_U1star= register_static_field(mod_name,'LC2_U1star',     (/id_lat/),         'initial surface LC1 wind', 'm/sec')
id_LC2_U2star= register_static_field(mod_name,'LC2_U2star',     (/id_lat/),         'initial surface LC2 wind', 'm/sec')
id_LC2_Ustar = register_static_field(mod_name,'LC2_Ustar',      (/id_lat/),         'initial surface wind',     'm/sec')
id_LC2_lapse00=register_static_field(mod_name,'LC2_lapse00',    (/id_lat/),         'initial surface lapse rate','K/m')

call compute_LC1
call compute_LC2
call compute_perturbation

do j=1,lat_max
  do k=1,num_levels+1
    if(k == num_levels+1) then
      LC1_pot_t(j,k) = LC1_t(j,k)
    else
      LC1_pot_t(j,k) = LC1_t(j,k)*(P00/LC1_p_full(k))**KAPPA
    endif
  enddo
enddo
if(id_LC1_u     > 0) used = send_data(id_LC1_u, LC1_u)
if(id_LC1_t     > 0) used = send_data(id_LC1_t, LC1_t)
if(id_LC1_pot_t > 0) used = send_data(id_LC1_pot_t, LC1_pot_t)
if(id_perturb   > 0) used = send_data(id_perturb,  perturbation)

call get_grid_domain(is, ie, js, je)

if(trim(type_of_init) == 'LC1') then
  do k=1,num_levels
    do i=1,size(ug,1)
      ug(i,:,k) = LC1_u(js:je,k)
      tg(i,:,k) = LC1_t(js:je,k)
    enddo
  enddo
  do j=1,size(ug,2)
    psg(:,j) = LC1_psurf(js+j-1)
  enddo
  call compute_p2007_tr(LC1_t,LC1_psurf)
else if(trim(type_of_init) == 'LC2') then
  do k=1,num_levels
    do i=1,size(ug,1)
      ug(i,:,k) = LC2_u(js:je,k)
      tg(i,:,k) = LC2_t(js:je,k)
    enddo
  enddo
  do j=1,size(ug,2)
    psg(:,j) = LC2_psurf(js+j-1)
  enddo
  call compute_p2007_tr(LC2_t,LC2_psurf)
else
  call mpp_error(FATAL,trim(type_of_init)//' is not a valid type_of_init for polvani_2007')
endif

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)
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)

return
end subroutine polvani_2007
!=================================================================================================================================
subroutine compute_LC1

real :: Tr

do k=1,num_levels+1
  if(k == num_levels+1) then
    z(k) = 0.0
  else
    z(k) = H*(ln_slp - LC1_ln_p_full(k))
  endif
  ztmp = z(k)/zt
  fz(k) = ztmp*exp(-0.5*(ztmp**2 - 1))
  dfz_dz(k) = ((1 - ztmp**2)/zt)*exp(-0.5*(ztmp**2 - 1))
enddo

do j=1,lat_max
  if(sin_lat(j) > 0.0) then
    fy1(j) = (sin(PI*sin_lat(j)**2))**3
  else
    fy1(j) = 0.0
  endif
enddo

do k=1,num_levels+1
  do j=1,lat_max
    LC1_u(j,k) = U0*fy1(j)*fz(k)
    du_dz   = U0*fy1(j)*dfz_dz(k)
    dTdy(j,k) = -(H/RDGAS) * (RADIUS*coriolis(j) + 2*LC1_u(j,k)*tan_lat(j)) * du_dz
  enddo
enddo

! The k loop below computes the integral of equation A4 in the paper.
! The inner loop over latitude does the actual integration.
do k=1,num_levels+1
  do j=1,lat_max
    if(j == 1) then
      LC1_t(j,k) = dTdy(j,k)*(rad_lat(j)-latb(j))
    else
      LC1_t(j,k) = LC1_t(j-1,k) + dTdy(j-1,k)*(latb(j)-rad_lat(j-1)) + dTdy(j,k)*(rad_lat(j)-latb(j))
    endif
  enddo
enddo
do k=1,num_levels+1
  if(k == num_levels+1) then
    Tr = T0
  else
    Tr = T0 + lapse/(zt**-alpha + z(k)**-alpha)**(1./alpha)
  endif
  LC1_t(:,k) = Tr + LC1_t(:,k)
enddo

LC1_psurf = P00

end subroutine compute_LC1
!=================================================================================================================================
subroutine compute_perturbation
    
do i=1,lon_max
  lon_factor(i) = cos(m*rad_lon(i))
enddo
ztmp = PI*theta_hat/180
do j=1,lat_max
  lat_factor(j) = 1.0/cosh(m*(rad_lat(j)-ztmp))**2
enddo
do j=1,lat_max 
  do i=1,lon_max
    perturbation(i,j) = T_hat*lon_factor(i)*lat_factor(j)
  enddo
enddo

end subroutine compute_perturbation
!=================================================================================================================================
subroutine compute_LC2

do k=1,num_levels+1
  fz(k) = exp(-z(k)/zs)
  dfz_dz(k) = -fz(k)/zs
enddo

do j=1,lat_max
  fy2(j) = (sin(2*rad_lat(j))**2)*((deg_lat(j)-theta_s)/delta_s)*exp(-((deg_lat(j)-theta_s)/delta_s)**2)
enddo

do k=1,num_levels+1
  do j=1,lat_max
    Uss(j,k) = -Us*fy2(j)*fz(k)
    du_dz    = -Us*fy2(j)*dfz_dz(k)
    dTdy(j,k) = -(H/RDGAS) * (RADIUS*coriolis(j) + 2*Uss(j,k)*tan_lat(j)) * du_dz
  enddo
enddo

! The k loop below computes the integral of equation A4 in the paper.
! The inner loop over latitude does the actual integration.
do k=1,num_levels+1
  do j=1,lat_max
    if(j == 1) then
      Tss(j,k) = dTdy(j,k)*(rad_lat(j)-latb(j))
    else
      Tss(j,k) = Tss(j-1,k) + dTdy(j-1,k)*(latb(j)-rad_lat(j-1)) + dTdy(j,k)*(rad_lat(j)-latb(j))
    endif
  enddo
enddo

if(id_Uss > 0) used = send_data(id_Uss, Uss)
if(id_Tss > 0) used = send_data(id_Tss, Tss)

LC2_u = LC1_u + Uss
LC2_t = LC1_t + Tss

! Surface pressure for LC2 is not constant.
! The paper does not desribe the details of how it is computed, so I concocted my own scheme to do it.
! I have isolated this calculation by putting it in a separate subroutine.
call compute_surf_press

call pressure_variables(LC2_p_half, LC2_ln_p_half, LC2_p_full, LC2_ln_p_full, LC2_psurf)

do j=1,lat_max
  do k=1,num_levels+1
    if(k == num_levels+1) then
      LC2_pot_t(j,k) = LC2_t(j,k)
    else
      LC2_pot_t(j,k) = LC2_t(j,k)*(P00/LC2_p_full(j,k))**KAPPA
    endif
  enddo
enddo
if(id_LC2_u     > 0) used = send_data(id_LC2_u, LC2_u)
if(id_LC2_t     > 0) used = send_data(id_LC2_t, LC2_t)
if(id_LC2_pot_t > 0) used = send_data(id_LC2_pot_t, LC2_pot_t)

return
end subroutine compute_LC2
!=================================================================================================================================
subroutine compute_surf_press

real, allocatable, dimension(:,:) :: Zstar     ! Zstar = H*alog(P00/psurf)
real, allocatable, dimension(:) :: af          ! af = RADIUS*coriolis
real, allocatable, dimension(:) :: dzstar_dy   ! derivative of Zstar with respect to latitude
real, allocatable, dimension(:) :: LC2_lapse00 ! lapse rate at pressure level P00
real, allocatable, dimension(:) :: dlapse00_dy ! derivative of LC2_lapse00 with respect to latitude
real :: rms, c1, c2, e, sqrt_e
integer :: j, itr

allocate(Zstar(lat_max,0:num_iter), af(lat_max), dzstar_dy(lat_max), dlapse00_dy(lat_max), LC2_lapse00(lat_max))
af = RADIUS*coriolis
e = exp(1.0)
sqrt_e = sqrt(e)
c1 = 2*e*(U0/zt)**2
c2 = Us/zs**2

do j=1,lat_max
  if(sin_lat(j) > 0.0) then
    dlapse00_dy(j) = c1*tan_lat(j)*fy1(j)**2 - (af(j) - 2*Us*fy2(j)*tan_lat(j))*c2*fy2(j)
  else
    dlapse00_dy(j) = 0.0
  endif
enddo
dlapse00_dy = -(H/RDGAS)*dlapse00_dy
do j=1,lat_max
  if(j == 1) then
    LC2_lapse00(j) = dlapse00_dy(j)*(rad_lat(j)-latb(j))
  else
    LC2_lapse00(j) = LC2_lapse00(j-1) + dlapse00_dy(j-1)*(latb(j)-rad_lat(j-1)) + dlapse00_dy(j)*(rad_lat(j)-latb(j))
  endif
enddo
LC2_lapse00 = LC2_lapse00 + lapse

! Zstar is computed by an iterative method. The zonal momentum equation for the
! surface winds is integrated from pole to pole each interation. Temperature and
! zonal wind are expressed as a Taylor series expansion about z=0, i.e., Pressure = P00.
Zstar(:,0) = 0.0
do itr=1,num_iter
  do j=1,lat_max
    LC2_Tstar(j)  = LC2_t(j,num_levels+1) + LC2_lapse00(j)*Zstar(j,itr-1)
    LC2_U1star(j) = U0*sqrt_e*fy1(j)*Zstar(j,itr-1)/zt
    LC2_U2star(j) = (Zstar(j,itr-1)/zs - 1)*Us*fy2(j)
    LC2_Ustar(j)  = LC2_U1star(j) + LC2_U2star(j)
    if(sin_lat(j) > 0.0) then
      dzstar_dy(j) = H*LC2_Ustar(j)*(af(j) + LC2_Ustar(j)*tan_lat(j)) / (RDGAS*LC2_Tstar(j))
    else
      dzstar_dy(j) = 0.0
    endif
  enddo
  do j=1,lat_max
    if(j == 1) then
      Zstar(j,itr) = dzstar_dy(j)*(rad_lat(j)-latb(j))
    else
      Zstar(j,itr) = Zstar(j-1,itr) + dzstar_dy(j-1)*(latb(j)-rad_lat(j-1)) + dzstar_dy(j)*(rad_lat(j)-latb(j))
    endif
  enddo
  rms = 0.0
  do j=1,lat_max
    rms = rms + (Zstar(j,itr) - Zstar(j,itr-1))**2
  enddo
  rms = sqrt(rms/lat_max)
  write(stdout(),'("Convergence of surface pressure for polvani_2007 LC2: iteration=",i2,"  rms=",1pe16.8)') itr,rms
enddo

LC2_Zstar = Zstar(:,num_iter)

do j=1,lat_max
  LC2_psurf(j) = P00*exp(-Zstar(j,num_iter)/H)
enddo

if(id_LC2_Pstar  > 0) used = send_data(id_LC2_Pstar,  LC2_psurf)
if(id_LC2_Tstar  > 0) used = send_data(id_LC2_Tstar,  LC2_Tstar )
if(id_LC2_U1star > 0) used = send_data(id_LC2_U1star, LC2_U1star)
if(id_LC2_U2star > 0) used = send_data(id_LC2_U2star, LC2_U2star)
if(id_LC2_Ustar  > 0) used = send_data(id_LC2_Ustar,  LC2_Ustar )
if(id_LC2_lapse00> 0) used = send_data(id_LC2_lapse00,LC2_lapse00)
deallocate(Zstar, af, dzstar_dy, dlapse00_dy, LC2_lapse00)

end subroutine compute_surf_press
!=================================================================================================================================
subroutine polvani_2007_tracer_init(tg, psg, surface_geopotential, tracers)

real, dimension(:,:,:),     intent(in)  :: tg
real, dimension(:,:),       intent(in)  :: psg, surface_geopotential
real, dimension(:,:,:,:,:), intent(out) :: tracers

real,    dimension(is:ie,js:je,num_levels)   :: z_full, p_full, pot_t
real,    dimension(is:ie,js:je,num_levels+1) :: z_half, p_half
integer, dimension(is:ie,js:je,num_levels)   :: world ! overworld=0, high middleworld=1, center middleworld=2, low middleworld=3, underworld=4
logical, dimension(is:ie,js:je,num_levels)   :: strat ! .true. if above tropopause, .false. if below. 

integer :: ntr, num_tracers, ntr_o, ntr_s1, ntr_s2, ntr_s3, ntr_t1, ntr_t2, ntr_t3, ntr_u
character(len=128) :: tname, longname, units

call compute_pressures_and_heights(tg, psg, surface_geopotential, z_full, z_half, p_full, p_half)
do k=1,num_levels
  do j=js,je
    do i=is,je
      pot_t(i,j,k) = tg(i,j,k)*(P00/p_full(i,j,k))**KAPPA
      if(pot_t(i,j,k) >= 380.) then
        world(i,j,k) = 0
      else if(pot_t(i,j,k) >= 350. .and. pot_t(i,j,k) < 380.) then
        world(i,j,k) = 1
      else if(pot_t(i,j,k) > 320. .and. pot_t(i,j,k) < 350.) then
        world(i,j,k) = 2
      else if(pot_t(i,j,k) > 290. .and. pot_t(i,j,k) <= 320.) then
        world(i,j,k) = 3
      else
        world(i,j,k) = 4
      endif
      if(z_full(i,j,k) > trop_strat_boundary(j)) then
        strat(i,j,k) = .true.
      else
        strat(i,j,k) = .false.
      endif
    enddo
  enddo
enddo

ntr_o=0; ntr_s1=0; ntr_s2=0; ntr_s3=0; ntr_t1=0; ntr_t2=0; ntr_t3=0; ntr_u=0
call get_number_tracers(MODEL_ATMOS, num_prog=num_tracers)
do ntr=1,num_tracers
  call get_tracer_names(MODEL_ATMOS, ntr, tname, longname, units)
  if(trim(lowercase(tname)) == 'polvani_2007_o') then
    ntr_o = ntr
    tracers(:,:,:,:,ntr) = 0.0
  else if(trim(lowercase(tname)) == 'polvani_2007_s1') then
    ntr_s1 = ntr
    tracers(:,:,:,:,ntr) = 0.0
  else if(trim(lowercase(tname)) == 'polvani_2007_s2') then
    ntr_s2 = ntr
    tracers(:,:,:,:,ntr) = 0.0
  else if(trim(lowercase(tname)) == 'polvani_2007_s3') then
    ntr_s3 = ntr
    tracers(:,:,:,:,ntr) = 0.0
  else if(trim(lowercase(tname)) == 'polvani_2007_t1') then
    ntr_t1 = ntr
    tracers(:,:,:,:,ntr) = 0.0
  else if(trim(lowercase(tname)) == 'polvani_2007_t2') then
    ntr_t2 = ntr
    tracers(:,:,:,:,ntr) = 0.0
  else if(trim(lowercase(tname)) == 'polvani_2007_t3') then
    ntr_t3 = ntr
    tracers(:,:,:,:,ntr) = 0.0
  else if(trim(lowercase(tname)) == 'polvani_2007_u') then
    ntr_u = ntr
    tracers(:,:,:,:,ntr) = 0.0
  endif
enddo

do k=1,num_levels
  do j=js,je
    do i=is,je
      if(world(i,j,k) == 0 .and. ntr_o > 0) then
        tracers(i,j,k,:,ntr_o) = 1.0
      else if(world(i,j,k) == 1) then
        if(strat(i,j,k)) then
          if(ntr_s1 > 0) tracers(i,j,k,:,ntr_s1) = 1.0
        else
          if(ntr_t1 > 0) tracers(i,j,k,:,ntr_t1) = 1.0
        endif
      else if(world(i,j,k) == 2) then
        if(strat(i,j,k)) then
          if(ntr_s2 > 0) tracers(i,j,k,:,ntr_s2) = 1.0
        else
          if(ntr_t2 > 0) tracers(i,j,k,:,ntr_t2) = 1.0
        endif
      else if(world(i,j,k) == 3) then
        if(strat(i,j,k)) then
          if(ntr_s3 > 0) tracers(i,j,k,:,ntr_s3) = 1.0
        else
          if(ntr_t3 > 0) tracers(i,j,k,:,ntr_t3) = 1.0
        endif
      else if(world(i,j,k) == 4 .and. ntr_u > 0) then
        tracers(i,j,k,:,ntr_u) = 1.0
      endif
    enddo
  enddo
enddo

end subroutine polvani_2007_tracer_init
!=================================================================================================================================

subroutine compute_p2007_tr(temp, surfp)
real, dimension(lat_max,num_levels), intent(in) :: temp
real, dimension(lat_max),            intent(in) :: surfp

real, dimension(lat_max)              :: surf_geopotential
real, dimension(lat_max,num_levels  ) :: pot_t, z_full, p_full
real, dimension(lat_max,num_levels+1) :: z_half, p_half
integer :: world(lat_max,num_levels) ! overworld=0, high middleworld=1, center middleworld=2, low middleworld=3, underworld=4
logical :: strat(lat_max,num_levels) ! .true. if above tropopause, .false. if below
character(len=15) :: tname
integer :: ntr, id_tr, id_latitude, id_pfull, axes(2)

surf_geopotential = 0.0
call compute_pressures_and_heights(temp, surfp, surf_geopotential, z_full, z_half, p_full, p_half)
do k=1,num_levels
  do j=1,lat_max
    pot_t(j,k) = temp(j,k)*(P00/p_full(j,k))**KAPPA
    if(pot_t(j,k) >= 380.) then
      world(j,k) = 0
    else if(pot_t(j,k) >= 350. .and. pot_t(j,k) < 380.) then
      world(j,k) = 1
    else if(pot_t(j,k) > 320. .and. pot_t(j,k) < 350.) then
      world(j,k) = 2
    else if(pot_t(j,k) > 290. .and. pot_t(j,k) <= 320.) then
      world(j,k) = 3
    else
      world(j,k) = 4
    endif
    if(z_full(j,k) > trop_strat_boundary(j)) then
      strat(j,k) = .true.
    else
      strat(j,k) = .false.
    endif
  enddo
enddo

p2007_tr = 0.0
do ntr=1,num_p2007_tr
  tname = p2007_tr_name(ntr)
  do k=1,num_levels
    do j=1,lat_max

      if(world(j,k) == 0 .and. tname == 'polvani_2007_ov') then
        p2007_tr(j,k,ntr) = 1.0
      else if(world(j,k) == 1) then
        if(strat(j,k)) then
          if(tname == 'polvani_2007_s1') p2007_tr(j,k,ntr) = 1.0
        else
          if(tname == 'polvani_2007_t1') p2007_tr(j,k,ntr) = 1.0
        endif
      else if(world(j,k) == 2) then
        if(strat(j,k)) then
          if(tname == 'polvani_2007_s2') p2007_tr(j,k,ntr) = 1.0
        else
          if(tname == 'polvani_2007_t2') p2007_tr(j,k,ntr) = 1.0
        endif
      else if(world(j,k) == 3) then
        if(strat(j,k)) then
          if(tname == 'polvani_2007_s3') p2007_tr(j,k,ntr) = 1.0
        else
          if(tname == 'polvani_2007_t3') p2007_tr(j,k,ntr) = 1.0
        endif
      else if(world(j,k) == 4 .and. tname == 'polvani_2007_un') then
        p2007_tr(j,k,ntr) = 1.0
      endif

      if(z_full(j,k) < 2500. .and. tname == 'polvani_2007_bp') then
        p2007_tr(j,k,ntr) = 1.0
      endif
      if(z_full(j,k) < 1500. .and. tname == 'polvani_2007_b0') then
        p2007_tr(j,k,ntr) = 1.0
      endif
      if(z_full(j,k) < 500.  .and. tname == 'polvani_2007_bm') then
        p2007_tr(j,k,ntr) = 1.0
      endif

    enddo
  enddo
enddo

id_latitude  = diag_axis_init('latitude' , deg_lat, 'degrees_N', 'y', 'latitude' , set_name=mod_name)
id_pfull     = diag_axis_init('pfull', .01*LC1_p_full, 'hPa'    , 'z', 'pressure' , set_name=mod_name)
axes = (/id_latitude,id_pfull/)
do ntr=1,num_p2007_tr
  id_tr = register_static_field(mod_name, 'initial_'//p2007_tr_name(ntr), axes, p2007_tr_longname(ntr), 'none')
  used = send_data(id_tr, p2007_tr(:,:,ntr))
enddo

end subroutine compute_p2007_tr
!=================================================================================================================================
subroutine get_polvani_2007_tracers(tracers)
real, dimension(:,:,:,:), intent(out) :: tracers
integer :: num_tracers, ntr1, ntr2
character(len=128) :: tname, longname, units

call get_number_tracers(MODEL_ATMOS, num_prog=num_tracers)
do ntr1=1,num_tracers
  call get_tracer_names(MODEL_ATMOS, ntr1, tname, longname, units)
  do ntr2=1,num_p2007_tr
    if(trim(lowercase(tname)) == p2007_tr_name(ntr2)) then
      do k=1,num_levels
        do i=1,size(tracers,1)
          tracers(i,:,k,ntr1) = p2007_tr(js:je,k,ntr2)
        enddo
      enddo
    endif
  enddo
enddo

end subroutine get_polvani_2007_tracers
!=================================================================================================================================
end module polvani_2007_mod
