module mod_strong
!$$$ module documentation block
! . . . .
! module: mod_strong
! prgmmr: parrish org: np23 date: 2007-02-15
!
! abstract: high level module for carrying all parameters used for various
! flavors of strong constraint.
! Implement implicit normal mod initialization routines for
! use with analysis increment and analysis increment tendencies.
! reference: Temperton, C., 1989: "Implicit Normal Mode Initialization
! for Spectral Models", . MWR, 117, 436-451.
!
! program history log:
! 2007-02-15 parrish
! 2012-02-08 kleist - add option tlnmc_option to control how TLNMC is applied
! 2013-07-02 parrish - change tlnmc_type to reg_tlnmc_type. tlnmc_type no
! longer used for global application of tlnmc.
! 2014-12-03 derber - remove unused variables
!
! Subroutines Included:
! sub init_strongvars - set default namelist variable values
! sub gproj - project input u,v,mass variable to gravity modes for
! update
! sub gproj_diag - project input u,v,mass variable to gravity modes plus
! diagnostics
! sub gproj_diag_update- project input u,v,mass variable to gravity modes plus
! diagnostic and update
! sub gproj0 -
! sub gproj_ad -
! sub dinmi - obtain balance increment from input tendencies
! sub dinmi_ad - adjoint of dinmi
! sub dinmi0 - lower level--balance increment from input tendencies
! sub balm_1 - compute balance diagnostic variable
! sub getbcf - compute matrices B,C,F as defined in above reference
! sub scale_vars - scale variables as defined in reference
! sub scale_vars_ad - adjoint of scale variables
! sub unscale_vars - unscale variables
! sub unscale_vars_ad - adjoint of unscale variables
! sub f_mult - multiply by F matrix
! sub c_mult - multiply by C matrix
! sub i_mult - multiply by sqrt(-1)
! sub solve_f2c2 - solve (F*F+C*C)*x = y
!
! Variable Definitions:
! def l_tlnmc - Logical for TLNMC (set to true if namelist option tlnmc_option
! is 1, 2, or 3
! def reg_tlnmc_type - =1 for regional 1st version of strong constraint
! =2 for regional 2nd version of strong constraint
! def nstrong - number of iterations of strong constraint initialization
! def scheme - which scheme (B, C or D) is being used (see reference above)
! def period_max - max period (hours) of gravity modes to be balanced
! def period_width - width of smooth transition (hours, centered on period_max)
! from balanced to unbalanced gravity modes
! def baldiag_full - flag to toggle balance diagnostics for the full fields
! def baldiag_inc - flag to toggle balance diagnostics for the analysis increment
! def tlnmc_option - Integer option for Incremental Normal Mode Constraint (inmc) / TLNMC
! when in hybrid ensemble mode:
! =0: no constraint at all
! =1: TLNMC on static contribution to increment (or if non-hybrid)
! =2: TLNMC on total increment (single time level only, or 3D mode)
! =3: TLNMC on total increment over all nobs_bins (if 4D mode)
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds,only: r_kind,i_kind
use constants, only: zero,half,one,two,four,r3600,omega,pi,rearth,one_tenth
implicit none
! set default to private
private
! set subroutines to public
public :: init_strongvars
public :: gproj
public :: gproj_diag
public :: gproj_diag_update
public :: gproj0
public :: gproj_ad
public :: dinmi
public :: dinmi_ad
public :: dinmi0
public :: balm_1
public :: getbcf
public :: scale_vars
public :: scale_vars_ad
public :: unscale_vars
public :: unscale_vars_ad
public :: f_mult
public :: c_mult
public :: i_mult
public :: solve_f2c2
! set passed variables to public
public :: nstrong,baldiag_full,l_tlnmc,baldiag_inc,period_width,period_max,scheme
public :: reg_tlnmc_type
public :: tlnmc_option
integer(i_kind) nstrong
integer(i_kind) reg_tlnmc_type,tlnmc_option
real(r_kind) period_max,period_width
logical l_tlnmc,baldiag_full,baldiag_inc
character(1) scheme
contains
subroutine init_strongvars
!$$$ subprogram documentation block
! . . .
! subprogram: init_strongvars
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block
! 2013-07-02 parrish - change tlnmc_type to reg_tlnmc_type. Set default for reg_tlnmc_type = 1.
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
l_tlnmc=.false.
reg_tlnmc_type=1
tlnmc_option=0
nstrong=0
period_max=1000000._r_kind
period_width=one_tenth
scheme='B'
baldiag_full=.false.
baldiag_inc =.false.
end subroutine init_strongvars
subroutine gproj_diag_update(vort,div,phi,vort_g,div_g,phi_g,rmstend,rmstend_g,rmstend_f,rmstend_fg, &
m,mmax,gspeed)
!$$$ subprogram documentation block
! . . .
! subprogram: gproj
!
! prgrmmr:
!
! abstract:
! for gravity wave projection: vort, div, phi --> vort_g, div_g, phi_g
!
! scale: vort,div,phi --> vort_hat,div_hat,phi_hat
!
! solve: (F*F+C*C)*x = F*vort_hat + C*phi_hat
! then:
! phi_hat_g = C*x
! vort_hat_g = F*x
! div_hat_g = div_hat
!
! unscale: vort_hat_g, div_hat_g, phi_hat_g --> vort_g, div_g, phi_g
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort,div,phi
! rmstend,rmstend_g
! filtered
!
! output argument list:
! vort_g,div_g,phi_g
! rmstend,rmstend_g
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(2,m:mmax):: vort,div,phi
real(r_kind),intent( out),dimension(2,m:mmax):: vort_g,div_g,phi_g
real(r_kind),intent(in ) :: gspeed
real(r_kind),intent(inout) :: rmstend,rmstend_g,rmstend_f,rmstend_fg
integer(i_kind),intent(in) :: m,mmax
real(r_kind),dimension(2,m:mmax):: vort_hat,div_hat,phi_hat,vort_hat_g,phi_hat_g
real(r_kind),dimension(m:mmax):: b,c,f,c2,c3
call getbcf(b,c,f,c2,c3,m,mmax,gspeed)
call scale_vars(vort,div,phi,vort_hat,div_hat,phi_hat,.false.,c2,c3,m,mmax,gspeed)
call balm_1(vort_hat,div_hat,phi_hat,rmstend,m,mmax)
call gproj0(vort_hat,phi_hat,vort_hat_g,phi_hat_g,c,f,m,mmax)
call balm_1(vort_hat_g,div_hat,phi_hat_g,rmstend_g,m,mmax)
call scale_vars(vort,div,phi,vort_hat,div_hat,phi_hat,.true.,c2,c3,m,mmax,gspeed)
call balm_1(vort_hat,div_hat,phi_hat,rmstend_f,m,mmax)
call gproj0(vort_hat,phi_hat,vort_hat_g,phi_hat_g,c,f,m,mmax)
call balm_1(vort_hat_g,div_hat,phi_hat_g,rmstend_fg,m,mmax)
call unscale_vars(vort_hat_g,div_hat,phi_hat_g,vort_g,div_g,phi_g,c2,c3,m,mmax)
end subroutine gproj_diag_update
subroutine gproj_diag(vort,div,phi,rmstend,rmstend_g,rmstend_f,rmstend_fg,&
m,mmax,gspeed)
!$$$ subprogram documentation block
! . . .
! subprogram: gproj
!
! prgrmmr:
!
! abstract:
! for gravity wave projection: vort, div, phi --> vort_g, div_g, phi_g
!
! scale: vort,div,phi --> vort_hat,div_hat,phi_hat
!
! solve: (F*F+C*C)*x = F*vort_hat + C*phi_hat
! then:
! phi_hat_g = C*x
! vort_hat_g = F*x
! div_hat_g = div_hat
!
! unscale: vort_hat_g, div_hat_g, phi_hat_g --> vort_g, div_g, phi_g
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort,div,phi
! rmstend,rmstend_g
! filtered
!
! output argument list:
! vort_g,div_g,phi_g
! rmstend,rmstend_g
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(2,m:mmax):: vort,div,phi
real(r_kind),intent(in ) :: gspeed
real(r_kind),intent(inout) :: rmstend,rmstend_g,rmstend_f,rmstend_fg
integer(i_kind),intent(in) :: m,mmax
real(r_kind),dimension(2,m:mmax):: vort_hat,div_hat,phi_hat,vort_hat_g,phi_hat_g
real(r_kind),dimension(m:mmax):: b,c,f,c2,c3
call getbcf(b,c,f,c2,c3,m,mmax,gspeed)
call scale_vars(vort,div,phi,vort_hat,div_hat,phi_hat,.false.,c2,c3,m,mmax,gspeed)
call balm_1(vort_hat,div_hat,phi_hat,rmstend,m,mmax)
call gproj0(vort_hat,phi_hat,vort_hat_g,phi_hat_g,c,f,m,mmax)
call balm_1(vort_hat_g,div_hat,phi_hat_g,rmstend_g,m,mmax)
call scale_vars(vort,div,phi,vort_hat,div_hat,phi_hat,.true.,c2,c3,m,mmax,gspeed)
call balm_1(vort_hat,div_hat,phi_hat,rmstend_f,m,mmax)
call gproj0(vort_hat,phi_hat,vort_hat_g,phi_hat_g,c,f,m,mmax)
call balm_1(vort_hat_g,div_hat,phi_hat_g,rmstend_fg,m,mmax)
end subroutine gproj_diag
subroutine gproj(vort,div,phi,vort_g,div_g,phi_g,m,mmax,gspeed)
!$$$ subprogram documentation block
! . . .
! subprogram: gproj
!
! prgrmmr:
!
! abstract:
! for gravity wave projection: vort, div, phi --> vort_g, div_g, phi_g
!
! scale: vort,div,phi --> vort_hat,div_hat,phi_hat
!
! solve: (F*F+C*C)*x = F*vort_hat + C*phi_hat
! then:
! phi_hat_g = C*x
! vort_hat_g = F*x
! div_hat_g = div_hat
!
! unscale: vort_hat_g, div_hat_g, phi_hat_g --> vort_g, div_g, phi_g
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort,div,phi
! rmstend,rmstend_g
! filtered
!
! output argument list:
! vort_g,div_g,phi_g
! rmstend,rmstend_g
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(2,m:mmax):: vort,div,phi
real(r_kind),intent( out),dimension(2,m:mmax):: vort_g,div_g,phi_g
real(r_kind),intent(in ) :: gspeed
integer(i_kind),intent(in) :: m,mmax
real(r_kind),dimension(2,m:mmax):: vort_hat,div_hat,phi_hat,vort_hat_g,phi_hat_g
real(r_kind),dimension(m:mmax):: b,c,f,c2,c3
call getbcf(b,c,f,c2,c3,m,mmax,gspeed)
call scale_vars(vort,div,phi,vort_hat,div_hat,phi_hat,.true.,c2,c3,m,mmax,gspeed)
call gproj0(vort_hat,phi_hat,vort_hat_g,phi_hat_g,c,f,m,mmax)
call unscale_vars(vort_hat_g,div_hat,phi_hat_g,vort_g,div_g,phi_g,c2,c3,m,mmax)
end subroutine gproj
subroutine gproj0(vort_hat,phi_hat,vort_hat_g,phi_hat_g,c,f,m,mmax)
!$$$ subprogram documentation block
! . . .
! subprogram: gproj0
!
! prgrmmr:
!
! abstract:
! for gravity wave projection: vort,div,phi --> vort_g,div_g,phi_g
! -----------------------------------------------------------------------------
!
! solve: (F*F+C*C)*x = F*vort_hat + C*phi_hat
!
! then:
! phi_hat_g = C*x
! vort_hat_g = F*x
! div_hat_g = div_hat
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort_hat,div_hat,phi_hat
!
! output argument list:
! vort_hat_g,div_hat_g,phi_hat_g
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),dimension(2,m:mmax),intent(in ) :: vort_hat,phi_hat
real(r_kind),dimension( m:mmax),intent(in ) :: c,f
real(r_kind),dimension(2,m:mmax),intent( out) :: vort_hat_g,phi_hat_g
integer(i_kind),intent(in) :: m,mmax
real(r_kind),dimension(2,m:mmax):: x,y
call f_mult(y,vort_hat,f,m,mmax)
call c_mult(x,phi_hat,c,m,mmax)
y=y+x
call solve_f2c2(x,y,f,c,m,mmax)
call c_mult(phi_hat_g,x,c,m,mmax)
call f_mult(vort_hat_g,x,f,m,mmax)
end subroutine gproj0
subroutine gproj_ad(vort,div,phi,vort_g,div_g,phi_g,m,mmax,gspeed)
!$$$ subprogram documentation block
! . . .
! subprogram: gproj_ad
!
! prgrmmr:
!
! abstract:
! for gravity wave projection: vort,div,phi --> vort_g,div_g,phi_g
!
! scale: vort,div,phi --> vort_hat,div_hat,phi_hat
!
! solve: (F*F+C*C)*x = F*vort_hat + C*phi_hat
!
! then:
! phi_hat_g = C*x
! vort_hat_g = F*x
! div_hat_g = div_hat
!
! unscale: vort_hat_g, div_hat_g, phi_hat_g --> vort_g, div_g, phi_g
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort,div,phi
! vort_g,div_g,phi_g
!
! output argument list:
! vort,div,phi
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(inout),dimension(2,m:mmax) :: vort,div,phi
real(r_kind),intent(in ),dimension(2,m:mmax) :: vort_g,div_g,phi_g
real(r_kind),intent(in ) :: gspeed
integer(i_kind),intent(in) :: m,mmax
real(r_kind),dimension(2,m:mmax)::vort_hat,div_hat,phi_hat,vort_hat_g,phi_hat_g
real(r_kind),dimension(m:mmax):: b,c,f,c2,c3
call getbcf(b,c,f,c2,c3,m,mmax,gspeed)
call unscale_vars_ad(vort_hat_g,div_hat,phi_hat_g,vort_g,div_g,phi_g,c2,c3,m,mmax)
call gproj0(vort_hat_g,phi_hat_g,vort_hat,phi_hat,c,f,m,mmax)
call scale_vars_ad(vort,div,phi,vort_hat,div_hat,phi_hat,c2,c3,m,mmax,gspeed)
end subroutine gproj_ad
subroutine dinmi(vort_t,div_t,phi_t,del_vort,del_div,del_phi,m,mmax,gspeed)
!$$$ subprogram documentation block
! . . .
! subprogram: dinmi
!
! prgrmmr:
!
! abstract:
! for implicit nmi correction: vort_t,div_t,phi_t --> del_vort,del_div,del_phi
!
! scale: vort_t,div_t,phi_t --> vort_t_hat,div_t_hat,phi_t_hat
!
! solve: (F*F+C*C)*del_div_hat = sqrt(-1)*(F*vort_t_hat + C*phi_t_hat)
!
! solve: (F*F+C*C)*x = sqrt(-1)*div_t_hat - B*del_div_hat
!
! then:
! del_phi_hat = C*x
! del_vort_hat = F*x
!
! unscale: del_vort_hat,del_div_hat,del_phi_hat --> del_vort,del_div,del_phi
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort_t,div_t,phi_t
!
! output argument list:
! del_vort,del_div,del_phi
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(2,m:mmax):: vort_t,div_t,phi_t
real(r_kind),intent( out),dimension(2,m:mmax):: del_vort,del_div,del_phi
real(r_kind),intent(in ):: gspeed
integer(i_kind),intent(in) :: m,mmax
real(r_kind),dimension(2,m:mmax):: vort_t_hat,div_t_hat,phi_t_hat
real(r_kind),dimension(2,m:mmax):: del_vort_hat,del_div_hat,del_phi_hat
real(r_kind),dimension(m:mmax):: b,c,f,c2,c3
call getbcf(b,c,f,c2,c3,m,mmax,gspeed)
call scale_vars(vort_t,div_t,phi_t,vort_t_hat,div_t_hat,phi_t_hat,.true.,c2,c3,m,mmax,gspeed)
call dinmi0(vort_t_hat,div_t_hat,phi_t_hat,del_vort_hat,del_div_hat,del_phi_hat,b,c,f,m,mmax)
call unscale_vars(del_vort_hat,del_div_hat,del_phi_hat,del_vort,del_div,del_phi,c2,c3,m,mmax)
end subroutine dinmi
subroutine dinmi_ad(vort_t,div_t,phi_t,del_vort,del_div,del_phi,m,mmax,gspeed)
!$$$ subprogram documentation block
! . . .
! subprogram: dinmi_ad
!
! prgrmmr:
!
! abstract:
! for implicit nmi correction: vort_t,div_t,phi_t --> del_vort,del_div,del_phi
!
! scale: vort_t,div_t,phi_t --> vort_t_hat,div_t_hat,phi_t_hat
!
! solve: (F*F+C*C)*del_div_hat = sqrt(-1)*(F*vort_t_hat + C*phi_t_hat)
!
! solve: (F*F+C*C)*x = sqrt(-1)*div_t_hat - B*del_div_hat
!
! then:
! del_phi_hat = C*x
! del_vort_hat = F*x
!
! unscale: del_vort_hat,del_div_hat,del_phi_hat --> del_vort,del_div,del_phi
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! del_vort,del_div,del_phi
! vort_t,div_t,phi_t
!
! output argument list:
! vort_t,div_t,phi_t
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(inout),dimension(2,m:mmax):: vort_t,div_t,phi_t
real(r_kind),intent(in ),dimension(2,m:mmax):: del_vort,del_div,del_phi
real(r_kind),intent(in ):: gspeed
integer(i_kind),intent(in) :: m,mmax
real(r_kind),dimension(2,m:mmax):: vort_t_hat,div_t_hat,phi_t_hat
real(r_kind),dimension(2,m:mmax):: del_vort_hat,del_div_hat,del_phi_hat
real(r_kind),dimension(m:mmax):: b,c,f,c2,c3
integer(i_kind) n
call getbcf(b,c,f,c2,c3,m,mmax,gspeed)
call unscale_vars_ad(del_vort_hat,del_div_hat,del_phi_hat,del_vort,del_div,del_phi,c2,c3,m,mmax)
call dinmi0(del_vort_hat,del_div_hat,del_phi_hat,vort_t_hat,div_t_hat,phi_t_hat,b,c,f,m,mmax)
do n=m,mmax
vort_t_hat(1,n)=-vort_t_hat(1,n)
vort_t_hat(2,n)=-vort_t_hat(2,n)
div_t_hat(1,n)=-div_t_hat(1,n)
div_t_hat(2,n)=-div_t_hat(2,n)
phi_t_hat(1,n)=-phi_t_hat(1,n)
phi_t_hat(2,n)=-phi_t_hat(2,n)
end do
call scale_vars_ad(vort_t,div_t,phi_t,vort_t_hat,div_t_hat,phi_t_hat,c2,c3,m,mmax,gspeed)
end subroutine dinmi_ad
subroutine dinmi0(vort_t_hat,div_t_hat,phi_t_hat,del_vort_hat,del_div_hat,del_phi_hat,b,c,f,&
m,mmax)
!$$$ subprogram documentation block
! . . .
! subprogram: dinmi0
!
! prgrmmr:
!
! abstract:
! for implicit nmi correction: vort_t,div_t,phi_t --> del_vort,del_div,del_phi
!
! solve: (F*F+C*C)*del_div_hat = sqrt(-1)*(F*vort_t_hat + C*phi_t_hat)
!
! solve: (F*F+C*C)*x = sqrt(-1)*div_t_hat - B*del_div_hat
!
! then:
! del_phi_hat = C*x
! del_vort_hat = F*x
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort_t_hat,div_t_hat,phi_t_hat
!
! output argument list:
! del_vort_hat,del_div_hat,del_phi_hat
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),dimension(2,m:mmax),intent(in ) :: vort_t_hat,div_t_hat,phi_t_hat
real(r_kind),dimension( m:mmax),intent(in ) :: b,c,f
real(r_kind),dimension(2,m:mmax),intent( out) :: del_vort_hat,del_div_hat,del_phi_hat
integer(i_kind),intent(in)::m,mmax
real(r_kind),dimension(2,m:mmax):: x,y
call f_mult(y,vort_t_hat,f,m,mmax)
call c_mult(x,phi_t_hat,c,m,mmax)
x=y+x
call i_mult(y,x,m,mmax)
call solve_f2c2(del_div_hat,y,f,c,m,mmax)
call c_mult(x,del_div_hat,b,m,mmax) ! actually multiplying by b
call i_mult(y,div_t_hat,m,mmax)
y=y-x
call solve_f2c2(x,y,f,c,m,mmax)
call c_mult(del_phi_hat,x,c,m,mmax)
call f_mult(del_vort_hat,x,f,m,mmax)
end subroutine dinmi0
subroutine balm_1(vort_t_hat,div_t_hat,phi_t_hat,balnm1,m,mmax)
!$$$ subprogram documentation block
! . . .
! subprogram: gproj
!
! prgrmmr:
!
! abstract: obtain balance diagnostic for each wave number n,m using
! method 1 (eq 4.23 of Temperton,1989)
!
! balnm1 = abs(vort_t_hat)(n,m)**2 + abs(div_t_hat)(n,m)**2 +
! abs(phi_t_hat)(n,m)**2
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
! 2011-07-01 todling -- balm1 must not be zeroed out
!
! input argument list:
! vort_t_hat,div_t_hat,phi_t_hat
!
! output argument list:
! balnm1
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(2,m:mmax) :: vort_t_hat,div_t_hat,phi_t_hat
real(r_kind),intent(inout) :: balnm1
integer(i_kind),intent(in)::m,mmax
integer(i_kind) n
do n=m,mmax
balnm1=balnm1+vort_t_hat(1,n)**2+vort_t_hat(2,n)**2 &
+div_t_hat(1,n)**2+ div_t_hat(2,n)**2 &
+phi_t_hat(1,n)**2+ phi_t_hat(2,n)**2
end do
end subroutine balm_1
subroutine getbcf(b,c,f,c2,c3,m,mmax,gspeed)
!$$$ subprogram documentation block
! . . .
! subprogram: getbcf
!
! prgrmmr:
!
! abstract: compute operators needed to do gravity wave projection
! and implicit normal mode initialization in spectral space
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
!
! output argument list:
! b,c,f,c2,c3
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent( out),dimension(m:mmax) :: b,c,f,c2,c3
real(r_kind),intent(in ) ::gspeed
integer(i_kind),intent(in)::m,mmax
integer(i_kind) n,nstart
real(r_kind) eps,rn,rm,rn1
! scheme B: b = 2*omega*m/(n*(n+1))
! f = 2*omega*sqrt(n*n-1)*eps/n
! c = gspeed*sqrt(n*(n+1))/erad
! scheme C: b = 0
! f = 2*omega*sqrt(n*n-1)*eps/n
! c = gspeed*sqrt(n*(n+1))/erad
! scheme D: b = 0
! f = 2*omega*eps
! c = gspeed*sqrt(n*(n+1))/erad
! in the above, eps = sqrt((n*n-m*m)/(4*n*n-1))
nstart=max(m,1)
rm=m
do n=nstart,mmax
rn=n
eps=sqrt((rn*rn-rm*rm )/(four*rn*rn-one))
rn1=sqrt(rn*(rn+one))
if(scheme=='B') then
b(n)=two*omega*rm/(rn*(rn+one))
f(n)=two*omega*sqrt(rn*rn-one)*eps/rn
c(n)=gspeed*rn1/rearth
c2(n)=rearth/rn1
c3(n)=one/gspeed
else if(scheme=='C') then
b(n)=zero
f(n)=two*omega*sqrt(rn*rn-one)*eps/rn
c(n)=gspeed*rn1/rearth
c2(n)=rearth/rn1
c3(n)=one/gspeed
else if(scheme=='D') then
b(n)=zero
f(n)=two*omega*eps
c(n)=gspeed*rn1/rearth
c2(n)=rearth
c3(n)=rn1/gspeed
else
write(6,*)' scheme = ',scheme,' incorrect, must be = B, C, or D'
end if
end do
end subroutine getbcf
subroutine scale_vars(vort,div,phi,vort_hat,div_hat,phi_hat,filtered,c2,c3,&
m,mmax,gspeed)
!$$$ subprogram documentation block
! . . .
! subprogram: scale_vars
!
! prgrmmr:
!
! abstract:
! input scaling:
!
! for schemes B, C:
! vort_hat(n) = erad*vort(n)/sqrt(n*(n+1))
! div_hat(n) = sqrt(-1)*erad*div(n)/sqrt(n*(n+1))
! phi_hat(n) = phi(n)/gspeed
!
! for scheme D:
! vort_hat(n) = erad*vort(n)
! div_hat(n) = sqrt(-1)*erad*div(n)
! phi_hat(n) = phi(n)*sqrt(n*(n+1))/gspeed
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort,div,phi
!
! output argument list:
! vort_hat,div_hat,phi_hat
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(2,m:mmax):: vort,div,phi
real(r_kind),intent(in ),dimension(m:mmax):: c2,c3
real(r_kind),intent(in ) :: gspeed
real(r_kind),intent( out),dimension(2,m:mmax):: vort_hat,div_hat,phi_hat
integer(i_kind),intent(in)::m,mmax
logical,intent(in)::filtered
real(r_kind) pmask,c1,pmax
integer(i_kind) n,nstart
! following is to account for 0,0 term being zero
vort_hat(1,m)=zero
vort_hat(2,m)=zero
div_hat(1,m)=zero
div_hat(2,m)=zero
phi_hat(1,m)=zero
phi_hat(2,m)=zero
nstart=max(m,1)
c1=two*pi*rearth/(gspeed*r3600*period_width)
pmax=period_max/period_width
do n=nstart,mmax
pmask=one
if(filtered)pmask=half*(one-tanh(c1/real(n,r_kind)-pmax))
vort_hat(1,n)=pmask*c2(n)*vort(1,n)
vort_hat(2,n)=pmask*c2(n)*vort(2,n)
div_hat(2,n) =pmask*c2(n)*div (1,n)
div_hat(1,n) =-pmask*c2(n)*div(2,n)
phi_hat(1,n) =pmask*c3(n)*phi(1,n)
phi_hat(2,n) =pmask*c3(n)*phi(2,n)
end do
end subroutine scale_vars
subroutine scale_vars_ad(vort,div,phi,vort_hat,div_hat,phi_hat,c2,c3,&
m,mmax,gspeed)
!$$$ subprogram documentation block
! . . .
! subprogram: scale_vars_ad
!
! prgrmmr:
!
! abstract:
! input scaling:
! for schemes B, C:
! vort_hat(n) = erad*vort(n)/sqrt(n*(n+1))
! div_hat(n) = sqrt(-1)*erad*div(n)/sqrt(n*(n+1))
! phi_hat(n) = phi(n)/gspeed
!
! for scheme D:
! vort_hat(n) = erad*vort(n)
! div_hat(n) = sqrt(-1)*erad*div(n)
! phi_hat(n) = phi(n)*sqrt(n*(n+1))/gspeed
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort,div,phi
! vort_hat,div_hat,phi_hat
!
! output argument list:
! vort,div,phi
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(inout),dimension(2,m:mmax):: vort,div,phi
real(r_kind),intent(in ),dimension(m:mmax):: c2,c3
real(r_kind),intent(in ) :: gspeed
real(r_kind),intent(in ),dimension(2,m:mmax):: vort_hat,div_hat,phi_hat
integer(i_kind),intent(in)::m,mmax
integer(i_kind) n,nstart
real(r_kind) pmask,pmax,c1
nstart=max(m,1)
c1=two*pi*rearth/(gspeed*r3600*period_width)
pmax=period_max/period_width
do n=nstart,mmax
pmask=half*(one-tanh(c1/real(n,r_kind)-pmax))
vort(1,n)=vort(1,n)+pmask*c2(n)*vort_hat(1,n)
vort(2,n)=vort(2,n)+pmask*c2(n)*vort_hat(2,n)
div(1,n) =div(1,n) +pmask*c2(n)*div_hat (2,n)
div(2,n) =div(2,n) -pmask*c2(n)*div_hat (1,n)
phi(1,n) =phi(1,n) +pmask*c3(n)*phi_hat (1,n)
phi(2,n) =phi(2,n) +pmask*c3(n)*phi_hat (2,n)
end do
end subroutine scale_vars_ad
subroutine unscale_vars(vort_hat,div_hat,phi_hat,vort,div,phi,c2,c3,m,mmax)
!$$$ subprogram documentation block
! . . .
! subprogram: unscale_vars
!
! prgrmmr:
!
! abstract:
! output scaling:
! for schemes B, C:
! vort(n,m) = sqrt(n*(n+1))*vort_hat(n,m)/erad
! div(n,m) = -sqrt(-1)*sqrt(n*(n+1))*div_hat(n,m)/erad
! phi(n,m) = gspeed*phi_hat(n,m)
! for scheme C:
! vort(n,m) = vort_hat(n,m)/erad
! div(n,m) = -sqrt(-1)*div_hat(n,m)/erad
! phi(n,m) = gspeed*phi_hat(n,m)/sqrt(n*(n+1))
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort_hat,div_hat,phi_hat
!
! output argument list:
! vort,div,phi
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(2,m:mmax):: vort_hat,div_hat,phi_hat
real(r_kind),intent(in ),dimension(m:mmax):: c2,c3
real(r_kind),intent( out),dimension(2,m:mmax):: vort,div,phi
integer(i_kind),intent(in)::m,mmax
integer(i_kind) n,nstart
! following is to account for 0,0 term being zero
vort(1,m)=zero
vort(2,m)=zero
div(1,m)=zero
div(2,m)=zero
phi(1,m)=zero
phi(2,m)=zero
nstart=max(m,1)
do n=nstart,mmax
vort(1,n)=vort_hat(1,n)/c2(n)
vort(2,n)=vort_hat(2,n)/c2(n)
div (1,n)=div_hat (2,n)/c2(n)
div (2,n)=-div_hat(1,n)/c2(n)
phi (1,n)=phi_hat (1,n)/c3(n)
phi (2,n)=phi_hat (2,n)/c3(n)
end do
end subroutine unscale_vars
subroutine unscale_vars_ad(vort_hat,div_hat,phi_hat,vort,div,phi,c2,c3,m,mmax)
!$$$ subprogram documentation block
! . . .
! subprogram: unscale_vars_ad
!
! prgrmmr:
!
! abstract:
! output scaling:
! for schemes B, C:
! vort(n,m) = sqrt(n*(n+1))*vort_hat(n,m)/erad
! div(n,m) = -sqrt(-1)*sqrt(n*(n+1))*div_hat(n,m)/erad
! phi(n,m) = gspeed*phi_hat(n,m)
! for scheme C:
! vort(n,m) = vort_hat(n,m)/erad
! div(n,m) = -sqrt(-1)*div_hat(n,m)/erad
! phi(n,m) = gspeed*phi_hat(n,m)/sqrt(n*(n+1))
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! vort,div,phi
! vort_hat,div_hat,phi_hat
!
! output argument list:
! vort_hat,div_hat,phi_hat
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(inout),dimension(2,m:mmax):: vort_hat,div_hat,phi_hat
real(r_kind),intent(in ),dimension(2,m:mmax):: vort,div,phi
real(r_kind),intent(in ),dimension(m:mmax):: c2,c3
integer(i_kind),intent(in)::m,mmax
integer(i_kind) n,nstart
! following is to account for 0,0 term being zero
vort_hat(1,m)=zero
vort_hat(2,m)=zero
div_hat(1,m)=zero
div_hat(2,m)=zero
phi_hat(1,m)=zero
phi_hat(2,m)=zero
nstart=max(m,1)
do n=nstart,mmax
vort_hat(1,n)=vort(1,n)/c2(n)
vort_hat(2,n)=vort(2,n)/c2(n)
div_hat (1,n)=-div (2,n)/c2(n)
div_hat (2,n)=div (1,n)/c2(n)
phi_hat (1,n)=phi (1,n)/c3(n)
phi_hat (2,n)=phi (2,n)/c3(n)
end do
end subroutine unscale_vars_ad
subroutine f_mult(x,y,f,m,mmax)
!$$$ subprogram documentation block
! . . .
! subprogram: f_mult
!
! prgrmmr:
!
! abstract: x = F*y
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! y,f
!
! output argument list:
! x
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(m:mmax):: f
real(r_kind),intent(in ),dimension(2,m:mmax):: y
real(r_kind),intent( out),dimension(2,m:mmax):: x
integer(i_kind),intent(in)::m,mmax
integer(i_kind) n,nstart
if(m==mmax) then
x=zero
return
end if
! following is to account for 0,0 term being zero
x(1,m)=zero
x(2,m)=zero
nstart=max(m,1)
x(1,mmax)=zero
x(2,mmax)=zero
if(nstart<mmax) then
do n=nstart,mmax-1
x(1,n)=f(n+1)*y(1,n+1)
x(2,n)=f(n+1)*y(2,n+1)
end do
do n=nstart+1,mmax
x(1,n)=x(1,n)+f(n)*y(1,n-1)
x(2,n)=x(2,n)+f(n)*y(2,n-1)
end do
end if
end subroutine f_mult
subroutine c_mult(x,y,c,m,mmax)
!$$$ subprogram documentation block
! . . .
! subprogram: c_mult
!
! prgrmmr:
!
! abstract: x = C*y
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! y,c
!
! output argument list:
! x
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(m:mmax):: c
real(r_kind),intent(in ),dimension(2,m:mmax):: y
real(r_kind),intent( out),dimension(2,m:mmax):: x
integer(i_kind),intent(in)::m,mmax
integer(i_kind) n,nstart
! following is to account for 0,0 term being zero
x(1,m)=zero
x(2,m)=zero
nstart=max(m,1)
do n=nstart,mmax
x(1,n)=c(n)*y(1,n)
x(2,n)=c(n)*y(2,n)
end do
end subroutine c_mult
subroutine i_mult(x,y,m,mmax)
!$$$ subprogram documentation block
! . . .
! subprogram: i_mult
!
! prgrmmr:
!
! abstract: x = sqrt(-1)*y
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! y
!
! output argument list:
! x
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(2,m:mmax):: y
real(r_kind),intent( out),dimension(2,m:mmax):: x
integer(i_kind),intent(in)::m,mmax
integer(i_kind) n,nstart
! following is to account for 0,0 term being zero
x(1,m)=zero
x(2,m)=zero
nstart=max(m,1)
do n=nstart,mmax
x(1,n)=-y(2,n)
x(2,n)=y(1,n)
end do
end subroutine i_mult
subroutine solve_f2c2(x,y,f,c,m,mmax)
!$$$ subprogram documentation block
! . . .
! subprogram: solve_f2c2
!
! prgrmmr:
!
! abstract: solve (F*F+C*C)*x = y
!
! program history log:
! 2008-05-05 safford -- add subprogram doc block, rm unused uses
!
! input argument list:
! y,f,c
!
! output argument list:
! x
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ),dimension(m:mmax):: f,c
real(r_kind),intent(in ),dimension(2,m:mmax):: y
real(r_kind),intent( out),dimension(2,m:mmax):: x
integer(i_kind),intent(in)::m,mmax
integer(i_kind) n,nstart
real(r_kind),dimension(m:mmax):: a,b
real(r_kind),dimension(2,m:mmax):: z
! following is to account for 0,0 term being zero
x(1,m)=zero
x(2,m)=zero
nstart=max(m,1)
! copy forcing y to internal array
do n=nstart,mmax
z(1,n)=y(1,n)
z(2,n)=y(2,n)
end do
! if nstart==mmax, then trivial solution
if(nstart==mmax) then
a( nstart)=c( nstart)*c(nstart)
x(1,nstart)=z(1,nstart)/a(nstart)
x(2,nstart)=z(2,nstart)/a(nstart)
else
! compute main diagonal of F*F + C*C
a(nstart)=f(nstart+1)*f(nstart+1)+c(nstart)*c(nstart)
if(nstart+1 < mmax) then
do n=nstart+1,mmax-1
a(n)=f(n)*f(n)+f(n+1)*f(n+1)+c(n)*c(n)
end do
end if
a(mmax)=f(mmax)*f(mmax)+c(mmax)*c(mmax)
! compute only non-zero off-diagonal of F*F + C*C
if(nstart+2<=mmax) then
do n=nstart+2,mmax
b(n)=f(n-1)*f(n)
a(n)=a(n)-b(n)*b(n)/a(n-2)
end do
! forward elimination:
do n=nstart,mmax-2
z(1,n+2)=z(1,n+2)-b(n+2)*z(1,n)/a(n)
z(2,n+2)=z(2,n+2)-b(n+2)*z(2,n)/a(n)
end do
end if
! backward substitution:
x(1,mmax )=z(1,mmax )/a(mmax )
x(2,mmax )=z(2,mmax )/a(mmax )
x(1,mmax-1)=z(1,mmax-1)/a(mmax-1)
x(2,mmax-1)=z(2,mmax-1)/a(mmax-1)
if(nstart+2 <= mmax) then
do n=mmax-2,nstart,-1
x(1,n)=(z(1,n) - b(n+2)*x(1,n+2))/a(n)
x(2,n)=(z(2,n) - b(n+2)*x(2,n+2))/a(n)
end do
end if
end if
end subroutine solve_f2c2
end module mod_strong