/[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 143 by guez, Tue Jun 9 14:32:46 2015 UTC trunk/filtrez/inifgn.f revision 265 by guez, Tue Mar 20 09:35:59 2018 UTC
# Line 1  Line 1 
1  module inifgn_m  module inifgn_m
2    
3    use dimens_m, only: iim    use dimensions, only: iim
4    
5    IMPLICIT NONE    IMPLICIT NONE
6    
7    private iim    private iim
8    
9    real sddu(iim), sddv(iim) ! SQRT(dx)    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      use acc_m, only: acc      ! Computes the eigenvalues and eigenvectors of the discrete analog
24      USE dimens_m, ONLY: iim      ! of the second derivative with respect to longitude.
25    
26        USE dimensions, 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)      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    
48        print *, "Call sequence information: inifgn"
49    
50      sddv = sqrt(xprimv(:iim))      sddv = sqrt(xprimv(:iim))
51      sddu = sqrt(xprimu(:iim))      sddu = sqrt(xprimu(:iim))
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.143  
changed lines
  Added in v.265
  ViewVC Help
Powered by ViewVC 1.1.21