dyn_em/module_polarfft.F

   1
   2 MODULE module_polarfft
   3
   4   USE module_model_constants
   5   USE module_wrf_error
   6   CHARACTER (LEN=256) , PRIVATE :: a_message
   7
   8 CONTAINS
   9
  10 SUBROUTINE couple_scalars_for_filter ( field    &
  11                  ,mu,mub,c1,c2                  &
  12                  ,ids,ide,jds,jde,kds,kde       &
  13                  ,ims,ime,jms,jme,kms,kme       &
  14                  ,ips,ipe,jps,jpe,kps,kpe       )
  15    IMPLICIT NONE
  16    INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde       &
  17                          ,ims,ime,jms,jme,kms,kme       &
  18                          ,ips,ipe,jps,jpe,kps,kpe
  19    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: field
  20    REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub
  21    REAL , DIMENSION(kms:kme) , INTENT(IN) :: c1,c2
  22
  23    INTEGER :: i , j , k
  24
  25    DO j = jps, MIN(jpe,jde-1)
  26    DO k = kps, kpe-1
  27    DO i = ips, MIN(ipe,ide-1)
  28       field(i,k,j)=field(i,k,j)*(  (c1(k)*mu(i,j))+(c1(k)*mub(i,j)+c2(k)))
  29    END DO
  30    END DO
  31    END DO
  32
  33 END SUBROUTINE couple_scalars_for_filter
  34
  35 SUBROUTINE uncouple_scalars_for_filter ( field    &
  36                  ,mu,mub,c1,c2                  &
  37                  ,ids,ide,jds,jde,kds,kde       &
  38                  ,ims,ime,jms,jme,kms,kme       &
  39                  ,ips,ipe,jps,jpe,kps,kpe       )
  40    IMPLICIT NONE
  41    INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde       &
  42                          ,ims,ime,jms,jme,kms,kme       &
  43                          ,ips,ipe,jps,jpe,kps,kpe
  44    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: field
  45    REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: mu,mub
  46    REAL , DIMENSION(kms:kme) , INTENT(IN) :: c1,c2
  47
  48    INTEGER :: i , j , k
  49
  50    DO j = jps, MIN(jpe,jde-1)
  51    DO k = kps, kpe-1
  52    DO i = ips, MIN(ipe,ide-1)
  53       field(i,k,j)=field(i,k,j)/(  (c1(k)*mu(i,j))+(c1(k)*mub(i,j)+c2(k)))
  54    END DO
  55    END DO
  56    END DO
  57
  58 END SUBROUTINE uncouple_scalars_for_filter
  59
  60 SUBROUTINE pxft ( grid                          &
  61                  ,lineno       &
  62                  ,flag_uv,flag_rurv             &
  63                  ,flag_wph,flag_ww              &
  64                  ,flag_t                        &
  65                  ,flag_mu,flag_mut              &
  66                  ,flag_moist                    &
  67                  ,flag_chem                     &
  68                  ,flag_tracer                   &
  69                  ,flag_scalar                   &
  70                  ,fft_filter_lat, dclat         &
  71                  ,actual_distance_average       &
  72                  ,pos_def                       &
  73                  ,swap_pole_with_next_j         &
  74                  ,moist,chem,tracer,scalar      &
  75                  ,ids,ide,jds,jde,kds,kde       &
  76                  ,ims,ime,jms,jme,kms,kme       &
  77                  ,ips,ipe,jps,jpe,kps,kpe       &
  78                  ,imsx,imex,jmsx,jmex,kmsx,kmex &
  79                  ,ipsx,ipex,jpsx,jpex,kpsx,kpex )
  80    USE module_state_description
  81    USE module_domain, ONLY : domain
  82 #ifdef DM_PARALLEL
  83    USE module_dm, ONLY : local_communicator, mytask, ntasks, ntasks_x, ntasks_y &
  84                        , local_communicator_periodic, itrace                    &
  85                        , local_communicator_x
  86    USE module_driver_constants
  87 #if 0
  88    USE module_comm_dm, ONLY : &
  89                  XPOSE_POLAR_FILTER_V_z2x_sub           &
  90                 ,XPOSE_POLAR_FILTER_V_x2z_sub           &
  91                 ,XPOSE_POLAR_FILTER_U_z2x_sub           &
  92                 ,XPOSE_POLAR_FILTER_U_x2z_sub           &
  93                 ,XPOSE_POLAR_FILTER_T_z2x_sub           &
  94                 ,XPOSE_POLAR_FILTER_T_x2z_sub           &
  95                 ,XPOSE_POLAR_FILTER_W_z2x_sub           &
  96                 ,XPOSE_POLAR_FILTER_W_x2z_sub           &
  97                 ,XPOSE_POLAR_FILTER_PH_z2x_sub          &
  98                 ,XPOSE_POLAR_FILTER_PH_x2z_sub          &
  99                 ,XPOSE_POLAR_FILTER_WW_z2x_sub          &
 100                 ,XPOSE_POLAR_FILTER_WW_x2z_sub          &
 101                 ,XPOSE_POLAR_FILTER_RV_z2x_sub          &
 102                 ,XPOSE_POLAR_FILTER_RV_x2z_sub          &
 103                 ,XPOSE_POLAR_FILTER_RU_z2x_sub          &
 104                 ,XPOSE_POLAR_FILTER_RU_x2z_sub          &
 105                 ,XPOSE_POLAR_FILTER_MOIST_z2x_sub       &
 106                 ,XPOSE_POLAR_FILTER_MOIST_x2z_sub       &
 107                 ,XPOSE_POLAR_FILTER_CHEM_z2x_sub        &
 108                 ,XPOSE_POLAR_FILTER_MOIST_x2z_sub       &
 109                 ,XPOSE_POLAR_FILTER_TRACER_z2x_sub      &
 110                 ,XPOSE_POLAR_FILTER_TRACER_x2z_sub      &
 111                 ,XPOSE_POLAR_FILTER_SCALAR_z2x_sub      &
 112                 ,XPOSE_POLAR_FILTER_SCALAR_x2z_sub
 113 #endif
 114 #else
 115    USE module_dm
 116 #endif
 117    IMPLICIT NONE
 118    !  Input data.
 119    TYPE(domain) , TARGET          :: grid
 120 integer, intent(in) :: lineno
 121 integer myproc, i, j, k
 122    LOGICAL, INTENT(IN) :: actual_distance_average
 123    LOGICAL, INTENT(IN) :: pos_def
 124    LOGICAL, INTENT(IN) :: swap_pole_with_next_j
 125    INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde       &
 126                          ,ims,ime,jms,jme,kms,kme       &
 127                          ,ips,ipe,jps,jpe,kps,kpe       &
 128                          ,imsx,imex,jmsx,jmex,kmsx,kmex &
 129                          ,ipsx,ipex,jpsx,jpex,kpsx,kpex
 130    REAL  , INTENT(IN) :: fft_filter_lat
 131    REAL,    INTENT(IN) :: dclat
 132    INTEGER, INTENT(IN) :: flag_uv                       &
 133                          ,flag_rurv                     &
 134                          ,flag_ww                       &
 135                          ,flag_t,flag_wph               &
 136                          ,flag_mu,flag_mut              &
 137                          ,flag_moist                    &
 138                          ,flag_chem                     &
 139                          ,flag_tracer                   &
 140                          ,flag_scalar
 141     REAL, DIMENSION(ims:ime,kms:kme,jms:jme,*) , INTENT(INOUT) :: moist, chem, scalar,tracer
 142
 143    ! Local
 144    LOGICAL piggyback_mu, piggyback_mut
 145    INTEGER ij, k_end
 146
 147 #ifdef DM_PARALLEL
 148 #else
 149    INTEGER itrace
 150 #endif
 151
 152
 153    piggyback_mu  = flag_mu .EQ. 1
 154    piggyback_mut = flag_mut .EQ. 1
 155
 156 !<DESCRIPTION>
 157 !
 158 ! The idea is that this parallel transpose fft routine can be
 159 ! called at various points in the solver (solve_em) and it will filter
 160 ! the correct fields based on the flag arguments.  There are two 2d
 161 ! fields mu_2 and mut  that may need to be filtered too. Since a two-d
 162 ! parallel transpose would be inefficient and because the fields that are
 163 ! not staggered in z have an extra layer anyway, these can be
 164 ! piggybacked.  This is handled using latches to makes sure that *somebody*
 165 ! carries one or both of these on their back through the filtering and then
 166 ! copies them back afterwards. IMPORTANT NOTE: for simplicity, this routine
 167 ! is not completely general.  It makes the following assumptions:
 168 !
 169 ! 1) If both flag_mu and flag_mut are specified then flag_uv is also specified
 170 !
 171 ! 2) If flag_uv is not specified, then only flag_mu and not flag_mut can be
 172 !
 173 ! 3) If flag_mu is specified than either flag_uv or flag_t must be
 174 !
 175 ! This is designed to keep the clutter to a minimum in solve_em.
 176 ! This is not intended to be a general abstraction of the polar filtering
 177 ! calls in in WRF solvers or if the solve_em algorithms change.
 178 ! If the needs of the calling solver change, this logic may have to be
 179 ! rewritten.
 180 !
 181 !</DESCRIPTION>
 182 !write(0,*)"<stdin>",__LINE__,' short circuit '
 183 !return
 184
 185 !write(0,*)'pxft called from ',lineno
 186 call wrf_get_myproc(myproc)
 187 !write(20+myproc,*)ipex-ipsx+1,jpex-jpsx+1,' clat_xxx '
 188 !do j = jpsx, jpex
 189 !do i = ipsx, ipex
 190 !write(20+myproc,*)grid%clat_xxx(i,j)
 191 !enddo
 192 !enddo
 193
 194 !!!!!!!!!!!!!!!!!!!!!!!
 195 ! U & V
 196    IF ( flag_uv .EQ. 1 ) THEN
 197      IF ( piggyback_mu ) THEN
 198        grid%u_2(ips:ipe,kde,jps:jpe) = grid%mu_2(ips:ipe,jps:jpe)
 199      ENDIF
 200 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 201 # include "XPOSE_POLAR_FILTER_V_z2x.inc"
 202
 203      CALL polar_filter_3d( grid%v_xxx, grid%clat_xxx, .false.,     &
 204                                 fft_filter_lat, dclat,                 &
 205                                 ids, ide, jds, jde, kds, kde-1,         &
 206                                 imsx, imex, jmsx, jmex, kmsx, kmex,     &
 207                                 ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex ) )
 208
 209 # include "XPOSE_POLAR_FILTER_V_x2z.inc"
 210 # include "XPOSE_POLAR_FILTER_U_z2x.inc"
 211      k_end = MIN(kde-1,kpex)
 212      IF ( piggyback_mu ) k_end = MIN(kde,kpex)
 213
 214      CALL polar_filter_3d( grid%u_xxx, grid%clat_xxx, piggyback_mu,     &
 215                                 fft_filter_lat, 0.,                &
 216                                 ids, ide, jds, jde, kds, kde,       &
 217                                 imsx, imex, jmsx, jmex, kmsx, kmex, &
 218                                 ipsx, ipex, jpsx, jpex, kpsx, k_end )
 219
 220 # include "XPOSE_POLAR_FILTER_U_x2z.inc"
 221 #else
 222
 223      CALL polar_filter_3d( grid%v_2, grid%clat, .false.,     &
 224                                 fft_filter_lat, dclat,             &
 225                                 ids, ide, jds, jde, kds, kde,       &
 226                                 ims, ime, jms, jme, kms, kme,       &
 227                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
 228
 229      k_end = MIN(kde-1,kpe)
 230      IF ( piggyback_mu ) k_end = MIN(kde,kpe)
 231
 232      CALL polar_filter_3d( grid%u_2, grid%clat, piggyback_mu,     &
 233                                 fft_filter_lat, 0.,                &
 234                                 ids, ide, jds, jde, kds, kde-1,     &
 235                                 ims, ime, jms, jme, kms, kme,       &
 236                                 ips, ipe, jps, jpe, kps, k_end )
 237
 238 #endif
 239
 240      IF ( piggyback_mu ) THEN
 241        grid%mu_2(ips:ipe,jps:jpe) = grid%u_2(ips:ipe,kde,jps:jpe)
 242        piggyback_mu = .FALSE.
 243      ENDIF
 244    ENDIF
 245
 246 !!!!!!!!!!!!!!!!!!!!!!!
 247 ! T
 248    IF ( flag_t .EQ. 1 ) THEN
 249      IF ( piggyback_mu ) THEN
 250        grid%t_2(ips:ipe,kde,jps:jpe) = grid%mu_2(ips:ipe,jps:jpe)
 251      ENDIF
 252 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 253 # include "XPOSE_POLAR_FILTER_T_z2x.inc"
 254      k_end = MIN(kde-1,kpex)
 255      IF ( piggyback_mu ) k_end = MIN(kde,kpex)
 256
 257      CALL polar_filter_3d( grid%t_xxx, grid%clat_xxx,piggyback_mu,     &
 258                                 fft_filter_lat, 0.,                &
 259                                 ids, ide, jds, jde, kds, kde-1,     &
 260                                 imsx, imex, jmsx, jmex, kmsx, kmex, &
 261                                 ipsx, ipex, jpsx, jpex, kpsx, k_end )
 262
 263      IF ( actual_distance_average ) THEN
 264         CALL filter_tracer ( grid%t_xxx , grid%clat_xxx , grid%mf_xxx , &
 265                              grid%fft_filter_lat , grid%mf_fft , &
 266                              pos_def, swap_pole_with_next_j , &
 267                              ids, ide, jds, jde, kds, kde , &
 268                              imsx, imex, jmsx, jmex, kmsx, kmex, &
 269                              ipsx, ipex, jpsx, jpex, kpsx, kpex )
 270      END IF
 271
 272 # include "XPOSE_POLAR_FILTER_T_x2z.inc"
 273 #else
 274      k_end = MIN(kde-1,kpe)
 275      IF ( piggyback_mu ) k_end = MIN(kde,kpe)
 276
 277      CALL polar_filter_3d( grid%t_2, grid%clat, piggyback_mu,     &
 278                                 fft_filter_lat, 0.,                &
 279                                 ids, ide, jds, jde, kds, kde-1,     &
 280                                 ims, ime, jms, jme, kms, kme,       &
 281                                 ips, ipe, jps, jpe, kps, k_end )
 282
 283      IF ( actual_distance_average ) THEN
 284         CALL filter_tracer ( grid%t_2 , grid%clat , grid%msft , &
 285                              grid%fft_filter_lat , grid%mf_fft , &
 286                              pos_def, swap_pole_with_next_j , &
 287                              ids, ide, jds, jde, kds, kde , &
 288                              ims, ime, jms, jme, kms, kme, &
 289                              ips, ipe, jps, jpe, kps, k_end )
 290      END IF
 291
 292 #endif
 293      IF ( piggyback_mu ) THEN
 294        grid%mu_2(ips:ipe,jps:jpe) = grid%t_2(ips:ipe,kde,jps:jpe)
 295        piggyback_mu = .FALSE.
 296      ENDIF
 297    ENDIF
 298
 299 !!!!!!!!!!!!!!!!!!!!!!!
 300 ! W and PH
 301    IF ( flag_wph .EQ. 1 ) THEN
 302       ! W AND PH USE ALL LEVELS SO NEVER PIGGYBACK, MU IS OUT OF LUCK HERE
 303 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 304 # include "XPOSE_POLAR_FILTER_W_z2x.inc"
 305       CALL polar_filter_3d( grid%w_xxx, grid%clat_xxx, .false.,     &
 306                                  fft_filter_lat, 0.,                &
 307                                  ids, ide, jds, jde, kds, kde,       &
 308                                  imsx, imex, jmsx, jmex, kmsx, kmex, &
 309                                  ipsx, ipex, jpsx, jpex, kpsx, kpex )
 310
 311 # include "XPOSE_POLAR_FILTER_W_x2z.inc"
 312 # include "XPOSE_POLAR_FILTER_PH_z2x.inc"
 313
 314       CALL polar_filter_3d( grid%ph_xxx, grid%clat_xxx, .false.,     &
 315                                  fft_filter_lat, 0.,                &
 316                                  ids, ide, jds, jde, kds, kde,       &
 317                                  imsx, imex, jmsx, jmex, kmsx, kmex, &
 318                                  ipsx, ipex, jpsx, jpex, kpsx, kpex )
 319
 320 # include "XPOSE_POLAR_FILTER_PH_x2z.inc"
 321 #else
 322
 323       CALL polar_filter_3d( grid%w_2, grid%clat,  .false.,     &
 324                                  fft_filter_lat, 0.,                &
 325                                  ids, ide, jds, jde, kds, kde,       &
 326                                  ims, ime, jms, jme, kms, kme, &
 327                                  ips, ipe, jps, jpe, kps, kpe )
 328
 329       CALL polar_filter_3d( grid%ph_2, grid%clat, .false.,     &
 330                                  fft_filter_lat, 0.,                &
 331                                  ids, ide, jds, jde, kds, kde,       &
 332                                  ims, ime, jms, jme, kms, kme, &
 333                                  ips, ipe, jps, jpe, kps, kpe )
 334 #endif
 335    ENDIF
 336
 337 !!!!!!!!!!!!!!!!!!!!!!!
 338 ! WW
 339    IF ( flag_ww .EQ. 1 ) THEN
 340       ! WW USES ALL LEVELS SO NEVER PIGGYBACK, MU IS OUT OF LUCK HERE
 341 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 342 # include "XPOSE_POLAR_FILTER_WW_z2x.inc"
 343       CALL polar_filter_3d( grid%ww_xxx, grid%clat_xxx, .false.,     &
 344                                  fft_filter_lat, 0.,                &
 345                                  ids, ide, jds, jde, kds, kde,       &
 346                                  imsx, imex, jmsx, jmex, kmsx, kmex, &
 347                                  ipsx, ipex, jpsx, jpex, kpsx, kpex )
 348
 349 # include "XPOSE_POLAR_FILTER_WW_x2z.inc"
 350 #else
 351       CALL polar_filter_3d( grid%ww_m, grid%clat, .false.,     &
 352                                  fft_filter_lat, 0.,                &
 353                                  ids, ide, jds, jde, kds, kde,       &
 354                                  ims, ime, jms, jme, kms, kme, &
 355                                  ips, ipe, jps, jpe, kps, kpe )
 356 #endif
 357    ENDIF
 358
 359 !!!!!!!!!!!!!!!!!!!!!!!
 360 ! RU AND RV
 361    IF ( flag_rurv .EQ. 1 ) THEN
 362      IF ( piggyback_mut ) THEN
 363        grid%ru_m(ips:ipe,kde,jps:jpe) = grid%mut(ips:ipe,jps:jpe)
 364      ENDIF
 365 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 366 # include "XPOSE_POLAR_FILTER_RV_z2x.inc"
 367      CALL polar_filter_3d( grid%rv_xxx, grid%clat_xxx, .false.,     &
 368                                 fft_filter_lat, dclat,             &
 369                                 ids, ide, jds, jde, kds, kde,       &
 370                                 imsx, imex, jmsx, jmex, kmsx, kmex, &
 371                                 ipsx, ipex, jpsx, jpex, kpsx, MIN(kpex,kde-1) )
 372
 373 # include "XPOSE_POLAR_FILTER_RV_x2z.inc"
 374 # include "XPOSE_POLAR_FILTER_RU_z2x.inc"
 375      k_end = MIN(kde-1,kpex)
 376      IF ( piggyback_mut ) k_end = MIN(kde,kpex)
 377
 378      CALL polar_filter_3d( grid%ru_xxx, grid%clat_xxx, piggyback_mut,     &
 379                                 fft_filter_lat, 0.,                &
 380                                 ids, ide, jds, jde, kds, kde,       &
 381                                 imsx, imex, jmsx, jmex, kmsx, kmex, &
 382                                 ipsx, ipex, jpsx, jpex, kpsx, k_end )
 383 #include "XPOSE_POLAR_FILTER_RU_x2z.inc"
 384 #else
 385
 386      CALL polar_filter_3d( grid%rv_m, grid%clat, .false.,     &
 387                                 fft_filter_lat, dclat,             &
 388                                 ids, ide, jds, jde, kds, kde,       &
 389                                 ims, ime, jms, jme, kms, kme, &
 390                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
 391
 392      k_end = MIN(kde-1,kpe)
 393      IF ( piggyback_mut ) k_end = MIN(kde,kpe)
 394
 395      CALL polar_filter_3d( grid%ru_m, grid%clat, piggyback_mut,     &
 396                                 fft_filter_lat, 0.,                &
 397                                 ids, ide, jds, jde, kds, kde-1,       &
 398                                 ims, ime, jms, jme, kms, kme, &
 399                                 ips, ipe, jps, jpe, kps, k_end )
 400 #endif
 401      IF ( piggyback_mut ) THEN
 402        grid%mut(ips:ipe,jps:jpe) = grid%ru_m(ips:ipe,kde,jps:jpe)
 403        piggyback_mut = .FALSE.
 404      ENDIF
 405    ENDIF
 406
 407 !!!!!!!!!!!!!!!!!!!!!!!
 408 ! MOIST
 409    IF ( flag_moist .GE. PARAM_FIRST_SCALAR ) THEN
 410      itrace = flag_moist
 411 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 412 # include "XPOSE_POLAR_FILTER_MOIST_z2x.inc"
 413      CALL polar_filter_3d( grid%fourd_xxx, grid%clat_xxx, .false. ,     &
 414                            fft_filter_lat, 0.,                &
 415                            ids, ide, jds, jde, kds, kde,       &
 416                            imsx, imex, jmsx, jmex, kmsx, kmex, &
 417                            ipsx, ipex, jpsx, jpex, kpsx, MIN(kpex,kde-1) )
 418
 419      IF ( actual_distance_average ) THEN
 420         CALL filter_tracer ( grid%fourd_xxx , grid%clat_xxx , grid%mf_xxx , &
 421                              grid%fft_filter_lat , grid%mf_fft , &
 422                              pos_def, swap_pole_with_next_j , &
 423                              ids, ide, jds, jde, kds, kde , &
 424                              imsx, imex, jmsx, jmex, kmsx, kmex, &
 425                              ipsx, ipex, jpsx, jpex, kpsx, kpex )
 426      END IF
 427 # include "XPOSE_POLAR_FILTER_MOIST_x2z.inc"
 428 #else
 429      CALL polar_filter_3d( moist(ims,kms,jms,itrace), grid%clat, .false.,     &
 430                            fft_filter_lat, 0.,                &
 431                            ids, ide, jds, jde, kds, kde,       &
 432                            ims, ime, jms, jme, kms, kme, &
 433                            ips, ipe, jps, jpe, kps, MIN(kpe,kde-1) )
 434
 435      IF ( actual_distance_average ) THEN
 436         CALL filter_tracer ( moist(ims,kms,jms,itrace) , grid%clat , grid%msft , &
 437                              grid%fft_filter_lat , grid%mf_fft , &
 438                              pos_def, swap_pole_with_next_j , &
 439                              ids, ide, jds, jde, kds, kde , &
 440                              ims, ime, jms, jme, kms, kme, &
 441                              ips, ipe, jps, jpe, kps, MIN(kpe,kde-1) )
 442      END IF
 443 #endif
 444    ENDIF
 445
 446 !!!!!!!!!!!!!!!!!!!!!!!
 447 ! CHEM
 448    IF ( flag_chem .GE. PARAM_FIRST_SCALAR ) THEN
 449      itrace = flag_chem
 450 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 451 # include "XPOSE_POLAR_FILTER_CHEM_z2x.inc"
 452      CALL polar_filter_3d( grid%fourd_xxx, grid%clat_xxx, .false. ,     &
 453                            fft_filter_lat, 0.,                &
 454                            ids, ide, jds, jde, kds, kde,       &
 455                            imsx, imex, jmsx, jmex, kmsx, kmex, &
 456                            ipsx, ipex, jpsx, jpex, kpsx, MIN(kpex,kde-1) )
 457
 458      IF ( actual_distance_average ) THEN
 459         CALL filter_tracer ( grid%fourd_xxx , grid%clat_xxx , grid%mf_xxx , &
 460                              grid%fft_filter_lat , grid%mf_fft , &
 461                              pos_def, swap_pole_with_next_j , &
 462                              ids, ide, jds, jde, kds, kde , &
 463                              imsx, imex, jmsx, jmex, kmsx, kmex, &
 464                              ipsx, ipex, jpsx, jpex, kpsx, kpex )
 465      END IF
 466 # include "XPOSE_POLAR_FILTER_CHEM_x2z.inc"
 467 #else
 468      CALL polar_filter_3d( chem(ims,kms,jms,itrace), grid%clat, .false. ,     &
 469                            fft_filter_lat, 0.,                &
 470                            ids, ide, jds, jde, kds, kde,       &
 471                            ims, ime, jms, jme, kms, kme, &
 472                            ips, ipe, jps, jpe, kps, MIN(kpe,kde-1) )
 473
 474      IF ( actual_distance_average ) THEN
 475         CALL filter_tracer ( chem(ims,kms,jms,itrace) , grid%clat , grid%msft , &
 476                              grid%fft_filter_lat , grid%mf_fft , &
 477                              pos_def, swap_pole_with_next_j , &
 478                              ids, ide, jds, jde, kds, kde , &
 479                              ims, ime, jms, jme, kms, kme, &
 480                              ips, ipe, jps, jpe, kps, MIN(kpe,kde-1) )
 481      END IF
 482 #endif
 483    ENDIF
 484
 485 !!!!!!!!!!!!!!!!!!!!!!!
 486 ! TRACER
 487    IF ( flag_tracer .GE. PARAM_FIRST_SCALAR ) THEN
 488      itrace = flag_tracer
 489 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 490 # include "XPOSE_POLAR_FILTER_TRACER_z2x.inc"
 491
 492      CALL polar_filter_3d( grid%fourd_xxx, grid%clat_xxx, .false. ,     &
 493                            fft_filter_lat, 0.,                &
 494                            ids, ide, jds, jde, kds, kde,       &
 495                            imsx, imex, jmsx, jmex, kmsx, kmex, &
 496                            ipsx, ipex, jpsx, jpex, kpsx, MIN(kpex,kde-1) )
 497
 498      IF ( actual_distance_average ) THEN
 499         CALL filter_tracer ( grid%fourd_xxx , grid%clat_xxx , grid%mf_xxx , &
 500                              grid%fft_filter_lat , grid%mf_fft , &
 501                              pos_def, swap_pole_with_next_j , &
 502                              ids, ide, jds, jde, kds, kde , &
 503                              imsx, imex, jmsx, jmex, kmsx, kmex, &
 504                              ipsx, ipex, jpsx, jpex, kpsx, kpex )
 505      END IF
 506
 507 # include "XPOSE_POLAR_FILTER_TRACER_x2z.inc"
 508 #else
 509      CALL polar_filter_3d( tracer(ims,kms,jms,itrace), grid%clat, .false. ,     &
 510                            fft_filter_lat, 0.,                &
 511                            ids, ide, jds, jde, kds, kde,       &
 512                            ims, ime, jms, jme, kms, kme, &
 513                            ips, ipe, jps, jpe, kps, MIN(kpe,kde-1) )
 514
 515      IF ( actual_distance_average ) THEN
 516         CALL filter_tracer ( tracer(ims,kms,jms,itrace) , grid%clat , grid%msft , &
 517                              grid%fft_filter_lat , grid%mf_fft , &
 518                              pos_def, swap_pole_with_next_j , &
 519                              ids, ide, jds, jde, kds, kde , &
 520                              ims, ime, jms, jme, kms, kme, &
 521                              ips, ipe, jps, jpe, kps, MIN(kpe,kde-1) )
 522      END IF
 523 #endif
 524    ENDIF
 525
 526 !!!!!!!!!!!!!!!!!!!!!!!
 527 ! SCALAR
 528    IF ( flag_scalar .GE. PARAM_FIRST_SCALAR ) THEN
 529      itrace = flag_scalar
 530 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 531 # include "XPOSE_POLAR_FILTER_SCALAR_z2x.inc"
 532      CALL polar_filter_3d( grid%fourd_xxx , grid%clat_xxx, .false. ,     &
 533                            fft_filter_lat, 0.,                &
 534                            ids, ide, jds, jde, kds, kde,       &
 535                            imsx, imex, jmsx, jmex, kmsx, kmex, &
 536                            ipsx, ipex, jpsx, jpex, kpsx, MIN(kpex,kde-1) )
 537
 538      IF ( actual_distance_average ) THEN
 539         CALL filter_tracer ( grid%fourd_xxx , grid%clat_xxx , grid%mf_xxx , &
 540                              grid%fft_filter_lat , grid%mf_fft , &
 541                              pos_def, swap_pole_with_next_j , &
 542                              ids, ide, jds, jde, kds, kde , &
 543                              imsx, imex, jmsx, jmex, kmsx, kmex, &
 544                              ipsx, ipex, jpsx, jpex, kpsx, kpex )
 545      END IF
 546 # include "XPOSE_POLAR_FILTER_SCALAR_x2z.inc"
 547 #else
 548      CALL polar_filter_3d( scalar(ims,kms,jms,itrace) , grid%clat, .false. ,     &
 549                            fft_filter_lat, 0.,                &
 550                            ids, ide, jds, jde, kds, kde,       &
 551                            ims, ime, jms, jme, kms, kme, &
 552                            ips, ipe, jps, jpe, kps, MIN(kpe,kde-1) )
 553
 554      IF ( actual_distance_average ) THEN
 555         CALL filter_tracer ( scalar(ims,kms,jms,itrace) , grid%clat , grid%msft , &
 556                              grid%fft_filter_lat , grid%mf_fft , &
 557                              pos_def, swap_pole_with_next_j , &
 558                              ids, ide, jds, jde, kds, kde , &
 559                              ims, ime, jms, jme, kms, kme, &
 560                              ips, ipe, jps, jpe, kps, MIN(kpe,kde-1) )
 561      END IF
 562 #endif
 563    ENDIF
 564
 565    IF ( flag_mu .EQ. 1 .AND. piggyback_mu ) THEN
 566       CALL wrf_error_fatal("mu needed to get piggybacked on a transpose and did not")
 567    ENDIF
 568    IF ( flag_mut .EQ. 1 .AND. piggyback_mut ) THEN
 569       CALL wrf_error_fatal("mut needed to get piggybacked on a transpose and did not")
 570    ENDIF
 571
 572 !write(0,*)'pxft back to ',lineno
 573    RETURN
 574 END SUBROUTINE pxft
 575
 576 SUBROUTINE polar_filter_3d( f, xlat, piggyback, fft_filter_lat, dvlat, &
 577                             ids, ide, jds, jde, kds, kde,    &
 578                             ims, ime, jms, jme, kms, kme,    &
 579                             its, ite, jts, jte, kts, kte     )
 580
 581   IMPLICIT NONE
 582
 583   INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 584                                    ims, ime, jms, jme, kms, kme, &
 585                                    its, ite, jts, jte, kts, kte
 586   REAL   ,       INTENT(IN   ) :: fft_filter_lat
 587
 588   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) ::  f
 589   REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) ::  xlat
 590   REAL , INTENT(IN) ::  dvlat
 591   LOGICAL , INTENT(IN) :: piggyback
 592
 593   REAL , DIMENSION(1:ide-ids,1:kte-kts+1) :: sheet
 594   REAL , DIMENSION(1:kte-kts+1) :: sheet_total
 595   REAL :: lat, avg, rnboxw
 596   INTEGER :: ig, jg, i, j, j_end, nx, ny, nmax, kw
 597   INTEGER :: k, nboxw, nbox2, istart, iend, overlap
 598   INTEGER, DIMENSION(6) :: wavenumber = (/ 1, 3, 7, 10, 13, 16 /)
 599
 600   INTEGER :: fftflag
 601
 602   ! Variables will stay in domain form since this routine is meaningless
 603   ! unless tile extent is the same as domain extent in E/W direction, i.e.,
 604   ! the processor has access to all grid points in E/W direction.
 605   ! There may be other ways of doing FFTs, but we have not learned them yet...
 606
 607   ! Check to make sure we have full access to all E/W points
 608   IF ((its /= ids) .OR. (ite /= ide)) THEN
 609      WRITE ( wrf_err_message , * ) 'module_polarfft: 3d: (its /= ids) or (ite /= ide)',its,ids,ite,ide
 610      CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
 611   END IF
 612
 613   fftflag= 1  ! call double-precision fft
 614 ! fftflag= 0  ! call single-precision fft
 615
 616   nx = ide-ids ! "U" stagger variables will be repeated by periodic BCs
 617   ny = kte-kts+1 ! we can filter extra level for variables that are non-Z-staggered
 618   lat = 0.
 619   j_end = MIN(jte, jde-1)
 620   IF (dvlat /= 0. .and. j_end == jde-1) j_end = jde
 621   DO j = jts, j_end
 622      ! jg is the J index in the global (domain) span of the array.
 623      jg = j-jds+1
 624
 625      ! determine whether or not to filter the data
 626
 627      if(xlat(ids,j).gt.0.) then
 628         lat = xlat(ids,j)-dvlat
 629      else
 630         lat = xlat(ids,j)+dvlat
 631      endif
 632
 633      IF (abs(lat) >= fft_filter_lat) THEN
 634
 635         DO k=kts,kte
 636         DO i=ids,ide-1
 637            sheet(i-ids+1,k-kts+1) = f(i,k,j)
 638         END DO
 639         END DO
 640
 641         call polar_filter_fft_2d_ncar(nx,ny,sheet,lat,fft_filter_lat,piggyback,fftflag)
 642
 643         DO k=kts,kte
 644            DO i=ids,ide-1
 645               f(i,k,j) = sheet(i-ids+1,k-kts+1)
 646            END DO
 647            ! setting up ims-ime with x periodicity:
 648            ! enforce periodicity as in set_physical_bc3d
 649
 650            DO i=1,ids-ims
 651               f(ids-i,k,j)=f(ide-i,k,j)
 652            END DO
 653            DO i=1,ime-ide+1
 654               f(ide+i-1,k,j)=f(ids+i-1,k,j)
 655            END DO
 656         END DO
 657
 658      END IF
 659   END DO ! outer j (latitude) loop
 660
 661 END SUBROUTINE polar_filter_3d
 662
 663 !------------------------------------------------------------------------------
 664
 665 SUBROUTINE polar_filter_fft_2d_ncar(nx,ny,fin,lat,filter_latitude,piggyback,fftflag)
 666   IMPLICIT NONE
 667   INTEGER , INTENT(IN) :: nx, ny
 668   REAL , DIMENSION(nx,ny), INTENT(INOUT) :: fin
 669   REAL , INTENT(IN) :: lat, filter_latitude
 670   LOGICAL, INTENT(IN) :: piggyback
 671
 672   REAL :: pi, rcosref, freq, c, cf
 673
 674   INTEGER :: k, fftflag
 675   REAL, DIMENSION(NX,NY) :: work
 676   REAL, DIMENSION(NX+15) :: wsave
 677   REAL(KIND=8), DIMENSION(NX,NY) :: fin8, work8
 678   REAL(KIND=8), DIMENSION(NX+15) :: wsave8
 679
 680   INTEGER :: i, j
 681   REAL, dimension(nx,ny) :: fp
 682
 683   INTEGER :: lensave, ier, nh, n1
 684   INTEGER :: lot, jump, n, inc, lenr, lensav, lenwrk
 685
 686   REAL, PARAMETER :: alpha = 0.0
 687   REAL :: factor_k
 688
 689   INTEGER :: ntop
 690
 691   pi = ACOS(-1.)
 692   rcosref = 1./COS(filter_latitude*pi/180.)
 693
 694 !  we are following the naming convention of the fftpack5 routines
 695
 696   n = nx
 697   lot = ny
 698   lensav = n+15
 699   inc = 1
 700   lenr = nx*ny
 701   jump = nx
 702   lenwrk = lenr
 703   ntop = ny
 704   IF(piggyback) ntop = ny-1
 705
 706 !  forward transform
 707 !  initialize coefficients, place in wsave
 708 !   (should place this in init and save wsave at program start)
 709
 710   if(fftflag.eq.0) then
 711      call rfftmi(n,wsave,lensav,ier)
 712   else
 713      call dfft1i(n,wsave8,lensav,ier)
 714   endif
 715
 716   IF(ier /= 0) THEN
 717     write(a_message,*) ' error in rfftmi ',ier
 718     CALL wrf_message ( a_message )
 719   END IF
 720
 721 !  do the forward transform
 722
 723   if(fftflag.eq.0) then
 724      call rfftmf(lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier)
 725   else
 726      fin8 = fin
 727      do k=1,ny
 728         call dfft1f(n, inc, fin8(1,k), lenr, wsave8, lensav, work8, lenwrk, ier)
 729      enddo
 730   endif
 731
 732   IF(ier /= 0) THEN
 733     write(a_message,*) ' error in rfftmf ',ier
 734     CALL wrf_message ( a_message )
 735   END IF
 736
 737   if(MOD(n,2) == 0) then
 738     nh = n/2 - 1
 739   else
 740     nh = (n-1)/2
 741   end if
 742
 743   DO j=1,ny
 744    fp(1,j) = 1.
 745   ENDDO
 746
 747   DO i=2,nh+1
 748     freq=REAL(i-1)/REAL(n)
 749     c = (rcosref*COS(lat*pi/180.)/SIN(freq*pi))**2
 750
 751 !    c = MAX(0.,MIN(1.,c))
 752     do j=1,ntop
 753       factor_k = (1.-alpha)+alpha*min(1.,float(ntop - j)/10.)
 754       cf = c*factor_k*factor_k
 755       cf = MAX(0.,MIN(1.,cf))
 756       fp(2*(i-1),j) = cf
 757       fp(2*(i-1)+1,j) = cf
 758     enddo
 759     if(piggyback) then
 760       cf = MAX(0.,MIN(1.,c))
 761       fp(2*(i-1),ny) = cf
 762       fp(2*(i-1)+1,ny) = cf
 763     endif
 764   END DO
 765
 766   IF(MOD(n,2) == 0) THEN
 767     c = (rcosref*COS(lat*pi/180.))**2
 768 !    c = MAX(0.,MIN(1.,c))
 769     do j=1,ntop
 770       factor_k = (1.-alpha)+alpha*min(1.,float(ntop - j)/10.)
 771       cf = c*factor_k*factor_k
 772       cf = MAX(0.,MIN(1.,cf))
 773       fp(n,j) = cf
 774     enddo
 775     if(piggyback) then
 776       cf = MAX(0.,MIN(1.,c))
 777       fp(n,ny) = cf
 778     endif
 779   END IF
 780
 781   if(fftflag.eq.0) then
 782      do j=1,ny
 783        do i=1,nx
 784           fin(i,j) = fp(i,j)*fin(i,j)
 785        enddo
 786      enddo
 787   else
 788      do j=1,ny
 789        do i=1,nx
 790           fin8(i,j) = fp(i,j)*fin8(i,j)
 791        enddo
 792      enddo
 793   endif
 794
 795 !  do the backward transform
 796
 797   if(fftflag.eq.0) then
 798      call rfftmb(lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier)
 799   else
 800      do k=1,ny
 801         call dfft1b(n, inc, fin8(1,k), lenr, wsave8, lensav, work8, lenwrk, ier)
 802      enddo
 803      fin= fin8
 804   endif
 805
 806   IF(ier /= 0) THEN
 807     write(a_message,*) ' error in rfftmb ',ier
 808     CALL wrf_message ( a_message )
 809   END IF
 810
 811 END SUBROUTINE polar_filter_fft_2d_ncar
 812
 813 !---------------------------------------------------------------------
 814
 815    SUBROUTINE filter_tracer ( tr3d_in , xlat , msftx , &
 816                               fft_filter_lat , mf_fft , &
 817                               pos_def , swap_pole_with_next_j , &
 818                               ids , ide , jds , jde , kds , kde , &
 819                               ims , ime , jms , jme , kms , kme , &
 820                               its , ite , jts , jte , kts , kte )
 821
 822       IMPLICIT NONE
 823
 824       INTEGER , INTENT(IN) ::  ids , ide , jds , jde , kds , kde , &
 825                                ims , ime , jms , jme , kms , kme , &
 826                                its , ite , jts , jte , kts , kte
 827
 828       REAL , INTENT(IN) :: fft_filter_lat , mf_fft
 829       REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: tr3d_in
 830       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: xlat , msftx
 831       LOGICAL , INTENT(IN) :: pos_def , swap_pole_with_next_j
 832
 833
 834       !  Local vars
 835
 836       INTEGER :: i , j , j_lat_pos , j_lat_neg , k
 837       INTEGER :: i_kicker , ik , i1, i2, i3, i4
 838       INTEGER :: i_left , i_right , ii , count
 839       REAL :: length_scale , sum
 840       REAL , DIMENSION(its:ite,jts:jte) :: tr_in, tr_out
 841       CHARACTER (LEN=256) :: message
 842
 843       !  The filtering is a simple average on a latitude loop.  Possibly a LONG list of
 844       !  numbers.  We assume that ALL of the 2d arrays have been transposed so that
 845       !  each patch has the entire domain size of the i-dimension available.
 846
 847       IF ( ( its .NE. ids ) .OR. ( ite .NE. ide ) ) THEN
 848          CALL wrf_error_fatal ( 'filtering assumes all values on X' )
 849       END IF
 850
 851       !  Starting at the south pole, we find where the
 852       !  grid distance is big enough, then go back a point.  Continuing to the
 853       !  north pole, we find the first small grid distance.  These are the
 854       !  computational latitude loops and the associated computational poles.
 855
 856       j_lat_neg = 0
 857       j_lat_pos = jde + 1
 858       loop_neg : DO j = MIN(jde-1,jte) , jts , -1
 859          IF ( xlat(its,j) .LT. 0.0 ) THEN
 860             IF ( ABS(xlat(its,j)) .GE. fft_filter_lat ) THEN
 861                j_lat_neg = j
 862                EXIT loop_neg
 863             END IF
 864          END IF
 865       END DO loop_neg
 866
 867       loop_pos : DO j = jts , MIN(jde-1,jte)
 868          IF ( xlat(its,j) .GT. 0.0 ) THEN
 869             IF ( xlat(its,j) .GE. fft_filter_lat ) THEN
 870                j_lat_pos = j
 871                EXIT loop_pos
 872             END IF
 873          END IF
 874       END DO loop_pos
 875
 876       !  Initialize the starting values for the averages.
 877
 878       DO k = kts, kte
 879          DO j = jts , MIN(jde-1,jte)
 880             DO i = its , MIN(ide-1,ite)
 881                tr_in(i,j) = tr3d_in(i,k,j)
 882                tr_out(i,j) = tr_in(i,j)
 883             END DO
 884          END DO
 885
 886          !  Filter the fields at the negative lats.
 887
 888          DO j = MIN(j_lat_neg,jte) , jts , -1
 889 !           i_kicker = MIN( MAX ( NINT(msftx(its,j)/2) , 1 ) , (ide - ids) / 2 )
 890             i_kicker = MIN( MAX ( NINT(msftx(its,j)/mf_fft/2) , 1 ) , (ide - ids) / 2 )
 891 !           WRITE (message,FMT='(a,i4,a,i4,f6.2,1x,f6.1,f6.1)') 'SOUTH j = ' , j, ', kicker = ',i_kicker,xlat(its,j),msftx(its,j),mf_fft
 892 !           CALL wrf_debug ( 0 , TRIM(message) )
 893             DO i = its , MIN(ide-1,ite)
 894                sum = 0.
 895                count = 0
 896                DO ik = 1 , i_kicker/2
 897                   ii = i-ik
 898                   IF ( ii .GE. ids ) THEN
 899                      i_left = ii
 900                   ELSE
 901                      i_left = ( ii - ids ) + (ide-1)+1
 902                   END IF
 903                   ii = i+ik
 904                   IF ( ii .LE. ide-1 ) THEN
 905                      i_right = ii
 906                   ELSE
 907                      i_right = ( ii - (ide-1) ) + its-1
 908                   END IF
 909                   sum = sum + tr_in(i_left,j) + tr_in(i_right,j)
 910                   count = count + 1
 911                END DO
 912                tr_out(i,j) = ( tr_in(i,j) + sum ) / REAL ( 2 * count + 1 )
 913             END DO
 914          END DO
 915
 916          !  Filter the fields at the positive lats.
 917
 918          DO j = MAX(j_lat_pos,jts) , MIN(jde-1,jte)
 919 !           i_kicker = MIN( MAX ( NINT(msftx(its,j)/2) , 1 ) , (ide - ids) / 2 )
 920             i_kicker = MIN( MAX ( NINT(msftx(its,j)/mf_fft/2) , 1 ) , (ide - ids) / 2 )
 921 !           WRITE (message,FMT='(a,i4,a,i4,f6.2,1x,f6.1,f6.1)') 'NORTH j = ' , j, ', kicker = ',i_kicker,xlat(its,j),msftx(its,j),mf_fft
 922 !           CALL wrf_debug ( 0 , TRIM(message) )
 923             DO i = its , MIN(ide-1,ite)
 924                count = 0
 925                sum = 0.
 926                DO ik = 1 , i_kicker/2
 927                   ii = i-ik
 928                   IF ( ii .GE. ids ) THEN
 929                      i_left = ii
 930                   ELSE
 931                      i_left = ( ii - ids ) + (ide-1)+1
 932                   END IF
 933                   ii = i+ik
 934                   IF ( ii .LE. ide-1 ) THEN
 935                      i_right = ii
 936                   ELSE
 937                      i_right = ( ii - (ide-1) ) + its-1
 938                   END IF
 939                   sum = sum + tr_in(i_left,j) + tr_in(i_right,j)
 940                   count = count + 1
 941                END DO
 942                tr_out(i,j) = ( tr_in(i,j) + sum ) / REAL ( 2 * count + 1 )
 943             END DO
 944          END DO
 945
 946          !  Set output values for whole patch.
 947
 948          DO j = jts , MIN(jde-1,jte)
 949             DO i = its , MIN(ide-1,ite)
 950                tr3d_in(i,k,j) = tr_out(i,j)
 951             END DO
 952          END DO
 953
 954          !  Positive definite on scalars?
 955
 956          IF ( pos_def ) THEN
 957             DO j = jts , MIN(jde-1,jte)
 958                DO i = its , MIN(ide-1,ite)
 959                   tr3d_in(i,k,j) = MAX( tr3d_in(i,k,j) , 0. )
 960                END DO
 961             END DO
 962          END IF
 963
 964          !  Remove values at j=1 and j=jde-1 locations, set them to the rows just next to them.
 965
 966          IF ( swap_pole_with_next_j ) THEN
 967             IF  ( jts .EQ. jds ) THEN
 968                DO i = its , MIN(ide-1,ite)
 969                   tr3d_in(i,k,jts) = tr3d_in(i,k,jts+1)
 970                END DO
 971             END IF
 972
 973             IF  ( jte .EQ. jde ) THEN
 974                DO i = its , MIN(ide-1,ite)
 975                   tr3d_in(i,k,jte-1) = tr3d_in(i,k,jte-2)
 976                END DO
 977             END IF
 978          END IF
 979
 980       END DO ! k-loop
 981
 982    END SUBROUTINE filter_tracer
 983
 984 !---------------------------------------------------------------------
 985
 986    SUBROUTINE filter_tracer_old ( tr3d_in , xlat , msftx , fft_filter_lat , &
 987                               ids , ide , jds , jde , kds , kde , &
 988                               ims , ime , jms , jme , kms , kme , &
 989                               its , ite , jts , jte , kts , kte )
 990
 991       IMPLICIT NONE
 992
 993       INTEGER , INTENT(IN) ::  ids , ide , jds , jde , kds , kde , &
 994                                ims , ime , jms , jme , kms , kme , &
 995                                its , ite , jts , jte , kts , kte
 996
 997       REAL , INTENT(IN) :: fft_filter_lat
 998       REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: tr3d_in
 999       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: xlat , msftx
