source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/drv/general_soa_check.f90 @ 3458

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   Driver to test the Second Order Adjoint (SOA) model
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

PROGRAM KPP_ROOT_Driver

  USE KPP_ROOT_Model
  USE KPP_ROOT_Initialize, ONLY: Initialize

      KPP_REAL :: T, DVAL(NSPEC), DV(NVAR)
      KPP_REAL :: ADJ1(NVAR), ADJ2(NVAR)
      KPP_REAL :: FinalRef(NVAR), FinalRef_tlm(NVAR)
      KPP_REAL :: FinalPert(NVAR), FinalPert_tlm(NVAR)
      KPP_REAL :: Hs(NVAR,NVAR), AdjTlm(NVAR,NVAR), SOA(NVAR)
      INTEGER :: i, j, ind_1 = 1, ind_2 = 2, ind_COST
      LOGICAL :: LHessian, LAdjTlm
 
!~~~>  Number of second order adjoints
!      i.e., number of vectors U_i s.t. Sigma_i = Hess*U_i
!      1 <= i <= NSOA

      INTEGER, PARAMETER :: NSOA = 1

      KPP_REAL :: Y_tlm(NVAR,NSOA)
      KPP_REAL :: Y_adj(NVAR,1)
      KPP_REAL :: Y_soa(NVAR,NSOA)
  
      STEPMIN = 0.0d0
      STEPMAX = 0.0d0

      DO i=1,NVAR
        RTOL(i) = 1.0d-4
        ATOL(i) = 1.0d-3
      END DO
     
     
!~~~>  The cost function is 0.5*VAR(ind_COST, tF)**2
      ind_COST = ind_O3
     
     
!~~~>  Compute the full Hessian NVARxNVAR (or not)
      LHessian = .TRUE.
     
     
!~~~>  Compute the full Adjoint*Tlm NVARxNVAR (or not)
      LAdjTlm = .TRUE.
     
!~~~~~~~~~~~~~~~~~~~~~~~~
      PRINT*,'FIRST RUN, Reference'
      CALL Initialize()
      DV = 0.0d0
      DV(ind_NO2) =  0.02*VAR(ind_NO2)
      DV(ind_NO)  =  0.04*VAR(ind_NO)
!      DV(ind_PAN) = -0.2*VAR(ind_PAN)
!      DV(ind_O3)  =  0.1*VAR(ind_O3)
!      DO i=1,NVAR
!        DV(i) = rand()*VAR(i)
!      END DO   
      Y_tlm(1:NVAR,1) = DV
      Y_adj(1:NVAR,1) = 0.0d0
      Y_soa(1:NVAR,1:NSOA) = 0.0d0
      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
                            ATOL, RTOL)
      FinalRef     = VAR
      FinalRef_tlm(1:NVAR) = Y_tlm(1:NVAR,1)
!~~~~~~~~~~~~~~~~~~~~~~~~   
     
!~~~~~~~~~~~~~~~~~~~~~~~~
      PRINT*,'FIRST RUN, Perturbed'
      CALL Initialize()
      VAR = VAR + DV
      Y_tlm(1:NVAR,1) = DV
      Y_adj(1:NVAR,1) = 0.0d0
      Y_soa(1:NVAR,1:NSOA) = 0.0d0
      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
                            ATOL, RTOL)
      FinalPert     = VAR
      FinalPert_tlm(1:NVAR) = Y_tlm(1:NVAR,1)
!~~~~~~~~~~~~~~~~~~~~~~~~
    
!~~~~~~~~~~~~~~~~~~~~~~~~
      PRINT*,'PERTURBED RUN'
      CALL Initialize()
      VAR = VAR + DV
      Y_tlm(1:NVAR,1) = DV
      Y_adj(1:NVAR,1) = 0.0d0
      Y_adj(ind_COST,1) = FinalPert(ind_COST)
      Y_soa(1:NVAR,1) = 0.0d0
      Y_soa(ind_COST,1) = FinalPert_tlm(ind_COST)
      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
                            ATOL, RTOL)
      F2 = VAR(ind_COST)**2/2
      ADJ2 = Y_adj(:,1)
      PRINT*,'PERTURBED FUN = ',F2
!~~~~~~~~~~~~~~~~~~~~~~~~      
      
!~~~~~~~~~~~~~~~~~~~~~~~~
      PRINT*,'REFERENCE RUN'
      CALL Initialize()
      Y_tlm(1:NVAR,1) = DV
      Y_adj(1:NVAR,1) = 0.0d0
      Y_adj(ind_COST,1) = FinalRef(ind_COST)
      Y_soa(1:NVAR,1) = 0.0d0
      Y_soa(ind_COST,1) = FinalRef_tlm(ind_COST)
      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
                            ATOL, RTOL)
      F1 = VAR(ind_COST)**2/2
      ADJ1 = Y_adj(:,1)
      SOA  = Y_soa(:,1)
      PRINT*,'REFERENCE FUN = ',F1
!~~~~~~~~~~~~~~~~~~~~~~~~   
      
