MODULE OPT

      USE PARAMS
      USE DATAFILES
      USE FRWDMDL
      USE SYNSPEC
      USE MATRIX
      USE WRITEOUT
      USE BANDPARAM

      IMPLICIT NONE

      LOGICAL                             :: Tflg
      LOGICAL                             :: ALL_SPEC_OUT

      REAL(DOUBLE), DIMENSION(:),   ALLOCATABLE :: SE
      REAL(DOUBLE), DIMENSION(:,:), ALLOCATABLE :: SA, SAINV, SHAT
      REAL(DOUBLE), DIMENSION(:,:), ALLOCATABLE :: KS, KSK, SPKSKINV
      REAL(DOUBLE), DIMENSION(:,:), ALLOCATABLE :: KHAT
      REAL(DOUBLE), DIMENSION(:,:), ALLOCATABLE :: G_LM

      LOGICAL        :: LM ! TRUE FOR LEVENBERG MARQUARDT
      INTEGER        :: ILEVENBERG
      REAL (DOUBLE)  :: GAMMA_START, GAMMA_INC, GAMMA_DEC, STOP_CRITERION
      REAL (DOUBLE)  :: CONVERGENCE

      CONTAINS

      !SUBROUTINE OPT_3(Y, XA, XHAT, YHAT, M, N, CONVERGE, MAXITER, TOL, DEBUG, RETFLG, DIVWARN, ITER, SNR, ISMIX, NLEV  )

      SUBROUTINE OPT_3(Y, XA, XHAT, YHAT, M, N, CONVERGE, MAXITER, TOL, DEBUG, RETFLG, DIVWARN, ITER, ISMIX, NLEV  )

! 17SEP02
!    - output format of kfile 2000e16.8 - compatible w/ idl
!      32676 chars per line

