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

Diff of /trunk/filtrez/inifilr_hemisph.f

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

revision 167 by guez, Mon Aug 24 16:30:33 2015 UTC revision 169 by guez, Mon Sep 14 17:13:16 2015 UTC
# Line 4  module inifilr_hemisph_m Line 4  module inifilr_hemisph_m
4    
5  contains  contains
6    
7    subroutine inifilr_hemisph(rlat, colat0, rlamda, unit, eignfn, jfilt, &    subroutine inifilr_hemisph(rlat, rlamda, unit, eignfn, jfilt, matrice, &
8         matrice, matrinv)         matrinv)
9    
10      ! See notes.      ! See notes.
11    
12      USE dimens_m, ONLY : iim      USE dimens_m, ONLY : iim
13      use nr_util, only: pi, ifirstloc      use nr_util, only: pi, ifirstloc
14    
15      real, intent(in):: rlat(:) ! (n)      real, intent(in):: rlat(:) ! (n_lat)
16      ! latitudes, in rad, in [0, pi / 2[, in strictly ascending order      ! latitudes, in rad, in [0, pi / 2[, in strictly ascending order
17    
     REAL, intent(in):: colat0 ! > 0  
18      REAL, intent(in):: rlamda(2:) ! (2:iim) > 0, in descending order      REAL, intent(in):: rlamda(2:) ! (2:iim) > 0, in descending order
19      integer, intent(in):: unit      integer, intent(in):: unit
20    
# Line 23  contains Line 22  contains
22      ! eigenvectors of the discrete second derivative with respect to longitude      ! eigenvectors of the discrete second derivative with respect to longitude
23    
24      integer, intent(out):: jfilt      integer, intent(out):: jfilt
     real, pointer:: matrice(:, :, :) ! (iim, iim, n - jfilt + 1) matrice filtre  
25    
26      real, pointer, optional:: matrinv(:, :, :) ! (iim, iim, n - jfilt + 1)      real, pointer:: matrice(:, :, :) ! (iim, iim, n_lat - jfilt + 1)
27        ! matrice filtre
28    
29        real, pointer, optional:: matrinv(:, :, :) ! (iim, iim, n_lat - jfilt + 1)
30      ! matrice pour le filtre inverse      ! matrice pour le filtre inverse
31    
32      ! Local:      ! Local:
33    
34      integer n, i, j      integer n_lat, i, j
35      REAL eignft(iim, iim)      REAL eignft(iim, iim)
36    
37      ! Index of the mode from where modes are filtered:      ! Index of the mode from where modes are filtered:
38      integer, allocatable:: modfrst(:) ! (jfilt:n) in {2, ..., iim}      integer modfrst ! in {2, ..., iim}
39    
40      ! Filtering coefficients (lamda_max * cos(rlat) / lamda):      ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
41      real coefil(2:iim)      real coefil(2:iim)
42    
43      !-----------------------------------------------------------      !-----------------------------------------------------------
44    
45      n = size(rlat)      n_lat = size(rlat)
46      jfilt = ifirstloc(cos(rlat) < min(colat0, 1. / rlamda(iim)))      jfilt = ifirstloc(cos(rlat) < 1. / rlamda(iim))
47      allocate(modfrst(jfilt:n))      allocate(matrice(iim, iim, n_lat - jfilt + 1))
48        if (present(matrinv)) allocate(matrinv(iim, iim, n_lat - jfilt + 1))
49      DO j = jfilt, n  
50         modfrst(j) = ifirstloc(rlamda < 1. / cos(rlat(j)), my_lbound = 2)      DO j = jfilt, n_lat
51         write(unit, fmt = *) rlat(j) / pi * 180., modfrst(j)         modfrst = ifirstloc(rlamda < 1. / cos(rlat(j)), my_lbound = 2)
52      END DO         write(unit, fmt = *) rlat(j) / pi * 180., modfrst
53           coefil(modfrst:) = rlamda(modfrst:) * cos(rlat(j)) - 1.
54      allocate(matrice(iim, iim, n - jfilt + 1))         eignft(:modfrst - 1, :) = 0.
     if (present(matrinv)) allocate(matrinv(iim, iim, n - jfilt + 1))  
   
     DO j = jfilt, n  
        coefil(modfrst(j):) = rlamda(modfrst(j):) * cos(rlat(j)) - 1.  
        eignft(:modfrst(j) - 1, :) = 0.  
55    
56         forall (i = modfrst(j):iim) eignft(i, :) = eignfn(:, i) * coefil(i)         forall (i = modfrst:iim) eignft(i, :) = eignfn(:, i) * coefil(i)
57         matrice(:, :, j - jfilt + 1) = matmul(eignfn, eignft)         matrice(:, :, j - jfilt + 1) = matmul(eignfn, eignft)
58    
59         if (present(matrinv)) then         if (present(matrinv)) then
60            forall (i = modfrst(j):iim) eignft(i, :) = eignfn(:, i) * coefil(i) &            forall (i = modfrst:iim) eignft(i, :) = eignfn(:, i) * coefil(i) &
61                 / (1. + coefil(i))                 / (1. + coefil(i))
62            matrinv(:, :, j - jfilt + 1) = matmul(eignfn, eignft)            matrinv(:, :, j - jfilt + 1) = matmul(eignfn, eignft)
63         end if         end if
Legend:
Removed from v.167  
changed lines
  Added in v.169
  ViewVC Help
Powered by ViewVC 1.1.21