| 3 |
michaesp |
1 |
! $RCSfile: utilities.f90,v $
|
|
|
2 |
! $Revision: 4.11 $ $Date: 2009/11/30 14:29:09 $
|
|
|
3 |
!+ Source module for utility routines
|
|
|
4 |
!==============================================================================
|
|
|
5 |
|
|
|
6 |
MODULE utilities
|
|
|
7 |
|
|
|
8 |
!==============================================================================
|
|
|
9 |
!
|
|
|
10 |
! Description:
|
|
|
11 |
! This module provides service utilities for the model. All routines are
|
|
|
12 |
! written in a manner that also other models can use it. That means:
|
|
|
13 |
! - no routine uses other modules, except the declarations for the
|
|
|
14 |
! KIND-type parameter; the data access is by parameter list only
|
|
|
15 |
! - no routine allocates dynamic memory; work space needed is
|
|
|
16 |
! provided via the parameter list
|
|
|
17 |
! - no derived data types are used
|
|
|
18 |
!
|
|
|
19 |
! Routines (module procedures) currently contained:
|
|
|
20 |
!
|
|
|
21 |
! - convert_month:
|
|
|
22 |
! Converts a 3-character string abbreviation of a month into the number
|
|
|
23 |
! of the month or vice versa.
|
|
|
24 |
!
|
|
|
25 |
! - dfilt4:
|
|
|
26 |
! Digital filter of length 4
|
|
|
27 |
!
|
|
|
28 |
! - dfilt8:
|
|
|
29 |
! Digital filter of length 8
|
|
|
30 |
!
|
|
|
31 |
! - dolph:
|
|
|
32 |
! Calculates the Dolph-Chebyshev window for the initialization
|
|
|
33 |
!
|
|
|
34 |
! - elapsed_time:
|
|
|
35 |
! Returns the elapsed wall-clock time in seconds since the last call.
|
|
|
36 |
! On the first call the variables are only initialized. If no system
|
|
|
37 |
! clock is present, an error-value will be returned
|
|
|
38 |
!
|
|
|
39 |
! - get_utc_date:
|
|
|
40 |
! Calculates the actual date using the date of the forecast-start and
|
|
|
41 |
! the number of timesteps performed.
|
|
|
42 |
!
|
|
|
43 |
! - horizontal_filtering
|
|
|
44 |
! horizontal filtering (at the moment especially for the pressure deviation)
|
|
|
45 |
!
|
|
|
46 |
! - phirot2phi:
|
|
|
47 |
! Converts phi from the rotated system to phi in the real
|
|
|
48 |
! geographical system.
|
|
|
49 |
!
|
|
|
50 |
! - phi2phirot:
|
|
|
51 |
! Converts phi from the real geographical system to phi
|
|
|
52 |
! in the rotated system.
|
|
|
53 |
!
|
|
|
54 |
! - rlarot2rla:
|
|
|
55 |
! Converts lambda from the rotated system to lambda in the real
|
|
|
56 |
! geographical system.
|
|
|
57 |
!
|
|
|
58 |
! - rla2rlarot:
|
|
|
59 |
! Converts lambda from the real geographical system to lambda
|
|
|
60 |
! in the rotated system.
|
|
|
61 |
!
|
|
|
62 |
! - sleve_split_oro
|
|
|
63 |
! Decomposes a given topography field in a large-scale and a small-scale
|
|
|
64 |
! part according to the definition of the SLEVE coordinate
|
|
|
65 |
!
|
|
|
66 |
! - smoother:
|
|
|
67 |
! Smoothes a 2-D field by applying digital filters
|
|
|
68 |
!
|
|
|
69 |
! - tautsp:
|
|
|
70 |
! Computes tension splines
|
|
|
71 |
!
|
|
|
72 |
! - tautsp2D:
|
|
|
73 |
! Computes tension splines for several columns
|
|
|
74 |
!
|
|
|
75 |
! - to_upper:
|
|
|
76 |
! Converts alphabetic characters from lower to upper case.
|
|
|
77 |
!
|
|
|
78 |
! - uvrot2uv:
|
|
|
79 |
! Converts the wind components u and v from the rotated system
|
|
|
80 |
! to the real geographical system.
|
|
|
81 |
!
|
|
|
82 |
! - uvrot2uv_vec:
|
|
|
83 |
! the same as above, but for a whole 2D field (in vectorized form).
|
|
|
84 |
!
|
|
|
85 |
! - uv2uvrot:
|
|
|
86 |
! Converts the wind components u and v from the real geographical
|
|
|
87 |
! system to the rotated system.
|
|
|
88 |
!
|
|
|
89 |
! - uv2uvrot_vec:
|
|
|
90 |
! the same as above, but for a whole 2D field (in vectorized form).
|
|
|
91 |
!
|
|
|
92 |
! - uv2df:
|
|
|
93 |
! Converts the wind components u and v to wind direction and speed.
|
|
|
94 |
!
|
|
|
95 |
! - uv2df_vec:
|
|
|
96 |
! the same as above, but for a whole 2D field (in vectorized form).
|
|
|
97 |
!
|
|
|
98 |
!
|
|
|
99 |
! Current Code Owner: DWD, Ulrich Schaettler
|
|
|
100 |
! phone: +49 69 8062 2739
|
|
|
101 |
! fax: +49 69 8062 3721
|
|
|
102 |
! email: ulrich.schaettler@dwd.de
|
|
|
103 |
!
|
|
|
104 |
! History:
|
|
|
105 |
! Version Date Name
|
|
|
106 |
! ---------- ---------- ----
|
|
|
107 |
! 1.1 1998/03/11 Ulrich Schaettler
|
|
|
108 |
! Initial release
|
|
|
109 |
! 1.2 1998/03/30 Ulrich Schaettler
|
|
|
110 |
! Introduction of subroutine dolph used during the initialization
|
|
|
111 |
! 1.9 1998/09/16 Guenther Doms
|
|
|
112 |
! Introduction of a smoothing routine 'smoother' which uses digital
|
|
|
113 |
! filters 'dfilt4' and 'dfilt8'.
|
|
|
114 |
! 1.10 1998/09/29 Ulrich Schaettler
|
|
|
115 |
! Routine remark eliminated and put to parallel_utilities.
|
|
|
116 |
! Routines uv2uvrot and uv2df introduced
|
|
|
117 |
! 1.16 1998/11/02 Guenther Doms
|
|
|
118 |
! Correction of filter processing in routine 'smoother'.
|
|
|
119 |
! 1.29 1999/05/11 Ulrich Schaettler
|
|
|
120 |
! Adaptations to use this module also in GME2LM
|
|
|
121 |
! 1.32 1999/08/24 Guenther Doms
|
|
|
122 |
! some _ireals declarations added.
|
|
|
123 |
! 2.8 2001/07/06 Ulrich Schaettler
|
|
|
124 |
! Added new subroutines tautsp2D, uv2uvrot_vec and uvrot2uv_vec for
|
|
|
125 |
! vectorization
|
|
|
126 |
! 2.14 2002/02/15 Ulrich Schaettler
|
|
|
127 |
! Correction and adaptations in tautsp2D (analogous to GME2LM)
|
|
|
128 |
! Added new subroutine dc_topo for the SLEVE coordinate
|
|
|
129 |
! 2.17 2002/05/08 Ulrich Schaettler
|
|
|
130 |
! Modifications for performing the filtering in irealgrib-format
|
|
|
131 |
! 2.18 2002/07/16 Guenther Doms
|
|
|
132 |
! Corrections for the rotation of the wind components from or to the
|
|
|
133 |
! geographical coordinate system.
|
|
|
134 |
! 3.3 2003/04/22 Christoph Schraff
|
|
|
135 |
! Introduction of subroutines 'convert_month' and 'to_upper' (for GPS data).
|
|
|
136 |
! 3.6 2003/12/11 Ulrich Schaettler
|
|
|
137 |
! Eliminated Subroutine istringlen (use F90 intrinsic LEN_TRIM instead)
|
|
|
138 |
! 3.13 2004/12/03 Ulrich Schaettler
|
|
|
139 |
! Eliminated dependency on data_io (put irealgrib to data_parameters)
|
|
|
140 |
! New SR horizontal_filtering (from INT2LM);
|
|
|
141 |
! Renamed SR dc_topo to sleve_split_oro
|
|
|
142 |
! 3.14 2005/01/25 Ulrich Schaettler
|
|
|
143 |
! New filter routine smooth9 for new type of Rayleigh damping (Lucio Torrisi)
|
|
|
144 |
! Changes in horizontal_filtering (Jochen Foerstner)
|
|
|
145 |
! 3.15 2005/03/03 Ulrich Schaettler
|
|
|
146 |
! Replaced FLOAT by REAL
|
|
|
147 |
! 3.16 2005/07/22 Ulrich Schaettler
|
|
|
148 |
! Bug correction in the call to intrinsic function REAL
|
|
|
149 |
! 3.18 2006/03/03 Ulrich Schaettler
|
|
|
150 |
! Introduced idouble/isingle as KIND parameters instead of ireals/irealgrib
|
|
|
151 |
! in the generic formulation of some routines (dfilt4, dfilt8, smoother)
|
|
|
152 |
! Changed get_utc_date to include also a climatological year with 360 days
|
|
|
153 |
! 3.21 2006/12/04 Burkhardt Rockel, Lucio Torrisi, Jochen Foerstner
|
|
|
154 |
! Added polgam in transformation function rla2rlarot
|
|
|
155 |
! polgam is not used as optional parameter any more
|
|
|
156 |
! Some adaptations in smooth9 for itype_spubc=2
|
|
|
157 |
! Some modifications in horizontal_filtering
|
|
|
158 |
! Function uv2df_vec introduced. (C. Schraff)
|
|
|
159 |
! V3_23 2007/03/30 Ulrich Schaettler
|
|
|
160 |
! Declared some constant variables as parameters to allow inlining on
|
|
|
161 |
! some platforms
|
|
|
162 |
! Changed computation of acthour in get_utc_date
|
|
|
163 |
! V3_24 2007/04/26 Ulrich Schaettler
|
|
|
164 |
! Bug correction in computation of acthour in get_utc_date
|
|
|
165 |
! V4_1 2007/12/04 Ulrich Schaettler
|
|
|
166 |
! Introduced parameter myid to sleve_split_oro (is called from all PEs in INT2LM)
|
|
|
167 |
! V4_4 2008/07/16 Ulrich Schaettler
|
|
|
168 |
! Adapted a debug printout in SR tautsp2D
|
|
|
169 |
! Changed NL parameter lyear_360 to itype_calendar, to have several options
|
|
|
170 |
! Vectorized SR horizontal_filtering by changing some loops
|
|
|
171 |
! Treatment of very small values for spline interpolation in tautsp2D
|
|
|
172 |
! V4_8 2009/02/16 Ulrich Schaettler (Andy Dobler)
|
|
|
173 |
! Corrected leap year calculation for centuries in Gregorian calendar
|
|
|
174 |
! @VERSION@ @DATE@ <Your name>
|
|
|
175 |
! <Modification comments>
|
|
|
176 |
!
|
|
|
177 |
! Code Description:
|
|
|
178 |
! Language: Fortran 90.
|
|
|
179 |
! Software Standards: "European Standards for Writing and
|
|
|
180 |
! Documenting Exchangeable Fortran 90 Code".
|
|
|
181 |
!==============================================================================
|
|
|
182 |
!
|
|
|
183 |
! Declarations:
|
|
|
184 |
!
|
|
|
185 |
! Modules used:
|
|
|
186 |
USE data_parameters , ONLY : &
|
|
|
187 |
ireals, & ! KIND-type parameter for real variables
|
|
|
188 |
iintegers, & ! KIND-type parameter for standard integer variables
|
|
|
189 |
irealgrib, & ! KIND-type parameter for real variables in the grib library
|
|
|
190 |
idouble, & ! KIND-type parameter for double precision real variables
|
|
|
191 |
isingle ! KIND-type parameter for single precision real variables
|
|
|
192 |
|
|
|
193 |
!==============================================================================
|
|
|
194 |
|
|
|
195 |
IMPLICIT NONE
|
|
|
196 |
|
|
|
197 |
!==============================================================================
|
|
|
198 |
|
|
|
199 |
! Interface Blocks
|
|
|
200 |
INTERFACE smoother
|
|
|
201 |
MODULE PROCEDURE &
|
|
|
202 |
smoother_double, &
|
|
|
203 |
smoother_single
|
|
|
204 |
END INTERFACE
|
|
|
205 |
|
|
|
206 |
INTERFACE dfilt4
|
|
|
207 |
MODULE PROCEDURE &
|
|
|
208 |
dfilt4_double, &
|
|
|
209 |
dfilt4_single
|
|
|
210 |
END INTERFACE
|
|
|
211 |
|
|
|
212 |
INTERFACE dfilt8
|
|
|
213 |
MODULE PROCEDURE &
|
|
|
214 |
dfilt8_double, &
|
|
|
215 |
dfilt8_single
|
|
|
216 |
END INTERFACE
|
|
|
217 |
|
|
|
218 |
!==============================================================================
|
|
|
219 |
|
|
|
220 |
CONTAINS
|
|
|
221 |
|
|
|
222 |
!==============================================================================
|
|
|
223 |
|
|
|
224 |
SUBROUTINE convert_month ( MonthString, MonthNumber, ind )
|
|
|
225 |
|
|
|
226 |
!-------------------------------------------------------------------------------
|
|
|
227 |
!
|
|
|
228 |
! Description:
|
|
|
229 |
! Convert 3-chr Month string to number (ind > 0) or vice-versa (ind <= 0).
|
|
|
230 |
! If ind > 0 and input string not valid, MonthNumber will be 0.
|
|
|
231 |
! If ind <=0 and input month number invalid, MonthString will be 'XXX'.
|
|
|
232 |
! Method:
|
|
|
233 |
! Uses subroutine 'to_upper'.
|
|
|
234 |
!-------------------------------------------------------------------------------
|
|
|
235 |
|
|
|
236 |
IMPLICIT NONE
|
|
|
237 |
|
|
|
238 |
! Subroutine arguments:
|
|
|
239 |
! --------------------
|
|
|
240 |
CHARACTER (LEN=3) , INTENT(INOUT) :: MonthString ! 3-chr Month name
|
|
|
241 |
INTEGER (KIND=iintegers), INTENT(INOUT) :: MonthNumber ! Month number (1-12)
|
|
|
242 |
INTEGER (KIND=iintegers), INTENT(IN) :: ind ! > 0 : chr --> no.
|
|
|
243 |
! <=0 : no --> chr
|
|
|
244 |
|
|
|
245 |
! Local parameters:
|
|
|
246 |
! ----------------
|
|
|
247 |
CHARACTER (LEN=36), PARAMETER :: &
|
|
|
248 |
MonthNames = "JANFEBMARAPRMAYJUNJULAUGSEPOCTNOVDEC"
|
|
|
249 |
|
|
|
250 |
! Local variables
|
|
|
251 |
! ---------------
|
|
|
252 |
CHARACTER (LEN=3) :: Month
|
|
|
253 |
INTEGER (KIND=iintegers) :: idx
|
|
|
254 |
!
|
|
|
255 |
!------------ End of header ----------------------------------------------------
|
|
|
256 |
|
|
|
257 |
IF ( ind > 0 ) THEN
|
|
|
258 |
|
|
|
259 |
! ----- String to number -----
|
|
|
260 |
|
|
|
261 |
Month = MonthString
|
|
|
262 |
CALL to_upper ( Month )
|
|
|
263 |
idx = INDEX ( MonthNames, Month )
|
|
|
264 |
IF ( MOD ( idx-1, 3 ) == 0 ) THEN
|
|
|
265 |
MonthNumber = ( idx + 2 ) / 3
|
|
|
266 |
ELSE
|
|
|
267 |
MonthNumber = 0
|
|
|
268 |
END IF
|
|
|
269 |
|
|
|
270 |
ELSE
|
|
|
271 |
|
|
|
272 |
! ----- Number to string -----
|
|
|
273 |
|
|
|
274 |
IF ( MonthNumber >= 1 .AND. &
|
|
|
275 |
MonthNumber <= 12 ) THEN
|
|
|
276 |
idx = MonthNumber * 3 - 2
|
|
|
277 |
MonthString = MonthNames(idx:idx+2)
|
|
|
278 |
ELSE
|
|
|
279 |
MonthString = "XXX"
|
|
|
280 |
END IF
|
|
|
281 |
|
|
|
282 |
END IF
|
|
|
283 |
|
|
|
284 |
END SUBROUTINE convert_month
|
|
|
285 |
|
|
|
286 |
!------------------------------------------------------------------------------
|
|
|
287 |
|
|
|
288 |
!==============================================================================
|
|
|
289 |
!==============================================================================
|
|
|
290 |
!+ Defines all subroutines for the generic routine dfilt4
|
|
|
291 |
!------------------------------------------------------------------------------
|
|
|
292 |
!
|
|
|
293 |
! SUBROUTINE dfilt4 (fin, idim, fhelp, fout, nfilt)
|
|
|
294 |
!
|
|
|
295 |
!------------------------------------------------------------------------------
|
|
|
296 |
!
|
|
|
297 |
! Description:
|
|
|
298 |
! This routine smoothes an arbitrary field (fin) of length idim by applying
|
|
|
299 |
! a digital filters of length nlength 4 nfilt times. The filterd field
|
|
|
300 |
! is written on fout.
|
|
|
301 |
!
|
|
|
302 |
! Method:
|
|
|
303 |
! Digital filter according to Shapiro
|
|
|
304 |
!
|
|
|
305 |
!------------------------------------------------------------------------------
|
|
|
306 |
!+ Implementation for double precision
|
|
|
307 |
!------------------------------------------------------------------------------
|
|
|
308 |
|
|
|
309 |
SUBROUTINE dfilt4_double (fin, idim, fhelp, fout, nfilt)
|
|
|
310 |
|
|
|
311 |
!------------------------------------------------------------------------------
|
|
|
312 |
|
|
|
313 |
! Parameter list:
|
|
|
314 |
INTEGER (KIND=iintegers), INTENT (IN) :: &
|
|
|
315 |
idim, & ! Dimension of the field
|
|
|
316 |
nfilt ! Number of iterative filerings
|
|
|
317 |
REAL (KIND=idouble), INTENT (IN) :: &
|
|
|
318 |
fin (idim) ! input field (unfilterd)
|
|
|
319 |
REAL (KIND=idouble), INTENT (OUT) :: &
|
|
|
320 |
fout (idim) ! smoothed output field (filtered)
|
|
|
321 |
REAL (KIND=idouble), INTENT (INOUT) :: &
|
|
|
322 |
fhelp(idim) ! additional storage supplied by the calling routine
|
|
|
323 |
|
|
|
324 |
! Local variables
|
|
|
325 |
INTEGER (KIND=iintegers) :: &
|
|
|
326 |
i,m, & ! loop indicees
|
|
|
327 |
nf_o2 ! nfilt/2
|
|
|
328 |
|
|
|
329 |
REAL (KIND=idouble) :: &
|
|
|
330 |
fw(5) ! filter weights
|
|
|
331 |
|
|
|
332 |
!------------------------------------------------------------------------------
|
|
|
333 |
DATA fw / -0.00390625_idouble, 0.03125_idouble, -0.109375_idouble, &
|
|
|
334 |
0.21875_idouble, 0.7265625_idouble /
|
|
|
335 |
|
|
|
336 |
|
|
|
337 |
! begin subroutine dfilt4_double
|
|
|
338 |
|
|
|
339 |
nf_o2 = (nfilt+1)/2
|
|
|
340 |
|
|
|
341 |
fout (:) = fin(:)
|
|
|
342 |
fhelp(:) = fin(:)
|
|
|
343 |
|
|
|
344 |
DO i = 2, idim-1
|
|
|
345 |
fhelp(i) = 0.15_idouble*fout (i-1) + 0.7_idouble*fout (i) &
|
|
|
346 |
+ 0.15_idouble*fout (i+1)
|
|
|
347 |
ENDDO
|
|
|
348 |
DO i = 2, idim-1
|
|
|
349 |
fout (i) = 0.15_idouble*fhelp(i-1) + 0.7_idouble*fhelp(i) &
|
|
|
350 |
+ 0.15_idouble*fhelp(i+1)
|
|
|
351 |
ENDDO
|
|
|
352 |
|
|
|
353 |
DO m = 1, nf_o2
|
|
|
354 |
DO i = 5, idim-4
|
|
|
355 |
fhelp(i) = fw(5)*fout(i) &
|
|
|
356 |
+ fw(4)*(fout(i-1)+fout(i+1)) + fw(3)*(fout(i-2)+fout(i+2)) &
|
|
|
357 |
+ fw(2)*(fout(i-3)+fout(i+3)) + fw(1)*(fout(i-4)+fout(i+4))
|
|
|
358 |
ENDDO
|
|
|
359 |
DO i = 5, idim-4
|
|
|
360 |
fout(i) = fw(5)*fhelp(i) &
|
|
|
361 |
+ fw(4)*(fhelp(i-1)+fhelp(i+1)) + fw(3)*(fhelp(i-2)+fhelp(i+2)) &
|
|
|
362 |
+ fw(2)*(fhelp(i-3)+fhelp(i+3)) + fw(1)*(fhelp(i-4)+fhelp(i+4))
|
|
|
363 |
ENDDO
|
|
|
364 |
ENDDO
|
|
|
365 |
|
|
|
366 |
DO i = 2, idim-1
|
|
|
367 |
fhelp(i) = 0.15_idouble*fout (i-1) + 0.7_idouble*fout (i) &
|
|
|
368 |
+ 0.15_idouble*fout (i+1)
|
|
|
369 |
ENDDO
|
|
|
370 |
DO i = 2, idim-1
|
|
|
371 |
fout (i) = 0.15_idouble*fhelp(i-1) + 0.7_idouble*fhelp(i) &
|
|
|
372 |
+ 0.15_idouble*fhelp(i+1)
|
|
|
373 |
ENDDO
|
|
|
374 |
|
|
|
375 |
END SUBROUTINE dfilt4_double
|
|
|
376 |
|
|
|
377 |
!------------------------------------------------------------------------------
|
|
|
378 |
!+ Implementation for single precision
|
|
|
379 |
!------------------------------------------------------------------------------
|
|
|
380 |
|
|
|
381 |
SUBROUTINE dfilt4_single (fin, idim, fhelp, fout, nfilt)
|
|
|
382 |
|
|
|
383 |
!------------------------------------------------------------------------------
|
|
|
384 |
|
|
|
385 |
! Parameter list:
|
|
|
386 |
INTEGER (KIND=iintegers), INTENT (IN) :: &
|
|
|
387 |
idim, & ! Dimension of the field
|
|
|
388 |
nfilt ! Number of iterative filerings
|
|
|
389 |
REAL (KIND=isingle), INTENT (IN) :: &
|
|
|
390 |
fin (idim) ! input field (unfilterd)
|
|
|
391 |
REAL (KIND=isingle), INTENT (OUT) :: &
|
|
|
392 |
fout (idim) ! smoothed output field (filtered)
|
|
|
393 |
REAL (KIND=isingle), INTENT (INOUT) :: &
|
|
|
394 |
fhelp(idim) ! additional storage supplied by the calling routine
|
|
|
395 |
|
|
|
396 |
! Local variables
|
|
|
397 |
INTEGER (KIND=iintegers) :: &
|
|
|
398 |
i,m, & ! loop indicees
|
|
|
399 |
nf_o2 ! nfilt/2
|
|
|
400 |
|
|
|
401 |
REAL (KIND=isingle) :: &
|
|
|
402 |
fw(5) ! filter weights
|
|
|
403 |
|
|
|
404 |
!------------------------------------------------------------------------------
|
|
|
405 |
DATA fw / -0.00390625_isingle, 0.03125_isingle, -0.109375_isingle, &
|
|
|
406 |
0.21875_isingle, 0.7265625_isingle /
|
|
|
407 |
|
|
|
408 |
|
|
|
409 |
! begin subroutine dfilt4_single
|
|
|
410 |
|
|
|
411 |
nf_o2 = (nfilt+1)/2
|
|
|
412 |
|
|
|
413 |
fout (:) = fin(:)
|
|
|
414 |
fhelp(:) = fin(:)
|
|
|
415 |
|
|
|
416 |
DO i = 2, idim-1
|
|
|
417 |
fhelp(i) = 0.15_isingle*fout (i-1) + 0.7_isingle*fout (i) &
|
|
|
418 |
+ 0.15_isingle*fout (i+1)
|
|
|
419 |
ENDDO
|
|
|
420 |
DO i = 2, idim-1
|
|
|
421 |
fout (i) = 0.15_isingle*fhelp(i-1) + 0.7_isingle*fhelp(i) &
|
|
|
422 |
+ 0.15_isingle*fhelp(i+1)
|
|
|
423 |
ENDDO
|
|
|
424 |
|
|
|
425 |
DO m = 1, nf_o2
|
|
|
426 |
DO i = 5, idim-4
|
|
|
427 |
fhelp(i) = fw(5)*fout(i) &
|
|
|
428 |
+ fw(4)*(fout(i-1)+fout(i+1)) + fw(3)*(fout(i-2)+fout(i+2)) &
|
|
|
429 |
+ fw(2)*(fout(i-3)+fout(i+3)) + fw(1)*(fout(i-4)+fout(i+4))
|
|
|
430 |
ENDDO
|
|
|
431 |
DO i = 5, idim-4
|
|
|
432 |
fout(i) = fw(5)*fhelp(i) &
|
|
|
433 |
+ fw(4)*(fhelp(i-1)+fhelp(i+1)) + fw(3)*(fhelp(i-2)+fhelp(i+2)) &
|
|
|
434 |
+ fw(2)*(fhelp(i-3)+fhelp(i+3)) + fw(1)*(fhelp(i-4)+fhelp(i+4))
|
|
|
435 |
ENDDO
|
|
|
436 |
ENDDO
|
|
|
437 |
|
|
|
438 |
DO i = 2, idim-1
|
|
|
439 |
fhelp(i) = 0.15_isingle*fout (i-1) + 0.7_isingle*fout (i) &
|
|
|
440 |
+ 0.15_isingle*fout (i+1)
|
|
|
441 |
ENDDO
|
|
|
442 |
DO i = 2, idim-1
|
|
|
443 |
fout (i) = 0.15_isingle*fhelp(i-1) + 0.7_isingle*fhelp(i) &
|
|
|
444 |
+ 0.15_isingle*fhelp(i+1)
|
|
|
445 |
ENDDO
|
|
|
446 |
|
|
|
447 |
END SUBROUTINE dfilt4_single
|
|
|
448 |
|
|
|
449 |
!------------------------------------------------------------------------------
|
|
|
450 |
|
|
|
451 |
!==============================================================================
|
|
|
452 |
!==============================================================================
|
|
|
453 |
!+ Defines all subroutines for the generic routine dfilt8
|
|
|
454 |
!------------------------------------------------------------------------------
|
|
|
455 |
!
|
|
|
456 |
! SUBROUTINE dfilt8 (fin, idim, fhelp, fout, nfilt)
|
|
|
457 |
!
|
|
|
458 |
!------------------------------------------------------------------------------
|
|
|
459 |
!
|
|
|
460 |
! Description:
|
|
|
461 |
! This routine smoothes an arbitrary field (fin) of length idim by applying
|
|
|
462 |
! a digital filters of length nlength 8 nfilt times. The filterd field
|
|
|
463 |
! is written on fout.
|
|
|
464 |
!
|
|
|
465 |
! Method:
|
|
|
466 |
! Digital filter according to Shapiro
|
|
|
467 |
!
|
|
|
468 |
!------------------------------------------------------------------------------
|
|
|
469 |
!+ Implementation for double precision
|
|
|
470 |
!------------------------------------------------------------------------------
|
|
|
471 |
|
|
|
472 |
SUBROUTINE dfilt8_double (fin, idim, fhelp, fout, nfilt)
|
|
|
473 |
|
|
|
474 |
!------------------------------------------------------------------------------
|
|
|
475 |
|
|
|
476 |
! Parameter list:
|
|
|
477 |
INTEGER (KIND=iintegers), INTENT (IN) :: &
|
|
|
478 |
idim, & ! Dimension of the field
|
|
|
479 |
nfilt ! Number of iterative filerings
|
|
|
480 |
REAL (KIND=idouble), INTENT (IN) :: &
|
|
|
481 |
fin (idim) ! input field (unfilterd)
|
|
|
482 |
REAL (KIND=idouble), INTENT (OUT) :: &
|
|
|
483 |
fout (idim) ! smoothed output field (filtered)
|
|
|
484 |
REAL (KIND=idouble), INTENT (INOUT) :: &
|
|
|
485 |
fhelp(idim) ! additional storage supplied by the calling routine
|
|
|
486 |
|
|
|
487 |
! Local variables
|
|
|
488 |
INTEGER (KIND=iintegers) :: &
|
|
|
489 |
i,m, & ! loop indicees
|
|
|
490 |
nf_o2 ! nfilt/2
|
|
|
491 |
|
|
|
492 |
REAL (KIND=idouble) :: &
|
|
|
493 |
fw(9) ! filter weights
|
|
|
494 |
|
|
|
495 |
!------------------------------------------------------------------------------
|
|
|
496 |
DATA fw /-0.0000152590_idouble, 0.0002441406_idouble, -0.0018310546_idouble, &
|
|
|
497 |
0.0085449218_idouble, -0.0277709960_idouble, 0.0666503906_idouble, &
|
|
|
498 |
-0.1221923828_idouble, 0.1745605469_idouble, 0.8036193848_idouble /
|
|
|
499 |
|
|
|
500 |
! begin subroutine dfilt8_double
|
|
|
501 |
|
|
|
502 |
nf_o2 = (nfilt+1)/2
|
|
|
503 |
|
|
|
504 |
fout (:) = fin(:)
|
|
|
505 |
fhelp(:) = fin(:)
|
|
|
506 |
|
|
|
507 |
DO i = 2, idim-1
|
|
|
508 |
fhelp(i) = 0.25_idouble*fout (i-1) + 0.5_idouble*fout (i) &
|
|
|
509 |
+ 0.25_idouble*fout (i+1)
|
|
|
510 |
ENDDO
|
|
|
511 |
DO i = 2, idim-1
|
|
|
512 |
fout (i) = 0.25_idouble*fhelp(i-1) + 0.5_idouble*fhelp(i) &
|
|
|
513 |
+ 0.25_idouble*fhelp(i+1)
|
|
|
514 |
ENDDO
|
|
|
515 |
|
|
|
516 |
DO m = 1, nf_o2
|
|
|
517 |
DO i = 9, idim-8
|
|
|
518 |
fhelp(i) = fw(9)*fout(i) &
|
|
|
519 |
+ fw(8)*(fout(i-1)+fout(i+1)) + fw(7)*(fout(i-2)+fout(i+2)) &
|
|
|
520 |
+ fw(6)*(fout(i-3)+fout(i+3)) + fw(5)*(fout(i-4)+fout(i+4)) &
|
|
|
521 |
+ fw(4)*(fout(i-5)+fout(i+5)) + fw(3)*(fout(i-6)+fout(i+6)) &
|
|
|
522 |
+ fw(2)*(fout(i-7)+fout(i+7)) + fw(1)*(fout(i-8)+fout(i+8))
|
|
|
523 |
ENDDO
|
|
|
524 |
DO i = 9, idim-8
|
|
|
525 |
fout(i) = fw(9)*fhelp(i) &
|
|
|
526 |
+ fw(8)*(fhelp(i-1)+fhelp(i+1)) + fw(7)*(fhelp(i-2)+fhelp(i+2)) &
|
|
|
527 |
+ fw(6)*(fhelp(i-3)+fhelp(i+3)) + fw(5)*(fhelp(i-4)+fhelp(i+4)) &
|
|
|
528 |
+ fw(4)*(fhelp(i-5)+fhelp(i+5)) + fw(3)*(fhelp(i-6)+fhelp(i+6)) &
|
|
|
529 |
+ fw(2)*(fhelp(i-7)+fhelp(i+7)) + fw(1)*(fhelp(i-8)+fhelp(i+8))
|
|
|
530 |
ENDDO
|
|
|
531 |
ENDDO
|
|
|
532 |
|
|
|
533 |
DO i = 2, idim-1
|
|
|
534 |
fhelp(i) = 0.25_idouble*fout (i-1) + 0.5_idouble*fout (i) &
|
|
|
535 |
+ 0.25_idouble*fout (i+1)
|
|
|
536 |
ENDDO
|
|
|
537 |
DO i = 2, idim-1
|
|
|
538 |
fout (i) = 0.25_idouble*fhelp(i-1) + 0.5_idouble*fhelp(i) &
|
|
|
539 |
+ 0.25_idouble*fhelp(i+1)
|
|
|
540 |
ENDDO
|
|
|
541 |
|
|
|
542 |
END SUBROUTINE dfilt8_double
|
|
|
543 |
|
|
|
544 |
!------------------------------------------------------------------------------
|
|
|
545 |
!+ Implementation for single precision
|
|
|
546 |
!------------------------------------------------------------------------------
|
|
|
547 |
|
|
|
548 |
SUBROUTINE dfilt8_single (fin, idim, fhelp, fout, nfilt)
|
|
|
549 |
|
|
|
550 |
!------------------------------------------------------------------------------
|
|
|
551 |
|
|
|
552 |
! Parameter list:
|
|
|
553 |
INTEGER (KIND=iintegers), INTENT (IN) :: &
|
|
|
554 |
idim, & ! Dimension of the field
|
|
|
555 |
nfilt ! Number of iterative filerings
|
|
|
556 |
REAL (KIND=isingle), INTENT (IN) :: &
|
|
|
557 |
fin (idim) ! input field (unfilterd)
|
|
|
558 |
REAL (KIND=isingle), INTENT (OUT) :: &
|
|
|
559 |
fout (idim) ! smoothed output field (filtered)
|
|
|
560 |
REAL (KIND=isingle), INTENT (INOUT) :: &
|
|
|
561 |
fhelp(idim) ! additional storage supplied by the calling routine
|
|
|
562 |
|
|
|
563 |
! Local variables
|
|
|
564 |
INTEGER (KIND=iintegers) :: &
|
|
|
565 |
i,m, & ! loop indicees
|
|
|
566 |
nf_o2 ! nfilt/2
|
|
|
567 |
|
|
|
568 |
REAL (KIND=isingle) :: &
|
|
|
569 |
fw(9) ! filter weights
|
|
|
570 |
|
|
|
571 |
!------------------------------------------------------------------------------
|
|
|
572 |
DATA fw /-0.0000152590_isingle, 0.0002441406_isingle, -0.0018310546_isingle, &
|
|
|
573 |
0.0085449218_isingle, -0.0277709960_isingle, 0.0666503906_isingle, &
|
|
|
574 |
-0.1221923828_isingle, 0.1745605469_isingle, 0.8036193848_isingle /
|
|
|
575 |
|
|
|
576 |
! begin subroutine dfilt8_single
|
|
|
577 |
|
|
|
578 |
nf_o2 = (nfilt+1)/2
|
|
|
579 |
|
|
|
580 |
fout (:) = fin(:)
|
|
|
581 |
fhelp(:) = fin(:)
|
|
|
582 |
|
|
|
583 |
DO i = 2, idim-1
|
|
|
584 |
fhelp(i) = 0.25_isingle*fout (i-1) + 0.5_isingle*fout (i) &
|
|
|
585 |
+ 0.25_isingle*fout (i+1)
|
|
|
586 |
ENDDO
|
|
|
587 |
DO i = 2, idim-1
|
|
|
588 |
fout (i) = 0.25_isingle*fhelp(i-1) + 0.5_isingle*fhelp(i) &
|
|
|
589 |
+ 0.25_isingle*fhelp(i+1)
|
|
|
590 |
ENDDO
|
|
|
591 |
|
|
|
592 |
DO m = 1, nf_o2
|
|
|
593 |
DO i = 9, idim-8
|
|
|
594 |
fhelp(i) = fw(9)*fout(i) &
|
|
|
595 |
+ fw(8)*(fout(i-1)+fout(i+1)) + fw(7)*(fout(i-2)+fout(i+2)) &
|
|
|
596 |
+ fw(6)*(fout(i-3)+fout(i+3)) + fw(5)*(fout(i-4)+fout(i+4)) &
|
|
|
597 |
+ fw(4)*(fout(i-5)+fout(i+5)) + fw(3)*(fout(i-6)+fout(i+6)) &
|
|
|
598 |
+ fw(2)*(fout(i-7)+fout(i+7)) + fw(1)*(fout(i-8)+fout(i+8))
|
|
|
599 |
ENDDO
|
|
|
600 |
DO i = 9, idim-8
|
|
|
601 |
fout(i) = fw(9)*fhelp(i) &
|
|
|
602 |
+ fw(8)*(fhelp(i-1)+fhelp(i+1)) + fw(7)*(fhelp(i-2)+fhelp(i+2)) &
|
|
|
603 |
+ fw(6)*(fhelp(i-3)+fhelp(i+3)) + fw(5)*(fhelp(i-4)+fhelp(i+4)) &
|
|
|
604 |
+ fw(4)*(fhelp(i-5)+fhelp(i+5)) + fw(3)*(fhelp(i-6)+fhelp(i+6)) &
|
|
|
605 |
+ fw(2)*(fhelp(i-7)+fhelp(i+7)) + fw(1)*(fhelp(i-8)+fhelp(i+8))
|
|
|
606 |
ENDDO
|
|
|
607 |
ENDDO
|
|
|
608 |
|
|
|
609 |
DO i = 2, idim-1
|
|
|
610 |
fhelp(i) = 0.25_isingle*fout (i-1) + 0.5_isingle*fout (i) &
|
|
|
611 |
+ 0.25_isingle*fout (i+1)
|
|
|
612 |
ENDDO
|
|
|
613 |
DO i = 2, idim-1
|
|
|
614 |
fout (i) = 0.25_isingle*fhelp(i-1) + 0.5_isingle*fhelp(i) &
|
|
|
615 |
+ 0.25_isingle*fhelp(i+1)
|
|
|
616 |
ENDDO
|
|
|
617 |
|
|
|
618 |
END SUBROUTINE dfilt8_single
|
|
|
619 |
|
|
|
620 |
!==============================================================================
|
|
|
621 |
!==============================================================================
|
|
|
622 |
!------------------------------------------------------------------------------
|
|
|
623 |
|
|
|
624 |
SUBROUTINE dolph (deltat, taus, m, window, t, time, time2, w, w2)
|
|
|
625 |
|
|
|
626 |
!------------------------------------------------------------------------------
|
|
|
627 |
!
|
|
|
628 |
! Description:
|
|
|
629 |
! Calculation of Dolph-Chebyshev window or, for short, Dolph Window, using
|
|
|
630 |
! the expression in the reference:
|
|
|
631 |
! Antoniou, Andreas, 1993: Digital Filters: Analysis,
|
|
|
632 |
! Design and Applications. McGraw-Hill, Inc., 689pp.
|
|
|
633 |
!
|
|
|
634 |
! The Dolph window is optimal in the following sense:
|
|
|
635 |
! For a given main-lobe width, the stop-band attenuation is minimal;
|
|
|
636 |
! for a given stop-band level, the main-lobe width is minimal.
|
|
|
637 |
!
|
|
|
638 |
! Method:
|
|
|
639 |
!
|
|
|
640 |
! Modules used: NONE
|
|
|
641 |
!
|
|
|
642 |
!------------------------------------------------------------------------------
|
|
|
643 |
|
|
|
644 |
! Parameter List:
|
|
|
645 |
! ---------------
|
|
|
646 |
|
|
|
647 |
INTEGER (KIND=iintegers), INTENT (IN) :: &
|
|
|
648 |
m ! for dimensioning the work arrays
|
|
|
649 |
|
|
|
650 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
651 |
deltat, taus ! time step and cutoff period for filtering
|
|
|
652 |
|
|
|
653 |
REAL (KIND=ireals), INTENT (OUT) :: &
|
|
|
654 |
window(0:2*m) ! result
|
|
|
655 |
|
|
|
656 |
! The following variables are only used for work space
|
|
|
657 |
REAL (KIND=ireals), INTENT (OUT) :: &
|
|
|
658 |
t(0:2*m), time(0:2*m), time2(0:2*m), w(0:2*m), w2(0:2*m)
|
|
|
659 |
|
|
|
660 |
! Local Variables:
|
|
|
661 |
! ----------------
|
|
|
662 |
|
|
|
663 |
INTEGER (KIND=iintegers) :: nt, i, n, nm1, nn
|
|
|
664 |
REAL (KIND=ireals) :: zpi, zthetas, zx0, zarg, zterm1, zterm2, zrr, &
|
|
|
665 |
zr, zdb, zsum, zsumw
|
|
|
666 |
|
|
|
667 |
!------------ End of header ---------------------------------------------------
|
|
|
668 |
|
|
|
669 |
! Begin subroutine dolph
|
|
|
670 |
|
|
|
671 |
zpi = 4.0_ireals * ATAN(1.0_ireals)
|
|
|
672 |
|
|
|
673 |
n = 2*m+1
|
|
|
674 |
nm1 = n-1
|
|
|
675 |
zthetas = 2.0_ireals*zpi*deltat/taus
|
|
|
676 |
zx0 = 1.0_ireals / COS(zthetas/2.0_ireals)
|
|
|
677 |
zterm1 = (zx0 + SQRT(zx0**2-1))**(REAL (N-1, ireals))
|
|
|
678 |
zterm2 = (zx0 - SQRT(zx0**2-1))**(REAL (N-1, ireals))
|
|
|
679 |
zrr = 0.5*(zterm1 + zterm2)
|
|
|
680 |
zr = 1/zrr
|
|
|
681 |
zdb = 20.0_ireals * LOG10(zr)
|
|
|
682 |
|
|
|
683 |
!------------------------------------------------------------
|
|
|
684 |
|
|
|
685 |
DO nt = 0, M
|
|
|
686 |
zsum = 1
|
|
|
687 |
DO i = 1, M
|
|
|
688 |
zarg = zx0 * cos(i*zpi/N)
|
|
|
689 |
! Calculate the Chebyshev polynomials
|
|
|
690 |
! Reference: Numerical Recipes, Page 184, recurrence
|
|
|
691 |
! T_n(x) = 2xT_{n-1}(x) - T_{n-2}(x) , n>=2.
|
|
|
692 |
T(0) = 1
|
|
|
693 |
T(1) = zarg
|
|
|
694 |
DO nn=2,nm1
|
|
|
695 |
T(nn) = 2*zarg*T(nn-1) - T(nn-2)
|
|
|
696 |
ENDDO
|
|
|
697 |
zterm1 = T(nm1)
|
|
|
698 |
zterm2 = cos(2*nt*zpi*i/n)
|
|
|
699 |
zsum = zsum + zr*2 * zterm1 * zterm2
|
|
|
700 |
ENDDO
|
|
|
701 |
w(nt) = zsum / n
|
|
|
702 |
TIME(nt) = nt
|
|
|
703 |
ENDDO
|
|
|
704 |
|
|
|
705 |
! Fill in the negative-time values by symmetry.
|
|
|
706 |
DO nt = 0, m
|
|
|
707 |
w2(m+nt) = w(nt)
|
|
|
708 |
w2(m-nt) = w(nt)
|
|
|
709 |
time2(m+nt) = time(nT)
|
|
|
710 |
time2(m-nt) = -time(nT)
|
|
|
711 |
ENDDO
|
|
|
712 |
|
|
|
713 |
! Fill up the array for return
|
|
|
714 |
zsumw = 0.0_ireals
|
|
|
715 |
DO nt = 0, 2*m
|
|
|
716 |
zsumw = zsumw + w2(nt)
|
|
|
717 |
ENDDO
|
|
|
718 |
|
|
|
719 |
DO nt=0,2*m
|
|
|
720 |
WINDOW(nt) = w2(nt)
|
|
|
721 |
ENDDO
|
|
|
722 |
!
|
|
|
723 |
!
|
|
|
724 |
!----------------------------------------------------------
|
|
|
725 |
! PRINT *, (w2(nT), nT=0,2*M)
|
|
|
726 |
!----------------------------------------------------------
|
|
|
727 |
!
|
|
|
728 |
|
|
|
729 |
END SUBROUTINE dolph
|
|
|
730 |
|
|
|
731 |
!==============================================================================
|
|
|
732 |
!==============================================================================
|
|
|
733 |
|
|
|
734 |
!------------------------------------------------------------------------------
|
|
|
735 |
|
|
|
736 |
SUBROUTINE elapsed_time (realtimedif, istat)
|
|
|
737 |
|
|
|
738 |
!------------------------------------------------------------------------------
|
|
|
739 |
!
|
|
|
740 |
! Description:
|
|
|
741 |
! Returns the elapsed wall-clock time in seconds since the last call. On
|
|
|
742 |
! the first call the variables are only initialized. If no system clock is
|
|
|
743 |
! present, an error value of istat=1 will be returned, if the optional
|
|
|
744 |
! argument istat was passed from the calling routine.
|
|
|
745 |
! realtimedif is set to 0 then.
|
|
|
746 |
!
|
|
|
747 |
! Method:
|
|
|
748 |
! The intrinsic function SYSTEM_CLOCK is used, that returns the number of
|
|
|
749 |
! clock counts since some system dependent event in the past (e.g. midnight
|
|
|
750 |
! for a 24-hour system clock). The difference of clock counts since the last
|
|
|
751 |
! call is determined and converted into seconds. The variables "lfirst"
|
|
|
752 |
! and "icountsold" (see below) have to be SAVEd for the next call.
|
|
|
753 |
!
|
|
|
754 |
! Modules used: NONE
|
|
|
755 |
!
|
|
|
756 |
!------------------------------------------------------------------------------
|
|
|
757 |
!
|
|
|
758 |
! Parameter List:
|
|
|
759 |
! ---------------
|
|
|
760 |
|
|
|
761 |
REAL (KIND=ireals), INTENT (OUT) :: &
|
|
|
762 |
realtimedif ! wall-clock time since the last call in seconds
|
|
|
763 |
! (0 if no system-clock is available)
|
|
|
764 |
|
|
|
765 |
INTEGER (KIND=iintegers), INTENT (OUT), OPTIONAL :: &
|
|
|
766 |
istat ! optional argument for error value
|
|
|
767 |
|
|
|
768 |
|
|
|
769 |
! Local Variables:
|
|
|
770 |
! ----------------
|
|
|
771 |
|
|
|
772 |
LOGICAL, SAVE :: lfirst = .TRUE. ! determine whether first call or not
|
|
|
773 |
|
|
|
774 |
INTEGER, SAVE :: icountsold ! number of counts in the last call
|
|
|
775 |
|
|
|
776 |
INTEGER :: icountsnew, & ! number of counts in this call
|
|
|
777 |
ir, im ! other arguments to SYSTEM_CLOCK
|
|
|
778 |
|
|
|
779 |
LOGICAL :: lpres ! if optional argument is present
|
|
|
780 |
|
|
|
781 |
!------------ End of header ---------------------------------------------------
|
|
|
782 |
|
|
|
783 |
! Begin subroutine elapsed_time
|
|
|
784 |
|
|
|
785 |
lpres = PRESENT (istat)
|
|
|
786 |
|
|
|
787 |
CALL SYSTEM_CLOCK ( COUNT=icountsnew, COUNT_RATE=ir, COUNT_MAX=im )
|
|
|
788 |
|
|
|
789 |
IF ( ir /= 0 ) THEN
|
|
|
790 |
! system clock is present
|
|
|
791 |
IF (lpres) THEN
|
|
|
792 |
istat = 0
|
|
|
793 |
ENDIF
|
|
|
794 |
|
|
|
795 |
IF (lfirst) THEN
|
|
|
796 |
! first call: store value for the number of clock counts
|
|
|
797 |
icountsold = icountsnew
|
|
|
798 |
lfirst = .FALSE.
|
|
|
799 |
ELSE
|
|
|
800 |
! convert the clock counts to seconds
|
|
|
801 |
IF ( icountsnew >= icountsold ) THEN
|
|
|
802 |
realtimedif = ( REAL (icountsnew - icountsold, ireals) ) &
|
|
|
803 |
/ REAL (ir,ireals)
|
|
|
804 |
ELSE
|
|
|
805 |
realtimedif = REAL (im- (icountsold-icountsnew ), ireals) &
|
|
|
806 |
/ REAL (ir, ireals)
|
|
|
807 |
ENDIF
|
|
|
808 |
icountsold = icountsnew
|
|
|
809 |
ENDIF
|
|
|
810 |
ELSE
|
|
|
811 |
! no system clock present: set error value
|
|
|
812 |
realtimedif = 0.0
|
|
|
813 |
IF ( lpres ) THEN
|
|
|
814 |
istat = 1
|
|
|
815 |
ENDIF
|
|
|
816 |
ENDIF
|
|
|
817 |
|
|
|
818 |
END SUBROUTINE elapsed_time
|
|
|
819 |
|
|
|
820 |
!==============================================================================
|
|
|
821 |
!==============================================================================
|
|
|
822 |
|
|
|
823 |
!------------------------------------------------------------------------------
|
|
|
824 |
|
|
|
825 |
SUBROUTINE get_utc_date (ntsteps, ystartdate, dt, itype_calendar, &
|
|
|
826 |
yactdate1, yactdate2, nactday, acthour)
|
|
|
827 |
|
|
|
828 |
!------------------------------------------------------------------------------
|
|
|
829 |
!
|
|
|
830 |
! Description:
|
|
|
831 |
! This routine determines the actual date of this forecast step.
|
|
|
832 |
!
|
|
|
833 |
! Method:
|
|
|
834 |
! Using the date of the forecast-start, the number of time steps
|
|
|
835 |
! already performed and the length of the time steps, the actual
|
|
|
836 |
! date is calculated taking leap-years into consideration.
|
|
|
837 |
! The date is given in three different formats.
|
|
|
838 |
!
|
|
|
839 |
! Modules used: NONE
|
|
|
840 |
!
|
|
|
841 |
!------------------------------------------------------------------------------
|
|
|
842 |
!
|
|
|
843 |
! Input Parameter list:
|
|
|
844 |
! ---------------------
|
|
|
845 |
|
|
|
846 |
INTEGER (KIND=iintegers), INTENT(IN) :: &
|
|
|
847 |
itype_calendar, & ! for specifying the calendar used
|
|
|
848 |
ntsteps ! number of actual performed time-steps
|
|
|
849 |
|
|
|
850 |
REAL (KIND=ireals), INTENT(IN) :: &
|
|
|
851 |
dt ! time step in seconds
|
|
|
852 |
|
|
|
853 |
CHARACTER (LEN=10), INTENT(IN) :: &
|
|
|
854 |
ystartdate ! start date of the forecast
|
|
|
855 |
|
|
|
856 |
! Output Parameter list:
|
|
|
857 |
! ----------------------
|
|
|
858 |
|
|
|
859 |
CHARACTER (LEN=10), INTENT(OUT) :: &
|
|
|
860 |
yactdate1 ! actual date in the form yyyymmddhh
|
|
|
861 |
|
|
|
862 |
CHARACTER (LEN=22), INTENT(OUT) :: &
|
|
|
863 |
yactdate2 ! actual date in the form wd dd.mm.yy hh UTC
|
|
|
864 |
|
|
|
865 |
|
|
|
866 |
INTEGER (KIND=iintegers), INTENT(OUT) :: &
|
|
|
867 |
nactday ! day of the year
|
|
|
868 |
|
|
|
869 |
REAL (KIND=ireals), INTENT(OUT) :: &
|
|
|
870 |
acthour ! actual hour of the day
|
|
|
871 |
|
|
|
872 |
! Local variables:
|
|
|
873 |
INTEGER (KIND=iintegers) :: &
|
|
|
874 |
month(12), monthsum(13), ileap, iweek, iy, m, &
|
|
|
875 |
idd, imm, iyy, ihh, iday, imonth, iyear, ihour, immhours, iyyhours, &
|
|
|
876 |
iyear_hours
|
|
|
877 |
|
|
|
878 |
CHARACTER (LEN=3) :: yweek(7)
|
|
|
879 |
|
|
|
880 |
! And for computing the amount of seconds of the whole forecast time,
|
|
|
881 |
! an 8-Byte INTEGER has to be used. Otherwise the computation fails after
|
|
|
882 |
! approx. 68 years!!
|
|
|
883 |
|
|
|
884 |
INTEGER, PARAMETER :: int_dp = KIND(1_8)
|
|
|
885 |
REAL(KIND=ireals) :: zseconds
|
|
|
886 |
|
|
|
887 |
!------------ End of header ---------------------------------------------------
|
|
|
888 |
|
|
|
889 |
! Begin subroutine get_utc_date
|
|
|
890 |
|
|
|
891 |
DATA month / 31 , 28 , 31 , 30 , 31 , 30 , &
|
|
|
892 |
31 , 31 , 30 , 31 , 30 , 31 /
|
|
|
893 |
DATA yweek /'MON', 'TUE', 'WED', 'THU', 'FRI', 'SAT', 'SUN' /
|
|
|
894 |
|
|
|
895 |
|
|
|
896 |
! Statementfunction: ileap(yy) = 0: no leap year,
|
|
|
897 |
! ileap(yy) = 1: leap year
|
|
|
898 |
! corrected version for Gregorian / Proleptic calendar
|
|
|
899 |
! by A. Dobler, CLM Community
|
|
|
900 |
ileap (iy) = IABS( MOD(iy, 4) - 4) / 4 & ! every 4 years is a leapyear
|
|
|
901 |
-IABS( MOD(iy,100) - 100) / 100 & ! every 100 years is no leapyear
|
|
|
902 |
+IABS( MOD(iy,400) - 400) / 400 ! but every 400 years is a leapyear
|
|
|
903 |
|
|
|
904 |
! Divide ystartdate in day, month, year and hour
|
|
|
905 |
! and calculate the sums of days from the beginning of the year to the
|
|
|
906 |
! end of the months
|
|
|
907 |
READ ( ystartdate, '(I4,3I2)' ) iyy, imm, idd, ihh
|
|
|
908 |
|
|
|
909 |
IF (itype_calendar == 0) THEN
|
|
|
910 |
month (2) = 28 + ileap (iyy)
|
|
|
911 |
monthsum(1) = 0
|
|
|
912 |
DO m = 2 , 13
|
|
|
913 |
monthsum(m) = monthsum(m-1) + month(m-1)
|
|
|
914 |
enddo
|
|
|
915 |
ELSEIF (itype_calendar == 1) THEN
|
|
|
916 |
monthsum(1) = 0
|
|
|
917 |
DO m = 2 , 13
|
|
|
918 |
monthsum(m) = monthsum(m-1) + 30
|
|
|
919 |
enddo
|
|
|
920 |
ENDIF
|
|
|
921 |
|
|
|
922 |
! Determine how many hours have passed in this year
|
|
|
923 |
iyyhours = (idd*24) + monthsum(imm)*24 + (ihh-24)
|
|
|
924 |
iyyhours = iyyhours + &
|
|
|
925 |
INT (NINT (ntsteps*dt, int_dp)/3600_int_dp, iintegers)
|
|
|
926 |
|
|
|
927 |
! Take turning of the year into account
|
|
|
928 |
IF (itype_calendar == 0) THEN
|
|
|
929 |
iyear_hours = 8760 + ileap(iyy)*24.0_ireals
|
|
|
930 |
ELSEIF (itype_calendar == 1) THEN
|
|
|
931 |
iyear_hours = 8640
|
|
|
932 |
ENDIF
|
|
|
933 |
|
|
|
934 |
IF (iyyhours < 0) THEN
|
|
|
935 |
iyear = iyy-1
|
|
|
936 |
IF (itype_calendar == 0) THEN
|
|
|
937 |
iyyhours = 8760 + ileap(iyear)*24.0_ireals + iyyhours
|
|
|
938 |
ELSEIF (itype_calendar == 1) THEN
|
|
|
939 |
iyyhours = 8640 + iyyhours
|
|
|
940 |
ENDIF
|
|
|
941 |
ELSE IF (iyyhours >= iyear_hours) THEN
|
|
|
942 |
! Take also into account if the run lasts for several years
|
|
|
943 |
iyear = iyy
|
|
|
944 |
IF (itype_calendar == 0) THEN
|
|
|
945 |
iyear_hours = 8760 + ileap(iyear)*24.0_ireals
|
|
|
946 |
ELSEIF (itype_calendar == 1) THEN
|
|
|
947 |
iyear_hours = 8640
|
|
|
948 |
ENDIF
|
|
|
949 |
|
|
|
950 |
DO WHILE (iyyhours >= iyear_hours)
|
|
|
951 |
iyyhours = iyyhours - iyear_hours
|
|
|
952 |
iyear=iyear+1
|
|
|
953 |
IF (itype_calendar == 0) THEN
|
|
|
954 |
iyear_hours = 8760 + ileap(iyear)*24.0_ireals
|
|
|
955 |
ELSEIF (itype_calendar == 1) THEN
|
|
|
956 |
iyear_hours = 8640
|
|
|
957 |
ENDIF
|
|
|
958 |
ENDDO
|
|
|
959 |
ELSE
|
|
|
960 |
iyear = iyy
|
|
|
961 |
ENDIF
|
|
|
962 |
|
|
|
963 |
! calculate monthsum for actual year
|
|
|
964 |
IF (itype_calendar == 0) THEN
|
|
|
965 |
month (2) = 28 + ileap (iyear)
|
|
|
966 |
monthsum(1) = 0
|
|
|
967 |
DO m = 2 , 13
|
|
|
968 |
monthsum(m) = monthsum(m-1) + month(m-1)
|
|
|
969 |
enddo
|
|
|
970 |
ELSEIF (itype_calendar == 1) THEN
|
|
|
971 |
monthsum(1) = 0
|
|
|
972 |
DO m = 2 , 13
|
|
|
973 |
monthsum(m) = monthsum(m-1) + 30
|
|
|
974 |
enddo
|
|
|
975 |
ENDIF
|
|
|
976 |
|
|
|
977 |
! Determine the actual date from iyyhours
|
|
|
978 |
m = 1
|
|
|
979 |
immhours = iyyhours
|
|
|
980 |
DO WHILE (immhours >= 0)
|
|
|
981 |
m = m+1
|
|
|
982 |
immhours = iyyhours - monthsum(m) * 24
|
|
|
983 |
ENDDO
|
|
|
984 |
imonth = m-1
|
|
|
985 |
|
|
|
986 |
immhours = iyyhours - monthsum(imonth)*24
|
|
|
987 |
iday = immhours/24 + 1
|
|
|
988 |
ihour = MOD(immhours,24)
|
|
|
989 |
|
|
|
990 |
|
|
|
991 |
!US: This was before, when 3600.0/dt was an integer value
|
|
|
992 |
! acthour = REAL (ihour, ireals) + dt/3600._ireals* &
|
|
|
993 |
! MOD(ntsteps,INT(3600./dt+0.01)) + 0.0001_ireals
|
|
|
994 |
|
|
|
995 |
! this is the more accurate computation
|
|
|
996 |
zseconds = REAL (ntsteps, ireals) * dt / 3600.0_ireals
|
|
|
997 |
acthour = REAL (ihour, ireals) + &
|
|
|
998 |
(zseconds - REAL(INT(zseconds, int_dp), ireals))
|
|
|
999 |
|
|
|
1000 |
ihour = INT(acthour)
|
|
|
1001 |
nactday = monthsum(imonth) + iday + INT(acthour/24. + 0.0001)
|
|
|
1002 |
iweek = MOD(monthsum(imonth) + iday + (iyear-1901) + (iyear-1901)/4, 7)+1
|
|
|
1003 |
|
|
|
1004 |
WRITE ( yactdate1(1:4) , '(I4.4)' ) iyear
|
|
|
1005 |
WRITE ( yactdate1(5:6) , '(I2.2)' ) imonth
|
|
|
1006 |
WRITE ( yactdate1(7:8) , '(I2.2)' ) iday
|
|
|
1007 |
WRITE ( yactdate1(9:10), '(I2.2)' ) ihour
|
|
|
1008 |
|
|
|
1009 |
IF (itype_calendar == 0) THEN
|
|
|
1010 |
yactdate2 = yweek(iweek)//' '//yactdate1(7:8)//'.'// yactdate1(5:6)//'.' &
|
|
|
1011 |
//yactdate1(1:4)//' '//yactdate1(9:10)//' UTC'
|
|
|
1012 |
ELSEIF (itype_calendar == 1) THEN
|
|
|
1013 |
yactdate2 = ' '//yactdate1(7:8)//'.'// yactdate1(5:6)//'.' &
|
|
|
1014 |
//yactdate1(1:4)//' '//yactdate1(9:10)//' UTC'
|
|
|
1015 |
ENDIF
|
|
|
1016 |
|
|
|
1017 |
END SUBROUTINE get_utc_date
|
|
|
1018 |
|
|
|
1019 |
!==============================================================================
|
|
|
1020 |
!==============================================================================
|
|
|
1021 |
!+ filter routines with 17-, 13-, 9- and 3-points stencil, resp.
|
|
|
1022 |
!------------------------------------------------------------------------------
|
|
|
1023 |
|
|
|
1024 |
SUBROUTINE horizontal_filtering( field_flt, ie_in, je_in, kedim, &
|
|
|
1025 |
nbdlines, nflt_width, ncutoff, &
|
|
|
1026 |
neighbors, hfx_mask, hfy_mask )
|
|
|
1027 |
|
|
|
1028 |
!------------------------------------------------------------------------------
|
|
|
1029 |
!
|
|
|
1030 |
! Description:
|
|
|
1031 |
!
|
|
|
1032 |
! Method:
|
|
|
1033 |
!
|
|
|
1034 |
!------------------------------------------------------------------------------
|
|
|
1035 |
|
|
|
1036 |
! Subroutine arguments:
|
|
|
1037 |
! ---------------------
|
|
|
1038 |
INTEGER (KIND=iintegers), INTENT(IN) :: &
|
|
|
1039 |
ie_in, je_in, kedim, & ! Dimensions of the field to be filtered
|
|
|
1040 |
nbdlines, & ! number of boundary lines from decomposition
|
|
|
1041 |
nflt_width, & ! width of field extension for filtering
|
|
|
1042 |
ncutoff, & ! filter-value for cutoff
|
|
|
1043 |
neighbors(4) ! process-id's of the neighbors in the grid
|
|
|
1044 |
|
|
|
1045 |
|
|
|
1046 |
REAL (KIND=ireals ), INTENT(INOUT) :: &
|
|
|
1047 |
field_flt(ie_in, je_in, kedim)
|
|
|
1048 |
|
|
|
1049 |
LOGICAL, INTENT(in), OPTIONAL :: &
|
|
|
1050 |
hfx_mask(ie_in, je_in), hfy_mask(ie_in, je_in)
|
|
|
1051 |
|
|
|
1052 |
! Local scalars:
|
|
|
1053 |
! -------------
|
|
|
1054 |
INTEGER (KIND=iintegers) :: &
|
|
|
1055 |
ilow, iup, & !
|
|
|
1056 |
jlow, jup, & !
|
|
|
1057 |
izstata, & ! error status at allocation
|
|
|
1058 |
izstatd, & ! error status at deallocation
|
|
|
1059 |
i, j, k, l ! Loop indices
|
|
|
1060 |
|
|
|
1061 |
INTEGER (KIND=iintegers) :: &
|
|
|
1062 |
istart, iend, jstart, jend, nfw_m_nb
|
|
|
1063 |
|
|
|
1064 |
! Local (automatic) arrays:
|
|
|
1065 |
! -------------------------
|
|
|
1066 |
REAL (KIND=ireals ) :: &
|
|
|
1067 |
field_tmp (ie_in, je_in, kedim), &
|
|
|
1068 |
field_tmp2(ie_in, je_in, kedim), &
|
|
|
1069 |
zfwnp(-nflt_width:nflt_width), & ! filter weights for n-point filter
|
|
|
1070 |
zfw3p(-1:1) ! filter weights for 3-point filter
|
|
|
1071 |
|
|
|
1072 |
!------------------------------------------------------------------------------
|
|
|
1073 |
|
|
|
1074 |
nfw_m_nb = nflt_width - nbdlines
|
|
|
1075 |
istart = 1 + nbdlines
|
|
|
1076 |
iend = ie_in - 2*nfw_m_nb - nbdlines
|
|
|
1077 |
jstart = 1 + nbdlines
|
|
|
1078 |
jend = je_in - 2*nfw_m_nb - nbdlines
|
|
|
1079 |
|
|
|
1080 |
! filter weights for n-point filter
|
|
|
1081 |
IF (ncutoff == 3 .AND. nflt_width == 4) THEN
|
|
|
1082 |
! --> dfilt4
|
|
|
1083 |
! filter weights for 9-point filter (approx. cutoff = 3)
|
|
|
1084 |
zfwnp = (/ -0.390625E-02_ireals, &
|
|
|
1085 |
+0.3125E-01_ireals, &
|
|
|
1086 |
-0.109375_ireals, &
|
|
|
1087 |
+0.21875_ireals, &
|
|
|
1088 |
+0.7265625_ireals, &
|
|
|
1089 |
+0.21875_ireals, &
|
|
|
1090 |
-0.109375_ireals, &
|
|
|
1091 |
+0.3125E-01_ireals, &
|
|
|
1092 |
-0.390625E-02_ireals /)
|
|
|
1093 |
ELSEIF (ncutoff == 3 .AND. nflt_width == 8) THEN
|
|
|
1094 |
! --> dfilt8
|
|
|
1095 |
! filter weights for 17-point filter (approx. cutoff = 3)
|
|
|
1096 |
zfwnp = (/ -0.15259E-04_ireals, &
|
|
|
1097 |
+0.2441406E-03_ireals, &
|
|
|
1098 |
-0.18310546E-02_ireals, &
|
|
|
1099 |
+0.85449218E-02_ireals, &
|
|
|
1100 |
-0.27770996E-01_ireals, &
|
|
|
1101 |
+0.666503906E-01_ireals, &
|
|
|
1102 |
-0.1221923828_ireals, &
|
|
|
1103 |
+0.1745605469_ireals, &
|
|
|
1104 |
+0.8036193848_ireals, &
|
|
|
1105 |
+0.1745605469_ireals, &
|
|
|
1106 |
-0.1221923828_ireals, &
|
|
|
1107 |
+0.666503906E-01_ireals, &
|
|
|
1108 |
-0.27770996E-01_ireals, &
|
|
|
1109 |
+0.85449218E-02_ireals, &
|
|
|
1110 |
-0.18310546E-02_ireals, &
|
|
|
1111 |
+0.2441406E-03_ireals, &
|
|
|
1112 |
-0.15259E-04_ireals /)
|
|
|
1113 |
ELSEIF (ncutoff == 4 .AND. nflt_width == 4) THEN
|
|
|
1114 |
! filter weights for 9-point filter (approx. cutoff = 4)
|
|
|
1115 |
zfwnp = (/ +0.1171875E-01_ireals, &
|
|
|
1116 |
-0.3125E-01_ireals, &
|
|
|
1117 |
-0.46875E-01_ireals, &
|
|
|
1118 |
+0.28125_ireals, &
|
|
|
1119 |
+0.5703125_ireals, &
|
|
|
1120 |
+0.28125_ireals, &
|
|
|
1121 |
-0.46875E-01_ireals, &
|
|
|
1122 |
-0.3125E-01_ireals, &
|
|
|
1123 |
+0.1171875E-01_ireals /)
|
|
|
1124 |
ELSEIF (ncutoff == 5 .AND. nflt_width == 6) THEN
|
|
|
1125 |
! filter weights for 13-point filter (approx. cutoff = 5)
|
|
|
1126 |
zfwnp = (/ +0.44023278E-02_ireals, &
|
|
|
1127 |
+0.13175894E-01_ireals, &
|
|
|
1128 |
-0.477203075E-01_ireals, &
|
|
|
1129 |
-0.435555245E-01_ireals, &
|
|
|
1130 |
+0.94700467E-01_ireals, &
|
|
|
1131 |
+0.2888298641_ireals, &
|
|
|
1132 |
+0.3803345582_ireals, &
|
|
|
1133 |
+0.2888298641_ireals, &
|
|
|
1134 |
+0.94700467E-01_ireals, &
|
|
|
1135 |
-0.435555245E-01_ireals, &
|
|
|
1136 |
-0.477203075E-01_ireals, &
|
|
|
1137 |
+0.13175894E-01_ireals, &
|
|
|
1138 |
+0.44023278E-02_ireals /)
|
|
|
1139 |
ELSEIF (ncutoff == 6 .AND. nflt_width == 4) THEN
|
|
|
1140 |
! filter weights for 9-point filter (approx. cutoff = 6)
|
|
|
1141 |
zfwnp = (/ -0.4694126E-01_ireals, &
|
|
|
1142 |
-0.50095541E-02_ireals, &
|
|
|
1143 |
+0.13528415_ireals, &
|
|
|
1144 |
+0.25500955_ireals, &
|
|
|
1145 |
+0.32331423_ireals, &
|
|
|
1146 |
+0.25500955_ireals, &
|
|
|
1147 |
+0.13528415_ireals, &
|
|
|
1148 |
-0.50095541E-02_ireals, &
|
|
|
1149 |
-0.4694126E-01_ireals /)
|
|
|
1150 |
ELSEIF (ncutoff == 8 .AND. nflt_width == 6) THEN
|
|
|
1151 |
! filter weights for 13-point filter (approx. cutoff = 8)
|
|
|
1152 |
zfwnp = (/ -0.16638111E-01_ireals, &
|
|
|
1153 |
-0.30753028E-01_ireals, &
|
|
|
1154 |
-0.17361869E-02_ireals, &
|
|
|
1155 |
+0.65428931E-01_ireals, &
|
|
|
1156 |
+0.14784805_ireals, &
|
|
|
1157 |
+0.2153241_ireals, &
|
|
|
1158 |
+0.2410525_ireals, &
|
|
|
1159 |
+0.2153241_ireals, &
|
|
|
1160 |
+0.14784805_ireals, &
|
|
|
1161 |
+0.65428931E-01_ireals, &
|
|
|
1162 |
-0.17361869E-02_ireals, &
|
|
|
1163 |
-0.30753028E-01_ireals, &
|
|
|
1164 |
-0.16638111E-01_ireals /)
|
|
|
1165 |
ELSE
|
|
|
1166 |
PRINT *, ' ERROR *** Wrong cutoff value for filtering or *** '
|
|
|
1167 |
PRINT *, ' ERROR *** wrong value for filter/field extension. *** '
|
|
|
1168 |
ENDIF
|
|
|
1169 |
|
|
|
1170 |
! filter weights for 3-point filter (approx. cutoff = 4)
|
|
|
1171 |
zfw3p = (/ 0.25_ireals, 0.5_ireals, 0.25_ireals /)
|
|
|
1172 |
|
|
|
1173 |
! west
|
|
|
1174 |
IF (neighbors(1) == -1) THEN
|
|
|
1175 |
ilow = 1 + 2*nflt_width
|
|
|
1176 |
ELSE
|
|
|
1177 |
ilow = istart + nfw_m_nb
|
|
|
1178 |
END IF
|
|
|
1179 |
! east
|
|
|
1180 |
IF (neighbors(3) == -1) THEN
|
|
|
1181 |
iup = iend - nbdlines
|
|
|
1182 |
ELSE
|
|
|
1183 |
iup = iend + nfw_m_nb
|
|
|
1184 |
END IF
|
|
|
1185 |
! south
|
|
|
1186 |
IF (neighbors(4) == -1) THEN
|
|
|
1187 |
jlow = 1 + 2*nflt_width
|
|
|
1188 |
ELSE
|
|
|
1189 |
jlow = jstart + nfw_m_nb
|
|
|
1190 |
END IF
|
|
|
1191 |
! north
|
|
|
1192 |
IF (neighbors(2) == -1) THEN
|
|
|
1193 |
jup = jend - nbdlines
|
|
|
1194 |
ELSE
|
|
|
1195 |
jup = jend + nfw_m_nb
|
|
|
1196 |
END IF
|
|
|
1197 |
|
|
|
1198 |
! init working array
|
|
|
1199 |
field_tmp (:,:,:) = field_flt(:,:,:)
|
|
|
1200 |
|
|
|
1201 |
|
|
|
1202 |
IF ( PRESENT( hfx_mask ) ) THEN
|
|
|
1203 |
|
|
|
1204 |
! apply n-point-filter in x-direction
|
|
|
1205 |
DO k = 1, kedim
|
|
|
1206 |
DO j = 1, je_in
|
|
|
1207 |
DO i = ilow, iup
|
|
|
1208 |
IF ( hfx_mask(i,j) ) THEN
|
|
|
1209 |
field_tmp(i,j,k) = 0.0_ireals
|
|
|
1210 |
ENDIF
|
|
|
1211 |
ENDDO
|
|
|
1212 |
DO l = -nflt_width, nflt_width
|
|
|
1213 |
DO i = ilow, iup
|
|
|
1214 |
IF ( hfx_mask(i,j) ) THEN
|
|
|
1215 |
field_tmp(i,j,k) = field_tmp(i,j,k) &
|
|
|
1216 |
+ zfwnp(l)*field_flt(i+l,j,k)
|
|
|
1217 |
END IF
|
|
|
1218 |
END DO
|
|
|
1219 |
END DO
|
|
|
1220 |
END DO
|
|
|
1221 |
END DO
|
|
|
1222 |
|
|
|
1223 |
! apply 3-point-filter in x-direction at west boundary
|
|
|
1224 |
IF (neighbors(1) == -1) THEN
|
|
|
1225 |
DO k = 1, kedim
|
|
|
1226 |
DO j = 1, je_in
|
|
|
1227 |
DO i = nfw_m_nb+1, ilow-1
|
|
|
1228 |
IF ( hfx_mask(i,j) ) THEN
|
|
|
1229 |
field_tmp(i,j,k) = 0.0_ireals
|
|
|
1230 |
ENDIF
|
|
|
1231 |
ENDDO
|
|
|
1232 |
DO l = -1, 1
|
|
|
1233 |
DO i = nfw_m_nb+1, ilow-1
|
|
|
1234 |
IF ( hfx_mask(i,j) ) THEN
|
|
|
1235 |
field_tmp(i,j,k) = field_tmp(i,j,k) &
|
|
|
1236 |
+ zfw3p(l)*field_flt(i+l,j,k)
|
|
|
1237 |
END IF
|
|
|
1238 |
END DO
|
|
|
1239 |
END DO
|
|
|
1240 |
END DO
|
|
|
1241 |
END DO
|
|
|
1242 |
END IF
|
|
|
1243 |
|
|
|
1244 |
! apply 3-point-filter in x-direction at east boundary
|
|
|
1245 |
IF (neighbors(3) == -1) THEN
|
|
|
1246 |
DO k = 1, kedim
|
|
|
1247 |
DO j = 1, je_in
|
|
|
1248 |
DO i = iup+1, ie_in-nfw_m_nb
|
|
|
1249 |
IF ( hfx_mask(i,j) ) THEN
|
|
|
1250 |
field_tmp(i,j,k) = 0.0_ireals
|
|
|
1251 |
ENDIF
|
|
|
1252 |
ENDDO
|
|
|
1253 |
DO l = -1, 1
|
|
|
1254 |
DO i = iup+1, ie_in-nfw_m_nb
|
|
|
1255 |
IF ( hfx_mask(i,j) ) THEN
|
|
|
1256 |
field_tmp(i,j,k) = field_tmp(i,j,k) &
|
|
|
1257 |
+ zfw3p(l)*field_flt(i+l,j,k)
|
|
|
1258 |
END IF
|
|
|
1259 |
END DO
|
|
|
1260 |
END DO
|
|
|
1261 |
END DO
|
|
|
1262 |
END DO
|
|
|
1263 |
END IF
|
|
|
1264 |
|
|
|
1265 |
ELSE
|
|
|
1266 |
|
|
|
1267 |
!
|
|
|
1268 |
! apply n-point-filter in x-direction
|
|
|
1269 |
!
|
|
|
1270 |
DO k = 1, kedim
|
|
|
1271 |
DO j = 1, je_in
|
|
|
1272 |
DO i = ilow, iup
|
|
|
1273 |
field_tmp(i,j,k) = 0.0_ireals
|
|
|
1274 |
END DO
|
|
|
1275 |
DO l = -nflt_width, nflt_width
|
|
|
1276 |
DO i = ilow, iup
|
|
|
1277 |
field_tmp(i,j,k) = field_tmp(i,j,k) &
|
|
|
1278 |
+ zfwnp(l)*field_flt(i+l,j,k)
|
|
|
1279 |
END DO
|
|
|
1280 |
END DO
|
|
|
1281 |
END DO
|
|
|
1282 |
END DO
|
|
|
1283 |
|
|
|
1284 |
! apply 3-point-filter in x-direction at west boundary
|
|
|
1285 |
IF (neighbors(1) == -1) THEN
|
|
|
1286 |
DO k = 1, kedim
|
|
|
1287 |
DO j = 1, je_in
|
|
|
1288 |
DO i = nfw_m_nb+1, ilow-1
|
|
|
1289 |
field_tmp(i,j,k) = 0.0_ireals
|
|
|
1290 |
END DO
|
|
|
1291 |
DO l = -1, 1
|
|
|
1292 |
DO i = nfw_m_nb+1, ilow-1
|
|
|
1293 |
field_tmp(i,j,k) = field_tmp(i,j,k) &
|
|
|
1294 |
+ zfw3p(l)*field_flt(i+l,j,k)
|
|
|
1295 |
END DO
|
|
|
1296 |
END DO
|
|
|
1297 |
END DO
|
|
|
1298 |
END DO
|
|
|
1299 |
END IF
|
|
|
1300 |
|
|
|
1301 |
! apply 3-point-filter in x-direction at east boundary
|
|
|
1302 |
IF (neighbors(3) == -1) THEN
|
|
|
1303 |
DO k = 1, kedim
|
|
|
1304 |
DO j = 1, je_in
|
|
|
1305 |
DO i = iup+1, ie_in-nfw_m_nb
|
|
|
1306 |
field_tmp(i,j,k) = 0.0_ireals
|
|
|
1307 |
END DO
|
|
|
1308 |
DO l = -1, 1
|
|
|
1309 |
DO i = iup+1, ie_in-nfw_m_nb
|
|
|
1310 |
field_tmp(i,j,k) = field_tmp(i,j,k) &
|
|
|
1311 |
+ zfw3p(l)*field_flt(i+l,j,k)
|
|
|
1312 |
END DO
|
|
|
1313 |
END DO
|
|
|
1314 |
END DO
|
|
|
1315 |
END DO
|
|
|
1316 |
END IF
|
|
|
1317 |
|
|
|
1318 |
END IF
|
|
|
1319 |
|
|
|
1320 |
|
|
|
1321 |
IF ( PRESENT( hfy_mask ) ) THEN
|
|
|
1322 |
|
|
|
1323 |
! apply n-point-filter in y-direction
|
|
|
1324 |
DO k = 1, kedim
|
|
|
1325 |
DO j = jlow, jup
|
|
|
1326 |
DO i = 1, ie_in
|
|
|
1327 |
IF ( hfy_mask(i,j) ) THEN
|
|
|
1328 |
field_flt(i,j,k) = 0.0_ireals
|
|
|
1329 |
ELSE
|
|
|
1330 |
field_flt(i,j,k) = field_tmp(i,j,k)
|
|
|
1331 |
ENDIF
|
|
|
1332 |
ENDDO
|
|
|
1333 |
DO l = -nflt_width, nflt_width
|
|
|
1334 |
DO i = 1, ie_in
|
|
|
1335 |
IF ( hfy_mask(i,j) ) THEN
|
|
|
1336 |
field_flt(i,j,k) = field_flt(i,j,k) + zfwnp(l)*field_tmp(i,j+l,k)
|
|
|
1337 |
END IF
|
|
|
1338 |
END DO
|
|
|
1339 |
END DO
|
|
|
1340 |
END DO
|
|
|
1341 |
END DO
|
|
|
1342 |
|
|
|
1343 |
! apply 3-point-filter in y-direction at south boundary
|
|
|
1344 |
IF (neighbors(4) == -1) THEN
|
|
|
1345 |
DO k = 1, kedim
|
|
|
1346 |
DO j = nfw_m_nb+1, jlow-1
|
|
|
1347 |
DO i = 1, ie_in
|
|
|
1348 |
IF ( hfy_mask(i,j) ) THEN
|
|
|
1349 |
field_flt(i,j,k) = 0.0_ireals
|
|
|
1350 |
ELSE
|
|
|
1351 |
field_flt(i,j,k) = field_tmp(i,j,k)
|
|
|
1352 |
ENDIF
|
|
|
1353 |
ENDDO
|
|
|
1354 |
DO l = -1, 1
|
|
|
1355 |
DO i = 1, ie_in
|
|
|
1356 |
IF ( hfy_mask(i,j) ) THEN
|
|
|
1357 |
field_flt(i,j,k) = field_flt(i,j,k)+zfw3p(l)*field_tmp(i,j+l,k)
|
|
|
1358 |
END IF
|
|
|
1359 |
END DO
|
|
|
1360 |
END DO
|
|
|
1361 |
END DO
|
|
|
1362 |
END DO
|
|
|
1363 |
END IF
|
|
|
1364 |
|
|
|
1365 |
! apply 3-point-filter in y-direction at north boundary
|
|
|
1366 |
IF (neighbors(2) == -1) THEN
|
|
|
1367 |
DO k = 1, kedim
|
|
|
1368 |
DO j = jup+1, je_in-nfw_m_nb
|
|
|
1369 |
DO i = 1, ie_in
|
|
|
1370 |
IF ( hfy_mask(i,j) ) THEN
|
|
|
1371 |
field_flt(i,j,k) = 0.0_ireals
|
|
|
1372 |
ELSE
|
|
|
1373 |
field_flt(i,j,k) = field_tmp(i,j,k)
|
|
|
1374 |
ENDIF
|
|
|
1375 |
ENDDO
|
|
|
1376 |
DO l = -1, 1
|
|
|
1377 |
DO i = 1, ie_in
|
|
|
1378 |
IF ( hfy_mask(i,j) ) THEN
|
|
|
1379 |
field_flt(i,j,k) = field_flt(i,j,k)+zfw3p(l)*field_tmp(i,j+l,k)
|
|
|
1380 |
END IF
|
|
|
1381 |
END DO
|
|
|
1382 |
END DO
|
|
|
1383 |
END DO
|
|
|
1384 |
END DO
|
|
|
1385 |
END IF
|
|
|
1386 |
|
|
|
1387 |
ELSE
|
|
|
1388 |
|
|
|
1389 |
!
|
|
|
1390 |
! apply n-point-filter in y-direction
|
|
|
1391 |
!
|
|
|
1392 |
DO k = 1, kedim
|
|
|
1393 |
DO j = jlow, jup
|
|
|
1394 |
DO i = 1, ie_in
|
|
|
1395 |
field_flt(i,j,k) = 0.0_ireals
|
|
|
1396 |
ENDDO
|
|
|
1397 |
DO l = -nflt_width, nflt_width
|
|
|
1398 |
DO i = 1, ie_in
|
|
|
1399 |
field_flt(i,j,k) = field_flt(i,j,k)+zfwnp(l)*field_tmp(i,j+l,k)
|
|
|
1400 |
END DO
|
|
|
1401 |
END DO
|
|
|
1402 |
END DO
|
|
|
1403 |
END DO
|
|
|
1404 |
|
|
|
1405 |
! apply 3-point-filter in y-direction at south boundary
|
|
|
1406 |
IF (neighbors(4) == -1) THEN
|
|
|
1407 |
DO k = 1, kedim
|
|
|
1408 |
DO j = nfw_m_nb+1, jlow-1
|
|
|
1409 |
DO i = 1, ie_in
|
|
|
1410 |
field_flt(i,j,k) = 0.0_ireals
|
|
|
1411 |
ENDDO
|
|
|
1412 |
DO l = -1, 1
|
|
|
1413 |
DO i = 1, ie_in
|
|
|
1414 |
field_flt(i,j,k) = field_flt(i,j,k)+zfw3p(l)*field_tmp(i,j+l,k)
|
|
|
1415 |
END DO
|
|
|
1416 |
END DO
|
|
|
1417 |
END DO
|
|
|
1418 |
END DO
|
|
|
1419 |
END IF
|
|
|
1420 |
|
|
|
1421 |
! apply 3-point-filter in y-direction at north boundary
|
|
|
1422 |
IF (neighbors(2) == -1) THEN
|
|
|
1423 |
DO k = 1, kedim
|
|
|
1424 |
DO j = jup+1, je_in-nfw_m_nb
|
|
|
1425 |
DO i = 1, ie_in
|
|
|
1426 |
field_flt(i,j,k) = 0.0_ireals
|
|
|
1427 |
ENDDO
|
|
|
1428 |
DO l = -1, 1
|
|
|
1429 |
DO i = 1, ie_in
|
|
|
1430 |
field_flt(i,j,k) = field_flt(i,j,k)+zfw3p(l)*field_tmp(i,j+l,k)
|
|
|
1431 |
END DO
|
|
|
1432 |
END DO
|
|
|
1433 |
END DO
|
|
|
1434 |
END DO
|
|
|
1435 |
END IF
|
|
|
1436 |
|
|
|
1437 |
END IF
|
|
|
1438 |
|
|
|
1439 |
!------------------------------------------------------------------------------
|
|
|
1440 |
|
|
|
1441 |
!------------------------------------------------------------------------------
|
|
|
1442 |
! End of subroutine horizontal_filtering
|
|
|
1443 |
!------------------------------------------------------------------------------
|
|
|
1444 |
|
|
|
1445 |
END SUBROUTINE horizontal_filtering
|
|
|
1446 |
|
|
|
1447 |
!==============================================================================
|
|
|
1448 |
!==============================================================================
|
|
|
1449 |
!+ Function for rotation of geographical coordinates
|
|
|
1450 |
!------------------------------------------------------------------------------
|
|
|
1451 |
|
|
|
1452 |
FUNCTION phirot2phi ( phirot, rlarot, polphi, pollam, polgam )
|
|
|
1453 |
|
|
|
1454 |
!------------------------------------------------------------------------------
|
|
|
1455 |
!
|
|
|
1456 |
! Description:
|
|
|
1457 |
! This function converts phi from one rotated system to phi in another
|
|
|
1458 |
! system. If the optional argument polgam is present, the other system
|
|
|
1459 |
! can also be a rotated one, where polgam is the angle between the two
|
|
|
1460 |
! north poles.
|
|
|
1461 |
! If polgam is not present, the other system is the real geographical
|
|
|
1462 |
! system.
|
|
|
1463 |
!
|
|
|
1464 |
! Method:
|
|
|
1465 |
! Transformation formulas for converting between these two systems.
|
|
|
1466 |
!
|
|
|
1467 |
!------------------------------------------------------------------------------
|
|
|
1468 |
|
|
|
1469 |
!------------------------------------------------------------------------------
|
|
|
1470 |
!
|
|
|
1471 |
! Declarations:
|
|
|
1472 |
!
|
|
|
1473 |
!------------------------------------------------------------------------------
|
|
|
1474 |
|
|
|
1475 |
! Parameter list:
|
|
|
1476 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
1477 |
polphi, & ! latitude of the rotated north pole
|
|
|
1478 |
pollam, & ! longitude of the rotated north pole
|
|
|
1479 |
phirot, & ! latitude in the rotated system
|
|
|
1480 |
rlarot ! longitude in the rotated system
|
|
|
1481 |
|
|
|
1482 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
1483 |
polgam ! angle between the north poles of the systems
|
|
|
1484 |
|
|
|
1485 |
REAL (KIND=ireals) :: &
|
|
|
1486 |
phirot2phi ! latitude in the geographical system
|
|
|
1487 |
|
|
|
1488 |
! Local variables
|
|
|
1489 |
REAL (KIND=ireals) :: &
|
|
|
1490 |
zsinpol, zcospol, zphis, zrlas, zarg, zgam
|
|
|
1491 |
|
|
|
1492 |
REAL (KIND=ireals), PARAMETER :: &
|
|
|
1493 |
zrpi18 = 57.2957795_ireals, &
|
|
|
1494 |
zpir18 = 0.0174532925_ireals
|
|
|
1495 |
|
|
|
1496 |
!------------------------------------------------------------------------------
|
|
|
1497 |
|
|
|
1498 |
! Begin function phirot2phi
|
|
|
1499 |
|
|
|
1500 |
zsinpol = SIN (zpir18 * polphi)
|
|
|
1501 |
zcospol = COS (zpir18 * polphi)
|
|
|
1502 |
|
|
|
1503 |
zphis = zpir18 * phirot
|
|
|
1504 |
IF (rlarot > 180.0_ireals) THEN
|
|
|
1505 |
zrlas = rlarot - 360.0_ireals
|
|
|
1506 |
ELSE
|
|
|
1507 |
zrlas = rlarot
|
|
|
1508 |
ENDIF
|
|
|
1509 |
zrlas = zpir18 * zrlas
|
|
|
1510 |
|
|
|
1511 |
IF (polgam /= 0.0_ireals) THEN
|
|
|
1512 |
zgam = zpir18 * polgam
|
|
|
1513 |
zarg = zsinpol*SIN (zphis) + &
|
|
|
1514 |
zcospol*COS(zphis) * ( COS(zrlas)*COS(zgam) - SIN(zgam)*SIN(zrlas) )
|
|
|
1515 |
ELSE
|
|
|
1516 |
zarg = zcospol * COS (zphis) * COS (zrlas) + zsinpol * SIN (zphis)
|
|
|
1517 |
ENDIF
|
|
|
1518 |
|
|
|
1519 |
phirot2phi = zrpi18 * ASIN (zarg)
|
|
|
1520 |
|
|
|
1521 |
END FUNCTION phirot2phi
|
|
|
1522 |
|
|
|
1523 |
!==============================================================================
|
|
|
1524 |
!==============================================================================
|
|
|
1525 |
|
|
|
1526 |
!------------------------------------------------------------------------------
|
|
|
1527 |
|
|
|
1528 |
FUNCTION phi2phirot ( phi, rla, polphi, pollam )
|
|
|
1529 |
|
|
|
1530 |
!------------------------------------------------------------------------------
|
|
|
1531 |
! Description:
|
|
|
1532 |
! This routine converts phi from the real geographical system to phi
|
|
|
1533 |
! in the rotated system.
|
|
|
1534 |
!
|
|
|
1535 |
! Method:
|
|
|
1536 |
! Transformation formulas for converting between these two systems.
|
|
|
1537 |
!
|
|
|
1538 |
!------------------------------------------------------------------------------
|
|
|
1539 |
! Parameter list:
|
|
|
1540 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
1541 |
polphi, & ! latitude of the rotated north pole
|
|
|
1542 |
pollam, & ! longitude of the rotated north pole
|
|
|
1543 |
phi, & ! latitude in the geographical system
|
|
|
1544 |
rla ! longitude in the geographical system
|
|
|
1545 |
|
|
|
1546 |
REAL (KIND=ireals) :: &
|
|
|
1547 |
phi2phirot ! longitude in the rotated system
|
|
|
1548 |
|
|
|
1549 |
! Local variables
|
|
|
1550 |
REAL (KIND=ireals) :: &
|
|
|
1551 |
zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
|
|
|
1552 |
|
|
|
1553 |
REAL (KIND=ireals), PARAMETER :: &
|
|
|
1554 |
zrpi18 = 57.2957795_ireals, & !
|
|
|
1555 |
zpir18 = 0.0174532925_ireals
|
|
|
1556 |
|
|
|
1557 |
!------------------------------------------------------------------------------
|
|
|
1558 |
|
|
|
1559 |
! Begin function phi2phirot
|
|
|
1560 |
|
|
|
1561 |
zsinpol = SIN (zpir18 * polphi)
|
|
|
1562 |
zcospol = COS (zpir18 * polphi)
|
|
|
1563 |
zlampol = zpir18 * pollam
|
|
|
1564 |
zphi = zpir18 * phi
|
|
|
1565 |
IF (rla > 180.0_ireals) THEN
|
|
|
1566 |
zrla1 = rla - 360.0_ireals
|
|
|
1567 |
ELSE
|
|
|
1568 |
zrla1 = rla
|
|
|
1569 |
ENDIF
|
|
|
1570 |
zrla = zpir18 * zrla1
|
|
|
1571 |
|
|
|
1572 |
zarg1 = SIN (zphi) * zsinpol
|
|
|
1573 |
zarg2 = COS (zphi) * zcospol * COS (zrla - zlampol)
|
|
|
1574 |
|
|
|
1575 |
phi2phirot = zrpi18 * ASIN (zarg1 + zarg2)
|
|
|
1576 |
|
|
|
1577 |
END FUNCTION phi2phirot
|
|
|
1578 |
|
|
|
1579 |
!==============================================================================
|
|
|
1580 |
!==============================================================================
|
|
|
1581 |
|
|
|
1582 |
!------------------------------------------------------------------------------
|
|
|
1583 |
|
|
|
1584 |
FUNCTION rlarot2rla (phirot, rlarot, polphi, pollam, polgam)
|
|
|
1585 |
|
|
|
1586 |
!------------------------------------------------------------------------------
|
|
|
1587 |
!
|
|
|
1588 |
! Description:
|
|
|
1589 |
! This function converts lambda from one rotated system to lambda in another
|
|
|
1590 |
! system. If the optional argument polgam is present, the other system
|
|
|
1591 |
! can also be a rotated one, where polgam is the angle between the two
|
|
|
1592 |
! north poles.
|
|
|
1593 |
! If polgam is not present, the other system is the real geographical
|
|
|
1594 |
! system.
|
|
|
1595 |
!
|
|
|
1596 |
! Method:
|
|
|
1597 |
! Transformation formulas for converting between these two systems.
|
|
|
1598 |
!
|
|
|
1599 |
! Modules used: NONE
|
|
|
1600 |
!
|
|
|
1601 |
!------------------------------------------------------------------------------
|
|
|
1602 |
|
|
|
1603 |
!------------------------------------------------------------------------------
|
|
|
1604 |
!
|
|
|
1605 |
! Declarations:
|
|
|
1606 |
!
|
|
|
1607 |
!------------------------------------------------------------------------------
|
|
|
1608 |
|
|
|
1609 |
! Parameter list:
|
|
|
1610 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
1611 |
polphi, & ! latitude of the rotated north pole
|
|
|
1612 |
pollam, & ! longitude of the rotated north pole
|
|
|
1613 |
phirot, & ! latitude in the rotated system
|
|
|
1614 |
rlarot ! longitude in the rotated system
|
|
|
1615 |
|
|
|
1616 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
1617 |
polgam ! angle between the north poles of the systems
|
|
|
1618 |
|
|
|
1619 |
REAL (KIND=ireals) :: &
|
|
|
1620 |
rlarot2rla ! latitude in the geographical system
|
|
|
1621 |
|
|
|
1622 |
! Local variables
|
|
|
1623 |
REAL (KIND=ireals) :: &
|
|
|
1624 |
zsinpol, zcospol, zlampol, zphis, zrlas, zarg1, zarg2, zgam
|
|
|
1625 |
|
|
|
1626 |
REAL (KIND=ireals), PARAMETER :: &
|
|
|
1627 |
zrpi18 = 57.2957795_ireals, & !
|
|
|
1628 |
zpir18 = 0.0174532925_ireals
|
|
|
1629 |
|
|
|
1630 |
!------------------------------------------------------------------------------
|
|
|
1631 |
|
|
|
1632 |
! Begin function rlarot2rla
|
|
|
1633 |
|
|
|
1634 |
zsinpol = SIN (zpir18 * polphi)
|
|
|
1635 |
zcospol = COS (zpir18 * polphi)
|
|
|
1636 |
|
|
|
1637 |
zlampol = zpir18 * pollam
|
|
|
1638 |
zphis = zpir18 * phirot
|
|
|
1639 |
IF (rlarot > 180.0_ireals) THEN
|
|
|
1640 |
zrlas = rlarot - 360.0_ireals
|
|
|
1641 |
ELSE
|
|
|
1642 |
zrlas = rlarot
|
|
|
1643 |
ENDIF
|
|
|
1644 |
zrlas = zpir18 * zrlas
|
|
|
1645 |
|
|
|
1646 |
IF (polgam /= 0.0_ireals) THEN
|
|
|
1647 |
zgam = zpir18 * polgam
|
|
|
1648 |
zarg1 = SIN (zlampol) * &
|
|
|
1649 |
(- zsinpol*COS(zphis) * (COS(zrlas)*COS(zgam) - SIN(zrlas)*SIN(zgam)) &
|
|
|
1650 |
+ zcospol * SIN(zphis)) &
|
|
|
1651 |
- COS (zlampol)*COS(zphis) * (SIN(zrlas)*COS(zgam) + COS(zrlas)*SIN(zgam))
|
|
|
1652 |
|
|
|
1653 |
zarg2 = COS (zlampol) * &
|
|
|
1654 |
(- zsinpol*COS(zphis) * (COS(zrlas)*COS(zgam) - SIN(zrlas)*SIN(zgam)) &
|
|
|
1655 |
+ zcospol * SIN(zphis)) &
|
|
|
1656 |
+ SIN (zlampol)*COS(zphis) * (SIN(zrlas)*COS(zgam) + COS(zrlas)*SIN(zgam))
|
|
|
1657 |
ELSE
|
|
|
1658 |
zarg1 = SIN (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + &
|
|
|
1659 |
zcospol * SIN(zphis)) - &
|
|
|
1660 |
COS (zlampol) * SIN(zrlas) * COS(zphis)
|
|
|
1661 |
zarg2 = COS (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + &
|
|
|
1662 |
zcospol * SIN(zphis)) + &
|
|
|
1663 |
SIN (zlampol) * SIN(zrlas) * COS(zphis)
|
|
|
1664 |
ENDIF
|
|
|
1665 |
|
|
|
1666 |
IF (zarg2 == 0.0) zarg2 = 1.0E-20_ireals
|
|
|
1667 |
|
|
|
1668 |
rlarot2rla = zrpi18 * ATAN2(zarg1,zarg2)
|
|
|
1669 |
|
|
|
1670 |
END FUNCTION rlarot2rla
|
|
|
1671 |
|
|
|
1672 |
!==============================================================================
|
|
|
1673 |
!==============================================================================
|
|
|
1674 |
|
|
|
1675 |
!------------------------------------------------------------------------------
|
|
|
1676 |
|
|
|
1677 |
FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam )
|
|
|
1678 |
|
|
|
1679 |
!------------------------------------------------------------------------------
|
|
|
1680 |
!
|
|
|
1681 |
! Description:
|
|
|
1682 |
! This routine converts lambda from the real geographical system to lambda
|
|
|
1683 |
! in the rotated system.
|
|
|
1684 |
!
|
|
|
1685 |
! Method:
|
|
|
1686 |
! Transformation formulas for converting between these two systems.
|
|
|
1687 |
!
|
|
|
1688 |
!------------------------------------------------------------------------------
|
|
|
1689 |
!
|
|
|
1690 |
! Parameter list:
|
|
|
1691 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
1692 |
polphi, & ! latitude of the rotated north pole
|
|
|
1693 |
pollam, & ! longitude of the rotated north pole
|
|
|
1694 |
phi, & ! latitude in geographical system
|
|
|
1695 |
rla ! longitude in geographical system
|
|
|
1696 |
|
|
|
1697 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
1698 |
polgam ! angle between the north poles of the systems
|
|
|
1699 |
|
|
|
1700 |
REAL (KIND=ireals) :: &
|
|
|
1701 |
rla2rlarot ! latitude in the the rotated system
|
|
|
1702 |
|
|
|
1703 |
! Local variables
|
|
|
1704 |
REAL (KIND=ireals) :: &
|
|
|
1705 |
zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
|
|
|
1706 |
|
|
|
1707 |
REAL (KIND=ireals), PARAMETER :: &
|
|
|
1708 |
zrpi18 = 57.2957795_ireals, & !
|
|
|
1709 |
zpir18 = 0.0174532925_ireals
|
|
|
1710 |
|
|
|
1711 |
!------------------------------------------------------------------------------
|
|
|
1712 |
|
|
|
1713 |
! Begin function rla2rlarot
|
|
|
1714 |
|
|
|
1715 |
zsinpol = SIN (zpir18 * polphi)
|
|
|
1716 |
zcospol = COS (zpir18 * polphi)
|
|
|
1717 |
zlampol = zpir18 * pollam
|
|
|
1718 |
zphi = zpir18 * phi
|
|
|
1719 |
IF (rla > 180.0_ireals) THEN
|
|
|
1720 |
zrla1 = rla - 360.0_ireals
|
|
|
1721 |
ELSE
|
|
|
1722 |
zrla1 = rla
|
|
|
1723 |
ENDIF
|
|
|
1724 |
zrla = zpir18 * zrla1
|
|
|
1725 |
|
|
|
1726 |
zarg1 = - SIN (zrla-zlampol) * COS(zphi)
|
|
|
1727 |
zarg2 = - zsinpol * COS(zphi) * COS(zrla-zlampol) + zcospol * SIN(zphi)
|
|
|
1728 |
|
|
|
1729 |
IF (zarg2 == 0.0) zarg2 = 1.0E-20_ireals
|
|
|
1730 |
|
|
|
1731 |
rla2rlarot = zrpi18 * ATAN2 (zarg1,zarg2)
|
|
|
1732 |
|
|
|
1733 |
IF (polgam /= 0.0_ireals ) THEN
|
|
|
1734 |
rla2rlarot = polgam + rla2rlarot
|
|
|
1735 |
IF (rla2rlarot > 180._ireals) rla2rlarot = rla2rlarot -360._ireals
|
|
|
1736 |
ENDIF
|
|
|
1737 |
|
|
|
1738 |
END FUNCTION rla2rlarot
|
|
|
1739 |
|
|
|
1740 |
!==============================================================================
|
|
|
1741 |
!==============================================================================
|
|
|
1742 |
|
|
|
1743 |
SUBROUTINE sleve_split_oro (hsurf, hsurfs, idim, jdim, nflt, nextralines, &
|
|
|
1744 |
svc1, svc2, vcflat, noutunit, myid, ierror, yerror)
|
|
|
1745 |
|
|
|
1746 |
!------------------------------------------------------------------------------
|
|
|
1747 |
!
|
|
|
1748 |
! Description:
|
|
|
1749 |
! decomposes a given topography field hsurf in a
|
|
|
1750 |
! large-scale (hsurfs(:,:,1)) and a small-scale (hsurfs(:,:,2)) part, where
|
|
|
1751 |
! hsurf(:,:) = hsurfs(:,:,1) + hsurfs(:,:,2)
|
|
|
1752 |
!
|
|
|
1753 |
! Method:
|
|
|
1754 |
! - a digital filter is applied for the computation of
|
|
|
1755 |
! the large scale part hsurfs(:,:,1).
|
|
|
1756 |
! - the boundary values are treated seperately to assure, that
|
|
|
1757 |
! also these points are smoothed:
|
|
|
1758 |
! i.e. at the i=1 boundary: A(1,j) = A(2,j) for all j
|
|
|
1759 |
! i=idim boundary: A(idim,j) = A(idim-1,j) for all j
|
|
|
1760 |
! j=1 boundary: A(i,1) = A(i,2) for all i
|
|
|
1761 |
! j=jdim boundary: A(i,jdim) = A(i,jdim-1) for all i
|
|
|
1762 |
! - nflt determines, how often the filter is applied
|
|
|
1763 |
! - Additionally, the maxima of hsurf, hsurfs(:,:,1) and hsurfs(:,:,2) are
|
|
|
1764 |
! computed and written to noutunit.
|
|
|
1765 |
!
|
|
|
1766 |
! written by Daniel Leuenberger, 03.10.2001
|
|
|
1767 |
!------------------------------------------------------------------------------
|
|
|
1768 |
|
|
|
1769 |
! Subroutine Arguments:
|
|
|
1770 |
|
|
|
1771 |
INTEGER (KIND=iintegers), INTENT(IN) :: &
|
|
|
1772 |
idim, jdim, & ! dimensions of hsurf
|
|
|
1773 |
nextralines, & ! number of extra lines around filtered
|
|
|
1774 |
! field (for interpolation program)
|
|
|
1775 |
nflt ! number of filter applications
|
|
|
1776 |
|
|
|
1777 |
REAL (KIND=ireals), INTENT(IN) :: &
|
|
|
1778 |
svc1, svc2, & ! decay rates for large and small scale
|
|
|
1779 |
vcflat ! vertical coordinate where the
|
|
|
1780 |
! terrain following system changes back
|
|
|
1781 |
! to an orthogonal z-system
|
|
|
1782 |
REAL (KIND=ireals), INTENT(IN) :: &
|
|
|
1783 |
hsurf(idim,jdim) ! height of full topography
|
|
|
1784 |
|
|
|
1785 |
INTEGER (KIND=iintegers), INTENT(IN) :: &
|
|
|
1786 |
noutunit ! unitnumber where output is written to
|
|
|
1787 |
|
|
|
1788 |
REAL (KIND=ireals), INTENT(OUT) :: &
|
|
|
1789 |
hsurfs(idim,jdim,2) ! height of splitted topography parts
|
|
|
1790 |
|
|
|
1791 |
INTEGER (KIND=iintegers), INTENT(IN) :: &
|
|
|
1792 |
myid ! PE number
|
|
|
1793 |
|
|
|
1794 |
INTEGER (KIND=iintegers), INTENT(OUT) :: &
|
|
|
1795 |
ierror ! error value
|
|
|
1796 |
|
|
|
1797 |
CHARACTER(LEN=*) :: &
|
|
|
1798 |
yerror ! error message
|
|
|
1799 |
|
|
|
1800 |
! Local variables
|
|
|
1801 |
REAL (KIND=ireals) :: &
|
|
|
1802 |
maxhsurf, & ! maximum of hsurf
|
|
|
1803 |
maxhsurf1, & ! maximum of hsurfs(:,:,1)
|
|
|
1804 |
maxhsurf2, & ! maximum of hsurfs(:,:,2)
|
|
|
1805 |
gammavc ! invertibility parameter
|
|
|
1806 |
|
|
|
1807 |
INTEGER (KIND=iintegers) :: &
|
|
|
1808 |
i,j,n,old,new,temp, istart, iend, jstart, jend
|
|
|
1809 |
|
|
|
1810 |
!------------------------------------------------------------------------------
|
|
|
1811 |
!- End of header -
|
|
|
1812 |
!------------------------------------------------------------------------------
|
|
|
1813 |
|
|
|
1814 |
!------------------------------------------------------------------------------
|
|
|
1815 |
!- Begin SUBROUTINE sleve_split_oro
|
|
|
1816 |
!------------------------------------------------------------------------------
|
|
|
1817 |
|
|
|
1818 |
ierror = 0_iintegers
|
|
|
1819 |
yerror = ' '
|
|
|
1820 |
|
|
|
1821 |
IF (myid == 0) THEN
|
|
|
1822 |
WRITE (noutunit,'(A)') ' '
|
|
|
1823 |
WRITE (noutunit,'(A)') ' Splitting of Topography for SLEVE coordinate:'
|
|
|
1824 |
ENDIF
|
|
|
1825 |
|
|
|
1826 |
IF ( (nextralines < 0) .OR. (nextralines > 1) ) THEN
|
|
|
1827 |
ierror = 1
|
|
|
1828 |
yerror = 'ERROR: nextralines outside range: 0 <= nextralines <= 1'
|
|
|
1829 |
RETURN
|
|
|
1830 |
ENDIF
|
|
|
1831 |
|
|
|
1832 |
! In order to obtain the same splitting in LM and in INT2LM, only the domain
|
|
|
1833 |
! without the extra boundary lines is splitted. These extra boundary lines
|
|
|
1834 |
! are indicated by nextralines
|
|
|
1835 |
|
|
|
1836 |
! compute index boundaries of LM-domain
|
|
|
1837 |
istart = 1 + nextralines
|
|
|
1838 |
iend = idim - nextralines
|
|
|
1839 |
jstart = 1 + nextralines
|
|
|
1840 |
jend = jdim - nextralines
|
|
|
1841 |
|
|
|
1842 |
maxhsurf = 0.0_ireals
|
|
|
1843 |
maxhsurf1 = 0.0_ireals
|
|
|
1844 |
maxhsurf2 = 0.0_ireals
|
|
|
1845 |
|
|
|
1846 |
old = 1
|
|
|
1847 |
new = 2
|
|
|
1848 |
|
|
|
1849 |
hsurfs(:,:,old) = hsurf(:,:)
|
|
|
1850 |
|
|
|
1851 |
! apply nflt times an ideal 2d filter to compute the
|
|
|
1852 |
! large-scale part hsurfs(:,:,1) of the topography
|
|
|
1853 |
|
|
|
1854 |
DO n = 1, nflt
|
|
|
1855 |
! treat inner points
|
|
|
1856 |
DO j=jstart+1, jend-1
|
|
|
1857 |
DO i=istart+1, iend-1
|
|
|
1858 |
hsurfs(i,j,new) = 0.25_ireals * hsurfs(i,j,old) &
|
|
|
1859 |
+ 0.125_ireals * (hsurfs(i-1,j ,old) + hsurfs(i+1,j ,old) + &
|
|
|
1860 |
hsurfs(i ,j-1,old) + hsurfs(i ,j+1,old)) &
|
|
|
1861 |
+ 0.0625_ireals * (hsurfs(i-1,j-1,old) + hsurfs(i+1,j-1,old) + &
|
|
|
1862 |
hsurfs(i-1,j+1,old) + hsurfs(i+1,j+1,old))
|
|
|
1863 |
ENDDO
|
|
|
1864 |
ENDDO
|
|
|
1865 |
|
|
|
1866 |
! treat corner points
|
|
|
1867 |
hsurfs(istart,jstart,new) = 0.25_ireals * &
|
|
|
1868 |
(hsurfs(istart, jstart ,old) + hsurfs(istart+1, jstart ,old) &
|
|
|
1869 |
+ hsurfs(istart, jstart+1,old) + hsurfs(istart+1, jstart+1,old))
|
|
|
1870 |
|
|
|
1871 |
hsurfs(istart, jend, new) = 0.25_ireals * &
|
|
|
1872 |
(hsurfs(istart, jend ,old) + hsurfs(istart+1, jend ,old) &
|
|
|
1873 |
+ hsurfs(istart, jend-1,old) + hsurfs(istart+1, jend-1,old))
|
|
|
1874 |
|
|
|
1875 |
hsurfs(iend, jstart, new) = 0.25_ireals * &
|
|
|
1876 |
(hsurfs(iend, jstart ,old) + hsurfs(iend-1,jstart ,old) &
|
|
|
1877 |
+ hsurfs(iend, jstart+1,old) + hsurfs(iend-1,jstart+1,old))
|
|
|
1878 |
|
|
|
1879 |
hsurfs(iend,jend,new) = 0.25_ireals * &
|
|
|
1880 |
(hsurfs(iend,jend ,old) + hsurfs(iend-1,jend ,old)&
|
|
|
1881 |
+ hsurfs(iend,jend-1,old) + hsurfs(iend-1,jend-1,old))
|
|
|
1882 |
|
|
|
1883 |
! treat edge points
|
|
|
1884 |
DO j = jstart+1,jend-1
|
|
|
1885 |
hsurfs(istart,j,new) = &
|
|
|
1886 |
0.25_ireals * (hsurfs(istart ,j ,old) + hsurfs(istart+1,j ,old)) &
|
|
|
1887 |
+ 0.125_ireals * (hsurfs(istart ,j-1,old) + hsurfs(istart ,j+1,old) +&
|
|
|
1888 |
hsurfs(istart+1,j-1,old) + hsurfs(istart+1,j+1,old) )
|
|
|
1889 |
|
|
|
1890 |
|
|
|
1891 |
hsurfs(iend,j,new) = &
|
|
|
1892 |
0.25_ireals * (hsurfs(iend ,j ,old) + hsurfs(iend-1,j ,old)) &
|
|
|
1893 |
+ 0.125_ireals * (hsurfs(iend ,j-1,old) + hsurfs(iend ,j+1,old) + &
|
|
|
1894 |
hsurfs(iend-1,j-1,old) + hsurfs(iend-1,j+1,old) )
|
|
|
1895 |
ENDDO
|
|
|
1896 |
|
|
|
1897 |
DO i = istart+1,iend-1
|
|
|
1898 |
hsurfs(i,jstart,new) = &
|
|
|
1899 |
0.25_ireals * (hsurfs(i ,jstart ,old) + hsurfs(i ,jstart+1,old)) &
|
|
|
1900 |
+ 0.125_ireals * (hsurfs(i-1,jstart ,old) + hsurfs(i+1,jstart ,old) +&
|
|
|
1901 |
hsurfs(i-1,jstart+1,old) + hsurfs(i+1,jstart+1,old) )
|
|
|
1902 |
|
|
|
1903 |
hsurfs(i,jend,new) = &
|
|
|
1904 |
0.25_ireals * (hsurfs(i ,jend ,old) + hsurfs(i ,jend-1,old)) &
|
|
|
1905 |
+ 0.125_ireals * (hsurfs(i-1,jend ,old) + hsurfs(i+1,jend ,old) + &
|
|
|
1906 |
hsurfs(i-1,jend-1,old) + hsurfs(i+1,jend-1,old) )
|
|
|
1907 |
ENDDO
|
|
|
1908 |
|
|
|
1909 |
temp = old
|
|
|
1910 |
old = new
|
|
|
1911 |
new = temp
|
|
|
1912 |
|
|
|
1913 |
ENDDO
|
|
|
1914 |
|
|
|
1915 |
! compute the large-scale part hsurfs(:,:,1) of the topo
|
|
|
1916 |
hsurfs(istart:iend,jstart:jend,1) = hsurfs(istart:iend,jstart:jend,old)
|
|
|
1917 |
|
|
|
1918 |
! compute the small-scale part hsurfs(:,:,2) of the topo
|
|
|
1919 |
hsurfs(istart:iend,jstart:jend,2) = hsurf (istart:iend,jstart:jend) - &
|
|
|
1920 |
hsurfs(istart:iend,jstart:jend,1)
|
|
|
1921 |
|
|
|
1922 |
! compute maxima of topographies
|
|
|
1923 |
maxhsurf = MAXVAL (hsurf (istart:iend,jstart:jend) )
|
|
|
1924 |
maxhsurf1 = MAXVAL (hsurfs(istart:iend,jstart:jend,1))
|
|
|
1925 |
maxhsurf2 = MAXVAL (hsurfs(istart:iend,jstart:jend,2))
|
|
|
1926 |
|
|
|
1927 |
IF (myid == 0) THEN
|
|
|
1928 |
WRITE(noutunit,'(A,I5,A)' ) ' nflt = ',nflt,' Applications of Filter'
|
|
|
1929 |
WRITE(noutunit,'(A)' ) ' Maxima of Topography Parts:'
|
|
|
1930 |
WRITE(noutunit,'(A,F10.5)') &
|
|
|
1931 |
' Max of Full Topography hsurf : ',maxhsurf
|
|
|
1932 |
WRITE(noutunit,'(A,F10.5)') &
|
|
|
1933 |
' Max of Large-Scale Topography hsurfs(:,:,1) : ',maxhsurf1
|
|
|
1934 |
WRITE(noutunit,'(A,F10.5)') &
|
|
|
1935 |
' Max of Small-Scale Topography hsurfs(:,:,2) : ',maxhsurf2
|
|
|
1936 |
ENDIF
|
|
|
1937 |
|
|
|
1938 |
! calculate SLEVE invertibility parameter gammavc
|
|
|
1939 |
gammavc = 1.0_ireals - &
|
|
|
1940 |
((MAXVAL(hsurfs(istart:iend,jstart:jend,1)) / svc1) / TANH(vcflat/svc1))- &
|
|
|
1941 |
((MAXVAL(hsurfs(istart:iend,jstart:jend,2)) / svc2) / TANH(vcflat/svc2))
|
|
|
1942 |
|
|
|
1943 |
IF (myid == 0) THEN
|
|
|
1944 |
WRITE (noutunit,'(A)') ' '
|
|
|
1945 |
WRITE (noutunit,'(A,F10.5)') &
|
|
|
1946 |
' Invertibility parameter for SLEVE coordinate: gammavc = ',gammavc
|
|
|
1947 |
WRITE (noutunit,'(A)') ' '
|
|
|
1948 |
ENDIF
|
|
|
1949 |
|
|
|
1950 |
! check if invertibility condition is fulfilled
|
|
|
1951 |
IF ( gammavc <= 0.0 ) THEN
|
|
|
1952 |
PRINT *, 'Invertibility parameter for SLEVE coordinate: gammavc = ',&
|
|
|
1953 |
gammavc
|
|
|
1954 |
PRINT *, 'vcflat = ',vcflat
|
|
|
1955 |
PRINT *, 'svc1 = ',svc1
|
|
|
1956 |
PRINT *, 'svc2 = ',svc2
|
|
|
1957 |
ierror = 2
|
|
|
1958 |
yerror = 'Invertibility condition of SLEVE coordinate not '// &
|
|
|
1959 |
'fulfilled, check values of svc1, svc2 and vcflat'
|
|
|
1960 |
RETURN
|
|
|
1961 |
ELSEIF ( gammavc < 0.05 ) THEN
|
|
|
1962 |
PRINT *, 'Invertibility parameter for SLEVE coordinate: gammavc = ',&
|
|
|
1963 |
gammavc
|
|
|
1964 |
PRINT *, 'vcflat = ',vcflat
|
|
|
1965 |
PRINT *, 'svc1 = ',svc1
|
|
|
1966 |
PRINT *, 'svc2 = ',svc2
|
|
|
1967 |
PRINT *, 'WARNING !!! SLEVE Invertibility parameter close to ', &
|
|
|
1968 |
'zero, check values of svc1, svc2 and vcflat'
|
|
|
1969 |
ENDIF
|
|
|
1970 |
|
|
|
1971 |
IF (nextralines > 0) THEN
|
|
|
1972 |
! The values of hsurfs outside the LM-domain are determined as follows:
|
|
|
1973 |
! First, the large-scale topo hsurfs(:,:,1) is linearly extrapolated
|
|
|
1974 |
! (from the point at the boundary of the LM-domain and the
|
|
|
1975 |
! first point inside the LM-domain),
|
|
|
1976 |
! then the small-scale topo hsurfs(:,:,2) is calculated from the
|
|
|
1977 |
! relationship h2 = h - h1
|
|
|
1978 |
|
|
|
1979 |
! ** Attention: This extrapolation works only in the case of
|
|
|
1980 |
! ** nextralines = 1 !!!
|
|
|
1981 |
! ** For nextralines > 1 the extrapolation has yet to be implemented !!!
|
|
|
1982 |
|
|
|
1983 |
DO i = istart, iend
|
|
|
1984 |
! extrapolation of south edge
|
|
|
1985 |
hsurfs(i,1,1) = 2 * hsurfs(i,2,1) - hsurfs(i,3,1)
|
|
|
1986 |
hsurfs(i,1,2) = hsurf (i,1) - hsurfs(i,1,1)
|
|
|
1987 |
|
|
|
1988 |
! extrapolation of north edge
|
|
|
1989 |
hsurfs(i,jdim,1) = 2 * hsurfs(i,jdim-1,1) - hsurfs(i,jdim-2,1)
|
|
|
1990 |
hsurfs(i,jdim,2) = hsurf (i,jdim) - hsurfs(i,jdim ,1)
|
|
|
1991 |
ENDDO
|
|
|
1992 |
|
|
|
1993 |
DO j = jstart, jend
|
|
|
1994 |
! extrapolation of west edge
|
|
|
1995 |
hsurfs(1,j,1) = 2 * hsurfs(2,j,1) - hsurfs(3,j,1)
|
|
|
1996 |
hsurfs(1,j,2) = hsurf (1,j) - hsurfs(1,j,1)
|
|
|
1997 |
|
|
|
1998 |
! extrapolation of east edge
|
|
|
1999 |
hsurfs(idim,j,1) = 2 * hsurfs(idim-1,j,1) - hsurfs(idim-2,j,1)
|
|
|
2000 |
hsurfs(idim,j,2) = hsurf (idim ,j) - hsurfs(idim ,j,1)
|
|
|
2001 |
ENDDO
|
|
|
2002 |
|
|
|
2003 |
! extrapolation of SW point
|
|
|
2004 |
hsurfs(1,1,1) = 2 * hsurfs(2,2,1) - hsurfs(3,3,1)
|
|
|
2005 |
hsurfs(1,1,2) = hsurf(1,1) - hsurfs(1,1,1)
|
|
|
2006 |
|
|
|
2007 |
! extrapolation of SE point
|
|
|
2008 |
hsurfs(idim,1,1) = 2 * hsurfs(idim-1,2,1) - hsurfs(idim-2,3,1)
|
|
|
2009 |
hsurfs(idim,1,2) = hsurf (idim ,1) - hsurfs(idim ,1,1)
|
|
|
2010 |
|
|
|
2011 |
! extrapolation of NW point
|
|
|
2012 |
hsurfs(1,jdim,1) = 2 * hsurfs(2,jdim-1,1) - hsurfs(3,jdim-2,1)
|
|
|
2013 |
hsurfs(1,jdim,2) = hsurf (1,jdim ) - hsurfs(1,jdim ,1)
|
|
|
2014 |
|
|
|
2015 |
! extrapolation of NE point
|
|
|
2016 |
hsurfs(idim,jdim,1) = 2*hsurfs(idim-1,jdim-1,1)-hsurfs(idim-2,jdim-2,1)
|
|
|
2017 |
hsurfs(idim,jdim,2) = hsurf (idim ,jdim) -hsurfs(idim ,jdim,1)
|
|
|
2018 |
ENDIF
|
|
|
2019 |
|
|
|
2020 |
!------------------------------------------------------------------------------
|
|
|
2021 |
! End of the Subroutine
|
|
|
2022 |
!------------------------------------------------------------------------------
|
|
|
2023 |
|
|
|
2024 |
END SUBROUTINE sleve_split_oro
|
|
|
2025 |
|
|
|
2026 |
!==============================================================================
|
|
|
2027 |
!==============================================================================
|
|
|
2028 |
!+ Defines all subroutines for the generic routine smoother
|
|
|
2029 |
!------------------------------------------------------------------------------
|
|
|
2030 |
!
|
|
|
2031 |
! SUBROUTINE smoother (finout, ie, je, nlength, nfilt)
|
|
|
2032 |
!
|
|
|
2033 |
!------------------------------------------------------------------------------
|
|
|
2034 |
!
|
|
|
2035 |
! Description:
|
|
|
2036 |
! This routine smoothes an arbitrary two-dimensional field (fin) by applying
|
|
|
2037 |
! digital filters of length nlength (4 or 8) nfilt times. The filterd field
|
|
|
2038 |
! is written on fout.
|
|
|
2039 |
!
|
|
|
2040 |
! Method:
|
|
|
2041 |
! Call of digital filters (dfilt4 or dfilt8) in each direction.
|
|
|
2042 |
!
|
|
|
2043 |
!------------------------------------------------------------------------------
|
|
|
2044 |
!+ Subroutine for double precision
|
|
|
2045 |
!------------------------------------------------------------------------------
|
|
|
2046 |
|
|
|
2047 |
SUBROUTINE smoother_double (finout, ie, je, nlength, nfilt)
|
|
|
2048 |
|
|
|
2049 |
!------------------------------------------------------------------------------
|
|
|
2050 |
!
|
|
|
2051 |
! Parameter list:
|
|
|
2052 |
INTEGER (KIND=iintegers), INTENT (IN) :: &
|
|
|
2053 |
ie, je, & ! Dimension of the field
|
|
|
2054 |
nlength, & ! Filter lenght
|
|
|
2055 |
nfilt ! Number of iterative filerings
|
|
|
2056 |
|
|
|
2057 |
REAL (KIND=idouble), INTENT (INOUT) :: &
|
|
|
2058 |
finout (ie*je) ! 2-d field: unfiltered at input, filtered at output
|
|
|
2059 |
|
|
|
2060 |
! Local variables
|
|
|
2061 |
INTEGER (KIND=iintegers) :: &
|
|
|
2062 |
i,j ! loop indicees
|
|
|
2063 |
|
|
|
2064 |
REAL (KIND=idouble) :: &
|
|
|
2065 |
f_2d_field(ie,je), & !
|
|
|
2066 |
sxin(ie), sxh(ie), sxout(ie), & ! local storage
|
|
|
2067 |
syin(je), syh(je), syout(je) ! local storage
|
|
|
2068 |
|
|
|
2069 |
!------------------------------------------------------------------------------
|
|
|
2070 |
! begin subroutine smoother_double
|
|
|
2071 |
|
|
|
2072 |
f_2d_field = RESHAPE (finout, (/ie,je/))
|
|
|
2073 |
|
|
|
2074 |
IF ( nlength /= 4 .AND. nlength /= 8 ) THEN
|
|
|
2075 |
PRINT*, ' CAUTION: Filterlength =',nlength,' not implemented'
|
|
|
2076 |
PRINT*, ' No filtering of output field done'
|
|
|
2077 |
RETURN
|
|
|
2078 |
ENDIF
|
|
|
2079 |
|
|
|
2080 |
DO j = 1, je
|
|
|
2081 |
sxin(:) = f_2d_field(:,j)
|
|
|
2082 |
IF(nlength==4) CALL dfilt4 ( sxin, ie, sxh, sxout, nfilt )
|
|
|
2083 |
IF(nlength==8) CALL dfilt8 ( sxin, ie, sxh, sxout, nfilt )
|
|
|
2084 |
f_2d_field(:,j) = sxout(:)
|
|
|
2085 |
ENDDO
|
|
|
2086 |
DO i = 1, ie
|
|
|
2087 |
syin(:) = f_2d_field(i,:)
|
|
|
2088 |
IF(nlength==4) CALL dfilt4 ( syin, je, syh, syout, nfilt )
|
|
|
2089 |
IF(nlength==8) CALL dfilt8 ( syin, je, syh, syout, nfilt )
|
|
|
2090 |
f_2d_field(i,:) = syout(:)
|
|
|
2091 |
ENDDO
|
|
|
2092 |
|
|
|
2093 |
finout = RESHAPE (f_2d_field, (/ie*je/))
|
|
|
2094 |
|
|
|
2095 |
END SUBROUTINE smoother_double
|
|
|
2096 |
|
|
|
2097 |
!------------------------------------------------------------------------------
|
|
|
2098 |
!+ Subroutine for single precision
|
|
|
2099 |
!------------------------------------------------------------------------------
|
|
|
2100 |
|
|
|
2101 |
SUBROUTINE smoother_single (finout, ie, je, nlength, nfilt)
|
|
|
2102 |
|
|
|
2103 |
!------------------------------------------------------------------------------
|
|
|
2104 |
!
|
|
|
2105 |
! Parameter list:
|
|
|
2106 |
INTEGER (KIND=iintegers), INTENT (IN) :: &
|
|
|
2107 |
ie, je, & ! Dimension of the field
|
|
|
2108 |
nlength, & ! Filter lenght
|
|
|
2109 |
nfilt ! Number of iterative filerings
|
|
|
2110 |
|
|
|
2111 |
REAL (KIND=isingle), INTENT (INOUT) :: &
|
|
|
2112 |
finout (ie*je) ! 2-d field: unfiltered at input, filtered at output
|
|
|
2113 |
|
|
|
2114 |
! Local variables
|
|
|
2115 |
INTEGER (KIND=iintegers) :: &
|
|
|
2116 |
i,j ! loop indicees
|
|
|
2117 |
|
|
|
2118 |
REAL (KIND=isingle) :: &
|
|
|
2119 |
f_2d_field(ie,je), & !
|
|
|
2120 |
sxin(ie), sxh(ie), sxout(ie), & ! local storage
|
|
|
2121 |
syin(je), syh(je), syout(je) ! local storage
|
|
|
2122 |
|
|
|
2123 |
!------------------------------------------------------------------------------
|
|
|
2124 |
! begin subroutine smoother_single
|
|
|
2125 |
|
|
|
2126 |
f_2d_field = RESHAPE (finout, (/ie,je/))
|
|
|
2127 |
|
|
|
2128 |
IF ( nlength /= 4 .AND. nlength /= 8 ) THEN
|
|
|
2129 |
PRINT*, ' CAUTION: Filterlength =',nlength,' not implemented'
|
|
|
2130 |
PRINT*, ' No filtering of output field done'
|
|
|
2131 |
RETURN
|
|
|
2132 |
ENDIF
|
|
|
2133 |
|
|
|
2134 |
DO j = 1, je
|
|
|
2135 |
sxin(:) = f_2d_field(:,j)
|
|
|
2136 |
IF(nlength==4) CALL dfilt4 ( sxin, ie, sxh, sxout, nfilt )
|
|
|
2137 |
IF(nlength==8) CALL dfilt8 ( sxin, ie, sxh, sxout, nfilt )
|
|
|
2138 |
f_2d_field(:,j) = sxout(:)
|
|
|
2139 |
ENDDO
|
|
|
2140 |
DO i = 1, ie
|
|
|
2141 |
syin(:) = f_2d_field(i,:)
|
|
|
2142 |
IF(nlength==4) CALL dfilt4 ( syin, je, syh, syout, nfilt )
|
|
|
2143 |
IF(nlength==8) CALL dfilt8 ( syin, je, syh, syout, nfilt )
|
|
|
2144 |
f_2d_field(i,:) = syout(:)
|
|
|
2145 |
ENDDO
|
|
|
2146 |
|
|
|
2147 |
finout = RESHAPE (f_2d_field, (/ie*je/))
|
|
|
2148 |
|
|
|
2149 |
END SUBROUTINE smoother_single
|
|
|
2150 |
|
|
|
2151 |
!==============================================================================
|
|
|
2152 |
!==============================================================================
|
|
|
2153 |
|
|
|
2154 |
!------------------------------------------------------------------------------
|
|
|
2155 |
|
|
|
2156 |
SUBROUTINE smooth9 (finout, imin,imaxx,jmin,jmaxx,ie,je,ke)
|
|
|
2157 |
|
|
|
2158 |
!------------------------------------------------------------------------------
|
|
|
2159 |
!
|
|
|
2160 |
! Description:
|
|
|
2161 |
! This routine smoothes an arbitrary two-dimensional field (finout) by applying
|
|
|
2162 |
! a 9 points smoother. The filtered field is written on finout.
|
|
|
2163 |
!
|
|
|
2164 |
! Method:
|
|
|
2165 |
! A 9 points smoother is applied for the computation of the
|
|
|
2166 |
! large scale part of finout. The boundary values are treated
|
|
|
2167 |
! separately to assure, that also these points are smoothed.
|
|
|
2168 |
!
|
|
|
2169 |
!------------------------------------------------------------------------------
|
|
|
2170 |
|
|
|
2171 |
! Parameter list:
|
|
|
2172 |
INTEGER (KIND=iintegers), INTENT (IN) :: &
|
|
|
2173 |
imin,jmin,imaxx,jmaxx, & ! Local Dimension of the field
|
|
|
2174 |
ie, je, ke ! Dimension of the field
|
|
|
2175 |
|
|
|
2176 |
REAL (KIND=ireals), INTENT (INOUT) :: &
|
|
|
2177 |
finout (ie,je,ke) ! 3-d field: unfiltered at input, filtered at output
|
|
|
2178 |
|
|
|
2179 |
! Local variables
|
|
|
2180 |
INTEGER (KIND=iintegers) :: &
|
|
|
2181 |
i,j,k, &! loop indicees
|
|
|
2182 |
imin1, imaxx1,jmin1, jmaxx1
|
|
|
2183 |
|
|
|
2184 |
LOGICAL lbord12,lbord13,lbord24,lbord34,lcorn1,lcorn2,lcorn3,lcorn4
|
|
|
2185 |
|
|
|
2186 |
REAL (KIND=ireals) :: fhelp (ie,je) ! local storage
|
|
|
2187 |
|
|
|
2188 |
!------------------------------------------------------------------------------
|
|
|
2189 |
! begin subroutine smooth9
|
|
|
2190 |
!
|
|
|
2191 |
lcorn1=.false.
|
|
|
2192 |
lcorn2=.false.
|
|
|
2193 |
lcorn3=.false.
|
|
|
2194 |
lcorn4=.false.
|
|
|
2195 |
lbord13=.false.
|
|
|
2196 |
lbord12=.false.
|
|
|
2197 |
lbord34=.false.
|
|
|
2198 |
lbord24=.false.
|
|
|
2199 |
|
|
|
2200 |
! If applied to a subdomain, smooth also points of a grid line halo
|
|
|
2201 |
|
|
|
2202 |
imin1=imin
|
|
|
2203 |
jmin1=jmin
|
|
|
2204 |
imaxx1=imaxx
|
|
|
2205 |
jmaxx1=jmaxx
|
|
|
2206 |
|
|
|
2207 |
! Adapt dimensions to the subdomain type (if it has edge and/or corner points)
|
|
|
2208 |
|
|
|
2209 |
IF( imin == 1 ) THEN
|
|
|
2210 |
imin1=imin+1
|
|
|
2211 |
lbord12=.true.
|
|
|
2212 |
IF(jmin == 1) THEN
|
|
|
2213 |
lcorn1=.true.
|
|
|
2214 |
ENDIF
|
|
|
2215 |
IF(jmaxx == je) THEN
|
|
|
2216 |
lcorn2=.true.
|
|
|
2217 |
ENDIF
|
|
|
2218 |
ENDIF
|
|
|
2219 |
IF( jmin == 1 ) THEN
|
|
|
2220 |
jmin1=jmin+1
|
|
|
2221 |
lbord13=.true.
|
|
|
2222 |
ENDIF
|
|
|
2223 |
IF(imaxx == ie )THEN
|
|
|
2224 |
imaxx1=imaxx-1
|
|
|
2225 |
lbord34=.true.
|
|
|
2226 |
IF (jmin == 1) THEN
|
|
|
2227 |
lcorn3=.true.
|
|
|
2228 |
ENDIF
|
|
|
2229 |
IF(jmaxx == je) THEN
|
|
|
2230 |
lcorn4=.true.
|
|
|
2231 |
ENDIF
|
|
|
2232 |
ENDIF
|
|
|
2233 |
IF(jmaxx == je) THEN
|
|
|
2234 |
jmaxx1=jmaxx-1
|
|
|
2235 |
lbord24=.true.
|
|
|
2236 |
ENDIF
|
|
|
2237 |
|
|
|
2238 |
DO k = 1, ke
|
|
|
2239 |
|
|
|
2240 |
fhelp(:,:) = finout(:,:,k)
|
|
|
2241 |
|
|
|
2242 |
! Treat inner points
|
|
|
2243 |
|
|
|
2244 |
DO i=imin1,imaxx1
|
|
|
2245 |
DO j=jmin1,jmaxx1
|
|
|
2246 |
finout(i,j,k) = 0.25_ireals * fhelp(i,j) &
|
|
|
2247 |
+ 0.125_ireals * (fhelp(i-1,j ) + fhelp(i+1,j ) + &
|
|
|
2248 |
fhelp(i ,j-1) + fhelp(i ,j+1)) &
|
|
|
2249 |
+ 0.0625_ireals * (fhelp(i-1,j-1) + fhelp(i+1,j-1) + &
|
|
|
2250 |
fhelp(i-1,j+1) + fhelp(i+1,j+1))
|
|
|
2251 |
ENDDO
|
|
|
2252 |
ENDDO
|
|
|
2253 |
|
|
|
2254 |
! Treat corner points
|
|
|
2255 |
|
|
|
2256 |
IF (lcorn1) finout(1,1,k) = 0.25_ireals * &
|
|
|
2257 |
(fhelp(1 ,1 ) + fhelp(2 ,1 ) &
|
|
|
2258 |
+ fhelp(1 ,2 ) + fhelp(2 ,2 ))
|
|
|
2259 |
|
|
|
2260 |
IF (lcorn2) finout(1 ,jmaxx,k) = 0.25_ireals * &
|
|
|
2261 |
(fhelp(1 ,jmaxx ) + fhelp( 2,jmaxx ) &
|
|
|
2262 |
+ fhelp(1 ,jmaxx-1) + fhelp( 2,jmaxx-1))
|
|
|
2263 |
|
|
|
2264 |
IF (lcorn3) finout(imaxx,1 ,k) = 0.25_ireals * &
|
|
|
2265 |
(fhelp(imaxx,1 ) + fhelp(imaxx-1,1 ) &
|
|
|
2266 |
+ fhelp(imaxx,2 ) + fhelp(imaxx-1,2 ))
|
|
|
2267 |
|
|
|
2268 |
IF (lcorn4) finout(imaxx,jmaxx,k) = 0.25_ireals * &
|
|
|
2269 |
(fhelp(imaxx,jmaxx ) + fhelp(imaxx-1,jmaxx )&
|
|
|
2270 |
+ fhelp(imaxx,jmaxx-1) + fhelp(imaxx-1,jmaxx-1))
|
|
|
2271 |
! Treat edge points
|
|
|
2272 |
|
|
|
2273 |
IF(lbord12)THEN
|
|
|
2274 |
DO j = jmin1,jmaxx1
|
|
|
2275 |
finout(1,j,k) = &
|
|
|
2276 |
0.25_ireals * (fhelp(1 ,j ) + fhelp(2 ,j )) &
|
|
|
2277 |
+ 0.125_ireals * (fhelp(1 ,j -1) + fhelp(1 ,j +1) + &
|
|
|
2278 |
fhelp(2 ,j -1) + fhelp(2 ,j +1) )
|
|
|
2279 |
ENDDO
|
|
|
2280 |
ENDIF
|
|
|
2281 |
|
|
|
2282 |
IF(lbord34)THEN
|
|
|
2283 |
DO j = jmin1,jmaxx1
|
|
|
2284 |
finout(imaxx,j,k) = &
|
|
|
2285 |
0.25_ireals * (fhelp(imaxx ,j ) + fhelp(imaxx-1,j )) &
|
|
|
2286 |
+ 0.125_ireals * (fhelp(imaxx ,j -1) + fhelp(imaxx ,j +1) + &
|
|
|
2287 |
fhelp(imaxx-1,j -1) + fhelp(imaxx-1,j +1) )
|
|
|
2288 |
ENDDO
|
|
|
2289 |
ENDIF
|
|
|
2290 |
|
|
|
2291 |
IF(lbord13)THEN
|
|
|
2292 |
DO i = imin1,imaxx1
|
|
|
2293 |
finout(i,1,k) = &
|
|
|
2294 |
0.25_ireals * (fhelp(i ,1 ) + fhelp(i ,2 )) &
|
|
|
2295 |
+ 0.125_ireals * (fhelp(i -1,1 ) + fhelp(i +1,1 ) + &
|
|
|
2296 |
fhelp(i -1,2 ) + fhelp(i +1,2 ) )
|
|
|
2297 |
ENDDO
|
|
|
2298 |
ENDIF
|
|
|
2299 |
|
|
|
2300 |
IF(lbord24)THEN
|
|
|
2301 |
DO i = imin1,imaxx1
|
|
|
2302 |
finout(i,jmaxx,k) = &
|
|
|
2303 |
0.25_ireals * (fhelp(i ,jmaxx ) + fhelp(i ,jmaxx-1)) &
|
|
|
2304 |
+ 0.125_ireals * (fhelp(i -1,jmaxx ) + fhelp(i +1,jmaxx ) + &
|
|
|
2305 |
fhelp(i -1,jmaxx-1) + fhelp(i +1,jmaxx-1) )
|
|
|
2306 |
ENDDO
|
|
|
2307 |
ENDIF
|
|
|
2308 |
|
|
|
2309 |
ENDDO
|
|
|
2310 |
|
|
|
2311 |
|
|
|
2312 |
END SUBROUTINE smooth9
|
|
|
2313 |
|
|
|
2314 |
!==============================================================================
|
|
|
2315 |
!==============================================================================
|
|
|
2316 |
|
|
|
2317 |
!------------------------------------------------------------------------------
|
|
|
2318 |
|
|
|
2319 |
SUBROUTINE TAUTSP(TAU,GTAU,NTAU,GAMMA,S,BREAK,COEF,L,IFLAG)
|
|
|
2320 |
|
|
|
2321 |
!------------------------------------------------------------------------------
|
|
|
2322 |
!
|
|
|
2323 |
! Description:
|
|
|
2324 |
! BERECHNET DEN TAUT-SPLINE FUER DIE DATEN: TAU(I),GTAU(I),I=1,.,NTAU
|
|
|
2325 |
!
|
|
|
2326 |
! Method:
|
|
|
2327 |
! WENN GAMMA GT.0 WERDEN ZUSAETZLICHE KNOTEN BERECHNET
|
|
|
2328 |
! GAMMA I.A. 2.5 BZW. 5.5
|
|
|
2329 |
!
|
|
|
2330 |
! BREAK,COEF,L,K GEBEN DIE PP-DARSTELLUNG DER INTERPOLATIONSFUNKTION
|
|
|
2331 |
!
|
|
|
2332 |
! FUER BREAK(I).LE.X.LE.BREAK(I+1) HAT DIE INTERPOLATIONSFUNKTION
|
|
|
2333 |
! FOLGENDE FORM
|
|
|
2334 |
!
|
|
|
2335 |
! F(X)=COEF(1,I)+DX(COEF(2,I)+DX/2(COEF(3,I)+DX/3(COEF(4,I)))
|
|
|
2336 |
! MIT DX=X-BREAK(I) UND I=1,...,L
|
|
|
2337 |
!
|
|
|
2338 |
! IFLAG=0 KEIN FEHLER IN TAUTSP
|
|
|
2339 |
! IFLAG=2 FALSCHER INPUT
|
|
|
2340 |
!
|
|
|
2341 |
! S(NTAU,6) WORK-ARRAY
|
|
|
2342 |
!
|
|
|
2343 |
!==============================================================================
|
|
|
2344 |
|
|
|
2345 |
INTEGER (KIND=iintegers) IFLAG,L,NTAU,I,METHOD,NTAUM1
|
|
|
2346 |
REAL (KIND=ireals) &
|
|
|
2347 |
BREAK(L),COEF(4,L),GAMMA,GTAU(NTAU),S(NTAU,6),TAU(NTAU), &
|
|
|
2348 |
ALPHA,C,D,DEL,DENOM,DIVDIF,ENTRY,ENTRY3,FACTOR,FACTR2,GAM, &
|
|
|
2349 |
ONEMG3,ONEMZT,RATIO,SIXTH,TEMP,X,Z,ZETA,ZT2,ALPH
|
|
|
2350 |
|
|
|
2351 |
!==============================================================================
|
|
|
2352 |
!
|
|
|
2353 |
ALPH(X)= MIN (1.0_ireals,ONEMG3/X)
|
|
|
2354 |
!
|
|
|
2355 |
! TEST DER INPUTPARAMETER
|
|
|
2356 |
!
|
|
|
2357 |
IF (NTAU .LT. 4) THEN
|
|
|
2358 |
PRINT 600,NTAU
|
|
|
2359 |
600 FORMAT(' NTAU =',I4,' NTAU SOLL GROESSER ALS 4 SEIN')
|
|
|
2360 |
GO TO 999
|
|
|
2361 |
ENDIF
|
|
|
2362 |
!
|
|
|
2363 |
! BERECHNUNG VON DELTA TAU UND DER 1. UND 2.ABLEITUNG DER DATEN
|
|
|
2364 |
!
|
|
|
2365 |
NTAUM1=NTAU-1
|
|
|
2366 |
DO I=1,NTAUM1
|
|
|
2367 |
S(I,1)=TAU(I+1)-TAU(I)
|
|
|
2368 |
IF (S(I,1) .LE. 0.0) THEN
|
|
|
2369 |
PRINT 610,I,TAU(I),TAU(I+1)
|
|
|
2370 |
610 FORMAT(' PUNKT ',I3,' UND DIE NAECHSTEN',2E15.6,'SIND IN&
|
|
|
2371 |
& FALSCHER REIHENFOLGE')
|
|
|
2372 |
GO TO 999
|
|
|
2373 |
ENDIF
|
|
|
2374 |
S(I+1,4) = (GTAU(I+1) - GTAU(I))/S(I,1)
|
|
|
2375 |
ENDDO
|
|
|
2376 |
!
|
|
|
2377 |
DO I=2,NTAUM1
|
|
|
2378 |
S(I,4) = S(I+1,4) - S(I,4)
|
|
|
2379 |
ENDDO
|
|
|
2380 |
!
|
|
|
2381 |
! 2.ABLEITUNG VON GTAU AN DEN PUNKTEN TAU
|
|
|
2382 |
!
|
|
|
2383 |
I=2
|
|
|
2384 |
S(2,2) = S(1,1)/3.0
|
|
|
2385 |
SIXTH = 1.0/6.0
|
|
|
2386 |
METHOD = 2
|
|
|
2387 |
GAM = GAMMA
|
|
|
2388 |
IF(GAM .LE. 0.0) METHOD = 1
|
|
|
2389 |
IF(GAM .GT. 3.0) THEN
|
|
|
2390 |
METHOD = 3
|
|
|
2391 |
GAM = GAM - 3.0
|
|
|
2392 |
ENDIF
|
|
|
2393 |
ONEMG3=1.0 - GAM/3.0
|
|
|
2394 |
!
|
|
|
2395 |
! SCHLEIFE UEBER I
|
|
|
2396 |
!
|
|
|
2397 |
10 CONTINUE
|
|
|
2398 |
Z=0.5
|
|
|
2399 |
IF (METHOD .EQ. 1) THEN
|
|
|
2400 |
GO TO 19
|
|
|
2401 |
ELSE IF (METHOD .EQ. 2) THEN
|
|
|
2402 |
GO TO 11
|
|
|
2403 |
ELSE IF (METHOD .EQ. 3) THEN
|
|
|
2404 |
GO TO 12
|
|
|
2405 |
ENDIF
|
|
|
2406 |
11 CONTINUE
|
|
|
2407 |
IF(S(I,4)*S(I+1,4).LT.0.0) GO TO 19
|
|
|
2408 |
12 CONTINUE
|
|
|
2409 |
TEMP = ABS(S(I+1,4))
|
|
|
2410 |
DENOM = ABS(S(I,4)) + TEMP
|
|
|
2411 |
IF(DENOM.EQ.0.0) GO TO 19
|
|
|
2412 |
Z = TEMP/DENOM
|
|
|
2413 |
IF(ABS(Z-0.5).LE.SIXTH) Z=0.5
|
|
|
2414 |
19 CONTINUE
|
|
|
2415 |
S(I,5) = Z
|
|
|
2416 |
!
|
|
|
2417 |
! ERSTELLEN EINES TEILES DER I-TEN GLEICHUNG
|
|
|
2418 |
!
|
|
|
2419 |
IF (Z-0.5 .LT. 0.) THEN
|
|
|
2420 |
ZETA = GAM*Z
|
|
|
2421 |
ONEMZT = 1.0 - ZETA
|
|
|
2422 |
ZT2 = ZETA**2
|
|
|
2423 |
ALPHA = ALPH(ONEMZT)
|
|
|
2424 |
FACTOR = ZETA/(ALPHA*(ZT2 - 1.0) + 1.0)
|
|
|
2425 |
S(I,6) = ZETA*FACTOR/6.0
|
|
|
2426 |
S(I,2) = S(I,2) + S(I,1)*((1.0-ALPHA*ONEMZT)*FACTOR*0.5-S(I,6))
|
|
|
2427 |
IF(S(I,2).LE.0.0) S(I,2) = 1.0
|
|
|
2428 |
S(I,3) = S(I,1)/6.0
|
|
|
2429 |
!
|
|
|
2430 |
ELSE IF (Z-0.5 .EQ. 0.) THEN
|
|
|
2431 |
!
|
|
|
2432 |
S(I,2) = S(I,2) + S(I,1)/3.0
|
|
|
2433 |
S(I,3) = S(I,1)/6.0
|
|
|
2434 |
!
|
|
|
2435 |
ELSE
|
|
|
2436 |
!
|
|
|
2437 |
ONEMZT = GAM*(1.0 - Z)
|
|
|
2438 |
ZETA = 1.0 - ONEMZT
|
|
|
2439 |
ALPHA = ALPH(ZETA)
|
|
|
2440 |
FACTOR = ONEMZT/(1.0 - ALPHA*ZETA*(1.0 + ONEMZT))
|
|
|
2441 |
S(I,6) = ONEMZT*FACTOR/6.0
|
|
|
2442 |
S(I,2) = S(I,2) + S(I,1)/3.0
|
|
|
2443 |
S(I,3) = S(I,6) * S(I,1)
|
|
|
2444 |
ENDIF
|
|
|
2445 |
!
|
|
|
2446 |
IF (I .GT. 2) GO TO 30
|
|
|
2447 |
S(1,5) = 0.5
|
|
|
2448 |
!
|
|
|
2449 |
! DIE ERSTEN BEIDEN GLEICHUNGEN ERZWINGEN STETIGKEIT DER 1. UND 3.AB-
|
|
|
2450 |
! LEITUNG IN TAU(I)
|
|
|
2451 |
!
|
|
|
2452 |
S(1,2) = S(1,1)/6.0
|
|
|
2453 |
S(1,3) = S(2,2)
|
|
|
2454 |
ENTRY3 = S(2,3)
|
|
|
2455 |
|
|
|
2456 |
IF (Z-0.5 .LT. 0.) THEN
|
|
|
2457 |
FACTR2 = ZETA*(ALPHA*(ZT2-1.0)+1.0)/(ALPHA*(ZETA*ZT2-1.0) + 1.0)
|
|
|
2458 |
RATIO = FACTR2*S(2,1)/S(1,2)
|
|
|
2459 |
S(2,2) = FACTR2*S(2,1) + S(1,1)
|
|
|
2460 |
S(2,3) = -FACTR2 * S(1,1)
|
|
|
2461 |
!
|
|
|
2462 |
ELSE IF (Z-0.5 .EQ. 0.) THEN
|
|
|
2463 |
!
|
|
|
2464 |
RATIO = S(2,1)/S(1,2)
|
|
|
2465 |
S(2,2) = S(2,1) + S(1,1)
|
|
|
2466 |
S(2,3) = -S(1,1)
|
|
|
2467 |
!
|
|
|
2468 |
ELSE
|
|
|
2469 |
!
|
|
|
2470 |
RATIO = S(2,1)/S(1,2)
|
|
|
2471 |
S(2,2) = S(2,1) + S(1,1)
|
|
|
2472 |
S(2,3) = -S(1,1)*6.0*ALPHA*S(2,6)
|
|
|
2473 |
ENDIF
|
|
|
2474 |
!
|
|
|
2475 |
! ELIMINATION DER 1.UNBEKANNTEN AUS DER 2.GLEICHUNG
|
|
|
2476 |
!
|
|
|
2477 |
S(2,2) = RATIO*S(1,3) + S(2,2)
|
|
|
2478 |
S(2,3) = RATIO*ENTRY3 + S(2,3)
|
|
|
2479 |
S(1,4) = S(2,4)
|
|
|
2480 |
S(2,4) = RATIO*S(1,4)
|
|
|
2481 |
GO TO 35
|
|
|
2482 |
!
|
|
|
2483 |
!
|
|
|
2484 |
30 CONTINUE
|
|
|
2485 |
S(I,2) = RATIO*S(I-1,3) + S(I,2)
|
|
|
2486 |
S(I,4) = RATIO*S(I-1,4) + S(I,4)
|
|
|
2487 |
!
|
|
|
2488 |
! AUFSTELLEN DES TEILES DER NAECHSTEN GLEICHUNG,DER VOM I-TEN INTERVAL
|
|
|
2489 |
! ABHAENGT
|
|
|
2490 |
!
|
|
|
2491 |
35 CONTINUE
|
|
|
2492 |
IF (Z-0.5 .LT. 0.) THEN
|
|
|
2493 |
RATIO = -S(I,6)*S(I,1)/S(I,2)
|
|
|
2494 |
S(I+1,2) = S(I,1)/3.0
|
|
|
2495 |
!
|
|
|
2496 |
ELSE IF (Z-0.5 .EQ. 0.) THEN
|
|
|
2497 |
!
|
|
|
2498 |
RATIO = -S(I,1)/(6.0*S(I,2))
|
|
|
2499 |
S(I+1,2) = S(I,1)/3.0
|
|
|
2500 |
!
|
|
|
2501 |
ELSE
|
|
|
2502 |
!
|
|
|
2503 |
RATIO = -(S(I,1)/6.0)/S(I,2)
|
|
|
2504 |
S(I+1,2) = S(I,1)*((1.0 - ZETA*ALPHA)*0.5*FACTOR-S(I,6))
|
|
|
2505 |
ENDIF
|
|
|
2506 |
!
|
|
|
2507 |
! ENDE DER SCHLEIFE UEBER I
|
|
|
2508 |
!
|
|
|
2509 |
I = I + 1
|
|
|
2510 |
IF(I.LT.NTAUM1) GO TO 10
|
|
|
2511 |
S(I,5) = 0.5
|
|
|
2512 |
!
|
|
|
2513 |
! DIE BEIDEN LETZTEN GLEICHUNGEN ERZWINGEN STETIGKEIT DER
|
|
|
2514 |
! 1. UND 3. ABLEITUNG IN TAU(NTAU-1)
|
|
|
2515 |
!
|
|
|
2516 |
ENTRY = RATIO*S(I-1,3) + S(I,2) + S(I,1)/3.0
|
|
|
2517 |
S(I+1,2) = S(I,1)/6.0
|
|
|
2518 |
S(I+1,4) = RATIO*S(I-1,4) + S(I,4)
|
|
|
2519 |
IF (Z-0.5 .LT. 0.) THEN
|
|
|
2520 |
RATIO = S(I,1)*6.0*S(I-1,6)*ALPHA/S(I-1,2)
|
|
|
2521 |
S(I,2) = RATIO*S(I-1,3) + S(I,1) + S(I-1,1)
|
|
|
2522 |
S(I,3) = -S(I-1,1)
|
|
|
2523 |
!
|
|
|
2524 |
ELSE IF (Z-0.5 .EQ. 0.) THEN
|
|
|
2525 |
!
|
|
|
2526 |
RATIO = S(I,1)/S(I-1,2)
|
|
|
2527 |
S(I,2) = RATIO*S(I-1,3) + S(I,1) + S(I-1,1)
|
|
|
2528 |
S(I,3) = -S(I-1,1)
|
|
|
2529 |
!
|
|
|
2530 |
ELSE
|
|
|
2531 |
!
|
|
|
2532 |
FACTR2 = ONEMZT*(ALPHA*(ONEMZT**2-1.0)+1.0)/ &
|
|
|
2533 |
(ALPHA*(ONEMZT**3-1.0)+1.0)
|
|
|
2534 |
RATIO = FACTR2*S(I,1)/S(I-1,2)
|
|
|
2535 |
S(I,2) = RATIO*S(I-1,3) + FACTR2*S(I-1,1) + S(I,1)
|
|
|
2536 |
S(I,3) = -FACTR2*S(I-1,1)
|
|
|
2537 |
ENDIF
|
|
|
2538 |
!
|
|
|
2539 |
! ELIMINATION VON XI
|
|
|
2540 |
!
|
|
|
2541 |
S(I,4) = RATIO*S(I-1,4)
|
|
|
2542 |
RATIO = -ENTRY/S(I,2)
|
|
|
2543 |
S(I+1,2) = RATIO*S(I,3) + S(I+1,2)
|
|
|
2544 |
S(I+1,4) = RATIO*S(I,4) + S(I+1,4)
|
|
|
2545 |
!
|
|
|
2546 |
! RUECKSUBSTITUTION
|
|
|
2547 |
!
|
|
|
2548 |
S(NTAU,4) = S(NTAU,4)/S(NTAU,2)
|
|
|
2549 |
50 CONTINUE
|
|
|
2550 |
S(I,4) = (S(I,4) - S(I,3)*S(I+1,4))/S(I,2)
|
|
|
2551 |
I = I - 1
|
|
|
2552 |
IF(I.GT.1) GO TO 50
|
|
|
2553 |
S(1,4) = (S(1,4) -S(1,3)*S(2,4)-ENTRY3*S(3,4))/S(1,2)
|
|
|
2554 |
!
|
|
|
2555 |
! ERZEUGEN DER POLYNOM-TEILE
|
|
|
2556 |
!
|
|
|
2557 |
BREAK(1) = TAU(1)
|
|
|
2558 |
L = 1
|
|
|
2559 |
DO 70 I = 1,NTAUM1
|
|
|
2560 |
COEF(1,L) = GTAU(I)
|
|
|
2561 |
COEF(3,L) = S(I,4)
|
|
|
2562 |
DIVDIF = (GTAU(I+1) - GTAU(I))/S(I,1)
|
|
|
2563 |
Z = S(I,5)
|
|
|
2564 |
!
|
|
|
2565 |
IF (Z-0.5 .LT. 0.) THEN
|
|
|
2566 |
IF(Z.EQ.0.0) GO TO 65
|
|
|
2567 |
ZETA = GAM*Z
|
|
|
2568 |
ONEMZT = 1.0 -ZETA
|
|
|
2569 |
C = S(I+1,4)/6.0
|
|
|
2570 |
D = S(I,4)*S(I,6)
|
|
|
2571 |
L = L + 1
|
|
|
2572 |
DEL = ZETA*S(I,1)
|
|
|
2573 |
BREAK(L) = TAU(I) + DEL
|
|
|
2574 |
ZT2 = ZETA**2
|
|
|
2575 |
ALPHA = ALPH(ONEMZT)
|
|
|
2576 |
FACTOR = ONEMZT**2*ALPHA
|
|
|
2577 |
COEF(1,L) = GTAU(I) + DIVDIF*DEL+S(I,1)**2*(D*ONEMZT*(FACTOR- &
|
|
|
2578 |
1.0) + C*ZETA*(ZT2 - 1.0))
|
|
|
2579 |
COEF(2,L) = DIVDIF + S(I,1)*(D*(1.0-3.0*FACTOR) + C*(3.0*ZT2- &
|
|
|
2580 |
1.0))
|
|
|
2581 |
COEF(3,L) = 6.0*(D*ALPHA*ONEMZT + C*ZETA)
|
|
|
2582 |
COEF(4,L) = 6.0*(C - D*ALPHA)/S(I,1)
|
|
|
2583 |
COEF(4,L-1) = COEF(4,L) -6.0*D*(1.0-ALPHA)/(DEL*ZT2)
|
|
|
2584 |
COEF(2,L-1) = COEF(2,L) - DEL*(COEF(3,L)-DEL*0.5*COEF(4,L-1))
|
|
|
2585 |
GO TO 68
|
|
|
2586 |
!
|
|
|
2587 |
ELSE IF (Z-0.5 .EQ. 0.) THEN
|
|
|
2588 |
!
|
|
|
2589 |
COEF(2,L) = DIVDIF - S(I,1)*(2.0*S(I,4) + S(I+1,4))/6.0
|
|
|
2590 |
COEF(4,L) = (S(I+1,4) - S(I,4))/S(I,1)
|
|
|
2591 |
GO TO 68
|
|
|
2592 |
!
|
|
|
2593 |
ELSE
|
|
|
2594 |
!
|
|
|
2595 |
ONEMZT = GAM*(1.0 - Z)
|
|
|
2596 |
IF(ONEMZT.EQ.0.0) GO TO 65
|
|
|
2597 |
ZETA = 1.0 - ONEMZT
|
|
|
2598 |
ALPHA = ALPH(ZETA)
|
|
|
2599 |
C = S(I+1,4)*S(I,6)
|
|
|
2600 |
D = S(I,4)/6.0
|
|
|
2601 |
DEL = ZETA*S(I,1)
|
|
|
2602 |
BREAK(L+1) = TAU(I) + DEL
|
|
|
2603 |
COEF(2,L) = DIVDIF -S(I,1)*(2.0*D + C)
|
|
|
2604 |
COEF(4,L) = 6.0*(C*ALPHA - D)/S(I,1)
|
|
|
2605 |
L = L + 1
|
|
|
2606 |
COEF(4,L) = COEF(4,L-1) + 6.0*(1.0-ALPHA)*C/(S(I,1)*ONEMZT**3)
|
|
|
2607 |
COEF(3,L) = COEF(3,L-1) + DEL* COEF(4,L-1)
|
|
|
2608 |
COEF(2,L) = COEF(2,L-1) + DEL*(COEF(3,L-1)+DEL*0.5*COEF(4,L-1))
|
|
|
2609 |
COEF(1,L) = COEF(1,L-1) + DEL*(COEF(2,L-1)+DEL*0.5* &
|
|
|
2610 |
(COEF(3,L-1) + (DEL/3.0)*COEF(4,L-1)))
|
|
|
2611 |
GO TO 68
|
|
|
2612 |
ENDIF
|
|
|
2613 |
!
|
|
|
2614 |
!
|
|
|
2615 |
65 CONTINUE
|
|
|
2616 |
COEF(2,L) = DIVDIF
|
|
|
2617 |
COEF(3,L) = 0.0
|
|
|
2618 |
COEF(4,L) = 0.0
|
|
|
2619 |
68 CONTINUE
|
|
|
2620 |
!
|
|
|
2621 |
L = L + 1
|
|
|
2622 |
BREAK(L) = TAU(I+1)
|
|
|
2623 |
70 CONTINUE
|
|
|
2624 |
|
|
|
2625 |
IFLAG = 0
|
|
|
2626 |
RETURN
|
|
|
2627 |
!
|
|
|
2628 |
999 IFLAG = 2
|
|
|
2629 |
|
|
|
2630 |
END SUBROUTINE tautsp
|
|
|
2631 |
|
|
|
2632 |
!==============================================================================
|
|
|
2633 |
!==============================================================================
|
|
|
2634 |
|
|
|
2635 |
!------------------------------------------------------------------------------
|
|
|
2636 |
|
|
|
2637 |
SUBROUTINE tautsp2D (TAU, GTAU, NTAU, NI, IMIN, IMAX, NTAUMAX, GAMMA, &
|
|
|
2638 |
S, BREAK, COEF, L_VEC, IFLAG)
|
|
|
2639 |
|
|
|
2640 |
!------------------------------------------------------------------------------
|
|
|
2641 |
!
|
|
|
2642 |
! Description:
|
|
|
2643 |
! BERECHNET DEN TAUT-SPLINE FUER DIE DATEN: TAU(I),GTAU(I),I=1,.,NTAU
|
|
|
2644 |
!
|
|
|
2645 |
! Method:
|
|
|
2646 |
! WENN GAMMA GT.0 WERDEN ZUSAETZLICHE KNOTEN BERECHNET
|
|
|
2647 |
! GAMMA I.A. 2.5 BZW. 5.5
|
|
|
2648 |
!
|
|
|
2649 |
! BREAK,COEF,L,K GEBEN DIE PP-DARSTELLUNG DER INTERPOLATIONSFUNKTION
|
|
|
2650 |
!
|
|
|
2651 |
! FUER BREAK(I).LE.X.LE.BREAK(I+1) HAT DIE INTERPOLATIONSFUNKTION
|
|
|
2652 |
! FOLGENDE FORM
|
|
|
2653 |
!
|
|
|
2654 |
! F(X)=COEF(1,I)+DX(COEF(2,I)+DX/2(COEF(3,I)+DX/3(COEF(4,I)))
|
|
|
2655 |
! MIT DX=X-BREAK(I) UND I=1,...,L
|
|
|
2656 |
!
|
|
|
2657 |
! IFLAG=0 KEIN FEHLER IN TAUTSP
|
|
|
2658 |
! IFLAG=2 FALSCHER INPUT
|
|
|
2659 |
!
|
|
|
2660 |
! S(NTAU,6) WORK-ARRAY
|
|
|
2661 |
!
|
|
|
2662 |
!==============================================================================
|
|
|
2663 |
|
|
|
2664 |
INTEGER(KIND=iintegers), INTENT(IN) :: &
|
|
|
2665 |
NI, IMIN, IMAX, NTAUMAX
|
|
|
2666 |
|
|
|
2667 |
INTEGER(KIND=iintegers), INTENT(IN) :: &
|
|
|
2668 |
NTAU(NI)
|
|
|
2669 |
|
|
|
2670 |
REAL (KIND=ireals), INTENT(IN) :: &
|
|
|
2671 |
GTAU(NI,NTAUMAX), &
|
|
|
2672 |
TAU (NI,NTAUMAX), &
|
|
|
2673 |
GAMMA
|
|
|
2674 |
|
|
|
2675 |
REAL (KIND=ireals), INTENT(OUT) :: &
|
|
|
2676 |
BREAK(NI,*), &
|
|
|
2677 |
COEF (NI,4,*), &
|
|
|
2678 |
S (NI,NTAUMAX,6)
|
|
|
2679 |
|
|
|
2680 |
INTEGER(KIND=iintegers), INTENT(OUT) :: &
|
|
|
2681 |
L_VEC(NI), &
|
|
|
2682 |
IFLAG
|
|
|
2683 |
|
|
|
2684 |
! Local Variables
|
|
|
2685 |
INTEGER(KIND=iintegers) I,J,K,L,METHOD,NTAUM1, mb_err_indx_i,mb_err_indx_k
|
|
|
2686 |
|
|
|
2687 |
REAL(KIND=ireals) C,D,DEL,DENOM,DIVDIF,ENTRY,FACTR2,GAM, &
|
|
|
2688 |
ONEMG3,SIXTH,TEMP,X,ALPH
|
|
|
2689 |
|
|
|
2690 |
REAL(KIND=ireals) :: &
|
|
|
2691 |
RATIO_VEC (NI), &
|
|
|
2692 |
Z_VEC (NI), &
|
|
|
2693 |
ZETA_VEC (NI), &
|
|
|
2694 |
ZT2_VEC (NI), &
|
|
|
2695 |
ALPHA_VEC (NI), &
|
|
|
2696 |
FACTOR_VEC (NI), &
|
|
|
2697 |
ONEMZT_VEC (NI), &
|
|
|
2698 |
ENTRY3 (NI)
|
|
|
2699 |
|
|
|
2700 |
!==============================================================================
|
|
|
2701 |
!
|
|
|
2702 |
ALPH(X)= MIN (1.0_ireals,ONEMG3/X)
|
|
|
2703 |
!
|
|
|
2704 |
! TEST DER INPUTPARAMETER
|
|
|
2705 |
!
|
|
|
2706 |
mb_err_indx_i = -1
|
|
|
2707 |
DO i = IMIN, IMAX
|
|
|
2708 |
IF (NTAU(i) .LT. 4) mb_err_indx_i = i
|
|
|
2709 |
ENDDO
|
|
|
2710 |
|
|
|
2711 |
IF ( mb_err_indx_i /= -1 ) THEN
|
|
|
2712 |
PRINT *, ' NTAU =', NTAU(mb_err_indx_i), mb_err_indx_i, &
|
|
|
2713 |
': MUST BE BIGGER THAN 4'
|
|
|
2714 |
! PRINT 600,NTAU(i)
|
|
|
2715 |
!600 FORMAT(' NTAU =',I4,' NTAU SOLL GROESSER ALS 4 SEIN')
|
|
|
2716 |
GO TO 999
|
|
|
2717 |
ENDIF
|
|
|
2718 |
|
|
|
2719 |
!
|
|
|
2720 |
! BERECHNUNG VON DELTA TAU UND DER 1. UND 2.ABLEITUNG DER DATEN
|
|
|
2721 |
!
|
|
|
2722 |
mb_err_indx_i = -1
|
|
|
2723 |
mb_err_indx_k = -1
|
|
|
2724 |
|
|
|
2725 |
DO k = 1, NTAUMAX
|
|
|
2726 |
DO i = IMIN, IMAX
|
|
|
2727 |
IF (k <= NTAU(i)-1) THEN
|
|
|
2728 |
S(i,k,1)=TAU(i,k+1)-TAU(i,k)
|
|
|
2729 |
IF (S(i,k,1) .LE. 0.0) THEN
|
|
|
2730 |
mb_err_indx_i = i
|
|
|
2731 |
mb_err_indx_k = k
|
|
|
2732 |
ELSE
|
|
|
2733 |
S(i,k+1,4) = (GTAU(i,k+1) - GTAU(i,k))/S(i,k,1)
|
|
|
2734 |
ENDIF
|
|
|
2735 |
ENDIF
|
|
|
2736 |
ENDDO
|
|
|
2737 |
|
|
|
2738 |
IF ( mb_err_indx_i /= -1 .OR. mb_err_indx_k /= -1) THEN
|
|
|
2739 |
PRINT 610,mb_err_indx_i,mb_err_indx_k, &
|
|
|
2740 |
& TAU(mb_err_indx_i,mb_err_indx_k),&
|
|
|
2741 |
& TAU(mb_err_indx_i,mb_err_indx_k+1)
|
|
|
2742 |
610 FORMAT(' PUNKT ',2I3,' UND DIE NAECHSTEN',2E15.6,'SIND IN&
|
|
|
2743 |
& FALSCHER REIHENFOLGE')
|
|
|
2744 |
GO TO 999
|
|
|
2745 |
ENDIF
|
|
|
2746 |
ENDDO
|
|
|
2747 |
|
|
|
2748 |
DO k = 1, NTAUMAX
|
|
|
2749 |
DO i = IMIN, IMAX
|
|
|
2750 |
IF (k >= 2 .AND. k <= NTAU(i)-1) THEN
|
|
|
2751 |
S(i,k,4) = S(i,k+1,4) - S(i,k,4)
|
|
|
2752 |
ENDIF
|
|
|
2753 |
ENDDO
|
|
|
2754 |
ENDDO
|
|
|
2755 |
!
|
|
|
2756 |
! 2.ABLEITUNG VON GTAU AN DEN PUNKTEN TAU
|
|
|
2757 |
!
|
|
|
2758 |
DO i = IMIN, IMAX
|
|
|
2759 |
S(i,2,2) = S(i,1,1)/3.0
|
|
|
2760 |
ENDDO
|
|
|
2761 |
|
|
|
2762 |
SIXTH = 1.0/6.0
|
|
|
2763 |
METHOD = 2
|
|
|
2764 |
GAM = GAMMA
|
|
|
2765 |
IF(GAM .LE. 0.0) METHOD = 1
|
|
|
2766 |
IF(GAM .GT. 3.0) THEN
|
|
|
2767 |
METHOD = 3
|
|
|
2768 |
GAM = GAM - 3.0
|
|
|
2769 |
ENDIF
|
|
|
2770 |
ONEMG3=1.0 - GAM/3.0
|
|
|
2771 |
!
|
|
|
2772 |
! SCHLEIFE UEBER K
|
|
|
2773 |
!
|
|
|
2774 |
DO k = 2, NTAUMAX
|
|
|
2775 |
DO i = IMIN, IMAX
|
|
|
2776 |
IF ( k <= NTAU(i)-2 ) THEN
|
|
|
2777 |
|
|
|
2778 |
Z_VEC(i)=0.5
|
|
|
2779 |
IF (METHOD /= 1) THEN
|
|
|
2780 |
IF ( ((METHOD == 2) .AND. &
|
|
|
2781 |
(S(i,k,4)*S(i,k+1,4) >= 0.0)) .OR. (METHOD == 3) ) THEN
|
|
|
2782 |
TEMP = ABS(S(i,k+1,4))
|
|
|
2783 |
DENOM = ABS(S(i,k,4)) + TEMP
|
|
|
2784 |
IF (DENOM /= 0.0) THEN
|
|
|
2785 |
Z_VEC(i) = TEMP/DENOM
|
|
|
2786 |
IF(ABS(Z_VEC(i)-0.5).LE.SIXTH) Z_VEC(i)=0.5
|
|
|
2787 |
ENDIF
|
|
|
2788 |
ENDIF
|
|
|
2789 |
ENDIF
|
|
|
2790 |
S(i,k,5) = Z_VEC(i)
|
|
|
2791 |
!
|
|
|
2792 |
! ERSTELLEN EINES TEILES DER I-TEN GLEICHUNG
|
|
|
2793 |
!
|
|
|
2794 |
IF (Z_VEC(i)-0.5 .LT. 0.) THEN
|
|
|
2795 |
ZETA_VEC(i) = GAM*Z_VEC(i)
|
|
|
2796 |
ONEMZT_VEC(i) = 1.0 - ZETA_VEC(i)
|
|
|
2797 |
ZT2_VEC(i) = ZETA_VEC(i)**2
|
|
|
2798 |
ALPHA_VEC(i) = ALPH(ONEMZT_VEC(i))
|
|
|
2799 |
FACTOR_VEC(i) = ZETA_VEC(i) / &
|
|
|
2800 |
(ALPHA_VEC(i)*(ZT2_VEC(i) - 1.0) + 1.0)
|
|
|
2801 |
S(i,k,6) = ZETA_VEC(i)*FACTOR_VEC(i)/6.0
|
|
|
2802 |
S(i,k,2) = S(i,k,2) + S(i,k,1) * &
|
|
|
2803 |
((1.0-ALPHA_VEC(i)*ONEMZT_VEC(i)) &
|
|
|
2804 |
*FACTOR_VEC(i)*0.5-S(i,k,6))
|
|
|
2805 |
IF(S(i,k,2).LE.0.0) S(i,k,2) = 1.0
|
|
|
2806 |
S(i,k,3) = S(i,k,1)/6.0
|
|
|
2807 |
!
|
|
|
2808 |
ELSE IF (Z_VEC(i)-0.5 .EQ. 0.) THEN
|
|
|
2809 |
!
|
|
|
2810 |
S(i,k,2) = S(i,k,2) + S(i,k,1)/3.0
|
|
|
2811 |
S(i,k,3) = S(i,k,1)/6.0
|
|
|
2812 |
!
|
|
|
2813 |
ELSE
|
|
|
2814 |
!
|
|
|
2815 |
ONEMZT_VEC(i) = GAM*(1.0 - Z_VEC(i))
|
|
|
2816 |
ZETA_VEC(i) = 1.0 - ONEMZT_VEC(i)
|
|
|
2817 |
ALPHA_VEC(i) = ALPH(ZETA_VEC(i))
|
|
|
2818 |
FACTOR_VEC(i) = ONEMZT_VEC(i) / &
|
|
|
2819 |
(1.0 - ALPHA_VEC(i)*ZETA_VEC(i)*(1.0 + ONEMZT_VEC(i)))
|
|
|
2820 |
S(i,k,6) = ONEMZT_VEC(i)*FACTOR_VEC(i)/6.0
|
|
|
2821 |
S(i,k,2) = S(i,k,2) + S(i,k,1)/3.0
|
|
|
2822 |
S(i,k,3) = S(i,k,6) * S(i,k,1)
|
|
|
2823 |
ENDIF
|
|
|
2824 |
!
|
|
|
2825 |
IF (k == 2) THEN
|
|
|
2826 |
S(i,1,5) = 0.5
|
|
|
2827 |
!
|
|
|
2828 |
! DIE ERSTEN BEIDEN GLEICHUNGEN ERZWINGEN STETIGKEIT DER 1. UND 3.AB-
|
|
|
2829 |
! LEITUNG IN TAU(i,k)
|
|
|
2830 |
!
|
|
|
2831 |
S(i,1,2) = S(i,1,1)/6.0
|
|
|
2832 |
S(i,1,3) = S(i,2,2)
|
|
|
2833 |
ENTRY3(i) = S(i,2,3)
|
|
|
2834 |
|
|
|
2835 |
IF (Z_VEC(i)-0.5 .LT. 0.) THEN
|
|
|
2836 |
FACTR2 = ZETA_VEC(i)*(ALPHA_VEC(i) &
|
|
|
2837 |
*(ZT2_VEC(i)-1.0)+1.0)/ &
|
|
|
2838 |
(ALPHA_VEC(i)*(ZETA_VEC(i)*ZT2_VEC(i)-1.0) + 1.0)
|
|
|
2839 |
RATIO_VEC(i) = FACTR2*S(i,2,1)/S(i,1,2)
|
|
|
2840 |
S(i,2,2) = FACTR2*S(i,2,1) + S(i,1,1)
|
|
|
2841 |
S(i,2,3) = -FACTR2 * S(i,1,1)
|
|
|
2842 |
!
|
|
|
2843 |
ELSE IF (Z_VEC(i)-0.5 .EQ. 0.) THEN
|
|
|
2844 |
!
|
|
|
2845 |
RATIO_VEC(i) = S(i,2,1)/S(i,1,2)
|
|
|
2846 |
S(i,2,2) = S(i,2,1) + S(i,1,1)
|
|
|
2847 |
S(i,2,3) = -S(i,1,1)
|
|
|
2848 |
!
|
|
|
2849 |
ELSE
|
|
|
2850 |
!
|
|
|
2851 |
RATIO_VEC(i) = S(i,2,1)/S(i,1,2)
|
|
|
2852 |
S(i,2,2) = S(i,2,1) + S(i,1,1)
|
|
|
2853 |
S(i,2,3) = -S(i,1,1)*6.0*ALPHA_VEC(i)*S(i,2,6)
|
|
|
2854 |
ENDIF
|
|
|
2855 |
!
|
|
|
2856 |
! ELIMINATION DER 1.UNBEKANNTEN AUS DER 2.GLEICHUNG
|
|
|
2857 |
!
|
|
|
2858 |
S(i,2,2) = RATIO_VEC(i)*S(i,1,3) + S(i,2,2)
|
|
|
2859 |
S(i,2,3) = RATIO_VEC(i)*ENTRY3(i) + S(i,2,3)
|
|
|
2860 |
S(i,1,4) = S(i,2,4)
|
|
|
2861 |
S(i,2,4) = RATIO_VEC(i)*S(i,1,4)
|
|
|
2862 |
!
|
|
|
2863 |
|
|
|
2864 |
ELSE ! k > 2
|
|
|
2865 |
!
|
|
|
2866 |
S(i,k,2) = RATIO_VEC(i)*S(i,k-1,3) + S(i,k,2)
|
|
|
2867 |
S(i,k,4) = RATIO_VEC(i)*S(i,k-1,4) + S(i,k,4)
|
|
|
2868 |
ENDIF ! k == 2
|
|
|
2869 |
!
|
|
|
2870 |
! AUFSTELLEN DES TEILES DER NAECHSTEN GLEICHUNG,DER VOM I-TEN INTERVAL
|
|
|
2871 |
! ABHAENGT
|
|
|
2872 |
!
|
|
|
2873 |
IF (Z_VEC(i)-0.5 .LT. 0.) THEN
|
|
|
2874 |
RATIO_VEC(i) = -S(i,k,6)*S(i,k,1)/S(i,k,2)
|
|
|
2875 |
S(i,k+1,2) = S(i,k,1)/3.0
|
|
|
2876 |
!
|
|
|
2877 |
ELSE IF (Z_VEC(i)-0.5 .EQ. 0.) THEN
|
|
|
2878 |
!
|
|
|
2879 |
RATIO_VEC(i) = -S(i,k,1)/(6.0*S(i,k,2))
|
|
|
2880 |
S(i,k+1,2) = S(i,k,1)/3.0
|
|
|
2881 |
!
|
|
|
2882 |
ELSE
|
|
|
2883 |
!
|
|
|
2884 |
RATIO_VEC(i) = -(S(i,k,1)/6.0)/S(i,k,2)
|
|
|
2885 |
S(i,k+1,2) = S(i,k,1) * &
|
|
|
2886 |
((1.0 - ZETA_VEC(i)*ALPHA_VEC(i))* &
|
|
|
2887 |
0.5*FACTOR_VEC(i)-S(i,k,6))
|
|
|
2888 |
ENDIF
|
|
|
2889 |
!
|
|
|
2890 |
! ENDE DER SCHLEIFE UEBER k
|
|
|
2891 |
!
|
|
|
2892 |
ENDIF ! k <= NTAU(i)-2
|
|
|
2893 |
ENDDO ! i = IMIN, IMAX
|
|
|
2894 |
ENDDO ! k = 2, NTAUMAX
|
|
|
2895 |
|
|
|
2896 |
DO i = IMIN, IMAX
|
|
|
2897 |
k = NTAU(i)-1
|
|
|
2898 |
S(i,k,5) = 0.5
|
|
|
2899 |
|
|
|
2900 |
!
|
|
|
2901 |
! DIE BEIDEN LETZTEN GLEICHUNGEN ERZWINGEN STETIGKEIT DER
|
|
|
2902 |
! 1. UND 3. ABLEITUNG IN TAU(NTAU-1)
|
|
|
2903 |
!
|
|
|
2904 |
ENTRY = RATIO_VEC(i)*S(i,k-1,3) + S(i,k,2) + S(i,k,1)/3.0
|
|
|
2905 |
S(i,k+1,2) = S(i,k,1)/6.0
|
|
|
2906 |
S(i,k+1,4) = RATIO_VEC(i)*S(i,k-1,4) + S(i,k,4)
|
|
|
2907 |
IF (Z_VEC(i)-0.5 .LT. 0.) THEN
|
|
|
2908 |
RATIO_VEC(i) = S(i,k,1) * 6.0 * S(i,k-1,6) * &
|
|
|
2909 |
ALPHA_VEC(i)/S(i,k-1,2)
|
|
|
2910 |
S(i,k,2) = RATIO_VEC(i)*S(i,k-1,3) +S(i,k,1) + S(i,k-1,1)
|
|
|
2911 |
S(i,k,3) = -S(i,k-1,1)
|
|
|
2912 |
!
|
|
|
2913 |
ELSE IF (Z_VEC(i)-0.5 .EQ. 0.) THEN
|
|
|
2914 |
!
|
|
|
2915 |
RATIO_VEC(i) = S(i,k,1)/S(i,k-1,2)
|
|
|
2916 |
S(i,k,2) = RATIO_VEC(i)*S(i,k-1,3) + S(i,k,1) + S(i,k-1,1)
|
|
|
2917 |
S(i,k,3) = -S(i,k-1,1)
|
|
|
2918 |
!
|
|
|
2919 |
ELSE
|
|
|
2920 |
!
|
|
|
2921 |
FACTR2 = ONEMZT_VEC(i) * (ALPHA_VEC(i) * &
|
|
|
2922 |
(ONEMZT_VEC(i)**2-1.0)+1.0) / &
|
|
|
2923 |
(ALPHA_VEC(i)*(ONEMZT_VEC(i)**3-1.0)+1.0)
|
|
|
2924 |
RATIO_VEC(i) = FACTR2*S(i,k,1)/S(i,k-1,2)
|
|
|
2925 |
S(i,k,2) = RATIO_VEC(i)*S(i,k-1,3) + FACTR2*S(i,k-1,1) &
|
|
|
2926 |
+ S(i,k,1)
|
|
|
2927 |
S(i,k,3) = -FACTR2*S(i,k-1,1)
|
|
|
2928 |
ENDIF
|
|
|
2929 |
!
|
|
|
2930 |
! ELIMINATION VON XI
|
|
|
2931 |
!
|
|
|
2932 |
S(i,k,4) = RATIO_VEC(i)*S(i,k-1,4)
|
|
|
2933 |
RATIO_VEC(i) = -ENTRY/S(i,k,2)
|
|
|
2934 |
S(i,k+1,2) = RATIO_VEC(i)*S(i,k,3) + S(i,k+1,2)
|
|
|
2935 |
S(i,k+1,4) = RATIO_VEC(i)*S(i,k,4) + S(i,k+1,4)
|
|
|
2936 |
ENDDO ! i = IMIN, IMAX
|
|
|
2937 |
|
|
|
2938 |
!
|
|
|
2939 |
! RUECKSUBSTITUTION
|
|
|
2940 |
!
|
|
|
2941 |
DO i = IMIN, IMAX
|
|
|
2942 |
S(i,NTAU(i),4) = S(i,NTAU(i),4)/S(i,NTAU(i),2)
|
|
|
2943 |
ENDDO
|
|
|
2944 |
|
|
|
2945 |
DO k = NTAUMAX,2,-1
|
|
|
2946 |
DO i = IMIN, IMAX
|
|
|
2947 |
IF (k <= NTAU(i)-1) THEN
|
|
|
2948 |
S(i,k,4) = (S(i,k,4) - S(i,k,3)*S(i,k+1,4))/S(i,k,2)
|
|
|
2949 |
ENDIF
|
|
|
2950 |
ENDDO
|
|
|
2951 |
ENDDO
|
|
|
2952 |
|
|
|
2953 |
DO i = IMIN, IMAX
|
|
|
2954 |
S(i,1,4) = (S(i,1,4) - S(i,1,3)*S(i,2,4) - &
|
|
|
2955 |
ENTRY3(i)*S(i,3,4))/S(i,1,2)
|
|
|
2956 |
ENDDO
|
|
|
2957 |
!
|
|
|
2958 |
! ERZEUGEN DER POLYNOM-TEILE
|
|
|
2959 |
!
|
|
|
2960 |
DO i = IMIN, IMAX
|
|
|
2961 |
BREAK(i,1) = TAU(i,1)
|
|
|
2962 |
L_VEC(i) = 1
|
|
|
2963 |
ENDDO
|
|
|
2964 |
|
|
|
2965 |
DO k = 1, NTAUMAX
|
|
|
2966 |
DO i = IMIN, IMAX
|
|
|
2967 |
IF ( k <= NTAU(i)-1) THEN
|
|
|
2968 |
L = L_VEC(i)
|
|
|
2969 |
COEF(i,1,L) = GTAU(i,k)
|
|
|
2970 |
COEF(i,3,L) = S(i,k,4)
|
|
|
2971 |
DIVDIF = (GTAU(i,k+1) - GTAU(i,k))/S(i,k,1)
|
|
|
2972 |
Z_VEC(i) = S(i,k,5)
|
|
|
2973 |
!
|
|
|
2974 |
IF (Z_VEC(i)-0.5 .LT. 0.) THEN
|
|
|
2975 |
! US avoid division by 0, if Z_VEC is veeeery small
|
|
|
2976 |
! by treating it as 0
|
|
|
2977 |
IF(ABS(Z_VEC(i)) < 1E-50_ireals) THEN
|
|
|
2978 |
COEF(i,2,L) = DIVDIF
|
|
|
2979 |
COEF(i,3,L) = 0.0
|
|
|
2980 |
COEF(i,4,L) = 0.0
|
|
|
2981 |
ELSE
|
|
|
2982 |
ZETA_VEC(i) = GAM*Z_VEC(i)
|
|
|
2983 |
ONEMZT_VEC(i) = 1.0 -ZETA_VEC(i)
|
|
|
2984 |
C = S(i,k+1,4)/6.0
|
|
|
2985 |
D = S(i,k,4)*S(i,k,6)
|
|
|
2986 |
L = L + 1
|
|
|
2987 |
DEL = ZETA_VEC(i)*S(i,k,1)
|
|
|
2988 |
BREAK(i,L) = TAU(i,k) + DEL
|
|
|
2989 |
ZT2_VEC(i) = ZETA_VEC(i)**2
|
|
|
2990 |
! ALPHA_VEC(i) = ALPH(ONEMZT_VEC(i))
|
|
|
2991 |
ALPHA_VEC(i) = MIN(1.0_ireals,ONEMG3/ONEMZT_VEC(i))
|
|
|
2992 |
FACTOR_VEC(i) = ONEMZT_VEC(i)**2*ALPHA_VEC(i)
|
|
|
2993 |
COEF(i,1,L) = GTAU(i,k) + DIVDIF*DEL + &
|
|
|
2994 |
S(i,k,1)**2*(D*ONEMZT_VEC(i)*(FACTOR_VEC(i)- &
|
|
|
2995 |
1.0) + C*ZETA_VEC(i)*(ZT2_VEC(i) - 1.0))
|
|
|
2996 |
COEF(i,2,L) = DIVDIF + S(i,k,1) * &
|
|
|
2997 |
(D*(1.0-3.0*FACTOR_VEC(i)) + C*(3.0*ZT2_VEC(i)- &
|
|
|
2998 |
1.0))
|
|
|
2999 |
COEF(i,3,L) = 6.0*(D*ALPHA_VEC(i)*ONEMZT_VEC(i) &
|
|
|
3000 |
+ C*ZETA_VEC(i))
|
|
|
3001 |
COEF(i,4,L) = 6.0*(C - D*ALPHA_VEC(i))/S(i,k,1)
|
|
|
3002 |
COEF(i,4,L-1) = COEF(i,4,L) -6.0*D* &
|
|
|
3003 |
(1.0-ALPHA_VEC(i))/(DEL*ZT2_VEC(i))
|
|
|
3004 |
COEF(i,2,L-1) = COEF(i,2,L) - &
|
|
|
3005 |
DEL*(COEF(i,3,L)-DEL*0.5*COEF(i,4,L-1))
|
|
|
3006 |
ENDIF
|
|
|
3007 |
!
|
|
|
3008 |
ELSE IF (Z_VEC(i)-0.5 .EQ. 0.) THEN
|
|
|
3009 |
!
|
|
|
3010 |
COEF(i,2,L) = DIVDIF - S(i,k,1) * &
|
|
|
3011 |
(2.0*S(i,k,4) + S(i,k+1,4))/6.0
|
|
|
3012 |
COEF(i,4,L) = (S(i,k+1,4) - S(i,k,4))/S(i,k,1)
|
|
|
3013 |
!
|
|
|
3014 |
ELSE
|
|
|
3015 |
!
|
|
|
3016 |
ONEMZT_VEC(i) = GAM*(1.0 - Z_VEC(i))
|
|
|
3017 |
IF(ONEMZT_VEC(i).EQ.0.0) THEN
|
|
|
3018 |
COEF(i,2,L) = DIVDIF
|
|
|
3019 |
COEF(i,3,L) = 0.0
|
|
|
3020 |
COEF(i,4,L) = 0.0
|
|
|
3021 |
ELSE
|
|
|
3022 |
ZETA_VEC(i) = 1.0 - ONEMZT_VEC(i)
|
|
|
3023 |
! ALPHA_VEC(i) = ALPH(ZETA_VEC(i))
|
|
|
3024 |
ALPHA_VEC(i) = MIN(1.0_ireals,ONEMG3/ZETA_VEC(i))
|
|
|
3025 |
C = S(i,k+1,4)*S(i,k,6)
|
|
|
3026 |
D = S(i,k,4)/6.0
|
|
|
3027 |
DEL = ZETA_VEC(i)*S(i,k,1)
|
|
|
3028 |
BREAK(i,L+1) = TAU(i,k) + DEL
|
|
|
3029 |
COEF(i,2,L) = DIVDIF -S(i,k,1)*(2.0*D + C)
|
|
|
3030 |
COEF(i,4,L) = 6.0*(C*ALPHA_VEC(i) - D)/S(i,k,1)
|
|
|
3031 |
L = L + 1
|
|
|
3032 |
COEF(i,4,L) = COEF(i,4,L-1) + 6.0 * &
|
|
|
3033 |
(1.0-ALPHA_VEC(i))*C/(S(i,k,1)*ONEMZT_VEC(i)**3)
|
|
|
3034 |
COEF(i,3,L) = COEF(i,3,L-1) + DEL* COEF(i,4,L-1)
|
|
|
3035 |
COEF(i,2,L) = COEF(i,2,L-1) + DEL*(COEF(i,3,L-1) &
|
|
|
3036 |
+DEL*0.5*COEF(i,4,L-1))
|
|
|
3037 |
COEF(i,1,L) = COEF(i,1,L-1) + DEL*(COEF(i,2,L-1) &
|
|
|
3038 |
+DEL*0.5* &
|
|
|
3039 |
(COEF(i,3,L-1) + (DEL/3.0)*COEF(i,4,L-1)))
|
|
|
3040 |
ENDIF
|
|
|
3041 |
ENDIF
|
|
|
3042 |
!
|
|
|
3043 |
L = L + 1
|
|
|
3044 |
BREAK(i,L) = TAU(i,k+1)
|
|
|
3045 |
L_VEC(i) = L
|
|
|
3046 |
ENDIF
|
|
|
3047 |
ENDDO ! i = IMIN, IMAX
|
|
|
3048 |
ENDDO ! k = 1, NTAUMAX
|
|
|
3049 |
|
|
|
3050 |
IFLAG = 0
|
|
|
3051 |
|
|
|
3052 |
RETURN
|
|
|
3053 |
!
|
|
|
3054 |
999 IFLAG = 2
|
|
|
3055 |
|
|
|
3056 |
END SUBROUTINE tautsp2D
|
|
|
3057 |
|
|
|
3058 |
!==============================================================================
|
|
|
3059 |
!==============================================================================
|
|
|
3060 |
|
|
|
3061 |
!------------------------------------------------------------------------------
|
|
|
3062 |
|
|
|
3063 |
SUBROUTINE to_upper ( string )
|
|
|
3064 |
|
|
|
3065 |
!-------------------------------------------------------------------------------
|
|
|
3066 |
!
|
|
|
3067 |
! Description:
|
|
|
3068 |
! Convert alphabetic characters in 'string' from lower to upper case
|
|
|
3069 |
!-------------------------------------------------------------------------------
|
|
|
3070 |
|
|
|
3071 |
IMPLICIT NONE
|
|
|
3072 |
|
|
|
3073 |
! Subroutine arguments:
|
|
|
3074 |
! --------------------
|
|
|
3075 |
CHARACTER (LEN=*), INTENT(INOUT) :: string
|
|
|
3076 |
|
|
|
3077 |
! Local parameters:
|
|
|
3078 |
! ----------------
|
|
|
3079 |
CHARACTER (LEN=26), PARAMETER :: UPPER="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
|
|
|
3080 |
CHARACTER (LEN=26), PARAMETER :: lower="abcdefghijklmnopqrstuvwxyz"
|
|
|
3081 |
|
|
|
3082 |
! Local variables:
|
|
|
3083 |
! ---------------
|
|
|
3084 |
INTEGER (KIND=iintegers) :: i, j
|
|
|
3085 |
!
|
|
|
3086 |
!------------ End of header ----------------------------------------------------
|
|
|
3087 |
|
|
|
3088 |
DO i = 1, LEN_TRIM(string)
|
|
|
3089 |
j = INDEX ( lower, string(i:i) )
|
|
|
3090 |
IF ( j > 0 ) string(i:i) = UPPER(j:j)
|
|
|
3091 |
END DO
|
|
|
3092 |
|
|
|
3093 |
END SUBROUTINE to_upper
|
|
|
3094 |
|
|
|
3095 |
|
|
|
3096 |
!==============================================================================
|
|
|
3097 |
!==============================================================================
|
|
|
3098 |
|
|
|
3099 |
!------------------------------------------------------------------------------
|
|
|
3100 |
|
|
|
3101 |
SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v)
|
|
|
3102 |
|
|
|
3103 |
|
|
|
3104 |
!------------------------------------------------------------------------------
|
|
|
3105 |
!
|
|
|
3106 |
! Description:
|
|
|
3107 |
! This routine converts the wind components u and v from the rotated system
|
|
|
3108 |
! to the real geographical system.
|
|
|
3109 |
!
|
|
|
3110 |
! Method:
|
|
|
3111 |
! Transformation formulas for converting between these two systems.
|
|
|
3112 |
!
|
|
|
3113 |
!------------------------------------------------------------------------------
|
|
|
3114 |
|
|
|
3115 |
! Parameter list:
|
|
|
3116 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
3117 |
urot, vrot, & ! wind components in the rotated grid
|
|
|
3118 |
rlat, rlon, & ! latitude and longitude in the true geographical system
|
|
|
3119 |
pollat, pollon ! latitude and longitude of the north pole of the
|
|
|
3120 |
! rotated grid
|
|
|
3121 |
|
|
|
3122 |
REAL (KIND=ireals), INTENT (OUT) :: &
|
|
|
3123 |
u, v ! wind components in the true geographical system
|
|
|
3124 |
|
|
|
3125 |
! Local variables
|
|
|
3126 |
|
|
|
3127 |
REAL (KIND=ireals) :: &
|
|
|
3128 |
zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
|
|
|
3129 |
|
|
|
3130 |
REAL (KIND=ireals) :: &
|
|
|
3131 |
zrpi18 = 57.2957795_ireals, & !
|
|
|
3132 |
zpir18 = 0.0174532925_ireals
|
|
|
3133 |
|
|
|
3134 |
!------------------------------------------------------------------------------
|
|
|
3135 |
! Begin subroutine uvrot2uv
|
|
|
3136 |
!------------------------------------------------------------------------------
|
|
|
3137 |
|
|
|
3138 |
! Converting from degree to radians
|
|
|
3139 |
zsinpol = SIN(pollat * zpir18)
|
|
|
3140 |
zcospol = COS(pollat * zpir18)
|
|
|
3141 |
zlonp = (pollon-rlon) * zpir18
|
|
|
3142 |
zlat = rlat * zpir18
|
|
|
3143 |
|
|
|
3144 |
zarg1 = zcospol*SIN(zlonp)
|
|
|
3145 |
zarg2 = zsinpol*COS(zlat) - zcospol*SIN(zlat)*COS(zlonp)
|
|
|
3146 |
znorm = 1.0/SQRT(zarg1**2 + zarg2**2)
|
|
|
3147 |
|
|
|
3148 |
! Convert the u- and v-components
|
|
|
3149 |
u = urot*zarg2*znorm + vrot*zarg1*znorm
|
|
|
3150 |
v = - urot*zarg1*znorm + vrot*zarg2*znorm
|
|
|
3151 |
|
|
|
3152 |
END SUBROUTINE uvrot2uv
|
|
|
3153 |
|
|
|
3154 |
!==============================================================================
|
|
|
3155 |
!==============================================================================
|
|
|
3156 |
|
|
|
3157 |
!------------------------------------------------------------------------------
|
|
|
3158 |
|
|
|
3159 |
SUBROUTINE uvrot2uv_vec(u, v, rlat, rlon, pollat, pollon, idim, jdim)
|
|
|
3160 |
|
|
|
3161 |
!------------------------------------------------------------------------------
|
|
|
3162 |
!
|
|
|
3163 |
! Description:
|
|
|
3164 |
! This routine converts the wind components u and v from the rotated
|
|
|
3165 |
! system to the real geographical system. This is the vectorized form
|
|
|
3166 |
! of the routine above, i.e. the computation is for a whole 2D field.
|
|
|
3167 |
!
|
|
|
3168 |
! Method:
|
|
|
3169 |
! Transformation formulas for converting between these two systems.
|
|
|
3170 |
!
|
|
|
3171 |
!------------------------------------------------------------------------------
|
|
|
3172 |
|
|
|
3173 |
!------------------------------------------------------------------------------
|
|
|
3174 |
! Parameter list:
|
|
|
3175 |
INTEGER (KIND=iintegers), INTENT(IN) :: &
|
|
|
3176 |
idim, jdim ! dimensions of the field
|
|
|
3177 |
|
|
|
3178 |
REAL (KIND=ireals), INTENT (INOUT) :: &
|
|
|
3179 |
u (idim,jdim), & ! wind components in the true geographical system
|
|
|
3180 |
v (idim,jdim) !
|
|
|
3181 |
|
|
|
3182 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
3183 |
rlat(idim,jdim),& ! coordinates in the true geographical system
|
|
|
3184 |
rlon(idim,jdim),& !
|
|
|
3185 |
pollat, pollon ! latitude and longitude of the north pole of the
|
|
|
3186 |
! rotated grid
|
|
|
3187 |
|
|
|
3188 |
! Local variables
|
|
|
3189 |
REAL (KIND=ireals) :: &
|
|
|
3190 |
zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm, zugeo, zvgeo
|
|
|
3191 |
|
|
|
3192 |
INTEGER (KIND=iintegers) :: i, j
|
|
|
3193 |
REAL (KIND=ireals) :: &
|
|
|
3194 |
zrpi18 = 57.2957795_ireals, & !
|
|
|
3195 |
zpir18 = 0.0174532925_ireals
|
|
|
3196 |
|
|
|
3197 |
!------------------------------------------------------------------------------
|
|
|
3198 |
! Begin Subroutine uvrot2uv_vec
|
|
|
3199 |
!------------------------------------------------------------------------------
|
|
|
3200 |
|
|
|
3201 |
! Converting from degree to radians
|
|
|
3202 |
zsinpol = SIN(pollat * zpir18)
|
|
|
3203 |
zcospol = COS(pollat * zpir18)
|
|
|
3204 |
|
|
|
3205 |
DO j = 1, jdim
|
|
|
3206 |
DO i = 1, idim
|
|
|
3207 |
|
|
|
3208 |
zlonp = (pollon-rlon(i,j)) * zpir18
|
|
|
3209 |
zlat = rlat(i,j) * zpir18
|
|
|
3210 |
|
|
|
3211 |
zarg1 = zcospol*SIN(zlonp)
|
|
|
3212 |
zarg2 = zsinpol*COS(zlat) - zcospol*SIN(zlat)*COS(zlonp)
|
|
|
3213 |
znorm = 1.0/SQRT(zarg1**2 + zarg2**2)
|
|
|
3214 |
|
|
|
3215 |
! Convert the u- and v-components
|
|
|
3216 |
zugeo = u(i,j)*zarg2*znorm + v(i,j)*zarg1*znorm
|
|
|
3217 |
zvgeo = -u(i,j)*zarg1*znorm + v(i,j)*zarg2*znorm
|
|
|
3218 |
u(i,j) = zugeo
|
|
|
3219 |
v(i,j) = zvgeo
|
|
|
3220 |
|
|
|
3221 |
ENDDO
|
|
|
3222 |
ENDDO
|
|
|
3223 |
|
|
|
3224 |
END SUBROUTINE uvrot2uv_vec
|
|
|
3225 |
|
|
|
3226 |
!==============================================================================
|
|
|
3227 |
!==============================================================================
|
|
|
3228 |
|
|
|
3229 |
!------------------------------------------------------------------------------
|
|
|
3230 |
|
|
|
3231 |
SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot)
|
|
|
3232 |
|
|
|
3233 |
!------------------------------------------------------------------------------
|
|
|
3234 |
!
|
|
|
3235 |
! Description:
|
|
|
3236 |
! This routine converts the wind components u and v from the real
|
|
|
3237 |
! geographical system to the rotated system.
|
|
|
3238 |
!
|
|
|
3239 |
! Method:
|
|
|
3240 |
! Transformation formulas for converting between these two systems.
|
|
|
3241 |
!
|
|
|
3242 |
!------------------------------------------------------------------------------
|
|
|
3243 |
|
|
|
3244 |
!------------------------------------------------------------------------------
|
|
|
3245 |
! Parameter list:
|
|
|
3246 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
3247 |
u , v , & ! wind components in the true geographical system
|
|
|
3248 |
rlat, rlon, & ! coordinates in the true geographical system
|
|
|
3249 |
pollat, pollon ! latitude and longitude of the north pole of the
|
|
|
3250 |
! rotated grid
|
|
|
3251 |
|
|
|
3252 |
REAL (KIND=ireals), INTENT (OUT) :: &
|
|
|
3253 |
urot, vrot ! wind components in the rotated grid
|
|
|
3254 |
|
|
|
3255 |
! Local variables
|
|
|
3256 |
|
|
|
3257 |
REAL (KIND=ireals) :: &
|
|
|
3258 |
zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
|
|
|
3259 |
|
|
|
3260 |
REAL (KIND=ireals) :: &
|
|
|
3261 |
zrpi18 = 57.2957795_ireals, & !
|
|
|
3262 |
zpir18 = 0.0174532925_ireals
|
|
|
3263 |
|
|
|
3264 |
!------------------------------------------------------------------------------
|
|
|
3265 |
! Begin Subroutine uv2uvrot
|
|
|
3266 |
!------------------------------------------------------------------------------
|
|
|
3267 |
|
|
|
3268 |
zsinpol = SIN(pollat * zpir18)
|
|
|
3269 |
zcospol = COS(pollat * zpir18)
|
|
|
3270 |
zlonp = (pollon-rlon) * zpir18
|
|
|
3271 |
zlat = rlat * zpir18
|
|
|
3272 |
|
|
|
3273 |
zarg1 = zcospol*SIN(zlonp)
|
|
|
3274 |
zarg2 = zsinpol*COS(zlat) - zcospol*SIN(zlat)*COS(zlonp)
|
|
|
3275 |
znorm = 1.0_ireals/SQRT( zarg1**2 + zarg2**2 )
|
|
|
3276 |
|
|
|
3277 |
! Transform the u and v wind components
|
|
|
3278 |
urot = u*zarg2*znorm - v*zarg1*znorm
|
|
|
3279 |
vrot = u*zarg1*znorm + v*zarg2*znorm
|
|
|
3280 |
|
|
|
3281 |
END SUBROUTINE uv2uvrot
|
|
|
3282 |
|
|
|
3283 |
!==============================================================================
|
|
|
3284 |
!==============================================================================
|
|
|
3285 |
|
|
|
3286 |
!------------------------------------------------------------------------------
|
|
|
3287 |
|
|
|
3288 |
SUBROUTINE uv2uvrot_vec(u, v, rlat, rlon, pollat, pollon, idim, jdim)
|
|
|
3289 |
|
|
|
3290 |
!------------------------------------------------------------------------------
|
|
|
3291 |
!
|
|
|
3292 |
! Description:
|
|
|
3293 |
! This routine converts the wind components u and v from the real
|
|
|
3294 |
! geographical system to the rotated system. This is the vectorized form
|
|
|
3295 |
! of the routine above, i.e. the computation is for a whole 2D field.
|
|
|
3296 |
!
|
|
|
3297 |
! Method:
|
|
|
3298 |
! Transformation formulas for converting between these two systems.
|
|
|
3299 |
!
|
|
|
3300 |
!------------------------------------------------------------------------------
|
|
|
3301 |
|
|
|
3302 |
!------------------------------------------------------------------------------
|
|
|
3303 |
! Parameter list:
|
|
|
3304 |
INTEGER (KIND=iintegers), INTENT(IN) :: &
|
|
|
3305 |
idim, jdim ! dimensions of the field
|
|
|
3306 |
|
|
|
3307 |
REAL (KIND=ireals), INTENT (INOUT) :: &
|
|
|
3308 |
u (idim,jdim), & ! wind components in the true geographical system
|
|
|
3309 |
v (idim,jdim) !
|
|
|
3310 |
|
|
|
3311 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
3312 |
rlat(idim,jdim),& ! coordinates in the true geographical system
|
|
|
3313 |
rlon(idim,jdim),& !
|
|
|
3314 |
pollat, pollon ! latitude and longitude of the north pole of the
|
|
|
3315 |
! rotated grid
|
|
|
3316 |
|
|
|
3317 |
! Local variables
|
|
|
3318 |
REAL (KIND=ireals) :: &
|
|
|
3319 |
zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm, zurot, zvrot
|
|
|
3320 |
|
|
|
3321 |
INTEGER (KIND=iintegers) :: i, j
|
|
|
3322 |
REAL (KIND=ireals) :: &
|
|
|
3323 |
zrpi18 = 57.2957795_ireals, & !
|
|
|
3324 |
zpir18 = 0.0174532925_ireals
|
|
|
3325 |
|
|
|
3326 |
!------------------------------------------------------------------------------
|
|
|
3327 |
! Begin Subroutine uv2uvrot_vec
|
|
|
3328 |
!------------------------------------------------------------------------------
|
|
|
3329 |
|
|
|
3330 |
zsinpol = SIN ( pollat * zpir18 )
|
|
|
3331 |
zcospol = COS ( pollat * zpir18 )
|
|
|
3332 |
|
|
|
3333 |
DO j = 1, jdim
|
|
|
3334 |
DO i = 1, idim
|
|
|
3335 |
|
|
|
3336 |
zlonp = ( pollon - rlon(i,j) ) * zpir18
|
|
|
3337 |
zlat = rlat(i,j) * zpir18
|
|
|
3338 |
|
|
|
3339 |
zarg1 = zcospol*SIN(zlonp)
|
|
|
3340 |
zarg2 = zsinpol*COS(zlat) - zcospol*SIN(zlat)*COS(zlonp)
|
|
|
3341 |
znorm = 1.0_ireals/SQRT( zarg1**2 + zarg2**2 )
|
|
|
3342 |
|
|
|
3343 |
! Transform the u and v wind components
|
|
|
3344 |
zurot = u(i,j)*zarg2*znorm - v(i,j)*zarg1*znorm
|
|
|
3345 |
zvrot = u(i,j)*zarg1*znorm + v(i,j)*zarg2*znorm
|
|
|
3346 |
u(i,j) = zurot
|
|
|
3347 |
v(i,j) = zvrot
|
|
|
3348 |
|
|
|
3349 |
ENDDO
|
|
|
3350 |
ENDDO
|
|
|
3351 |
|
|
|
3352 |
END SUBROUTINE uv2uvrot_vec
|
|
|
3353 |
|
|
|
3354 |
!==============================================================================
|
|
|
3355 |
!==============================================================================
|
|
|
3356 |
|
|
|
3357 |
!------------------------------------------------------------------------------
|
|
|
3358 |
|
|
|
3359 |
SUBROUTINE uv2df (u, v, d, f)
|
|
|
3360 |
|
|
|
3361 |
!------------------------------------------------------------------------------
|
|
|
3362 |
!
|
|
|
3363 |
! Description:
|
|
|
3364 |
! This routine computes wind speed and wind direction from the wind
|
|
|
3365 |
! components.
|
|
|
3366 |
!
|
|
|
3367 |
! Method:
|
|
|
3368 |
! Straightforward.
|
|
|
3369 |
!
|
|
|
3370 |
!------------------------------------------------------------------------------
|
|
|
3371 |
!
|
|
|
3372 |
! Parameter list:
|
|
|
3373 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
3374 |
u , v ! wind components in the true geographical system
|
|
|
3375 |
|
|
|
3376 |
REAL (KIND=ireals), INTENT (OUT) :: &
|
|
|
3377 |
f , d ! wind speed and wind direction
|
|
|
3378 |
|
|
|
3379 |
! Local variables
|
|
|
3380 |
|
|
|
3381 |
REAL (KIND=ireals) :: &
|
|
|
3382 |
zrpi18 = 57.2957795_ireals, & ! conversion from radians to degrees
|
|
|
3383 |
zsmall = 0.001_ireals
|
|
|
3384 |
|
|
|
3385 |
!------------------------------------------------------------------------------
|
|
|
3386 |
! Begin Subroutine uv2df
|
|
|
3387 |
!------------------------------------------------------------------------------
|
|
|
3388 |
|
|
|
3389 |
IF (ABS(u) > zsmall) THEN
|
|
|
3390 |
f = SQRT( u*u + v*v )
|
|
|
3391 |
d = v / u
|
|
|
3392 |
d = 180.0_ireals + SIGN( 90.0_ireals , u ) - ATAN( d ) *zrpi18
|
|
|
3393 |
ELSEIF (ABS(v) > zsmall) THEN
|
|
|
3394 |
f = ABS( v )
|
|
|
3395 |
d = 270.0_ireals - SIGN( 90.0_ireals , v )
|
|
|
3396 |
ELSE
|
|
|
3397 |
f = 0.0_ireals
|
|
|
3398 |
d = 0.0_ireals
|
|
|
3399 |
ENDIF
|
|
|
3400 |
|
|
|
3401 |
END SUBROUTINE uv2df
|
|
|
3402 |
|
|
|
3403 |
!==============================================================================
|
|
|
3404 |
!==============================================================================
|
|
|
3405 |
|
|
|
3406 |
!------------------------------------------------------------------------------
|
|
|
3407 |
|
|
|
3408 |
SUBROUTINE uv2df_vec (u, v, d, f, idim, jdim)
|
|
|
3409 |
|
|
|
3410 |
!------------------------------------------------------------------------------
|
|
|
3411 |
!
|
|
|
3412 |
! Description:
|
|
|
3413 |
! This routine computes wind speed and wind direction from the wind
|
|
|
3414 |
! components. This is the vectorized form of the routine above,
|
|
|
3415 |
! i.e. the computation is for a whole 2D field.
|
|
|
3416 |
!
|
|
|
3417 |
! Method:
|
|
|
3418 |
! Straightforward.
|
|
|
3419 |
!
|
|
|
3420 |
!------------------------------------------------------------------------------
|
|
|
3421 |
!
|
|
|
3422 |
! Parameter list:
|
|
|
3423 |
INTEGER (KIND=iintegers), INTENT(IN) :: &
|
|
|
3424 |
idim, jdim ! dimensions of the field
|
|
|
3425 |
|
|
|
3426 |
REAL (KIND=ireals), INTENT (IN) :: &
|
|
|
3427 |
u (idim,jdim) ,& ! wind components in the true geographical system
|
|
|
3428 |
v (idim,jdim) !
|
|
|
3429 |
|
|
|
3430 |
REAL (KIND=ireals), INTENT (OUT) :: &
|
|
|
3431 |
f (idim,jdim) ,& ! wind speed
|
|
|
3432 |
d (idim,jdim) ! wind direction
|
|
|
3433 |
|
|
|
3434 |
! Local variables
|
|
|
3435 |
|
|
|
3436 |
INTEGER (KIND=iintegers) :: i, j
|
|
|
3437 |
REAL (KIND=ireals) :: &
|
|
|
3438 |
zrpi18 = 57.2957795_ireals, & ! conversion from radians to degrees
|
|
|
3439 |
zsmall = 0.001_ireals
|
|
|
3440 |
|
|
|
3441 |
!------------------------------------------------------------------------------
|
|
|
3442 |
! Begin Subroutine uv2df_vec
|
|
|
3443 |
!------------------------------------------------------------------------------
|
|
|
3444 |
|
|
|
3445 |
DO j = 1, jdim
|
|
|
3446 |
DO i = 1, idim
|
|
|
3447 |
|
|
|
3448 |
IF (ABS(u(i,j)) > zsmall) THEN
|
|
|
3449 |
f (i,j) = SQRT( u(i,j)*u(i,j) + v(i,j)*v(i,j) )
|
|
|
3450 |
d (i,j) = 180.0_ireals + SIGN( 90.0_ireals , u(i,j) ) &
|
|
|
3451 |
- ATAN( v(i,j) / u(i,j) ) *zrpi18
|
|
|
3452 |
ELSEIF (ABS(v(i,j)) > zsmall) THEN
|
|
|
3453 |
f (i,j) = ABS( v(i,j) )
|
|
|
3454 |
d (i,j) = 270.0_ireals - SIGN( 90.0_ireals , v(i,j) )
|
|
|
3455 |
ELSE
|
|
|
3456 |
f (i,j) = 0.0_ireals
|
|
|
3457 |
d (i,j) = 0.0_ireals
|
|
|
3458 |
ENDIF
|
|
|
3459 |
|
|
|
3460 |
ENDDO
|
|
|
3461 |
ENDDO
|
|
|
3462 |
|
|
|
3463 |
END SUBROUTINE uv2df_vec
|
|
|
3464 |
|
|
|
3465 |
!==============================================================================
|
|
|
3466 |
!==============================================================================
|
|
|
3467 |
|
|
|
3468 |
END MODULE utilities
|