!
!
! Version 2; 1/17/91; bjc
! version x; 6/2/95; bjc
! Version 3; 10/12/95; bjc
! Retrieval/forward model switch (RetFlg) added
! Version 4; 8/13/98; psm
! DivWarn and Iter passed out, & WrtK option added
!
! KFlg tells FM whether to calculate weighting functions - BJC 12/96
!
! change from version 1: Se is now a vector; a new matrix arithmetic
! subroutine (MultDiag) has been added
!
!

      LOGICAL       :: KFLG, FILOPEN
      LOGICAL       :: CONVERGE, DEBUG, RETFLG, DIVWARN
      INTEGER       :: N, M, ISMIX, MAXITER, ITER, NAERR
      REAL(DOUBLE)  :: TOL !, SNR
      REAL(DOUBLE)  :: Y(M)
      REAL(DOUBLE)  :: XA(N) ! PARM on input
      REAL(DOUBLE)  :: XHAT(NMAX)
      REAL(DOUBLE)  :: YHAT(MMAX)

      REAL(DOUBLE)   :: YN_OLD(M), SEINVDY(M), SEINVDY_OLD_SE(M)
      REAL(DOUBLE)   :: XN_OLD(N), SAINVDX(N)
      REAL(DOUBLE)   :: GAMMA, RED_GAMMA, INC_GAMMA
      REAL(DOUBLE)   :: KN_OLD(M*N), CHI_2, CHI_2_X, CHI_2_Y, CHI_2_Y_OLD_SE
      REAL(DOUBLE)   :: CHI_2_OLD, D_CHI_2, CHI_2_OLD_SE, D_CHI_2_OLD_SE

      INTEGER      :: I, J, ONE, SGN, KK, NS, IYDX1, IYDX2, JSCAN, IBAND, NLEV
      REAL(DOUBLE) :: RMSDELY, RMSDY, RMSIM1, CHGY, UNCY, SQRMS, VQRMS, RMSKDX !, SNR

      REAL(DOUBLE), DIMENSION(M) :: SEINV, SEINV_OLD, DY, DELY, DYMKDX
      REAL(DOUBLE), DIMENSION(N) :: XNP1, XN, DX, DX_OLD, KSDYMKDX, DELX

      ! need two extra vectors of size N for Levenberg Marquardt -- mp
      REAL(DOUBLE), DIMENSION(N)       :: KSDYMKDX_LM, GSAINVDX, G
      REAL(DOUBLE), DIMENSION(N,N)     :: GSAINV

      ! for calculating AVK when LM
      REAL(DOUBLE), DIMENSION(N,N)     :: IDNN, T2, T3, T4
      REAL(DOUBLE), DIMENSION(N,M)     :: T1, T5

      REAL(DOUBLE), DIMENSION(MMAX)          :: YN, KDX
      REAL(DOUBLE), DIMENSION(NMAX*NMAX)     :: SPKSK, SHATINV
      REAL(DOUBLE), DIMENSION(MMAX*NMAX)     :: KN
      REAL(DOUBLE), DIMENSION(NMAX*MMAX)     :: KNT, KHATT

      EQUIVALENCE (KNT, KHATT), (SPKSK, SHATINV), (YN, KDX)
      ! removed (DX, KSDYMKDX), because need both at the same time for
      ! Levenberg MArquardt, (DX, KSDYMKDX)

      DATA ONE/ 1/

      ALLOCATE( KS(N,M), KSK(N,N), SPKSKINV(N,N), STAT=NAERR )

      IF (NAERR /= 0) THEN
         WRITE (6, *) 'COULD NOT ALLOCATE KS, KSK or SPKSKINV ARRAY'
         WRITE (6, *) 'ERROR NUMBER = ', NAERR
         STOP 'OPT ALLOCATION'
      ENDIF

      ALLOCATE( KHAT(M,N), G_LM(N,M), STAT=NAERR )
      IF (NAERR /= 0) THEN
         WRITE (6, *) 'COULD NOT ALLOCATE KHAT or G_LM ARRAY'
         WRITE (6, *) 'ERROR NUMBER = ', NAERR
         STOP 'OPT ALLOCATION'
      ENDIF

      YN = 0.0D0
      DY = 0.0D0
      XN = 0.0D0
      DX = 0.0D0
      KN = 0.0D0
      DELX = 0.0D0
      XNP1 = 0.0D0
      XN_OLD = 0.0D0
      FILOPEN   = .FALSE.
      DEBUG     = .FALSE.
      DIVWARN   = .FALSE.
      CONVERGE  = .FALSE.

      G_LM = 0.0
      IDNN = 0.0
      DO I=1,N
         IDNN(I,I) = 1.0
      ENDDO

! --- LEVENBERG MARQUARDT METHOD MP
      IF( ILEVENBERG .NE. 0 )THEN
         LM = .TRUE.
         GAMMA = GAMMA_START
         RED_GAMMA = GAMMA_DEC
         INC_GAMMA = GAMMA_INC
      ELSE
         LM = .FALSE.
         GAMMA = 0.0
         RED_GAMMA = 0.0
         INC_GAMMA = 0.0
      ENDIF

      IF( DEBUG )WRITE (*, *) ' m,n:', M, N
      KFLG = .TRUE.

      RMSDELY = 0.D0
      DO I = 1, M
         DELY(I) = SQRT(SE(I))
         RMSDELY = RMSDELY + SE(I)
      END DO

      RMSDELY = SQRT(RMSDELY/M)
      IF (DEBUG) WRITE (0, *) 'rmsdely= ', RMSDELY

! --- XN WORKING VERSION OF XA (PARM)
      XN(:N) = XA(:N)
      CALL INVRT( SA, SAINV, N )

! --- SUBSTITUTE SAINV VALUES FOR ANY PROFILE RETRIEVE GAS FROM FILE "SA.INPUT"
! --- LOOP OVER RETRIEVAL GASES BUT CALL ONCE IF NEEDED
      DO KK = 1, NRET
! --- PICK OUT GASES WITH SET FLAG (IFOFF=5)
         IF (( IFPRF(KK) ) .AND. ( IFOFF(KK) == 5 )) THEN
            CALL GETSAINV( ISMIX )
            EXIT
         ENDIF
      ENDDO

      SEINV(:M) = 1.D0/SE(:M)
      ITER = 0
      RMSDY = 0.D0