!~~~~~~~~~~~~~~~~~~~~~~~~
      IF (LHessian) THEN
      DO i = 1,NVAR
        PRINT*,'.... FULL HESSIAN column ',i,' of ',NVAR
        CALL Initialize()
        Y_tlm(1:NVAR,1) = 0.0d0
        Y_tlm(i,1)      = 1.0d0
        Y_adj(1:NVAR,1) = 0.0d0
        Y_adj(ind_COST,1) = FinalRef(ind_COST)
        Y_soa(1:NVAR,1) = 0.0d0
        Y_soa(ind_COST,1) = FinalRef_tlm(ind_COST)
        CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
                            ATOL, RTOL)
        CALL Initialize()
        Y_soa(1:NVAR,1) = 0.0d0
        Y_soa(ind_COST,1) = Y_tlm(ind_COST,1)
        Y_tlm(1:NVAR,1) = 0.0d0
        Y_tlm(i,1)      = 1.0d0
        Y_adj(1:NVAR,1) = 0.0d0
        Y_adj(ind_COST,1) = FinalRef(ind_COST)
        CALL INTEGRATE_soa(Nsoa, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
                            ATOL, RTOL)
        Hs(1:NVAR,i) = Y_soa(1:NVAR,1)
      END DO
      ! Write Hessian in file
      OPEN(150,FILE='KPP_ROOT_Hess.m')
      DO i=1,NVAR
        WRITE(150,224) (Hs(i,j),j=1,NVAR)
      END DO  
      CLOSE(150)      
      END IF
!~~~~~~~~~~~~~~~~~~~~~~~~
    
!~~~~~~~~~~~~~~~~~~~~~~~~
      IF (LAdjTlm) THEN
      DO i = 1,NVAR
        PRINT*,'.... FULL ADJ*TLM column ',i,' of ',NVAR
        CALL Initialize()
        Y_tlm(1:NVAR,1) = 0.0d0
        Y_tlm(i,1)      = 1.0d0
        Y_adj(1:NVAR,1) = 0.0d0
        Y_adj(ind_COST,1) = FinalRef(ind_COST)
        Y_soa(1:NVAR,1) = 0.0d0
        Y_soa(ind_COST,1) = FinalRef_tlm(ind_COST)
        CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
                            ATOL, RTOL)
        CALL Initialize()
        Y_adj(1:NVAR,1) = Y_tlm(1:NVAR)
        Y_soa(1:NVAR,1) = 0.0d0
        Y_tlm(1:NVAR,1) = 0.0d0
        CALL INTEGRATE_soa(Nsoa, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
                            ATOL, RTOL)
        AdjTlm(1:NVAR,i) = Y_adj(1:NVAR,1)
      END DO
      ! Write Hessian in file
      OPEN(150,FILE='KPP_ROOT_AdjTlm.m')
      DO i=1,NVAR
        WRITE(150,224) (AdjTlm(i,j),j=1,NVAR)
      END DO  
      CLOSE(150)      
      END IF
!~~~~~~~~~~~~~~~~~~~~~~~~

      WRITE(6,*) 'FD versus SOA:'      
      DO i=1,NVAR
        WRITE(6,222) ADJ2(i)-ADJ1(i), SOA(i),           &
          ABS(ADJ2(i)-ADJ1(i)-SOA(i))/                  &
          MAX(ABS(ADJ2(i)-ADJ1(i)),ABS(SOA(i)),1.d-16), &
          TRIM(SPC_NAMES(i))
      END DO  
222   FORMAT('FD=',E12.5,'   soa=',E12.5,'  Diff=',E12.5,'  (',A,')')      

      WRITE(6,FMT="(A14,1X,E14.7)") 'F2-F1 = ',F2-F1
      WRITE(6,FMT="(A14,1X,E14.7)") 'First = ',DOT_PRODUCT(ADJ1,DV)      
      WRITE(6,FMT="(A14,1X,E14.7)") 'F2-F1-First = ',F2-F1-DOT_PRODUCT(ADJ1,DV)      
      WRITE(6,FMT="(A14,1X,E14.7)") 'Second = ',DOT_PRODUCT(SOA,DV)/2      

      IF (LHessian) THEN
      WRITE(6,*) 'The Hessian:'      
      DO i=1,NVAR
        WRITE(6,223) (Hs(i,j),j=1,NVAR)
      END DO  
223   FORMAT(100(E12.5,1X))
      END IF


      OPEN(150,FILE='KPP_ROOT_Cost.m')
      WRITE(150,FMT="(A14,1X,E14.7)") 'F2-F1 = ',F2-F1
      WRITE(150,FMT="(A14,1X,E14.7)") 'First = ',DOT_PRODUCT(ADJ1,DV)      
      WRITE(150,FMT="(A14,1X,E14.7)") 'F2-F1-First = ',F2-F1-DOT_PRODUCT(ADJ1,DV)      
      WRITE(150,FMT="(A14,1X,E14.7)") 'Second = ',DOT_PRODUCT(SOA,DV)/2      
      CLOSE(150)      

      OPEN(150,FILE='KPP_ROOT_FD.m')
      DO i=1,NVAR
        WRITE(150,224) ADJ2(i)-ADJ1(i), SOA(i),        &
          ABS(ADJ2(i)-ADJ1(i)-SOA(i))/                 &
          MAX(ABS(ADJ2(i)-ADJ1(i)),ABS(SOA(i)),1.d-16)
      END DO  
      CLOSE(150)      

      OPEN(150,FILE='KPP_ROOT_Sol.m')
      DO i=1,NVAR
        WRITE(150,224) FinalRef(i)
      END DO  
      CLOSE(150)      

      OPEN(150,FILE='KPP_ROOT_Adj.m')
      DO i=1,NVAR
        WRITE(150,224) ADJ1(i)
      END DO  
      CLOSE(150)      

224   FORMAT(10000(E24.14,1X))      


END PROGRAM KPP_ROOT_Driver