module mod_inmi !$$$ module documentation block ! . . . . ! module: obsmod ! prgmmr: derber org: np23 date: 2003-09-25 ! ! abstract: 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: ! 2006-08-15 parrish ! ! Subroutines Included: ! sub gproj - project input u,v,mass variable to gravity modes ! 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 get_periodmask - create mask to only balance gravity modes with periods ! less than period_max ! 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 ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds,only: r_kind,i_kind implicit none ! set default to private private ! set subroutines to public public :: gproj public :: gproj0 public :: gproj_ad public :: dinmi public :: dinmi_ad public :: dinmi0 public :: balm_1 public :: get_periodmask 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 :: mmax,gspeed,m integer(i_kind) mmax integer(i_kind) m real(r_kind) gspeed contains subroutine gproj(vort,div,phi,vort_g,div_g,phi_g,rmstend,rmstend_g,filtered) !$$$ 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 ) :: vort(2,m:mmax),div(2,m:mmax),phi(2,m:mmax) real(r_kind),intent( out) :: vort_g(2,m:mmax),div_g(2,m:mmax),phi_g(2,m:mmax) real(r_kind),intent(inout) :: rmstend,rmstend_g logical ,intent(in ) :: filtered real(r_kind) pmask(m:mmax) real(r_kind) vort_hat(2,m:mmax),div_hat(2,m:mmax),phi_hat(2,m:mmax) real(r_kind) vort_hat_g(2,m:mmax),div_hat_g(2,m:mmax),phi_hat_g(2,m:mmax) real(r_kind) rmstendnm(m:mmax) integer(i_kind) n call scale_vars(vort,div,phi,vort_hat,div_hat,phi_hat) if(filtered) then call get_periodmask(pmask) do n=m,mmax vort_hat(1,n)=pmask(n)*vort_hat(1,n) vort_hat(2,n)=pmask(n)*vort_hat(2,n) div_hat(1,n)=pmask(n)*div_hat(1,n) div_hat(2,n)=pmask(n)*div_hat(2,n) phi_hat(1,n)=pmask(n)*phi_hat(1,n) phi_hat(2,n)=pmask(n)*phi_hat(2,n) end do end if call balm_1(vort_hat,div_hat,phi_hat,rmstendnm) do n=m,mmax rmstend=rmstend+rmstendnm(n) end do call gproj0(vort_hat,div_hat,phi_hat,vort_hat_g,div_hat_g,phi_hat_g) call balm_1(vort_hat_g,div_hat_g,phi_hat_g,rmstendnm) do n=m,mmax rmstend_g=rmstend_g+rmstendnm(n) end do call unscale_vars(vort_hat_g,div_hat_g,phi_hat_g,vort_g,div_g,phi_g) end subroutine gproj subroutine gproj0(vort_hat,div_hat,phi_hat,vort_hat_g,div_hat_g,phi_hat_g) !$$$ 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,div_hat,phi_hat real(r_kind),dimension(2,m:mmax),intent( out) :: vort_hat_g,div_hat_g,phi_hat_g real(r_kind) b(m:mmax),c(m:mmax),f(m:mmax) real(r_kind) x(2,m:mmax),y(2,m:mmax) call getbcf(b,c,f) call f_mult(y,vort_hat,f) call c_mult(x,phi_hat,c) y=y+x call solve_f2c2(x,y,f,c) call c_mult(phi_hat_g,x,c) call f_mult(vort_hat_g,x,f) div_hat_g=div_hat end subroutine gproj0 subroutine gproj_ad(vort,div,phi,vort_g,div_g,phi_g) !$$$ 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 use constants, only: zero implicit none real(r_kind),intent(inout) :: vort(2,m:mmax),div(2,m:mmax),phi(2,m:mmax) real(r_kind),intent(in ) :: vort_g(2,m:mmax),div_g(2,m:mmax),phi_g(2,m:mmax) real(r_kind) pmask(m:mmax) real(r_kind) vort_hat(2,m:mmax),div_hat(2,m:mmax),phi_hat(2,m:mmax) real(r_kind) vort_hat_g(2,m:mmax),div_hat_g(2,m:mmax),phi_hat_g(2,m:mmax) integer(i_kind) n vort_hat_g=zero div_hat_g=zero phi_hat_g=zero call unscale_vars_ad(vort_hat_g,div_hat_g,phi_hat_g,vort_g,div_g,phi_g) call gproj0(vort_hat_g,div_hat_g,phi_hat_g,vort_hat,div_hat,phi_hat) call get_periodmask(pmask) do n=m,mmax vort_hat(1,n)=pmask(n)*vort_hat(1,n) vort_hat(2,n)=pmask(n)*vort_hat(2,n) div_hat(1,n)=pmask(n)*div_hat(1,n) div_hat(2,n)=pmask(n)*div_hat(2,n) phi_hat(1,n)=pmask(n)*phi_hat(1,n) phi_hat(2,n)=pmask(n)*phi_hat(2,n) end do call scale_vars_ad(vort,div,phi,vort_hat,div_hat,phi_hat) end subroutine gproj_ad subroutine dinmi(vort_t,div_t,phi_t,del_vort,del_div,del_phi) !$$$ 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 ) :: vort_t(2,m:mmax),div_t(2,m:mmax),phi_t(2,m:mmax) real(r_kind),intent( out) :: del_vort(2,m:mmax),del_div(2,m:mmax),del_phi(2,m:mmax) real(r_kind) pmask(m:mmax) real(r_kind) vort_t_hat(2,m:mmax),div_t_hat(2,m:mmax),phi_t_hat(2,m:mmax) real(r_kind) del_vort_hat(2,m:mmax),del_div_hat(2,m:mmax),del_phi_hat(2,m:mmax) integer(i_kind) n call scale_vars(vort_t,div_t,phi_t,vort_t_hat,div_t_hat,phi_t_hat) call get_periodmask(pmask) do n=m,mmax vort_t_hat(1,n)=pmask(n)*vort_t_hat(1,n) vort_t_hat(2,n)=pmask(n)*vort_t_hat(2,n) div_t_hat(1,n)=pmask(n)*div_t_hat(1,n) div_t_hat(2,n)=pmask(n)*div_t_hat(2,n) phi_t_hat(1,n)=pmask(n)*phi_t_hat(1,n) phi_t_hat(2,n)=pmask(n)*phi_t_hat(2,n) end do call dinmi0(vort_t_hat,div_t_hat,phi_t_hat,del_vort_hat,del_div_hat,del_phi_hat) call unscale_vars(del_vort_hat,del_div_hat,del_phi_hat,del_vort,del_div,del_phi) end subroutine dinmi subroutine dinmi_ad(vort_t,div_t,phi_t,del_vort,del_div,del_phi) !$$$ 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 use constants, only: zero implicit none real(r_kind),intent(inout) :: vort_t(2,m:mmax),div_t(2,m:mmax),phi_t(2,m:mmax) real(r_kind),intent(in ) :: del_vort(2,m:mmax),del_div(2,m:mmax),del_phi(2,m:mmax) real(r_kind) pmask(m:mmax) real(r_kind) vort_t_hat(2,m:mmax),div_t_hat(2,m:mmax),phi_t_hat(2,m:mmax) real(r_kind) del_vort_hat(2,m:mmax),del_div_hat(2,m:mmax),del_phi_hat(2,m:mmax) integer(i_kind) n del_vort_hat=zero ; del_div_hat=zero ; del_phi_hat=zero call unscale_vars_ad(del_vort_hat,del_div_hat,del_phi_hat,del_vort,del_div,del_phi) call dinmi0(del_vort_hat,del_div_hat,del_phi_hat,vort_t_hat,div_t_hat,phi_t_hat) call get_periodmask(pmask) do n=m,mmax vort_t_hat(1,n)=-pmask(n)*vort_t_hat(1,n) vort_t_hat(2,n)=-pmask(n)*vort_t_hat(2,n) div_t_hat(1,n)=-pmask(n)*div_t_hat(1,n) div_t_hat(2,n)=-pmask(n)*div_t_hat(2,n) phi_t_hat(1,n)=-pmask(n)*phi_t_hat(1,n) phi_t_hat(2,n)=-pmask(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) end subroutine dinmi_ad subroutine dinmi0(vort_t_hat,div_t_hat,phi_t_hat,del_vort_hat,del_div_hat,del_phi_hat) !$$$ 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(2,m:mmax),intent( out) :: del_vort_hat,del_div_hat,del_phi_hat real(r_kind) b(m:mmax),c(m:mmax),f(m:mmax) real(r_kind) x(2,m:mmax),y(2,m:mmax) call getbcf(b,c,f) call f_mult(y,vort_t_hat,f) call c_mult(x,phi_t_hat,c) x=y+x call i_mult(y,x) call solve_f2c2(del_div_hat,y,f,c) call c_mult(x,del_div_hat,b) ! actually multiplying by b call i_mult(y,div_t_hat) y=y-x call solve_f2c2(x,y,f,c) call c_mult(del_phi_hat,x,c) call f_mult(del_vort_hat,x,f) end subroutine dinmi0 subroutine balm_1(vort_t_hat,div_t_hat,phi_t_hat,balnm1) !$$$ 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 ! ! 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 ) :: vort_t_hat(2,m:mmax),div_t_hat(2,m:mmax),phi_t_hat(2,m:mmax) real(r_kind),intent( out) :: balnm1(m:mmax) integer(i_kind) n do n=m,mmax balnm1(n)=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 get_periodmask(pmask) !$$$ subprogram documentation block ! . . . ! subprogram: gproj ! ! prgrmmr: ! ! abstract: create mask to zero out components with estimated periods longer ! than period_max, where period_max is in hours ! ! c = phase speed for given vertical mode ! ! L/c = period, where L is wavelength ! ! L = 2*pi*erad/n where n is wavenumber and 2*pi*erad is ! circumference of earth ! ! program history log: ! 2008-05-05 safford -- add subprogram doc block, rm unused uses ! ! input argument list: ! ! output argument list: ! pmask ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use constants, only: izero,zero,half,one,two,four,rearth,r3600 use mod_strong, only: period_max,period_width implicit none real(r_kind),intent( out) :: pmask(m:mmax) real(r_kind) pi,thislength,thisperiod integer(i_kind) n pi=four*atan(one) do n=m,mmax pmask(n)=zero if(n==izero) cycle thislength=two*pi*rearth/n thisperiod=thislength/(gspeed*r3600) pmask(n)=half*(one-tanh((thisperiod-period_max)/period_width)) end do end subroutine get_periodmask subroutine getbcf(b,c,f) !$$$ 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 ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use constants, only: ione,zero,one,two,four,omega,rearth use mod_strong, only: scheme implicit none real(r_kind),intent( out) :: b(m:mmax),c(m:mmax),f(m:mmax) integer(i_kind) n,nstart real(r_kind) eps,rn,rm ! 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,ione) rm=m do n=nstart,mmax rn=n eps=sqrt((rn*rn-rm*rm )/(four*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*sqrt(rn*(rn+one))/rearth else if(scheme=='C') then b(n)=zero f(n)=two*omega*sqrt(rn*rn-one)*eps/rn c(n)=gspeed*sqrt(rn*(rn+one))/rearth else if(scheme=='D') then b(n)=zero f(n)=two*omega*eps c(n)=gspeed*sqrt(rn*(rn+one))/rearth 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) !$$$ 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 use constants, only: ione,zero,one,rearth use mod_strong, only: scheme implicit none real(r_kind),intent(in ) :: vort(2,m:mmax),div(2,m:mmax),phi(2,m:mmax) real(r_kind),intent( out) :: vort_hat(2,m:mmax),div_hat(2,m:mmax),phi_hat(2,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,ione) if(scheme/='D') then do n=nstart,mmax vort_hat(1,n)=rearth*vort(1,n)/sqrt(n*(n+one)) vort_hat(2,n)=rearth*vort(2,n)/sqrt(n*(n+one)) div_hat(2,n) =rearth*div (1,n)/sqrt(n*(n+one)) div_hat(1,n) =-rearth*div(2,n)/sqrt(n*(n+one)) phi_hat(1,n) =phi(1,n)/gspeed phi_hat(2,n) =phi(2,n)/gspeed end do else do n=nstart,mmax vort_hat(1,n)=rearth*vort(1,n) vort_hat(2,n)=rearth*vort(2,n) div_hat(2,n) =rearth *div(1,n) div_hat(1,n) =-rearth*div(2,n) phi_hat(1,n) =sqrt(n*(n+one))*phi(1,n)/gspeed phi_hat(2,n) =sqrt(n*(n+one))*phi(2,n)/gspeed end do end if end subroutine scale_vars subroutine scale_vars_ad(vort,div,phi,vort_hat,div_hat,phi_hat) !$$$ 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 use constants, only: ione,one,rearth use mod_strong, only: scheme implicit none real(r_kind),intent(inout) :: vort(2,m:mmax),div(2,m:mmax),phi(2,m:mmax) real(r_kind),intent(in ) :: vort_hat(2,m:mmax),div_hat(2,m:mmax),phi_hat(2,m:mmax) integer(i_kind) n,nstart nstart=max(m,ione) if(scheme/='D') then do n=nstart,mmax vort(1,n)=vort(1,n)+vort_hat(1,n)*rearth/sqrt(n*(n+one)) vort(2,n)=vort(2,n)+vort_hat(2,n)*rearth/sqrt(n*(n+one)) div(1,n) =div(1,n) +div_hat (2,n)*rearth/sqrt(n*(n+one)) div(2,n) =div(2,n) -div_hat (1,n)*rearth/sqrt(n*(n+one)) phi(1,n) =phi(1,n) +phi_hat (1,n)/gspeed phi(2,n) =phi(2,n) +phi_hat (2,n)/gspeed end do else do n=nstart,mmax vort(1,n)=vort(1,n)+rearth*vort_hat(1,n) vort(2,n)=vort(2,n)+rearth*vort_hat(2,n) div(1,n) =div (1,n)+div_hat(2,n)*rearth div(2,n) =div (2,n)-div_hat(1,n)*rearth phi(1,n) =phi (1,n)+sqrt(n*(n+one))*phi_hat(1,n)/gspeed phi(2,n) =phi (2,n)+sqrt(n*(n+one))*phi_hat(2,n)/gspeed end do end if end subroutine scale_vars_ad subroutine unscale_vars(vort_hat,div_hat,phi_hat,vort,div,phi) !$$$ 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 use constants, only: ione,zero,one,rearth use mod_strong, only: scheme implicit none real(r_kind),intent(in ) :: vort_hat(2,m:mmax),div_hat(2,m:mmax),phi_hat(2,m:mmax) real(r_kind),intent( out) :: vort(2,m:mmax),div(2,m:mmax),phi(2,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,ione) if(scheme/='D') then do n=nstart,mmax vort(1,n)=vort_hat(1,n)*sqrt(n*(n+one))/rearth vort(2,n)=vort_hat(2,n)*sqrt(n*(n+one))/rearth div (1,n)=div_hat (2,n)*sqrt(n*(n+one))/rearth div (2,n)=-div_hat(1,n)*sqrt(n*(n+one))/rearth phi (1,n)=phi_hat (1,n)*gspeed phi (2,n)=phi_hat (2,n)*gspeed end do else do n=nstart,mmax vort(1,n)=vort_hat(1,n)/rearth vort(2,n)=vort_hat(2,n)/rearth div (1,n)=div_hat (2,n)/rearth div (2,n)=-div_hat(1,n)/rearth phi (1,n)=phi_hat (1,n)*gspeed/sqrt(n*(n+one)) phi (2,n)=phi_hat (2,n)*gspeed/sqrt(n*(n+one)) end do end if end subroutine unscale_vars subroutine unscale_vars_ad(vort_hat,div_hat,phi_hat,vort,div,phi) !$$$ 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 use constants, only: ione,one,rearth use mod_strong, only: scheme implicit none real(r_kind),intent(inout) :: vort_hat(2,m:mmax),div_hat(2,m:mmax),phi_hat(2,m:mmax) real(r_kind),intent(in ) :: vort(2,m:mmax),div(2,m:mmax),phi(2,m:mmax) integer(i_kind) n,nstart nstart=max(m,ione) if(scheme/='D') then do n=nstart,mmax vort_hat(1,n)=vort_hat(1,n)+vort(1,n)*sqrt(n*(n+one))/rearth vort_hat(2,n)=vort_hat(2,n)+vort(2,n)*sqrt(n*(n+one))/rearth div_hat (2,n)=div_hat (2,n)+div (1,n)*sqrt(n*(n+one))/rearth div_hat (1,n)=div_hat (1,n)-div (2,n)*sqrt(n*(n+one))/rearth phi_hat (1,n)=phi_hat (1,n)+phi (1,n)*gspeed phi_hat (2,n)=phi_hat (2,n)+phi (2,n)*gspeed end do else do n=nstart,mmax vort_hat(1,n)=vort_hat(1,n)+vort(1,n)/rearth vort_hat(2,n)=vort_hat(2,n)+vort(2,n)/rearth div_hat (2,n)=div_hat (2,n)+div (1,n)/rearth div_hat (1,n)=div_hat (1,n)-div (2,n)/rearth phi_hat (1,n)=phi_hat (1,n)+phi (1,n)*gspeed/sqrt(n*(n+one)) phi_hat (2,n)=phi_hat (2,n)+phi (2,n)*gspeed/sqrt(n*(n+one)) end do end if end subroutine unscale_vars_ad subroutine f_mult(x,y,f) !$$$ 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 use constants, only: ione,zero implicit none real(r_kind),intent(in ) :: y(2,m:mmax),f(m:mmax) real(r_kind),intent( out) :: x(2,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,ione) x(1,mmax)=zero x(2,mmax)=zero if(nstart