[2696] | 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 | ATOL, RTOL) |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | CALL GetMass( C, DVAL ) |
---|
| 67 | WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T, & |
---|
| 68 | (TRIM(SPC_NAMES(MONITOR(i))), & |
---|
| 69 | C(MONITOR(i))/CFACTOR, i=1,NMONITOR),& |
---|
| 70 | (TRIM(SMASS(i)),DVAL(i)/CFACTOR, i=1,NMASS) |
---|
| 71 | |
---|
| 72 | CALL SaveData() |
---|
| 73 | |
---|
| 74 | !~~~~~~~~~~~ END Time LOOP ~~~~~~~~~~~~~ |
---|
| 75 | |
---|
| 76 | WRITE(6,*) '**************************************************' |
---|
| 77 | WRITE(6,*) ' Results were written in the files' |
---|
| 78 | WRITE(6,*) ' KPP_ROOT_[TLM|ADJ|SOA].m' |
---|
| 79 | WRITE(6,*) '**************************************************' |
---|
| 80 | |
---|
| 81 | PRINT 995,TRIM(SPC_NAMES(ind_1)), & |
---|
| 82 | TRIM(SPC_NAMES(ind_1)), & |
---|
| 83 | Y_tlm(ind_1,1), Y_adj(ind_1) |
---|
| 84 | PRINT 995,TRIM(SPC_NAMES(ind_1)), & |
---|
| 85 | TRIM(SPC_NAMES(ind_2)), & |
---|
| 86 | Y_tlm(ind_1,2), Y_adj(ind_2) |
---|
| 87 | |
---|
| 88 | DO j=1,NSOA |
---|
| 89 | PRINT 997, j,(TRIM(SPC_NAMES(i)),Y_soa(i,j),i=1,NVAR) |
---|
| 90 | END DO |
---|
| 91 | |
---|
| 92 | OPEN(53,FILE='KPP_ROOT_TLM.m') |
---|
| 93 | DO j=1, NSOA |
---|
| 94 | WRITE(53,993) (Y_tlm(i,j),i=1,NVAR) |
---|
| 95 | END DO |
---|
| 96 | CLOSE(53) |
---|
| 97 | |
---|
| 98 | OPEN(54,FILE='KPP_ROOT_ADJ.m') |
---|
| 99 | WRITE(54,993) (Y_adj(i),i=1,NVAR) |
---|
| 100 | CLOSE(54) |
---|
| 101 | |
---|
| 102 | OPEN(55,FILE='KPP_ROOT_SOA.m') |
---|
| 103 | DO j=1, NSOA |
---|
| 104 | WRITE(55,993) (Y_soa(i,j),i=1,NVAR) |
---|
| 105 | END DO |
---|
| 106 | CLOSE(55) |
---|
| 107 | |
---|
| 108 | 991 FORMAT(F6.1,'%. T=',E10.3,3X,20(A,'=',E10.4,';',1X)) |
---|
| 109 | 993 FORMAT(1000(E24.16,2X)) |
---|
| 110 | 995 FORMAT('d[',A,'](tf)/d[',A,'](t0). TLM=',E14.7,' ADJ=',E14.7) |
---|
| 111 | 997 FORMAT('2nd ADJ (',I3,'): ',200(A,'=',E10.3,'; ')) |
---|
| 112 | |
---|
| 113 | CALL CloseSaveData() |
---|
| 114 | |
---|
| 115 | END PROGRAM KPP_ROOT_SOA_Driver |
---|
| 116 | |
---|