1000
1001
1002       !  Local vars
1003
1004       INTEGER :: i , j , j_lat_pos , j_lat_neg , k
1005       INTEGER :: i_kicker , ik , i1, i2, i3, i4
1006       REAL :: length_scale , sum
1007       REAL , DIMENSION(its:ite,jts:jte) :: tr_in, tr_out
1008
1009       !  The filtering is a simple average on a latitude loop.  Possibly a LONG list of
1010       !  numbers.  We assume that ALL of the 2d arrays have been transposed so that
1011       !  each patch has the entire domain size of the i-dim local.
1012
1013       IF ( ( its .NE. ids ) .OR. ( ite .NE. ide ) ) THEN
1014          CALL wrf_error_fatal ( 'filtering assumes all values on X' )
1015       END IF
1016
1017       !  Starting at the south pole, we find where the
1018       !  grid distance is big enough, then go back a point.  Continuing to the
1019       !  north pole, we find the first small grid distance.  These are the
1020       !  computational latitude loops and the associated computational poles.
1021
1022       j_lat_neg = 0
1023       j_lat_pos = jde + 1
1024       loop_neg : DO j = jts , MIN(jde-1,jte)
1025          IF ( xlat(its,j) .LT. 0.0 ) THEN
1026             IF ( ABS(xlat(its,j)) .LT. fft_filter_lat ) THEN
1027                j_lat_neg = j - 1
1028                EXIT loop_neg
1029             END IF
1030          END IF
1031       END DO loop_neg
1032
1033       loop_pos : DO j = jts , MIN(jde-1,jte)
1034          IF ( xlat(its,j) .GT. 0.0 ) THEN
1035             IF ( xlat(its,j) .GE. fft_filter_lat ) THEN
1036                j_lat_pos = j
1037                EXIT loop_pos
1038             END IF
1039          END IF
1040       END DO loop_pos
1041
1042       !  Set output values to initial input topo values for whole patch.
1043
1044       DO k = kts, kte
1045          DO j = jts , MIN(jde-1,jte)
1046             DO i = its , MIN(ide-1,ite)
1047                tr_in(i,j) = tr3d_in(i,k,j)
1048                tr_out(i,j) = tr_in(i,j)
1049             END DO
1050          END DO
1051
1052          !  Filter the topo at the negative lats.
1053
1054          DO j = j_lat_neg , jts , -1
1055             i_kicker = MIN( MAX ( NINT(msftx(its,j)) , 1 ) , (ide - ids) / 2 )
1056 !           print *,'j = ' , j, ', kicker = ',i_kicker
1057             DO i = its , MIN(ide-1,ite)
1058                IF      ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
1059                   sum = 0.0
1060                   DO ik = 1 , i_kicker
1061                      sum = sum + tr_in(i+ik,j) + tr_in(i-ik,j)
1062                   END DO
1063                   tr_out(i,j) = ( tr_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
1064                ELSE IF ( ( i - i_kicker .LT. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
1065                   sum = 0.0
1066                   DO ik = 1 , i_kicker
1067                      sum = sum + tr_in(i+ik,j)
1068                   END DO
1069                   i1 = i - i_kicker + ide -1
1070                   i2 = ide-1
1071                   i3 = ids
1072                   i4 = i-1
1073                   DO ik = i1 , i2
1074                      sum = sum + tr_in(ik,j)
1075                   END DO
1076                   DO ik = i3 , i4
1077                      sum = sum + tr_in(ik,j)
1078                   END DO
1079                   tr_out(i,j) = ( tr_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
1080                ELSE IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .GT. ide-1 ) ) THEN
1081                   sum = 0.0
1082                   DO ik = 1 , i_kicker
1083                      sum = sum + tr_in(i-ik,j)
1084                   END DO
1085                   i1 = i+1
1086                   i2 = ide-1
1087                   i3 = ids
1088                   i4 = ids + ( i_kicker+i ) - ide
1089                   DO ik = i1 , i2
1090                      sum = sum + tr_in(ik,j)
1091                   END DO
1092                   DO ik = i3 , i4
1093                      sum = sum + tr_in(ik,j)
1094                   END DO
1095                   tr_out(i,j) = ( tr_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
1096                END IF
1097             END DO
1098          END DO
1099
1100          !  Filter the topo at the positive lats.
1101
1102          DO j = j_lat_pos , MIN(jde-1,jte)
1103             i_kicker = MIN( MAX ( NINT(msftx(its,j)) , 1 ) , (ide - ids) / 2 )
1104 !           print *,'j = ' , j, ', kicker = ',i_kicker
1105             DO i = its , MIN(ide-1,ite)
1106                IF      ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
1107                   sum = 0.0
1108                   DO ik = 1 , i_kicker
1109                      sum = sum + tr_in(i+ik,j) + tr_in(i-ik,j)
1110                   END DO
1111                   tr_out(i,j) = ( tr_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
1112                ELSE IF ( ( i - i_kicker .LT. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
1113                   sum = 0.0
1114                   DO ik = 1 , i_kicker
1115                      sum = sum + tr_in(i+ik,j)
1116                   END DO
1117                   i1 = i - i_kicker + ide -1
1118                   i2 = ide-1
1119                   i3 = ids
1120                   i4 = i-1
1121                   DO ik = i1 , i2
1122                      sum = sum + tr_in(ik,j)
1123                   END DO
1124                   DO ik = i3 , i4
1125                      sum = sum + tr_in(ik,j)
1126                   END DO
1127                   tr_out(i,j) = ( tr_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
1128                ELSE IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .GT. ide-1 ) ) THEN
1129                   sum = 0.0
1130                   DO ik = 1 , i_kicker
1131                      sum = sum + tr_in(i-ik,j)
1132                   END DO
1133                   i1 = i+1
1134                   i2 = ide-1
1135                   i3 = ids
1136                   i4 = ids + ( i_kicker+i ) - ide
1137                   DO ik = i1 , i2
1138                      sum = sum + tr_in(ik,j)
1139                   END DO
1140                   DO ik = i3 , i4
1141                      sum = sum + tr_in(ik,j)
1142                   END DO
1143                   tr_out(i,j) = ( tr_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
1144                END IF
1145             END DO
1146          END DO
1147
1148          !  Set output values to initial input topo values for whole patch.
1149
1150          DO j = jts , MIN(jde-1,jte)
1151             DO i = its , MIN(ide-1,ite)
1152                tr3d_in(i,k,j) = tr_out(i,j)
1153             END DO
1154          END DO
1155       END DO ! k-loop
1156
1157    END SUBROUTINE filter_tracer_old
1158
1159 !---------------------------------------------------------------------
1160 END MODULE module_polarfft
1161