! --- FOR FORWARD MODEL CALC. ONLY:
      IF(( .NOT. RETFLG ) .OR. MAXITER .EQ. 0 )THEN
         IF( .NOT. F_WRTK )THEN
            KFLG = .FALSE.
         ENDIF
         XHAT(:N) = XA(:N)
         SHAT(:N, :N) = SA(:N, :N)
         CALL FM (XHAT, YHAT, KHAT, M, N, KFLG, 0, Tflg )
         WRITE (0, *) 'OPT RETURNING EARLY...'
         IF ( ALL_SPEC_OUT ) THEN
            WRITE(91,*) Y(:M)
            WRITE(91,*) YHAT(:M)
            WRITE(91,*)
         END IF
         IF (F_WRTK) THEN
            CALL FILEOPEN( 66, 2 )
            WRITE (66, *) 'K by columns m x n'
            WRITE (66, *) M, N, ISMIX, NLEV

            DO I = 1, N
               WRITE (66, '(2000ES26.18)') (KHAT(J,I),J=1,M)
            END DO

            WRITE (66, *) 'SeInv (diagonal)'
            WRITE (66, '(2000ES26.18)') (SEINV(I),I=1,M)

            WRITE (66, *) 'SaInv n x n (block diagonal)'
            DO I = 1, N
               WRITE (66, '(2000ES26.18)') (SAINV(I,J), J=1, N)
            END DO

            WRITE (66, *) 'Sa n x n (block diagonal)'
            DO I=1,N
               WRITE(66, '(2000ES26.18)') (SA(I,J), J=1, N)
            END DO

            CALL TRNSPS (KHAT, KHATT, N, M)
            CALL MULTDIAG (KHATT, SEINV, KS(:N,:M), N, M)
            CALL MULT (KS(:N,:M), KHAT, KSK, N, M, N)
            SGN = 1
            CALL ADDSUB( SAINV(:N,:N), KSK, SHATINV, N, N, SGN)
            CALL INVRT( SHATINV, SHAT(:N,:N), N)

            WRITE (66, *) 'Shat n x n'
            DO I=1,N
               WRITE(66, '(2000ES26.18)') (SHAT(I,J), J=1, N)
            END DO

            CALL FILECLOSE( 66, 1 )

         ENDIF
         RETURN
      ENDIF

      IF( OUTPUT_LM .AND. LM )THEN
         CALL FILEOPEN( 70, 2 )
         WRITE( 70,* ) "GAMMA                 = ", GAMMA
         WRITE( 70,* ) "RED GAMMA             = ", RED_GAMMA
         WRITE( 70,* ) "INC GAMMA             = ", INC_GAMMA
         WRITE(70, '(7(A10,1X))') 'CHI_2_X', 'CHI_2_Y', 'CHI_2', 'D_CHI_2', &
                    'CHI_2_LIN', 'DIFF_LIN_EXACT', 'RATIO_D_CHI'
      END IF

      CHI_2_OLD = HUGE(CHI_2_OLD)
      SEINV_OLD = SEINV
!
! ITERATION LOOP
!
   10 CONTINUE
      ITER = ITER + 1
      IF (DEBUG) WRITE (*, *) ' iteration:', ITER

! --- CALC YN, KN
      CALL FM (XN, YN, KN, M, N, KFLG, ITER, TFLG )
      !print *,'rtn fm'
      IF ( ALL_SPEC_OUT ) THEN
         WRITE(91,*) Y(:M)
         WRITE(91,*) YN(:M)
         WRITE(91,*)
      END IF

      IF (TFLG) RETURN

      IF (DEBUG) THEN
         WRITE (0, *) 'b- ret from FM'
         WRITE (0, '(6(ES12.4))') (YN(I),I=1,M)
      ENDIF

! --- CALC X(N+1)
      SGN = -1
