source: trunk/LMDZ.MARS/libf/phymars/elimin_mz1d.F @ 426

c       ************************************************************************
        subroutine elimin_mz1d (c,vc, ilayer,nan,itblout, nw)

c       Eliminate anomalous negative numbers in c(nl,nl) according to "nan":

c        nan = 0 -> no eliminate 
c              @       -> eliminate all numbers with absol.value<abs(max(c(n,r)))/300.
c              2 -> eliminate all anomalous negative numbers in c(n,r).
c              3 -> eliminate all anomalous negative numbers far from the main
c                   diagonal. 
c              8 -> eliminate all non-zero numbers outside the main diagonal,
c                   and the contibution of lower boundary. 
c              9 -> eliminate all non-zero numbers outside the main diagonal. 
c              4 -> hace un smoothing cuando la distancia de separacion entre
c                   el valor maximo y el minimo de cf > 50 capas.
c              5 -> elimina valores menores que 1.0d-19
c              6 -> incluye los dos casos 4 y 5 
c              7 -> llama a lisa: smooth con width=nw & elimina mejorado
c              78-> incluye los dos casos 7 y 8
c              79-> incluye los dos casos 7 y 9

c       itblout (itableout in calling program) is the option for writing 
c       out or not the purged c(n,r) matrix: 
c       itblout = 0 -> no write  
c               = 1 -> write out in curtis***.out according to ilayer

c       ilayer is the index for the layer selected to write out the matrix:
c       ilayer = 0  => matrix elements written out cover all the altitudes
c                                                     with 5 layers steps
c              > 0  =>   "        "      "      "  are  c(ilayer,*)
c       NOTA: 
c         EXISTE LA POSIBILIDAD DE SACAR TODAS LAS CAPAS (TODA LA MATRIZ)
c         UTILIZANDO itableout=30 EN MZTUD

c       jul 2011        malv+fgg       adapted to LMD-MGCM
c       Sep-04          FGG+MALV        Correct include and call parameters
c       cristina        25-sept-1996   y  27-ene-1997
c       JAN 98          MALV            Version for mz1d
c       ************************************************************************

        implicit none

        include 'nltedefs.h'

        integer   nan,j,i,itblout,kk,k,ir,in
        integer   ilayer,jmin, jmax,np,nw,ntimes,ntimes2
!*      real*8    c(nl,nl), vc(nl), amax, cmax, cmin, cs(nl,nl), mini
        real*8    c(nl,nl), vc(nl), amax, cmax, cmin, mini
        real*8 aux(nl), auxs(nl)
        character layercode*3

        ntimes=0
        ntimes2=0
!       type *,'from elimin_mz4: nan, nw',nan,nw

        if (nan .eq. 0) goto 200 

        if(nan.eq.1)then
                do i=1,nl
                        amax=1.0d-36
                        do j=1,nl
                          if(abs(c(i,j)).gt.amax)amax=abs(c(i,j))
                        end do
                        do j=1,nl
                          if(abs(c(i,j)).lt.amax/300.0d0)c(i,j)=0.0d0
                        end do
                enddo
        elseif(nan.eq.2)then
                do i=1,nl
                        do j=1,nl
                          if( ( j.le.(i-2) .or. j.gt.(i+2) ).and.
     @                    ( c(i,j).lt.0.0d0 ) ) c(i,j)=0.0d0
                        end do
                enddo
        elseif(nan.eq.3)then
                do i=1,nl
                        do j=1,nl
                            if (abs(i-j).ge.10) c(i,j)=0.0d0
                        end do
                enddo
        elseif(nan.eq.8)then
                do i=1,nl
                        do j=1,i-1 
                                c(i,j)=0.0d0
                        enddo
                        do j=i+1,nl 
                                c(i,j)=0.0d0
                        enddo
                        vc(i)= 0.d0
                enddo
        elseif(nan.eq.9)then
                do i=1,nl
                        do j=1,i-1 
                                c(i,j)=0.0d0
                        enddo
                        do j=i+1,nl 
                                c(i,j)=0.0d0
                        enddo
                enddo
!       elseif(nan.eq.7.or.nan.eq.78.or.nan.eq.79)then
!               call lisa(c, vc, nl, nw)
        end if
        if(nan.eq.78)then
                do i=1,nl
                        do j=1,i-1 
                                c(i,j)=0.0d0
                        enddo
                        do j=i+1,nl 
                                c(i,j)=0.0d0
                        enddo
                        vc(i)= 0.d0
                enddo
        endif
        if(nan.eq.79)then
                do i=1,nl
                        do j=1,i-1 
                                c(i,j)=0.0d0
                        enddo
                        do j=i+1,nl 
                                c(i,j)=0.0d0
                        enddo
                enddo
        endif

        if(nan.eq.5.or.nan.eq.6)then
                do i=1,nl
                        mini = 1.0d-19
                        do j=1,nl
                         if(abs(c(i,j)).le.mini.and.c(i,j).ne.0.d0) then
                            ntimes2=ntimes2+1
                         end if
                           if ( abs(c(i,j)).le.mini) c(i,j)=0.d0
                        end do
                enddo
        end if

        if(nan.eq.4.or.nan.eq.6)then
                do i=1,nl
                        do j=1,nl
                          aux(j)=c(i,j)
                          auxs(j)=c(i,j)
                        end do
                        call maxdp_2(aux,nl,cmax,jmax)
                           if(abs(jmax-i).ge.50) then
                                call smooth_cf(aux,auxs,i,nl,3)
                                !!!call smooth_cf(aux,auxs,i,nl,5)
                                ntimes=ntimes+1
                           end if
                        do j=1,nl
                          c(i,j)=auxs(j)
                        end do
                end do
        end if

