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