source: branches/2011/dev_MERCATOR_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DIA/surdetermine.F90 @ 3047

MODULE surdetermine 

   USE in_out_manager  ! I/O units

   IMPLICIT NONE
   PUBLIC
   
   INTEGER, PARAMETER :: jpincomax = 18
   INTEGER, PARAMETER :: jpdimsparse = jpincomax*300*24

   INTEGER, PUBLIC :: nsparse, ninco
   REAL(wp), PUBLIC, DIMENSION(jpdimsparse) :: valuesparse
   INTEGER , PUBLIC, DIMENSION(jpdimsparse) :: njsparse, nisparse

   INTEGER, SAVE, DIMENSION(jpincomax) :: ipos1
   REAL(wp), DIMENSION(jpincomax) :: tmp4, tmp7
   REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3, zpilier
   REAL(wp), SAVE, DIMENSION(jpincomax) :: zpivot

   !!---------------------------------------------------------------------------------
   !!
   !!---------------------------------------------------------------------------------

CONTAINS

      SUBROUTINE SUR_DETERMINE(init)
 
      INTEGER, INTENT(in) :: init

      INTEGER  :: &
         ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
      REAL(wp) :: zval1, zval2, zx1
      REAL(wp), DIMENSION(jpincomax) :: ztmpx, zcol1, zcol2
      INTEGER, DIMENSION(jpincomax) :: ipos2, ipivot      
   !---------------------------------------------------------------------------------

      IF (init==1) THEN
        IF(nsparse .gt. jpdimsparse) STOP 'surdetermine erreur1'
        IF(ninco .gt. jpincomax)THEN
          IF (lwp) WRITE(numout,*)'ninco   =',ninco
          IF (lwp) WRITE(numout,*)'jpincomax=',jpincomax
          STOP 'DONC dans surdetermine erreur2'
        END IF

        ztmp3(:,:)=0.e0

        DO jk1_sd=1,nsparse
          DO jk2_sd=1,nsparse 
            nisparse(jk2_sd)=nisparse(jk2_sd)
            njsparse(jk2_sd)=njsparse(jk2_sd)
            IF(nisparse(jk2_sd) .eq. nisparse(jk1_sd)) &
            ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))= &
            ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  &
           +valuesparse(jk1_sd)*valuesparse(jk2_sd)
          END DO
        END DO

        DO jj_sd=1,ninco
          ipos1(jj_sd)=jj_sd
          ipos2(jj_sd)=jj_sd
        ENDDO

        DO ji_sd=1,ninco
! recherche du plus grand pivot non nul
          zval1=ABS(ztmp3(ji_sd,ji_sd))

          ipivot(ji_sd)=ji_sd
          DO jj_sd=ji_sd,ninco
             zval2=ABS(ztmp3(ji_sd,jj_sd))
             IF(zval2.GE.zval1) THEN
              ipivot(ji_sd)=jj_sd
              zval1=zval2
             ENDIF
          END DO

          DO ji1_sd=1,ninco
             zcol1(ji1_sd)=ztmp3(ji1_sd,ji_sd)
             zcol2(ji1_sd)=ztmp3(ji1_sd,ipivot(ji_sd))
             ztmp3(ji1_sd,ji_sd)=zcol2(ji1_sd)
             ztmp3(ji1_sd,ipivot(ji_sd))=zcol1(ji1_sd)
          END DO

          ipos2(ji_sd)=ipos1(ipivot(ji_sd))
          ipos2(ipivot(ji_sd))=ipos1(ji_sd)
          ipos1(ji_sd)=ipos2(ji_sd)
          ipos1(ipivot(ji_sd))=ipos2(ipivot(ji_sd))


!-------------------------------
          zpivot(ji_sd)=ztmp3(ji_sd,ji_sd)
          DO jj_sd=1,ninco
             ztmp3(ji_sd,jj_sd)=ztmp3(ji_sd,jj_sd)/zpivot(ji_sd)
          END DO
!-------------------------------

!-------------------------------
          DO ji2_sd=ji_sd+1,ninco
             zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
             DO jj_sd=1,ninco
               ztmp3(ji2_sd,jj_sd)= &
               ztmp3(ji2_sd,jj_sd)-ztmp3(ji_sd,jj_sd)*zpilier(ji2_sd,ji_sd)
             END DO 
          END DO

        END DO

      ENDIF ! End init==1

!
      DO ji_sd=1,ninco
        tmp4(ji_sd)=tmp4(ji_sd)/zpivot(ji_sd)
        DO ji2_sd=ji_sd+1,ninco
           tmp4(ji2_sd)=tmp4(ji2_sd)-tmp4(ji_sd)*zpilier(ji2_sd,ji_sd)
        END DO
      END DO

!  resolution du systeme:
      ztmpx(ninco)=tmp4(ninco)/ztmp3(ninco,ninco)
      ji_sd=ninco
      DO ji_sd=ninco-1,1,-1
        zx1=0.
        DO jj_sd=ji_sd+1,ninco
          zx1=zx1+ztmpx(jj_sd)*ztmp3(ji_sd,jj_sd)
        END DO 
        ztmpx(ji_sd)=tmp4(ji_sd)-zx1
      END DO 

      DO jj_sd=1,ninco
        tmp7(ipos1(jj_sd))=ztmpx(jj_sd)
      END DO


  END SUBROUTINE SUR_DETERMINE 

END MODULE surdetermine