program alfanew

! **********************************************************************
! *              .      .    .
! * PROGRAM:    alfanew
! * PRGMMR:     M. ZUPANSKI     ORG: CIRA/CSU   DATE: 2001-06-07
! *
! * ABSTRACT:   THIS PROGRAM UPDATES THE STEP-LENGTH
! *
! * PROGRAM LOG:
! *
! * 06/07/2001    ..... M. ZUPANSKI:
! * 09/10/2003    ..... M. ZUPANSKI: Upgrade for use with ensemble DA
! *
! **********************************************************************

  integer,parameter::ialfa_old=21    ! old step-length file
  integer,parameter::icost_b=22      ! background cost function 
  integer,parameter::icost_obs=23    ! observation cost function 
  integer,parameter::icost_pen=24    ! GW penalty cost function 
  integer,parameter::ialfa=51        ! new step-length file

!-----
 real ALFA,ALFA_OLD
 real COST_b_top,COST_obs_top,COST_pen_top
 real COST_b_bot,COST_obs_bot,COST_pen_bot
 real COST_top,COST_bot

 real RATIO_max
 integer iout,ioutm1,index

 character*7 alfaname

!
!      DECLARE NAMELIST
!
      NAMELIST /CURRENT_ITERATION/ iout,ioutm1,index
      NAMELIST /STEP/ RATIO_max

!=========================================================
!==============start calculation===================
        write(3,*) "start alfanew"
        write(*,*) "start alfanew"

!-- initialize step-length sum contributions 

       COST_top=0.0
       COST_bot=0.0

!-- read current iteration

       REWIND 13
       READ(13,CURRENT_ITERATION)

       write(*,*) "alfanew: iout,ioutm1,index=",iout,ioutm1,index

!-- read RATIO_max

       REWIND 14
       READ(14,STEP)

       write(*,*) "alfanew: RATIO_max=",RATIO_max

!    read old step-length 
!  (Note: the index is already changed to current in 'alfa_initlze')

      write(alfaname,100) iout
 100  format('alfa_',i2.2)

      close(ialfa)
      open(UNIT=ialfa,FILE=alfaname,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' ialfa OPEN UNIT ERROR IER=',IER
      read(ialfa) ALFA_OLD
      CLOSE(ialfa,status='KEEP')

      WRITE(3,*)"OLD STEP-LENGTH  ",ALFA_OLD
      WRITE(*,*)"OLD STEP-LENGTH  ",ALFA_OLD

!-- background cost function

      close(icost_b)
      open(icost_b,file='backsum',status='unknown',form='unformatted')
      IF(IER.NE.0) WRITE(3,*)' icost_b OPEN UNIT ERROR IER=',IER
      read(icost_b) COST_b_top,COST_b_bot
      close(icost_b,status='KEEP')

      write(*,*) "Background: top, bot =",COST_b_top,COST_b_bot
      write(3,*) "Background: top, bot =",COST_b_top,COST_b_bot
!--------------------------------------------------------------
!-- this is an old fix from 4DVAR, however not needed with EnKF
!!!    if(COST_b_top.eq.0.) COST_b_bot=0.0
!---------------------------------------------------------------

      COST_top = COST_top + COST_b_top 
      COST_bot = COST_bot + COST_b_bot 

!-- observation cost function

      close(icost_obs)
      open(icost_obs,file='obssum',status='unknown',form='unformatted')
      IF(IER.NE.0) WRITE(3,*)' icost_obs OPEN UNIT ERROR IER=',IER
      read(icost_obs) COST_obs_top,COST_obs_bot
      close(icost_obs,status='KEEP')
       if(COST_obs_top.eq.0.) COST_obs_bot=0.0

      write(*,*) "Observation: top, bot =",COST_obs_top,COST_obs_bot
      write(3,*) "Observation: top, bot =",COST_obs_top,COST_obs_bot

      COST_top = COST_top + COST_obs_top 
      COST_bot = COST_bot + COST_obs_bot 

!-- gravity-wave penalty (high-frewuency) cost function

!-- (no gravity wave penalty at this time)

      COST_pen_top=0.0
      COST_pen_bot=0.0

      write(*,*) "GW penalty: top, bot =",COST_pen_top,COST_pen_bot
      write(3,*) "GW penalty: top, bot =",COST_pen_top,COST_pen_bot
       if(COST_pen_top.eq.0.) COST_pen_bot=0.0

      COST_top = COST_top + COST_pen_top 
      COST_bot = COST_bot + COST_pen_bot 

        WRITE(3,*)"ALFANEW: SSBOT = ",COST_bot
        WRITE(3,*)"ALFANEW: SSTOP = ",COST_top
        WRITE(*,*)"ALFANEW: SSBOT = ",COST_bot
        WRITE(*,*)"ALFANEW: SSTOP = ",COST_top

!---   calculate new alfa

      IF(COST_bot.LE.0.0) THEN
        WRITE(3,*)"PROBLEM: SSBOT = ",COST_bot
        WRITE(*,*)"PROBLEM: SSBOT = ",COST_bot
        RATIO = 1.
      ELSE
        RATIO = COST_top / COST_bot
      END IF

!-- safeguard step length

        write(*,*) "before update: RATIO=",RATIO

        if(abs(RATIO).gt.RATIO_max) then
         RATIO=sign(1.,RATIO)*RATIO_max
        end if

        write(*,*) "after update: RATIO=",RATIO

!-- new step length

        ALFA = ALFA_OLD * RATIO


      WRITE(3,*)"NEW STEP-LENGTH  ",ALFA
      WRITE(*,*)"NEW STEP-LENGTH  ",ALFA
!---

!--    write new step-length

      write(alfaname,100) iout

      close(ialfa)
      open(UNIT=ialfa,FILE=alfaname,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' ialfa OPEN UNIT ERROR IER=',IER
      write(ialfa) ALFA
      CLOSE(ialfa,status='KEEP')

      close(11)
      open(UNIT=11,FILE='alfa.out',FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' 11 OPEN UNIT ERROR IER=',IER
      write(11,200) iout,ALFA_OLD,ALFA
      CLOSE(11,status='KEEP')

 200  format('iter=',I4,' alfa_old=',E15.5,' alfa=',E15.5)
!====================================================================

 stop
 end