Ticket #461: ropp_fm_refrac_1d_new.f90

File ropp_fm_refrac_1d_new.f90, 9.0 KB (added by Ian Culverwell, 8 years ago)

ropp_fm_refrac_1d_new.f90_above_highest_lev

Line 
1! $Id: ropp_fm_refrac_1d_new.f90 4060 2014-02-03 11:13:00Z idculv $
2
3SUBROUTINE ropp_fm_refrac_1d_new(x, y)
4
5!****s* Refractivity/ropp_fm_refrac_1d_new *
6!
7! NAME
8! ropp_fm_refrac_1d_new - Forward model calculating a one dimensional
9! refractivity profile from the state vector
10! assuming hydrostatic variation of refractivity
11! between model levels.
12!
13! SYNOPSIS
14! call ropp_fm_refrac_1d_new(x, y)
15!
16! DESCRIPTION
17! This routine is a forward model calculating a vertical profile of
18! refractivity from profiles of temperature, humidity and pressure.
19! Refractivity values are calculated for the geopotential height levels
20! given in the observation vector.
21!
22! INPUTS
23! type(State1dFM) :: x ! State vector structure
24! type(Obs1dRefrac) :: y ! Observation vector (levels required)
25!
26! OUTPUT
27! type(Obs1dRefrac) :: y ! Obs vector with forward modelled refrac
28!
29! NOTES
30! The forward model assumes that the state vector structure contains
31! temperature, humidity and pressure values on common geopotential height
32! levels, independent of the source of those data. Model-dependent
33! conversion routines can be used to accomplish this within the
34! ropp_fm_roprof2state() subroutine.
35!
36! The forward model assumes that the observation vector contains
37! geopotential height levels onto which the forward simulated
38! observations are interpolated.
39!
40! The interpolation of refractivity calculated at the state vector's
41! geopotential height levels to the observation vector's geopotential height
42! levels is carried out assuming that the temperature varies linearly within
43! the layer and the pressures maintain hydrostatic balance.
44!
45! SEE ALSO
46! ropp_fm_types
47! ropp_fm_refrac_1d_new_ad
48! ropp_fm_refrac_1d_new_tl
49!
50! AUTHOR
51! Met Office, Exeter, UK.
52! Any comments on this software should be given via the ROM SAF
53! Helpdesk at http://www.romsaf.org
54!
55! COPYRIGHT
56! (c) EUMETSAT. All rights reserved.
57! For further details please refer to the file COPYRIGHT
58! which you should have received as part of this distribution.
59!
60!****
61
62!-------------------------------------------------------------------------------
63! 1. Declarations
64!-------------------------------------------------------------------------------
65
66 USE typesizes, ONLY: wp => EightByteReal
67 USE ropp_utils, ONLY: ropp_ZDTV, ropp_MDFV
68 USE ropp_fm, not_this => ropp_fm_refrac_1d_new
69 USE ropp_fm_types
70 USE ropp_fm_constants
71
72 IMPLICIT NONE
73
74 TYPE(State1dFM), INTENT(in) :: x ! State vector
75 TYPE(Obs1dRefrac), INTENT(inout) :: y ! Observation vector
76
77! Local variables
78 REAL(wp), DIMENSION(x%n_lev) :: z_geop ! Geopotential height of model levels
79 REAL(wp), DIMENSION(x%n_lev) :: zcomp_dry_inv ! Dry compressibility
80 REAL(wp), DIMENSION(x%n_lev) :: zcomp_wet_inv ! Wet compressibility
81
82 REAL(wp) :: kap1,kap2,kap3 ! Refractivity coefficients used in routine
83
84 REAL(wp) :: q_above_temp, q_below_temp ! temporary humidity values
85
86 LOGICAL :: below_lowest_lev ! .TRUE. if ob is below lowest model level
87 LOGICAL :: above_highest_lev ! .TRUE. if ob is above highest model level
88
89 INTEGER :: i,j
90
91 REAL(wp) :: beta,t_ob,alpha,q_ob,gamma ! [Temporary variables for new
92 REAL(wp) :: p_ob,ptemp_wet,ptemp_dry,pgrad ! refractivity interpolation]
93
94 REAL(wp) :: zcomp_dry_inv_ob,zcomp_wet_inv_ob ! Comp. factors on ob level
95
96integer, parameter :: i_debug=776
97!-------------------------------------------------------------------------------
98! 2. Non ideal gas options
99!-------------------------------------------------------------------------------
100print*,'********** q_min = ', q_min, ' *****************'
101 DO i=1, SIZE(x%geop)
102 write(*, '(a,i6,es15.5)')'i, x%geop = ', i, x%geop(i)
103 enddo
104! set inverse of compressibilities
105
106 zcomp_dry_inv(:) = 1.0_wp
107 zcomp_wet_inv(:) = 1.0_wp
108
109! initialise geopotential heights
110
111 z_geop(:) = x%geop(:)
112
113 IF (x%non_ideal) THEN
114
115! if non ideal gas calculation, use adjusted coefficients
116
117 kap1 = kappa1_comp
118 kap2 = kappa2_comp
119 kap3 = kappa3_comp
120
121! calculate compressibilty and adjust geopotential heights in z_geop
122
123 CALL ropp_fm_compress(x,z_geop,zcomp_dry_inv,zcomp_wet_inv)
124
125 ELSE
126
127 kap1 = kappa1
128 kap2 = kappa2
129 kap3 = kappa3
130
131 ENDIF
132
133!-------------------------------------------------------------------------------
134! 3. New non-linear refractivity interpolation option - see RSR 15
135!-------------------------------------------------------------------------------
136
137! 3.1 Initialise refractivity to ropp_MDFV
138! ----------------------------------------
139 y%refrac(:) = ropp_MDFV
140
141! 3.2 Cycle over each observation level
142! -------------------------------------
143
144 DO i=1, SIZE(y%geop)
145
146 below_lowest_lev=.FALSE.
147 above_highest_lev=.FALSE.
148
149! 3.3 Initialise compressibilities for ideal case
150! -----------------------------------------------
151
152 zcomp_dry_inv_ob = 1.0_wp
153
154 zcomp_wet_inv_ob = 1.0_wp
155
156! 3.4 Find model layer (j to j+1) in which ob lies
157! ------------------------------------------------
158if (abs(i-i_debug) <= 2) print*,'x%n_lev, i, y%geop(i) = ', x%n_lev, i, y%geop(i)
159
160 DO j=1, x%n_lev-1
161 IF (y%geop(i) < x%geop(1)) THEN
162 below_lowest_lev=.TRUE.
163 EXIT
164 END IF
165 IF (y%geop(i) > x%geop(x%n_lev)) THEN
166 above_highest_lev=.TRUE.
167 EXIT
168 END IF
169if (abs(i-i_debug) <= 2) write(*,'(a,i6,2es15.5)') 'j, x%geop(j), x%geop(j+1) = ', j, x%geop(j), x%geop(j+1)
170 IF (y%geop(i) > x%geop(j) .AND. y%geop(i) < x%geop(j+1)) THEN
171 EXIT
172 END IF
173 END DO
174
175!if (j == x%n_lev) j = x%n_lev - 1 ! So that we use the gradients etc at the top level.
176
177if (abs(i-i_debug) <= 2) print*,'i, j = ', i, j
178
179! 3.5 If ob is below lowest model level or above highest model level, cycle
180!--------------------------------------------------------------------------
181
182! IF (below_lowest_lev) THEN
183 IF (below_lowest_lev .OR. above_highest_lev) THEN
184 y%weights(i) = 0.0_wp
185 CYCLE
186 END IF
187
188! 3.6 Calculate temperature gradient
189! ----------------------------------
190
191 beta = (x%temp(j+1)-x%temp(j))/(x%geop(j+1)-x%geop(j))
192
193 t_ob = x%temp(j) + beta * (y%geop(i) - x%geop(j))
194
195!if (t_ob < ropp_ZDTV) t_ob = ropp_ZDTV
196! 3.7 Calculate T,q,P on observation level (N(z) is forced to be continuous)
197! --------------------------------------------------------------------------
198
199 q_above_temp = x%shum(j+1)
200
201 q_below_temp = x%shum(j)
202
203 IF (x%shum(j+1)<0.0_wp) q_above_temp = q_min
204! IF (x%shum(j+1)<ropp_ZDTV) q_above_temp = q_min
205
206 IF (x%shum(j)<0.0_wp) q_below_temp = q_min
207! IF (x%shum(j)<ropp_ZDTV) q_below_temp = q_min
208
209if (abs(i-i_debug) <= 2) print*,' '
210
211if (abs(i-i_debug) <= 2) write(*,'(a,2i6,4es15.5)')'i, j, x%shum(j+1), x%geop(j+1), x%shum(j), x%geop(j) = ', &
212 i, j, x%shum(j+1), x%geop(j+1), x%shum(j), x%geop(j)
213
214 alpha = LOG(q_below_temp / q_above_temp) / (x%geop(j+1) - x%geop(j))
215
216if (abs(i-i_debug) <= 2) write(*,'(a,i6,5es15.5)')'i, q_below_temp, q_above_temp, alpha, (y%geop(i) - x%geop(j)), (x%geop(j+1) - x%geop(j)) = ', &
217 i, q_below_temp, q_above_temp, alpha, (y%geop(i) - x%geop(j)), (x%geop(j+1) - x%geop(j))
218
219 q_ob = q_below_temp * EXP(-alpha * (y%geop(i) - x%geop(j)))
220
221! 3.8 Assume exponential pressure if T(j)~=T(j+1) to avoid P_ob=infinity
222! ----------------------------------------------------------------------
223
224 IF (ABS(x%temp(j+1) - x%temp(j)) < ropp_ZDTV) THEN
225 IF (MIN(x%pres(j+1),x%pres(j)) <= 0.0_wp) THEN
226 pgrad = (x%pres(j+1)-x%pres(j))/(x%geop(j+1)-x%geop(j))
227 p_ob = x%pres(j) + pgrad * (y%geop(i) - x%geop(j))
228 ELSE
229 p_ob=x%pres(j) * EXP(LOG(x%pres(j+1)/x%pres(j))* &
230 ((y%geop(i) - x%geop(j))/(x%geop(j+1)-x%geop(j))))
231 END IF
232 ELSE
233 gamma = -g_wmo/r_dry * (LOG(x%temp(j+1)/x%temp(j)))/ &
234 (LOG(x%pres(j+1)/x%pres(j)))
235 p_ob = x%pres(j) * (t_ob / x%temp(j))**(-g_wmo / (r_dry * gamma))
236 END IF
237
238! 3.9 Non-ideal gas compressibility factors
239! -----------------------------------------
240
241 IF (x%non_ideal) THEN
242 CALL ropp_fm_compress_single( &
243 t_ob, p_ob, q_ob, zcomp_dry_inv_ob, zcomp_wet_inv_ob)
244 END IF
245
246! 3.10 Calculate refractivity
247! --------------------------
248
249 ptemp_wet = p_ob * q_ob / &
250 (epsilon_water + (1.0_wp - epsilon_water) * q_ob)
251
252 ptemp_dry = p_ob - ptemp_wet
253
254 y%refrac(i) = kap1 * ptemp_dry * zcomp_dry_inv_ob/ t_ob + &
255 kap2 * ptemp_wet * zcomp_wet_inv_ob/ t_ob**2 + &
256 kap3 * ptemp_wet * zcomp_wet_inv_ob/ t_ob
257
258if (abs(i-i_debug) <= 2) write(*,'(A,I6,8ES15.5)')'i, y%refrac(i), t_ob, p_ob, q_ob, ptemp_dry, zcomp_dry_inv_ob, ptemp_wet, zcomp_wet_inv_ob', &
259 i, y%refrac(i), t_ob, p_ob, q_ob, ptemp_dry, zcomp_dry_inv_ob, ptemp_wet, zcomp_wet_inv_ob
260
261 END DO
262
263END SUBROUTINE ropp_fm_refrac_1d_new