! --- NEEDED SOMEWHERE ELSE, CALCULATED HERE TO SAVE SOME TIME
      CALL ADDSUB( Y, YN, DY, M, ONE, SGN )
      CALL TRNSPS( KN, KNT, N, M )
      CALL ADDSUB( XN, XN_OLD, DX_OLD, N, ONE, SGN )
      CALL ADDSUB( XA, XN, DX, N, ONE, SGN )

      CALL MULT( SAINV(:N,:N), DX, SAINVDX, N, N, ONE )
      !CALL MULT( DX, SAINVDX, CHI_2_X, ONE, N, ONE )
      CHI_2_X = DOT_PRODUCT( DX(:N), SAINVDX(:N) )
      ! The expected value of the cost function is m at the minimum
      ! (Rodgers, 2000, p89)
      CHI_2_X = CHI_2_X / M

      SEINVDY(:M) = SEINV(:M)*DY(:M)
      !CALL MULT( DY, SEINVDY, CHI_2_Y, ONE, M, ONE)
      CHI_2_Y = DOT_PRODUCT( DY(:M), SEINVDY(:M) )
      CHI_2_Y = CHI_2_Y / M

      SEINVDY_OLD_SE(:M) = SEINV_OLD(:M)*DY(:M)
      !CALL MULT( DY, SEINVDY, CHI_2_Y, ONE, M, ONE)
      CHI_2_Y_OLD_SE = DOT_PRODUCT( DY(:M), SEINVDY_OLD_SE(:M) )
      CHI_2_Y_OLD_SE = CHI_2_Y_OLD_SE / M

      CHI_2 = CHI_2_X + CHI_2_Y
      CHI_2_OLD_SE = CHI_2_X + CHI_2_Y_OLD_SE

      WRITE(*,306)' GAMMA = ', GAMMA
      WRITE(*, 306) ' CHI_2_X = ' , CHI_2_X

      IF (ITER.GT.1)THEN
         WRITE(*, 306) ' CHI_2_Y = ' , CHI_2_Y , CHI_2_Y_OLD_SE
         WRITE(*, 306) ' CHI2 = '    , CHI_2, CHI_2_OLD_SE
         WRITE(*, 306) ' CHI2_OLD = ', CHI_2_OLD
         D_CHI_2 = CHI_2_OLD - CHI_2
         D_CHI_2_OLD_SE = CHI_2_OLD - CHI_2_OLD_SE
         WRITE(*, 306) ' DELTA CHI_2 = ', D_CHI_2, D_CHI_2_OLD_SE
      ELSE
         WRITE(*, 306) ' CHI_2_Y = ' , CHI_2_Y
      END IF


      WRITE(16, '(A)') 'COST FUNCTION'
      WRITE(16, '(4(A,1X))') 'CHI_2_X', 'CHI_2_Y', 'CHI_2', 'CHI_2_OLD-CHI_2'
      WRITE(16, '(4ES11.3)') CHI_2_X, CHI_2_Y, CHI_2, D_CHI_2

      IF(ITER.gt.1 .and. (CONVERGENCE .GT. 0.0) &
           & .AND.   ( D_CHI_2_OLD_SE .LT. CONVERGENCE ) &
           & .AND.  (( CHI_2_OLD      .GT. CHI_2_OLD_SE ) &
           & .OR. ABS(D_CHI_2_OLD_SE) .LT. 1.0E-5 &
           & .OR. ABS(D_CHI_2)        .LT. 1.0E-5)) THEN
         ! CONVERGED
         ! CONVERGE = .TRUE.
         GOTO 20
      END IF

      IF( LM .AND. ( ITER .GT. 1 ))THEN

         IF( OUTPUT_LM )THEN
            WRITE(70, '(4ES11.3)') CHI_2_X, CHI_2_Y, CHI_2, D_CHI_2
         END IF


         IF( CHI_2_OLD .GT. CHI_2_OLD_SE )THEN
            GAMMA = GAMMA / RED_GAMMA
            SEINV_OLD(:M) = SEINV(:M)
         ELSE
            GAMMA    = GAMMA * INC_GAMMA
            XN(:N)   = XN_OLD(:N)
            YN(:M)   = YN_OLD(:M)
            KN(:M*N) = KN_OLD(:M*N)
            CHI_2    = CHI_2_OLD
            SEINV(:M) = SEINV_OLD(:M)
            ! CHI_2_LIN = CHI_2_LIN_OLD;
