dox/html/legsym_8f90_source.html

00001 ! ************************************************************************
00002 ! * module legsym - direct and indirect Legendre transformation routines *
00003 ! * using symmetric and antisymmetric fourier coefficients               *
00004 ! * E. Kirk 08-Sep-2010 tested for T21 - T682 resolutions 32 & 64 bit    *
00005 ! ************************************************************************
00006
00007 module legsym
00008 ! ************************
00009 ! * Legendre Polynomials *
00010 ! ************************
00011
00012 integer :: ntru
00013 integer :: ntp1
00014 integer :: ncsp
00015 integer :: nesp
00016 integer :: nlon
00017 integer :: nlpp
00018 integer :: nhpp
00019 integer :: nlat
00020 integer :: nlev
00021 real    :: plavor
00022
00023 real   , allocatable :: qi(:,:) ! P(m,n) = Associated Legendre Polynomials
00024 real   , allocatable :: qj(:,:) ! Q(m,n) = Used for d/d(mu)
00025 real   , allocatable :: qc(:,:) ! P(m,n) * gwd              used in fc2sp
00026 real   , allocatable :: qe(:,:) ! Q(mn,) * gwd / cos2       used in mktend
00027 real   , allocatable :: qq(:,:) ! P(m,n) * gwd / cos2 * n * (n+1) / 2  "
00028 real   , allocatable :: qu(:,:) ! P(m,n) / (n*(n+1)) * m    used in dv2uv
00029 real   , allocatable :: qv(:,:) ! Q(m,n) / (n*(n+1))        used in dv2uv
00030 complex, allocatable :: qx(:,:) ! P(m,n) * gwd / cos2 * m   used in mktend
00031 end module legsym
00032
00033 ! =================
00034 ! SUBROUTINE LEGINI
00035 ! =================
00036
00037 subroutine legini(klat,klpp,kesp,klev,vorpla,sid,gwd)
00038 use legsym
00039 implicit none
00040
00041 integer :: klat
00042 integer :: klpp
00043 integer :: kesp
00044 integer :: klev
00045
00046 real          :: vorpla   ! planetary vorticity
00047 real (kind=8) :: sid(*)   ! sin(phi)
00048 real (kind=8) :: gwd(*)   ! Gaussian weight (phi)
00049
00050 integer :: jlat ! Latitude
00051 integer :: lm
00052 integer :: m
00053 integer :: n
00054
00055 real (kind=8) :: amsq
00056 real (kind=8) :: z1
00057 real (kind=8) :: z2
00058 real (kind=8) :: z3
00059 real (kind=8) :: f1m
00060 real (kind=8) :: f2m
00061 real (kind=8) :: znn1
00062 real (kind=8) :: zsin    ! sin
00063 real (kind=8) :: zcsq    ! cos2
00064 real (kind=8) :: zcos    ! cos
00065 real (kind=8) :: zgwd    ! gw
00066 real (kind=8) :: zgwdcsq ! gw / cos2
00067 real (kind=8) :: zpli(kesp)
00068 real (kind=8) :: zpld(kesp)
00069
00070 nlat = klat
00071 nlpp = klpp
00072 nhpp = klpp / 2
00073 nlon = klat + klat
00074 ntru = (nlon - 1) /3
00075 ntp1 = ntru + 1
00076 ncsp = ((ntru + 1) * (ntru + 2)) / 2
00077 nesp = kesp
00078 nlev = klev
00079
00080 plavor = vorpla
00081
00082 allocate(qi(ncsp,nhpp))
00083 allocate(qj(ncsp,nhpp))
00084 allocate(qc(ncsp,nhpp))
00085 allocate(qe(ncsp,nhpp))
00086 allocate(qx(ncsp,nhpp))
00087 allocate(qq(ncsp,nhpp))
00088 allocate(qu(ncsp,nhpp))
00089 allocate(qv(ncsp,nhpp))
00090
00091 do jlat = 1 , nhpp
00092
00093 ! set p(0,0) and p(0,1)
00094
00095    zgwd    = gwd(jlat)            ! gaussian weight - from inigau
00096    zsin    = sid(jlat)            ! sin(phi) - from inigau
00097    zcsq    = 1.0_8 - zsin * zsin  ! cos(phi) squared
00098    zgwdcsq = zgwd / zcsq          ! weight / cos squared
00099    zcos    = sqrt(zcsq)           ! cos(phi)
00100    f1m     = sqrt(1.5_8)
00101    zpli(1) = sqrt(0.5_8)
00102    zpli(2) = f1m * zsin
00103    zpld(1) = 0.0
00104    lm      = 2
00105
00106 ! loop over wavenumbers
00107
00108    do m = 0 , ntru
00109       if (m > 0) then
00110          lm  = lm + 1
00111          f2m = -f1m * sqrt(zcsq / (m+m))
00112          f1m =  f2m * sqrt(m+m + 3.0_8)
00113          zpli(lm) = f2m
00114          if (lm < ncsp) then
00115             lm = lm + 1
00116             zpli(lm  ) =       f1m * zsin
00117             zpld(lm-1) =  -m * f2m * zsin
00118          endif ! (lm < ncsp)
00119       endif ! (m > 0)
00120
00121       amsq = m * m
00122
00123       do n = m+2 , ntru
00124          lm = lm + 1
00125          z1 = sqrt(((n-1)*(n-1) - amsq) / (4*(n-1)*(n-1)-1))
00126          z2 = zsin * zpli(lm-1) - z1 * zpli(lm-2)
00127          zpli(lm  ) = z2 * sqrt((4*n*n-1) / (n*n-amsq))
00128          zpld(lm-1) = (1-n) * z2 + n * z1 * zpli(lm-2)
00129       enddo ! n
00130
00131       if (lm < ncsp) then ! mode (m,ntru)
00132          z3 = sqrt((ntru*ntru-amsq) / (4*ntru*ntru-1))
00133          zpld(lm)=-ntru*zsin*zpli(lm) + (ntru+ntru+1)*zpli(lm-1)*z3
00134       else                ! mode (ntru,ntru)
00135          zpld(lm)=-ntru*zsin*zpli(lm)
00136       endif
00137    enddo ! m
00138
00139    lm = 0
00140    do m = 0 , ntru
00141       do n = m , ntru
00142            lm = lm + 1
00143            znn1 = 0.0
00144            if (n > 0) znn1 = 1.0_8 / (n*(n+1))
00145            qi(lm,jlat) = zpli(lm)
00146            qj(lm,jlat) = zpld(lm)
00147            qc(lm,jlat) = zpli(lm) * zgwd
00148            qu(lm,jlat) = zpli(lm) * znn1 * m
00149            qv(lm,jlat) = zpld(lm) * znn1
00150            qe(lm,jlat) = zpld(lm) * zgwdcsq
00151            qq(lm,jlat) = zpli(lm) * zgwdcsq * n * (n+1) * 0.5_8
00152            qx(lm,jlat) = zpli(lm) * zgwdcsq * m * (0.0,1.0)
00153       enddo ! n
00154    enddo ! m
00155 enddo ! jlat
00156 return
00157 end
00158
00159
00160 ! ================
00161 ! SUBROUTINE FC2SP
00162 ! ================
00163
00164 subroutine fc2sp(fc,sp)
00165 use legsym
00166 implicit none
00167 complex, intent(in ) :: fc(nlon,nhpp)
00168 complex, intent(out) :: sp(nesp/2)
00169
00170 integer :: l ! Index for latitude
00171 integer :: m ! Index for zonal wavenumber
00172 integer :: w ! Index for spherical harmonic
00173 integer :: e ! Index for last wavenumber
00174
00175 sp(:) = (0.0,0.0)
00176
00177 do l = 1 , nhpp
00178    w = 1
00179    do m = 1 , ntp1
00180       e = w + ntp1 - m
00181       sp(w  :e:2) = sp(w  :e:2) + qc(w  :e:2,l) * (fc(m,l) + fc(m+nlat,l))
00182       sp(w+1:e:2) = sp(w+1:e:2) + qc(w+1:e:2,l) * (fc(m,l) - fc(m+nlat,l))
00183       w = e + 1
00184    enddo ! m
00185 enddo ! l
00186 return
00187 end
00188
00189
00190 ! ================
00191 ! SUBROUTINE SP2FC
00192 ! ================
00193
00194 subroutine sp2fc(sp,fc) ! Spectral to Fourier
00195 use legsym
00196 implicit none
00197
00198 complex :: sp(ncsp)      ! Coefficients of spherical harmonics
00199 complex :: fc(nlon,nhpp) ! Fourier coefficients
00200
00201 integer :: l ! Loop index for latitude
00202 integer :: m ! Loop index for zonal wavenumber m
00203 integer :: w ! Index for spectral mode
00204 integer :: e ! Index for last wavenumber
00205
00206 complex :: fs,fa
00207
00208 fc(:,:) = (0.0,0.0)
00209
00210 do l = 1 , nhpp
00211   w = 1
00212   do m = 1 ,ntp1
00213     e = w + ntp1 - m
00214     fs = dot_product(qi(w  :e:2,l),sp(w  :e:2))
00215     fa = dot_product(qi(w+1:e:2,l),sp(w+1:e:2))
00216     fc(m     ,l) = fs + fa
00217     fc(m+nlat,l) = fs - fa
00218     w = e + 1
00219   enddo ! m
00220 enddo ! l
00221 return
00222 end
00223
00224
00225 ! ===================
00226 ! SUBROUTINE SP2FCDMU
00227 ! ===================
00228
00229 subroutine sp2fcdmu(sp,fc) ! Spectral to Fourier d/dmu
00230 use legsym
00231 implicit none
00232
00233 complex :: sp(ncsp)      ! Coefficients of spherical harmonics
00234 complex :: fc(nlon,nhpp) ! Fourier coefficients
00235
00236 integer :: l ! Loop index for latitude
00237 integer :: m ! Loop index for zonal wavenumber m
00238 integer :: w ! Index for spectral mode
00239 integer :: e ! Index for last wavenumber
00240
00241 complex :: fs,fa
00242
00243 fc(:,:) = (0.0,0.0)
00244
00245 do l = 1 , nhpp
00246   w = 1
00247   do m = 1 , ntp1
00248     e = w + ntp1 - m
00249     fs = dot_product(qj(w  :e:2,l),sp(w  :e:2))
00250     fa = dot_product(qj(w+1:e:2,l),sp(w+1:e:2))
00251     fc(m     ,l) = fa + fs
00252     fc(m+nlat,l) = fa - fs
00253     w = e + 1
00254   enddo ! m
00255 enddo ! l
00256 return
00257 end
00258
00259
00260 ! ================
00261 ! SUBROUTINE DV2UV
00262 ! ================
00263
00264 ! This is an alternative subroutine for computing U and V from Div and Vor.
00265 ! It looks much prettier than the regular one, but unfortunately it is slower
00266 ! if compiling with "gfortran".
00267 ! I leave it (unused) in this module for educational purposes.
00268
00269 subroutine dv2uv_alt(pd,pz,pu,pv)
00270 use legsym
00271 implicit none
00272 complex, parameter :: i = (0.0,1.0)
00273 complex :: pd(nesp/2)    ! Spherical harmonics  of divergence
00274 complex :: pz(nesp/2)    ! Spherical harmonics  of vorticity
00275 complex :: pu(nlon,nhpp) ! Fourier coefficients of u
00276 complex :: pv(nlon,nhpp) ! Fourier coefficients of v
00277 complex :: zsave,uds,vds,uzs,vzs,uda,vda,uza,vza
00278
00279 integer :: l ! Loop index for latitude
00280 integer :: m ! Loop index for zonal wavenumber m
00281 integer :: w ! Loop index for spectral mode
00282 integer :: e ! End index
00283
00284 pu(:,:) = (0.0,0.0)
00285 pv(:,:) = (0.0,0.0)
00286
00287 zsave = pz(2)                      ! Save mode(0,1) of vorticity
00288 pz(2) = zsave - cmplx(plavor,0.0)  ! Convert pz from absolute to relative vorticity
00289
00290 do l = 1 , nhpp
00291   w = 1
00292   do m = 1 ,ntp1
00293     e = w + ntp1 - m
00294     uds = i * dot_product(qu(w  :e:2,l),pd(w:  e:2))
00295     vds =     dot_product(qv(w  :e:2,l),pd(w:  e:2))
00296     uzs = i * dot_product(qu(w  :e:2,l),pz(w:  e:2))
00297     vzs =     dot_product(qv(w  :e:2,l),pz(w:  e:2))
00298     uda = i * dot_product(qu(w+1:e:2,l),pd(w+1:e:2))
00299     vda =     dot_product(qv(w+1:e:2,l),pd(w+1:e:2))
00300     uza = i * dot_product(qu(w+1:e:2,l),pz(w+1:e:2))
00301     vza =     dot_product(qv(w+1:e:2,l),pz(w+1:e:2))
00302     pu(m     ,l) =  vzs - uds + vza - uda
00303     pu(m+nlat,l) = -vzs - uds + vza + uda
00304     pv(m     ,l) = -vds - uzs - vda - uza
00305     pv(m+nlat,l) =  vds - uzs - vda + uza
00306     w = e + 1
00307   enddo ! m
00308 enddo ! l
00309 pz(2) = zsave
00310 return
00311 end
00312
00313
00314 ! ================
00315 ! SUBROUTINE DV2UV
00316 ! ================
00317
00318 subroutine dv2uv(pd,pz,pu,pv)
00319 use legsym
00320 implicit none
00321
00322 real :: pd(2,nesp/2)    ! Spherical harmonics  of divergence
00323 real :: pz(2,nesp/2)    ! Spherical harmonics  of vorticity
00324 real :: pu(2,nlon,nhpp) ! Fourier coefficients of u
00325 real :: pv(2,nlon,nhpp) ! Fourier coefficients of v
00326 real :: zsave
00327 real :: unr,uni,usr,usi,vnr,vni,vsr,vsi
00328 real :: zdr,zdi,zzr,zzi
00329
00330 integer :: l ! Loop index for latitude
00331 integer :: m ! Loop index for zonal wavenumber m
00332 integer :: n ! Loop index for total wavenumber n
00333 integer :: w ! Loop index for spectral mode
00334
00335 pu(:,:,:) = 0.0
00336 pv(:,:,:) = 0.0
00337
00338 zsave = pz(1,2)           ! Save mode(0,1) of vorticity
00339 pz(1,2) = zsave - plavor  ! Convert pz from absolute to relative vorticity
00340
00341 do l = 1 , nhpp
00342    w = 1
00343    do m = 1 , ntp1
00344       unr = 0.0 ! u - north - real
00345       uni = 0.0 ! u - north - imag
00346       usr = 0.0 ! u - south - real
00347       usi = 0.0 ! u - south - imag
00348       vnr = 0.0 ! v - north - real
00349       vni = 0.0 ! v - north - imag
00350       vsr = 0.0 ! v - south - real
00351       vsi = 0.0 ! v - south - imag
00352
00353 !     process two modes per iteration, one symmetric (m+n = even) and one anti
00354 !     we start the loop with (n=m), so the starting mode is always symmetric
00355
00356       do n = m , ntru , 2
00357          zdr = qu(w,l) * pd(1,w)  ! symmetric mode
00358          zdi = qu(w,l) * pd(2,w)
00359          zzr = qv(w,l) * pz(1,w)
00360          zzi = qv(w,l) * pz(2,w)
00361          unr = unr + zzr + zdi
00362          uni = uni + zzi - zdr
00363          usr = usr - zzr + zdi
00364          usi = usi - zzi - zdr
00365          zzr = qu(w,l) * pz(1,w)
00366          zzi = qu(w,l) * pz(2,w)
00367          zdr = qv(w,l) * pd(1,w)
00368          zdi = qv(w,l) * pd(2,w)
00369          vnr = vnr + zzi - zdr
00370          vni = vni - zzr - zdi
00371          vsr = vsr + zzi + zdr
00372          vsi = vsi - zzr + zdi
00373          w = w + 1
00374          zdr = qu(w,l) * pd(1,w)  ! antisymmetric mode
00375          zdi = qu(w,l) * pd(2,w)
00376          zzr = qv(w,l) * pz(1,w)
00377          zzi = qv(w,l) * pz(2,w)
00378          unr = unr + zzr + zdi
00379          uni = uni + zzi - zdr
00380          usr = usr + zzr - zdi
00381          usi = usi + zzi + zdr
00382          zzr = qu(w,l) * pz(1,w)
00383          zzi = qu(w,l) * pz(2,w)
00384          zdr = qv(w,l) * pd(1,w)
00385          zdi = qv(w,l) * pd(2,w)
00386          vnr = vnr + zzi - zdr
00387          vni = vni - zzr - zdi
00388          vsr = vsr - zzi - zdr
00389          vsi = vsi + zzr - zdi
00390          w = w + 1
00391       enddo
00392       if (n == ntp1) then         ! additional symmetric mode
00393          zdr = qu(w,l) * pd(1,w)  ! if (ntp1-m) is even
00394          zdi = qu(w,l) * pd(2,w)
00395          zzr = qv(w,l) * pz(1,w)
00396          zzi = qv(w,l) * pz(2,w)
00397          unr = unr + zzr + zdi
00398          uni = uni + zzi - zdr
00399          usr = usr - zzr + zdi
00400          usi = usi - zzi - zdr
00401          zzr = qu(w,l) * pz(1,w)
00402          zzi = qu(w,l) * pz(2,w)
00403          zdr = qv(w,l) * pd(1,w)
00404          zdi = qv(w,l) * pd(2,w)
00405          vnr = vnr + zzi - zdr
00406          vni = vni - zzr - zdi
00407          vsr = vsr + zzi + zdr
00408          vsi = vsi - zzr + zdi
00409          w = w + 1
00410       endif
00411       pu(1,m     ,l) = unr
00412       pu(2,m     ,l) = uni
00413       pu(1,m+nlat,l) = usr
00414       pu(2,m+nlat,l) = usi
00415       pv(1,m     ,l) = vnr
00416       pv(2,m     ,l) = vni
00417       pv(1,m+nlat,l) = vsr
00418       pv(2,m+nlat,l) = vsi
00419    enddo ! m
00420 enddo ! l
00421 pz(1,2) = zsave ! Restore pz to absolute vorticity
00422 return
00423 end
00424
00425
00426 ! =================
00427 ! SUBROUTINE MKTEND
00428 ! =================
00429
00430 subroutine mktend(d,t,z,tn,fu,fv,ke,ut,vt)
00431 use legsym
00432 implicit none
00433
00434 complex, intent(in) :: tn(nlon,nhpp)
00435 complex, intent(in) :: fu(nlon,nhpp)
00436 complex, intent(in) :: fv(nlon,nhpp)
00437 complex, intent(in) :: ke(nlon,nhpp)
00438 complex, intent(in) :: ut(nlon,nhpp)
00439 complex, intent(in) :: vt(nlon,nhpp)
00440
00441 complex, intent(out) :: d(nesp/2)
00442 complex, intent(out) :: t(nesp/2)
00443 complex, intent(out) :: z(nesp/2)
00444
00445 integer :: l ! Loop index for latitude
00446 integer :: m ! Loop index for zonal wavenumber m
00447 integer :: w ! Loop index for spectral mode
00448 integer :: e ! End index for w
00449
00450 complex :: fus,fua,fvs,fva,kes,kea,tns,tna,uts,uta,vts,vta
00451
00452 d(:) = (0.0,0.0) ! divergence
00453 t(:) = (0.0,0.0) ! temperature
00454 z(:) = (0.0,0.0) ! vorticity
00455
00456 do l = 1 , nhpp  ! process pairs of Nort-South latitudes
00457    w = 1
00458    do m = 1 , ntp1
00459       kes = ke(m,l) + ke(m+nlat,l) ; kea = ke(m,l) - ke(m+nlat,l)
00460       fvs = fv(m,l) + fv(m+nlat,l) ; fva = fv(m,l) - fv(m+nlat,l)
00461       fus = fu(m,l) + fu(m+nlat,l) ; fua = fu(m,l) - fu(m+nlat,l)
00462       uts = ut(m,l) + ut(m+nlat,l) ; uta = ut(m,l) - ut(m+nlat,l)
00463       vts = vt(m,l) + vt(m+nlat,l) ; vta = vt(m,l) - vt(m+nlat,l)
00464       tns = tn(m,l) + tn(m+nlat,l) ; tna = tn(m,l) - tn(m+nlat,l)
00465       e = w + ntp1 - m    ! vector of symmetric modes
00466       d(w:e:2) = d(w:e:2) + qq(w:e:2,l) * kes - qe(w:e:2,l) * fva + qx(w:e:2,l) * fus
00467       t(w:e:2) = t(w:e:2) + qe(w:e:2,l) * vta + qc(w:e:2,l) * tns - qx(w:e:2,l) * uts
00468       z(w:e:2) = z(w:e:2) + qe(w:e:2,l) * fua + qx(w:e:2,l) * fvs
00469       w = w + 1           ! vector of antisymmetric modes
00470       d(w:e:2) = d(w:e:2) + qq(w:e:2,l) * kea - qe(w:e:2,l) * fvs + qx(w:e:2,l) * fua
00471       t(w:e:2) = t(w:e:2) + qe(w:e:2,l) * vts + qc(w:e:2,l) * tna - qx(w:e:2,l) * uta
00472       z(w:e:2) = z(w:e:2) + qe(w:e:2,l) * fus + qx(w:e:2,l) * fva
00473       w = e + 1
00474    enddo ! m
00475 enddo ! l
00476 return
00477 end
00478
00479
00480 ! ==================
00481 ! SUBROUTINE REG2ALT
00482 ! ==================
00483
00484 subroutine reg2alt(pr,klev)
00485 use legsym
00486 implicit none
00487
00488 real :: pr(nlon,nlat,klev)
00489 real :: pa(nlon,nlat,klev)
00490
00491 integer :: jlat
00492 integer :: klev
00493
00494 do jlat = 1 , nlat / 2
00495   pa(:,2*jlat-1,:) = pr(:,jlat       ,:)
00496   pa(:,2*jlat  ,:) = pr(:,nlat-jlat+1,:)
00497 enddo
00498
00499 pr = pa
00500
00501 return
00502 end
00503
00504
00505 ! ==================
00506 ! SUBROUTINE ALT2REG
00507 ! ==================
00508
00509 subroutine alt2reg(pa,klev)
00510 use legsym
00511 implicit none
00512
00513 real :: pa(nlon,nlat,klev)
00514 real :: pr(nlon,nlat,klev)
00515
00516 integer :: jlat
00517 integer :: klev
00518
00519 do jlat = 1 , nlat / 2
00520   pr(:,jlat       ,:) = pa(:,2*jlat-1,:)
00521   pr(:,nlat-jlat+1,:) = pa(:,2*jlat  ,:)
00522 enddo
00523
00524 pa = pr
00525
00526 return
00527 end
00528
00529
00530 ! ================
00531 ! SUBROUTINE ALTCS
00532 ! ================
00533
00534 subroutine altcs(pcs)
00535 use legsym
00536 implicit none
00537 real :: pcs(nlat,nlev)
00538 real :: pal(nlat,nlev)
00539 integer :: jlat
00540
00541 do jlat = 1 , nlat / 2
00542   pal(jlat       ,:) = pcs(2*jlat-1,:)
00543   pal(nlat-jlat+1,:) = pcs(2*jlat  ,:)
00544 enddo
00545
00546 pcs = pal
00547
00548 return
00549 end
00550
00551
00552 ! =================
00553 ! SUBROUTINE ALTLAT
00554 ! =================
00555
00556 subroutine altlat(pr,klat)
00557 implicit none
00558 integer :: jlat
00559 integer :: klat
00560 real    :: pr(klat)  ! regular     grid
00561 real    :: pa(klat)  ! alternating grid
00562
00563 do jlat = 1 , klat / 2
00564   pa(2*jlat-1) = pr(jlat       )
00565   pa(2*jlat  ) = pr(klat-jlat+1)
00566 enddo
00567
00568 pr(:) = pa(:)
00569
00570 return
00571 end