chem/KPP/kpp/kpp-2.1/drv/general_soa.f90

   1 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   2 !   Driver for the Second Order Adjoint (SOA) model
   3 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   4
   5 PROGRAM KPP_ROOT_SOA_Driver
   6
   7   USE KPP_ROOT_Model
   8   USE KPP_ROOT_Initialize, ONLY: Initialize
   9
  10       KPP_REAL :: T, DVAL(NSPEC)
  11       INTEGER :: i, j, ind_1 = 1, ind_2 = 2
  12       ! INTEGER :: ind_1 = ind_NO2, ind_2 = ind_O3
  13
  14 !~~~>  Number of second order adjoints
  15 !      i.e., number of vectors U_i s.t. SOA_i = Hess*U_i
  16 !      1 <= i <= NSOA
  17
  18       INTEGER, PARAMETER :: NSOA = 2
  19
  20       KPP_REAL :: Y_tlm(NVAR,NSOA)
  21       KPP_REAL :: Y_adj(NVAR)
  22       KPP_REAL :: Y_soa(NVAR,NSOA)
  23
  24       STEPMIN = 0.0d0
  25       STEPMAX = 0.0d0
  26
  27       DO i=1,NVAR
  28         RTOL(i) = 1.0d-4
  29         ATOL(i) = 1.0d-3
  30       END DO
  31
  32       CALL Initialize()
  33 !~~~>  Tangent linear variable values at the initial time
  34       Y_tlm(1:NVAR,1:NSOA) = 0.0d0
  35       Y_tlm(ind_1,1) = 1.0d0
  36       Y_tlm(ind_2,2) = 1.0d0
  37 !~~~>  Adjoint values at the final time
  38       Y_adj(1:NVAR) = 0.0d0
  39       Y_adj(ind_1)  = 1.0d0
  40 !~~~>  2nd order adjoint values at the final time
  41       Y_soa(1:NVAR,1:NSOA) = 0.0d0
  42       Y_soa(ind_1,1) = 1.0d0
  43       Y_soa(ind_2,2) = 1.0d0
  44
  45 !~~~~~~~~~~~ Time LOOP ~~~~~~~~~~~~~
  46
  47       CALL InitSaveData()
  48
  49       T = TSTART
  50
  51         CALL GetMass( C, DVAL )
  52         WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T,      &
  53                     (TRIM(SPC_NAMES(MONITOR(i))),          &
  54                       C(MONITOR(i))/CFACTOR, i=1,NMONITOR),&
  55                     (TRIM(SMASS(i)),DVAL(i)/CFACTOR, i=1,NMASS)
  56
  57         CALL SaveData()
  58
  59         CALL Update_SUN()
  60         CALL Update_RCONST()
  61
  62         CALL INTEGRATE_SOA(NSOA, VAR, Y_tlm, Y_adj, Y_soa, T, TEND)
  63
  64
  65       CALL GetMass( C, DVAL )
  66       WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T,        &
  67                     (TRIM(SPC_NAMES(MONITOR(i))),          &
  68                       C(MONITOR(i))/CFACTOR, i=1,NMONITOR),&
  69                     (TRIM(SMASS(i)),DVAL(i)/CFACTOR, i=1,NMASS)
  70
  71       CALL SaveData()
  72
  73 !~~~~~~~~~~~ END Time LOOP ~~~~~~~~~~~~~
  74
  75       WRITE(6,*) '**************************************************'
  76       WRITE(6,*) ' Results were written in the files'
  77       WRITE(6,*) ' KPP_ROOT_[TLM|ADJ|SOA].m'
  78       WRITE(6,*) '**************************************************'
  79
  80       PRINT 995,TRIM(SPC_NAMES(ind_1)),       &
  81             TRIM(SPC_NAMES(ind_1)),           &
  82             Y_tlm(ind_1,1), Y_adj(ind_1)
  83       PRINT 995,TRIM(SPC_NAMES(ind_1)),       &
  84             TRIM(SPC_NAMES(ind_2)),           &
  85             Y_tlm(ind_1,2), Y_adj(ind_2)
  86
  87       DO j=1,NSOA
  88          PRINT 997, j,(TRIM(SPC_NAMES(i)),Y_soa(i,j),i=1,NVAR)
  89       END DO
  90
  91       OPEN(53,FILE='KPP_ROOT_TLM.m')
  92       DO j=1, NSOA
  93          WRITE(53,993), (Y_tlm(i,j),i=1,NVAR)
  94       END DO
  95       CLOSE(53)
  96
  97       OPEN(54,FILE='KPP_ROOT_ADJ.m')
  98       WRITE(54,993), (Y_adj(i),i=1,NVAR)
  99       CLOSE(54)
 100
 101       OPEN(55,FILE='KPP_ROOT_SOA.m')
 102       DO j=1, NSOA
 103          WRITE(55,993), (Y_soa(i,j),i=1,NVAR)
 104       END DO
 105       CLOSE(55)
 106
 107  991  FORMAT(F6.1,'%. T=',E10.3,3X,20(A,'=',E10.4,';',1X))
 108  993  FORMAT(1000(E24.16,2X))
 109  995  FORMAT('d[',A,'](tf)/d[',A,'](t0). TLM=',E14.7,'  ADJ=',E14.7)
 110  997  FORMAT('2nd ADJ (',I3,'): ',200(A,'=',E10.3,'; '))
 111
 112       CALL CloseSaveData()
 113
 114 END PROGRAM KPP_ROOT_SOA_Driver
 115