! --- THOSE HAVE TO BE RECALULATED, CAN SURELY BE SHORTCUT A LITTLE BIT.
            SGN = -1
            CALL ADDSUB (Y, YN, DY, M, ONE, SGN)
! --- HAS BEEN SOMEWHERE ELSE, JUST SHIFTED, NEED KNT LATER
            CALL TRNSPS (KN, KNT, N, M)
            CALL ADDSUB (XA, XN, DX, N, ONE, SGN)
         END IF
      ELSE
         SGN = -1
         SEINV_OLD(:M) = SEINV(:M)
      END IF

! --- END CALCULATION COST FUNCTION
! --- KEEP OLD STATE INFORMATION -- MP
      XN_OLD(:N)   = XN(:N)
      YN_OLD(:M)   = YN(:M)
      KN_OLD(:M*N) = KN(:M*N)
      CHI_2_OLD    = CHI_2;
     ! CHI_2_LIN_OLD = CHI_2_LIN;

      RMSIM1 = RMSDY
      RMSDY = DOT_PRODUCT(DY(:M),DY(:M))
      RMSDY = SQRT(RMSDY/M)
      IF(( ITER .GT. 1 ) .AND.(( RMSDY-RMSIM1 ) > 2*RMSDELY ))THEN
         WRITE (*, *) RMSDY, RMSIM1, RMSDELY
         WRITE (*, *) 'Threat of divergence after', ITER, ' iterations'
         DIVWARN = .TRUE.
      ENDIF

      IF( DEBUG )THEN
         WRITE (0, *) 'xa:'
         WRITE (0, '(6(ES12.4))') XA
         WRITE (0, *) 'xn:'
         WRITE (0, '(6(ES12.4))') XN
         WRITE (0, *) 'dx:'
         WRITE (0, '(6(ES12.4))') DX
         WRITE (0, *) 'yn:'
         WRITE (0, '(6(ES12.4))') YN
         WRITE (0, *) 'dy:'
         WRITE (0, '(6(ES12.4))') DY
      ENDIF

      CALL MULT (KN, DX, KDX, M, N, ONE)

      CALL ADDSUB (DY, KDX, DYMKDX, M, ONE, SGN)

      BND: DO IBAND = 1, NBAND
         NS = NSCAN(IBAND)
         IF (NS == 0) CYCLE
         IF( BSNR(IBAND).lt.0)THEN
            SPC: DO JSCAN = 1, NS
               IYDX1 = ISCNDX(1,IBAND,JSCAN)
               IYDX2 = ISCNDX(2,IBAND,JSCAN)
               SQRMS = (IYDX2-IYDX1) &
                    & /DOT_PRODUCT(DY(IYDX1:IYDX2),DY(IYDX1:IYDX2))
               VQRMS = SQRT(1.D0/SQRMS)
               SEINV(IYDX1:IYDX2) = SQRMS
               DELY(IYDX1:IYDX2)  = VQRMS
               WRITE(*,305) IBAND, JSCAN, SQRT( SQRMS)
            ENDDO SPC
         ENDIF
      ENDDO BND
      !RMSDELY = 1.D0/SUM(SEINV(:M))
      !RMSDELY = SQRT(RMSDELY/M)
      !WRITE(*,300)' RMS:DELY = ',RMSDELY, 'avgSNR = ',SNR,' TOL = ',TOL

!203   CONTINUE

      CALL MULTDIAG (KNT, SEINV, KS(:N,:M), N, M)
      CALL MULT (KS(:N,:M), DYMKDX, KSDYMKDX, N, M, ONE)
      CALL MULT (KS(:N,:M), KN, KSK, N, M, N)

      IF( DEBUG )THEN
