/[lmdze]/trunk/filtrez/inifgn.f
ViewVC logotype

Diff of /trunk/filtrez/inifgn.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/Sources/filtrez/inifgn.f revision 151 by guez, Tue Jun 23 15:14:20 2015 UTC trunk/filtrez/inifgn.f revision 254 by guez, Mon Feb 5 10:39:38 2018 UTC
# Line 6  module inifgn_m Line 6  module inifgn_m
6    
7    private iim    private iim
8    
9    real sddu(iim), sddv(iim) ! SQRT(dx / di)    real sddu(iim), sddv(iim)
10    real unsddu(iim), unsddv(iim)    ! sdd[uv] = sqrt(2 pi / iim * (derivative of the longitudinal zoom
11      ! function)(rlon[uv]))
12    
13    real eignfnu(iim, iim), eignfnv(iim, iim)    real unsddu(iim), unsddv(iim)
   ! eigenfunctions of the discrete laplacian  
14    
15  contains  contains
16    
17    SUBROUTINE inifgn(dv)    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      ! From LMDZ4/libf/filtrez/inifgn.F, v 1.1.1.1 2004/05/19 12:53:09
20    
21      ! H. Upadyaya, O. Sharma      ! 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    
     use acc_m, only: acc  
26      USE dimens_m, ONLY: iim      USE dimens_m, ONLY: iim
27      USE dynetat0_m, ONLY: xprimu, xprimv      USE dynetat0_m, ONLY: xprimu, xprimv
28      use nr_util, only: pi      use jumble, only: new_unit
29      use numer_rec_95, only: jacobi, eigsrt      use numer_rec_95, only: jacobi, eigsrt
30    
31      real, intent(out):: dv(:) ! (iim) eigenvalues sorted in descending order      real, intent(out):: eignval_v(:) ! (iim)
32        ! eigenvalues sorted in descending order
33    
34        real, intent(out):: eignfnu(:, :), eignfnv(:, :) ! (iim, iim) eigenvectors
35    
36      ! Local:      ! Local:
37      REAL vec(iim, iim), vec1(iim, iim)  
38      REAL du(iim)      REAL delta(iim, iim) ! second derivative, symmetric, elements are angle^{-2}
39      INTEGER i, j, k, nrot  
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        INTEGER i, unit
45    
46      !----------------------------------------------------------------      !----------------------------------------------------------------
47    
# Line 42  contains Line 52  contains
52      unsddu = 1. / sddu      unsddu = 1. / sddu
53      unsddv = 1. / sddv      unsddv = 1. / sddv
54    
55      DO j = 1, iim      deriv_u = 0.
56         DO i = 1, iim      deriv_u(iim, 1) = unsddu(iim) * unsddv(1)
57            vec(i, j) = 0.      forall (i = 1:iim) deriv_u(i, i) = - unsddu(i) * unsddv(i)
58            vec1(i, j) = 0.      forall (i = 1:iim - 1) deriv_u(i, i + 1) = unsddu(i) * unsddv(i + 1)
59            eignfnv(i, j) = 0.  
60            eignfnu(i, j) = 0.      deriv_v = - transpose(deriv_u)
61         END DO  
62      END DO      delta = matmul(deriv_v, deriv_u) ! second derivative at v longitudes
63        CALL jacobi(delta, eignval_v, eignfnv)
64      eignfnv(1, 1) = - 1.      CALL eigsrt(eignval_v, eignfnv)
65      eignfnv(iim, 1) = 1.  
66      DO i = 1, iim - 1      delta = matmul(deriv_u, deriv_v) ! second derivative at u longitudes
67         eignfnv(i+1, i+1) = - 1.      CALL jacobi(delta, eignval_u, eignfnu)
68         eignfnv(i, i+1) = 1.      CALL eigsrt(eignval_u, eignfnu)
69      END DO  
70        call new_unit(unit)
71      DO j = 1, iim      open(unit, file = "inifgn_out.txt", status = "replace", action = "write")
72         DO i = 1, iim      write(unit, fmt = *) '"eignval_v"', eignval_v
73            eignfnv(i, j) = eignfnv(i, j) / (sddu(i) * sddv(j))      close(unit)
        END DO  
     END DO  
   
     DO j = 1, iim  
        DO i = 1, iim  
           eignfnu(i, j) = - eignfnv(j, i)  
        END DO  
     END DO  
   
     DO j = 1, iim  
        DO i = 1, iim  
           vec(i, j) = 0.0  
           vec1(i, j) = 0.0  
           DO k = 1, iim  
              vec(i, j) = vec(i, j) + eignfnu(i, k) * eignfnv(k, j)  
              vec1(i, j) = vec1(i, j) + eignfnv(i, k) * eignfnu(k, j)  
           END DO  
        END DO  
     END DO  
   
     CALL jacobi(vec, dv, eignfnv, nrot)  
     CALL acc(eignfnv)  
     CALL eigsrt(dv, eignfnv)  
   
     CALL jacobi(vec1, du, eignfnu, nrot)  
     CALL acc(eignfnu)  
     CALL eigsrt(du, eignfnu)  
74    
75    END SUBROUTINE inifgn    END SUBROUTINE inifgn
76    
Legend:
Removed from v.151  
changed lines
  Added in v.254
  ViewVC Help
Powered by ViewVC 1.1.21