trunk/filtrez/inifgn.f

module inifgn_m

  use dimens_m, only: iim

  IMPLICIT NONE

  private iim

  real sddu(iim), sddv(iim) 
  ! sdd[uv] = sqrt(2 pi / iim * (derivative of the longitudinal zoom
  ! function)(rlon[uv]))

  real unsddu(iim), unsddv(iim)

contains

  SUBROUTINE inifgn(eignval_v, eignfnu, eignfnv)

    ! From LMDZ4/libf/filtrez/inifgn.F, v 1.1.1.1 2004/05/19 12:53:09

    ! Authors: H. Upadyaya, O. Sharma

    ! Computes the eigenvalues and eigenvectors of the discrete analog
    ! of the second derivative with respect to longitude.

    use acc_m, only: acc 
    USE dimens_m, ONLY: iim
    USE dynetat0_m, ONLY: xprimu, xprimv
    use numer_rec_95, only: jacobi, eigsrt

    real, intent(out):: eignval_v(:) ! (iim)
    ! eigenvalues sorted in descending order

    real, intent(out):: eignfnu(:, :), eignfnv(:, :) ! (iim, iim) eigenvectors

    ! Local:

    REAL delta(iim, iim) ! second derivative, symmetric, elements are angle^{-2}

    REAL deriv_u(iim, iim), deriv_v(iim, iim) 
    ! first derivative at u and v longitudes, elements are angle^{-1}

    REAL eignval_u(iim)
    INTEGER i

    !----------------------------------------------------------------

    print *, "Call sequence information: inifgn"

    sddv = sqrt(xprimv(:iim))
    sddu = sqrt(xprimu(:iim))
    unsddu = 1. / sddu
    unsddv = 1. / sddv

    deriv_u = 0.
    deriv_u(iim, 1) = unsddu(iim) * unsddv(1)
    forall (i = 1:iim) deriv_u(i, i) = - unsddu(i) * unsddv(i)
    forall (i = 1:iim - 1) deriv_u(i, i + 1) = unsddu(i) * unsddv(i + 1)

    deriv_v = - transpose(deriv_u)

    delta = matmul(deriv_v, deriv_u) ! second derivative at v longitudes
    CALL jacobi(delta, eignval_v, eignfnv)
    CALL acc(eignfnv)
    CALL eigsrt(eignval_v, eignfnv)

    delta = matmul(deriv_u, deriv_v) ! second derivative at u longitudes
    CALL jacobi(delta, eignval_u, eignfnu)
    CALL acc(eignfnu)
    CALL eigsrt(eignval_u, eignfnu)

  END SUBROUTINE inifgn

end module inifgn_m
1	guez	113	module inifgn_m
2	guez	3
3	guez	140	use dimens_m, only: iim
4
5	guez	113	IMPLICIT NONE
6	guez	3
7	guez	140	private iim
8
9	guez	154	real sddu(iim), sddv(iim)
10			! sdd[uv] = sqrt(2 pi / iim * (derivative of the longitudinal zoom
11			! function)(rlon[uv]))
12
13	guez	140	real unsddu(iim), unsddv(iim)
14
15	guez	113	contains
16	guez	3
17	guez	154	SUBROUTINE inifgn(eignval_v, eignfnu, eignfnv)
18	guez	3
19	guez	113	! From LMDZ4/libf/filtrez/inifgn.F, v 1.1.1.1 2004/05/19 12:53:09
20	guez	3
21	guez	154	! Authors: H. Upadyaya, O. Sharma
22	guez	3
23	guez	154	! Computes the eigenvalues and eigenvectors of the discrete analog
24			! of the second derivative with respect to longitude.
25
26	guez	143	use acc_m, only: acc
27	guez	113	USE dimens_m, ONLY: iim
28	guez	139	USE dynetat0_m, ONLY: xprimu, xprimv
29	guez	143	use numer_rec_95, only: jacobi, eigsrt
30	guez	3
31	guez	154	real, intent(out):: eignval_v(:) ! (iim)
32			! eigenvalues sorted in descending order
33	guez	81
34	guez	154	real, intent(out):: eignfnu(:, :), eignfnv(:, :) ! (iim, iim) eigenvectors
35
36	guez	113	! Local:
37	guez	154
38	guez	161	REAL delta(iim, iim) ! second derivative, symmetric, elements are angle^{-2}
39	guez	154
40			REAL deriv_u(iim, iim), deriv_v(iim, iim)
41			! first derivative at u and v longitudes, elements are angle^{-1}
42
43			REAL eignval_u(iim)
44	guez	152	INTEGER i
45	guez	81
46	guez	113	!----------------------------------------------------------------
47	guez	81
48	guez	151	print *, "Call sequence information: inifgn"
49
50	guez	140	sddv = sqrt(xprimv(:iim))
51			sddu = sqrt(xprimu(:iim))
52			unsddu = 1. / sddu
53			unsddv = 1. / sddv
54	guez	81
55	guez	154	deriv_u = 0.
56			deriv_u(iim, 1) = unsddu(iim) * unsddv(1)
57			forall (i = 1:iim) deriv_u(i, i) = - unsddu(i) * unsddv(i)
58			forall (i = 1:iim - 1) deriv_u(i, i + 1) = unsddu(i) * unsddv(i + 1)
59	guez	81
60	guez	154	deriv_v = - transpose(deriv_u)
61	guez	140
62	guez	161	delta = matmul(deriv_v, deriv_u) ! second derivative at v longitudes
63			CALL jacobi(delta, eignval_v, eignfnv)
64	guez	143	CALL acc(eignfnv)
65	guez	154	CALL eigsrt(eignval_v, eignfnv)
66	guez	81
67	guez	161	delta = matmul(deriv_u, deriv_v) ! second derivative at u longitudes
68			CALL jacobi(delta, eignval_u, eignfnu)
69	guez	143	CALL acc(eignfnu)
70	guez	154	CALL eigsrt(eignval_u, eignfnu)
71	guez	81
72	guez	113	END SUBROUTINE inifgn
73	guez	81
74	guez	113	end module inifgn_m