!          write(*,'(6(es12.4))')'KS',KS
!          write(*,'(6(es12.4))')'KSdymKdx:',KSdymKdx
!          write(*,*)'KSK:',KSK
      ENDIF

      SGN = 1
! --- USE LEVENBERG MARQUARDT METHOD -- MP
      IF( LM )THEN
         G(:N) = GAMMA
         CALL MULTDIAG ( SAINV(:N,:N), G, GSAINV(:N,:N), N, N )
         CALL MULT ( GSAINV(:N,:N), DX, GSAINVDX, N, N, ONE )
         CALL ADDSUB ( KSDYMKDX, GSAINVDX, KSDYMKDX_LM, N, ONE, -1 )
         KSDYMKDX(:N) = KSDYMKDX_LM(:N)

         G(:N) = 1.0D0 + GAMMA
         CALL MULTDIAG ( SAINV(:N,:N), G, GSAINV(:N,:N), N, N )
         CALL ADDSUB ( GSAINV(:N,:N), KSK, SPKSK, N, N, SGN )
         CALL INVRT (SPKSK, SPKSKINV, N)

         ! Calculate T for AVK
         ! (Ceccherini & Ridolfi, ACP, 10, 3131-3139, 2009)
         CALL MULT(SPKSKINV, KS, T1,N,N,M)
         CALL MULT(SPKSKINV, KSK, T2,N,N,N)
         CALL MULT(SPKSKINV, SAINV, T3,N,N,N)
         CALL ADDSUB(IDNN,T2,T4,N,N,-1)
         CALL ADDSUB(T4,T3,T2,N,N,-1)         ! T2 is recycled
         CALL MULT(T2,G_LM,T5,N,N,M)
         CALL ADDSUB(T1,T5,G_LM,N,M,1)           ! T is updated
      ELSE
         CALL ADDSUB ( SAINV(:N,:N), KSK, SPKSK, N, N, SGN )
         CALL INVRT (SPKSK, SPKSKINV, N)
      END IF

      CALL MULT (SPKSKINV, KSDYMKDX, DELX, N, N, ONE)
      CALL ADDSUB (XA, DELX, XNP1, N, ONE, SGN)

      IF( DEBUG )THEN
!          write(*,*)'SpKSK:',SpKSK
!          write(*,*)'SpKSKInv',SpKSKInv
         WRITE (0, *) 'delx:'
         WRITE (0, '(6(ES12.4))') DELX
         WRITE (0, *) ' xnp1: '
         WRITE (0, '(6(ES12.4))') XNP1(:N)
      ENDIF

! --- CHECK CONVERGENCE
      SGN = -1
      IF( CONVERGENCE .gt. 0.0 )THEN
         CONVERGE = .FALSE.
      ELSE
         CALL ADDSUB (XNP1, XN, DX, N, ONE, SGN)
         IF( DEBUG )THEN
            WRITE (0, *) ' dx: '
            WRITE (0, '(6(ES12.4))') DX(:N)
         ENDIF
         CALL MULT (KN, DX, KDX, M, N, ONE)
         IF( DEBUG )THEN
            WRITE (0, *) 'Tol:'
            WRITE (0, '(6(ES12.4))') TOL
         ENDIF

         RMSKDX = 0.0D0
         DO I=1,M
            RMSKDX = RMSKDX + KDX(I)*KDX(I)
         ENDDO
         RMSKDX = SQRT(RMSKDX/M)
         WRITE(*,300)' RMS:Kdx = ', RMSKDX

         IF( DEBUG )THEN
            DO I = 1, M
               WRITE (0, *) 'Kdx(i),dely(i):', KDX(I), DELY(I)
               CHGY = ABS(KDX(I))
               UNCY = ABS(TOL*DELY(I))
               IF( CHGY .GT. UNCY )THEN
                  CONVERGE = .FALSE.
                  EXIT
               ENDIF
               CONVERGE = .TRUE.
            END DO
         ELSE
            DO I = 1, M
               CHGY = ABS(KDX(I))
               UNCY = ABS(TOL*DELY(I))
               IF( CHGY .GT. UNCY )THEN
                  CONVERGE = .FALSE.
                  EXIT
               ENDIF
               CONVERGE = .TRUE.
            END DO
         ENDIF
      END IF

