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 jumble, only: new_unit
    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, unit

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

    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)

    call new_unit(unit)
    open(unit, file = "inifgn_out.txt", status = "replace", action = "write") 
    write(unit, fmt = *) '"eignval_v"', eignval_v
    close(unit)

  END SUBROUTINE inifgn

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