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

   1 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   2 !   Driver for the tangent linear model
   3 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   4
   5 PROGRAM KPP_ROOT_ADJ_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 = ind_NO2, ind_2 = ind_O3
  12
  13 ! ---  Number of functional for which sensitivities are computed
  14 ! ---  Note: this value is set for sensitivities w.r.t. all initial values
  15 ! ---       it may have to be changed for other applications
  16       INTEGER NADJ
  17       PARAMETER (NADJ = 2)
  18       KPP_REAL  Y_ADJ(NVAR,NADJ)
  19       REAL(kind=dble_p)  R1(NVAR), R2(NVAR), V1, V2
  20
  21 ! ---- TIME VARIABLES ------------------
  22
  23       STEPMIN = 0.0d0
  24       STEPMAX = 0.0d0
  25
  26       CALL SRAND(89)
  27       RTOLS = 1.0d-3
  28       DO i=1,NVAR
  29         RTOL(i) = RTOLS
  30         ATOL(i) = 1.0d-2
  31         R1(i) = 10*(RAND()-0.5d0)
  32         R2(i) = 10*(RAND()-0.5d0)
  33       END DO
  34
  35       CALL Initialize()
  36 ! ---  Note: the initial values below are adjoint values at the final time
  37      Y_ADJ(1:NVAR,1) = R1(1:NVAR)
  38      Y_ADJ(1:NVAR,2) = R2(1:NVAR)
  39
  40 ! ********** T LOOP *************************
  41
  42       CALL InitSaveData()
  43
  44       WRITE(6,990) (SPC_NAMES(MONITOR(i)), i=1,NMONITOR)
  45 990   FORMAT('DOne[%] Time[h] ',20(4X,A12))
  46
  47       T = TSTART
  48
  49         CALL GetMass( C, DVAL )
  50         WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T/3600.,   &
  51                     (C(MONITOR(i))/CFACTOR, i=1,NMONITOR),  &
  52                    (DVAL(i)/CFACTOR, i=1,NMASS)
  53 991     FORMAT(F6.1,'% ',F7.2,3X,20(E10.4,2X))
  54
  55         CALL SaveData()
  56
  57         CALL Update_SUN()
  58         CALL Update_RCONST()
  59
  60         CALL INTEGRATE_ADJ( NADJ, Y_ADJ, T, TEND )
  61
  62       V1 = 0.0d0
  63       V2 = 0.0d0
  64       DO i=1,NVAR
  65         V1 = V1 + Y_ADJ(i,1)*R2(i)
  66         V2 = V2 + Y_ADJ(i,2)*R1(i)
  67       END DO
  68
  69       PRINT*,'**************************************************'
  70       WRITE(6,887) V1
  71       WRITE(6,888) V2
  72       WRITE(6,889) ABS(V1-V2)/MAX(ABS(V1),ABS(V2))
  73 887   FORMAT('u.M''*M''.v = ',E24.14 )
  74 888   FORMAT('v.M''*M''.u = ',E24.14 )
  75 889   FORMAT('RelativeErr=',E10.3 )
  76       PRINT*,'**************************************************'
  77
  78       CALL GetMass( C, DVAL )
  79       WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T/3600.,    &
  80                     (C(MONITOR(i))/CFACTOR, i=1,NMONITOR), &
  81                     (DVAL(i)/CFACTOR, i=1,NMASS)
  82
  83       CALL SaveData()
  84
  85 ! *********** END TIME LOOP ********
  86       OPEN(20, FILE='KPP_ROOT_ADJ_results.m')
  87       WRITE(6,*) '**************************************************'
  88       WRITE(6,*) ' Concentrations and Sensitivities at final time '
  89       WRITE(6,*) ' were written in the file KPP_ROOT_ADJ_results.m'
  90       WRITE(6,*) '**************************************************'
  91       DO j=1,NADJ
  92         WRITE(20,993) ( Y_ADJ(i,j), i=1,NVAR )
  93       END DO
  94  993  FORMAT(1000(E24.16,2X))
  95
  96       PRINT*,'ADJ: d[',TRIM(SPC_NAMES(ind_1)),'](tf) / d[', &
  97             TRIM(SPC_NAMES(ind_1)),'](t0)=', Y_ADJ(ind_1,1)
  98       PRINT*,'ADJ: d[',TRIM(SPC_NAMES(ind_2)),'](tf) / d[', &
  99             TRIM(SPC_NAMES(ind_2)),'](t0)=', Y_ADJ(ind_2,2)
 100       PRINT*,'ADJ: d[',TRIM(SPC_NAMES(ind_1)),'](tf) / d[', &
 101             TRIM(SPC_NAMES(ind_2)),'](t0)=', Y_ADJ(ind_1,2)
 102       PRINT*,'ADJ: d[',TRIM(SPC_NAMES(ind_2)),'](tf) / d[', &
 103             TRIM(SPC_NAMES(ind_1)),'](t0)=', Y_ADJ(ind_2,1)
 104
 105       PRINT*,'TLM: d[NO2](tf) / d[NO2](t0)=  1.714961808143527E-002'
 106       PRINT*,'TLM: d[ O3](tf) / d[ O3](t0)= -4.447774183920545E-003'
 107       PRINT*,'TLM: d[NO2](tf) / d[ O3](t0)=  0.897512294491540'
 108       PRINT*,'TLM: d[ O3](tf) / d[NO2](t0)= -5.543729901774693E-005'
 109
 110
 111       CALL CloseSaveData()
 112
 113 END PROGRAM KPP_ROOT_ADJ_Driver
 114