! --- IN LEVENBERG MARQUARDT MODE, THIS FLAG IS ALWAYS FALSE,
!     CONVERGANCE IS DETERMINED SOMEWHERE ELSE.
      IF( CONVERGE )THEN
!         GO TO 20
!         SEINV(:M) = SEINV(:M) * 0.75
!         SNR = SQRT(SUM(SEINV(:M))/M)
!         PRINT *, ""
!         PRINT *, "  Reducing final SNR to : ", SNR
!         GO TO 10
      ELSE
         IF( ITER .GE. MAXITER )THEN
            IF( LM )THEN
               XNP1(:N) = XN_OLD(:N)
            END IF
            WRITE(0, 307) ITER
!            SEINV(:M) = SEINV(:M) * 0.75
            SNR = SQRT(SUM(SEINV(:M))/M)
            GO TO 20
         ENDIF
         XN(:N) = XNP1(:N)
         GOTO 10
      ENDIF

   20 CONTINUE
      CONVERGE = .TRUE.
      XHAT(:N) = XNP1(:N)
      CALL FM (XHAT, YHAT, KHAT, M, N, KFLG, -1, TFLG )
      IF( DEBUG )WRITE (0, *) 'a- ret from FM'

      IF( F_WRTK )THEN
!          write out K's - BJC 8/97
         CALL FILEOPEN( 66, 2 )
         WRITE (66, *) 'K by columns m x n'
         WRITE (66, *) M, N, ISMIX, NLEV

         DO I = 1, N
            WRITE (66, '(2000ES26.18)') (KHAT(J,I),J=1,M)
         END DO

         WRITE (66, *) 'SeInv (diagonal)'
         WRITE (66, '(2000ES26.18)') (SEINV(I),I=1,M)

         WRITE (66, *) 'SaInv n x n (block diagonal)'
         DO I = 1, N
            WRITE (66, '(2000ES26.18)') (SAINV(I,J), J=1, N)
         END DO

         WRITE (66, *) 'Sa n x n (block diagonal)'
         DO I=1,N
            WRITE(66, '(2000ES26.18)') (SA(I,J), J=1, N)
         END DO

      ENDIF

      IF( DEBUG )THEN
         WRITE (0, *) 'yhat:'
         WRITE (0, '(6(ES12.4))') YHAT(1:M)
         WRITE (0, *) 'xhat:'
         WRITE (0, '(6(ES12.4))') XHAT(1:N)
      ENDIF

      CALL TRNSPS (KHAT, KHATT, N, M)
      CALL MULTDIAG (KHATT, SEINV, KS(:N,:M), N, M)
      CALL MULT (KS(:N,:M), KHAT, KSK, N, M, N)
      SGN = 1

      IF (LM) THEN
         ! Calculate T for AVK iteratively
         ! (Ceccherini & Ridolfi, ACP, 10, 3131-3139, 2009)
         CALL MULTDIAG ( SAINV(:N,:N), G, GSAINV(:N,:N), N, N )
         CALL ADDSUB ( GSAINV(:N,:N), KSK, SPKSK, N, N, SGN )
         CALL INVRT (SPKSK, SPKSKINV, N)
         CALL MULT(SPKSKINV, KS, T1,N,N,M)
         CALL MULT(SPKSKINV, KSK, T2,N,N,N)
         CALL MULT(SPKSKINV, SAINV, T3,N,N,N)
         CALL ADDSUB(IDNN,T2,T4,N,N,-1)
         CALL ADDSUB(T4,T3,T2,N,N,-1)         ! T2 is recycled
         CALL MULT(T2,G_LM,T5,N,N,M)
         CALL ADDSUB(T1,T5,G_LM,N,M,1)           ! T is updated
      end if
      CALL ADDSUB( SAINV(:N,:N), KSK, SHATINV, N, N, SGN)
      CALL INVRT( SHATINV, SHAT(:N,:N), N)

      IF( F_WRTK )THEN
         WRITE (66, *) 'Shat n x n'
         DO I=1,N
            !J = (I-1)*N
            !WRITE (61, '(255ES14.6)') (SA(K),K=J+1,J+N)
            !WRITE(66, *) (SA(i,K),K=1,N)
            WRITE(66, '(2000ES26.18)') (SHAT(I,J), J=1, N)
         END DO
         CALL FILECLOSE( 66, 1 )
      ENDIF

      IF( OUTPUT_LM .AND. LM )THEN
         CALL FILECLOSE( 70, 1 )
      END IF

