Ticket #200: gpsro_opTL.f90

File gpsro_opTL.f90, 9.7 KB (added by Huw Lewis, 14 years ago)

tangent linear subroutines

Line 
1
2!
3
4! calculate the bending angle on impact params.
5
6SUBROUTINE alpha_optl(nlev, & ! no. of model levs (=38)
7 nobs, & ! no. of bending angles in profile
8 roc, & ! radius of curv.
9 undul, & ! undulation (set to 0.0)
10 lat, & ! latitude of ob. location (degrees)
11 pres, & ! pressure on mod levels (hPa)
12 pres_prime, &
13 temp, & ! temp on model levels
14 temp_prime, &
15 q, & ! specific humidity (g/kg)
16 q_prime, &
17 zg, & ! geopotential height of model levels (m)
18 a, & ! impact parameters
19 alpha, &
20 alpha_prime) ! bending angles
21
22
23
24
25
26
27IMPLICIT NONE
28
29!
30! subroutine args.
31!
32
33INTEGER, INTENT(IN) :: nlev ! no. of p levels in state vec.
34INTEGER, INTENT(IN) :: nobs
35REAL, INTENT(IN) :: roc
36REAL, INTENT(IN) :: undul
37REAL, INTENT(IN) :: lat
38REAL, INTENT(IN) :: pres(nlev),temp(nlev),q(nlev)
39REAL, INTENT(IN) :: pres_prime(nlev),temp_prime(nlev),q_prime(nlev)
40REAL, INTENT(IN) :: zg(nlev)
41REAL, INTENT(IN) :: a(nobs)
42REAL, INTENT(OUT) :: alpha(nobs),alpha_prime(nobs)
43
44!
45! local variables
46!
47INTEGER :: i
48REAL :: refrac(nlev)
49REAL :: nr(nlev),zg_prime(nlev)
50REAL :: refrac_prime(nlev),nr_prime(nlev)
51REAL :: roc2
52
53!
54! using fixed geopotential heights
55!
56
57zg_prime = 0.0
58
59!
60! calculate refractivity on model levels
61!
62
63CALL refrac_levstl(nlev, &
64 pres, &
65 pres_prime, &
66 temp, &
67 temp_prime, &
68 q, &
69 q_prime, &
70 refrac, &
71 refrac_prime)
72
73
74!
75! calculate the refractive index * radius on model levels
76!
77
78roc2 = roc + undul
79
80CALL calc_nrtl(nlev, &
81 roc2, &
82 lat, &
83 zg, &
84 zg_prime, &
85 refrac, &
86 refrac_prime, &
87 nr, &
88 nr_prime)
89
90
91
92
93!
94! calculate the bending angle for the derived profile from profile
95!
96
97
98CALL calc_alphatl(nobs, &
99 nlev, &
100 a, &
101 refrac, &
102 refrac_prime, &
103 nr, &
104 nr_prime, &
105 alpha, &
106 alpha_prime)
107
108
109RETURN
110
111END SUBROUTINE alpha_optl
112
113! calculate the refractivity on observation heights from state vector x
114
115
116SUBROUTINE refrac_levsTL(nlev, &
117 pres, &
118 pres_prime, &
119 temp, &
120 temp_prime,&
121 q, &
122 q_prime, &
123 refrac, &
124 refrac_prime)
125
126
127
128USE refrac_info, ONLY: &
129 epsilon,&
130 aval, &
131 bval, &
132 RMDI
133
134
135IMPLICIT NONE
136
137!
138! subroutine args.
139!
140
141INTEGER, INTENT(IN) :: nlev ! no. of p levels in state vec.
142REAL, INTENT(IN) :: pres(nlev) ! in hPa
143REAL, INTENT(IN) :: pres_prime(nlev) ! in hPa
144
145REAL, INTENT(IN) :: temp(nlev)
146REAL, INTENT(IN) :: temp_prime(nlev)
147REAL, INTENT(IN) :: q(nlev)
148REAL, INTENT(IN) :: q_prime(nlev)
149REAL, INTENT(OUT) :: refrac(nlev)
150REAL, INTENT(OUT) :: refrac_prime(nlev)
151
152!
153! local variables
154!
155
156INTEGER :: i
157
158REAL :: Ndry,Nhum
159REAL :: Ndry_prime,Nhum_prime
160
161
162!
163! now calculate the refractivity on pressure levels
164!
165
166
167refrac(:) = RMDI
168
169refrac_prime(:) = 0.0
170
171DO i = 1, nlev
172
173 Ndry = aval * pres(i)/ temp(i)
174
175 Ndry_prime = Ndry/pres(i)*pres_prime(i) - Ndry/temp(i) * temp_prime(i)
176
177 Nhum = 0.0
178 Nhum_prime = 0.0
179
180
181 Nhum = &
182 1.0E-3*bval* pres(i) * q(i)/(Epsilon*temp(i)**2)
183
184 Nhum_prime = &
185 Nhum/pres(i)*pres_prime(i)-2.0*Nhum/temp(i)*temp_prime(i) + Nhum/q(i)*q_prime(i)
186
187
188!
189! refractivity on ith pressure level
190!
191
192 refrac(i) = Ndry + Nhum
193
194 refrac_prime(i) = Ndry_prime + Nhum_prime
195
196
197ENDDO
198
199RETURN
200
201END SUBROUTINE refrac_levsTL
202
203! calculate the nr product
204
205SUBROUTINE calc_nrTL(nlev, &
206 roc, &
207 lat, &
208 zg, &
209 zg_prime, &
210 refrac, &
211 refrac_prime, &
212 nr, &
213 nr_prime)
214
215
216
217
218USE refrac_info, ONLY: &
219 g, &
220 RMDI
221
222
223IMPLICIT NONE
224
225!
226! subroutine args.
227!
228
229INTEGER, INTENT(IN) :: nlev ! no. of p levels in state vec.
230REAL, INTENT(IN) :: roc
231REAL, INTENT(IN) :: lat
232REAL, INTENT(IN) :: zg(nlev)
233REAL, INTENT(IN) :: zg_prime(nlev)
234REAL, INTENT(IN) :: refrac(nlev)
235REAL, INTENT(IN) :: refrac_prime(nlev)
236REAL, INTENT(OUT) :: nr(nlev)
237REAL, INTENT(OUT) :: nr_prime(nlev)
238
239!
240! local variables
241!
242
243INTEGER :: i
244REAL :: zed(nlev)
245REAL :: rad(nlev)
246REAL :: zed_prime(nlev)
247REAL :: rad_prime(nlev)
248REAL :: grat
249REAL :: radius
250REAL :: E_rad
251REAL :: g_lat
252
253!
254! calculate the radius and g values used in the geopotential/geometric height
255! conversion.
256!
257
258grat = g_lat(lat)/g
259radius = E_rad(lat)
260
261!
262! calculate the geometric heights
263!
264
265rad(:) = RMDI
266nr(:) = RMDI
267
268rad_prime(:) = 0.0
269nr_prime(:) = 0.0
270
271
272DO i=1,nlev
273
274 IF (zg(i) > 0.0 .AND. refrac(i) > 0.0 ) THEN
275
276 zed(i)=zg(i)/(grat - zg(i)/radius)
277
278 zed_prime(i) = grat/(grat-zg(i)/radius)**2 * zg_prime(i)
279
280!
281! calculate radius value
282!
283
284 rad(i) = roc + zed(i)
285
286 rad_prime(i) = zed_prime(i)
287
288!
289! calculate the radius times refractive index product.
290!
291
292 nr(i) = rad(i) * (1.0+1.0E-6*refrac(i))
293
294 nr_prime(i) = rad_prime(i) * (1.0+1.0E-6*refrac(i)) &
295 + 1.0E-6* rad(i)* refrac_prime(i)
296
297
298
299 ENDIF
300
301
302ENDDO
303
304
305END SUBROUTINE calc_nrTL
306
307! calculate the bending angle "alpha" for impact parameters a
308
309SUBROUTINE calc_alphaTL(nobs, &
310 nlev, &
311 a, &
312 refrac, &
313 refrac_prime, &
314 nr, &
315 nr_prime, &
316 alpha, &
317 alpha_prime)
318
319
320USE refrac_info, ONLY: &
321 RMDI, &
322 pi
323
324
325
326IMPLICIT NONE
327
328!
329! subroutine args.
330!
331
332INTEGER, INTENT(IN) :: nobs ! size of ob. vector
333INTEGER, INTENT(IN) :: nlev ! no. of refractivity levels
334REAL, INTENT(IN) :: a(nobs) ! impact parameter
335REAL, INTENT(IN) :: refrac(nlev) ! refractivity values on levels
336REAL, INTENT(IN) :: refrac_prime(nlev)
337REAL, INTENT(IN) :: nr(nlev) ! index * radius product
338REAL, INTENT(IN) :: nr_prime(nlev)
339REAL, INTENT(OUT) :: alpha(nobs)
340REAL, INTENT(OUT) :: alpha_prime(nobs) ! bending angle
341
342!
343! local variables
344!
345
346INTEGER :: i,n,ibot,jbot
347REAL :: kval(nlev-1)
348REAL :: kval_prime(nlev-1)
349REAL :: root_2pia
350REAL :: ref_low
351REAL :: ref_low_prime
352REAL :: nr_low
353REAL :: nr_low_prime
354REAL :: tup,tlow
355REAL :: tup_prime,tlow_prime
356REAL :: dalpha
357REAL :: dalpha_prime
358REAL :: erf
359REAL :: diff_erf
360REAL :: diff_erf_prime
361
362
363alpha(:) = RMDI
364
365alpha_prime(:) = 0.0
366
367jbot = 1
368
369DO
370
371 IF (refrac(jbot) > 0.0 .AND. nr(jbot) > 0.0) EXIT
372
373 jbot = jbot + 1
374
375ENDDO
376
377
378!
379! calculate the exponential decay rate between levels
380!
381
382
383DO i=jbot,nlev-1
384
385 kval(i) = LOG(refrac(i)/refrac(i+1)) / &
386 MAX(1.0,(nr(i+1)-nr(i)))
387
388 Kval(i) = MAX(1.0E-6,kval(i))
389 kval_prime(i) = 0.0
390
391
392 IF (kval(i) > 1.0E-6) &
393 kval_prime(i) = ((kval(i)*(nr_prime(i)-nr_prime(i+1))) + &
394 (refrac_prime(i)/refrac(i)- &
395 refrac_prime(i+1)/refrac(i+1)))/ &
396 MAX(1.0,(nr(i+1)-nr(i)))
397
398ENDDO
399
400!
401! now calculate the bending angle values
402!
403
404
405
406DO n=1,nobs
407
408
409 IF (a(n) < nr(jbot) .OR. a(n) > nr(nlev)) CYCLE
410
411 Root_2PIa = SQRT(2.0*pi*a(n))
412
413 ibot = jbot
414
415 DO
416
417 IF (a(n) < nr(ibot+1)) EXIT
418
419 ibot=ibot+1
420
421 ENDDO
422
423!
424! set bending angle value
425!
426
427 alpha(n) = 0.0
428
429 alpha_prime(n) = 0.0
430
431 DO i = ibot, nlev-1
432
433 IF ( i == ibot) THEN
434
435 ref_low = refrac(i)*EXP(-kval(i)*(a(n)-nr(i)))
436
437
438 ref_low_prime = ref_low* &
439 (refrac_prime(i)/refrac(i) - &
440 kval_prime(i)*(a(n)-nr(i)) + &
441 kval(i)*nr_prime(i))
442
443
444 nr_low = a(n)
445
446 nr_low_prime = 0.0
447
448
449
450 ELSE
451
452 ref_low = refrac(i)
453
454 ref_low_prime = refrac_prime(i)
455
456 nr_low = nr(i)
457
458 nr_low_prime = nr_prime(i)
459
460 ENDIF
461 tup = SQRT(kval(i)*(nr(i+1)-a(n)))
462
463 tup_prime = 0.5*(kval_prime(i)*(nr(i+1)-a(n)) + &
464 kval(i)*nr_prime(i+1))/tup
465
466 tlow = 0.0
467
468 tlow_prime = 0.0
469
470 IF (i > ibot) THEN
471
472 tlow = SQRT(kval(i)*(nr(i) -a(n)))
473
474 tlow_prime = 0.5*(kval_prime(i)*(nr(i)-a(n)) + &
475 kval(i)*nr_prime(i))/tlow
476
477 ENDIF
478
479
480 IF (i < nlev-1) THEN
481
482
483 diff_erf = erf(tup) - erf(tlow)
484
485 diff_erf_prime = 2.0/SQRT(pi)*(EXP(-tup**2)*tup_prime - &
486 EXP(-tlow**2)*tlow_prime )
487
488 ELSE
489
490
491 diff_erf = 1.0 - erf(tlow)
492
493 diff_erf_prime = -2.0/SQRT(pi)*EXP(-tlow**2)*tlow_prime
494
495 ENDIF
496
497
498 dalpha = &
499 + 1.0E-6 * Root_2PIa * SQRT(kval(i)) &
500 * ref_low*EXP(kval(i)*(nr_low-a(n)))*diff_erf
501
502
503 dalpha_prime = dalpha * ( &
504 ref_low_prime/MAX(1.0E-10,ref_low) + &
505 diff_erf_prime/MAX(1.0E-10,diff_erf) + &
506 kval(i)*nr_low_prime + &
507 (nr_low-a(n) + 0.5/kval(i))* kval_prime(i))
508
509
510 alpha(n) = alpha(n) + dalpha
511
512
513
514 alpha_prime(n) = alpha_prime(n) + dalpha_prime
515
516
517 ENDDO
518
519
520ENDDO
521
522
523END SUBROUTINE calc_alphaTL
524