wrfv2_fire/share/dfi.F

   1    SUBROUTINE dfi_accumulate( grid )
   2
   3       USE module_domain
   4       USE module_configure
   5       USE module_driver_constants
   6       USE module_machine
   7       USE module_dm
   8       USE module_model_constants
   9       USE module_state_description
  10
  11       IMPLICIT NONE
  12
  13       REAL hn
  14
  15       !  Input data.
  16
  17       TYPE(domain) , POINTER          :: grid
  18
  19 #if (EM_CORE == 1)
  20
  21       IF ( grid%dfi_opt .EQ. DFI_NODFI .OR. grid%dfi_stage .EQ. DFI_FST ) RETURN
  22
  23       hn = grid%hcoeff(grid%itimestep+1)
  24
  25       ! accumulate dynamic variables
  26        grid%dfi_mu(:,:)    = grid%dfi_mu(:,:)    + grid%mu_2(:,:)    * hn
  27        grid%dfi_u(:,:,:)   = grid%dfi_u(:,:,:)   + grid%u_2(:,:,:)   * hn
  28        grid%dfi_v(:,:,:)   = grid%dfi_v(:,:,:)   + grid%v_2(:,:,:)   * hn
  29        grid%dfi_w(:,:,:)   = grid%dfi_w(:,:,:)   + grid%w_2(:,:,:)   * hn
  30        grid%dfi_ww(:,:,:)  = grid%dfi_ww(:,:,:)  + grid%ww(:,:,:)    * hn
  31        grid%dfi_t(:,:,:)   = grid%dfi_t(:,:,:)   + grid%t_2(:,:,:)   * hn
  32        grid%dfi_phb(:,:,:) = grid%dfi_phb(:,:,:) + grid%phb(:,:,:)   * hn
  33        grid%dfi_ph0(:,:,:) = grid%dfi_ph0(:,:,:) + grid%ph0(:,:,:)   * hn
  34        grid%dfi_php(:,:,:) = grid%dfi_php(:,:,:) + grid%php(:,:,:)   * hn
  35        grid%dfi_p(:,:,:)   = grid%dfi_p(:,:,:)   + grid%p(:,:,:)     * hn
  36        grid%dfi_ph(:,:,:)  = grid%dfi_ph(:,:,:)  + grid%ph_2(:,:,:)  * hn
  37        grid%dfi_tke(:,:,:) = grid%dfi_tke(:,:,:) + grid%tke_2(:,:,:) * hn
  38        grid%dfi_al(:,:,:)  = grid%dfi_al(:,:,:)  + grid%al(:,:,:)    * hn
  39        grid%dfi_alt(:,:,:) = grid%dfi_alt(:,:,:) + grid%alt(:,:,:)   * hn
  40        grid%dfi_pb(:,:,:)  = grid%dfi_pb(:,:,:)  + grid%pb(:,:,:)    * hn
  41        grid%dfi_moist(:,:,:,:)  = grid%dfi_moist(:,:,:,:)  + grid%moist(:,:,:,:)  * hn
  42        grid%dfi_scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) + grid%scalar(:,:,:,:) * hn
  43
  44       ! accumulate DFI coefficient
  45       grid%hcoeff_tot = grid%hcoeff_tot + hn
  46 #endif
  47
  48    END SUBROUTINE dfi_accumulate
  49
  50
  51 #if (EM_CORE == 1)
  52    SUBROUTINE wrf_dfi_bck_init ( )
  53
  54       USE module_domain
  55       USE module_utility
  56
  57       IMPLICIT NONE
  58
  59       INTEGER rc
  60
  61       INTERFACE
  62          SUBROUTINE Setup_Timekeeping(grid)
  63             USE module_domain
  64             TYPE (domain), POINTER :: grid
  65          END SUBROUTINE Setup_Timekeeping
  66
  67          SUBROUTINE dfi_save_arrays(grid)
  68             USE module_domain
  69             TYPE (domain), POINTER :: grid
  70          END SUBROUTINE dfi_save_arrays
  71
  72          SUBROUTINE dfi_clear_accumulation(grid)
  73             USE module_domain
  74             TYPE (domain), POINTER :: grid
  75          END SUBROUTINE dfi_clear_accumulation
  76
  77          SUBROUTINE optfil_driver(grid)
  78             USE module_domain
  79             TYPE (domain), POINTER :: grid
  80          END SUBROUTINE optfil_driver
  81
  82          SUBROUTINE start_domain(grid, allowed_to_read)
  83             USE module_domain
  84             TYPE (domain)          :: grid
  85             LOGICAL, INTENT(IN)    :: allowed_to_read
  86          END SUBROUTINE start_domain
  87       END INTERFACE
  88
  89       head_grid%dfi_stage = DFI_BCK
  90
  91       ! Negate time step
  92       CALL nl_set_time_step ( 1, -head_grid%time_step )
  93
  94       CALL Setup_Timekeeping (head_grid)
  95
  96       ! set physics options to zero
  97       CALL nl_set_mp_physics( 1, 0 )
  98       CALL nl_set_ra_lw_physics( 1, 0 )
  99       CALL nl_set_ra_sw_physics( 1, 0 )
 100       CALL nl_set_sf_surface_physics( 1, 0 )
 101       CALL nl_set_sf_sfclay_physics( 1, 0 )
 102       CALL nl_set_bl_pbl_physics( 1, 0 )
 103       CALL nl_set_cu_physics( 1, 0 )
 104
 105       ! set diffusion to zero for backward integration
 106       CALL nl_set_km_opt( 1, head_grid%km_opt_dfi)
 107       CALL nl_set_pd_moist( 1, head_grid%pd_moist_dfi)
 108
 109       head_grid%start_subtime = domain_get_start_time ( head_grid )
 110       head_grid%stop_subtime = domain_get_stop_time ( head_grid )
 111
 112       CALL WRFU_ClockSet(head_grid%domain_clock, currTime=head_grid%start_subtime, rc=rc)
 113
 114       CALL dfi_save_arrays ( head_grid )
 115       CALL dfi_clear_accumulation( head_grid )
 116       CALL optfil_driver(head_grid)
 117
 118       !tgs need to call start_domain here to reset bc initialization for negative dt
 119       CALL start_domain ( head_grid , .TRUE. )
 120
 121    END SUBROUTINE wrf_dfi_bck_init
 122
 123
 124    SUBROUTINE wrf_dfi_fwd_init ( )
 125
 126       USE module_domain
 127       USE module_utility
 128
 129       IMPLICIT NONE
 130
 131       INTEGER rc
 132
 133       INTERFACE
 134          SUBROUTINE Setup_Timekeeping(grid)
 135             USE module_domain
 136             TYPE (domain), POINTER :: grid
 137          END SUBROUTINE Setup_Timekeeping
 138
 139          SUBROUTINE dfi_save_arrays(grid)
 140             USE module_domain
 141             TYPE (domain), POINTER :: grid
 142          END SUBROUTINE dfi_save_arrays
 143
 144          SUBROUTINE dfi_clear_accumulation(grid)
 145             USE module_domain
 146             TYPE (domain), POINTER :: grid
 147          END SUBROUTINE dfi_clear_accumulation
 148
 149          SUBROUTINE optfil_driver(grid)
 150             USE module_domain
 151             TYPE (domain), POINTER :: grid
 152          END SUBROUTINE optfil_driver
 153
 154          SUBROUTINE start_domain(grid, allowed_to_read)
 155             USE module_domain
 156             TYPE (domain)          :: grid
 157             LOGICAL, INTENT(IN)    :: allowed_to_read
 158          END SUBROUTINE start_domain
 159       END INTERFACE
 160
 161       head_grid%dfi_stage = DFI_FWD
 162
 163       ! get the negative time step from the namelist and store
 164       ! it as positive again.
 165       ! for Setup_Timekeeping to use when setting the clock
 166       ! note that this ignores fractional parts of time step
 167       CALL nl_get_time_step( 1, head_grid%time_step )
 168       head_grid%time_step = abs(head_grid%time_step)
 169       CALL nl_set_time_step( 1, head_grid%time_step )
 170
 171       head_grid%itimestep=0
 172       head_grid%xtime=0.
 173
 174       ! reset physics options to normal
 175       CALL nl_set_mp_physics( 1, head_grid%mp_physics)
 176       CALL nl_set_ra_lw_physics( 1, head_grid%ra_lw_physics)
 177       CALL nl_set_ra_sw_physics( 1, head_grid%ra_sw_physics)
 178       CALL nl_set_sf_surface_physics( 1, head_grid%sf_surface_physics)
 179       CALL nl_set_sf_sfclay_physics( 1, head_grid%sf_sfclay_physics)
 180       CALL nl_set_bl_pbl_physics( 1, head_grid%bl_pbl_physics)
 181       CALL nl_set_cu_physics( 1, head_grid%cu_physics)
 182
 183       ! reset km_opt to norlmal
 184       CALL nl_set_km_opt( 1, head_grid%km_opt)
 185       CALL nl_set_pd_moist( 1, head_grid%pd_moist)
 186
 187       CALL Setup_Timekeeping (head_grid)
 188       head_grid%start_subtime = domain_get_start_time ( head_grid )
 189       head_grid%stop_subtime = domain_get_stop_time ( head_grid )
 190
 191       CALL WRFU_ClockSet(head_grid%domain_clock, currTime=head_grid%start_subtime, rc=rc)
 192
 193       IF ( head_grid%dfi_opt .EQ. DFI_DFL ) THEN
 194          CALL dfi_save_arrays ( head_grid )
 195       END IF
 196       CALL dfi_clear_accumulation( head_grid )
 197       CALL optfil_driver(head_grid)
 198
 199       !tgs need to call it here to reset bc initialization for positive time_step
 200       CALL start_domain ( head_grid , .TRUE. )
 201
 202    END SUBROUTINE wrf_dfi_fwd_init
 203
 204
 205    SUBROUTINE wrf_dfi_fst_init ( )
 206
 207       USE module_domain
 208
 209       IMPLICIT NONE
 210
 211       CHARACTER (LEN=80)      :: wrf_error_message
 212
 213       INTERFACE
 214          SUBROUTINE Setup_Timekeeping(grid)
 215             USE module_domain
 216             TYPE (domain), POINTER :: grid
 217          END SUBROUTINE Setup_Timekeeping
 218
 219          SUBROUTINE dfi_save_arrays(grid)
 220             USE module_domain
 221             TYPE (domain), POINTER :: grid
 222          END SUBROUTINE dfi_save_arrays
 223
 224          SUBROUTINE dfi_clear_accumulation(grid)
 225             USE module_domain
 226             TYPE (domain), POINTER :: grid
 227          END SUBROUTINE dfi_clear_accumulation
 228
 229          SUBROUTINE optfil_driver(grid)
 230             USE module_domain
 231             TYPE (domain), POINTER :: grid
 232          END SUBROUTINE optfil_driver
 233
 234          SUBROUTINE start_domain(grid, allowed_to_read)
 235             USE module_domain
 236             TYPE (domain)          :: grid
 237             LOGICAL, INTENT(IN)    :: allowed_to_read
 238          END SUBROUTINE start_domain
 239       END INTERFACE
 240
 241       head_grid%dfi_stage = DFI_FST
 242
 243       head_grid%itimestep=0
 244       head_grid%xtime=0.         ! BUG: This will probably not work for all DFI options
 245
 246       CALL Setup_Timekeeping (head_grid)
 247       head_grid%start_subtime = domain_get_start_time ( head_grid )
 248       head_grid%stop_subtime = domain_get_stop_time ( head_grid )
 249
 250       CALL start_domain ( head_grid , .TRUE. )
 251
 252    END SUBROUTINE wrf_dfi_fst_init
 253
 254
 255    SUBROUTINE wrf_dfi_write_initialized_state( )
 256
 257       ! Driver layer
 258       USE module_domain
 259       USE module_io_domain
 260       USE module_configure
 261
 262       IMPLICIT NONE
 263
 264       INTEGER             :: fid, ierr
 265       CHARACTER (LEN=80)  :: wrf_error_message
 266       CHARACTER (LEN=80)  :: rstname
 267       CHARACTER (LEN=132) :: message
 268
 269       TYPE (grid_config_rec_type) :: config_flags
 270
 271       CALL model_to_grid_config_rec ( head_grid%id , model_config_rec , config_flags )
 272
 273       WRITE (wrf_err_message,'(A,I4)') 'Writing out initialized model state'
 274       CALL wrf_message(TRIM(wrf_err_message))
 275
 276       rstname = 'wrfinput_initialized_d01'
 277       CALL open_w_dataset ( fid, TRIM(rstname), head_grid, config_flags, output_model_input, "DATASET=INPUT", ierr )
 278       IF ( ierr .NE. 0 ) THEN
 279          WRITE( message , '("program wrf: error opening ",A," for writing")') TRIM(rstname)
 280          CALL WRF_ERROR_FATAL ( message )
 281       END IF
 282       CALL output_model_input ( fid,   head_grid, config_flags, ierr )
 283       CALL close_dataset ( fid, config_flags, "DATASET=INPUT" )
 284
 285    END SUBROUTINE wrf_dfi_write_initialized_state
 286
 287
 288    SUBROUTINE wrf_dfi_array_reset ( )
 289
 290       USE module_domain
 291
 292       IMPLICIT NONE
 293
 294       INTERFACE
 295          SUBROUTINE dfi_array_reset(grid)
 296             USE module_domain
 297             TYPE (domain), POINTER :: grid
 298          END SUBROUTINE dfi_array_reset
 299       END INTERFACE
 300
 301       ! Copy filtered arrays back into state arrays in grid structure, and
 302       !   restore original land surface fields
 303       CALL dfi_array_reset( head_grid )
 304
 305    END SUBROUTINE wrf_dfi_array_reset
 306
 307
 308    SUBROUTINE optfil_driver( grid )
 309
 310       USE module_domain
 311       USE module_wrf_error
 312       USE module_timing
 313       USE module_date_time
 314       USE module_configure
 315       USE module_state_description
 316       USE module_model_constants
 317
 318       IMPLICIT NONE
 319
 320       TYPE (domain) , POINTER                 :: grid
 321
 322       ! Local variables
 323       integer :: nstep2, nstepmax, rundfi, i, rc
 324       real :: timestep, tauc
 325       TYPE(WRFU_TimeInterval) :: run_interval
 326
 327       timestep=abs(grid%dt)
 328       run_interval = grid%stop_subtime - grid%start_subtime
 329
 330       CALL WRFU_TimeIntervalGet( run_interval, S=rundfi, rc=rc )
 331       rundfi = abs(rundfi)
 332
 333       nstep2= ceiling((1.0 + real(rundfi)/timestep) / 2.0)
 334
 335       ! nstep2 is equal to a half of timesteps per initialization period,
 336       ! should not exceed nstepmax
 337
 338       tauc = real(grid%dfi_cutoff_seconds)
 339
 340       ! Get DFI coefficient
 341       grid%hcoeff(:) = 0.0
 342       IF ( grid%dfi_nfilter < 0 .OR. grid%dfi_nfilter > 8 ) THEN
 343 write(0,*) 'Invalid filter specified in namelist.'
 344       ELSE
 345          call dfcoef(nstep2-1, grid%dt, tauc, grid%dfi_nfilter, grid%hcoeff)
 346       END IF
 347
 348       IF ( MOD(int(1.0 + real(rundfi)/timestep),2) /= 0 ) THEN
 349          DO i=1,nstep2-1
 350             grid%hcoeff(2*nstep2-i) = grid%hcoeff(i)
 351          END DO
 352       ELSE
 353          DO i=1,nstep2
 354             grid%hcoeff(2*nstep2-i+1) = grid%hcoeff(i)
 355          END DO
 356       END IF
 357
 358    END SUBROUTINE optfil_driver
 359
 360
 361    SUBROUTINE dfi_clear_accumulation( grid )
 362
 363       USE module_domain
 364       USE module_configure
 365       USE module_driver_constants
 366       USE module_machine
 367       USE module_dm
 368       USE module_model_constants
 369       USE module_state_description
 370
 371       IMPLICIT NONE
 372
 373       !  Input data.
 374       TYPE(domain) , POINTER          :: grid
 375
 376       grid%dfi_mu(:,:)       = 0.
 377       grid%dfi_u(:,:,:)      = 0.
 378       grid%dfi_v(:,:,:)      = 0.
 379       grid%dfi_w(:,:,:)      = 0.
 380       grid%dfi_ww(:,:,:)     = 0.
 381       grid%dfi_t(:,:,:)      = 0.
 382       grid%dfi_phb(:,:,:)    = 0.
 383       grid%dfi_ph0(:,:,:)    = 0.
 384       grid%dfi_php(:,:,:)    = 0.
 385       grid%dfi_p(:,:,:)      = 0.
 386       grid%dfi_ph(:,:,:)     = 0.
 387       grid%dfi_tke(:,:,:)    = 0.
 388       grid%dfi_al(:,:,:)     = 0.
 389       grid%dfi_alt(:,:,:)    = 0.
 390       grid%dfi_pb(:,:,:)     = 0.
 391       grid%dfi_moist(:,:,:,:)  = 0.
 392       grid%dfi_scalar(:,:,:,:) = 0.
 393
 394       grid%hcoeff_tot  = 0.0
 395
 396    END SUBROUTINE dfi_clear_accumulation
 397
 398
 399    SUBROUTINE dfi_save_arrays( grid )
 400
 401       USE module_domain
 402       USE module_configure
 403       USE module_driver_constants
 404       USE module_machine
 405       USE module_dm
 406       USE module_model_constants
 407       USE module_state_description
 408
 409       IMPLICIT NONE
 410
 411       !  Input data.
 412       TYPE(domain) , POINTER          :: grid
 413
 414       ! save surface 2-D fields
 415       grid%dfi_SNOW(:,:)   = grid%SNOW(:,:)
 416       grid%dfi_SNOWH(:,:)  = grid%SNOWH(:,:)
 417       grid%dfi_SNOWC(:,:)  = grid%SNOWC(:,:)
 418       grid%dfi_CANWAT(:,:) = grid%CANWAT(:,:)
 419       grid%dfi_TSK(:,:)    = grid%TSK(:,:)
 420       grid%dfi_QVG(:,:)    = grid%QVG(:,:)
 421       grid%dfi_SOILT1(:,:) = grid%SOILT1(:,:)
 422       grid%dfi_TSNAV(:,:)  = grid%TSNAV(:,:)
 423
 424       ! save soil fields
 425       grid%dfi_TSLB(:,:,:)         = grid%TSLB(:,:,:)
 426       grid%dfi_SMOIS(:,:,:)        = grid%SMOIS(:,:,:)
 427       grid%dfi_SMFR3D(:,:,:)       = grid%SMFR3D(:,:,:)
 428       grid%dfi_KEEPFR3DFLAG(:,:,:) = grid%KEEPFR3DFLAG(:,:,:)
 429
 430    END SUBROUTINE dfi_save_arrays
 431
 432
 433    SUBROUTINE dfi_array_reset( grid )
 434
 435       USE module_domain
 436       USE module_configure
 437       USE module_driver_constants
 438       USE module_machine
 439       USE module_dm
 440       USE module_model_constants
 441       USE module_state_description
 442
 443       IMPLICIT NONE
 444
 445       !  Input data.
 446       TYPE(domain) , POINTER          :: grid
 447
 448       IF ( grid%dfi_opt .EQ. DFI_NODFI ) RETURN
 449
 450
 451       ! Set dynamic variables
 452       ! divide by total DFI coefficient
 453       grid%mu_2(:,:)    = grid%dfi_mu(:,:)      / grid%hcoeff_tot
 454       grid%u_2(:,:,:)   = grid%dfi_u(:,:,:)     / grid%hcoeff_tot
 455       grid%v_2(:,:,:)   = grid%dfi_v(:,:,:)     / grid%hcoeff_tot
 456       grid%w_2(:,:,:)   = grid%dfi_w(:,:,:)     / grid%hcoeff_tot
 457       grid%ww(:,:,:)    = grid%dfi_ww(:,:,:)    / grid%hcoeff_tot
 458       grid%t_2(:,:,:)   = grid%dfi_t(:,:,:)     / grid%hcoeff_tot
 459       grid%phb(:,:,:)   = grid%dfi_phb(:,:,:)   / grid%hcoeff_tot
 460       grid%ph0(:,:,:)   = grid%dfi_ph0(:,:,:)   / grid%hcoeff_tot
 461       grid%php(:,:,:)   = grid%dfi_php(:,:,:)   / grid%hcoeff_tot
 462       grid%p(:,:,:)     = grid%dfi_p(:,:,:)     / grid%hcoeff_tot
 463       grid%ph_2(:,:,:)  = grid%dfi_ph(:,:,:)    / grid%hcoeff_tot
 464       grid%tke_2(:,:,:) = grid%dfi_tke(:,:,:)   / grid%hcoeff_tot
 465       grid%al(:,:,:)    = grid%dfi_al(:,:,:)    / grid%hcoeff_tot
 466       grid%alt(:,:,:)   = grid%dfi_alt(:,:,:)   / grid%hcoeff_tot
 467       grid%pb(:,:,:)    = grid%dfi_pb(:,:,:)    / grid%hcoeff_tot
 468       grid%moist(:,:,:,:)  = grid%dfi_moist(:,:,:,:)  / grid%hcoeff_tot
 469       grid%scalar(:,:,:,:) = grid%dfi_scalar(:,:,:,:) / grid%hcoeff_tot
 470
 471       ! restore initial fields
 472       grid%SNOW  (:,:) = grid%dfi_SNOW  (:,:)
 473       grid%SNOWH (:,:) = grid%dfi_SNOWH (:,:)
 474       grid%SNOWC (:,:) = grid%dfi_SNOWC (:,:)
 475       grid%CANWAT(:,:) = grid%dfi_CANWAT(:,:)
 476       grid%TSK   (:,:) = grid%dfi_TSK   (:,:)
 477       grid%QVG   (:,:) = grid%dfi_QVG   (:,:)
 478       grid%SOILT1(:,:) = grid%dfi_SOILT1(:,:)
 479       grid%TSNAV (:,:) = grid%dfi_TSNAV (:,:)
 480
 481       grid%TSLB  (:,:,:)       = grid%dfi_TSLB   (:,:,:)
 482       grid%SMOIS (:,:,:)       = grid%dfi_SMOIS  (:,:,:)
 483       grid%SMFR3D(:,:,:)       = grid%dfi_SMFR3D (:,:,:)
 484       grid%KEEPFR3DFLAG(:,:,:) = grid%dfi_KEEPFR3DFLAG(:,:,:)
 485
 486    END SUBROUTINE dfi_array_reset
 487
 488
 489    SUBROUTINE dfcoef (NSTEPS,DT,TAUC,NFILT,H)
 490 !
 491 !     calculate filter weights with selected window.
 492 !
 493 !       peter lynch and xiang-yu huang
 494 !
 495 !       ref: see hamming, r.w., 1989: digital filters,
 496 !                prentice-hall international. 3rd edition.
 497 !
 498 !       input:      nsteps  -  number of timesteps
 499 !                              forward or backward.
 500 !                   dt      -  time step in seconds.
 501 !                   tauc    -  cut-off period in seconds.
 502 !                   nfilt   -  indicator for selected filter.
 503 !
 504 !       output:     h       -  array(0:nsteps) with the
 505 !                              required filter weights
 506 !
 507 !------------------------------------------------------------
 508
 509       implicit none
 510
 511       integer, intent(in)   :: nsteps, nfilt
 512       real   , intent(in)   :: dt, tauc
 513       real, intent(out)  :: h(1:nsteps+1)
 514
 515       ! Local data
 516
 517       integer               :: n
 518       real                  :: pi, omegac, x, window, deltat
 519       real                  :: hh(0:nsteps)
 520
 521       pi=4*ATAN(1.)
 522       deltat=ABS(dt)
 523
 524       ! windows are defined by a call and held in hh.
 525
 526       if ( nfilt .eq. -1) call debug    (nsteps,h)
 527
 528       IF ( NFILT .EQ. 0 ) CALL UNIFORM  (NSTEPS,HH)
 529       IF ( NFILT .EQ. 1 ) CALL LANCZOS  (NSTEPS,HH)
 530       IF ( NFILT .EQ. 2 ) CALL HAMMING  (NSTEPS,HH)
 531       IF ( NFILT .EQ. 3 ) CALL BLACKMAN (NSTEPS,HH)
 532       IF ( NFILT .EQ. 4 ) CALL KAISER   (NSTEPS,HH)
 533       IF ( NFILT .EQ. 5 ) CALL POTTER2  (NSTEPS,HH)
 534       IF ( NFILT .EQ. 6 ) CALL DOLPHWIN (NSTEPS,HH)
 535
 536       IF ( NFILT .LE. 6 ) THEN     ! sinc-windowed filters
 537
 538          ! calculate the cutoff frequency
 539          OMEGAC = 2.*PI/TAUC
 540
 541          DO N=0,NSTEPS
 542             WINDOW = HH(N)
 543             IF ( N .EQ. 0 ) THEN
 544               X = (OMEGAC*DELTAT/PI)
 545             ELSE
 546               X = SIN(N*OMEGAC*DELTAT)/(N*PI)
 547             END IF
 548             HH(N) = X*WINDOW
 549          END DO
 550
 551          ! normalize the sums to be unity
 552          CALL NORMLZ(HH,NSTEPS)
 553
 554          DO N=0,NSTEPS
 555             H(N+1) = HH(NSTEPS-N)
 556          END DO
 557
 558       ELSE IF ( NFILT .EQ. 7 ) THEN     ! dolph filter
 559
 560          CALL DOLPH(DT,TAUC,NSTEPS,H)
 561
 562       ELSE IF ( NFILT .EQ. 8 ) THEN     ! 2nd order, 2nd type quick start filter (RHO)
 563
 564          CALL RHOFIL(DT,TAUC,2,NSTEPS*2,2,H,NSTEPS)
 565
 566       END IF
 567
 568       RETURN
 569
 570    END SUBROUTINE dfcoef
 571
 572
 573    SUBROUTINE NORMLZ(HH,NMAX)
 574
 575    ! normalize the sum of hh to be unity
 576
 577       implicit none
 578
 579       integer, intent(in)                       :: nmax
 580       real   , dimension(0:nmax), intent(out)   :: hh
 581
 582       ! local data
 583       real     :: sumhh
 584       integer  :: n
 585
 586       SUMHH = HH(0)
 587       DO N=1,NMAX
 588         SUMHH = SUMHH + 2*HH(N)
 589       ENDDO
 590       DO N=0,NMAX
 591         HH(N)  = HH(N)/SUMHH
 592       ENDDO
 593
 594       RETURN
 595
 596    END subroutine normlz
 597
 598
 599    subroutine debug(nsteps, ww)
 600
 601       implicit none
 602
 603       integer, intent(in)                        :: nsteps
 604       real   , dimension(0:nsteps), intent(out)  :: ww
 605       integer :: n
 606
 607       do n=0,nsteps
 608          ww(n)=0
 609       end do
 610
 611       ww(int(nsteps/2))=1
 612
 613       return
 614
 615    end subroutine debug
 616
 617
 618    SUBROUTINE UNIFORM(NSTEPS,WW)
 619
 620    !  define uniform or rectangular window function.
 621
 622       implicit none
 623
 624       integer, intent(in)                        :: nsteps
 625       real   , dimension(0:nsteps), intent(out)  :: ww
 626
 627       integer          :: n
 628
 629       DO N=0,NSTEPS
 630         WW(N) = 1.
 631       ENDDO
 632
 633       RETURN
 634
 635    END subroutine uniform
 636
 637
 638    SUBROUTINE LANCZOS(NSTEPS,WW)
 639
 640    ! define (genaralised) lanczos window function.
 641
 642       implicit none
 643
 644       integer,  parameter                      :: nmax = 1000
 645       integer,  intent(in)                     :: nsteps
 646       real   ,  dimension(0:nmax), intent(out) :: ww
 647       integer  ::  n
 648       real     :: power, pi, w
 649
 650       ! (for the usual lanczos window, power = 1 )
 651       POWER = 1
 652
 653       PI=4*ATAN(1.)
 654       DO N=0,NSTEPS
 655          IF ( N .EQ. 0 ) THEN
 656             W = 1.0
 657          ELSE
 658             W = SIN(N*PI/(NSTEPS+1)) / ( N*PI/(NSTEPS+1))
 659          ENDIF
 660          WW(N) = W**POWER
 661       ENDDO
 662
 663       RETURN
 664
 665    END SUBROUTINE lanczos
 666
 667
 668    SUBROUTINE HAMMING(NSTEPS,WW)
 669
 670    ! define (genaralised) hamming window function.
 671
 672       implicit none
 673
 674       integer, intent(in)           :: nsteps
 675       real, dimension(0:nsteps)    :: ww
 676       integer   ::   n
 677       real      :: alpha, pi, w
 678
 679       ! (for the usual hamming window, alpha=0.54,
 680       !      for the hann window, alpha=0.50).
 681       ALPHA=0.54
 682
 683       PI=4*ATAN(1.)
 684       DO N=0,NSTEPS
 685          IF ( N .EQ. 0 ) THEN
 686             W = 1.0
 687          ELSE
 688             W = ALPHA + (1-ALPHA)*COS(N*PI/(NSTEPS))
 689          ENDIF
 690          WW(N) = W
 691       ENDDO
 692
 693       RETURN
 694
 695    END SUBROUTINE hamming
 696
 697
 698    SUBROUTINE BLACKMAN(NSTEPS,WW)
 699
 700    ! define blackman window function.
 701
 702       implicit none
 703
 704       integer, intent(in)           :: nsteps
 705       real, dimension(0:nsteps)    :: ww
 706       integer  :: n
 707
 708       real :: pi, w
 709
 710       PI=4*ATAN(1.)
 711       DO N=0,NSTEPS
 712          IF ( N .EQ. 0 ) THEN
 713             W = 1.0
 714          ELSE
 715             W = 0.42 + 0.50*COS(  N*PI/(NSTEPS))   &
 716                      + 0.08*COS(2*N*PI/(NSTEPS))
 717          ENDIF
 718          WW(N) = W
 719       ENDDO
 720
 721       RETURN
 722
 723    END SUBROUTINE blackman
 724
 725
 726    SUBROUTINE KAISER(NSTEPS,WW)
 727
 728    ! define kaiser window function.
 729
 730       implicit none
 731
 732       real, external :: bessi0
 733
 734       integer, intent(in)           :: nsteps
 735       real, dimension(0:nsteps)    :: ww
 736       integer   :: n
 737       real      :: alpha, xi0a, xn, as
 738
 739       ALPHA = 1
 740
 741       XI0A =  BESSI0(ALPHA)
 742       DO N=0,NSTEPS
 743          XN = N
 744          AS = ALPHA*SQRT(1.-(XN/NSTEPS)**2)
 745          WW(N) = BESSI0(AS) / XI0A
 746       ENDDO
 747
 748       RETURN
 749
 750    END SUBROUTINE kaiser
 751
 752
 753    REAL FUNCTION BESSI0(X)
 754
 755    ! from numerical recipes (press, et al.)
 756
 757       implicit none
 758
 759       real(8)   :: Y
 760       real(8)   :: P1 = 1.0d0
 761       real(8)   :: P2 = 3.5156230D0
 762       real(8)   :: P3 = 3.0899424D0
 763       real(8)   :: P4 = 1.2067492D0
 764       real(8)   :: P5 = 0.2659732D0
 765       real(8)   :: P6 = 0.360768D-1
 766       real(8)   :: P7 = 0.45813D-2
 767
 768       real*8    :: Q1 = 0.39894228D0
 769       real*8    :: Q2 = 0.1328592D-1
 770       real*8    :: Q3 = 0.225319D-2
 771       real*8    :: Q4 = -0.157565D-2
 772       real*8    :: Q5 = 0.916281D-2
 773       real*8    :: Q6 = -0.2057706D-1
 774       real*8    :: Q7 = 0.2635537D-1
 775       real*8    :: Q8 = -0.1647633D-1
 776       real*8    :: Q9 = 0.392377D-2
 777
 778       real            :: x, ax
 779
 780
 781       IF (ABS(X).LT.3.75) THEN
 782         Y=(X/3.75)**2
 783         BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
 784       ELSE
 785         AX=ABS(X)
 786         Y=3.75/AX
 787         BESSI0=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4    &
 788            +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
 789       ENDIF
 790       RETURN
 791
 792    END FUNCTION bessi0
 793
 794
 795    SUBROUTINE POTTER2(NSTEPS,WW)
 796
 797       ! define potter window function.
 798       ! modified to fall off over twice the range.
 799
 800       implicit none
 801
 802       integer, intent(in)                       :: nsteps
 803       real, dimension(0:nsteps),intent(out)     ::  ww
 804       integer  :: n
 805       real     :: ck, sum, arg
 806
 807       ! local data
 808       real, dimension(0:3)   :: d
 809       real                   :: pi
 810       integer                :: ip
 811
 812       d(0) = 0.35577019
 813       d(1) = 0.2436983
 814       d(2) = 0.07211497
 815       d(3) = 0.00630165
 816
 817       PI=4*ATAN(1.)
 818
 819       CK = 1.0
 820       DO N=0,NSTEPS
 821          IF (N.EQ.NSTEPS) CK = 0.5
 822          ARG = PI*FLOAT(N)/FLOAT(NSTEPS)
 823 !min---        modification in next statement
 824          ARG = ARG/2.
 825 !min---        end of modification
 826          SUM = D(0)
 827          DO IP=1,3
 828             SUM = SUM + 2.*D(IP)*COS(ARG*FLOAT(IP))
 829          END DO
 830          WW(N) = CK*SUM
 831       END DO
 832
 833       RETURN
 834
 835    END SUBROUTINE potter2
 836
 837
 838    SUBROUTINE dolphwin(m, window)
 839
 840 !     calculation of dolph-chebyshev window or, for short,
 841 !     dolph window, using the expression in the reference:
 842 !
 843 !     antoniou, andreas, 1993: digital filters: analysis,
 844 !     design and applications. mcgraw-hill, inc., 689pp.
 845 !
 846 !     the dolph window is optimal in the following sense:
 847 !     for a given main-lobe width, the stop-band attenuation
 848 !     is minimal; for a given stop-band level, the main-lobe
 849 !     width is minimal.
 850 !
 851 !     it is possible to specify either the ripple-ratio r
 852 !     or the stop-band edge thetas.
 853
 854       IMPLICIT NONE
 855
 856       ! Arguments
 857       INTEGER, INTENT(IN)                  ::  m
 858       REAL, DIMENSION(0:M), INTENT(OUT)    ::  window
 859
 860       ! local data
 861       REAL, DIMENSION(0:2*M)               :: t
 862       REAL, DIMENSION(0:M)                 :: w, time
 863       REAL    :: pi, thetas, x0, term1, term2, rr, r, db, sum, arg
 864       INTEGER :: n, nm1, nt, i
 865
 866       PI = 4*ATAN(1.D0)
 867       THETAS = 2*PI/M
 868
 869       N = 2*M+1
 870       NM1 = N-1
 871       X0 = 1/COS(THETAS/2)
 872
 873       TERM1 = (X0 + SQRT(X0**2-1))**(FLOAT(N-1))
 874       TERM2 = (X0 - SQRT(X0**2-1))**(FLOAT(N-1))
 875       RR = 0.5*(TERM1+TERM2)
 876       R = 1/RR
 877       DB = 20*LOG10(R)
 878       WRITE(*,'(1X,''DOLPH: M,N='',2I8)')M,N
 879       WRITE(*,'(1X,''DOLPH: THETAS (STOP-BAND EDGE)='',F10.3)')THETAS
 880       WRITE(*,'(1X,''DOLPH: R,DB='',2F10.3)')R, DB
 881
 882       DO NT=0,M
 883         SUM = RR
 884         DO I=1,M
 885           ARG = X0*COS(I*PI/N)
 886           CALL CHEBY(T,NM1,ARG)
 887           TERM1 = T(NM1)
 888           TERM2 = COS(2*NT*PI*I/N)
 889           SUM = SUM + 2*TERM1*TERM2
 890         ENDDO
 891         W(NT) = SUM/N
 892         TIME(NT) = NT
 893       ENDDO
 894
 895 !     fill up the array for return
 896       DO NT=0,M
 897         WINDOW(NT) = W(NT)
 898       ENDDO
 899
 900       RETURN
 901
 902    END SUBROUTINE dolphwin
 903
 904
 905    SUBROUTINE dolph(deltat, taus, m, window)
 906
 907 !     calculation of dolph-chebyshev window or, for short,
 908 !     dolph window, using the expression in the reference:
 909 !
 910 !     antoniou, andreas, 1993: digital filters: analysis,
 911 !     design and applications. mcgraw-hill, inc., 689pp.
 912 !
 913 !     the dolph window is optimal in the following sense:
 914 !     for a given main-lobe width, the stop-band attenuation
 915 !     is minimal; for a given stop-band level, the main-lobe
 916 !     width is minimal.
 917
 918       IMPLICIT NONE
 919
 920       ! Arguments
 921       INTEGER, INTENT(IN)                  ::  m
 922       REAL, DIMENSION(0:M), INTENT(OUT)    ::  window
 923       REAL, INTENT(IN)                     :: deltat, taus
 924
 925       ! local data
 926       integer, PARAMETER        :: NMAX = 5000
 927       REAL, dimension(0:NMAX)   :: t, w, time
 928       real, dimension(0:2*nmax) :: w2
 929       INTEGER                   :: NPRPE=0        ! no of pe
 930       CHARACTER*80              :: MES
 931
 932       real    :: pi, thetas, x0, term1, term2, rr, r,db, sum, arg, sumw
 933       integer :: n, nm1, i, nt
 934
 935       PI = 4*ATAN(1.D0)
 936
 937       print *, 'in dfcoef, deltat = ', deltat, 'taus=',taus
 938
 939       N = 2*M+1
 940       NM1 = N-1
 941
 942       THETAS = 2*PI*ABS(DELTAT/TAUS)
 943       X0 = 1/COS(THETAS/2)
 944       TERM1 = (X0 + SQRT(X0**2-1))**(FLOAT(N-1))
 945       TERM2 = (X0 - SQRT(X0**2-1))**(FLOAT(N-1))
 946       RR = 0.5*(TERM1+TERM2)
 947       R = 1/RR
 948       DB = 20*LOG10(R)
 949
 950
 951       WRITE(*,'(1X,''DOLPH: M,N='',2I8)')M,N
 952       WRITE(*,'(1X,''DOLPH: THETAS (STOP-BAND EDGE)='',F10.3)')THETAS
 953       WRITE(*,'(1X,''DOLPH: R,DB='',2F10.3)')R, DB
 954
 955       DO NT=0,M
 956          SUM = 1
 957          DO I=1,M
 958             ARG = X0*COS(I*PI/N)
 959             CALL CHEBY(T,NM1,ARG)
 960             TERM1 = T(NM1)
 961             TERM2 = COS(2*NT*PI*I/N)
 962             SUM = SUM + R*2*TERM1*TERM2
 963          ENDDO
 964          W(NT) = SUM/N
 965          TIME(NT) = NT
 966          WRITE(*,'(1X,''DOLPH: TIME, W='',F10.6,2X,E17.7)')  &
 967            TIME(NT), W(NT)
 968       ENDDO
 969 !     fill in the negative-time values by symmetry.
 970       DO NT=0,M
 971          W2(M+NT) = W(NT)
 972          W2(M-NT) = W(NT)
 973       ENDDO
 974
 975 !     fill up the array for return
 976       SUMW = 0.
 977       DO NT=0,2*M
 978          SUMW = SUMW + W2(NT)
 979       ENDDO
 980       WRITE(*,'(1X,''DOLPH: SUM OF WEIGHTS W2='',F10.4)')SUMW
 981
 982       DO NT=0,2*M
 983          WINDOW(NT) = W2(NT)
 984       ENDDO
 985
 986       RETURN
 987
 988    END SUBROUTINE dolph
 989
 990
 991    SUBROUTINE cheby(t, n, x)
 992
 993 !     calculate all chebyshev polynomials up to order n
 994 !     for the argument value x.
 995
 996 !     reference: numerical recipes, page 184, recurrence
 997 !         t_n(x) = 2xt_{n-1}(x) - t_{n-2}(x) ,  n>=2.
 998
 999       IMPLICIT NONE
