Changeset 1591


Ignore:
Timestamp:
Nov 25, 2011, 10:12:33 AM (13 years ago)
Author:
Ehouarn Millour
Message:

Minor corrections in filter: some arrays were oversized and the computation of the indexes from which the filter is applied could go wrong (in extreme cases, e.g. with very little grid points combined to a zoom).
EM

Location:
LMDZ5/trunk/libf/filtrez
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/filtrez/coefils.h

    r1146 r1591  
    11!
    2 ! $Header$
     2! $Id $
    33!
    44      COMMON/coefils/jfiltnu,jfiltsu,jfiltnv,jfiltsv,sddu(iim),sddv(iim)&
    … …  
    77     & ,coefilu2(iim,jjm),coefilv2(iim,jjm)
    88!c
    9       INTEGER jfiltnu,jfiltsu,jfiltnv,jfiltsv,modfrstu,modfrstv
     9      INTEGER jfiltnu ! index of the last lat line filtered in NH (U grid)
     10      INTEGER jfiltsu ! index of the first lat line filtered in SH (U grid)
     11      INTEGER jfiltnv ! index of the last lat line filtered in NH (V grid)
     12      INTEGER jfiltsv ! index of the first lat line filtered in SH (V grid)
     13      INTEGER modfrstu ! number of retained (ie: unfiltered) modes on U grid
     14      INTEGER modfrstv ! number of retained (ie: unfiltered) modes on V grid
    1015      REAL    sddu,sddv,unsddu,unsddv,coefilu,coefilv,eignfnu,eignfnv
    1116      REAL    coefilu2,coefilv2

  • LMDZ5/trunk/libf/filtrez/filtreg_mod.F90

    r1279 r1591  
     1!
     2! $Id $
     3!
    14MODULE filtreg_mod
    25
    … …  
    4245    INTEGER ixmineq
    4346#endif
    44     EXTERNAL  inifgn
    4547    !
    4648    ! ------------------------------------------------------------
    … …  
    7173    CALL inifgn(eignvl)
    7274    !
    73     PRINT *,' EIGNVL '
     75    PRINT *,'inifilr: EIGNVL '
    7476    PRINT 250,eignvl
    75 250 FORMAT( 1x,5e13.6)
     77250 FORMAT( 1x,5e14.6)
    7678    !
    7779    ! compute eigenvalues and eigenfunctions
    … …  
    113115#endif
    114116    !
     117    ! For a regular grid, we want the filter to start at latitudes
     118    ! corresponding to lengths dx of the same size as dy (in terms
     119    ! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
     120    !  <=> latitude=60 degrees).
     121    ! Same idea for the zoomed grid: start filtering polewards as soon
     122    ! as length dx becomes of the same size as dy
    115123    !
    116124    colat0  =  MIN( 0.5, dymin/dxmin )
    … …  
    158166    imx  = iim
    159167    !
    160     PRINT *,' TRUNCATION AT ',imx
    161     !
     168    PRINT *,'inifilr: TRUNCATION AT ',imx
     169    !
     170! Ehouarn: set up some defaults
     171    jfiltnu=2 ! avoid north pole
     172    jfiltsu=jjm ! avoid south pole (which is at jjm+1)
     173    jfiltnv=1 ! NB: no poles on the V grid
     174    jfiltsv=jjm
     175
    162176    DO j = 2, jjm/2+1
    163177       cof = COS( rlatu(j) )/ colat0
    164178       IF ( cof .LT. 1. ) THEN
    165           IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) jfiltnu= j
     179          IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) THEN
     180            jfiltnu= j
     181          ENDIF
    166182       ENDIF
    167183
    168184       cof = COS( rlatu(jjp1-j+1) )/ colat0
    169185       IF ( cof .LT. 1. ) THEN
    170           IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) &
     186          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) THEN
    171187               jfiltsu= jjp1-j+1
     188          ENDIF
    172189       ENDIF
    173190    ENDDO
    … …  
    176193       cof = COS( rlatv(j) )/ colat0
    177194       IF ( cof .LT. 1. ) THEN
    178           IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) jfiltnv= j
     195          IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) THEN
     196            jfiltnv= j
     197          ENDIF
    179198       ENDIF
    180199
    181200       cof = COS( rlatv(jjm-j+1) )/ colat0
    182201       IF ( cof .LT. 1. ) THEN
    183           IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) &
     202          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) THEN
    184203               jfiltsv= jjm-j+1
     204          ENDIF
    185205       ENDIF
    186206    ENDDO
    187207    !                                 
    188208
    189     IF ( jfiltnu.LE.0 ) jfiltnu=1
    190209    IF( jfiltnu.GT. jjm/2 +1 )  THEN
    191210       PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
    … …  
    193212    ENDIF
    194213
    195     IF( jfiltsu.LE.0) jfiltsu=1
    196214    IF( jfiltsu.GT.  jjm  +1 )  THEN
    197215       PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
    … …  
    199217    ENDIF
    200218
    201     IF( jfiltnv.LE.0) jfiltnv=1
    202219    IF( jfiltnv.GT. jjm/2    )  THEN
    203220       PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
    … …  
    205222    ENDIF
    206223
    207     IF( jfiltsv.LE.0) jfiltsv=1
    208224    IF( jfiltsv.GT.     jjm  )  THEN
    209225       PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
    … …  
    211227    ENDIF
    212228
    213     PRINT *,' jfiltnv jfiltsv jfiltnu jfiltsu ' , &
     229    PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
    214230         jfiltnv,jfiltsv,jfiltnu,jfiltsu
    215231
    216232    IF(first_call_inifilr) THEN
    217233       ALLOCATE(matriceun(iim,iim,jfiltnu))
    218        ALLOCATE(matriceus(iim,iim,jfiltsu))
     234       ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
    219235       ALLOCATE(matricevn(iim,iim,jfiltnv))
    220        ALLOCATE(matricevs(iim,iim,jfiltsv))
     236       ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
    221237       ALLOCATE( matrinvn(iim,iim,jfiltnu))
    222        ALLOCATE( matrinvs(iim,iim,jfiltsu))
     238       ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
    223239       first_call_inifilr = .FALSE.
    224240    ENDIF
    … …  
    230246    !
    231247    DO j = 1,jjm
     248    !default initialization: all modes are retained (i.e. no filtering)
    232249       modfrstu( j ) = iim
    233250       modfrstv( j ) = iim
    … …  
    306323
    307324    IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
    308 
     325! Ehouarn: and what are these for??? Trying to handle a limit case
     326!          where filters extend to and meet at the equator?
    309327       IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
    310328       IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
    … …  
    334352             eignft(i,k) = eignfnv(k,i) * coff
    335353          ENDDO
    336        ENDDO
     354       ENDDO ! of DO i=1,iim
    337355#ifdef CRAY
    338356       CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
    … …  
    350368             ENDDO
    351369          ENDDO
    352        ENDDO
    353 #endif
    354 #endif
    355 
    356     ENDDO
     370       ENDDO ! of DO k = 1, iim
     371#endif
     372#endif
     373
     374    ENDDO ! of DO j = 2, jfiltnu
    357375
    358376    DO j = jfiltsu, jjm
    … …  
    364382             eignft(i,k) = eignfnv(k,i) * coff
    365383          ENDDO
    366        ENDDO
     384       ENDDO ! of DO i=1,iim
    367385#ifdef CRAY
    368386       CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
    … …  
    381399             ENDDO
    382400          ENDDO
    383        ENDDO
    384 #endif
    385 #endif
    386 
    387     ENDDO
     401       ENDDO ! of DO k = 1, iim
     402#endif
     403#endif
     404
     405    ENDDO ! of DO j = jfiltsu, jjm
    388406
    389407    !   ...................................................................
    … …  
    421439#endif
    422440
    423     ENDDO
     441    ENDDO ! of DO j = 1, jfiltnv
    424442
    425443    DO j = jfiltsv, jjm
    … …  
    452470#endif
    453471
    454     ENDDO
     472    ENDDO ! of DO j = jfiltsv, jjm
    455473
    456474    !   ...................................................................
    … …  
    488506#endif
    489507
    490     ENDDO
     508    ENDDO ! of DO j = 2, jfiltnu
    491509
    492510    DO j = jfiltsu, jjm
    … …  
    518536#endif
    519537
    520     ENDDO
     538    ENDDO ! of DO j = jfiltsu, jjm
    521539
    522540    IF (use_filtre_fft) THEN
Note: See TracChangeset for help on using the changeset viewer.