Changeset 213 in lmdz_wrf for trunk/tools


Ignore:
Timestamp:
Jan 9, 2015, 3:21:22 PM (10 years ago)
Author:
lfita
Message:

Version:

  • without filter functions
  • without python typos
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/iniaqua.py

    r212 r213  
    2121## g.e. # iniaqua.py -d 32,32,39 -p hybdrid -o WRF -t 19791201000000 -z param
    2222
    23 # from dyn3d/fxy_new.h
    24 #c
    25 #c....stretching in x...
    26 #c
    27 def fripx(ri,dx):
    28     fname = 'fripx'
    29 
    30     ripx = (ri-1.0) *2.*np.pi/np.float(dx)
    31 
    32     return ripx
    33 
    34 def ffx(ri):
    35     fname = 'ffx'
    36 
    37     fx = fripx(ri) + transx  + alphax * np.sin(fripx(ri)+transx-pxo ) - np.pi
    38 
    39     return fx
    40 
    41 def ffxprim(ri,dx):
    42     fname = 'ffxprim'
    43 
    44     fxprim = 2.*np.pi/REAL(dx)*( 1.+ alphax*np.cos( fripx(ri)+transx-pxo ) )
    45 
    46     return fxprim
    47 
    48 #c....stretching in y...
    49 #c
    50 def fbigy(rj,jjp1,jjm):
    51     fname = 'fbigy'
    52 
    53     bigy = 2.*(np.float(jjp1)-rj)*np.pi/jjm
    54 
    55     return bigy
    56 
    57 def fy(rj):
    58     fname = 'fy'
    59 
    60     fy = ( fbigy(rj) + transy + alphay * np.sin(fbigy(rj)+transy-pyo) )/2. - np.pi/2.
    61 
    62     return fy
    63 
    64 def ffyprim(rj,jjm):
    65     fname = 'ffyprim'
    66 
    67     fyprim = ( np.pi/jjm ) * ( 1. + alphay * np.cos( fbigy(rj)+transy-pyo ) )
    68 
    69     return fyprim
    70 
    71 def fxy(dx, dy):
    72 
    73     """!
    74        ! $Id: fxy.F 1403 2010-07-01 09:02:53Z fairhead $
    75        !
    76 
    77        c     Auteur  :  P. Le Van
    78        c
    79        c     Calcul  des longitudes et des latitudes  pour une fonction f(x,y)
    80        c           a tangente sinusoidale et eventuellement avec zoom  .
    81        c
    82        c
    83     """
    84 
    85     fname ='fxy'
    86 
    87 #c    ......  calcul  des  latitudes  et de y'   .....
    88 #c
    89     vrlatu = np.zeros((dy+1), dtype=np.float)
    90     vyprimu = np.zeros((dy+1), dtype=np.float)
    91 
    92     vrlatu = ffy(np.arange(dy+1)*1. + 1.)
    93     vyprimu = ffy(np.arange(dy+1)*1. + 1.)
    94 
    95     vrlatv = np.zeros((dy), dtype=np.float)
    96     vrlatu1 = np.zeros((dy), dtype=np.float)
    97     vrlatu2 = np.zeros((dy), dtype=np.float)
    98     yprimv = np.zeros((dy), dtype=np.float)
    99     vyprimu1 = np.zeros((dy), dtype=np.float)
    100     vyprimu2 = np.zeros((dy), dtype=np.float)
    101 
    102     vrlatv = ffy(np.arange(dy)*1. + 1. + 0.5)
    103     vrlatu1 = ffy(np.arange(dy)*1. + 1. + 0.25)
    104     vrlatu2 = ffy(np.arange(dy)*1. + 1. + 0.75)
    105 
    106     vyprimv = fyprim(np.arange(dy)*1. + 1. + 0.5)
    107     vyprimu1 = fyprim(np.arange(dy)*1. + 1. + 0.25)
    108     vyprimu2 = fyprim(np.arange(dy)*1. + 1. + 0.75)
    109 
    110 #c
    111 #c     .....  calcul   des  longitudes et de  x'   .....
    112 #c
    113     vrlonv = np.zeros((dx+1), dtype=np.float)
    114     vrlonu = np.zeros((dx+1), dtype=np.float)
    115     vrlonm025 = np.zeros((dx+1), dtype=np.float)
    116     vrlonp025 = np.zeros((dx+1), dtype=np.float)
    117     vxprimv = np.zeros((dx+1), dtype=np.float)
    118     vxprimu = np.zeros((dx+1), dtype=np.float)
    119     vxprimm025 = np.zeros((dx+1), dtype=np.float)
    120     vxprimo025 = np.zeros((dx+1), dtype=np.float)
    121 
    122     vrlonv = ffy(np.arange(dx+1)*1. + 1.)
    123     vrlonu = ffy(np.arange(dx+1)*1. + 1. + 0.5)
    124     vrlonm025 = ffy(np.arange(dx+1)*1. + 1. - 0.25)
    125     vrlonp025 = ffy(np.arange(dx+1)*1. + 1. + 0.25)
    126 
    127     vxprimv = fxprim(np.arange(dx+1)*1. + 1.)
    128     vxprimu = fxprim(np.arange(dx+1)*1. + 1. + 0.5)
    129     vxprimm025 = fxprim(np.arange(dx+1)*1. + 1. - 0.25)
    130     vxprimp025 = fxprim(np.arange(dx+1)*1. + 1. + 0.25)
    131 
    132     return vrlatu, vyprimu, vrlatv, vyprimv, vrlatu1, vyprimu1, vrlatu2, vyprimu2,   \
    133       vrlonu, vxprimu, vrlonv, vxprimv, vrlonm025, vxprimm025, vrlonp025, vxprimp025
    134 
    135 def jacobi(A,N,NP):
    136     """!
    137        ! $Id: jacobi.F90 1289 2009-12-18 14:51:22Z emillour $
    138        !
    139     """
    140     fname = 'jacobi'
    141 
    142     D = np.zeros((NP), dtype=np.float)
    143     V = np.zeros((NP,NP), dtype=np.float)
    144     B = np.zeros((NP), dtype=np.float)
    145     Z = np.zeros((NP), dtype=np.float)
    146 
    147     for IP in range(N):
    148         V[IP,IP]=1.
    149         B[IP] = A[IP,IP]
    150 
    151     D = B
    152     NROT = 0
    153     for I in range(50):
    154 #      DO I=1,50 ! 50? I suspect this should be NP
    155 #                !     but convergence is fast enough anyway
    156         SM = 0.
    157         for IP in range(N-1):
    158             for IQ in range(IP+1,N):
    159                 SM = SM + np.abs(A[IP,IQ])
    160 
    161         if SM  == 0.: return A, D, V, NROT
    162 
    163         if I < 4:
    164           TRESH = 0.2*SM/N**2
    165         else:
    166           TRESH = 0.
    167 
    168         for IP in range(N-1):
    169             for IQ in range(IP+1,N):
    170                 G = 100.*np.abs(A[IP,IQ])
    171                 if I > 4 and np.abs(D[I]) + G == np.abs(D[IP]) and                   \
    172                   np.abs(D[IQ])+G == np.abs(D[IQ]):
    173                     A[IP,IQ] = 0.
    174                 elif np.abs(A[IP,IQ]) > TRESH:
    175                     H = D[IQ]-D[IP]
    176                     if np.abs(H)+G == np.abs(H):
    177                         T = A[IP,IQ]/H
    178                     else:
    179                         THETA = 0.5*H/A[IP,IQ]
    180                         T = 1./(np.abs(THETA)+np.sqrt(1.+THETA**2))
    181                         if THETA < 0.: T = -T
    182 
    183                     C=1./SQRT(1+T**2)
    184                     S=T*C
    185                     TAU=S/(1.+C)
    186                     H=T*A[IP,IQ]
    187                     Z[IP]=Z[IP]-H
    188                     Z[IQ]=Z[IQ]+H
    189                     D[IP]=D[IP]-H
    190                     D[IQ]=D[IQ]+H
    191                     A[IP,IQ]=0.
    192                     for J in range(IP-1):
    193                         G=A[J,IP]
    194                         H=A[J,IQ]
    195                         A[J,IP]=G-S*(H+G*TAU)
    196                         A[J,IQ]=H+S*(G-H*TAU)
    197 
    198                     for J in range(IP+1,IQ-1):
    199                         G = A[IP,J]
    200                         H = A[J,IQ]
    201                         A[IP,J] = G-S*(H+G*TAU)
    202                         A[J,IQ] = H+S*(G-H*TAU)
    203 
    204                     for J in range(IQ+1,N):
    205                         G = A[IP,J]
    206                         H = A[IQ,J]
    207                         A[IP,J] = G-S*(H+G*TAU)
    208                         A[IQ,J] = H+S*(G-H*TAU)
    209 
    210                     for J in range(N):
    211                         G = V[J,IP]
    212                         H = V[J,IQ]
    213                         V[J,IP] = G-S*(H+G*TAU)
    214                         V[J,IQ] = H+S*(G-H*TAU)
    215 
    216                     NROT=NROT+1
    217 
    218           B = B + Z
    219           D = B
    220           Z = 0.
    221 
    222     print errormsg
    223     print '  ' + fname + ': 50 iterations should never happen !!'
    224     quit(-1)
    225 
    226     return
    227 
    228 def acc(vec,im):
    229     """!
    230        ! $Header$
    231        !
    232     """
    233     fname ='acc'
    234 
    235     d = np.zeros((im), dtype=np.float)
    236 
    237     for j in range(im):
    238         for i in range(im):
    239             d[i] = vec[i,j]*vec[i,j]
    240 
    241         vsum = np.sqrt(np.sum(d))
    242         for i in range(im):
    243             vec[i,j] = vec[i,j] / vsum
    244 
    245     return vec, d
    246 
    247 def eigen_sort(d,n,np):
    248     """!
    249        ! $Header$
    250        !
    251     """
    252     fname = 'eigen_sort'
    253 
    254     v = np.zeros((np,np), dtype=np.float)
    255 
    256     for i in range(n-1):     
    257         k=i
    258         p=d[i]
    259         for j in range(i+1,n):
    260             if d[j] >= p:
    261                 k=j
    262                 p=d[j]
    263          
    264         if k != i:
    265             d[k]=d[i]
    266             d[i]=p
    267             for j in range(n):
    268                 p=v[j,i]
    269                 v[j,i]=v[j,k]
    270                 v[j,k]=p
    271 
    272     return d, v
    273 
    274 def inifgn(dx, dy, dv):
    275     """! $Header: /home/cvsroot/LMDZ4/libf/filtrez/inifgn.F,v 1.1.1.1 2004-05-19 12:53:09 lmdzadmin Exp $
    276        !
    277        c 
    278        c    ...  H.Upadyaya , O.Sharma  ...
    279        c
    280     """
    281     fname = 'inifgn'
    282 
    283     imm1  = dx - 1
    284 
    285     rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu, xprimu,   \
    286       rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025 = fxy(dx,dy)
    287 
    288     dlonu = xprimu
    289     dlonv = xprimv
    290 
    291     sddv = np.sqrt(dlonv)
    292     sddu = np.sqrt(dlonu)
    293     unsddu = 1./sddu
    294     unsddv = 1./sddv
    295 
    296     vec = np.zeros((dx,dy), dtype=np.float)
    297     vec1 = np.zeros((dx,dy), dtype=np.float)
    298     eigenfnv = np.zeros((dx,dy), dtype=np.float)
    299     eigenfnu = np.zeros((dx,dy), dtype=np.float)
    300 c
    301 c
    302     eignfnv[0,0] = -1.
    303     eignfnv[dx-1,0] = 1.
    304 
    305     for i in range(imm1):
    306         eignfnv[i+1,i+1]= -1.
    307         eignfnv[i,i+1] = 1.
    308 
    309     for j in range(dx):
    310         for i in range(dx):
    311             eignfnv[i,j] = eignfnv[i,j]/(sddu[i]*sddv[j])
    312 
    313     for j in range(dx):
    314         for i in range(dx):
    315             eignfnu[i,j] = -eignfnv[j,i]
    316 
    317     for j in range(dx):
    318         for i in range(dx):
    319             vec[i,j] = 0.0
    320             vec1[i,j] = 0.0
    321             for k in range(dx):
    322                 vec[i,j] = vec[i,j]  + eignfnu[i,k]*eignfnv[k,j]
    323                 vec1[i,j] = vec1[i,j] + eignfnv[i,k]*eignfnu[k,j]
    324 
    325     vec, dv, eignfnv, nrot = jacobi(vec,iim,iim)
    326     eignfnv, d = acc(eignfnv,iim)
    327     dv, eignfnv = eigen_sort(dv,iim,iim)
    328 c
    329     vec1, du, eignfnu, nrot = jacobi(vec1,iim,iim)
    330     eignfnu, d = acc(eignfnu,iim)
    331     du, eignfnu = eigen_sort(du,iim,iim)
    332 
    333     return sddu, unsddu, sddv, unsddv
    334 
    335 
    336 def inifilr(iim, xprimu):
    337     """ from: filtrez/filtreg_mod.F90
    338         !    ... H. Upadhyaya, O.Sharma   ...
    339         !
    340         !     version 3 .....
    341 
    342         !     Correction  le 28/10/97    P. Le Van .
    343         !  -------------------------------------------------------------------
    344         !
    345         ! ------------------------------------------------------------
    346         !   This routine computes the eigenfunctions of the laplacien
    347         !   on the stretched grid, and the filtering coefficients
    348         !     
    349         !  We designate:
    350         !   eignfn   eigenfunctions of the discrete laplacien
    351         !   eigenvl  eigenvalues
    352         !   jfiltn   indexof the last scalar line filtered in NH
    353         !   jfilts   index of the first line filtered in SH
    354         !   modfrst  index of the mode from WHERE modes are filtered
    355         !   modemax  maximum number of modes ( im )
    356         !   coefil   filtering coefficients ( lamda_max*COS(rlat)/lamda )
    357         !   sdd      SQRT( dx )
    358         !     
    359         !     the modes are filtered from modfrst to modemax
    360         !     
    361         !-----------------------------------------------------------
    362         !
    363     """
    364     fname = 'inifilr'
    365 
    366     dlonu = xprimu
    367 
    368     !
    369     CALL inifgn(eignvl)
    370     !
    371     PRINT *,'inifilr: EIGNVL '
    372     PRINT 250,eignvl
    373 250 FORMAT( 1x,5e14.6)
    374     !
    375     ! compute eigenvalues and eigenfunctions
    376     !
    377     !
    378     !.................................................................
    379     !
    380     !  compute the filtering coefficients for scalar lines and
    381     !  meridional wind v-lines
    382     !
    383     !  we filter all those latitude lines WHERE coefil < 1
    384     !  NO FILTERING AT POLES
    385     !
    386     !  colat0 is to be used  when alpha (stretching coefficient)
    387     !  is set equal to zero for the regular grid CASE
    388     !
    389     !    .......   Calcul  de  colat0   .........
    390     !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...
    391     !
    392     !
    393     DO j = 1,jjm
    394        dlatu( j ) = rlatu( j ) - rlatu( j+1 )
    395     ENDDO
    396     !
    397 #ifdef CRAY
    398     iymin   = ISMIN( jjm, dlatu, 1 )
    399     ixmineq = ISMIN( iim, dlonu, 1 )
    400     dymin   = dlatu( iymin )
    401     dxmin   = dlonu( ixmineq )
    402 #else
    403     dxmin   =  dlonu(1)
    404     DO  i  = 2, iim
    405        dxmin = MIN( dxmin,dlonu(i) )
    406     ENDDO
    407     dymin  = dlatu(1)
    408     DO j  = 2, jjm
    409        dymin = MIN( dymin,dlatu(j) )
    410     ENDDO
    411 #endif
    412     !
    413     ! For a regular grid, we want the filter to start at latitudes
    414     ! corresponding to lengths dx of the same size as dy (in terms
    415     ! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
    416     !  <=> latitude=60 degrees).
    417     ! Same idea for the zoomed grid: start filtering polewards as soon
    418     ! as length dx becomes of the same size as dy
    419     !
    420     colat0  =  MIN( 0.5, dymin/dxmin )
    421     !
    422     IF( .NOT.fxyhypb.AND.ysinus )  THEN
    423        colat0 = 0.6
    424        !         ...... a revoir  pour  ysinus !   .......
    425        alphax = 0.
    426     ENDIF
    427     !
    428     PRINT 50, colat0,alphax
    429 50  FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
    430     !
    431     IF(alphax.EQ.1. )  THEN
    432        PRINT *,' Inifilr  alphax doit etre  <  a 1.  Corriger '
    433        STOP
    434     ENDIF
    435     !
    436     lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
    437 
    438     !                        ... Correction  le 28/10/97  ( P.Le Van ) ..
    439     !
    440     DO i = 2,iim
    441        rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
    442     ENDDO
    443     !
    444 
    445     DO j = 1,jjm
    446        DO i = 1,iim
    447           coefilu( i,j )  = 0.0
    448           coefilv( i,j )  = 0.0
    449           coefilu2( i,j ) = 0.0
    450           coefilv2( i,j ) = 0.0
    451        ENDDO
    452     ENDDO
    453 
    454     !
    455     !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
    456     !    .........................................................
    457     !
    458     modemax = iim
    459 
    460 !!!!    imx = modemax - 4 * (modemax/iim)
    461 
    462     imx  = iim
    463     !
    464     PRINT *,'inifilr: TRUNCATION AT ',imx
    465     !
    466 ! Ehouarn: set up some defaults
    467     jfiltnu=2 ! avoid north pole
    468     jfiltsu=jjm ! avoid south pole (which is at jjm+1)
    469     jfiltnv=1 ! NB: no poles on the V grid
    470     jfiltsv=jjm
    471 
    472     DO j = 2, jjm/2+1
    473        cof = COS( rlatu(j) )/ colat0
    474        IF ( cof .LT. 1. ) THEN
    475           IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) THEN
    476             jfiltnu= j
    477           ENDIF
    478        ENDIF
    479 
    480        cof = COS( rlatu(jjp1-j+1) )/ colat0
    481        IF ( cof .LT. 1. ) THEN
    482           IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) THEN
    483                jfiltsu= jjp1-j+1
    484           ENDIF
    485        ENDIF
    486     ENDDO
    487     !
    488     DO j = 1, jjm/2
    489        cof = COS( rlatv(j) )/ colat0
    490        IF ( cof .LT. 1. ) THEN
    491           IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) THEN
    492             jfiltnv= j
    493           ENDIF
    494        ENDIF
    495 
    496        cof = COS( rlatv(jjm-j+1) )/ colat0
    497        IF ( cof .LT. 1. ) THEN
    498           IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) THEN
    499                jfiltsv= jjm-j+1
    500           ENDIF
    501        ENDIF
    502     ENDDO
    503     !                                 
    504 
    505     IF( jfiltnu.GT. jjm/2 +1 )  THEN
    506        PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
    507        STOP
    508     ENDIF
    509 
    510     IF( jfiltsu.GT.  jjm  +1 )  THEN
    511        PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
    512        STOP
    513     ENDIF
    514 
    515     IF( jfiltnv.GT. jjm/2    )  THEN
    516        PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
    517        STOP
    518     ENDIF
    519 
    520     IF( jfiltsv.GT.     jjm  )  THEN
    521        PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
    522        STOP
    523     ENDIF
    524 
    525     PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
    526          jfiltnv,jfiltsv,jfiltnu,jfiltsu
    527 
    528     IF(first_call_inifilr) THEN
    529        ALLOCATE(matriceun(iim,iim,jfiltnu))
    530        ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
    531        ALLOCATE(matricevn(iim,iim,jfiltnv))
    532        ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
    533        ALLOCATE( matrinvn(iim,iim,jfiltnu))
    534        ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
    535        first_call_inifilr = .FALSE.
    536     ENDIF
    537 
    538     !                                 
    539     !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
    540     !................................................................
    541     !
    542     !
    543     DO j = 1,jjm
    544     !default initialization: all modes are retained (i.e. no filtering)
    545        modfrstu( j ) = iim
    546        modfrstv( j ) = iim
    547     ENDDO
    548     !
    549     DO j = 2,jfiltnu
    550        DO k = 2,modemax
    551           cof = rlamda(k) * COS( rlatu(j) )
    552           IF ( cof .LT. 1. ) GOTO 82
    553        ENDDO
    554        GOTO 84
    555 82     modfrstu( j ) = k
    556        !
    557        kf = modfrstu( j )
    558        DO k = kf , modemax
    559           cof = rlamda(k) * COS( rlatu(j) )
    560           coefilu(k,j) = cof - 1.
    561           coefilu2(k,j) = cof*cof - 1.
    562        ENDDO
    563 84     CONTINUE
    564     ENDDO
    565     !                                 
    566     !
    567     DO j = 1,jfiltnv
    568        !
    569        DO k = 2,modemax
    570           cof = rlamda(k) * COS( rlatv(j) )
    571           IF ( cof .LT. 1. ) GOTO 87
    572        ENDDO
    573        GOTO 89
    574 87     modfrstv( j ) = k
    575        !
    576        kf = modfrstv( j )
    577        DO k = kf , modemax
    578           cof = rlamda(k) * COS( rlatv(j) )
    579           coefilv(k,j) = cof - 1.
    580           coefilv2(k,j) = cof*cof - 1.
    581        ENDDO
    582 89     CONTINUE
    583     ENDDO
    584     !
    585     DO j = jfiltsu,jjm
    586        DO k = 2,modemax
    587           cof = rlamda(k) * COS( rlatu(j) )
    588           IF ( cof .LT. 1. ) GOTO 92
    589        ENDDO
    590        GOTO 94
    591 92     modfrstu( j ) = k
    592        !
    593        kf = modfrstu( j )
    594        DO k = kf , modemax
    595           cof = rlamda(k) * COS( rlatu(j) )
    596           coefilu(k,j) = cof - 1.
    597           coefilu2(k,j) = cof*cof - 1.
    598        ENDDO
    599 94     CONTINUE
    600     ENDDO
    601     !                                 
    602     DO j = jfiltsv,jjm
    603        DO k = 2,modemax
    604           cof = rlamda(k) * COS( rlatv(j) )
    605           IF ( cof .LT. 1. ) GOTO 97
    606        ENDDO
    607        GOTO 99
    608 97     modfrstv( j ) = k
    609        !
    610        kf = modfrstv( j )
    611        DO k = kf , modemax
    612           cof = rlamda(k) * COS( rlatv(j) )
    613           coefilv(k,j) = cof - 1.
    614           coefilv2(k,j) = cof*cof - 1.
    615        ENDDO
    616 99     CONTINUE
    617     ENDDO
    618     !
    619 
    620     IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
    621 ! Ehouarn: and what are these for??? Trying to handle a limit case
    622 !          where filters extend to and meet at the equator?
    623        IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
    624        IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
    625 
    626        PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' , &
    627             jfiltnv,jfiltsv,jfiltnu,jfiltsu
    628     ENDIF
    629 
    630     PRINT *,'   Modes premiers  v  '
    631     PRINT 334,modfrstv
    632     PRINT *,'   Modes premiers  u  '
    633     PRINT 334,modfrstu
    634 
    635     ! 
    636     !   ...................................................................
    637     !
    638     !   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes
    639     !                       sur la grille scalaire                 ........
    640     !   ...................................................................
    641     !
    642     DO j = 2, jfiltnu
    643 
    644        DO i=1,iim
    645           coff = coefilu(i,j)
    646           IF( i.LT.modfrstu(j) ) coff = 0.
    647           DO k=1,iim
    648              eignft(i,k) = eignfnv(k,i) * coff
    649           ENDDO
    650        ENDDO ! of DO i=1,iim
    651 #ifdef CRAY
    652        CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
    653 #else
    654 #ifdef BLAS
    655        CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
    656             eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
    657 #else
    658        DO k = 1, iim
    659           DO i = 1, iim
    660              matriceun(i,k,j) = 0.0
    661              DO ii = 1, iim
    662                 matriceun(i,k,j) = matriceun(i,k,j) &
    663                      + eignfnv(i,ii)*eignft(ii,k)
    664              ENDDO
    665           ENDDO
    666        ENDDO ! of DO k = 1, iim
    667 #endif
    668 #endif
    669 
    670     ENDDO ! of DO j = 2, jfiltnu
    671 
    672     DO j = jfiltsu, jjm
    673 
    674        DO i=1,iim
    675           coff = coefilu(i,j)
    676           IF( i.LT.modfrstu(j) ) coff = 0.
    677           DO k=1,iim
    678              eignft(i,k) = eignfnv(k,i) * coff
    679           ENDDO
    680        ENDDO ! of DO i=1,iim
    681 #ifdef CRAY
    682        CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
    683 #else
    684 #ifdef BLAS
    685        CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
    686             eignfnv, iim, eignft, iim, 0.0, &
    687             matriceus(1,1,j-jfiltsu+1), iim)
    688 #else
    689        DO k = 1, iim
    690           DO i = 1, iim
    691              matriceus(i,k,j-jfiltsu+1) = 0.0
    692              DO ii = 1, iim
    693                 matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) &
    694                      + eignfnv(i,ii)*eignft(ii,k)
    695              ENDDO
    696           ENDDO
    697        ENDDO ! of DO k = 1, iim
    698 #endif
    699 #endif
    700 
    701     ENDDO ! of DO j = jfiltsu, jjm
    702 
    703     !   ...................................................................
    704     !
    705     !   ... Calcul de la matrice filtre 'matricev'  pour les champs situes
    706     !                       sur la grille   de V ou de Z           ........
    707     !   ...................................................................
    708     !
    709     DO j = 1, jfiltnv
    710 
    711        DO i = 1, iim
    712           coff = coefilv(i,j)
    713           IF( i.LT.modfrstv(j) ) coff = 0.
    714           DO k = 1, iim
    715              eignft(i,k) = eignfnu(k,i) * coff
    716           ENDDO
    717        ENDDO
    718 #ifdef CRAY
    719        CALL MXM( eignfnu,iim,eignft,iim,matricevn(1,1,j),iim )
    720 #else
    721 #ifdef BLAS
    722        CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
    723             eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
    724 #else
    725        DO k = 1, iim
    726           DO i = 1, iim
    727              matricevn(i,k,j) = 0.0
    728              DO ii = 1, iim
    729                 matricevn(i,k,j) = matricevn(i,k,j) &
    730                      + eignfnu(i,ii)*eignft(ii,k)
    731              ENDDO
    732           ENDDO
    733        ENDDO
    734 #endif
    735 #endif
    736 
    737     ENDDO ! of DO j = 1, jfiltnv
    738 
    739     DO j = jfiltsv, jjm
    740 
    741        DO i = 1, iim
    742           coff = coefilv(i,j)
    743           IF( i.LT.modfrstv(j) ) coff = 0.
    744           DO k = 1, iim
    745              eignft(i,k) = eignfnu(k,i) * coff
    746           ENDDO
    747        ENDDO
    748 #ifdef CRAY
    749        CALL MXM(eignfnu,iim,eignft,iim,matricevs(1,1,j-jfiltsv+1),iim)
    750 #else
    751 #ifdef BLAS
    752        CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
    753             eignfnu, iim, eignft, iim, 0.0, &
    754             matricevs(1,1,j-jfiltsv+1), iim)
    755 #else
    756        DO k = 1, iim
    757           DO i = 1, iim
    758              matricevs(i,k,j-jfiltsv+1) = 0.0
    759              DO ii = 1, iim
    760                 matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) &
    761                      + eignfnu(i,ii)*eignft(ii,k)
    762              ENDDO
    763           ENDDO
    764        ENDDO
    765 #endif
    766 #endif
    767 
    768     ENDDO ! of DO j = jfiltsv, jjm
    769 
    770     !   ...................................................................
    771     !
    772     !   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes
    773     !              sur la grille scalaire , pour le filtre inverse ........
    774     !   ...................................................................
    775     !
    776     DO j = 2, jfiltnu
    777 
    778        DO i = 1,iim
    779           coff = coefilu(i,j)/ ( 1. + coefilu(i,j) )
    780           IF( i.LT.modfrstu(j) ) coff = 0.
    781           DO k=1,iim
    782              eignft(i,k) = eignfnv(k,i) * coff
    783           ENDDO
    784        ENDDO
    785 #ifdef CRAY
    786        CALL MXM( eignfnv,iim,eignft,iim,matrinvn(1,1,j),iim )
    787 #else
    788 #ifdef BLAS
    789        CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
    790             eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
    791 #else
    792        DO k = 1, iim
    793           DO i = 1, iim
    794              matrinvn(i,k,j) = 0.0
    795              DO ii = 1, iim
    796                 matrinvn(i,k,j) = matrinvn(i,k,j) &
    797                      + eignfnv(i,ii)*eignft(ii,k)
    798              ENDDO
    799           ENDDO
    800        ENDDO
    801 #endif
    802 #endif
    803 
    804     ENDDO ! of DO j = 2, jfiltnu
    805 
    806     DO j = jfiltsu, jjm
    807 
    808        DO i = 1,iim
    809           coff = coefilu(i,j) / ( 1. + coefilu(i,j) )
    810           IF( i.LT.modfrstu(j) ) coff = 0.
    811           DO k=1,iim
    812              eignft(i,k) = eignfnv(k,i) * coff
    813           ENDDO
    814        ENDDO
    815 #ifdef CRAY
    816        CALL MXM(eignfnv,iim,eignft,iim,matrinvs(1,1,j-jfiltsu+1),iim)
    817 #else
    818 #ifdef BLAS
    819        CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
    820             eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim)
    821 #else
    822        DO k = 1, iim
    823           DO i = 1, iim
    824              matrinvs(i,k,j-jfiltsu+1) = 0.0
    825              DO ii = 1, iim
    826                 matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) &
    827                      + eignfnv(i,ii)*eignft(ii,k)
    828              ENDDO
    829           ENDDO
    830        ENDDO
    831 #endif
    832 #endif
    833 
    834     ENDDO ! of DO j = jfiltsu, jjm
    835 
    836     IF (use_filtre_fft) THEN
    837        CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu,  &
    838                            coefilv,modfrstv,jfiltnv,jfiltsv)
    839        CALL Init_filtre_fft_loc(coefilu,modfrstu,jfiltnu,jfiltsu,  &
    840                            coefilv,modfrstv,jfiltnv,jfiltsv)
    841     ENDIF
    842 
    843     !   ...................................................................
    844 
    845     !
    846 334 FORMAT(1x,24i3)
    847 755 FORMAT(1x,6f10.3,i3)
    848 
    849     RETURN
    850   END SUBROUTINE inifilr
    851 
    852 
    853 def filtreg(dx, dy, champorig, nlat, nbniv, ifiltre, iaire, griscal ,iterv):
    854     """c=======================================================================
    855        c
    856        c   Auteur: P. Le Van        07/10/97
    857        c   ------
    858        c
    859        c   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
    860        c                     pour l'operateur  Filtre    .
    861        c   ------
    862        c
    863        c   Arguments:
    864        c   ----------
    865        c
    866        c      nblat                 nombre de latitudes a filtrer
    867        c      nbniv                 nombre de niveaux verticaux a filtrer
    868        c      champ(iip1,nblat,nbniv)  en entree : champ a filtrer
    869        c                            en sortie : champ filtre
    870        c      ifiltre               +1  Transformee directe
    871        c                            -1  Transformee inverse
    872        c                            +2  Filtre directe
    873        c                            -2  Filtre inverse
    874        c
    875        c      iaire                 1   si champ intensif
    876        c                            2   si champ extensif (pondere par les aires)
    877        c
    878        c      iterv                 1   filtre simple
    879        c
    880        c=======================================================================
    881        c
    882        c
    883        c                      Variable Intensive
    884        c                ifiltre = 1     filtre directe
    885        c                ifiltre =-1     filtre inverse
    886        c
    887        c                      Variable Extensive
    888        c                ifiltre = 2     filtre directe
    889        c                ifiltre =-2     filtre inverse
    890        c
    891     """
    892 
    893     fname = 'filtreg'
    894     first = True
    895 
    896     sdd12 = np.zeros((dx,4), dtype = np.float)
    897 
    898     champ = fortran_mat(champorig)
    899 
    900     sddu, unsddu, sddv, unsddv = inifgn(dx, dy)
    901 
    902     if first:
    903          sdd12[:,type_sddu] = sddu
    904          sdd12[:,type_sddv] = sddv
    905          sdd12[:,type_unsddu] = unsddu
    906          sdd12[:,type_unsddv] = unsddv
    907 
    908          first=False
    909 
    910     if ifiltre  == 1 or ifiltre == -1:
    911         print '  ' + fname + ': Pas de transformee simple dans cette version'
    912         quit()
    913      
    914     if iter == 2:
    915         print errormsg
    916         print '  ' + fname + ': Pas d iteration du filtre dans cette version !'
    917         print ' Utiliser old_filtreg et repasser !'
    918         quit()
    919      
    920     if ifiltre == -2 and not griscal:
    921         print errormsg
    922         print '  ' + fname + ' Cette routine ne calcule le filtre inverse que ' +    \
    923           ' sur la grille des scalaires !'
    924         quit(-1)
    925      
    926     if ifiltre != 2 and ifiltre != - 2:
    927         print errormsg
    928         print '  ' + fname + ': Probleme dans filtreg car ifiltre NE 2 et NE -2' +   \
    929           ' corriger et repasser !'
    930         quit(-1)
    931      
    932     iim = dx
    933     jjm = dy
    934 
    935     iim2 = iim*iim
    936     immjm = iim*jjm
    937 
    938     if griscal:
    939         if nlat != jjp1:
    940             print errormsg
    941             print '  ' + fname  + ': 1111'
    942             quit(-1)
    943         else
    944             if iaire == 1:
    945                 sdd1_type = type_sddv
    946                 sdd2_type = type_unsddv
    947             else:
    948                 sdd1_type = type_unsddv
    949                 sdd2_type = type_sddv
    950 
    951             jdfil1 = 2
    952             jffil1 = jfiltnu
    953             jdfil2 = jfiltsu
    954             jffil2 = jjm
    955     else:
    956         if nlat != jjm:
    957             print errormsg
    958             print '  ' + fname  + ': 2222'
    959             quit(-1)
    960         else:           
    961             if iaire == 1:
    962                 sdd1_type = type_sddu
    963                 sdd2_type = type_unsddu
    964             else:
    965                 sdd1_type = type_unsddu
    966                 sdd2_type = type_sddu
    967 
    968            
    969             jdfil1 = 1
    970             jffil1 = jfiltnv
    971             jdfil2 = jfiltsv
    972             jffil2 = jjm
    973      
    974     for hemisph in range(1,3):
    975         if (hemisph == 1 ):
    976             jdfil = jdfil1
    977             jffil = jffil1
    978         else:
    979             jdfil = jdfil2
    980             jffil = jffil2
    981  
    982         for l in range(nbniv):
    983             for j in rangE(jdfil,jffil):
    984                 for i in range(iim):
    985                     champ[i,j,l] = champ[i,j,l]*sdd12[i,sdd1_type] # sdd1(i)
    986          
    987         if hemisph == 1:
    988             if ifiltre == -2:
    989                for j in range(jdfil,jffil):
    990                   eignq[:,j-jdfil+1,:] = np.matmul(matrinvn[:,:,j], champ[:,j,:])
    991 #endif
    992                END DO
    993                
    994             ELSE IF ( griscal )     THEN
    995                
    996                DO j = jdfil,jffil
    997 #ifdef BLAS
    998                   CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
    999      &                 matriceun(1,1,j),
    1000      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
    1001      &                 eignq(1,j-jdfil+1,1), iim*nlat)
    1002 #else
    1003                   eignq(:,j-jdfil+1,:)
    1004      $                 = matmul(matriceun(:,:,j), champ(:iim,j,:))
    1005 #endif
    1006                END DO
    1007                
    1008             ELSE
    1009                
    1010                DO j = jdfil,jffil
    1011 #ifdef BLAS
    1012                   CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
    1013      &                 matricevn(1,1,j),
    1014      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
    1015      &                 eignq(1,j-jdfil+1,1), iim*nlat)
    1016 #else
    1017                   eignq(:,j-jdfil+1,:)
    1018      $                 = matmul(matricevn(:,:,j), champ(:iim,j,:))
    1019 #endif
    1020                END DO
    1021                
    1022             ENDIF
    1023            
    1024          ELSE
    1025            
    1026             IF( ifiltre. EQ. -2 )   THEN
    1027                
    1028                DO j = jdfil,jffil
    1029 #ifdef BLAS
    1030                   CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
    1031      &                 matrinvs(1,1,j-jfiltsu+1),
    1032      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
    1033      &                 eignq(1,j-jdfil+1,1), iim*nlat)
    1034 #else
    1035                   eignq(:,j-jdfil+1,:)
    1036      $                 = matmul(matrinvs(:,:,j-jfiltsu+1),
    1037      $                 champ(:iim,j,:))
    1038 #endif
    1039                END DO
    1040                
    1041                
    1042             ELSE IF ( griscal )     THEN
    1043                
    1044                DO j = jdfil,jffil
    1045 #ifdef BLAS
    1046                   CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
    1047      &                 matriceus(1,1,j-jfiltsu+1),
    1048      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
    1049      &                 eignq(1,j-jdfil+1,1), iim*nlat)
    1050 #else
    1051                   eignq(:,j-jdfil+1,:)
    1052      $                 = matmul(matriceus(:,:,j-jfiltsu+1),
    1053      $                 champ(:iim,j,:))
    1054 #endif
    1055                END DO
    1056                              
    1057             ELSE
    1058                
    1059                DO j = jdfil,jffil
    1060 #ifdef BLAS
    1061                   CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
    1062      &                 matricevs(1,1,j-jfiltsv+1),
    1063      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
    1064      &                 eignq(1,j-jdfil+1,1), iim*nlat)
    1065 #else
    1066                   eignq(:,j-jdfil+1,:)
    1067      $                 = matmul(matricevs(:,:,j-jfiltsv+1),
    1068      $                 champ(:iim,j,:))
    1069 #endif
    1070                END DO
    1071                              
    1072             ENDIF
    1073            
    1074          ENDIF
    1075          
    1076          IF( ifiltre.EQ. 2 )  THEN
    1077            
    1078             DO l = 1, nbniv
    1079                DO j = jdfil,jffil
    1080                   DO i = 1, iim
    1081                      champ( i,j,l ) =
    1082      &                    (champ(i,j,l) + eignq(i,j-jdfil+1,l))
    1083      &                    * sdd12(i,sdd2_type) ! sdd2(i)
    1084                   END DO
    1085                END DO
    1086             END DO
    1087 
    1088          ELSE
    1089 
    1090             DO l = 1, nbniv
    1091                DO j = jdfil,jffil
    1092                   DO i = 1, iim
    1093                      champ( i,j,l ) =
    1094      &                    (champ(i,j,l) - eignq(i,j-jdfil+1,l))
    1095      &                    * sdd12(i,sdd2_type) ! sdd2(i)
    1096                   END DO
    1097                END DO
    1098             END DO
    1099 
    1100          ENDIF
    1101 
    1102          DO l = 1, nbniv
    1103             DO j = jdfil,jffil
    1104                champ( iip1,j,l ) = champ( 1,j,l )
    1105             END DO
    1106          END DO
    1107 
    1108      
    1109       ENDDO
    1110 
    1111 1111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
    1112      &     filtrer, sur la grille des scalaires'/)
    1113 2222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
    1114      &     ltrer, sur la grille de V ou de Z'/)
    1115       RETURN
    1116       END
    1117 
    111823def SSUM(n,sx,incx):
    111924    """ Obsolete version of sum for non Fortan 90 code
    … …  
    112530    ix = 0
    112631
    1127     for i in range(n)
     32    for i in range(n):
    112833        ssumv=ssumv+sx[ix]
    112934        ix=ix+incx
    … …  
    114550    return sy
    114651
    1147 def exner_hyb (dx, dy, dz, psv, pv, alphav, betav)
     52def exner_hyb (dx, dy, dz, psv, pv, alphav, betav):
    114853    """c
    114954       c     Auteurs :  P.Le Van  , Fr. Hourdin  .
    … …  
    119297        return pksv, pkv, pkfv
    119398
    1194 #!!!! General case:
     99#### General case:
    1195100     
    1196101    unpl2k = 1.+ 2.* kappa
    … …  
    1221126#
    1222127    for l in range (llm-1,1,-1):
    1223 c
     128
    1224129        for ij in range(ngrid):
    1225130            dellta = p[ij,l]* unpl2k + p[ij,l+1]* ( beta[ij,l+1]-unpl2k )
    … …  
    1280185    pkf = np.zeros((dz, dy, dx), dtype=np.float)
    1281186
    1282     ppn = np.zeros((iim), dtype=np.float))
    1283     ppn = np.zeros((iim), dtype=np.float))
     187    ppn = np.zeros((iim), dtype=np.float)
     188    ppn = np.zeros((iim), dtype=np.float)
    1284189
    1285190    firstcall = True
    … …  
    1304209
    1305210#### Specific behaviour for Shallow Water (1 vertical layer) case:
    1306       if llm == 1:
     211    if llm == 1:
    1307212     
    1308213#         Compute pks(:),pk(:),pkf(:)
    1309214       
    1310           for ij in range(ngrid):
    1311               pks[ij] = (cpp/preff) * ps[ij]
    1312               pk[ij,1] = .5*pks[ij]
     215        for ij in range(ngrid):
     216            pks[ij] = (cpp/preff) * ps[ij]
     217            pk[ij,1] = .5*pks[ij]
    1313218
    1314219         
    1315           pkf = SCOPY( ngrid * llm, pk, 1, pkf, 1 )
     220        pkf = SCOPY( ngrid * llm, pk, 1, pkf, 1 )
    1316221# We do not filter for iniaqua
    1317222#        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
    1318223       
    1319224#         our work is done, exit routine
    1320           return pksv, pkv, pkfv
    1321 
     225        return pksv, pkv, pkfv
    1322226
    1323227#### General case:
    … …  
    1330234        pks[ij] = cpp * ( ps[ij]/preff ) ** kappa
    1331235
    1332     for ij in range(iim)
     236    for ij in range(iim):
    1333237        ppn[ij] = aire[ij] * pks[ij]
    1334238        pps[ij] = aire[ij+ip1jm] * pks[ij+ip1jm]
    … …  
    1337241    xps = SSUM(iim,pps,1) /apols
    1338242
    1339     for ij in range(iip1)
     243    for ij in range(iip1):
    1340244        pks[ij] = xpn
    1341245        pks[ij+ip1jm] = xps
    … …  
    1373277    return pksv, pkv, pkfv
    1374278
    1375 def pression(dx, dy, dz, apv, bpv, psv)
     279def pression(dx, dy, dz, apv, bpv, psv):
    1376280    """c
    1377281
    … …  
    1881785            masse[ij,l] = airesurg[ij] * ( p[ij,l] - p[ij,l+1] )
    1882786
    1883         for ij in range(ip1jmp1,iip1)
     787        for ij in range(ip1jmp1,iip1):
    1884788            masse[ij+iim,l] = masse[ij,l]
    1885789       
    1886790    return masse
    1887791
    1888 def geopot(ngrid, teta, pk, pks)
     792def geopot(ngrid, teta, pk, pks):
    1889793    """c=======================================================================
    1890794       c
    … …  
    1975879    u = np.zeros((iip1,jjm,llm), dtype=np.float)
    1976880
    1977 
    1978     for j in range(jjm)
    1979 
     881    for j in range(jjm):
    1980882        if np.abs(np.sin(rlatv[j])) < 1.e-4:
    1981883            zlat = 1.e-4
    … …  
    2051953    IX3 = np.mod(IA3*IX3+IC3,M3)
    2052954    J = 1+(Nvals*IX3)/M3
    2053     if J > Nvals or J lt 1: quit()
     955    if J > Nvals or J < 1: quit()
    2054956    ran1=R[J]
    2055957    R[J]=(np.float(IX1)+np.float(IX2)*RM2)*RM1
    … …  
    2067969parser.add_option("-d", "--dimensions", dest="dims",
    2068970  help="dimensions of the initial conditions: dimx,dimy,dimz", metavar="VALUES")
    2069 parser.add_option("-p", "--pressure_exner", dest="pexner,
     971parser.add_option("-p", "--pressure_exner", dest="pexner",
    2070972  help="how as to b computed Exner pressure ('hybrid': hybrid coordinates, 'middle': middle layer)", 
    2071973  metavar="VALUE")
Note: See TracChangeset for help on using the changeset viewer.