module_dfi.F
References to this file elsewhere.
1 MODULE module_dfi
2
3 USE module_domain
4 USE module_configure
5 USE module_wrf_error
6
7 implicit none
8
9 real, allocatable,dimension(:),public :: H
10
11 real, private :: da_dot
12
13 #ifdef WRF_PLUS
14 CONTAINS
15
16 SUBROUTINE jcdfi_init_coef ( grid , config_flags )
17
18 IMPLICIT NONE
19
20 ! Input data.
21
22 TYPE(domain) , intent(in) :: grid
23 TYPE (grid_config_rec_type), intent(in) :: config_flags
24 integer :: nsteps
25
26 if (mod ( INT(config_flags%interval_seconds / config_flags%dt) , 2 ).ne.0) &
27 CALL wrf_error_fatal( 'Total # of time steps should be EVEN' )
28
29 nsteps = config_flags%interval_seconds / config_flags%dt / 2
30
31 CALL JCDFI_COEF ( 0 , NSTEPS , config_flags%dt , grid%jcdfi_tauc , 7 )
32
33 END SUBROUTINE jcdfi_init_coef
34 !
35 !
36 !
37 SUBROUTINE jcdfi_zero_ad ( grid )
38
39 IMPLICIT NONE
40
41 ! Input data.
42
43 TYPE(domain) , intent(inout) :: grid
44
45 grid%em_a_u_2 = 0.0
46 grid%em_a_v_2 = 0.0
47 grid%em_a_t_2 = 0.0
48 grid%em_a_mu_2 = 0.0
49 grid%a_moist_2 = 0.0
50
51 END SUBROUTINE jcdfi_zero_ad
52 !
53 !
54 !
55 SUBROUTINE jcdfi_tl ( grid ,i)
56
57 IMPLICIT NONE
58
59 ! Input data.
60
61 TYPE(domain) , intent(inout) :: grid
62
63 integer, intent(in) :: i
64
65 grid%em_a_u_2 = grid%em_a_u_2 + h(i)*grid%em_g_u_2
66 grid%em_a_v_2 = grid%em_a_v_2 + h(i)*grid%em_g_v_2
67 grid%em_a_t_2 = grid%em_a_t_2 + h(i)*grid%em_g_t_2
68 grid%em_a_mu_2 = grid%em_a_mu_2 + h(i)*grid%em_g_mu_2
69 grid%a_moist_2 = grid%a_moist_2 + h(i)*grid%g_moist_2
70
71 END SUBROUTINE jcdfi_tl
72 !
73 !
74 !
75 SUBROUTINE jcdfi_output_forcing ( grid , config_flags )
76
77 IMPLICIT NONE
78
79 TYPE ( domain ) , intent(inout) :: grid
80 type(grid_config_rec_type) , intent(in) :: config_flags
81 real :: jc
82
83 grid%em_g_u_1 = grid%em_g_u_2
84 grid%em_g_u_2 = sqrt(grid%jcdfi_gama) * 1.0/grid%jcdfi_error_wind * grid%em_a_u_2
85
86 grid%em_g_v_1 = grid%em_g_v_2
87 grid%em_g_v_2 = sqrt(grid%jcdfi_gama) * 1.0/grid%jcdfi_error_wind * grid%em_a_v_2
88
89 grid%em_g_t_1 = grid%em_g_t_2
90 grid%em_g_t_2 = sqrt(grid%jcdfi_gama) * 1.0/grid%jcdfi_error_t * grid%em_a_t_2
91
92 grid%em_g_mu_1 = grid%em_g_mu_2
93 grid%em_g_mu_2 = sqrt(grid%jcdfi_gama) * 1.0/grid%jcdfi_error_mu * grid%em_a_mu_2
94
95 grid%g_moist_1 = grid%g_moist_2
96 grid%g_moist_2 = sqrt(grid%jcdfi_gama) * 1.0/grid%jcdfi_error_q * grid%a_moist_2
97
98 call med_hist_out ( grid , 3, config_flags )
99
100 ! jc = 0.0
101 ! jc = jc + 0.5*da_dot( grid%em_g_u_2 , grid%em_g_u_2 , size(grid%em_g_u_2) )
102 ! print*,'jc(u) =' , 0.5*da_dot( grid%em_g_u_2 , grid%em_g_u_2 , size(grid%em_g_u_2) )
103 ! jc = jc + 0.5*da_dot( grid%em_g_v_2 , grid%em_g_v_2 , size(grid%em_g_v_2) )
104 ! print*,'jc(v) =' , 0.5**da_dot( grid%em_g_v_2 , grid%em_g_v_2 , size(grid%em_g_v_2) )
105 ! jc = jc + 0.5*da_dot( grid%em_g_w_2 , grid%em_g_w_2 , size(grid%em_g_w_2) )
106 ! print*,'jc(w) =' , 0.5*da_dot( grid%em_g_w_2 , grid%em_g_w_2 , size(grid%em_g_w_2) )
107 ! jc = jc + 0.5*da_dot( grid%em_g_t_2 , grid%em_g_t_2 , size(grid%em_g_t_2) )
108 ! print*,'jc(t) =' , 0.5*da_dot( grid%em_g_t_2 , grid%em_g_t_2 , size(grid%em_g_t_2) )
109 ! jc = jc + 0.5*da_dot( grid%em_g_mu_2 , grid%em_g_mu_2 , size(grid%em_g_mu_2) )
110 ! print*,'jc(mu) =' , 0.5*da_dot( grid%em_g_mu_2 , grid%em_g_mu_2 , size(grid%em_g_mu_2) )
111 ! jc = jc + 0.5*da_dot( grid%em_g_ph_2 , grid%em_g_ph_2 , size(grid%em_g_ph_2) )
112 ! print*,'jc(ph) =' ,0.5*da_dot( grid%em_g_ph_2 , grid%em_g_ph_2 , size(grid%em_g_ph_2) )
113 ! jc = jc + 0.5*da_dot( grid%g_moist_2 , grid%g_moist_2 , size(grid%g_moist_2) )
114 ! print*,'jc(q) =' , 0.5*da_dot( grid%g_moist_2 , grid%g_moist_2 , size(grid%g_moist_2) )
115 ! print*,'Total jc =', jc
116
117
118 grid%em_g_u_2 = grid%em_g_u_1
119 grid%em_g_v_2 = grid%em_g_v_1
120 grid%em_g_t_2 = grid%em_g_t_1
121 grid%em_g_mu_2 = grid%em_g_mu_1
122 grid%g_moist_2 = grid%g_moist_1
123
124 END SUBROUTINE jcdfi_output_forcing
125 !
126 !
127 !
128 SUBROUTINE jcdfi_add_forcing ( grid ,i )
129
130 IMPLICIT NONE
131
132 TYPE ( domain ) , intent(inout) :: grid
133 integer , intent(in) :: i
134 integer :: jcdfi_flag
135
136 jcdfi_flag = 0
137 if ( grid%jcdfi_use ) jcdfi_flag = 1
138 grid%em_a_u_2 = grid%em_a_u_2 + jcdfi_flag * sqrt(grid%jcdfi_gama) * 1.0/grid%jcdfi_error_wind * h(i) * grid%em_g_u_1
139 grid%em_a_v_2 = grid%em_a_v_2 + jcdfi_flag * sqrt(grid%jcdfi_gama) * 1.0/grid%jcdfi_error_wind * h(i) * grid%em_g_v_1
140 grid%em_a_t_2 = grid%em_a_t_2 + jcdfi_flag * sqrt(grid%jcdfi_gama) * 1.0/grid%jcdfi_error_t * h(i) * grid%em_g_t_1
141 grid%em_a_mu_2 = grid%em_a_mu_2 + jcdfi_flag * sqrt(grid%jcdfi_gama) * 1.0/grid%jcdfi_error_mu * h(i) * grid%em_g_mu_1
142 grid%a_moist_2 = grid%a_moist_2 + jcdfi_flag * sqrt(grid%jcdfi_gama) * 1.0/grid%jcdfi_error_q * h(i) * grid%g_moist_1
143
144
145 END SUBROUTINE jcdfi_add_forcing
146 !
147 !
148 !
149 SUBROUTINE jcdfi_input_forcing ( grid ,config_flags )
150
151 IMPLICIT NONE
152
153 TYPE ( domain ) , intent(inout) :: grid
154 type(grid_config_rec_type) , intent(in) :: config_flags
155
156 call med_auxinput3dfi_in ( grid , config_flags )
157
158 grid%em_g_u_1 = grid%em_g_u_2
159 grid%em_g_v_1 = grid%em_g_v_2
160 grid%em_g_t_1 = grid%em_g_t_2
161 grid%em_g_mu_1 = grid%em_g_mu_2
162 grid%g_moist_1 = grid%g_moist_2
163
164
165 grid%em_g_u_2 = 0.0
166 grid%em_g_v_2 = 0.0
167 grid%em_g_t_2 = 0.0
168 grid%em_g_mu_2 = 0.0
169 grid%g_moist_2 = 0.0
170
171 END SUBROUTINE jcdfi_input_forcing
172 !
173 !
174 SUBROUTINE jcdfi_finalize
175 deallocate (H)
176 END SUBROUTINE jcdfi_finalize
177 !
178 SUBROUTINE JCDFI_COEF (MYPE,NSTEPS,DT,TAUC,NFILT)
179 !min---C
180 !min---C calculate filter weights with selected window.
181 !min---C
182 !min---C peter lynch and xiang-yu huang
183 !min---C
184 !min---C ref: see hamming, r.w., 1989: digital filters,
185 !min---C prentice-hall international. 3rd edition.
186 !min---C
187 !min---C input: nsteps - number of timesteps
188 !min---C forward or backward.
189 !min---C dt - time step in seconds.
190 !min---C tauc - cut-off period in seconds.
191 !min---C nfilt - indicator for selected filter.
192 !min---C
193 !min---C output: h - array(0:nsteps) with the
194 !min---C required filter weights
195 !min---C
196 !min---C------------------------------------------------------------
197 !min---C
198
199 implicit none
200
201 integer, intent(in) :: mype, nsteps, nfilt
202 real , intent(in) :: dt, tauc
203
204 !min--- Local data
205
206 integer :: n
207 real :: pi, omegac, x, window, deltat
208 real :: hh(0:nsteps)
209
210 allocate ( h(0:nsteps*2) )
211 !min---
212
213 pi=4*ATAN(1.)
214 deltat=ABS(dt)
215 !min---
216 !min---------------------------------------------------------------
217 !min---
218 !min--- windows are defined by a call and held in hh.
219 !min---
220 if ( nfilt .eq. -1) call debug (nsteps,h)
221
222 IF ( NFILT .EQ. 0 ) CALL UNIFORM (NSTEPS,HH)
223 IF ( NFILT .EQ. 1 ) CALL LANCZOS (NSTEPS,HH)
224 IF ( NFILT .EQ. 2 ) CALL HAMMING (NSTEPS,HH)
225 IF ( NFILT .EQ. 3 ) CALL BLACKMAN (NSTEPS,HH)
226 IF ( NFILT .EQ. 4 ) CALL KAISER (NSTEPS,HH)
227 IF ( NFILT .EQ. 5 ) CALL POTTER2 (NSTEPS,HH)
228 IF ( NFILT .EQ. 6 ) CALL DOLPHWIN (NSTEPS,HH)
229 !min---
230 !min---------------------------------------------------------------
231 !min---
232 IF ( NFILT .LE. 6 ) THEN ! sinc-windowed filters
233 !min---
234 !min--- calculate the cutoff frequency
235 OMEGAC = 2.*PI/TAUC
236 !min---
237 DO N=0,NSTEPS
238 WINDOW = HH(N)
239 IF ( N .EQ. 0 ) THEN
240 X = (OMEGAC*DELTAT/PI)
241 ELSE
242 X = SIN(N*OMEGAC*DELTAT)/(N*PI)
243 ENDIF
244 HH(N) = X*WINDOW
245 ENDDO
246
247 !min---
248 !min--- normalize the sums to be unity
249 CALL NORMLZ(HH,NSTEPS)
250
251 DO N=0,NSTEPS
252 H(NSTEPS+N) = HH(N)
253 H(NSTEPS-N) = HH(N)
254 ENDDO
255 ! IF(MYPE.EQ.0) THEN
256 ! DO N=0,NSTEPS*2
257 ! WRITE(6,'(A,I3,1X,F10.7)') 'DF COEF: N, H =',N,H(N)
258 ! ENDDO
259 ! ENDIF
260
261 ELSEIF ( NFILT .EQ. 7 ) THEN ! dolph filter
262
263 CALL DOLPH(MYPE,DT,TAUC,NSTEPS,H)
264
265 ELSEIF ( NFILT .EQ. 8 ) THEN ! 2nd order, 2nd type quick start filter (RHO)
266
267 CALL RHOFIL(DT,TAUC,2,NSTEPS*2,2,H,NSTEPS*2)
268 IF(MYPE.EQ.0) THEN
269 DO N=0,NSTEPS*2
270 WRITE(6,'(A,I3,1X,F10.7)') 'RHO COEF: N, H =',N,H(N)
271 ENDDO
272 ENDIF
273 ENDIF
274
275 DO N=0,NSTEPS*2
276 H(N) = -1.0*H(N)
277 IF (N.EQ.NSTEPS) H(N) = 1.0 + H(N)
278 ENDDO
279
280 RETURN
281 END subroutine jcdfi_coef
282
283 !min---
284 SUBROUTINE NORMLZ(HH,NMAX)
285 !min---
286 !min--- normalize the sum of hh to be unity
287 implicit none
288
289 integer, intent(in) :: nmax
290 real , dimension(0:nmax), intent(out) :: hh
291
292 !min--- local data
293 real :: sumhh
294 integer :: n
295
296 SUMHH = HH(0)
297 DO N=1,NMAX
298 SUMHH = SUMHH + 2*HH(N)
299 ENDDO
300 DO N=0,NMAX
301 HH(N) = HH(N)/SUMHH
302 ENDDO
303 !min---
304 RETURN
305 END subroutine normlz
306
307
308
309 subroutine debug(nsteps, ww)
310
311 implicit none
312
313 integer, intent(in) :: nsteps
314 real , dimension(0:nsteps), intent(out) :: ww
315 integer :: n
316
317 do n=0,nsteps
318 ww(n)=0
319 enddo
320
321 ww(int(nsteps/2))=1
322
323 return
324 end subroutine debug
325
326
327 SUBROUTINE UNIFORM(NSTEPS,WW)
328 !min---
329 !min--- define uniform or rectangular window function.
330 !min---
331 implicit none
332
333 integer, intent(in) :: nsteps
334 real , dimension(0:nsteps), intent(out) :: ww
335
336 integer :: n
337
338 DO N=0,NSTEPS
339 WW(N) = 1.
340 ENDDO
341 RETURN
342 END subroutine uniform
343
344
345
346
347 SUBROUTINE LANCZOS(NSTEPS,WW)
348 !min---
349 !min--- define (genaralised) lanczos window function.
350 !min---
351 implicit none
352
353 integer, parameter :: nmax = 1000
354 integer, intent(in) :: nsteps
355 real , dimension(0:nmax), intent(out) :: ww
356 integer :: n
357 real :: power, pi, w
358
359 !min---
360 !min--- (for the usual lanczos window, power = 1 )
361 POWER = 1
362 !min---
363 PI=4*ATAN(1.)
364 DO N=0,NSTEPS
365 IF ( N .EQ. 0 ) THEN
366 W = 1.0
367 ELSE
368 W = SIN(N*PI/(NSTEPS+1)) / ( N*PI/(NSTEPS+1))
369 ENDIF
370 WW(N) = W**POWER
371 ENDDO
372 RETURN
373 END subroutine lanczos
374
375
376
377
378 SUBROUTINE HAMMING(NSTEPS,WW)
379 !min---
380 !min--- define (genaralised) hamming window function.
381 !min---
382 implicit none
383
384 integer, intent(in) :: nsteps
385 real, dimension(0:nsteps) :: ww
386 integer :: n
387 real :: alpha, pi, w
388
389 !min---
390 !min--- (for the usual hamming window, alpha=0.54,
391 !min--- for the hann window, alpha=0.50).
392 ALPHA=0.54
393 !min---
394 PI=4*ATAN(1.)
395 DO N=0,NSTEPS
396 IF ( N .EQ. 0 ) THEN
397 W = 1.0
398 ELSE
399 W = ALPHA + (1-ALPHA)*COS(N*PI/(NSTEPS))
400 ENDIF
401 WW(N) = W
402 ENDDO
403 RETURN
404 END subroutine hamming
405
406
407
408
409 SUBROUTINE BLACKMAN(NSTEPS,WW)
410 !min---
411 !min--- define blackman window function.
412 !min---
413 implicit none
414
415 integer, intent(in) :: nsteps
416 real, dimension(0:nsteps) :: ww
417 integer :: n
418 !min---
419 real :: pi, w
420
421 PI=4*ATAN(1.)
422 DO N=0,NSTEPS
423 IF ( N .EQ. 0 ) THEN
424 W = 1.0
425 ELSE
426 W = 0.42 + 0.50*COS( N*PI/(NSTEPS)) &
427 + 0.08*COS(2*N*PI/(NSTEPS))
428 ENDIF
429 WW(N) = W
430 ENDDO
431 RETURN
432 END subroutine blackman
433
434
435
436
437
438 SUBROUTINE KAISER(NSTEPS,WW)
439 !min---
440 !min--- define kaiser window function.
441 !min---
442 implicit none
443
444 integer, intent(in) :: nsteps
445 real, dimension(0:nsteps) :: ww
446 integer :: n
447 real :: alpha, xi0a, xn, as
448 real, external :: bessi0
449
450 ALPHA = 1
451 !min---
452 XI0A = BESSI0(ALPHA)
453 DO N=0,NSTEPS
454 XN = N
455 AS = ALPHA*SQRT(1.-(XN/NSTEPS)**2)
456 WW(N) = BESSI0(AS) / XI0A
457 ENDDO
458 RETURN
459 END subroutine kaiser
460
461
462
463
464 SUBROUTINE POTTER2(NSTEPS,WW)
465 !min---
466 !min--- define potter window function.
467 !min--- modified (by me) to fall off over twice the range.
468 !min---
469 implicit none
470 integer, intent(in) :: nsteps
471 real, dimension(0:nsteps),intent(out) :: ww
472 integer :: n
473 real :: ck, sum, arg
474 !min--- local data
475
476 real, dimension(0:3) :: d
477 real :: pi
478 integer :: ip
479
480 d(0) = 0.35577019
481 d(1) = 0.2436983
482 d(2) = 0.07211497
483 d(3) = 0.00630165
484
485 PI=4*ATAN(1.)
486 !min---
487 CK = 1.0
488 DO N=0,NSTEPS
489 IF(N.EQ.NSTEPS) CK = 0.5
490 ARG = PI*FLOAT(N)/FLOAT(NSTEPS)
491 !min--- modification in next statement
492 ARG = ARG/2.
493 !min--- end of modification
494 SUM = D(0)
495 DO IP=1,3
496 SUM = SUM + 2.*D(IP)*COS(ARG*FLOAT(IP))
497 ENDDO
498 WW(N) = CK*SUM
499 ENDDO
500 RETURN
501 END SUBROUTINE POTTER2
502
503
504 SUBROUTINE DOLPHWIN(M,WINDOW)
505 !min---
506 !min--- calculation of dolph-chebyshev window or, for short,
507 !min--- dolph window, using the expression in the reference:
508 !min---
509 !min---- antoniou, andreas, 1993: digital filters: analysis,
510 !min---- design and applications. mcgraw-hill, inc., 689pp.
511 !min----
512 !min---- the dolph window is optimal in the following sense:
513 !min---- for a given main-lobe width, the stop-band attenuation
514 !min---- is minimal; for a given stop-band level, the main-lobe
515 !min---- width is minimal.
516 !min----
517 !min---- it is possible to specify either the ripple-ratio r
518 !min---- or the stop-band edge thetas.
519 !min----
520 implicit none
521 integer, intent(in) :: m
522 real, dimension(0:m),intent(out) :: window
523
524 !min---- local data
525 real, dimension(0:2*m) :: t
526 real, dimension(0:m) :: w, time
527 real :: pi, thetas, x0, term1, term2, rr, r, db, sum, arg
528 integer :: n, nm1, nt, i
529
530 !min----
531 PI = 4*ATAN(1.D0)
532 !min----
533 THETAS = 2*PI/(M)
534
535 !min----
536 N = 2*M+1
537 NM1 = N-1
538 X0 = 1/COS(THETAS/2)
539 !min----
540 TERM1 = (X0 + SQRT(X0**2-1))**(FLOAT(N-1))
541 TERM2 = (X0 - SQRT(X0**2-1))**(FLOAT(N-1))
542 RR = 0.5*(TERM1+TERM2)
543 R = 1/RR
544 DB = 20*LOG10(R)
545 WRITE(*,'(1X,''DOLPH: M,N='',2I8)')M,N
546 !CIVY--- CALL MESAGE(MYPE,NPRPE,MES)
547 WRITE(*,'(1X,''DOLPH: THETAS (STOP-BAND EDGE)='',F10.3)')THETAS
548 !CIVY--- CALL MESAGE(MYPE,NPRPE,MES)
549 WRITE(*,'(1X,''DOLPH: R,DB='',2F10.3)')R, DB
550 !min----
551 DO NT=0,M
552 SUM = RR
553 DO I=1,M
554 ARG = X0*COS(I*PI/N)
555 CALL CHEBY(T,NM1,ARG)
556 TERM1 = T(NM1)
557 TERM2 = COS(2*NT*PI*I/N)
558 SUM = SUM + 2*TERM1*TERM2
559 ENDDO
560 W(NT) = SUM/N
561 TIME(NT) = NT
562 ENDDO
563 !min----
564 !min------ fill up the array for return
565 DO NT=0,M
566 WINDOW(NT) = W(NT)
567 ENDDO
568 !min----
569 RETURN
570 END subroutine dolphwin
571 !min----
572 !min----<A NAME="DOLPH">
573 !min----
574 SUBROUTINE DOLPH(MYPE, DELTAT, TAUS, M, WINDOW)
575 !min----
576 !min---- calculation of dolph-chebyshev window or, for short,
577 !min---- dolph window, using the expression in the reference:
578 !min----
579 !min---- antoniou, andreas, 1993: digital filters: analysis,
580 !min---- design and applications. mcgraw-hill, inc., 689pp.
581 !min----
582 !min---- the dolph window is optimal in the following sense:
583 !min---- for a given main-lobe width, the stop-band attenuation
584 !min---- is minimal; for a given stop-band level, the main-lobe
585 !min---- width is minimal.
586 !min----
587 !min---- mype ! current pe number (used only for message)
588 !min----
589 !min----------------------------------------------------------------
590 !min----
591 implicit none
592
593 integer, intent(in) :: m
594 real, dimension(0:m),intent(out) :: window
595 integer, intent(in) :: mype
596 real, intent(in) :: deltat, taus
597
598
599
600 !min---- local data
601 integer, PARAMETER :: NMAX = 5000
602 REAL, dimension(0:NMAX) :: t, w, time
603 real, dimension(0:2*nmax) :: w2
604 INTEGER :: NPRPE=0 ! no of pe
605 CHARACTER*80 :: MES
606 !min----
607 real :: pi, thetas, x0, term1, term2, rr, r,db, sum, arg, sumw
608 integer :: n, nm1, i, nt
609
610 PI = 4*ATAN(1.D0)
611
612 print *, 'in dfcoef, deltat = ', deltat, 'taus=',taus
613
614
615 !min----
616 N = 2*M+1
617 NM1 = N-1
618 !min----
619 THETAS = 2*PI*ABS(DELTAT/TAUS)
620 X0 = 1/COS(THETAS/2)
621 TERM1 = (X0 + SQRT(X0**2-1))**(FLOAT(N-1))
622 TERM2 = (X0 - SQRT(X0**2-1))**(FLOAT(N-1))
623 RR = 0.5*(TERM1+TERM2)
624 R = 1/RR
625 DB = 20*LOG10(R)
626
627
628 WRITE(*,'(1X,''DOLPH: M,N='',2I8)')M,N
629 !CIVY--- CALL MESAGE(MYPE,NPRPE,MES)
630 WRITE(*,'(1X,''DOLPH: THETAS (STOP-BAND EDGE)='',F10.3)')THETAS
631 !CIVY--- CALL MESAGE(MYPE,NPRPE,MES)
632 WRITE(*,'(1X,''DOLPH: R,DB='',2F10.3)')R, DB
633 !CIVY--- CALL MESAGE(MYPE,NPRPE,MES)
634 !min----
635 !min----------------------------------------------------------------
636 !min----
637 DO NT=0,M
638 SUM = 1
639 DO I=1,M
640 ARG = X0*COS(I*PI/N)
641 CALL CHEBY(T,NM1,ARG)
642 TERM1 = T(NM1)
643 TERM2 = COS(2*NT*PI*I/N)
644 SUM = SUM + R*2*TERM1*TERM2
645 ENDDO
646 W(NT) = SUM/N
647 TIME(NT) = NT
648 WRITE(*,'(1X,''DOLPH: TIME, W='',F10.6,2X,E17.7)') &
649 TIME(NT), W(NT)
650 !CIVY--- CALL MESAGE(MYPE,NPRPE,MES)
651 ENDDO
652 !min---- fill in the negative-time values by symmetry.
653 DO NT=0,M
654 W2(M+NT) = W(NT)
655 W2(M-NT) = W(NT)
656 ENDDO
657 !min----
658 !min------ fill up the array for return
659 SUMW = 0.
660 DO NT=0,2*M
661 SUMW = SUMW + W2(NT)
662 ENDDO
663 WRITE(*,'(1X,''DOLPH: SUM OF WEIGHTS W2='',F10.4)')SUMW
664 !CIVY--- CALL MESAGE(MYPE,NPRPE,MES)
665 !min----
666 DO NT=0,2*M
667 WINDOW(NT) = W2(NT)
668 ENDDO
669 !min----
670 RETURN
671 END subroutine dolph
672 !min----
673 !min----<A NAME="CHEBY">
674 !min----
675 SUBROUTINE CHEBY(T,N,X)
676 !min----
677 !min---- calculate all chebyshev polynomials up to order n
678 !min---- for the argument value x.
679 !min----
680 !min---- reference: numerical recipes, page 184, recurrence
681 !min---- t_n(x) = 2xt_{n-1}(x) - t_{n-2}(x) , n>=2.
682 !min----
683 implicit none
684
685 integer, intent(in) :: n
686 real , intent(in) :: x
687 real, dimension(0:n):: t
688
689 integer :: nn
690 !min----
691 T(0) = 1
692 T(1) = X
693 IF(N.LT.2) RETURN
694 DO NN=2,N
695 T(NN) = 2*X*T(NN-1) - T(NN-2)
696 ENDDO
697 !min----
698 RETURN
699 END subroutine cheby
700 !min----
701 !min-<A NAME="RHO">
702 !min-
703
704
705
706 SUBROUTINE RHOFIL(DT,TAUC,NORDER,NSTEP,ICTYPE,FROW,NOSIZE)
707 !min---c
708 !min---c RHO = recurssive high order.
709 !min---c
710 !min---- This routine calculates and returns the
711 !min---- Last Row, FROW, of the FILTER matrix.
712 !min----
713 !min---- Input Parameters:
714 !min---- DT : Time Step in seconds
715 !min---- TAUC : Cut-off period (hours)
716 !min---- NORDER : Order of QS Filter
717 !min---- NSTEP : Number of step/Size of row.
718 !min---- ICTYPE : Initial Conditions
719 !min---- NOSIZE : Max. side of FROW.
720 !min----
721 !min---- Working Fields:
722 !min---- ACOEF : X-coefficients of filter
723 !min---- BCOEF : Y-coefficients of filter
724 !min---- FILTER : Filter Matrix.
725 !min----
726 !min---- Output Parameters:
727 !min---- FROW : Last Row of Filter Matrix.
728 !min----
729 !min---- Note: Two types of initial conditions are permitted.
730 !min---- ICTYPE = 1 : Order increasing each row to NORDER.
731 !min---- ICTYPE = 2 : Order fixed at NORDER throughout.
732 !min----
733 !min---- DOUBLE PRECISION USED THROUGHOUT.
734 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
735 !min----
736 DOUBLE PRECISION MUC
737 !min----
738 !min---- N.B. Single Precision for List Parameters.
739 REAL, intent(in) :: DT,TAUC
740 !min---- Space for the last row of FILTER.
741 integer, intent(in) :: norder, nstep, ictype, nosize
742 REAL , dimension(0:nosize), intent(out):: FROW
743 !min----
744 !min---- Arrays for rho filter.
745 integer, PARAMETER :: NOMAX=100
746 real , dimension(0:NOMAX) :: acoef, bcoef
747 real , dimension(0:NOMAX,0:NOMAX) :: filter
748 !min---- Working space.
749 real , dimension(0:NOMAX) :: alpha, beta
750 !kjw----
751 real :: DTT
752 integer :: kk, k, NN, NORD, NR, NC, NROW, LL
753 COMPLEX :: IOTA
754 !min----
755 !min-------------------------------------------------------------------
756 DTT = ABS(DT)
757 ! DT = ABS(DT)
758 PI = 2*DASIN(1.D0)
759 IOTA = CMPLX(0.,1.)
760 !min-------------------------------------------------------------------
761 !min----
762 !min--- Filtering Parameters (derived).
763 THETAC = 2*PI*DTT/(TAUC)
764 ! THETAC = 2*PI*DT/(TAUC)
765 MUC = tan(THETAC/2)
766 FC = THETAC/(2*PI)
767 !min---
768 !min---------------------------------------------------------------
769 !min---
770 !min--- Clear the arrays.
771 DO NC=0,NOMAX
772 ACOEF(NC) = 0.
773 BCOEF(NC) = 0.
774 ALPHA(NC) = 0.
775 BETA (NC) = 0.
776 FROW (NC) = 0.
777 DO NR=0,NOMAX
778 FILTER(NR,NC) = 0.
779 ENDDO
780 ENDDO
781 !min---
782 !min---- Fill up the Filter Matrix.
783 FILTER(0,0) = 1.
784 !min---
785 !min--- Get the coefficients of the Filter.
786 IF ( ICTYPE.eq.2 ) THEN
787 CALL RHOCOF(NORDER,NOMAX,MUC, ACOEF,BCOEF)
788 ENDIF
789 !min----
790 DO 100 NROW=1,NSTEP
791 !min----
792 IF ( ICTYPE.eq.1 ) THEN
793 NORD = MIN(NROW,NORDER)
794 IF ( NORD.le.NORDER) THEN
795 CALL RHOCOF(NORD,NOMAX,MUC, ACOEF,BCOEF)
796 ENDIF
797 ENDIF
798 !min---
799 DO K=0,NROW
800 ALPHA(K) = ACOEF(NROW-K)
801 IF(K.lt.NROW) BETA(K) = BCOEF(NROW-K)
802 ENDDO
803 !min----
804 !min---- Correction for terms of negative index.
805 IF ( ICTYPE.eq.2 ) THEN
806 IF ( NROW.lt.NORDER ) THEN
807 CN = 0.
808 DO NN=NROW+1,NORDER
809 CN = CN + (ACOEF(NN)+BCOEF(NN))
810 ENDDO
811 !min----CC ACOEF(NROW) = ACOEF(NROW) + CN
812 ALPHA(0) = ALPHA(0) + CN
813 ENDIF
814 ENDIF
815 !min----
816 !min---- Check sum of ALPHAs and BETAs = 1
817 SUMAB = 0.
818 DO NN=0,NROW
819 SUMAB = SUMAB + ALPHA(NN)
820 IF(NN.lt.NROW) SUMAB = SUMAB + BETA(NN)
821 ENDDO
822 !min---
823 DO KK=0,NROW-1
824 SUMBF = 0.
825 DO LL=0,NROW-1
826 SUMBF = SUMBF + BETA(LL)*FILTER(LL,KK)
827 ENDDO
828 FILTER(NROW,KK) = ALPHA(KK)+SUMBF
829 ENDDO
830 FILTER(NROW,NROW) = ALPHA(NROW)
831 !min---
832 !min---- Check sum of row elements = 1
833 SUMROW = 0.
834 DO NN=0,NROW
835 SUMROW = SUMROW + FILTER(NROW,NN)
836 ENDDO
837 !min----
838 100 CONTINUE
839 !min----
840 DO NC=0,NSTEP
841 FROW(NC) = FILTER(NSTEP,NC)
842 ENDDO
843 !min----
844 RETURN
845 END subroutine RHOFIL
846 !min----
847 !min---***************************************************************
848 !min---
849 SUBROUTINE RHOCOF(NORD,NOMAX,MUC,CA,CB)
850 !min---
851 !min--- Get the coefficients of the RHO Filter
852 !min---
853 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
854 !min---
855 integer, intent(in) :: nord, nomax
856 real, dimension(0:nomax) :: ca, cb
857 integer :: NN
858 !min---
859 COMPLEX :: IOTA
860 !min--- REAL MUC
861 !min--- COMPLEX ZN
862 DOUBLE PRECISION :: MUC, ZN
863 !min---
864 DOUBLE PRECISION :: CNR
865 !min---
866 PI = 2*ASIN(1.)
867 ROOT2 = SQRT(2.)
868 IOTA = (0.,1.)
869 !min---
870 RN = 1./FLOAT(NORD)
871 SIGMA = 1./( SQRT(2.**RN-1.) )
872 !min---
873 GAIN = (MUC*SIGMA/(1+MUC*SIGMA))**NORD
874 ZN = (1-MUC*SIGMA)/(1+MUC*SIGMA)
875 !min---
876 DO NN=0,NORD
877 CA(NN) = CNR(NORD,NN)*GAIN
878 IF(NN.gt.0) CB(NN) = -CNR(NORD,NN)*(-ZN)**NN
879 ENDDO
880 !min---C Check sum of coefficients = 1
881 SUMCOF = 0.
882 DO NN=0,NORD
883 SUMCOF = SUMCOF + CA(NN)
884 IF(NN.gt.0) SUMCOF = SUMCOF + CB(NN)
885 ENDDO
886 !min---
887 RETURN
888 END subroutine RHOCOF
889
890 #endif
891 END MODULE module_dfi
892 !min---***************************************************************
893 DOUBLE PRECISION FUNCTION CNR(N,R)
894 !min---
895 !min--- Binomial Coefficient (n,r).
896 IMPLICIT DOUBLE PRECISION(C,X)
897 !min---
898 INTEGER , intent(in) :: n, R
899 IF ( R.eq.0 ) THEN
900 CNR = 1.0
901 RETURN
902 ENDIF
903 Coeff = 1.0
904 XN = DFLOAT(N)
905 XR = DFLOAT(R)
906 DO K=1,R
907 XK = DFLOAT(K)
908 COEFF = COEFF * ( (XN-XR+XK)/XK )
909 ENDDO
910 CNR = COEFF
911 RETURN
912 END FUNCTION
913 !min---***************************************************************
914
915
916
917 real FUNCTION BESSI0(X)
918 !min---
919 !min--- from numerical recipes (press, et al.)
920 !min---
921 implicit none
922
923 real(8) :: Y
924 real(8) :: P1 = 1.0d0
925 real(8) :: P2 = 3.5156230D0
926 real(8) :: P3 = 3.0899424D0
927 real(8) :: P4 = 1.2067492D0
928 real(8) :: P5 = 0.2659732D0
929 real(8) :: P6 = 0.360768D-1
930 real(8) :: P7 = 0.45813D-2
931
932 real*8 :: Q1 = 0.39894228D0
933 real*8 :: Q2 = 0.1328592D-1
934 real*8 :: Q3 = 0.225319D-2
935 real*8 :: Q4 = -0.157565D-2
936 real*8 :: Q5 = 0.916281D-2
937 real*8 :: Q6 = -0.2057706D-1
938 real*8 :: Q7 = 0.2635537D-1
939 real*8 :: Q8 = -0.1647633D-1
940 real*8 :: Q9 = 0.392377D-2
941
942 real :: x, ax
943
944
945 IF (ABS(X).LT.3.75) THEN
946 Y=(X/3.75)**2
947 BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
948 ELSE
949 AX=ABS(X)
950 Y=3.75/AX
951 BESSI0=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4 &
952 +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
953 ENDIF
954 RETURN
955 END function bessi0
956
957 real function da_dot(x,y,n)
958 !
959 ! forms the dot product of two vectors.
960 ! uses unrolled loops for increments equal to one.
961 ! jack dongarra, linpack, 3/11/78.
962 !
963 implicit none
964
965 integer, intent(in) :: n
966 real, dimension(n), intent(in) :: x,y
967
968 real dtemp1,dtemp1x
969 integer i,m,mp1
970 integer ierror ! MPI error code.
971 !
972 da_dot = 0.0
973 if(n <= 0)return
974 dtemp1 = 0.0
975 !
976 ! code for both increments equal to 1
977 !
978 if(n > 0 ) then
979 m = mod(n,5)
980 if( m /= 0 ) then
981 do i = 1,m
982 dtemp1 = dtemp1 + x(i)*y(i)
983 end do
984 end if
985 if( n >= 5 ) then
986 mp1 = m + 1
987 do i = mp1,n,5
988 dtemp1 = dtemp1 + x(i )*y(i ) + x(i + 1)*y(i + 1) + &
989 x(i + 2)*y(i + 2) + x(i + 3)*y(i + 3) + &
990 x(i + 4)*y(i + 4)
991 end do
992 end if
993 end if
994 da_dot = dtemp1
995 end function da_dot
996