wrfv2_fire/var/da/da_recursive_filter/da_mat_cv3.f90

   1 module da_mat_cv3
   2
   3 contains
   4
   5       SUBROUTINE ZERM(A,MI,MJ,NA) ! Set the elements of general matrix to zero
   6       INTEGER MI, MJ, NA
   7       DIMENSION A(NA,*)
   8       DO J=1,MJ
   9        CALL ZERV(A(1,J),MI)
  10       ENDDO
  11       RETURN
  12       END SUBROUTINE ZERM
  13
  14       SUBROUTINE ZERV(D,M)    ! set elements of a vector to zero
  15       INTEGER M, I
  16       DIMENSION D(M)
  17       DO I=1,M
  18        D(I)=0.
  19       ENDDO
  20       RETURN
  21       END SUBROUTINE ZERV
  22
  23       SUBROUTINE MULMV(A,D,E,MI,MJ,NA)
  24       INTEGER MI, MJ, NA, J
  25       DIMENSION A(NA,*),D(*),E(*)
  26       CALL ZERV(E,MJ)
  27       DO J=1,MJ
  28        CALL MADVS(A(1,J),D(J),E,MI)
  29       ENDDO
  30       RETURN
  31       END SUBROUTINE MULMV
  32
  33       SUBROUTINE MADVS(D,S,E,M)
  34       INTEGER M, I
  35       DIMENSION D(*),E(*)
  36       DO I=1,M
  37        E(I)=E(I)+D(I)*S
  38       ENDDO
  39       RETURN
  40       END SUBROUTINE MADVS
  41
  42      subroutine LINMM(a,b,m,mm,na,nb)
  43       DIMENSION A(NA,*),B(NB,*),ipiv(m)
  44       CALL LDUM(A,IPIV,D,M,NA)
  45       CALL UDLMM(A,B,IPIV,M,MM,NA,NB)
  46       RETURN
  47       END subroutine LINMM
  48
  49 !------------------------------------------------------------------------------
  50 !   R.J.Purser, NCEP, Washington D.C.   1996
  51 !                   SUBROUTINE  LDUM
  52 !  perform l-d-u decomposition of square matrix a in place with
  53 !  IMPLICIT pivoting
  54 !
  55 !  --> a    square matrix to be factorized
  56 !  <-- ipiv array encoding the pivoting sequence
  57 !  <-- d    indicator for possible sign change of determinant
  58 !  --> m    degree of (active part of) a
  59 !  --> na   first fortran dimension of a
  60 !
  61 !   LIMITATION:
  62 !       S is an array, internal to this routine, containing the
  63 !    scaling factors of each row used for pivoting decisions. It is given a
  64 !    fortran dimension of NN=500 in the parameter statement below.
  65 !    If the order of the linear system exceeds NN, increase NN.
  66 !------------------------------------------------------------------------------
  67       SUBROUTINE LDUM(A,IPIV,D,M,NA)
  68       PARAMETER(NN=500)
  69       DIMENSION A(NA,*),IPIV(*),S(NN)
  70 !      IF(M.GT.NN)STOP'MATRIX TOO LARGE FOR LDUM'
  71       IF(M.GT.NN)STOP
  72       DO I=1,M
  73        AAM=0.
  74        DO J=1,M
  75         AA=ABS(A(I,J))
  76         IF(AA.GT.AAM)AAM=AA
  77        ENDDO
  78        IF(AAM.EQ.0.)THEN
  79         PRINT'('' ROW '',I3,'' OF MATRIX IN LUFM VANISHES'')',I
  80         STOP
  81        ENDIF
  82        S(I)=1./AAM
  83       ENDDO
  84       D=1.
  85       IPIV(M)=M
  86       DO J=1,M-1
  87        JP=J+1
  88        ABIG=S(J)*ABS(A(J,J))
  89        IBIG=J
  90        DO I=JP,M
  91         AA=S(I)*ABS(A(I,J))
  92         IF(AA.GT.ABIG)THEN
  93          IBIG=I
  94          ABIG=AA
  95         ENDIF
  96        ENDDO
  97 !  swap rows, recording changed sign of determinant
  98        IPIV(J)=IBIG
  99        IF(IBIG.NE.J)THEN
 100         D=-D
 101         DO K=1,M
 102          T=A(J,K)
 103          A(J,K)=A(IBIG,K)
 104          A(IBIG,K)=T
 105         ENDDO
 106         S(IBIG)=S(J)
 107        ENDIF
 108        AJJ=A(J,J)
 109        IF(AJJ.EQ.0.)THEN
 110                            JM=J-1
 111   PRINT'('' FAILURE IN LDUM:''/'' MATRIX SINGULAR, RANK='',i3)',JM
 112                            STOP
 113        ENDIF
 114        AJJI=1./AJJ
 115        DO I=JP,M
 116         AIJ=AJJI*A(I,J)
 117         A(I,J)=AIJ
 118         DO K=JP,M
 119          A(I,K)=A(I,K)-AIJ*A(J,K)
 120         ENDDO
 121        ENDDO
 122       ENDDO
 123       RETURN
 124       END SUBROUTINE LDUM
 125
 126 !------------------------------------------------------------------------------
 127 !   R.J.Purser, National Meteorological Center, Washington D.C.  1993
 128 !                   SUBROUTINE UDLMM
 129 !  use l-u factors in a to back-substitute for mm rhs in b, using ipiv to
 130 !  define the pivoting permutation used in the l-u decomposition.
 131 !
 132 !  --> A    L-D-U factorization of linear system matrux
 133 !  <-> B    right-hand-sides on entry, corresponding matrix of solution
 134 !           vectors on return
 135 !  --> IPIV array encoding the pivoting sequence
 136 !  --> M    degree of (active part of) B and A
 137 !  --> MM   number of right-hand-side vectors (active columns of B)
 138 !  --> NA   first fortran dimension of A
 139 !  --> NB   first fortran dimension of B
 140 !------------------------------------------------------------------------------
 141       SUBROUTINE UDLMM(A,B,IPIV,M,MM,NA,NB)
 142       DIMENSION A(NA,*),B(NB,*),IPIV(*)
 143       DO K=1,MM !loop over columns of B
 144        DO I=1,M
 145         L=IPIV(I)
 146         S=B(L,K)
 147         B(L,K)=B(I,K)
 148         CALL DSBVR(B(1,K),A(I,1),S,I-1,NA)
 149         B(I,K)=S
 150        ENDDO
 151        B(M,K)=B(M,K)/A(M,M)
 152        DO I=M-1,1,-1
 153         AIII=1./A(I,I)
 154         CALL DSBVR(B(I+1,K),A(I,I+1),B(I,K),M-I,NA)
 155         B(I,K)=B(I,K)*AIII
 156        ENDDO
 157       ENDDO
 158       RETURN
 159       END SUBROUTINE UDLMM
 160
 161       subroutine DSBVR(D,A,S,M,NA)
 162       DIMENSION D(M),A(NA,*)
 163       DO I=1,M
 164        S=S-D(I)*A(1,I)
 165       ENDDO
 166       RETURN
 167       END subroutine DSBVR
 168
 169
 170 end module da_mat_cv3