Ticket #200: gpsro_op.f90

File gpsro_op.f90, 8.0 KB (added by Huw Lewis, 14 years ago)

forward operator subroutines

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