dox/html/gaussmod_8f90_source.html

00001 ! =================
00002 ! SUBROUTINE INIGAU
00003 ! =================
00004
00005 subroutine inigau(klat,pz0,pzw)        ! pz0 & pzw are (kind=8) reals !!!
00006 implicit none
00007 integer                  :: klat       ! Number of Gaussian latitudes
00008 real (kind=8)            :: pz0(klat)  ! Gaussian abscissas
00009 real (kind=8)            :: pzw(klat)  ! Gaussian weights
00010 integer                  :: jlat       ! Latitudinal loop index
00011 integer                  :: jiter      ! Iteration loop index
00012 integer      , parameter :: NITER = 50 ! Maximum # of iterations
00013 real (kind=8), parameter :: PI    =  3.14159265358979_8
00014 real (kind=8), parameter :: ZEPS  =  1.0e-16 ! Convergence criterion
00015 real (kind=8) :: z0,z1,z2,z3,z4,z5
00016 real (kind=8) :: ql,qld
00017
00018 ! Compute Gaussian abscissas & weights
00019
00020 z0 = PI / (2*klat+1)
00021 z1 = 1.0_8 / (klat*klat*8)
00022 z4 = 2.0_8 / (klat*klat)
00023
00024 do jlat = 1 , klat/2
00025    z2 = z0 * (2*jlat - 0.5_8)
00026    z2 = cos(z2 + z1 / tan(z2))
00027    do jiter = 1 , NITER
00028       z3 = ql(klat,z2) * qld(klat,z2)
00029       z2 = z2 - z3
00030       if (abs(z3) < ZEPS) exit ! converged
00031    enddo ! jiter
00032    z5 = ql(klat-1,z2) / sqrt(klat - 0.5_8)
00033    pz0(jlat) = z2
00034    pzw(jlat) = z4 * (1.0_8 - z2 * z2) / (z5 * z5)
00035    pz0(klat-jlat+1) = -z2
00036    pzw(klat-jlat+1) = pzw(jlat)
00037 enddo ! jlat
00038
00039 return
00040 end subroutine inigau
00041
00042 ! ===========
00043 ! FUNCTION QL
00044 ! ===========
00045
00046 real (kind=8) function ql(k,p)
00047 implicit none
00048 integer      , intent(IN) :: k
00049 real (kind=8), intent(IN) :: p
00050 real (kind=8) :: z0,z1,z2,z3,z4
00051 integer :: j
00052 z0 = acos(p)
00053 z1 = 1.0
00054 z2 = 0.0
00055 do j = k , 0 , -2
00056    z3 = z1 * cos(z0 * j)
00057    z2 = z2 + z3
00058    z4 = (k-j+1) * (k+j) * 0.5_8
00059    z1 = z1 * z4 / (z4 + (j-1))
00060 enddo ! j
00061 if (mod(k,2) == 0) z2 = z2 - 0.5_8 * z3
00062
00063 z0 = sqrt(2.0_8)
00064 do j = 1 ,k
00065    z0 = z0 * sqrt(1.0_8 - 0.25_8 / (j*j))
00066 enddo ! j
00067 ql = z0 * z2
00068 return
00069 end function ql
00070
00071 ! ============
00072 ! FUNCTION QLD
00073 ! ============
00074
00075 real (kind=8) function qld(k,p)
00076 implicit none
00077 integer      , intent(IN) :: k
00078 real (kind=8), intent(IN) :: p
00079 real (kind=8) :: z
00080 real (kind=8) :: ql
00081
00082 z = p * ql(k,p) - sqrt((k + k + 1.0_8) / (k + k - 1.0_8)) * ql(k-1,p)
00083 qld = (p * p - 1.0_8) / (k * z)
00084
00085 return
00086 end function qld