module general_specmod !$$$ module documentation block ! . . . . ! module: general_specmod ! prgmmr: treadon org: np23 date: 2003-11-24 ! ! abstract: copy of specmod, introducing structure variable spec_vars, so ! spectral code can be used for arbitrary resolutions. ! ! program history log: ! 2003-11-24 treadon ! 2004-04-28 d. kokron, updated SGI's fft to use scsl ! 2004-05-18 kleist, documentation ! 2004-08-27 treadon - add/initialize variables/arrays needed by ! splib routines for grid <---> spectral ! transforms ! 2007-04-26 yang - based on idrt value xxxx descriptionxxx ! 2010-02-18 parrish - copy specmod to general_specmod and add structure variable spec_vars. ! remove all *_b variables, since now not necessary to have two ! resolutions. any number of resolutions can be now contained in ! type(spec_vars) variables passed in through init_spec_vars. also ! remove init_spec, since not really necessary. ! 2011-05-01 rancic - add parameters jcap_trunc, nc_trunc ! 2014-12-03 derber - add alp0 array to spectral structures to eliminate ! repeated calculations ! ! subroutines included: ! sub general_init_spec_vars ! sub general_destroy_spec_vars ! ! remarks: variable definitions below ! def jcap - spectral (assumed triangular) truncation ! def nc - (N+1)*(N+2); N=truncation ! def nc1 - 2*(N+1); N=truncation ! def ncd2 - [(N+1)*(N+2)]/2; N=truncation ! def jnpe - (N+2)/2; N=truncation ! def factsml - factor to ensure proper scalar coefficients are zero ! def factvml - factor to ensure proper vector coefficients are zero ! def iromb - integer spectral domain shape ! (0 for triangular, 1 for rhomboidal) ! def idrt - integer grid identifier ! (idrt=4 for gaussian grid, ! idrt=0 for equally-spaced grid including poles, ! idrt=256 for equally-spaced grid excluding poles) ! def imax - integer even number of longitudes for transform ! def jmax - integer number of latitudes for transform ! def ijmax - integer imax*jmax ! def jn - integer skip number between n.h. latitudes from north ! def js - integer skip number between s.h. latitudes from south ! def kw - integer skip number between wave fields ! def jb - integer latitude index (from pole) to begin transform ! def je - integer latitude index (from pole) to end transform ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,r_double,i_kind implicit none ! set default to private private ! set subroutines to public public :: general_init_spec_vars public :: general_destroy_spec_vars ! set passed variables to public public :: spec_vars public :: spec_cut integer(i_kind) :: spec_cut type spec_vars integer(i_kind) jcap integer(i_kind) nc integer(i_kind) jcap_trunc integer(i_kind) nc_trunc integer(i_kind) ncd2 integer(i_kind) iromb integer(i_kind) idrt integer(i_kind) imax integer(i_kind) jmax integer(i_kind) ijmax integer(i_kind) jn integer(i_kind) js integer(i_kind) kw integer(i_kind) jb integer(i_kind) je integer(i_kind) ioffset logical, pointer :: factsml(:) => NULL() logical, pointer :: factvml(:) => NULL() real(r_kind), pointer :: eps(:) => NULL() real(r_kind), pointer :: epstop(:) => NULL() real(r_kind), pointer :: enn1(:) => NULL() real(r_kind), pointer :: elonn1(:) => NULL() real(r_kind), pointer :: eon(:) => NULL() real(r_kind), pointer :: eontop(:) => NULL() real(r_kind), pointer :: clat(:) => NULL() real(r_kind), pointer :: slat(:) => NULL() real(r_kind), pointer :: wlat(:) => NULL() real(r_kind), pointer :: pln(:,:) => NULL() real(r_kind), pointer :: plntop(:,:) => NULL() real(r_kind), pointer :: test_mask(:) => NULL() real(r_kind), pointer :: rlats(:) => NULL() real(r_kind), pointer :: slats(:) => NULL() real(r_kind), pointer :: clats(:) => NULL() real(r_kind), pointer :: rlons(:) => NULL() real(r_kind), pointer :: slons(:) => NULL() real(r_kind), pointer :: clons(:) => NULL() real(r_kind), pointer :: alp0(:) => NULL() real(r_double),pointer :: afft(:) => NULL() logical:: lallocated = .false. logical:: precalc_pln = .true. end type spec_vars contains subroutine general_init_spec_vars(sp,jcap,jcap_test,nlat_a,nlon_a,eqspace) !$$$ subprogram documentation block ! . . . . ! subprogram: init_spec_vars ! prgmmr: treadon org: np23 date: 2003-11-24 ! ! abstract: initialize spectral variables ! ! program history log: ! 2003-11-24 treadon ! 2004-05-18 kleist, new variables and documentation ! 2004-08-27 treadon - add call to sptranf0 and associated arrays, ! remove del21 and other unused arrays/variables ! 2006-04-06 middlecoff - remove jc=ncpus() since not used ! 2008-04-11 safford - rm unused vars ! 2010-02-18 parrish - substantial changes to simplify and introduce input/output variable ! type(spec_vars) sp ! 2010-04-01 treadon - remove mpimod and rad2deg constants (not used) ! 2013-10-23 el akkraoui - initialize lats to zero (otherwise point is undefined) ! ! input argument list: ! sp - type(spec_vars) variable ! jcap - target resolution ! jcap_test - test resolution, used to construct mask which will zero out coefs ! with total wavenumber n in range jcap_test < n <= jcap ! nlat_a - number of latitudes on target grid ! nlon_a - number of longitudes on target grid ! eqspace - optional flag for equal spaced grid instead of gaussian grid ! ! output argument list: ! sp - type(spec_vars) variable, ready to use for target spectral resolution/grid. ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ use constants, only: zero,half,one,two,pi,three implicit none ! Declare passed variables type(spec_vars) ,intent(inout) :: sp integer(i_kind) ,intent(in ) :: jcap,jcap_test,nlat_a,nlon_a logical,optional,intent(in ) :: eqspace ! Declare local variables integer(i_kind) i,ii1,j,l,m,jhe,n integer(i_kind) :: ldafft real(r_kind) :: dlon_a,half_pi,two_pi real(r_kind),dimension(nlat_a-2) :: wlatx,slatx real(r_kind) :: epsi0(0:jcap) ! epsilon factor for m=0 real(r_kind) :: fnum, fden ! Set constants used in transforms for analysis grid sp%jcap=jcap sp%nc=(jcap+1)*(jcap+2) sp%jcap_trunc=jcap sp%nc_trunc=(sp%jcap_trunc+1)*(sp%jcap_trunc+2) sp%ncd2=sp%nc/2 sp%iromb=0 sp%idrt=4 if(present(eqspace)) then if(eqspace) sp%idrt=256 endif sp%imax=nlon_a sp%jmax=nlat_a-2 sp%ijmax=sp%imax*sp%jmax sp%ioffset=sp%imax*(sp%jmax-1) sp%jn=sp%imax sp%js=-sp%jn sp%kw=2*sp%ncd2 sp%jb=1 sp%je=(sp%jmax+1)/2 ! Allocate and initialize fact arrays if(sp%lallocated) then deallocate(sp%factsml,sp%factvml,sp%eps,sp%epstop,sp%enn1,sp%elonn1,sp%eon,sp%eontop) deallocate(sp%clat,sp%slat,sp%wlat,sp%pln,sp%plntop,sp%test_mask,sp%afft) deallocate(sp%rlats,sp%clats,sp%slats,sp%alp0) deallocate(sp%rlons,sp%clons,sp%slons) end if allocate(sp%factsml(sp%nc),sp%factvml(sp%nc),sp%test_mask(sp%nc)) sp%factsml=.false. sp%factvml=.false. sp%test_mask=one ii1=0 do l=0,sp%jcap do m=0,sp%jcap-l ii1=ii1+2 if(l == 0)then sp%factsml(ii1)=.true. sp%factvml(ii1)=.true. end if if(l+m.gt.jcap_test) then sp%test_mask(ii1-1)=zero sp%test_mask(ii1)=zero end if end do end do sp%factvml(1)=.true. ! Allocate and initialize arrays used in transforms allocate( sp%eps(sp%ncd2) ) allocate( sp%epstop(sp%jcap+1) ) allocate( sp%enn1(sp%ncd2) ) allocate( sp%elonn1(sp%ncd2) ) allocate( sp%eon(sp%ncd2) ) allocate( sp%eontop(sp%jcap+1) ) ldafft=50000+4*sp%imax ! ldafft=256+imax would be sufficient at GMAO. allocate( sp%afft(ldafft)) allocate( sp%clat(sp%jb:sp%je) ) allocate( sp%slat(sp%jb:sp%je) ) allocate( sp%wlat(sp%jb:sp%je) ) call spwget(sp%iromb,sp%jcap,sp%eps,sp%epstop,sp%enn1, & sp%elonn1,sp%eon,sp%eontop) call spffte(sp%imax,(sp%imax+2)/2,sp%imax,2,0.,0.,0,sp%afft) call splat(sp%idrt,sp%jmax,slatx,wlatx) jhe=(sp%jmax+1)/2 if(jhe > sp%jmax/2)wlatx(jhe)=wlatx(jhe)/2 do j=sp%jb,sp%je sp%clat(j)=sqrt(1.-slatx(j)**2) sp%slat(j)=slatx(j) sp%wlat(j)=wlatx(j) end do if(sp%jcap < spec_cut)then sp%precalc_pln=.true. allocate( sp%pln(sp%ncd2,sp%jb:sp%je) ) allocate( sp%plntop(sp%jcap+1,sp%jb:sp%je) ) do j=sp%jb,sp%je call splegend(sp%iromb,sp%jcap,sp%slat(j),sp%clat(j),sp%eps, & sp%epstop,sp%pln(1,j),sp%plntop(1,j)) end do else sp%precalc_pln=.false. allocate( sp%pln(sp%ncd2,1) ) allocate( sp%plntop(sp%jcap+1,1) ) end if ! obtain rlats and rlons half_pi=half*pi two_pi=two*pi allocate(sp%rlats(nlat_a),sp%rlons(nlon_a)) allocate(sp%clats(nlat_a),sp%clons(nlon_a)) allocate(sp%slats(nlat_a),sp%slons(nlon_a)) allocate(sp%alp0(0:sp%jcap)) sp%rlats=zero sp%clats=zero sp%slats=zero sp%rlats(1)=-half_pi sp%clats(1)=zero sp%slats(1)=-one sp%rlats(nlat_a)=half_pi sp%clats(nlat_a)=zero sp%slats(nlat_a)=one do i=1,(nlat_a-2)/2 sp%rlats(nlat_a-i)= asin(sp%slat(i)) sp%clats(nlat_a-i)= sp%clat(i) sp%slats(nlat_a-i)= sp%slat(i) sp%rlats(1+i )=-asin(sp%slat(i)) sp%clats(1+i )= sp%clat(i) sp%slats(1+i )= -sp%slat(i) end do dlon_a=two_pi/nlon_a do j=1,nlon_a sp%rlons(j)=(j-one)*dlon_a sp%clons(j)=cos(sp%rlons(j)) sp%slons(j)=sin(sp%rlons(j)) end do ! Compute epsilon for m=0. epsi0(0)=zero do n=1,sp%jcap fnum=real(n**2) fden=real(4*n**2-1) epsi0(n)=sqrt(fnum/fden) enddo ! ! Compute Legendre polynomials for m=0 at North Pole sp%alp0(0)=sqrt(half) sp%alp0(1)=sqrt(three)*sp%alp0(0) do n=2,sp%jcap sp%alp0(n)=(sp%alp0(n-1)-epsi0(n-1)*sp%alp0(n-2))/epsi0(n) enddo sp%lallocated=.true. return end subroutine general_init_spec_vars subroutine general_destroy_spec_vars(sp) !$$$ subprogram documentation block ! . . . . ! subprogram: destroy_spec_vars ! prgmmr: treadon org: np23 date: 2003-11-24 ! ! abstract: deallocate memory from spectral variables ! ! program history log: ! 2003-11-24 treadon ! 2004-05-18 kleist, new variables and documentation ! 2010-02-18 parrish, copy destroy_spec_vars to general_destroy_spec_vars ! ! input argument list: ! sp - type(spec_vars) variable, which is no longer required ! ! output argument list: ! sp - arrays in sp deleted ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ implicit none type(spec_vars),intent(inout) :: sp if(sp%lallocated) then deallocate(sp%factsml,sp%factvml) deallocate(sp%eps,sp%epstop,sp%enn1,sp%elonn1,sp%eon,sp%eontop,sp%afft,& sp%clat,sp%slat,sp%wlat) deallocate(sp%pln) deallocate(sp%plntop) deallocate(sp%rlats,sp%rlons) deallocate(sp%clats,sp%clons) deallocate(sp%slats,sp%slons) deallocate(sp%alp0) sp%lallocated=.false. end if return end subroutine general_destroy_spec_vars end module general_specmod