!          type *, 'elimin_mz4: c(n,r) procesed for elimination. '
!          type *, ' '
!          if(nan.eq.4.or.nan.eq.6) type *, '    call smoothing:',ntimes
!          if(nan.eq.5.or.nan.eq.6) type *, '    call elimina:  ',ntimes2
!          if(nan.eq.7)   type *, '    from elimin: lisa w=',nw
!          type *, ' '


 200    continue

c       writting out of c(n,r) in ascii file 

!       if(itblout.eq.1) then

!         if (ilayer.eq.0) then

!          open (unit=2, status='new', 
!     @    file=dircurtis//'curtis_gnu.out', recl=1024)
!           write(2,'(a)') 
!     @    ' curtis matrix:     table with   1.e+7 * acf(n,r) '
!           write(2,114) 'n,r', ( i, i=nl,1,-5 )
!           do in=nl,1,-5
!             write(2,*)
!             write(2,115) in, ( c(in,ir)*1.d7, ir=nl,1,-5 )
!           end do 
!          close(2)


!          write (*,*)  ' '
!          write (*,*)  '  curtis.out has been created. '
!          write (*,*)  ' '

!         else

!            write (layercode,132) ilayer
!           open (2, status='new', 
!     @    file=dircurtis//'curtis'//layercode//'.out')
!           write(2,'(a)') 
!     @    ' curtis matrix:     table with   1.e+7 * acf(n,r) '
!           write(2,116) ' layer x       c(',layercode,
!     @    ',x)           c(x,', layercode,')' 
!           do in=nl,1,-1
!            if (c(ilayer,ilayer).ne.0.d0) then 
!             write(2,117) in, c(ilayer,in), c(in,ilayer), 
!     @        c(ilayer,in)/c(ilayer,ilayer),
!     @        c(in,ilayer)/c(ilayer,ilayer)
!            else
!             write(2,118) in, c(ilayer,in), c(in,ilayer)
!            end if
!           end do 
!           close(2)
!           write (*,*)  ' '
!           write (*,*)  dircurtis//'curtis'//layercode//'.out', 
!     @ ' has been created.'
!           write (*,*) ' '

!         end if

!       elseif(itblout.eq.0)then

!         continue

!       else

!         write (*,*) ' error from elimin: ', 
!     @      ' itblout should be 1 or 0;   itblout= ',itblout
!         stop

!       end if
        
        return

112     format(10x,10(i3,9x))
113     format(1x,i3,2x,9(1pe9.2,2x))

114     format(1x,a3, 11(8x,i3))
115     format( 1x,i3, 2x, 11(1pe10.3))
116     format( 1x,a17,a2,a18,a2,a1 ) 
117     format( 3x,i3, 4(8x,1pe10.3) )
118     format( 3x,i3, 2(8x,1pe10.3) )
120     format( 1x,i3, 1x,i3, 2x, 11(1pe10.3))

132     format(i3)

!  cambio: los formatos 114, 115 , 117 y 118
!  cambio: al cambia nl de 51 a 140 hay que cambiar el formato i2-->i3
!          y ahora en vez de 11 capas de 5 en 5, hay 28
!
        end
c**************************************************************************
        subroutine smooth_cf( c, cs, i, nl, w )
c       hace un smoothing de c(i,*), de la contribucion de todas las capas
c       menos de la capa en cuestion, la i.
c       opcion w (width): el tamanho de la ventana del smoothing.
c       output values: cs
c**************************************************************************

        implicit none

        integer  j,np,i,nl,w
        real*8   c(nl), cs(nl)

        if(w.eq.0) then
         do j=1,nl
                cs(j)=c(j)
         end do

        elseif(w.eq.3) then

!       write (*,*) 'smoothing w=3'
         do j=1,i-4
                   if(j.eq.1) then
                        cs(j)=c(j)
                   else
                        cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1))
                   end if
         end do
         do j=i+4,nl-1
                   if(j.eq.nl) then
                        cs(j)=c(j)
                   else
                        cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1))
                   end if
         end do
        elseif(w.eq.5) then

!       type *,'smoothing w=5'
         do j=3,i-4
                   if(j.eq.1) then
                        cs(j)=c(j)
                   else
                        cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2))
                   end if
         end do
         do j=i+4,nl-2
                   if(j.eq.nl) then
                        cs(j)=c(j)
                   else
                        cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2))
                   end if
         end do
        end if
        return
        end