| 1 | SUBROUTINE comp_slant_tec_tl(nstate,nobs,r_leo,r_gps,x,x_tl,impact,slant_tec,slant_tec_tl)
|
|---|
| 2 |
|
|---|
| 3 | IMPLICIT NONE
|
|---|
| 4 |
|
|---|
| 5 | INTEGER, INTENT(IN) :: nstate,nobs
|
|---|
| 6 | REAL, INTENT(IN) :: r_leo(nobs),r_gps(nobs)
|
|---|
| 7 | REAL, INTENT(IN) :: x(nstate)
|
|---|
| 8 | REAL, INTENT(IN) :: x_tl(nstate)
|
|---|
| 9 | REAL, INTENT(IN) :: impact(nobs)
|
|---|
| 10 | REAL, INTENT(OUT) :: slant_tec(nobs)
|
|---|
| 11 | REAL, INTENT(OUT) :: slant_tec_tl(nobs)
|
|---|
| 12 |
|
|---|
| 13 | ! local variables
|
|---|
| 14 |
|
|---|
| 15 | INTEGER :: nstep_leo,nstep_gps
|
|---|
| 16 |
|
|---|
| 17 | INTEGER :: i,j,k
|
|---|
| 18 |
|
|---|
| 19 | REAL :: ne_max,rmax,H0,H_grad
|
|---|
| 20 |
|
|---|
| 21 | REAL :: ne_max_tl,rmax_tl,H0_tl,H_grad_tl
|
|---|
| 22 |
|
|---|
| 23 | REAL :: dtheta,theta_leo,theta_gps
|
|---|
| 24 |
|
|---|
| 25 | REAL :: x_leo,x_gps
|
|---|
| 26 | REAL :: theta(3),h(3)
|
|---|
| 27 | REAL :: h_tl(3)
|
|---|
| 28 |
|
|---|
| 29 | REAL :: rad(3)
|
|---|
| 30 | REAL :: u(3)
|
|---|
| 31 | REAL :: ne(3)
|
|---|
| 32 |
|
|---|
| 33 | REAL :: u_tl(3)
|
|---|
| 34 | REAL :: ne_tl(3)
|
|---|
| 35 |
|
|---|
| 36 | REAL :: afac,S_term,S_term_tl
|
|---|
| 37 |
|
|---|
| 38 | ! for integral - quick test
|
|---|
| 39 |
|
|---|
| 40 | nstep_leo = 40
|
|---|
| 41 | nstep_gps = 20
|
|---|
| 42 |
|
|---|
| 43 | ! state vector contents
|
|---|
| 44 |
|
|---|
| 45 | ne_max = x(1)
|
|---|
| 46 | rmax = x(2)
|
|---|
| 47 | H0 = x(3)
|
|---|
| 48 | H_grad = x(4)
|
|---|
| 49 |
|
|---|
| 50 | ne_max_tl = x_tl(1)
|
|---|
| 51 | rmax_tl = x_tl(2)
|
|---|
| 52 | H0_tl = x_tl(3)
|
|---|
| 53 | H_grad_tl = x_tl(4)
|
|---|
| 54 |
|
|---|
| 55 |
|
|---|
| 56 | ! Maximum for the VaryCHAP
|
|---|
| 57 |
|
|---|
| 58 | !
|
|---|
| 59 | ! Compute slant TECs for impact parameters
|
|---|
| 60 | !
|
|---|
| 61 |
|
|---|
| 62 | slant_tec(:) = -9998765.0 ! missing
|
|---|
| 63 |
|
|---|
| 64 | slant_tec_tl(:) = -9998765.0 ! missing
|
|---|
| 65 |
|
|---|
| 66 | DO i = 1, nobs
|
|---|
| 67 |
|
|---|
| 68 | IF (r_leo(i) < impact(i) .OR. impact(i) < 0.0) CYCLE
|
|---|
| 69 |
|
|---|
| 70 | x_leo = r_leo(i)/impact(i)
|
|---|
| 71 |
|
|---|
| 72 | x_gps = r_gps(i)/impact(i)
|
|---|
| 73 |
|
|---|
| 74 | theta_leo = log(x_leo+SQRT(MAX(1.0E-20,(x_leo*x_leo-1.0))))
|
|---|
| 75 |
|
|---|
| 76 | theta_gps = log(x_gps+SQRT(MAX(1.0E-20,(x_gps*x_gps-1.0))))
|
|---|
| 77 |
|
|---|
| 78 | ! integrate to LEO
|
|---|
| 79 |
|
|---|
| 80 | slant_tec(i) = 0.0
|
|---|
| 81 |
|
|---|
| 82 | slant_tec_tl(i) = 0.0
|
|---|
| 83 |
|
|---|
| 84 | DO j = 1,(nstep_leo+nstep_gps)
|
|---|
| 85 |
|
|---|
| 86 | IF (j <= nstep_leo) THEN
|
|---|
| 87 |
|
|---|
| 88 | ! path from tangent point to LEO
|
|---|
| 89 |
|
|---|
| 90 | dtheta = theta_leo/REAL(nstep_leo)
|
|---|
| 91 |
|
|---|
| 92 | afac = 0.33333333
|
|---|
| 93 |
|
|---|
| 94 | ELSE
|
|---|
| 95 |
|
|---|
| 96 | ! path from LEO to GPS
|
|---|
| 97 |
|
|---|
| 98 | dtheta = (theta_gps-theta_leo)/nstep_gps
|
|---|
| 99 |
|
|---|
| 100 | afac = 0.16666667
|
|---|
| 101 |
|
|---|
| 102 | ENDIF
|
|---|
| 103 |
|
|---|
| 104 | DO k = 1,3
|
|---|
| 105 |
|
|---|
| 106 | IF (j <= nstep_leo) THEN
|
|---|
| 107 |
|
|---|
| 108 | theta(k) = REAL(j-1)*dtheta + REAL(k-1)*0.5*dtheta
|
|---|
| 109 |
|
|---|
| 110 | ELSE
|
|---|
| 111 |
|
|---|
| 112 | theta(k) = theta_leo + REAL(j-nstep_leo-1)*dtheta + REAL(k-1)*0.5*dtheta
|
|---|
| 113 |
|
|---|
| 114 | ENDIF
|
|---|
| 115 |
|
|---|
| 116 | rad(k) = impact(i)*COSH(theta(k))
|
|---|
| 117 |
|
|---|
| 118 | IF (h_grad > 1.0e-6 .AND. rad(k) > rmax ) THEN
|
|---|
| 119 |
|
|---|
| 120 | u(k) = LOG(1.0+H_grad*(rad(k)-rmax)/H0)/H_grad
|
|---|
| 121 |
|
|---|
| 122 | u_tl(k) = -(u(k)/H_grad)*H_grad_tl + &
|
|---|
| 123 | &((1.0/(1.0+H_grad*(rad(k)-rmax)/H0))*(H_grad_tl*(rad(k)-rmax)/H0 - &
|
|---|
| 124 | &(H_grad/H0)*rmax_tl - (H_grad*(rad(k)-rmax)/H0**2)*H0_tl ))/H_grad
|
|---|
| 125 |
|
|---|
| 126 | h(k) = H0 + h_grad*(rad(k)-rmax)
|
|---|
| 127 |
|
|---|
| 128 | h_tl(k) = H0_tl + h_grad_tl*(rad(k)-rmax) - h_grad*rmax_tl
|
|---|
| 129 |
|
|---|
| 130 | ELSE
|
|---|
| 131 |
|
|---|
| 132 | u(k) = (rad(k)-rmax)/H0
|
|---|
| 133 |
|
|---|
| 134 | u_tl(k) = -rmax_tl/H0 - (u(k)/H0)*H0_tl
|
|---|
| 135 |
|
|---|
| 136 | h(k) = H0
|
|---|
| 137 |
|
|---|
| 138 | h_tl(k) = H0_tl
|
|---|
| 139 |
|
|---|
| 140 | ENDIF
|
|---|
| 141 |
|
|---|
| 142 | ne(k) = ne_max*EXP(0.5*(1.0-u(k)-EXP(-u(k))))*SQRT(H0/h(k))
|
|---|
| 143 |
|
|---|
| 144 | ne_tl(k) = ne(k)*(ne_max_tl/ne_max + 0.5*(EXP(-u(k))-1.0)*u_tl(k) + &
|
|---|
| 145 | & 0.5*(H0_tl/H0 - h_tl(k)/h(k)))
|
|---|
| 146 |
|
|---|
| 147 |
|
|---|
| 148 | ENDDO ! k
|
|---|
| 149 |
|
|---|
| 150 | !!! write (6,*) j,'ne',2.0*(ne_up-ne_low)/(ne_low+ne_up),ne_up,ne_low
|
|---|
| 151 |
|
|---|
| 152 | !!! write (6,*) j,theta(1),dtheta,rad(1)-6.371e6,rad(3)-6.371e6,ne(1),ne(3)
|
|---|
| 153 |
|
|---|
| 154 | ! Simpson's rule
|
|---|
| 155 |
|
|---|
| 156 | S_term = (ne(1)*rad(1) + 4.0*ne(2)*rad(2) + ne(3)*rad(3))
|
|---|
| 157 |
|
|---|
| 158 | S_term_tl = (ne_tl(1)*rad(1) + 4.0*ne_tl(2)*rad(2) + ne_tl(3)*rad(3))
|
|---|
| 159 | !
|
|---|
| 160 | ! afac depends on which section of path
|
|---|
| 161 | !
|
|---|
| 162 | slant_tec(i) = slant_tec(i) + afac*S_term*dtheta
|
|---|
| 163 |
|
|---|
| 164 | slant_tec_tl(i) = slant_tec_tl(i) + afac*S_term_tl*dtheta
|
|---|
| 165 |
|
|---|
| 166 | ENDDO ! j
|
|---|
| 167 |
|
|---|
| 168 |
|
|---|
| 169 | !! write (6,*) i,impact(i),slant_tec(i),slant_tec_tl(i)
|
|---|
| 170 |
|
|---|
| 171 | ENDDO ! i
|
|---|
| 172 |
|
|---|
| 173 |
|
|---|
| 174 | END
|
|---|
| 175 |
|
|---|
| 176 |
|
|---|
