123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049 |
- ! NAME
- ! grid_type_sh - manipulate spherical harmonic field
- !
- ! may 2002, Arjo Segers
- !
- module grid_type_sh
- implicit none
-
- ! --- in/out ----------------------------
-
- private
-
- public :: TshGridInfo
- public :: TshGrid
- public :: Init, Done, Check
-
- public :: SpN
- public :: Set
- public :: Truncate
-
- public :: sh_Pnm
- public :: FourrierCoeff
- public :: Eval
- public :: EvalFast
- public :: Eval_Lons
- public :: Nabla
-
- public :: vod2uv
-
- ! --- const ------------------------------
-
- character(len=*), parameter :: mname = 'grid_type_sh'
-
- ! --- types -------------------------
-
- ! Coefficients of Associated Legendre Function
- ! * Triangular truncation: m=0,..,M=T n=m,..,N=T
- ! * Coefficients p(m,n), stored in order:
- ! (m,n) = (0,0) (0,1) .. (0,N) (1,1) .. (1,N) ..
- ! The storage order is similar to the order used in the EMOS library.
- ! Only coefficients for m>=0 are stored, since
- ! for real valued ALF the p(-m,n)=conj(p(m,n) .
-
- type TshGridInfo
- integer :: T ! triangular truncation
- integer :: np ! number of data points
- end type TshGridInfo
-
- type TshGrid
- integer :: T
- complex, pointer :: c(:)
- end type TshGrid
-
-
- ! --- interfaces ----------------------
-
- interface Init
- module procedure shi_Init
- module procedure shi_Init_copy
- module procedure shgrid_Init
- module procedure shgrid_Init_T
- module procedure shgrid_Init_1
- module procedure shgrid_Init_1_T
- end interface
-
- interface Done
- module procedure shi_Done
- module procedure shgrid_Done
- module procedure shgrid_Done_1
- end interface
-
- interface Check
- module procedure shgrid_Check
- module procedure shgrid_Check_T_1
- module procedure shgrid_Check_T
- end interface
-
- interface Set
- module procedure shi_Set_T
- module procedure shi_Set_shi
- module procedure Set_T
- module procedure Set_c
- module procedure Set_sh
- module procedure Set_T_sh
- end interface
-
- interface Truncate
- module procedure shgrid_Truncate
- end interface
-
- interface FourrierCoeff
- module procedure shgrid_FourrierCoeff_pnm
- module procedure shi_FourrierCoeff_pnm
- module procedure shgrid_FourrierCoeff_pnm2
- end interface
- interface Eval
- module procedure shgrid_Eval_Xm_lon
- module procedure shgrid_Eval_Pnm_lon
- module procedure shgrid_Eval_lat_lon
- end interface
-
- interface EvalFast
- module procedure shgrid_Eval_Pnm_lon_fast
- end interface
-
- interface Eval_Lons
- module procedure shi_Eval_Lons_Ulat
- module procedure shi_Eval_Lons_UPnm_1
- module procedure shgrid_Eval_Lons_Ulat
- module procedure shgrid_Eval_Lons_UPnm_1
- module procedure shgrid_Eval_Lons_Xm_1
- end interface
-
- interface Nabla
- module procedure shgrid_Nabla
- end interface
-
- interface vod2uv
- module procedure sh_vod2uv
- module procedure shi_vod2uv
- end interface
-
-
- contains
-
- ! ====================================================================
-
-
- subroutine shi_Init( shi, T, status )
-
- ! --- in/out ------------------------------------
-
- type(TshGridInfo), intent(out) :: shi
- integer, intent(in) :: T
- integer, intent(out) :: status
-
- ! --- begin ----------------------------------
- ! store triangular trunction:
- shi%T = T
-
- ! number of complex data points:
- shi%np = SpN( shi%T )
-
- ! ok
- status = 0
- end subroutine shi_Init
- ! ***
-
-
- subroutine shi_Init_copy( shi, shi2, status )
-
- ! --- in/out ------------------------------------
-
- type(TshGridInfo), intent(out) :: shi
- type(TshGridInfo), intent(in) :: shi2
- integer, intent(out) :: status
-
- ! --- begin ----------------------------------
- call Init( shi, shi2%T, status )
- end subroutine shi_Init_copy
- ! ***
-
- subroutine shi_Done( shi )
-
- ! --- in/out ------------------------------------
-
- type(TshGridInfo), intent(inout) :: shi
-
- ! --- begin ----------------------------------
- ! nothing to be done ...
- end subroutine shi_Done
- ! ***
-
-
- subroutine shi_Set_T( shi, T )
-
- ! --- in/out ------------------------------------
-
- type(TshGridInfo), intent(out) :: shi
- integer, intent(in) :: T
-
- ! --- begin ----------------------------------
- shi%T = T
- end subroutine shi_Set_T
- ! ***
-
-
- subroutine shi_Set_shi( shi, shi2 )
-
- ! --- in/out ------------------------------------
-
- type(TshGridInfo), intent(out) :: shi
- type(TshGridInfo), intent(in) :: shi2
-
- ! --- begin ----------------------------------
- shi%T = shi2%T
- end subroutine shi_Set_shi
- ! ====================================================================
-
-
- subroutine shgrid_Init( SH )
-
- ! --- in/out ------------------------------------
-
- type(TshGrid), intent(out) :: SH
-
- ! --- begin ----------------------------------
- ! ! check arguments ...
- ! if (associated(SH%c)) then
- ! print *, 'shgrid_Init : coefficient array already associated ...'
- ! stop
- ! end if
-
- SH%T = -1
- nullify( SH%c )
- end subroutine shgrid_Init
- ! ***
-
- subroutine shgrid_Init_T( SH, T )
-
- ! --- in/out ------------------------------------
-
- type(TshGrid), intent(out) :: SH
- integer, intent(in) :: T
-
- ! --- begin ----------------------------------
- ! ! check arguments ...
- ! if (associated(SH%c)) then
- ! print *, 'shgrid_Init : coefficient array already associated ...'
- ! stop
- ! end if
- ! set Triangular Truncation
- if ( T < 0 ) then
- print *, 'shgrid_Init : tried to set triangular truncation ', T
- stop
- end if
- SH%T = T
- ! allocate ...
- allocate( SH%c(SpN(T)) )
- ! set to zero:
- SH%c = (0.0,0.0)
-
- end subroutine shgrid_Init_T
- ! ***
-
- subroutine shgrid_Init_1( SH )
-
- ! --- in/out ------------------------------------
-
- type(TshGrid), intent(out) :: SH(:)
-
- ! --- local ----------------------------------
-
- integer :: l
-
- ! --- begin ----------------------------------
- do l = 1, size(SH)
-
- ! ! check arguments ...
- ! if ( associated(SH(l)%c) ) then
- ! print *, 'shgrid_Init : coefficient array already associated ...'
- ! stop
- ! end if
- SH(l)%T = -1
- nullify( SH(l)%c )
- end do
-
- end subroutine shgrid_Init_1
- ! ***
-
- subroutine shgrid_Init_1_T( SH, T )
-
- ! --- in/out ------------------------------------
-
- type(TshGrid), intent(out) :: SH(:)
- integer, intent(in) :: T
-
- ! --- local ----------------------------------
-
- integer :: l
-
- ! --- begin ----------------------------------
- do l = 1, size(SH)
-
- ! ! check arguments ...
- ! if (associated(SH(l)%c)) then
- ! print *, 'shgrid_Init : coefficient array already associated ...'
- ! stop
- ! end if
- ! set Triangular Truncation
- if ( T < 0 ) then
- print *, 'shgrid_Init : tried to set triangular truncation ', T
- stop
- end if
- SH(l)%T = T
- ! allocate ...
- allocate( SH(l)%c(SpN(T)) )
- ! set to zero:
- SH(l)%c = (0.0,0.0)
- end do
-
- end subroutine shgrid_Init_1_T
-
- ! ***
-
- subroutine shgrid_Done( SH )
-
- ! --- in/out ------------------------------------
-
- type(TshGrid), intent(out) :: SH
-
- ! --- begin ----------------------------------
-
- ! deallocate if necessary ...
- if (associated(SH%c)) deallocate(SH%c)
-
- end subroutine shgrid_Done
-
-
- ! ***
-
- subroutine shgrid_Done_1( SH )
-
- ! --- in/out ------------------------------------
-
- type(TshGrid), intent(out) :: SH(:)
-
- ! --- local ----------------------------------
-
- integer :: l
-
- ! --- begin ----------------------------------
- ! deallocate if necessary ...
- do l = 1, size(SH)
- if ( associated(SH(l)%c) ) deallocate( SH(l)%c )
- end do
-
- end subroutine shgrid_Done_1
-
-
- ! ====================================================================
- ! call Check( SH ) : check wether T and c are consistent
- ! call Check( SH, T ) : check wether SH is truncated at T
-
- subroutine shgrid_Check( SH )
-
- ! --- in/out ----------------------------
-
- type(TshGrid), intent(in) :: SH
-
- ! --- begin ----------------------------
-
- if (SH%T<0) then
- if (associated(SH%c)) then
- print *, 'shgrid_Check : error : coefficient array associated', &
- ' while T<0'
- stop
- end if
- else
- if (.not.associated(SH%c)) then
- print *, 'shgrid_Check : error : coefficient array not associated'
- stop
- end if
- if (size(SH%c) < SpN(SH%T) ) then
- print *, 'shgrid_Check : error : coefficient array has size ', size(SH%c), &
- 'while triangular truncation ', SH%T, &
- 'requires spectral number ', SpN(SH%T)
- stop
- end if
- end if
- end subroutine shgrid_Check
- ! ===
- subroutine shgrid_Check_T( SH, T )
-
- ! --- in/out ----------------------------
-
- type(TshGrid), intent(in) :: SH
- integer, intent(in) :: T
-
- ! --- begin ----------------------------
-
- if ( T < 0 ) then
- print *, 'shgrid_Check : tried to check for truncation ', T
- stop
- end if
- if ( SH%T /= T ) then
- print *, 'shgrid_Check : truncated at ', SH%T, ' instead of ', T
- stop
- end if
- if (.not.associated(SH%c)) then
- print *, 'shgrid_Check : error : coefficient array not associated'
- stop
- end if
- if (size(SH%c) < SpN(SH%T) ) then
- print *, 'shgrid_Check : error : coefficient array has size ', size(SH%c), &
- 'while triangular truncation ', SH%T, &
- 'requires spectral number ', SpN(SH%T)
- stop
- end if
- end subroutine shgrid_Check_T
- ! ===
- subroutine shgrid_Check_T_1( SH, T )
-
- ! --- in/out ----------------------------
-
- type(TshGrid), intent(in) :: SH(:)
- integer, intent(in) :: T
-
- ! --- local ---------------------------
-
- integer :: l
-
- ! --- begin ----------------------------
-
- do l = 1, size(SH)
- call Check( SH(l), T )
- end do
- end subroutine shgrid_Check_T_1
- ! ====================================================================
- ! Spectral number : SN = (T+1)*(T+2)/2
- ! (half triangle m=0,..,T , n=m,..,T ;
- ! we assume real restult thus only m>=0 required)
-
- integer function SpN( TT )
-
- ! --- in/out -----------------
-
- integer, intent(in) :: TT ! triangular truncation
-
- ! --- begin ------------------
-
- SpN = (TT+1)*(TT+2)/2
-
- end function SpN
- ! ====================================================================
-
- subroutine Set_T( SH, TT )
-
- ! --- in/out -------------------------
-
- type(TshGrid), intent(inout) :: SH
- integer, intent(in) :: TT
-
- ! --- local ---------------------------
-
- integer :: SN
-
- ! --- in/out -------------------------
- call check( SH )
-
- ! -- set triangular truncation
- ! check ...
- if ( TT < 0 ) then
- print *, 'Set : tried to set triangular truncation ', TT
- stop
- end if
- ! set:
- SH%T = TT
-
- ! -- set spectral coefficients
-
- ! number of coeff to be set:
- SN = SpN(SH%T)
- ! (de)allocate coefficient array if necessary ...
- if ( associated(SH%c) ) then
- if ( size(SH%c)/=SN ) deallocate(SH%c)
- end if
- if (.not. associated(SH%c) ) allocate(SH%c(SN))
- ! fill with zeros
- SH%c = (0.0,0.0)
-
- end subroutine Set_T
-
- ! ===
- subroutine Set_c( SH, TT, c )
-
- ! --- in/out -------------------------
-
- type(TshGrid), intent(inout) :: SH
- integer, intent(in) :: TT
- complex, intent(in) :: c(:)
-
- ! --- local ---------------------------
-
- integer :: SN
-
- ! --- in/out -------------------------
-
- ! init with zeros:
- call Set( SH, TT )
-
- ! number of coeff to be set:
- SN = SpN(SH%T)
- ! fill with c
- if ( size(c) /= SN ) then
- print *, 'Set_c : size of input array is ', size(c), &
- ' while triangular truncation ', TT, ' requires ', SN, &
- ' spectral coefficients'
- stop
- end if
- SH%c = c
-
- end subroutine Set_c
-
-
- ! *
-
-
- subroutine Set_T_sh( SH, TT, sh2 )
-
- ! --- in/out -------------------------
-
- type(TshGrid), intent(inout) :: SH
- integer, intent(in) :: TT
- type(TshGrid), intent(in) :: sh2
-
- ! --- in/out -------------------------
-
- ! init with full truncation:
- call Set( SH, sh2%T, sh2%c )
-
- ! truncate to requested number:
- call Truncate( SH, TT )
-
- end subroutine Set_T_sh
-
-
- ! *
-
-
- subroutine Set_sh( SH, sh2 )
-
- ! --- in/out -------------------------
-
- type(TshGrid), intent(inout) :: SH
- type(TshGrid), intent(in) :: sh2
-
- ! --- in/out -------------------------
-
- ! init with full truncation:
- call Set( SH, sh2%T, sh2%c )
-
- end subroutine Set_sh
-
-
- ! ====================================================================
-
- ! truncates spherical harmonic coeffiecients to lower
- ! triangular truncation.
-
- subroutine shgrid_Truncate( SH, T_new )
-
- ! --- in/out ------------------------------
- type(TshGrid), intent(inout) :: SH
- integer, intent(in) :: T_new
-
- ! --- local ----------------------------
-
- integer :: T, m, n
- integer :: k, k_new
- complex, allocatable :: c_new(:)
-
- ! ---- begin ----------------------------
-
- ! check arguments ...
- call Check( SH )
- if (T_new<0) then
- print *, 'shgrid_Truncate : error : tried to truncate to T=', T_new
- stop
- end if
- if (T_new > SH%T) then
- print *, 'shgrid_Truncate : error : tried to truncate to T=', T_new, &
- ' while input spherical harmonics are truncated at ', SH%T
- stop
- end if
-
- ! shorthands:
- T = SH%T
- ! truncate only if necessary ...
- if ( T_new < T ) then
- ! allocate temporary array for truncated coefficients:
- allocate( c_new(SpN(T_new)) )
- ! loop over original triangle for m <= T_new;
- k = 0
- k_new = 0
- do m = 0, T_new
- do n = m, T
- k = k + 1
- if (n <= T_new) then
- k_new = k_new + 1
- c_new(k_new) = SH%c(k)
- end if
- end do
- end do
-
- ! fill SH with truncated coefficients:
- call Set( SH, T_new, c_new )
-
- ! done
- deallocate( c_new )
-
- end if
-
- end subroutine shgrid_Truncate
-
- ! =======================================================
-
- !
- ! Evaluate associate Legendre functions at given latitude.
- ! Corner (T,T) is set to 0.0 .
- !
- ! Based on recurent formula:
- !
- !
- ! mu P(mu;m,n-1) = eps(m,n) P(mu;m,n) + eps(m,n-1) P(mu;m,n-2)
- !
- ! n^2 - m^2
- ! mu = sin(lat) , eps(m,n) = sqrt( --------- )
- ! 4n^2 - 1
- !
- ! P(m,n) = ( mu P(m,n-1) - eps(m,n-1) P(m,n-2) ) / eps(m,n)
- !
- !
- ! P(0,0) = 1
- !
- !
- ! P(0,1) = mu / eps(0,1) = mu sqrt(3)
- !
- !
- ! 2m+1 1/2 sqrt(1-mu^2)^m
- ! P(m,m) = ( ----- ) -------------- (2m)!
- ! (2m)! 2^m m!
- !
- ! sqrt(1-mu^2)^m
- ! = sqrt((2m+1)!) --------------
- ! 2^m m!
- !
- ! sqrt(1-mu^2)sqrt(1-mu^2)^(m-1)
- ! = sqrt((2m+1)2m(2(m-1))!) ------------------------------
- ! 2 2^(m-1) m (m-1)!
- !
- ! cos(lat)
- ! = P(m-1,m-1) sqrt((2m+1)2m) --------
- ! 2m
- !
- ! = P(m-1,m-1) sqrt((2m+1)) cos(lat) / sqrt(2m)
- !
- ! P(m,m+1) = mu P(m,m) / eps(m,m+1)
- !
- !
- subroutine sh_Pnm( Pnm, T, lat, status )
- ! --- in/out ----------------------------------
- real, intent(out) :: Pnm(:)
- integer, intent(in) :: T
- real, intent(in) :: lat ! rad
- integer, intent(out) :: status
- ! --- const --------------------------------------
-
- character(len=*), parameter :: rname = mname//'/sh_Pnm'
-
- ! --- local ------------------------------------
- integer :: m, n
- integer :: k
- real :: mu, rmu
- real :: fmm, fmmp
- real :: eps, eps1
-
- ! --- begin ---------------------------------------
- ! Use EMOS library:
- !Pnm = 0.0
- !call jspleg1( Pnm, lat, T-1 )
- !return
- if ( size(Pnm) /= SpN(T) ) then
- write (*,'("ERROR - wrong size of output array:")')
- write (*,'("ERROR - size(Pnm) : ",i6)') size(Pnm)
- write (*,'("ERROR - expected : ",i6)') SpN(T)
- write (*,'("ERROR - truncation : ",i6)') T
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
-
- mu = sin(lat)
- rmu = sqrt(1.0-mu*mu) ! cos(lat)
- do m = 0, T-1
- if ( m == 0 ) then
-
- fmmp = sqrt(3.0)
- k = 1; Pnm(k) = 1.0 ! (0,0)
- k = 2; Pnm(k) = fmmp*mu ! (0,1)
- else
- fmm = fmmp * rmu / sqrt( 2.0*m )
- fmmp = fmm * sqrt( 2*m + 3.0 )
- ! n = m and n = m+1
- k = k+1; Pnm(k) = fmm
- k = k+1; Pnm(k) = fmmp * mu
- endif
- ! leave if truncation is reached:
- if ( m == T-1 ) exit
- eps1 = 1.0 / sqrt( 2*m + 3.0 )
- do n = m+2, T
- eps = sqrt((n*n*1.0-m*m*1.0)/(4.0*n*n-1.0))
- k = k+1; Pnm(k) = ( mu*Pnm(k-1) - eps1*Pnm(k-2) ) / eps
- eps1 = eps
- end do
- end do
-
- ! set corner:
- k = k + 1; Pnm(k) = 0.0
-
- ! ok
- status = 0
-
- end subroutine sh_Pnm
- ! ***
-
-
- !! same, but without if statements in do loop
- !! (faster on some compilers)
- !
- !subroutine sh_Pnm_fast( Pnm, T, lat, status )
- !
- ! ! --- in/out ----------------------------------
- !
- ! real, intent(out) :: Pnm(:)
- ! integer, intent(in) :: T
- ! real, intent(in) :: lat ! rad
- ! integer, intent(out) :: status
- !
- ! ! --- const --------------------------------------
- !
- ! character(len=*), parameter :: rname = mname//'/sh_Pnm_fast'
- !
- ! ! --- local ------------------------------------
- !
- ! integer :: m, n
- ! integer :: k
- ! real :: mu, rmu
- ! real :: fmm, fmmp
- ! real :: eps, eps1
- !
- ! ! --- begin ---------------------------------------
- !
- ! ! Use EMOS library:
- ! !Pnm = 0.0
- ! !call jspleg1( Pnm, lat, T-1 )
- ! !return
- !
- ! if ( size(Pnm) /= SpN(T) ) then
- ! write (*,'("ERROR - wrong size of output array:")')
- ! write (*,'("ERROR - size(Pnm) : ",i6)') size(Pnm)
- ! write (*,'("ERROR - expected : ",i6)') SpN(T)
- ! write (*,'("ERROR - truncation : ",i6)') T
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end if
- !
- ! mu = sin(lat)
- ! rmu = sqrt(1.0-mu*mu) ! cos(lat)
- !
- ! ! -- first m
- !
- ! m = 0
- !
- ! fmmp = sqrt(3.0)
- ! k = 1; Pnm(k) = 1.0 ! (0,0)
- ! k = 2; Pnm(k) = fmmp*mu ! (0,1)
- !
- ! eps1 = 1.0 / sqrt( 2*m + 3.0 )
- !
- ! do n = m+2, T
- !
- ! eps = sqrt((n*n*1.0-m*m*1.0)/(4.0*n*n-1.0))
- ! k = k+1; Pnm(k) = ( mu*Pnm(k-1) - eps1*Pnm(k-2) ) / eps
- ! eps1 = eps
- !
- ! end do
- !
- ! ! -- 0 < m < T-1
- !
- ! do m = 1, T-2
- !
- ! fmm = fmmp * rmu / sqrt( 2.0*m )
- ! fmmp = fmm * sqrt( 2*m + 3.0 )
- !
- ! ! n = m and n = m+1
- !
- ! k = k+1; Pnm(k) = fmm
- ! k = k+1; Pnm(k) = fmmp * mu
- !
- ! eps1 = 1.0 / sqrt( 2*m + 3.0 )
- !
- ! do n = m+2, T
- !
- ! eps = sqrt((n*n*1.0-m*m*1.0)/(4.0*n*n-1.0))
- ! k = k+1; Pnm(k) = ( mu*Pnm(k-1) - eps1*Pnm(k-2) ) / eps
- ! eps1 = eps
- !
- ! end do
- !
- ! end do
- !
- ! ! -- m = T-1
- !
- ! m = T-1
- !
- ! fmm = fmmp * rmu / sqrt( 2.0*m )
- ! fmmp = fmm * sqrt( 2*m + 3.0 )
- !
- ! ! n = m and n = m+1
- !
- ! k = k+1; Pnm(k) = fmm
- ! k = k+1; Pnm(k) = fmmp * mu
- !
- ! eps1 = 1.0 / sqrt( 2*m + 3.0 )
- !
- ! do n = m+2, T
- !
- ! eps = sqrt((n*n*1.0-m*m*1.0)/(4.0*n*n-1.0))
- ! k = k+1; Pnm(k) = ( mu*Pnm(k-1) - eps1*Pnm(k-2) ) / eps
- ! eps1 = eps
- !
- ! end do
- !
- ! ! set corner:
- ! k = k + 1; Pnm(k) = 0.0
- !
- !
- !end subroutine sh_Pnm_fast
- ! ====================================================================
- ! M N i m lon
- ! Eval(p,lon,lat) = sum sum p(m,n) P(m,n,sin(lat)) e
- ! m=-M n=m
- !
- ! M
- ! = c(0) + sum 2[Re{c(m)}cos(m lon) - Im{c(m)}sin(m lon)] ,
- ! m=1
- !
- ! N
- ! c(m) = sum p(m,n) P(m,n,sin(lat))
- ! n=m
- !
-
- real function shgrid_Eval_Xm_lon( Xm, T, lon )
-
- ! --- in/out -------------------------
-
- integer, intent(in) :: T
- complex, intent(in) :: Xm(0:T)
- real, intent(in) :: lon
-
- ! --- local --------------------------
-
- integer :: m
- ! --- begin --------------------------
-
- ! summation over triangle
- shgrid_eval_Xm_lon = Xm(0)
- do m = 1, T
- shgrid_eval_Xm_lon = shgrid_eval_Xm_lon + 2.0*( real(Xm(m))*cos(m*lon) - aimag(Xm(m))*sin(m*lon) )
- end do
-
- end function shgrid_Eval_Xm_lon
-
-
- ! ***
-
-
- real function shgrid_Eval_Pnm_lon( SH, Pnm, lon )
-
- ! --- in/out -------------------------
-
- type(TshGrid), intent(in) :: SH
- real, intent(in) :: Pnm(:)
- real, intent(in) :: lon
-
- ! --- local --------------------------
-
- integer :: T, SN
- integer :: m, n, k
- complex :: cm
- ! --- begin --------------------------
-
- ! shorthands:
- T = SH%T
- SN = SpN(T)
-
- ! check arguments:
- call Check( SH )
- if ( size(Pnm) /= SN ) then
- print *, 'shgrid_FourrierCoeff_pnm : size(Pnm)=',size(Pnm), &
- 'while',SN,'expected for T',T
- stop
- end if
- ! summation over triangle
- k = 0
- do m = 0, T
- cm = (0.0,0.0)
- do n = m, T
- k = k + 1
- cm = cm + SH%c(k) * Pnm(k)
- end do
- if (m==0) then
- shgrid_eval_Pnm_lon = real(cm)
- else
- shgrid_eval_Pnm_lon = shgrid_eval_Pnm_lon + 2.0*( real(cm)*cos(m*lon) - aimag(cm)*sin(m*lon) )
- end if
- end do
- end function shgrid_Eval_Pnm_lon
-
-
- ! ***
-
-
- real function shgrid_Eval_Pnm_lon_fast( SH, Pnm, lon )
-
- ! --- in/out -------------------------
-
- type(TshGrid), intent(in) :: SH
- real, intent(in) :: Pnm(:)
- real, intent(in) :: lon
-
- ! --- local --------------------------
-
- integer :: T, SN
- integer :: m, n, k
- complex :: cm
- complex :: cp(size(Pnm))
- real :: res
- ! --- begin --------------------------
-
- ! shorthands:
- T = SH%T
- SN = SpN(T)
-
- ! check arguments:
- call Check( SH )
- if ( size(Pnm) /= SN ) then
- print *, 'shgrid_FourrierCoeff_pnm : size(Pnm)=',size(Pnm), &
- 'while',SN,'expected for T',T
- stop
- end if
- cp = SH%c * Pnm
- ! summation over triangle
- res = sum( cp(1:T+1) )
- k = T+2
- do m = 1, T
- cm = sum( cp(k:k+T-m) )
- k = k + T-m+1
- res = res + 2.0*( real(cm)*cos(m*lon) - aimag(cm)*sin(m*lon) )
- end do
-
- shgrid_eval_Pnm_lon_fast = res
-
- end function shgrid_Eval_Pnm_lon_fast
-
-
- ! ***
-
-
- real function shgrid_Eval_lat_lon( SH, lat, lon, status )
-
- ! --- in/out -------------------------
-
- type(TshGrid), intent(in) :: SH
- real, intent(in) :: lon, lat
- integer, intent(out) :: status
-
- ! --- const --------------------------------------
-
- character(len=*), parameter :: rname = mname//'/shgrid_Eval_lat_lon'
-
- ! --- local --------------------------
-
- integer :: T
- real, allocatable :: Pnm(:)
- ! --- begin --------------------------
-
- ! check arguments:
- call Check( SH )
- ! shorthands:
- T = SH%T
-
- ! allocate array for associated Legendre functions:
- allocate( Pnm(SpN(T)) )
- ! evaluate Legendre functions at given latitude:
- call sh_Pnm( Pnm, T, lat, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- ! expand summation:
- shgrid_Eval_lat_lon = Eval( SH, Pnm, lon )
- ! done
- deallocate( Pnm )
-
- ! ok
- status = 0
-
- end function shgrid_Eval_lat_lon
-
-
-
-
- ! ==============================================
-
- ! compute Fourrier coefficients:
- !
- ! T
- ! X(m;mu) = sum U(m,n) P(m,n;mu) = Xe(m;mu) + Xo(m;mu)
- ! n=|m|
- !
- ! T T
- ! = sum U(m,n) P(m,n;mu) + sum U(m,n) P(m,n;mu)
- ! n=|m| n=|m|
- ! (n-m) even (n-m) odd
- !
- ! and return coefficients for both +mu and -mu using that
- !
- ! X(m,1) = X(m;+mu) = Xe(m;mu) + Xo(m;mu) 'north'
- !
- ! X(m,2) = X(m;-mu) = Xe(m;mu) - Xo(m;mu) 'south'
- !
- ! (valid for real symmetric Pmn(mu) )
-
- ! subroutine shgrid_FourrierCoeff( Xnm, Pnm, T, X )
- !
- ! ! --- in/out -----------------------------------
- !
- ! type(TshGrid), intent(in) :: Xnm, Pnm
- ! integer, intent(in) :: T
- ! complex, intent(out) :: X(0:T,2)
- !
- ! ! --- local ---------------------------------
- !
- ! complex :: Xe, Xo
- ! integer :: m, n, k
- !
- ! ! --- begin ----------------------------------
- !
- ! ! check input ...
- ! call Check( Xnm, T )
- ! call Check( Pnm, T )
- !
- ! ! loop over all coeff:
- ! m = 0
- ! n = 0
- ! Xe = (0.0,0.0)
- ! Xo = (0.0,0.0)
- ! do k = 1, SpN(T)
- ! if ( mod(n+m,2)==0 ) then
- ! Xe = Xe + Xnm%c(k)*Pnm%c(k)
- ! else
- ! Xo = Xo + Xnm%c(k)*Pnm%c(k)
- ! end if
- ! ! next n
- ! n = n + 1
- ! ! end of current m ?
- ! if ( n > T ) then
- ! ! save fourrier coeff for current m :
- ! X(m,1) = Xe + Xo
- ! X(m,2) = Xe - Xo
- ! ! next m
- ! m = m + 1
- ! n = m
- ! Xe = (0.0,0.0)
- ! Xo = (0.0,0.0)
- ! end if
- ! end do
- !
- ! ! adhoc: ensure that X(0) is real:
- ! X(0,:) = real(X(0,:)) * (1.0,0.0)
- !
- ! end subroutine shgrid_FourrierCoeff
-
-
-
-
- subroutine shgrid_FourrierCoeff_pnm( Xnm, Pnm, T, X )
-
- ! --- in/out -----------------------------------
-
- type(TshGrid), intent(in) :: Xnm
- real, intent(in) :: Pnm(:)
- integer, intent(in) :: T
- complex, intent(out) :: X(0:T)
-
- ! --- local ---------------------------------
-
- complex :: Xe, Xo
- integer :: m, n, k
-
- ! --- begin ----------------------------------
-
- ! check input ...
- call Check( Xnm, T )
- if ( size(Pnm) /= SpN(T) ) then
- print *, 'shgrid_FourrierCoeff_pnm : size(Pnm)=',size(Pnm), &
- 'while',SpN(T),'expected for T',T
- stop
- end if
-
- ! loop over all coeff:
- k = 0
- do m = 0, T
- X(m) = (0.0,0.0)
- do n = m, T
- k = k + 1
- X(m) = X(m) + Xnm%c(k) * Pnm(k)
- end do
- end do
-
- ! adhoc: ensure that X(0) is real:
- X(0) = real(X(0)) * (1.0,0.0)
-
- end subroutine shgrid_FourrierCoeff_pnm
-
-
- ! ***
-
-
- pure subroutine shi_FourrierCoeff_pnm( shi_in, Xnm, shi, Pnm, X )
-
- ! --- in/out -----------------------------------
-
- type(TshGridInfo), intent(in) :: shi_in
- complex, intent(in) :: Xnm(1:shi_in%np)
- type(TshGridInfo), intent(in) :: shi
- real, intent(in) :: Pnm(1:shi%np)
- complex, intent(out) :: X(0:shi%T)
-
- ! --- const --------------------------------------
-
- character(len=*), parameter :: rname = mname//'/shi_FourrierCoeff_pnm'
-
- ! --- local ---------------------------------
-
- integer :: m, n, kx, kp
-
- ! --- begin ----------------------------------
- #ifdef check_all
- ! check input ...
- if ( shi_in%T < shi%T ) then
- write (*,'("ERROR - input truncation should not be less than output:")')
- write (*,'("ERROR - shi_in%T : ",i6)') shi_in%T
- write (*,'("ERROR - shi%T : ",i6)') shi%T
- write (*,'("ERROR in ",a)') rname; stop
- end if
- if ( size(Xnm) /= shi_in%np ) then
- write (*,'("ERROR - mismatch between size of Xnm and specified truncation:")')
- write (*,'("ERROR - shi_in%T : ",i6)') shi_in%T
- write (*,'("ERROR - shi_in%np : ",i6)') shi_in%np
- write (*,'("ERROR - Xnm : ",i6)') size(Xnm)
- write (*,'("ERROR in ",a)') rname; stop
- end if
- if ( size(Pnm) /= shi%np ) then
- write (*,'("ERROR - mismatch between size of Pnm and specified truncation:")')
- write (*,'("ERROR - shi%T : ",i6)') shi%T
- write (*,'("ERROR - shi%np : ",i6)') shi%np
- write (*,'("ERROR - Pnm : ",i6)') size(Pnm)
- write (*,'("ERROR in ",a)') rname; stop
- end if
- #endif
-
- ! loop over all coeff:
- kx = 0
- kp = 0
- do m = 0, shi%T
- X(m) = (0.0,0.0)
- do n = m, shi%T
- kx = kx + 1
- kp = kp + 1
- X(m) = X(m) + Xnm(kx) * Pnm(kp)
- end do
- kx = kx + shi_in%T - shi%T
- end do
-
- ! adhoc: ensure that X(0) is real:
- X(0) = real(X(0)) * (1.0,0.0)
-
- end subroutine shi_FourrierCoeff_pnm
-
-
- ! ***
-
-
- subroutine shgrid_FourrierCoeff_pnm2( Xnm, Pnm, T, X )
-
- ! --- in/out -----------------------------------
-
- type(TshGrid), intent(in) :: Xnm
- real, intent(in) :: Pnm(:)
- integer, intent(in) :: T
- complex, intent(out) :: X(0:T,2)
-
- ! --- local ---------------------------------
-
- complex :: Xe, Xo
- integer :: m, n, k
-
- ! --- begin ----------------------------------
-
- ! check input ...
- call Check( Xnm, T )
- if ( size(Pnm) /= SpN(T) ) then
- print *, 'shgrid_FourrierCoeff_pnm : size(Pnm)=',size(Pnm), &
- 'while',SpN(T),'expected for T',T
- stop
- end if
-
- ! loop over all coeff:
- m = 0
- n = 0
- Xe = (0.0,0.0)
- Xo = (0.0,0.0)
- do k = 1, SpN(T)
- if ( mod(n+m,2)==0 ) then
- Xe = Xe + Xnm%c(k)*Pnm(k)
- else
- Xo = Xo + Xnm%c(k)*Pnm(k)
- end if
- ! next n
- n = n + 1
- ! end of current m ?
- if ( n > T ) then
- ! save fourrier coeff for current m :
- X(m,1) = Xe + Xo
- X(m,2) = Xe - Xo
- ! next m
- m = m + 1
- n = m
- Xe = (0.0,0.0)
- Xo = (0.0,0.0)
- end if
- end do
-
- ! adhoc: ensure that X(0) is real:
- X(0,:) = real(X(0,:)) * (1.0,0.0)
-
- end subroutine shgrid_FourrierCoeff_pnm2
-
-
- ! ====
-
- !
- ! Evaluate x(lat) and x(-lat) given:
- ! Xnm : complex Associated Legendre Func. coefficients
- ! Pnm : Associated Legendre Functions evaluated at mu=sin(lat)
- ! K : number of lon values on a complete latitude circle;
- ! distance between to lon values is thus 2pi/K;
- ! lon0 : western boundary of lon grid (rad)
- ! Kout : actual number of lon points from lon0 to east .
- !
- ! Evaluation at oposite latitdues at the same time saves some
- ! computation time.
- !
- ! x(t0+k2pi/K) , k=0,1,..,K0-1
- !
- ! T i m [t0 + k 2pi/K]
- ! = sum X(m) e , k=0,1,..,K0-1
- ! m=-T
- !
- ! L such that KL >= (2T+1)
- ! k=(j+1)/L , j=kL-1 , J=KL
- !
- ! T i m t0 i m (j+1)2pi/KL
- ! = sum X(m) e e , j=-1,L-1,2L-1,...,(K0-1)L-1
- ! m=-T
- !
- ! T im[t0+2pi/J] i m j 2pi/J
- ! = sum [X(m)e ] e , j=-1,L-1,2L-1,...,(K0-1)L-1
- ! m=-T
- !
- ! _
- ! X(m) = X(m) exp(i phi) , phi=m[t0+2pi/J]
- ! = (a+ib)(cos(phi)+isin(phi))
- ! = [ a cos(phi) - b sin(phi) ] + i [ a sin(phi) + b cos(phi) ]
- !
- ! FFT99 uses that X(-m) = conjg(X(m)), and thus requires X(m) for m>=0 only;
- ! higher fourrier coefficients are zero:
- ! X(0) X(1) ... X(T) 0 0 0 ...
- ! |--- 1+floor(J/2) elements ---|
- !
- ! o o x x x 0 0 ... 0 0 x x
- ! m = 0 1 2 J-2 J-1
- !
- ! Number of K longitudes on complete latitude circle,
- ! while Kout evaluations are required.
- ! Note that Kout might be both smaller and larger than K.
- ! Routine FFT99 returns J+2 evaluations:
- ! x(j=-1) x(j=0) ... x(j=J-1=-1) x(j=J=0)
- ! thus if Kout exceeds K+2 then the result should be taken
- ! cyclic from the output of FFT99.
- !
-
- ! Evaluate given spherical harmonic coefficients and latitude
- subroutine shi_Eval_Lons_Ulat( shi, Xnm, lat, KK, lon_start, Kout, llgrid, status )
- ! --- in/out -------------------------
-
- type(TshGridInfo), intent(in) :: shi
- complex, intent(in) :: Xnm(1:shi%np)
- real, intent(in) :: lat
- integer, intent(in) :: KK
- real, intent(in) :: lon_start
- integer, intent(in) :: Kout
- real,intent(out) :: llgrid(Kout)
- integer, intent(out) :: status
-
- ! --- const --------------------------------------
-
- character(len=*), parameter :: rname = mname//'/shi_Eval_Lons_Ulat'
-
- ! --- local --------------------------
- complex :: X(0:shi%T)
- real :: Pnm(1:shi%np)
-
- ! --- begin --------------------------
-
- ! evaluate Legendre functions at given latitude:
- call sh_Pnm( Pnm, shi%T, lat, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- ! compute Fourrier coefficients
- ! (latitude is implicit defined in Pnm)
- call FourrierCoeff( shi, Xnm, shi, Pnm, X )
- ! evaluate at longitudes given Fourrier coefficients:
- call Eval_Lons( llgrid, X, shi%T, KK, lon_start, Kout )
-
- ! ok
- status = 0
-
- end subroutine shi_Eval_Lons_Ulat
-
- ! *
-
- subroutine shgrid_Eval_Lons_Ulat( llgrid, Xnm, lat, KK, lon_start, Kout, status )
- ! --- in/out -------------------------
-
- type(TshGrid), intent(in) :: Xnm
- real, intent(in) :: lat
- integer, intent(in) :: KK
- real, intent(in) :: lon_start
- integer, intent(in) :: Kout
- real,intent(out) :: llgrid(Kout)
- integer, intent(out) :: status
-
- ! --- const --------------------------------------
-
- character(len=*), parameter :: rname = mname//'/shgrid_Eval_Lons_Ulat'
-
- ! --- local --------------------------
- integer :: T
- complex, allocatable :: X(:,:)
- real, allocatable :: Pnm(:)
-
- ! --- begin --------------------------
-
- ! timing:
- ! spleg1 65 %
- ! fourrier 25 %
- ! eval 10 %
-
- ! extract triangular truncation
- T = Xnm%T
-
- ! allocate array for associated Legendre functions:
- allocate( Pnm(SpN(T)) )
- ! evaluate Legendre functions at given latitude:
- call sh_Pnm( Pnm, T, lat, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- ! compute Fourrier coefficients
- ! (latitude is implicit defined in Pnm)
- allocate( X(0:T,2) )
- call FourrierCoeff( Xnm, Pnm, T, X )
- ! evaluate at longitudes given Fourrier coefficients:
- call Eval_Lons( llgrid, X(:,1), T, KK, lon_start, Kout )
-
- ! done
- deallocate( X )
- deallocate( Pnm )
-
- ! ok
- status = 0
-
- end subroutine shgrid_Eval_Lons_Ulat
-
- ! ***
-
- ! Evaluate given spherical harmonic coefficients and Legendre functions
- subroutine shgrid_Eval_Lons_UPnm_1( llgrid, Xnm, Pnm, KK, lon_start, Kout, status )
- ! --- in/out -------------------------
-
- type(TshGrid), intent(in) :: Xnm
- real, intent(in) :: Pnm(:)
- integer, intent(in) :: KK
- real, intent(in) :: lon_start
- integer, intent(in) :: Kout
- real,intent(out) :: llgrid(Kout)
- integer, intent(out) :: status
-
- ! --- const --------------------------------------
-
- character(len=*), parameter :: rname = mname//'/shgrid_Eval_Lons_UPnm_1'
-
- ! --- local --------------------------
- integer :: T
- complex, allocatable :: X(:)
-
- ! --- begin --------------------------
-
- T = Xnm%T
-
- ! compute Fourrier coefficients at oposite latitudes
- ! (latitude is implicit defined in Pnm)
- allocate( X(0:T) )
- call FourrierCoeff( Xnm, Pnm, T, X )
- ! evaluate given Fourrier coefficients:
- call Eval_Lons( llgrid, X, T, KK, lon_start, Kout )
-
- ! done
- deallocate( X )
-
- ! ok
- status = 0
-
- end subroutine shgrid_Eval_Lons_UPnm_1
-
- ! ***
-
- ! Evaluate given spherical harmonic coefficients and Legendre functions
- ! o TXnm, Xnm : input truncation and sh coeff
- ! o T, Pnm : target truncation and evaluated Legende functions
- ! o KK : number of longitudes on full circle
- ! o lon_start : first longitude
- ! o Kout : number of output longitudes
- ! o llgrid : output array
- subroutine shi_Eval_Lons_UPnm_1( shi_in, Xnm, shi, Pnm, KK, lon_start, Kout, llgrid, status )
- ! --- in/out -------------------------
-
- type(TShGridInfo), intent(in) :: shi_in
- complex, intent(in) :: Xnm(shi_in%np)
- type(TShGridInfo), intent(in) :: shi
- real, intent(in) :: Pnm(shi%np)
- integer, intent(in) :: KK
- real, intent(in) :: lon_start
- integer, intent(in) :: Kout
- real,intent(out) :: llgrid(Kout)
- integer, intent(out) :: status
-
- ! --- const --------------------------------------
-
- character(len=*), parameter :: rname = mname//'/sh_Eval_Lons_UPnm_1'
-
- ! --- local --------------------------
- complex :: X(0:shi%T)
-
- ! --- begin --------------------------
- ! compute Fourrier coefficients at oposite latitudes
- ! (latitude is implicit defined in Pnm)
- call FourrierCoeff( shi_in, Xnm, shi, Pnm, X )
- ! evaluate given Fourrier coefficients:
- call Eval_Lons( llgrid, X, shi%T, KK, lon_start, Kout )
-
- ! ok
- status = 0
-
- end subroutine shi_Eval_Lons_UPnm_1
-
- ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
-
- ! Shift to first lon:
- !
- ! x(t0+k2pi/K) , k=0,1,..,KK-1
- !
- ! T i m [t0 + k 2pi/K]
- ! = sum X(m) e , k=0,1,..,K0-1
- ! m=-T
- !
- ! T i m t0 k 2pi/K
- ! = sum [X(m) e ] e , k=0,1,..,K0-1
- ! m=-T
- !
- ! PERFORMANCE
- ! This routine is about 7% slower than with emos fft ...
- !
- ! real(4), external :: etime
- ! real :: tarr(2), t0
- ! t0=etime(tarr)
- ! print *, 'time:',etime(tarr)-t0
-
-
- subroutine shgrid_Eval_Lons_Xm_1( llgrid, Xm, T, KK, lon_start, Kout )
- use grid_tools, only : pi
- use grid_singleton, only : fft
-
- ! --- in/out -------------------------
-
- integer, intent(in) :: T
- complex, intent(in) :: Xm(0:T)
- integer, intent(in) :: KK
- real, intent(in) :: lon_start
- integer, intent(in) :: Kout
- real,intent(out) :: llgrid(Kout)
-
- ! --- local --------------------------
- integer :: i, j, k
- integer :: m, n
-
- complex :: X(0:T)
- integer :: LL, JJ
- complex, allocatable :: C(:), F(:)
- real :: fac
- ! --- begin --------------------------
-
- ! shift over west most longitude:
- ! X(m) := X(m) exp( i m lon0 )
- do m = 0, T
- X(m) = Xm(m) * exp( (0.0,1.0)*m*lon_start )
- end do
- ! choose L such that JJ=KK*LL >= 2T+1
- LL = 1
- do
- if ( KK*LL >= 2*T+1 ) exit
- LL = LL * 2
- end do
-
- ! number of longitudes to be evaluated:
- JJ = KK*LL
- ! fft arrays:
- allocate( C(JJ) )
- allocate( F(JJ) )
-
- ! fill coeff array:
- C = (0.0,0.0)
- C(1) = X(0)
- do m = 1, T
- C(1+m) = X(m)
- C(JJ+1-m) = conjg(X(m))
- end do
- ! Apply fast fourrier transform.
- !
- ! 1 J-1 i m lambda(j)
- ! F(j) = ------- sum C(m) e
- ! sqrt(J) m=0
- !
- ! where
- !
- ! lambda(j) = j 2pi/J = 0, dlon, 2*dlon, ...
- ! ____ ____
- ! C = [ X(0), X(1), .., X(T), 0, ..., 0, X(T), .., X(1) ]
- !
- ! Coeff in C are complex conj, thus imag part of F is zero.returned array is real.
- !
- F = fft( C, inv=.true. )
- fac = sqrt(JJ*1.0)
- ! Extract result.
- ! Since 'only' KK lons are evaluated (covering complete circle),
- ! select values cyclic for k>KK.
- ! Select first each group of LL elements
- do k = 0, Kout-1
- j = mod(k,KK)*LL + 1 ! 1, LL+1, 2*LL+1, ...
- llgrid(k+1) = real(F(j))*fac
- end do
-
- ! done
- deallocate( C )
- deallocate( F )
-
- end subroutine shgrid_Eval_Lons_Xm_1
-
-
- ! ==============================================
-
-
- ! compute associated legendre coeff for
- !
- ! d X 2 d X dX dX
- ! ( -------- , (1-mu ) ---- ) = r cos(theta) ( -- , -- )
- ! d lambda d mu dx dy
- !
- ! given associated legendre coeff of X(lambda,mu)
- !
-
- subroutine shgrid_Nabla( Xnm, NablaX )
- ! --- in/out -----------------------------
-
- type(TshGrid), intent(in) :: Xnm
- type(TshGrid), intent(inout) :: NablaX(2)
-
- ! --- local ------------------------
-
- integer :: m, n, T
- integer :: k
-
- ! --- begin ----------------------------
- ! extract triangular truncation
- T = Xnm%T
-
- ! allocate output arrays
- call Set( NablaX(1), T )
- call Set( NablaX(2), T )
-
- ! loop over all coeff
- k = 0
- do m = 0, T
- do n = m, T
- k = k + 1
- ! d(ln ps) i m lon
- ! -------- = sum {i m X(m,n)} P(mu;m,n) e
- ! d lon m,n
- NablaX(1)%c(k) = (0.0,1.0) * m * Xnm%c(k)
- ! 2 d(ln ps) ~ i m lon
- ! (1-mu ) -------- = sum X(m,n) P(mu;m,n) e
- ! d mu m,n
- !
- ! ~ = (n+2) eps(m,n+1) X(m,n+1) , m<T, m=n
- ! X(m,n) = (n+2) eps(m,n+1) X(m,n+1) - (n-1) eps(m,n) X(m,n-1) , m<T, m<n<T
- ! = - (n-1) eps(m,n) X(m,n-1) , m<T, n=T
- ! = 0.0 , m=T, n=T
-
- if ( m < T ) then
- if ( n == m ) then
- NablaX(2)%c(k) = (n+2) * epsilon(m,n+1) * Xnm%c(k+1)
- else if ( n > m .and. n < T ) then
- NablaX(2)%c(k) = (n+2) * epsilon(m,n+1) * Xnm%c(k+1) - &
- (n-1) * epsilon(m,n) * Xnm%c(k-1)
- else ! n == T
- NablaX(2)%c(k) = (n-1) * epsilon(m,n) * Xnm%c(k-1)
- end if
- else
- NablaX(2)%c(k) = 0.0
- end if
- end do
- end do
-
- contains
- ! n*n - m*m 1/2
- ! epsilon = ( --------- )
- ! 4*n*n - 1
- real function epsilon( m, n )
- ! --- in/out --------------------------------
- integer, intent(in) :: m, n
- ! --- begin --------------------------------
- epsilon = sqrt((n*n-m*m)*1.0/(4*n*n-1.0))
- end function epsilon
-
- end subroutine shgrid_Nabla
- ! ==================================================================
-
- !
- ! Compute (U,V) = (u,v) cos(lat) from vorticity VO and divergence D.
- !
- ! Use relations:
- !
- ! U(m,n) = ( + delta(m,n) VO(m,n-1) + i sigma(m,n) D(m,n) - delta(m,n+1) VO(m,n+1) ) R
- !
- ! V(m,n) = ( - delta(m,n) D(m,n-1) + i sigma(m,n) VO(m,n) + delta(m,n+1) D(m,n+1) ) R
- !
- ! where
- !
- ! delta(m,n) = - eps(m,n)/n , eps(m,n) = sqrt( (n^2-m^2)/(4n^2-1) )
- !
- ! sigma(m,n) = - m / (n(n+1))
- !
- ! NOTE: In EMOS library, vor. and div. are assumed to be truncated at T-1;
- ! for consitency, this is used here too ...
- !
-
- subroutine sh_vod2uv( VO, DI, U, V )
-
- use Binas, only : R => ae
- ! --- in/out -----------------------------
- type(TshGrid), intent(in) :: VO, DI
- type(TshGrid), intent(inout) :: U, V
- ! --- local ------------------------------
- integer :: m, n, T
- integer :: k
- complex :: z
- ! --- begin ----------------------------
- ! Use EMOS library:
- !call jvod2uv( VO%c, DI%c, VO%T, U%c, V%c, VO%T )
- ! check arguments:
- call Check( VO )
- T = VO%T
- call Check( DI, T )
- call Check( U, T )
- call Check( V, T )
- ! complex i for multiplication:
- z = (0.0, 1.0)
- ! index of coeff in array; initialise:
- k = 0
-
- ! loop over all m:
- do m = 0, T-1
- ! (m,m)
- n = m
- k = k + 1
- if ( m == 0 ) then
- U%c(k) = ( - delta(m,n+1)*VO%c(k+1) )*R
- V%c(k) = ( delta(m,n+1)*DI%c(k+1) )*R
- else
- U%c(k) = ( z*sigma(m,n)*DI%c(k) - delta(m,n+1)*VO%c(k+1) )*R
- V%c(k) = ( z*sigma(m,n)*VO%c(k) + delta(m,n+1)*DI%c(k+1) )*R
- end if
- ! internal area:
- if ( m <= T-2 ) then
- do n = m+1, T-2
- k = k + 1
- U%c(k) = ( delta(m,n)*VO%c(k-1) + z*sigma(m,n)*DI%c(k) - delta(m,n+1)*VO%c(k+1) )*R
- V%c(k) = ( - delta(m,n)*DI%c(k-1) + z*sigma(m,n)*VO%c(k) + delta(m,n+1)*DI%c(k+1) )*R
- end do
- ! (m,T-1)
- n = T-1
- k = k + 1
- U%c(k) = ( delta(m,n)*VO%c(k-1) + z*sigma(m,n)*DI%c(k))*R
- V%c(k) = ( - delta(m,n)*DI%c(k-1) + z*sigma(m,n)*VO%c(k))*R
- end if
- ! (m,T)
- n = T
- k = k + 1
- U%c(k) = ( delta(m,n)*VO%c(k-1) )*R
- V%c(k) = ( - delta(m,n)*DI%c(k-1) )*R
- end do
- ! (T,T)
- k = k + 1
- U%c(k) = 0.0
- V%c(k) = 0.0
- contains
- real function delta( m, n )
- integer, intent(in) :: m, n
- delta = - sqrt( (n*n-m*m)*1.0/(4.0*n*n-1.0) ) / (n*1.0)
- end function delta
- real function sigma( m, n )
- integer, intent(in) :: m, n
- sigma = - m * 1.0/(n*(n+1.0))
- end function sigma
- end subroutine sh_vod2uv
-
- ! ***
-
-
- pure subroutine shi_vod2uv( shi_in, VO, DI, shi, U, V )
-
- use Binas, only : R => ae
- ! --- in/out -----------------------------
- type(TShGridInfo), intent(in) :: shi_in
- complex, intent(in) :: VO(shi_in%np), DI(shi_in%np)
- type(TShGridInfo), intent(in) :: shi
- complex, intent(out) :: U(shi%np), V(shi%np)
- ! --- local ------------------------------
- integer :: m, n
- integer :: kin, k
- complex :: z
- ! --- begin ----------------------------
- ! Use EMOS library:
- !call jvod2uv( VO%c, DI%c, VO%T, U%c, V%c, VO%T )
- ! check arguments:
- !call Check( VO )
- !T = VO%T
- !call Check( DI, T )
- !call Check( U, T )
- !call Check( V, T )
- ! complex i for multiplication:
- z = (0.0, 1.0)
- ! index of coeff in array; initialise:
- kin = 0
- k = 0
-
- ! loop over all m:
- do m = 0, shi%T-1
- ! (m,m)
- n = m
- kin = kin + 1
- k = k + 1
- if ( m == 0 ) then
- U(k) = ( - delta(m,n+1)*VO(kin+1) )*R
- V(k) = ( delta(m,n+1)*DI(kin+1) )*R
- else
- U(k) = ( z*sigma(m,n)*DI(kin) - delta(m,n+1)*VO(kin+1) )*R
- V(k) = ( z*sigma(m,n)*VO(kin) + delta(m,n+1)*DI(kin+1) )*R
- end if
- ! internal area:
- if ( m <= shi%T-2 ) then
- do n = m+1, shi%T-2
- kin = kin + 1
- k = k + 1
- U(k) = ( delta(m,n)*VO(kin-1) + z*sigma(m,n)*DI(kin) - delta(m,n+1)*VO(kin+1) )*R
- V(k) = ( - delta(m,n)*DI(kin-1) + z*sigma(m,n)*VO(kin) + delta(m,n+1)*DI(kin+1) )*R
- end do
- ! (m,T-1)
- n = shi%T-1
- kin = kin + 1
- k = k + 1
- U(k) = ( delta(m,n)*VO(kin-1) + z*sigma(m,n)*DI(kin))*R
- V(k) = ( - delta(m,n)*DI(kin-1) + z*sigma(m,n)*VO(kin))*R
- end if
- ! (m,T)
- n = shi%T
- kin = kin + 1
- k = k + 1
- U(k) = ( delta(m,n)*VO(kin-1) )*R
- V(k) = ( - delta(m,n)*DI(kin-1) )*R
- kin = kin + shi_in%T - shi%T
-
- end do
- ! (T,T)
- k = k + 1
- U(k) = (0.0,0.0)
- V(k) = (0.0,0.0)
- contains
- pure real function delta( m, n )
- integer, intent(in) :: m, n
- delta = - sqrt( (n*n-m*m)*1.0/(4.0*n*n-1.0) ) / (n*1.0)
- end function delta
- pure real function sigma( m, n )
- integer, intent(in) :: m, n
- sigma = - m * 1.0/(n*(n+1.0))
- end function sigma
- end subroutine shi_vod2uv
-
-
-
-
- end module grid_type_sh
-
|