1000
1001       ! Arguments
1002       INTEGER, INTENT(IN)  :: n
1003       REAL, INTENT(IN)     :: x
1004       REAL, DIMENSION(0:N) :: t
1005
1006       integer  :: nn
1007
1008       T(0) = 1
1009       T(1) = X
1010       IF(N.LT.2) RETURN
1011       DO NN=2,N
1012          T(NN) = 2*X*T(NN-1) - T(NN-2)
1013       ENDDO
1014
1015       RETURN
1016
1017    END SUBROUTINE cheby
1018
1019
1020    SUBROUTINE rhofil(dt, tauc, norder, nstep, ictype, frow, nosize)
1021
1022 !       RHO = recurssive high order.
1023 !
1024 !       This routine calculates and returns the
1025 !       Last Row, FROW, of the FILTER matrix.
1026 !
1027 !       Input Parameters:
1028 !              DT  :  Time Step in seconds
1029 !            TAUC  :  Cut-off period (hours)
1030 !          NORDER  :  Order of QS Filter
1031 !           NSTEP  :  Number of step/Size of row.
1032 !          ICTYPE  :  Initial Conditions
1033 !          NOSIZE  :  Max. side of FROW.
1034 !
1035 !       Working Fields:
1036 !           ACOEF  :  X-coefficients of filter
1037 !           BCOEF  :  Y-coefficients of filter
1038 !          FILTER  :  Filter Matrix.
1039 !
1040 !       Output Parameters:
1041 !            FROW  : Last Row of Filter Matrix.
1042 !
1043 !       Note: Two types of initial conditions are permitted.
1044 !          ICTYPE = 1 : Order increasing each row to NORDER.
1045 !          ICTYPE = 2 : Order fixed at NORDER throughout.
1046 !
1047 !       DOUBLE PRECISION USED THROUGHOUT.
1048
1049       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1050
1051       DOUBLE PRECISION MUC
1052
1053 !       N.B. Single Precision for List Parameters.
1054       REAL, intent(in)    ::  DT,TAUC
1055
1056 !       Space for the last row of FILTER.
1057       integer, intent(in) ::  norder, nstep, ictype, nosize
1058       REAL   , dimension(0:nosize), intent(out)::  FROW
1059
1060 !       Arrays for rho filter.
1061       integer, PARAMETER  :: NOMAX=100
1062       real   , dimension(0:NOMAX) :: acoef, bcoef
1063       real   , dimension(0:NOMAX,0:NOMAX) :: filter
1064 !       Working space.
1065       real   , dimension(0:NOMAX) :: alpha, beta
1066
1067       real   :: DTT
1068
1069       DTT = ABS(DT)
1070       PI = 2*DASIN(1.D0)
1071       IOTA = CMPLX(0.,1.)
1072
1073 !       Filtering Parameters (derived).
1074       THETAC = 2*PI*DTT/(TAUC)
1075       MUC = tan(THETAC/2)
1076       FC = THETAC/(2*PI)
1077
1078 !     Clear the arrays.
1079     DO NC=0,NOMAX
1080        ACOEF(NC) = 0.
1081        BCOEF(NC) = 0.
1082        ALPHA(NC) = 0.
1083        BETA (NC) = 0.
1084        FROW (NC) = 0.
1085        DO NR=0,NOMAX
1086           FILTER(NR,NC) = 0.
1087        ENDDO
1088     ENDDO
1089
1090 !     Fill up the Filter Matrix.
1091     FILTER(0,0) = 1.
1092
1093 !     Get the coefficients of the Filter.
1094     IF ( ICTYPE.eq.2 ) THEN
1095        CALL RHOCOF(NORDER,NOMAX,MUC, ACOEF,BCOEF)
1096     ENDIF
1097
1098     DO 100 NROW=1,NSTEP
1099
1100        IF ( ICTYPE.eq.1 ) THEN
1101           NORD = MIN(NROW,NORDER)
1102           IF ( NORD.le.NORDER) THEN
1103              CALL RHOCOF(NORD,NOMAX,MUC, ACOEF,BCOEF)
1104           ENDIF
1105        ENDIF
1106
1107        DO K=0,NROW
1108           ALPHA(K) = ACOEF(NROW-K)
1109           IF(K.lt.NROW) BETA(K) = BCOEF(NROW-K)
1110        ENDDO
1111
1112 !        Correction for terms of negative index.
1113        IF ( ICTYPE.eq.2 ) THEN
1114           IF ( NROW.lt.NORDER ) THEN
1115              CN = 0.
1116              DO NN=NROW+1,NORDER
1117                 CN = CN + (ACOEF(NN)+BCOEF(NN))
1118              ENDDO
1119              ALPHA(0) = ALPHA(0) + CN
1120           ENDIF
1121        ENDIF
1122
1123 !       Check sum of ALPHAs and BETAs = 1
1124       SUMAB = 0.
1125       DO NN=0,NROW
1126         SUMAB = SUMAB + ALPHA(NN)
1127           IF(NN.lt.NROW) SUMAB = SUMAB + BETA(NN)
1128       ENDDO
1129
1130       DO KK=0,NROW-1
1131          SUMBF = 0.
1132          DO LL=0,NROW-1
1133             SUMBF = SUMBF + BETA(LL)*FILTER(LL,KK)
1134          ENDDO
1135          FILTER(NROW,KK) = ALPHA(KK)+SUMBF
1136       ENDDO
1137       FILTER(NROW,NROW) = ALPHA(NROW)
1138
1139 !       Check sum of row elements = 1
1140       SUMROW = 0.
1141       DO NN=0,NROW
1142         SUMROW = SUMROW + FILTER(NROW,NN)
1143       ENDDO
1144
1145 100 CONTINUE
1146
1147       DO NC=0,NSTEP
1148         FROW(NC) = FILTER(NSTEP,NC)
1149       ENDDO
1150
1151       RETURN
1152
1153    END SUBROUTINE rhofil
1154
1155
1156    SUBROUTINE rhocof(nord, nomax, muc, ca, cb)
1157
1158 !     Get the coefficients of the RHO Filter
1159
1160 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1161       IMPLICIT NONE
1162
1163       ! Arguments
1164       integer, intent(in)      :: nord, nomax
1165       real, dimension(0:nomax) :: ca, cb
1166
1167       ! Functions
1168       double precision, external :: cnr
1169
1170       ! Local variables
1171       INTEGER                  :: nn
1172       COMPLEX                  :: IOTA
1173       DOUBLE PRECISION         :: MUC, ZN
1174       DOUBLE PRECISION         :: pi, root2, rn, sigma, gain, sumcof
1175
1176       PI = 2*ASIN(1.)
1177       ROOT2 = SQRT(2.)
1178       IOTA = (0.,1.)
1179
1180       RN = 1./FLOAT(NORD)
1181       SIGMA = 1./( SQRT(2.**RN-1.) )
1182
1183       GAIN  = (MUC*SIGMA/(1+MUC*SIGMA))**NORD
1184       ZN = (1-MUC*SIGMA)/(1+MUC*SIGMA)
1185
1186       DO NN=0,NORD
1187         CA(NN) = CNR(NORD,NN)*GAIN
1188         IF(NN.gt.0) CB(NN) = -CNR(NORD,NN)*(-ZN)**NN
1189       ENDDO
1190
1191 !     Check sum of coefficients = 1
1192       SUMCOF = 0.
1193       DO NN=0,NORD
1194         SUMCOF = SUMCOF + CA(NN)
1195         IF(NN.gt.0) SUMCOF = SUMCOF + CB(NN)
1196       ENDDO
1197
1198       RETURN
1199
1200    END SUBROUTINE RHOCOF
1201
1202
1203    DOUBLE PRECISION FUNCTION cnr(n,r)
1204
1205    ! Binomial Coefficient (n,r).
1206
1207 !      IMPLICIT DOUBLE PRECISION(C,X)
1208       IMPLICIT NONE
1209
1210       ! Arguments
1211       INTEGER , intent(in)    :: n, R
1212
1213       ! Local variables
1214       INTEGER          :: k
1215       DOUBLE PRECISION :: coeff, xn, xr, xk
1216
1217       IF ( R.eq.0 ) THEN
1218          CNR = 1.0
1219          RETURN
1220       ENDIF
1221       Coeff = 1.0
1222       XN = DFLOAT(N)
1223       XR = DFLOAT(R)
1224       DO K=1,R
1225         XK = DFLOAT(K)
1226         COEFF = COEFF * ( (XN-XR+XK)/XK )
1227       ENDDO
1228       CNR = COEFF
1229
1230       RETURN
1231
1232    END FUNCTION cnr
1233
1234
1235     SUBROUTINE optfil (grid,NH,DELTAT,NHMAX)
1236 !----------------------------------------------------------------------
1237
1238 !    SUBROUTINE optfil (NH,DELTAT,TAUP,TAUS,LPRINT,         &
1239 !                         H,NHMAX)
1240 !
1241 ! - Huang and Lynch optimal filter
1242 !   Monthly Weather Review, Feb 1993
1243 !----------------------------------------------------------
1244 !       Input Parameters in List:
1245 !             NH     :  Half-length of the Filter
1246 !             DELTAT :  Time-step (in seconds).
1247 !             TAUP   :  Period of pass-band edge (hours).
1248 !             TAUS   :  Period of stop-band edge (hours).
1249 !             LPRINT :  Logical switch for messages.
1250 !             NHMAX  :  Maximum permitted Half-length.
1251 !
1252 !       Output Parameters in List:
1253 !             H      :  Impulse Response.
1254 !             DP     :  Deviation in pass-band (db)
1255 !             DS     :  Deviation in stop-band (db)
1256 !----------------------------------------------------------
1257 !
1258      USE module_domain
1259
1260      TYPE(domain) , POINTER :: grid
1261
1262      REAL,DIMENSION( 20) :: EDGE
1263      REAL,DIMENSION( 10) :: FX, WTX, DEVIAT
1264      REAL,DIMENSION(2*NHMAX+1) :: H
1265      logical LPRINT
1266      REAL, INTENT (IN) :: DELTAT
1267      INTEGER, INTENT (IN) :: NH, NHMAX
1268 !
1269          TAUP = 3.
1270          TAUS = 1.5
1271          LPRINT = .true.
1272 !initialize H array
1273
1274          NL=2*NHMAX+1
1275         do 101 n=1,NL
1276           H(n)=0.
1277  101    continue
1278
1279         NFILT = 2*NH+1
1280         print *,' start optfil, NFILT=', nfilt
1281
1282 !
1283 ! 930325 PL & XYH : the upper limit is changed from 64 to 128.
1284         IF(NFILT.LE.0 .OR. NFILT.GT.128 ) THEN
1285            WRITE(6,*) 'NH=',NH
1286        CALL wrf_error_fatal (' Sorry, error 1 in call to OPTFIL ')
1287         ENDIF
1288 !
1289 !       The following four should always be the same.
1290         JTYPE = 1
1291         NBANDS = 2
1292 !CC     JPRINT = 0
1293         LGRID = 16
1294 !
1295 !       calculate transition frequencies.
1296           DT = ABS(DELTAT)
1297           FS = DT/(TAUS*3600.)
1298           FP = DT/(TAUP*3600.)
1299           IF(FS.GT.0.5) then
1300 !            print *,' FS too large in OPTFIL '
1301        CALL wrf_error_fatal (' FS too large in OPTFIL ')
1302 !            return
1303           end if
1304           IF(FP.LT.0.0) then
1305 !            print *, ' FP too small in OPTFIL '
1306        CALL wrf_error_fatal (' FP too small in OPTFIL ')
1307 !            return
1308           end if
1309 !
1310 !       Relative Weights in pass- and stop-bands.
1311           WTP = 1.0
1312           WTS = 1.0
1313 !
1314 !CC     NOTE: (FP,FC,FS) is an arithmetic progression, so
1315 !CC           (1/FS,1/FC,1/FP) is a harmonic one.
1316 !CC           TAUP = 1/( (1/TAUC)-(1/DTAU) )
1317 !CC           TAUS = 1/( (1/TAUC)+(1/DTAU) )
1318 !CC           TAUC   :  Cut-off Period (hours).
1319 !CC           DTAU   :  Transition half-width (hours).
1320 !CC           FC = 1/TAUC  ;  DF = 1/DTAU
1321 !CC           FP = FC - DF :  FS = FC + DF
1322 !
1323           IF ( LPRINT ) THEN
1324                TAUC = 2./((1/TAUS)+(1/TAUP))
1325                DTAU = 2./((1/TAUS)-(1/TAUP))
1326                FC = DT/(TAUC*3600.)
1327                DF = DT/(DTAU*3600.)
1328                WRITE(6,*) ' DT ' , dt
1329                WRITE(6,*) ' TAUS, TAUP ' , TAUS,TAUP
1330                WRITE(6,*) ' TAUC, DTAU ' , TAUC,DTAU
1331                WRITE(6,*) '   FP,   FS ' ,   FP,  FS
1332                WRITE(6,*) '   FC,   DF ' ,   FC,  DF
1333                WRITE(6,*) '  WTS,  WTP ' ,  WTS, WTP
1334           ENDIF
1335 !
1336 !       Fill the control vectors for MCCPAR
1337         EDGE(1) = 0.0
1338         EDGE(2) = FP
1339         EDGE(3) = FS
1340         EDGE(4) = 0.5
1341         FX(1) = 1.0
1342         FX(2) = 0.0
1343         WTX(1) = WTP
1344         WTX(2) = WTS
1345
1346         CALL MCCPAR(NFILT,JTYPE,NBANDS,LPRINT,LGRID,       &
1347                    EDGE,FX,WTX,DEVIAT, h )
1348 !
1349 !       Save the deviations in the pass- and stop-bands.
1350         DP = DEVIAT(1)
1351         DS = DEVIAT(2)
1352 !
1353 !     Fill out the array H (only first half filled in MCCPAR).
1354       IF(MOD(NFILT,2).EQ.0) THEN
1355          NHALF = ( NFILT )/2
1356       ELSE
1357          NHALF = (NFILT+1)/2
1358       ENDIF
1359       DO 100 nn=1,NHALF
1360          H(NFILT+1-nn) = h(nn)
1361   100 CONTINUE
1362
1363 !       normalize the sums to be unity
1364         sumh = 0
1365         do 150 n=1,NFILT
1366           sumh = sumh + H(n)
1367   150   continue
1368   print *,'SUMH =', sumh
1369
1370         do 200 n=1,NFILT
1371           H(n)  = H(n)/sumh
1372   200   continue
1373         do 201 n=1,NFILT
1374         grid%hcoeff(n)=H(n)
1375   201   continue
1376 !       print *,'HCOEFF(n) ', grid%hcoeff
1377 !
1378    END SUBROUTINE optfil
1379
1380
1381       SUBROUTINE MCCPAR (NFILT,JTYPE,NBANDS,LPRINT,LGRID,    &
1382                    EDGE,FX,WTX,DEVIAT,h )
1383
1384 !      PROGRAM FOR THE DESIGN OF LINEAR PHASE FINITE IMPULSE
1385 !      REPONSE (FIR) FILTERS USING THE REMEZ EXCHANGE ALGORITHM
1386 !
1387 !************************************************************
1388 !*  Reference: McClellan, J.H., T.W. Parks and L.R.Rabiner, *
1389 !*             1973:  A computer program for designing      *
1390 !*             optimum FIR linear phase digital filters.    *
1391 !*             IEEE Trans. on Audio and Electroacoustics,   *
1392 !*             Vol AU-21, No. 6, 506-526.                   *
1393 !************************************************************
1394 !
1395 !      THREE TYPES OF FILTERS ARE INCLUDED -- BANDPASS FILTERS
1396 !      DIFFERENTIATORS, AND HILBERT TRANSFORM FILTERS
1397 !
1398 !---------------------------------------------------------------
1399 !
1400 !      COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
1401       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
1402       DIMENSION H(66)
1403       DIMENSION DES(1045),GRID(1045),WT(1045)
1404       DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10)
1405       DOUBLE PRECISION PI2,PI
1406       DOUBLE PRECISION AD,DEV,X,Y
1407       LOGICAL LPRINT
1408
1409       PI  = 3.141592653589793
1410       PI2 = 6.283185307179586
1411
1412 !      ......
1413
1414       NFMAX = 128
1415 100      CONTINUE
1416
1417 !      PROGRAM INPUT SECTION
1418
1419 !CC     READ(5,*) NFILT,JTYPE,NBANDS,JPRINT,LGRID
1420
1421       IF(NFILT.GT.NFMAX.OR.NFILT.LT.3) THEN
1422          CALL wrf_error_fatal (' **** ERROR IN INPUT DATA ****' )
1423       END IF
1424       IF(NBANDS.LE.0) NBANDS = 1
1425
1426 !      ....
1427
1428       IF(LGRID.LE.0) LGRID = 16
1429       JB = 2*NBANDS
1430 !cc       READ(5,*) (EDGE(J),J=1,JB)
1431 !cc       READ(5,*) (FX(J),J=1,NBANDS)
1432 !cc       READ(5,*) (WTX(J),J=1,NBANDS)
1433       IF(JTYPE.EQ.0) THEN
1434          CALL wrf_error_fatal (' **** ERROR IN INPUT DATA ****' )
1435       END IF
1436       NEG = 1
1437       IF(JTYPE.EQ.1) NEG = 0
1438       NODD = NFILT/2
1439       NODD = NFILT-2*NODD
1440       NFCNS = NFILT/2
1441       IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS = NFCNS+1
1442
1443 !      ...
1444
1445       GRID(1) = EDGE(1)
1446       DELF = LGRID*NFCNS
1447       DELF = 0.5/DELF
1448       IF(NEG.EQ.0) GOTO 135
1449       IF(EDGE(1).LT.DELF) GRID(1) = DELF
1450 135      CONTINUE
1451       J = 1
1452       L = 1
1453       LBAND = 1
1454 140      FUP = EDGE(L+1)
1455 145      TEMP = GRID(J)
1456
1457 !      ....
1458
1459       DES(J) = EFF(TEMP,FX,WTX,LBAND,JTYPE)
1460       WT(J) = WATE(TEMP,FX,WTX,LBAND,JTYPE)
1461       J = J+1
1462       GRID(J) = TEMP+DELF
1463       IF(GRID(J).GT.FUP) GOTO 150
1464       GOTO 145
1465 150      GRID(J-1) = FUP
1466       DES(J-1) = EFF(FUP,FX,WTX,LBAND,JTYPE)
1467       WT(J-1) = WATE(FUP,FX,WTX,LBAND,JTYPE)
1468       LBAND = LBAND+1
1469       L = L+2
1470       IF(LBAND.GT.NBANDS) GOTO 160
1471       GRID(J) = EDGE(L)
1472       GOTO 140
1473 160      NGRID = J-1
1474       IF(NEG.NE.NODD) GOTO 165
1475       IF(GRID(NGRID).GT.(0.5-DELF)) NGRID = NGRID-1
1476 165      CONTINUE
1477
1478 !      ......
1479
1480       IF(NEG) 170,170,180
1481 170      IF(NODD.EQ.1) GOTO 200
1482       DO 175 J=1,NGRID
1483             CHANGE = DCOS(PI*GRID(J))
1484             DES(J) = DES(J)/CHANGE
1485             WT(J) = WT(J)*CHANGE
1486 175      CONTINUE
1487       GOTO 200
1488 180      IF(NODD.EQ.1) GOTO 190
1489       DO 185 J = 1,NGRID
1490             CHANGE = DSIN(PI*GRID(J))
1491             DES(J) = DES(J)/CHANGE
1492             WT(J)  = WT(J)*CHANGE
1493 185      CONTINUE
1494       GOTO 200
1495 190      DO 195 J =1,NGRID
1496             CHANGE = DSIN(PI2*GRID(J))
1497             DES(J) = DES(J)/CHANGE
1498             WT(J)  = WT(J)*CHANGE
1499 195      CONTINUE
1500
1501 !      ......
1502
1503 200      TEMP = FLOAT(NGRID-1)/FLOAT(NFCNS)
1504       DO 210 J = 1,NFCNS
1505             IEXT(J) = (J-1)*TEMP+1
1506 210      CONTINUE
1507       IEXT(NFCNS+1) = NGRID
1508       NM1 = NFCNS-1
1509       NZ  = NFCNS+1
1510
1511 ! CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION PROBLEM
1512
1513       CALL REMEZ(EDGE,NBANDS,PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID)
1514
1515 ! CALCULATE THE IMPULSE RESPONSE
1516
1517       IF(NEG) 300,300,320
1518 300      IF(NODD.EQ.0) GOTO 310
1519       DO 305 J=1,NM1
1520             H(J) = 0.5*ALPHA(NZ-J)
1521 305      CONTINUE
1522       H(NFCNS)=ALPHA(1)
1523       GOTO 350
1524 310      H(1) = 0.25*ALPHA(NFCNS)
1525       DO 315 J = 2,NM1
1526             H(J) = 0.25*(ALPHA(NZ-J)+ALPHA(NFCNS+2-J))
1527 315      CONTINUE
1528       H(NFCNS) = 0.5*ALPHA(1)+0.25*ALPHA(2)
1529       GOTO 350
1530 320      IF(NODD.EQ.0) GOTO 330
1531       H(1) = 0.25*ALPHA(NFCNS)
1532       H(2) = 0.25*ALPHA(NM1)
1533       DO 325 J = 3,NM1
1534             H(J) = 0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+3-J))
1535 325      CONTINUE
1536       H(NFCNS) = 0.5*ALPHA(1)-0.25*ALPHA(3)
1537       H(NZ) = 0.0
1538       GOTO 350
1539 330      H(1) = 0.25*ALPHA(NFCNS)
1540       DO 335 J =2,NM1
1541             H(J) = 0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+2-J))
1542 335      CONTINUE
1543       H(NFCNS) = 0.5*ALPHA(1)-0.25*ALPHA(2)
1544
1545 !   PROGRAM OUTPUT SECTION
1546
1547 350   CONTINUE
1548 !
1549       IF(LPRINT) THEN
1550
1551          print *, '****************************************************'
1552          print *, 'FINITE IMPULSE RESPONSE (FIR)'
1553          print *, 'LINEAR PHASE DIGITAL FILTER DESIGN'
1554          print *, 'REMEZ EXCHANGE ALGORITHM'
1555
1556       IF(JTYPE.EQ.1) WRITE(6,365)
1557 365      FORMAT(25X,'BANDPASS FILTER'/)
1558
1559       IF(JTYPE.EQ.2) WRITE(6,370)
1560 370      FORMAT(25X,'DIFFERENTIATOR '/)
1561
1562       IF(JTYPE.EQ.3) WRITE(6,375)
1563 375      FORMAT(25X,'HILBERT TRANSFORMER '/)
1564
1565       WRITE(6,378) NFILT
1566 378      FORMAT(15X,'FILTER LENGTH =',I3/)
1567
1568       WRITE(6,380)
1569 380      FORMAT(15X,'***** IMPULSE RESPONSE *****')
1570
1571       DO 381 J = 1,NFCNS
1572             K = NFILT+1-J
1573             IF(NEG.EQ.0) WRITE(6,382) J,H(J),K
1574             IF(NEG.EQ.1) WRITE(6,383) J,H(J),K
1575 381      CONTINUE
1576 382      FORMAT(20X,'H(',I3,') = ',E15.8,' = H(',I4,')')
1577 383      FORMAT(20X,'H(',I3,') = ',E15.8,' = -H(',I4,')')
1578
1579       IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(6,384) NZ
1580 384      FORMAT(20X,'H(',I3,') = 0.0')
1581
1582       DO 450 K=1,NBANDS,4
1583             KUP = K+3
1584             IF(KUP.GT.NBANDS) KUP = NBANDS
1585             print *
1586             WRITE(6,385) (J,J=K,KUP)
1587 385            FORMAT(24X,4('BAND',I3,8X))
1588             WRITE(6,390) (EDGE(2*J-1),J=K,KUP)
1589 390            FORMAT(2X,'LOWER BAND EDGE',5F15.8)
1590             WRITE(6,395) (EDGE(2*J),J=K,KUP)
1591 395            FORMAT(2X,'UPPER BAND EDGE',5F15.8)
1592             IF(JTYPE.NE.2) WRITE(6,400) (FX(J),J=K,KUP)
1593 400            FORMAT(2X,'DESIRED VALUE',2X,5F15.8)
1594             IF(JTYPE.EQ.2) WRITE(6,405) (FX(J),J=K,KUP)
1595 405            FORMAT(2X,'DESIRED SLOPE',2X,5F15.8)
1596             WRITE(6,410) (WTX(J),J=K,KUP)
1597 410            FORMAT(2X,'WEIGHTING',6X,5F15.8)
1598             DO 420 J = K,KUP
1599                   DEVIAT(J) = DEV/WTX(J)
1600 420            CONTINUE
1601             WRITE(6,425) (DEVIAT(J),J=K,KUP)
1602 425            FORMAT(2X,'DEVIATION',6X,5F15.8)
1603             IF(JTYPE.NE.1) GOTO 450
1604             DO 430 J = K,KUP
1605                   DEVIAT(J) = 20.0*ALOG10(DEVIAT(J))
1606 430            CONTINUE
1607             WRITE(6,435) (DEVIAT(J),J=K,KUP)
1608 435            FORMAT(2X,'DEVIATION IN DB',5F15.8)
1609 450      CONTINUE
1610       print *, 'EXTREMAL FREQUENCIES'
1611       WRITE(6,455) (GRID(IEXT(J)),J=1,NZ)
1612 455      FORMAT((2X,5F15.7))
1613       WRITE(6,460)
1614 460      FORMAT(1X,70(1H*))
1615 !
1616       ENDIF
1617 !
1618 !CC     IF(NFILT.NE.0) GOTO 100  ! removal of re-run loop.
1619 !
1620    END SUBROUTINE mccpar
1621
1622
1623       FUNCTION EFF(TEMP,FX,WTX,LBAND,JTYPE)
1624          DIMENSION FX(5),WTX(5)
1625          IF(JTYPE.EQ.2) GOTO 1
1626          EFF = FX(LBAND)
1627          RETURN
1628 1        EFF = FX(LBAND)*TEMP
1629       END FUNCTION eff
1630
1631
1632       FUNCTION WATE(TEMP,FX,WTX,LBAND,JTYPE)
1633          DIMENSION FX(5),WTX(5)
1634          IF(JTYPE.EQ.2) GOTO 1
1635          WATE = WTX(LBAND)
1636          RETURN
1637 1        IF(FX(LBAND).LT.0.0001) GOTO 2
1638          WATE = WTX(LBAND)/TEMP
1639          RETURN
1640 2        WATE = WTX(LBAND)
1641       END FUNCTION wate
1642
1643
1644 !      SUBROUTINE ERROR
1645 !!         WRITE(6,*)' **** ERROR IN INPUT DATA ****'
1646 !          CALL wrf_error_fatal (' **** ERROR IN INPUT DATA ****' )
1647 !      END SUBROUTINE error
1648
1649
1650       SUBROUTINE REMEZ(EDGE,NBANDS,PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID)
1651 !  THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM
1652 !  FOR THE WEIGHTED CHEBCHEV APPROXIMATION OF A CONTINUOUS
1653 !  FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE
1654 !  ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE
1655 !  DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE
1656 !  GRID, THE NUMBER OF COSINES, AND THE INITIAL GUESS OF THE
1657 !  EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHEBYSHEV
1658 !  ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL
1659 !  FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES
1660 !  THE COEFFICIENTS OF THE BEST APPROXIMATION.
1661 !
1662 !      COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
1663       DIMENSION EDGE(20)
1664       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
1665       DIMENSION DES(1045),GRID(1045),WT(1045)
1666       DIMENSION A(66),P(65),Q(65)
1667       DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q
1668       DOUBLE PRECISION AD,DEV,X,Y
1669       DOUBLE PRECISION, EXTERNAL :: D, GEE
1670 !
1671 !  THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25
1672 !
1673       ITRMAX=25
1674       DEVL=-1.0
1675       NZ=NFCNS+1
1676       NZZ=NFCNS+2
1677       NITER=0
1678   100 CONTINUE
1679       IEXT(NZZ)=NGRID+1
1680       NITER=NITER+1
1681       IF(NITER.GT.ITRMAX) GO TO 400
1682       DO 110 J=1,NZ
1683       DTEMP=GRID(IEXT(J))
1684       DTEMP=DCOS(DTEMP*PI2)
1685   110 X(J)=DTEMP
1686       JET=(NFCNS-1)/15+1
1687       DO 120 J=1,NZ
1688   120 AD(J)=D(J,NZ,JET,X)
1689       DNUM=0.0
1690       DDEN=0.0
1691       K=1
1692       DO 130 J=1,NZ
1693       L=IEXT(J)
1694       DTEMP=AD(J)*DES(L)
1695       DNUM=DNUM+DTEMP
1696       DTEMP=K*AD(J)/WT(L)
1697       DDEN=DDEN+DTEMP
1698   130 K=-K
1699       DEV=DNUM/DDEN
1700       NU=1
1701       IF(DEV.GT.0.0) NU=-1
1702       DEV=-NU*DEV
1703       K=NU
1704       DO 140 J=1,NZ
1705       L=IEXT(J)
1706       DTEMP=K*DEV/WT(L)
1707       Y(J)=DES(L)+DTEMP
1708   140 K=-K
1709       IF(DEV.GE.DEVL) GO TO 150
1710       WRITE(6,*) ' ******** FAILURE TO CONVERGE *********** '
1711       WRITE(6,*) ' PROBABLE CAUSE IS MACHINE ROUNDING ERROR '
1712       WRITE(6,*) ' THE IMPULSE RESPONSE MAY BE CORRECT '
1713       WRITE(6,*) ' CHECK WITH A FREQUENCY RESPONSE '
1714       WRITE(6,*) ' **************************************** '
1715       GO TO 400
1716   150 DEVL=DEV
1717       JCHNGE=0
1718       K1=IEXT(1)
1719       KNZ=IEXT(NZ)
1720       KLOW=0
1721       NUT=-NU
1722       J=1
1723 !
1724 !  SEARCH FOR THE EXTERMAL FREQUENCIES OF THE BEST
1725 !  APPROXIMATION.
1726
1727   200 IF(J.EQ.NZZ) YNZ=COMP
1728       IF(J.GE.NZZ) GO TO 300
1729       KUP=IEXT(J+1)
1730       L=IEXT(J)+1
1731       NUT=-NUT
1732       IF(J.EQ.2) Y1=COMP
1733       COMP=DEV
1734       IF(L.GE.KUP) GO TO 220
1735       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
1736       ERR=(ERR-DES(L))*WT(L)
1737       DTEMP=NUT*ERR-COMP
1738       IF(DTEMP.LE.0.0) GO TO 220
1739       COMP=NUT*ERR
1740   210 L=L+1
1741       IF(L.GE.KUP) GO TO 215
1742       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
1743       ERR=(ERR-DES(L))*WT(L)
1744       DTEMP=NUT*ERR-COMP
1745       IF(DTEMP.LE.0.0) GO TO 215
1746       COMP=NUT*ERR
1747       GO TO 210
1748   215 IEXT(J)=L-1
1749       J=J+1
1750       KLOW=L-1
1751       JCHNGE=JCHNGE+1
1752       GO TO 200
1753   220 L=L-1
1754   225 L=L-1
1755       IF(L.LE.KLOW) GO TO 250
1756       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
1757       ERR=(ERR-DES(L))*WT(L)
1758       DTEMP=NUT*ERR-COMP
1759       IF(DTEMP.GT.0.0) GO TO 230
1760       IF(JCHNGE.LE.0) GO TO 225
1761       GO TO 260
1762   230 COMP=NUT*ERR
1763   235 L=L-1
1764       IF(L.LE.KLOW) GO TO 240
1765       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
1766       ERR=(ERR-DES(L))*WT(L)
1767       DTEMP=NUT*ERR-COMP
1768       IF(DTEMP.LE.0.0) GO TO 240
1769       COMP=NUT*ERR
1770       GO TO 235
1771   240 KLOW=IEXT(J)
1772       IEXT(J)=L+1
1773       J=J+1
1774       JCHNGE=JCHNGE+1
1775       GO TO 200
1776   250 L=IEXT(J)+1
1777       IF(JCHNGE.GT.0) GO TO 215
1778   255 L=L+1
1779       IF(L.GE.KUP) GO TO 260
1780       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
1781       ERR=(ERR-DES(L))*WT(L)
1782       DTEMP=NUT*ERR-COMP
1783       IF(DTEMP.LE.0.0) GO TO 255
1784       COMP=NUT*ERR
1785       GO TO 210
1786   260 KLOW=IEXT(J)
1787       J=J+1
1788       GO TO 200
1789   300 IF(J.GT.NZZ) GO TO 320
1790       IF(K1.GT.IEXT(1)) K1=IEXT(1)
1791       IF(KNZ.LT.IEXT(NZ)) KNZ=IEXT(NZ)
1792       NUT1=NUT
1793       NUT=-NU
1794       L=0
1795       KUP=K1
1796       COMP=YNZ*(1.00001)
1797       LUCK=1
1798   310 L=L+1
1799       IF(L.GE.KUP) GO TO 315
1800       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
1801       ERR=(ERR-DES(L))*WT(L)
1802       DTEMP=NUT*ERR-COMP
1803       IF(DTEMP.LE.0.0) GO TO 310
1804       COMP=NUT*ERR
1805       J=NZZ
1806       GO TO 210
1807   315 LUCK=6
1808       GO TO 325
1809   320 IF(LUCK.GT.9) GO TO 350
1810       IF(COMP.GT.Y1) Y1=COMP
1811       K1=IEXT(NZZ)
1812   325 L=NGRID+1
1813       KLOW=KNZ
1814       NUT=-NUT1
1815       COMP=Y1*(1.00001)
1816   330 L=L-1
1817       IF(L.LE.KLOW) GO TO 340
1818       ERR=GEE(L,NZ,GRID,PI2,X,Y,AD)
1819       ERR=(ERR-DES(L))*WT(L)
1820       DTEMP=NUT*ERR-COMP
1821       IF(DTEMP.LE.0.0) GO TO 330
1822       J=NZZ
1823       COMP=NUT*ERR
1824       LUCK=LUCK+10
1825       GO TO 235
1826   340 IF(LUCK.EQ.6) GO TO 370
1827       DO 345 J=1,NFCNS
1828   345 IEXT(NZZ-J)=IEXT(NZ-J)
1829       IEXT(1)=K1
1830       GO TO 100
1831   350 KN=IEXT(NZZ)
1832       DO 360 J=1,NFCNS
1833   360 IEXT(J)=IEXT(J+1)
1834       IEXT(NZ)=KN
1835       GO TO 100
1836   370 IF(JCHNGE.GT.0) GO TO 100
1837 !
1838 !  CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
1839 !  USING THE INVERSE DISCRETE FOURIER TRANSFORM.
1840 !
1841   400 CONTINUE
1842       NM1=NFCNS-1
1843       FSH=1.0E-06
1844       GTEMP=GRID(1)
1845       X(NZZ)=-2.0
1846       CN=2*NFCNS-1
1847       DELF=1.0/CN
1848       L=1
1849       KKK=0
1850       IF(EDGE(1).EQ.0.0.AND.EDGE(2*NBANDS).EQ.0.5) KKK=1
1851       IF(NFCNS.LE.3) KKK=1
1852       IF(KKK.EQ.1) GO TO 405
1853       DTEMP=DCOS(PI2*GRID(1))
1854       DNUM=DCOS(PI2*GRID(NGRID))
1855       AA=2.0/(DTEMP-DNUM)
1856       BB=-(DTEMP+DNUM)/(DTEMP-DNUM)
1857   405 CONTINUE
1858       DO 430 J=1,NFCNS
1859       FT=(J-1)*DELF
1860       XT=DCOS(PI2*FT)
1861       IF(KKK.EQ.1) GO TO 410
1862       XT=(XT-BB)/AA
1863 ! original :      FT=ARCOS(XT)/PI2
1864       FT=ACOS(XT)/PI2
1865   410 XE=X(L)
1866       IF(XT.GT.XE) GO TO 420
1867       IF((XE-XT).LT.FSH) GO TO 415
1868       L=L+1
1869       GO TO 410
1870   415 A(J)=Y(L)
1871       GO TO 425
1872   420 IF((XT-XE).LT.FSH) GO TO 415
1873       GRID(1)=FT
1874       A(J)=GEE(1,NZ,GRID,PI2,X,Y,AD)
1875   425 CONTINUE
1876       IF(L.GT.1) L=L-1
1877   430 CONTINUE
1878       GRID(1)=GTEMP
1879       DDEN=PI2/CN
1880       DO 510 J=1,NFCNS
1881       DTEMP=0.0
1882       DNUM=(J-1)*DDEN
1883       IF(NM1.LT.1) GO TO 505
1884       DO 500 K=1,NM1
1885   500 DTEMP=DTEMP+A(K+1)*DCOS(DNUM*K)
1886   505 DTEMP=2.0*DTEMP+A(1)
1887   510 ALPHA(J)=DTEMP
1888       DO 550 J=2,NFCNS
1889   550 ALPHA(J)=2*ALPHA(J)/CN
1890       ALPHA(1)=ALPHA(1)/CN
1891       IF(KKK.EQ.1) GO TO 545
1892       P(1)=2.0*ALPHA(NFCNS)*BB+ALPHA(NM1)
1893       P(2)=2.0*AA*ALPHA(NFCNS)
1894       Q(1)=ALPHA(NFCNS-2)-ALPHA(NFCNS)
1895       DO 540 J=2,NM1
1896       IF(J.LT.NM1) GO TO 515
1897       AA=0.5*AA
1898       BB=0.5*BB
1899   515 CONTINUE
1900       P(J+1)=0.0
1901       DO 520 K=1,J
1902       A(K)=P(K)
1903   520 P(K)=2.0*BB*A(K)
1904       P(2)=P(2)+A(1)*2.0*AA
1905       JM1=J-1
1906       DO 525 K=1,JM1
1907   525 P(K)=P(K)+Q(K)+AA*A(K+1)
1908       JP1=J+1
1909       DO 530 K=3,JP1
1910   530 P(K)=P(K)+AA*A(K-1)
1911       IF(J.EQ.NM1) GO TO 540
1912       DO 535 K=1,J
1913   535 Q(K)=-A(K)
1914       Q(1)=Q(1)+ALPHA(NFCNS-1-J)
1915   540 CONTINUE
1916       DO 543 J=1,NFCNS
1917   543 ALPHA(J)=P(J)
1918   545 CONTINUE
1919       IF(NFCNS.GT.3) RETURN
1920       ALPHA(NFCNS+1)=0.0
1921       ALPHA(NFCNS+2)=0.0
1922    END SUBROUTINE remez
1923
1924    DOUBLE PRECISION FUNCTION D(K,N,M,X)
1925 !      COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
1926       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
1927       DIMENSION DES(1045),GRID(1045),WT(1045)
1928       DOUBLE PRECISION AD,DEV,X,Y
1929       DOUBLE PRECISION Q
1930       DOUBLE PRECISION PI2
1931       D = 1.0
1932       Q = X(K)
1933       DO 3 L = 1,M
1934       DO 2 J = L,N,M
1935             IF(J-K) 1,2,1
1936 1                  D = 2.0*D*(Q-X(J))
1937 2      CONTINUE
1938 3      CONTINUE
1939       D = 1.0/D
1940    END FUNCTION D
1941
1942
1943    DOUBLE PRECISION FUNCTION GEE(K,N,GRID,PI2,X,Y,AD)
1944 !      COMMON /x3x/ PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
1945       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
1946       DIMENSION DES(1045),GRID(1045),WT(1045)
1947       DOUBLE PRECISION AD,DEV,X,Y
1948       DOUBLE PRECISION P,C,D,XF
1949       DOUBLE PRECISION PI2
1950       P = 0.0
1951       XF = GRID(K)
1952       XF = DCOS(PI2*XF)
1953       D = 0.0
1954       DO 1 J =1,N
1955             C = XF-X(J)
1956             C = AD(J)/C
1957             D = D+C
1958             P = P+C*Y(J)
1959 1      CONTINUE
1960       GEE = P/D
1961    END FUNCTION GEE
1962
1963 #endif