! --- DEALLOCATE ARRAYS
      DEALLOCATE( SPKSKINV )

      RETURN

 300  FORMAT( 3(A16, ES11.4 ))
! 304  FORMAT( "    BAND    SCAN      SeINV         DelY         SNR" )
 305  FORMAT( "BAND: ",I2,1x,"SCAN: ",I2,1x,"SNR:",1x,F7.1 )
 306  FORMAT( A20, 2ES12.4 )
 307  FORMAT(/, ' NO CONVERGENCE AFTER', I4, ' ITERATIONS')

      END SUBROUTINE OPT_3


SUBROUTINE GETSAINV( ISMIX )

      REAL(DOUBLE), DIMENSION(:,:), ALLOCATABLE :: SAINP
      REAL(DOUBLE), DIMENSION(:),   ALLOCATABLE :: SAROOT
      LOGICAL ::FILOPEN
      INTEGER :: ISMIX, KK, I, J, INDXX, NAERR

      FILOPEN   = .FALSE.

      INDXX = ISMIX
      DO KK = 1, NRET
! --- PICK OUT GASES WITH SET FLAG (IFOFF=5)
         IF (( IFPRF(KK) ) .AND. ( IFOFF(KK) == 5 )) THEN
            ALLOCATE( SAINP(LAYMAX,LAYMAX), SAROOT(LAYMAX), STAT=NAERR )
            IF (NAERR /= 0) THEN
               WRITE (6, *) 'COULD NOT ALLOCATE SAINP ARRAY'
               WRITE (6, *) 'ERROR NUMBER = ', NAERR
               STOP 'OPT ALLOCATION'
            ENDIF
            IF ( .NOT. FILOPEN ) THEN
               CALL FILEOPEN( 62, 3 )
               WRITE(16,301) TRIM(TFILE(62))
               FILOPEN = .TRUE.
            ENDIF
            WRITE(16,302) KK, IFOFF(KK)
! --- GET DIAGONAL FRACTIONAL VARIANCES FROM BINPUT && ALREADY SQUARED
            DO I = 1, NLAYERS
               SAROOT(I) = SQRT( SA(I+INDXX,I+INDXX) )
            END DO
! --- READ IN PUT NEW SAINV VALUES
            READ(62,*) SAINP(:NLAYERS, :NLAYERS )
! --- SCALE AND STORE THEM IN SAINV
            DO I = 1, NLAYERS
               DO J = 1, NLAYERS
                  SAINV(I+INDXX,J+INDXX) = SAINP(I,J)*( 1.0D0/(SAROOT(I)*SAROOT(J)))
               END DO
            END DO
            DEALLOCATE( SAINP, SAROOT )
         ELSE
            WRITE(16,303) KK, IFOFF(KK)
         ENDIF
! --- BUMP UP INDEX OF MIXING RATIO BLOCK IN SA(INV) MATRIX
         INDXX = INDXX + NLAYERS
      ENDDO
      IF( FILOPEN ) CALL FILECLOSE( 62, 2 )

      RETURN

 301  FORMAT(/," FILE: ", A, " OPENED IN OPT" )
 302  FORMAT("  RETRIEVAL GAS # :",I3," HAS IFOFF FLAG:", I3, &
             " ...READING IN NEW SAINV VALUES." )
 303  FORMAT("  RETRIEVAL GAS # :",I3," HAS IFOFF FLAG:", I3, " ...SKIPPING." )

END SUBROUTINE GETSAINV


      END MODULE OPT