! ------------------------------------------------------------------------------- ! INPE/CPTEC, DMD, Modelling and Development Division ! ------------------------------------------------------------------------------- ! ! modULE: Mod_LegendreTransform ! ! REVISION HISTORY: ! - Jose P. Bonatti - version: 1 ! 01-08-2007 - Tomita - version: 1.1.1.1 ! 01-04-2018 - Daniel M. Lamosa - version: 2.0 ! Barbara A. G. P. Yamada - version: 2.0 ! 08-04-2019 - Eduardo Khamis !> @author !> Jose P. Bonatti \n !> Daniel M. Lamosa (last revision) \n !> Barbara A. G. P. Yamada (last revision) \n !> Eduardo Khamis \n !! !> @brief !! !! !! @version 2.0 !! @date 01-04-2018 !! !! @copyright Under GLP-3.0 !! @link: https://opensource.org/licenses/GPL-3.0 ! ------------------------------------------------------------------------------- module LegendreTransform use typeKinds, only: p_r8 => Double, p_i4 => Long use ModConstants, only: emRad use ModConstants, only: emRad1 use ModConstants, only: emRad2 use ModConstants, only: emRad12 implicit none private integer, parameter :: p_nferr = 0 public :: getxMaxyMax type, public :: legendre ! private logical :: created integer :: Mend integer :: Mend1 integer :: Mend2 integer :: Mend3 integer :: Mends integer :: MnWv0 integer :: MnWv1 integer :: MnWv2 integer :: MnWv3 integer :: xMax integer :: yMax integer :: yMaxHf integer :: xMx integer(p_i4), pointer, dimension(:,:) :: la0 => null() integer(p_i4), pointer, dimension(:,:) :: la1 => null() real(kind=p_r8), pointer, dimension(:) :: e0 => null() real(kind=p_r8), pointer, dimension(:) :: e1 => null() real(kind=p_r8), pointer, dimension(:) :: eps => null() real(kind=p_r8), pointer, dimension(:) :: colRad => null() real(kind=p_r8), pointer, dimension(:) :: rCs2 => null() real(kind=p_r8), pointer, dimension(:) :: wgt => null() real(kind=p_r8), pointer, dimension(:) :: gLats => null() real(kind=p_r8), pointer, dimension(:) :: snnp => null() real(kind=p_r8), pointer, dimension(:,:) :: legS2F => null() real(kind=p_r8), pointer, dimension(:,:) :: legExtS2F => null() real(kind=p_r8), pointer, dimension(:,:) :: legDerS2F => null() real(kind=p_r8), pointer, dimension(:,:) :: legDerNS => null() real(kind=p_r8), pointer, dimension(:,:) :: legDerEW => null() contains procedure, public :: initLegendre => initialization_ procedure, public :: destroyLegendre => destroylegendre procedure, public :: transp => specTransp_ procedure, public :: SymAsy => SymAsy_ procedure, public :: Four2Spec => Fourier2SpecCoef_ procedure, public :: Spec2Four => SpecCoef2Fourier_ procedure, public :: GetSize procedure, private :: createSpectralRep procedure, private :: createGaussRep procedure, private :: createLegPols procedure, private :: DestroySpectralRep procedure, private :: DestroyGaussRep procedure, private :: DestroyLegPols end type contains function initialization_(self, mend, imax, jmax) result(iret) class(legendre) :: self integer :: Mend integer, optional :: imax integer, optional :: jmax integer :: iret iret = 0 ! precisa implementar isso self%Mend = Mend self%Mend1 = Mend + 1 self%Mend2 = Mend + 2 self%Mend3 = Mend + 3 self%Mends = 2*self%Mend1 self%MnWv2 = (Mend + 1) * (Mend + 2) self%MnWv0 = self%MnWv2 / 2 self%MnWv3 = self%MnWv2 + 2 * (self%Mend1) self%MnWv1 = self%MnWv3 / 2 call getxMaxyMax(Mend, self%xMax, self%yMax) if (present(imax)) self%xMax = imax if (present(jmax)) self%yMax = jmax self%yMaxHf = self%yMax/2 self%xmx = self%xMax + 2 self%created = .false. call self%createSpectralRep( ) call self%createGaussRep( ) call self%createLegPols( ) end function function destroylegendre(self) result(iret) class(legendre) :: self integer :: iret iret = 0 ! precisa implementar isso self%Mend = -1 self%Mend1 = -1 self%Mend2 = -1 self%MnWv0 = -1 self%MnWv1 = -1 self%MnWv2 = -1 self%MnWv3 = -1 self%xMax = -1 self%yMax = -1 self%yMaxHf= -1 self%created = .false. call self%destroySpectralRep( ) call self%destroyGaussRep( ) call self%destroyLegPols( ) end function ! --------------------------------------------------------------------------- !> @brief createSpectralRep !! !! @details !! !! @author Jose P. Bonatti \n !> Daniel M. Lamosa (last revision) \n !> Barbara A. G. P. Yamada (last revision) !> Eduardo Khamis 08-04-2019 !! !! @date mar/2018 ! --------------------------------------------------------------------------- subroutine createSpectralRep(self) class(legendre) :: self integer :: k integer :: l !< integer :: m !< integer :: n !< real(kind=p_r8) :: am !< real(kind=p_r8) :: an !< real(kind=p_r8) :: sn !< allocate(self%la0(self%Mend1,self%Mend1)) allocate(self%la1(self%Mend1,self%Mend2)) allocate(self%eps(self%mnwv1)) allocate(self%e0(self%mnwv1)) allocate(self%e1(self%mnwv0)) allocate(self%snnp(self%mnwv2)) l = 0 k = 0 do n=1, self%mend1 l = l + 1 self%la1(n,1) = l self%eps(l) = 0.0_p_r8 self%e0(l) = 0.0_p_r8 self%e1(l) = emRad/REAL(l,p_r8) do m=1, self%mend2 - n k = k + 1 self%la0(m,n) = k sn=REAL(m+n-2,p_r8) if(sn.ne.0.0_p_r8)then sn = -EMRad2/(sn*(sn+1.0_p_r8)) self%snnp(2*l-1) = sn self%snnp(2*l ) = sn else self%snnp(2*l-1) = 0.0_p_r8 self%snnp(2*l ) = 0.0_p_r8 endif if(n.ge.2)then an=n+m-2 am=m-1 self%e1(k) = emRad*am/(an+an*an) endif end do end do self%e1(1) = 0.0_p_r8 l = self%mend1 do n=2, self%mend2 do m=1, self%mend1 + 2 - n l = l + 1 self%la1(m,n) = l am = m - 1 an = m + n - 2 self%eps(l) = sqrt((an * an - am * am) / (4.0_p_r8 * an * an - 1.0_p_r8)) self%e0(l) = emRad*self%eps(l)/REAL(n+m-2,p_r8) end do end do end subroutine createSpectralRep ! --------------------------------------------------------------------------- !> @brief createGaussRep !! !! @details !! !! @author Jose P. Bonatti \n !> Daniel M. Lamosa (last revision) \n !> Barbara A. G. P. Yamada (last revision) !> Eduardo Khamis 09-04-2019 !! !! @date mar/2018 ! --------------------------------------------------------------------------- subroutine createGaussRep(self) class(legendre) :: self allocate(self%colRad(self%yMaxHf)) allocate(self%rCs2(self%yMaxHf)) allocate(self%wgt(self%yMaxHf)) allocate(self%gLats(self%yMax)) call gaussianLatitudes(self) end subroutine createGaussRep subroutine CreateLegPols(self) class(legendre) :: self allocate(self%legS2F (self%mnwv0, self%yMaxHf)) allocate(self%legExtS2F(self%mnwv1, self%yMaxHf)) allocate(self%legDerS2F(self%mnwv0, self%yMaxHf)) allocate(self%legDerNS (self%mnwv0, self%yMaxHf)) allocate(self%legDerEW (self%mnwv0, self%yMaxHf)) call legPols(self) endsubroutine ! --------------------------------------------------------------------------- !> @brief gLats !! !! @details Calculates Gaussian Latitudes and Gaussian Weights !! for Use in Grid-Spectral and Spectral-Grid Transforms !! !! @author Jose P. Bonatti \n !> Daniel M. Lamosa (last revision) \n !> Barbara A. G. P. Yamada (last revision) !> Eduardo Khamis 09-04-2019 !! !! @date mar/2018 ! --------------------------------------------------------------------------- ! ! subroutine gaussianLatitudes(self) ! glats: calculates gaussian latitudes and ! gaussian weights for use in grid-spectral ! and spectral-grid transforms. ! ! glats calls the subroutine poly ! ! argument(dimensions) description ! ! colrad(Jmaxhf) output: co-latitudes for gaussian ! latitudes in one hemisphere. ! wgt(Jmaxhf) output: gaussian weights. ! rcs2(Jmaxhf) output: 1.0/cos(gaussian latitude)**2 ! Jmaxhf input: number of gaussian latitudes ! in one hemisphere. type(legendre), intent(inout) :: self real(kind=p_r8) :: epsil, scal, dgColIn integer :: j !< Loop iterator real(kind=p_r8) :: gCol !< real(kind=p_r8) :: dgCol !< real(kind=p_r8) :: p2 !< real(kind=p_r8) :: p1 !< real(kind=p_r8) :: rad !< epsil = epsilon(1.0_p_r8) * 100.0_p_r8 rad = 45.0_p_r8 / atan(1.0_p_r8) scal = 2.0_p_r8 / (real(self%yMax, p_r8) * real(self%yMax, p_r8)) dgColIn = atan(1.0_p_r8) / real(self%yMax, p_r8) gCol = 0.0_p_r8 do j=1, self%yMaxHf dgCol = dgColIn do call legendrePolynomial(self%yMax, gCol, p2) do p1 = p2 gCol = gCol + dgCol call legendrePolynomial(self%yMax, gCol, p2) if(sign(1.0_p_r8, p1) /= sign(1.0_p_r8, p2)) exit end do if(dgCol <= epsil) exit gCol = gCol - dgCol dgCol = dgCol * 0.25_p_r8 end do self%colRad(j) = gCol self%gLats(j) = 90.0_p_r8 - rad * gCol self%gLats(self%yMax - j + 1) = -self%gLats(j) call legendrePolynomial(self%yMax-1, gCol, p1) self%wgt(j) = scal * (1.0_p_r8 - cos(gCol) * cos(gCol)) / (p1 * p1) self%rCs2(j) = 1.0_p_r8 / (sin(gCol) * sin(gCol)) end do end subroutine gaussianLatitudes ! --------------------------------------------------------------------------- !> @brief legendrePolynomial !! !! @details Calculates the Value of the Ordinary Legendre !! Function of Given Order at a Specified Colatitude. !! Used to Determine Gaussian Latitudes. !! !! @author Jose P. Bonatti \n !> Daniel M. Lamosa (last revision) \n !> Barbara A. G. P. Yamada (last revision) !! !! @date mar/2018 ! --------------------------------------------------------------------------- subroutine legendrePolynomial(n, rad, pln) ! legendrePolynomial: ! calculates the value of the ordinary legendre function ! of given order at a specified latitude. used to ! determine gaussian latitudes. ! !*********************************************************************** ! ! This routine is called by the subroutine glats. ! ! this routine calls no subroutines. ! !*********************************************************************** ! ! argument(dimensions) description ! ! n input : order of the ordinary legendre ! function whose value is to be ! calculated. set in routine ! "glats". ! rad input : colatitude (in radians) at ! which the value of the ordinar ! legendre function is to be ! calculated. set in routine ! "glats". ! pln output : value of the ordinary legendre ! function of order "n" at ! colatitude "rad". ! !*********************************************************************** ! integer, intent(in) :: n !< real(kind=p_r8), intent(in) :: rad !< real(kind=p_r8), intent(out) :: pln !< integer :: i !< Loop iterator real(kind=p_r8) :: x !< real(kind=p_r8) :: y1 !< real(kind=p_r8) :: y2 !< real(kind=p_r8) :: y3 !< real(kind=p_r8) :: g !< x = cos(rad) y1 = 1.0_p_r8 y2 = x do i=2, n g = x * y2 y3 = g - y1 + g - (g - y1) / real(i, p_r8) y1 = y2 y2 = y3 end do pln = y3 end subroutine legendrePolynomial ! --------------------------------------------------------------------------- !> @brief legPols !! !! @details !! !! !! @author Jose P. Bonatti \n !> Daniel M. Lamosa (last revision) \n !> Barbara A. G. P. Yamada (last revision) !> Eduardo Khamis 19-04-2019 !! !! @date mar/2018 ! --------------------------------------------------------------------------- subroutine legPols(self) type(legendre), intent(inout) :: self integer :: j !< integer :: l !< integer :: nn !< integer :: mm !< integer :: mn !< integer :: lx !< real(kind=p_r8) :: pln(self%mnwv1) do j=1, self%yMaxHf call pln2(self, pln, j) l = 0 do nn=1, self%mend1 do mm=1, self%mend2 - nn l = l + 1 lx = self%la1(mm,nn) self%legS2F(l,j) = pln(lx) ! just half part end do end do do mn=1, self%mnwv1 self%legExtS2F(mn,j) = pln(mn) end do call plnder(self, j, pln, & self%legDerS2F(:,j), & self%legDerNS (:,j), & self%legDerEW (:,j) & ) end do end subroutine legPols subroutine pln2 (self, sln, plat) ! pln2: calculates the associated legendre functions ! at one specified latitude. ! ! plat input: gaussian latitude index. set ! by calling routine. ! sln (Mnwv1) output: values of associated legendre ! functions at one gaussian ! latitude specified by the ! argument "lat". ! uses from self: ! colrad(Jmaxhf) colatitudes of gaussian grid ! (in radians). calculated ! in routine "glats". ! eps factor that appears in recusio ! formula of a.l.f. ! Mend triangular truncation wave number ! Mnwv1 number of elements ! la1(Mend1,Mend1+1) numbering array of pln1 !-------------------------------------------------------! type(legendre), intent(in) :: self real(kind=p_r8), intent(out) :: sln(self%mnwv1) integer, intent(in) :: plat integer :: mm, nn, lx, ly, lz real(kind=p_r8) :: colr, sinlat, coslat, prod logical, SAVE :: first = .true. real(kind=p_r8), allocatable, dimension(:), SAVE :: x, y real(kind=p_r8), SAVE :: rthf if (first) then allocate (x(self%mend1)) allocate (y(self%mend1)) first=.false. do mm=1,self%mend1 x(mm)=sqrt(2.0_p_r8*mm+1.0_p_r8) y(mm)=sqrt(1.0_p_r8+0.5_p_r8/real(mm,p_r8)) end do rthf=sqrt(0.5_p_r8) endif colr=self%colrad(plat) sinlat=cos(colr) coslat=sin(colr) prod=1.0_p_r8 do mm=1,self%mend1 sln(mm)=rthf*prod ! line below should only be used where exponent range is limted ! if(prod < flim) prod=0.0_p_r8 prod=prod*coslat*y(mm) end do do mm=1,self%mend1 sln(mm+self%mend1)=x(mm)*sinlat*sln(mm) end do do nn=3,self%mend2 do mm=1,self%mend1+2-nn lx=self%la1(mm,nn) ly=self%la1(mm,nn-1) lz=self%la1(mm,nn-2) sln(lx)=(sinlat*sln(ly)-self%eps(ly)*sln(lz))/self%eps(lx) end do end do end subroutine pln2 subroutine plnder(self, plat, pln, dpln, der, plnwcs) ! calculates zonal and meridional pseudo-derivatives as ! well as laplacians of the associated legendre functions. ! ! argument(dimensions) description ! ! pln (Mnwv1) input : associated legendre function ! values at gaussian latitude ! output : pln(l,n)= ! pln(l,n)*(l+n-2)*(l+n-1)/a**2. ! ("a" denotes earth's radius) ! used in calculating the ! laplacian of global fields ! in spectral form. ! dpln (Mnwv0) output : derivatives of ! associated legendre functions ! at gaussian latitude ! der (Mnwv0) output : cosine-weighted derivatives of ! associated legendre functions ! at gaussian latitude ! plnwcs(Mnwv0) output : plnwcs(l,n)= ! pln(l,n)*(l-1)/sin(latitude)**2. ! !uses from self: ! ! rcs2l input : 1.0/sin(latitude)**2 at ! gaussian latitude ! wgtl input : gaussian weight, at gaussian ! latitude ! eps (Mnwv1) input : array of constants used to ! calculate "der" from "pln". ! computed in routine "epslon". ! la1(Mend1,Mend2) input : numbering array for pln type(legendre), intent(in ) :: self integer, intent(in ) :: plat real(kind=p_r8), intent(inout) :: pln(self%mnwv1) real(kind=p_r8), intent( out) :: dpln(self%mnwv0) real(kind=p_r8), intent( out) :: der(self%mnwv0) real(kind=p_r8), intent( out) :: plnwcs(self%mnwv0) integer :: n, l, nn, mm, mn, lm, l0, lp real(kind=p_r8) :: raa, wcsa real(kind=p_r8) :: x(self%mnwv1) logical, SAVE :: first=.true. real(kind=p_r8), allocatable, SAVE :: an(:) ! ! compute pln derivatives ! if (first) then allocate (an(self%mend2)) do n=1,self%mend2 an(n)=real(n-1,p_r8) end do first=.false. end if raa=self%wgt(plat)*emRad12 wcsa=self%rCs2(plat)*self%wgt(plat)*emRad1 l=0 do mm=1,self%mend1 l=l+1 x(l)=an(mm) end do do nn=2,self%mend2 do mm=1,self%mend1+2-nn l=l+1 x(l)=an(mm+nn-1) end do end do l=self%mend1 do nn=2,self%mend1 do mm=1,self%mend2-nn l=l+1 lm=self%la1(mm,nn-1) l0=self%la1(mm,nn) lp=self%la1(mm,nn+1) der(l)=x(lp)*self%eps(l0)*pln(lm)-x(l0)*self%eps(lp)*pln(lp) end do end do do mm=1,self%mend1 der(mm)=-x(mm)*self%eps(mm+self%Mend1)*pln(mm+self%Mend1) end do do mn=1,self%mnwv0 dpln(mn) = der(mn) der(mn) = wcsa*der(mn) end do l=0 do nn=1,self%mend1 do mm=1,self%mend2-nn l=l+1 l0=self%la1(mm,nn) plnwcs(l)=an(mm)*pln(l0) end do end do do mn=1,self%mnwv0 plnwcs(mn)=wcsa*plnwcs(mn) end do do nn=1,self%mend1 do mm=1,self%mend2-nn l0=self%la1(mm,nn) lp=self%la1(mm,nn+1) pln(l0)=x(l0)*x(lp)*raa*pln(l0) end do end do end subroutine plnder subroutine SpecTransp_ (self, itype, a, MnWv) ! transp: after input, transposes scalar arrays ! of spectral coefficients by swapping ! the order of the subscripts ! representing the degree and order ! of the associated legendre functions. ! ! argument(dimensions) description ! ! a(:) input: spectral representation of a ! global field. Size should be ! mnwv2 or mnwv3 ! signal=+1 diagonalwise storage ! signal=-1 coluMnwise storage ! output: spectral representation of a ! global field. ! signal=+1 coluMnwise storage ! signal=-1 diagonalwise storage class(legendre), intent(in ) :: self integer, intent(in ) :: itype integer, intent(in ) :: MnWv real (kind=p_r8), intent(inout) :: a(MnWv) integer :: l, lx, n, m, mn integer :: mMax integer, pointer :: idx(:,:) => null() real (kind=p_r8), allocatable :: qwork(:) if(MnWv == self%mnwv2)then idx => self%la0 mMax = self%Mend2 elseif(MnWv == self%mnwv3)then idx => self%la1 mMax = self%Mend2 + 1 else write(*,'(A,1x,I6)')'>> error: wrong size of spectral coefficents:',MnWv write(*,'(A,1x,I6,1x,A,1x,I6)')'Size should be:',self%MnWv2,'or',self%MnWv3 return endif allocate(qwork(MnWv)) qwork=0.0_p_r8 if (itype == +1) then l=0 do m=1,self%mend1 do n=1,mMax-m l=l+1 lx=idx(m,n) qwork( 2*l-1 ) = a( 2*lx-1 ) qwork( 2*l ) = a( 2*lx ) end do end do do mn=1,mnwv a(mn) = qwork(mn) end do else if(itype == -1)then l=0 do m=1,self%mend1 do n=1,mMax-m l=l+1 lx=idx(m,n) qwork( 2*lx-1 ) = a( 2*l-1 ) qwork( 2*lx ) = a( 2*l ) end do end do do mn=1,mnwv a(mn) = qwork(mn) end do end if deallocate(qwork) end subroutine SpecTransp_ subroutine Fourier2SpecCoef_ (self, fp, fm, lat, MnWv, fln) ! calculates spectral representations of ! global fields from fourier representations ! of symmetric and anti-symmetric portions ! of global fields. ! ! argument(dimensions) description ! ! fp(Imax+2) input: fourier representation of ! symmetric portion of a ! global field at one gaussian ! latitude. ! fm(Imax+2) input: fourier representation of ! anti-symmetric portion of a ! global field at one gaussian ! latitude. ! fln(Mnwv2) or input: spectral representation of ! fln(MnWv3) (the laplacian of) a global ! field. includes contributions ! from gaussian latitudes up to ! but not including current ! iteration of gaussian loop ! in calling routine. ! output: spectral representation of ! (the laplacian of) a global ! field. includes contributions ! from gaussian latitudes up to ! and including current ! iteration of gaussian loop ! in calling routine. ! lat input: current index of gaussian ! loop in calling routine. class(legendre), intent (in ) :: self real (kind=p_r8), intent (in ) :: fp(:) real (kind=p_r8), intent (in ) :: fm(:) integer, intent (in ) :: lat integer, intent (in ) :: MnWv real (kind=p_r8), intent (inout) :: fln(MnWv) integer :: k, l, nn, mmax, mm, mn real (kind=p_r8), dimension (self%mnwv2) :: s l=0 do nn=1,self%mend1 mmax=2*(self%mend2-nn) if (mod(nn-1,2) == 0) then do mm=1,mmax l=l+1 s(l)=fp(mm) end do else do mm=1,mmax l=l+1 s(l)=fm(mm) end do end if end do do mn=1,self%mnwv0 fln(2*mn-1)=fln(2*mn-1)+s(2*mn-1)*self%legS2F(mn,lat)*self%wgt(lat) fln(2*mn )=fln(2*mn )+s(2*mn )*self%legS2F(mn,lat)*self%wgt(lat) end do end subroutine Fourier2SpecCoef_ subroutine SpecCoef2Fourier_(self, fln, ap, am, lat) ! calculates the fourier representation of a field at a ! pair of latitudes symmetrically located about the ! equator. the calculation is made using the spectral ! representation of the field and the values of the ! associated legendre functions at that latitude. ! it is designed to triangular truncation only and ! for scalar fields. ! ! argument(dimensions) description ! ! fln(Mnwv2) input: spectral representation of a ! MnWv3 global field. ! ap(Imax+2) output: fourier representation of ! a global field at the ! latitude in the northern ! hemisphere at which the ! associated legendre functions ! have been defined. ! am(Imax+2) output: fourier representation of ! a global field at the ! latitude in the southern ! hemisphere at which the ! associated legendre functions ! have been defined. ! lat input: current index of gaussian ! loop in calling routine. character(len=17) :: myname_="**(Spec2Fourier_)**" class(legendre), intent (in ) :: self real (kind=p_r8), intent (in ) :: fln(:) integer, intent (in ) :: lat real (kind=p_r8), intent ( out) :: ap(:) real (kind=p_r8), intent ( out) :: am(:) integer :: l, k, mm, mn, nn, mstr, mmax1, mend1d integer :: nHarm integer :: MnWv integer :: nMax integer :: mMax real (kind=p_r8), allocatable :: s(:) real (kind=p_r8), pointer :: leg(:,:) => null() nHarm = size(fln) if(nHarm .eq. self%MnWv2)then nMax = self%Mend1 MnWv = self%MnWv0 mend1d = 2 * self%mend1 mmax1 = 2 * self%mend leg => self%legS2F allocate(s(self%MnWv2)) elseif (nHarm .eq. self%MnWv3)then nMax = self%Mend2 MnWv = self%MnWv1 mend1d = 2 * self%mend1 mmax1 = 2 * self%mend1 leg => self%legExtS2F allocate(s(self%MnWv3)) else write(unit=p_nferr, fmt='(2(A,1x),2(I6,1x,A,1x),I6)')& trim(myname_),& 'unknow size of spectral field ',mnwv,& 'should be',self%mnwv2,'or',self%mnwv3 stop 3022 endif l = mend1d + mmax1 do mn=1,MnWv s(2*mn-1) = leg(mn,lat)*fln(2*mn-1) s(2*mn ) = leg(mn,lat)*fln(2*mn ) end do do nn=3,nMax mMax = 2 * (nMax+1-nn) if (mod(nn-1,2) == 0) then mstr=0 else mstr=mend1d end if do mm=1,mMax l=l+1 s(mm+mstr)=s(mm+mstr)+s(l) end do end do do mm=1,mend1d ap(mm)=s(mm) am(mm)=s(mm) end do do mm=1,mmax1 ap(mm)=ap(mm)+s(mm+mend1d) am(mm)=am(mm)-s(mm+mend1d) end do deallocate(s) end subroutine subroutine SymAsy_ (self, a, b) ! converts the fourier representations of a field at two ! parallels at the same latitude in the northern and ! southern hemispheres into the fourier representations ! of the symmetric and anti-symmetric portions of that ! field at the same distance from the equator as the ! input latitude circles. ! ! argument(dimensions) description ! ! a(Imx,Ldim) input: fourier representation of one ! latitude circle of a field ! from the northern hemisphere ! at "n" levels in the vertical. ! output: fourier representation of the ! symmetric portion of a field ! at the same latitude as the ! input, at "n" levels in ! the vertical. ! b(Imx,Ldim) input: fourier representation of one ! latitude circle of a field ! from the southern hemisphere ! at "n" levels in the vertical. ! output: fourier representation of the ! anti-symmetric portion of a ! field at the same latitude as ! the input, at "n" levels in ! the vertical. ! t(Imx,Ldim) temporary storage ! class(legendre), intent (in ) :: self real (kind=p_r8), intent (inout) :: a(self%xMx) real (kind=p_r8), intent (inout) :: b(self%xMx) integer :: i, k real (kind=p_r8) :: t(self%xMx) do i=1,self%xMax ! rever essa dimensao, acredito que deveria se xMx t(i)=a(i) a(i)=t(i)+b(i) b(i)=t(i)-b(i) end do end subroutine SymAsy_ ! Eduardo Khamis subroutine getxMaxyMax (mend, xMax, yMax, linG) implicit none integer, intent(in) :: mend integer, intent(out) :: xMax, yMax logical, intent(in), optional :: linG ! Flag for linear (T) or quadratic (F) triangular truncation logical :: linearGrid integer :: nx, nm, n2m, n3m, n5m, n2, n3, n5, j, n, check, yfft integer, save :: lfft = 40000 integer, dimension(:), allocatable, save :: xfft if (present(linG)) then linearGrid = linG else linearGrid = .false. end if n2m = ceiling(log(real(lfft, p_r8)) / log(2.0_p_r8)) n3m = ceiling(log(real(lfft, p_r8)) / log(3.0_p_r8)) n5m = ceiling(log(real(lfft, p_r8)) / log(5.0_p_r8)) nx = n2m * (n3m + 1) * (n5m + 1) allocate(xfft (nx)) xfft = 0 n = 0 do n2 = 1, n2m yfft = (2**n2) if (yfft > lfft) exit do n3 = 0, n3m yfft = (2**n2) * (3**n3) if (yfft > lfft) exit do n5 = 0, n5m yfft = (2**n2) * (3**n3) * (5**n5) if (yfft > lfft) exit n = n + 1 xfft(n) = yfft end do end do end do nm = n n = 0 do check = 0 n = n + 1 do j = 1, nm - 1 if (xfft(j) > xfft(j + 1)) then yfft = xfft(j) xfft(j) = xfft(j + 1) xfft(j + 1) = yfft check = 1 end if end do if (check == 0) exit end do if (linearGrid) then yfft=2 else yfft=3 end if xMax = yfft * mend + 1 do n = 1, nm if (xfft(n) >= xMax) then xMax = xfft(n) exit end if end do yMax = xMax / 2 if (mod(yMax, 2) /= 0) yMax = yMax + 1 deallocate(xfft) end subroutine getxMaxyMax function GetSize(self,sname)result(sval) class(legendre) :: self character(len=*) :: sname integer :: sval ! ! Colocar aqui um lower case em sname! ! sname = lower(sname) ! select case(trim(sname)) case('mend') sval = self%Mend case('mend1') sval = self%Mend1 case('mend2') sval = self%Mend2 case('mnwv0') sval = self%MnWv0 case('mnwv1') sval = self%MnWv1 case('mnwv2') sval = self%MnWv2 case('mnwv3') sval = self%MnWv3 case('xmax') sval = self%xMax case('ymax') sval = self%yMax case('ymaxhf') sval = self%yMaxHf case default write(*,'(A,1x,A)')'unknown size var name:',trim(sname) sval = -1 end select end function subroutine DestroySpectralRep( self ) class(legendre), intent(inout) :: self deallocate( self%la0 ) deallocate( self%la1 ) deallocate( self%eps ) deallocate( self%e0 ) deallocate( self%e1 ) deallocate( self%snnp) end subroutine subroutine DestroyGaussRep( self ) class(legendre), intent(inout) :: self deallocate ( self%colRad ) deallocate ( self%rCs2 ) deallocate ( self%wgt ) deallocate ( self%gLats ) end subroutine subroutine DestroyLegPols( self ) class(legendre), intent(inout) :: self deallocate ( self%legS2F ) deallocate ( self%legDerS2F ) deallocate ( self%legExtS2F ) deallocate ( self%legDerNS ) deallocate ( self%legDerEW ) end subroutine subroutine print_size(label) !Cray-JNT! ! Esta subrotina reporta o Virtual Memory High Water Mark (VmHWM) ! presente no arquivo /proc/self/status do sistema Linux. ! character(len=256),intent(in) :: label integer :: iunit character(len=80) :: linha iunit = 7 open(iunit,file="/proc/self/status") do read(iunit,'(a)',END=999) linha if (INDEX(linha,'VmHWM').eq.1) then print *,trim(label), linha endif enddo 999 close(iunit) end subroutine print_size end module LegendreTransform