phys/module_fr_fire_util.F

   1 !
   2 ! This module contains general purpose utilities and WRF wrappers because some
   3 ! applications require this module to run standalone. No physics here.
   4 ! Some are dependent on WRF indexing scheme. Some violate WRF conventions but these
   5 ! are not called from the WRF dependent code. Some are not called at all.
   6 !
   7
   8 module module_fr_fire_util
   9
  10 ! various method selection parameters
  11 ! 1. add the parameter and its static default here
  12 ! optional:
  13 ! 2. add copy from config_flags in module_fr_fire_driver%%set_flags
  14 ! 3. add a line in Registry.EM to define the variable and set default value
  15 ! 4. add the variable and value in namelist.input
  16
  17
  18 integer,save::          &
  19  fire_print_msg=2,      &  ! print FIRE progress
  20  fire_print_file=1,     &  ! write many files by write_array_m; compile with DEBUG_OUT, do not run in parallel
  21  fuel_left_method=1,    &  ! 1=simple, 2=exact in linear case
  22  fuel_left_irl=2,       &  ! refinement for fuel calculation, must be even
  23  fuel_left_jrl=2,       &
  24  boundary_guard=-1,     &  ! crash if fire gets this many cells to domain boundary, -1=off
  25  fire_grows_only=1,     &  ! fire can spread out only (level set functions may not increase)
  26  fire_upwinding=9,      &  ! upwind normal spread: 1=standard, 2=godunov, 3=eno, 4=sethian, 5=2nd-order, 6=WENO3, 7=WENO5, 8=hybrid WENO3/ENO1, 9=hybrid WENO5/ENO1
  27  fire_upwind_split=0,   &  ! 1=upwind advection separately from normal direction spread
  28  fire_test_steps=0,     &  ! 0=no fire, 1=normal, >1 = do specified number of steps and terminate (testing only)
  29  fire_topo_from_atm=1,  &  ! 0 = expect ZSF set correctly on entry, 1 = populate by interploating from atmosphere
  30  fire_advection=0,      &  ! 0 = fire spread from normal wind/slope (CAWFE), 1 = full speed projected
  31  fire_fmc_read=1,       &  ! fuel moisture: 0 from wrfinput, 1 from namelist.fire, 2 read from file in ideal
  32  fire_upwinding_reinit=4, &  ! spatial discretization for reinitialization PDE: 1=WENO3, 2=WENO5, 3=hybrid WENO3-ENO1, 4=hybrid WENO5-ENO1
  33  fire_lsm_reinit_iter=1,  &  ! number of iterations for the reinitialization PDE
  34  fire_lsm_band_ngp=4,     &  ! number of grid points around lfn=0 that the higher-order scheme WENO5/WENO3 is used (ENO1 elsewhere), for ire_upwinding=8,9 and fire_upwinding_reinit=3,4 and
  35  fire_viscosity_ngp=4        ! number of grid points around lfn=0 where low artificial viscosity is used = fire_viscosity_bg
  36
  37 real, save::            &
  38  fire_const_time=-1,    &  ! time from ignition to start constant heat output  <0=never
  39  fire_const_grnhfx=-1., &  ! if both >=0, the constant heat flux to output then
  40  fire_const_grnqfx=-1., &  ! if both >=0, the constant heat flux to output then
  41  fire_atm_feedback=1. , &  ! 1 = normal, 0. = one way coupling atmosphere -> fire only
  42  fire_viscosity=0.4,    &  ! artificial viscosity
  43  fire_lfn_ext_up=1.,    &  ! 0.=extend level set function at boundary by reflection, 1.=always up
  44  fire_lsm_zcoupling_ref=50., &  ! Reference height from wich u at fire_wind_hegiht is calculated using a logarithmic profile
  45  fire_viscosity_bg=0.4,      &  ! artificial viscosity in the near-front region
  46  fire_viscosity_band=0.5,    &  ! number of times the hybrid advection band to transition from fire_viscosity_bg to fire_viscosity
  47  fire_slope_factor=1.0          ! slope correction factor
  48
  49 logical,save::              &
  50  fire_lsm_reinit=.true.,    &  ! flag to activate reinitialization of level set method
  51  fire_lsm_zcoupling=.false.     ! flag to activate reference velocity at a different height from fire_wind_height
  52
  53 integer, parameter:: REAL_SUM=10, REAL_MAX=20, RNRM_SUM=30, RNRM_MAX=40
  54
  55 contains
  56
  57
  58 subroutine interpolate_z2fire(id,                 & ! for debug output, <= 0 no output
  59     ids,ide, jds,jde,                    & ! atm grid dimensions
  60     ims,ime, jms,jme,                    &
  61     ips,ipe,jps,jpe,                              &
  62     its,ite,jts,jte,                              &
  63     ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
  64     ifms, ifme, jfms, jfme,                       &
  65     ifts,ifte,jfts,jfte,                          &
  66     ir,jr,                                        & ! atm/fire grid ratio
  67     zs,                                       & ! atm grid arrays in
  68     zsf,flag_z0)                                      ! fire grid arrays out
  69
  70 implicit none
  71 !*** purpose: interpolate height or any other 2d variable defined at mesh cell centers
  72
  73 !*** arguments
  74 integer, intent(in)::id,                          &
  75     ids,ide, jds,jde,                    & ! atm domain bounds
  76     ims,ime,jms,jme,                    & ! atm memory bounds
  77     ips,ipe,jps,jpe,                              &
  78     its,ite,jts,jte,                              & ! atm tile bounds
  79     ifds, ifde, jfds, jfde,                       & ! fire domain bounds
  80     ifms, ifme, jfms, jfme,                       & ! fire memory bounds
  81     ifts,ifte,jfts,jfte,                          & ! fire tile bounds
  82     ir,jr                                         ! atm/fire grid refinement ratio
  83 real, intent(in), dimension(ims:ime, jms:jme):: zs  ! terrain height at atm cell centers                                        & ! terrain height
  84 real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
  85     zsf                                             ! terrain height fire grid nodes
  86 integer,intent(in)::flag_z0
  87
  88
  89 !*** local
  90 real, dimension(its-2:ite+2,jts-2:jte+2):: za      ! terrain height
  91 integer:: i,j,jts1,jte1,its1,ite1,jfts1,jfte1,ifts1,ifte1,itso,jtso,iteo,jteo
  92
  93 ! terrain height
  94
  95     jts1=max(jts-1,jds) ! lower loop limit by one less when at end of domain
  96     its1=max(its-1,ids) ! ASSUMES THE HALO IS THERE if patch != domain
  97     jte1=min(jte+1,jde)
  98     ite1=min(ite+1,ide)
  99     do j = jts1,jte1
 100         do i = its1,ite1
 101             ! copy to local array
 102             za(i,j)=zs(i,j)
 103         enddo
 104     enddo
 105
 106     call continue_at_boundary(1,1,0., & ! do x direction or y direction
 107     its-2,ite+2,jts-2,jte+2,           &                ! memory dims
 108     ids,ide,jds,jde, &            ! domain dims - winds defined up to +1
 109     ips,ipe,jps,jpe, &            ! patch dims - winds defined up to +1
 110     its1,ite1,jts1,jte1, &                ! tile dims
 111     itso,jtso,iteo,jteo, &
 112     za)                               ! array
 113
 114     ! interpolate to tile plus strip along domain boundary if at boundary
 115     jfts1=snode(jfts,jfds,-1) ! lower loop limit by one less when at end of domain
 116     ifts1=snode(ifts,ifds,-1)
 117     jfte1=snode(jfte,jfde,+1)
 118     ifte1=snode(ifte,ifde,+1)
 119
 120     call interpolate_2d(  &
 121         its-2,ite+2,jts-2,jte+2, & ! memory dims atm grid tile
 122         its1-1,ite1+1,jts1-1,jte1+1, & ! where atm grid values set
 123         ifms,ifme,jfms,jfme,    & ! array dims fire grid
 124         ifts1,ifte1,jfts1,jfte1,  & ! dimensions fire grid tile
 125         ir,jr,                  & ! refinement ratio
 126         real(ids),real(jds),ifds+(ir-1)*0.5,jfds+(jr-1)*0.5, & ! line up by lower left corner of domain
 127         za,                     & ! in atm grid
 128         zsf)                      ! out fire grid
 129
 130 if (flag_z0 .eq. 1 ) then
 131   do j=jfts1,jfte1
 132     do i=ifts1,ifte1
 133       zsf(i,j)=max(zsf(i,j),0.001)
 134     enddo
 135   enddo
 136 endif
 137
 138 end subroutine interpolate_z2fire
 139
 140 !
 141 !****************
 142 !
 143
 144 subroutine crash(s)
 145 use module_wrf_error
 146 implicit none
 147 character(len=*), intent(in)::s
 148 character(len=128)msg
 149 msg='crash: '//s
 150 call message(msg,level=0)
 151 !$OMP CRITICAL(FIRE_MESSAGE_CRIT)
 152 call wrf_error_fatal(msg)
 153 !$OMP END CRITICAL(FIRE_MESSAGE_CRIT)
 154 end subroutine crash
 155
 156 !
 157 !****************
 158 !
 159
 160 subroutine message(s,level)
 161 use module_wrf_error
 162 #ifdef _OPENMP
 163 use OMP_LIB
 164 #endif
 165 implicit none
 166 !*** arguments
 167 character(len=*), intent(in)::s
 168 integer,intent(in),optional::level
 169 !*** local
 170 character(len=128)::msg
 171 character(len=118)::t
 172 integer m,mlevel
 173 logical op
 174 !*** executable
 175 if(present(level))then
 176     mlevel=level
 177 else
 178     mlevel=2  ! default message level
 179 endif
 180 if(fire_print_msg.ge.mlevel)then
 181       m=0
 182 !$OMP CRITICAL(FIRE_MESSAGE_CRIT)
 183 #ifdef _OPENMP
 184       m=omp_get_thread_num()
 185       t=s
 186       write(msg,'(a6,i3,a1,a118)')'FIRE:',m,':',t
 187 #else
 188       msg='FIRE:'//s
 189 #endif
 190       call wrf_message(msg)
 191 !      flush(6)
 192 !      flush(0)
 193 !$OMP END CRITICAL(FIRE_MESSAGE_CRIT)
 194 endif
 195 end subroutine message
 196
 197 !
 198 !****************
 199 !
 200
 201
 202 integer function open_text_file(filename,rw,allow_fail)
 203 implicit none
 204 character(len=*),intent(in):: filename,rw
 205 logical, intent(in), optional:: allow_fail
 206 #ifdef _OPENMP
 207 !integer, external:: OMP_GET_THREAD_NUM
 208 #endif
 209 character(len=128):: msg
 210 character(len=1)::act
 211 integer::iounit,ierr
 212 logical::op,ok
 213
 214 #ifdef _OPENMP
 215 !   if (OMP_GET_THREAD_NUM() .ne. 0)then
 216 !      call crash('open_input_text_file: called from parallel loop')
 217 !   endif
 218 #endif
 219 if(present(allow_fail))then
 220     ok=allow_fail
 221 else
 222     ok=.false.  ! default
 223 endif
 224
 225     do iounit=19,99
 226        inquire(iounit,opened=op)
 227        if(.not.op)goto 1
 228     enddo
 229     call crash('open_text_file: Cannot find any available I/O unit')
 230 1   continue
 231     act=rw(1:1)
 232     select case (act)
 233         case ('r','R')
 234             OPEN(iounit, FILE=filename,FORM='FORMATTED',STATUS='OLD',ACTION='READ',IOSTAT=ierr)
 235         case ('w','W')
 236             OPEN(iounit, FILE=filename,FORM='FORMATTED',STATUS='UNKNOWN',ACTION='WRITE',IOSTAT=ierr)
 237         case default
 238             write(msg,*)'open_text_file: bad mode ',trim(rw),' for file ',trim(filename)
 239     end select
 240
 241     if(ierr.ne.0)then
 242         write(msg,*)'Cannot open file ',filename
 243         if(ok)then
 244             call message(msg)
 245             open_text_file=-1
 246         else
 247             call crash(msg)
 248         endif
 249     else
 250         open_text_file=iounit
 251     endif
 252
 253 end function open_text_file
 254
 255 !
 256 !****************
 257 !
 258
 259
 260 subroutine set_ideal_coord( dxf,dyf, &
 261                 ifds,ifde,jfds,jfde,  &
 262                 ifms,ifme,jfms,jfme,  &
 263                 ifts,ifte,jfts,jfte,  &
 264                 fxlong,fxlat          &
 265             )
 266 implicit none
 267 ! arguments
 268 real, intent(in)::dxf,dyf
 269 integer, intent(in):: &
 270                 ifds,ifde,jfds,jfde,  &
 271                 ifms,ifme,jfms,jfme,  &
 272                 ifts,ifte,jfts,jfte
 273 real, intent(out),dimension(ifms:ifme,jfms:jfme)::fxlong,fxlat
 274 ! local
 275 integer::i,j
 276                 ! could we not just as easily get something that
 277                 ! that looks like a lat/lon
 278                 ! set fake  coordinates, in m
 279                 do j=jfts,jfte
 280                     do i=ifts,ifte
 281                         ! uniform mesh, lower left domain corner is (0,0)
 282                         fxlong(i,j)=(i-ifds+0.5)*dxf
 283                         fxlat (i,j)=(j-jfds+0.5)*dyf
 284                     enddo
 285                 enddo
 286 end subroutine set_ideal_coord
 287
 288 !
 289 !****************
 290 !
 291
 292
 293 subroutine continue_at_boundary(ix,iy,bias, & ! do x direction or y direction
 294     ims,ime,jms,jme, &                ! memory dims
 295     ids,ide,jds,jde, &                ! domain dims
 296     ips,ipe,jps,jpe, &                ! patch dims
 297     its,ite,jts,jte, &                ! tile dims
 298     itso,iteo,jtso,jteo, &            ! tile dims where set
 299     lfn)                              ! array
 300 implicit none
 301 !*** description
 302 ! extend array by one beyond the domain by linear continuation
 303 !*** arguments
 304 integer, intent(in)::ix,iy              ! not 0 = do x or y (1 or 2) direction
 305 real,intent(in)::bias                   ! 0=none, 1.=max
 306 integer, intent(in)::ims,ime,jms,jme, &                ! memory dims
 307     ids,ide,jds,jde, &                ! domain dims
 308     ips,ipe,jps,jpe, &                ! patch dims
 309     its,ite,jts,jte                   ! tile dims
 310 integer, intent(out)::itso,jtso,iteo,jteo ! where set
 311 real,intent(inout),dimension(ims:ime,jms:jme)::lfn
 312 !*** local
 313 integer i,j
 314 character(len=128)::msg
 315 integer::its1,ite1,jts1,jte1
 316 integer,parameter::halo=1           ! only 1 domain halo is needed since ENO1 is used near domain boundaries
 317 !*** executable
 318
 319 ! for dislay only
 320 itso=its
 321 jtso=jts
 322 iteo=ite
 323 jteo=jte
 324 ! go halo width beyond if at patch boundary but not at domain boudnary
 325 ! assume we have halo need to compute the value we do not have
 326 ! the next thread that would conveniently computer the extended values at patch corners
 327 ! besides halo may not transfer values outside of the domain
 328 !
 329 its1=its
 330 jts1=jts
 331 ite1=ite
 332 jte1=jte
 333 if(its.eq.ips.and..not.its.eq.ids)its1=its-halo
 334 if(jts.eq.jps.and..not.jts.eq.jds)jts1=jts-halo
 335 if(ite.eq.ipe.and..not.ite.eq.ide)ite1=ite+halo
 336 if(jte.eq.jpe.and..not.jte.eq.jde)jte1=jte+halo
 337 !$OMP CRITICAL(FIRE_UTIL_CRIT)
 338 write(msg,'(a,2i5,a,f5.2)')'continue_at_boundary: directions',ix,iy,' bias ',bias
 339 call message(msg)
 340 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
 341 if(ix.ne.0)then
 342     if(its.eq.ids)then
 343         do j=jts1,jte1
 344             lfn(ids-1,j)=EX(lfn(ids,j),lfn(ids+1,j))
 345         enddo
 346         itso=ids-1
 347     endif
 348     if(ite.eq.ide)then
 349         do j=jts1,jte1
 350             lfn(ide+1,j)=EX(lfn(ide,j),lfn(ide-1,j))
 351         enddo
 352         iteo=ide+1
 353     endif
 354 !$OMP CRITICAL(FIRE_UTIL_CRIT)
 355     write(msg,'(8(a,i5))')'continue_at_boundary: x:',its,':',ite,',',jts,':',jte,' ->',itso,':',iteo,',',jts1,':',jte1
 356     call message(msg)
 357 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
 358 endif
 359 if(iy.ne.0)then
 360     if(jts.eq.jds)then
 361         do i=its1,ite1
 362             lfn(i,jds-1)=EX(lfn(i,jds),lfn(i,jds+1))
 363         enddo
 364         jtso=jds-1
 365     endif
 366     if(jte.eq.jde)then
 367         do i=its1,ite1
 368             lfn(i,jde+1)=EX(lfn(i,jde),lfn(i,jde-1))
 369         enddo
 370         jteo=jde+1
 371     endif
 372 !$OMP CRITICAL(FIRE_UTIL_CRIT)
 373     write(msg,'(8(a,i5))')'continue_at_boundary: y:',its,':',ite,',',jts,':',jte,' ->',its1,':',ite1,',',jtso,':',jteo
 374 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
 375     call message(msg)
 376 endif
 377 ! corners of the domain
 378 if(ix.ne.0.and.iy.ne.0)then
 379     if(its.eq.ids.and.jts.eq.jds)lfn(ids-1,jds-1)=EX(lfn(ids,jds),lfn(ids+1,jds+1))
 380     if(its.eq.ids.and.jte.eq.jde)lfn(ids-1,jde+1)=EX(lfn(ids,jde),lfn(ids+1,jde-1))
 381     if(ite.eq.ide.and.jts.eq.jds)lfn(ide+1,jds-1)=EX(lfn(ide,jds),lfn(ide-1,jds+1))
 382     if(ite.eq.ide.and.jte.eq.jde)lfn(ide+1,jde+1)=EX(lfn(ide,jde),lfn(ide-1,jde-1))
 383 endif
 384 return
 385 contains
 386 real function EX(a,b)
 387 !*** statement function
 388 real a,b
 389 EX=(1.-bias)*(2.*a-b)+bias*max(2.*a-b,a,b)   ! extrapolation, max quarded
 390 end function EX
 391 end subroutine continue_at_boundary
 392
 393 !
 394 !*****************************
 395 !
 396
 397 subroutine check_mesh_2dim(ids,ide,jds,jde,ims,ime,jms,jme)
 398 implicit none
 399 integer, intent(in)::ids,ide,jds,jde,ims,ime,jms,jme
 400 character(len=128)msg
 401 if(ids<ims.or.ide>ime.or.jds<jms.or.jde>jme)then
 402 !$OMP CRITICAL(FIRE_UTIL_CRIT)
 403     write(msg,*)'mesh dimensions:  ',ids,ide,jds,jde
 404     call message(msg)
 405     write(msg,*)'memory dimensions:',ims,ime,jms,jme
 406 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
 407     call message(msg)
 408     call crash('check_mesh_2dim: memory dimensions too small')
 409 endif
 410 end subroutine check_mesh_2dim
 411
 412 !
 413 !****************
 414 !
 415
 416 subroutine check_mesh_3dim(ids,ide,kds,kde,jds,jde,ims,ime,kms,kme,jms,jme)
 417 integer, intent(in)::ids,ide,jds,jde,ims,ime,jms,jme,kds,kde,kms,kme
 418 if(ids<ims.or.ide>ime.or.jds<jms.or.jde>jme.or.kds<kms.or.kde>kme) then
 419     call crash('memory dimensions too small')
 420 endif
 421 end subroutine check_mesh_3dim
 422
 423 !
 424 !****************
 425 !
 426
 427 subroutine sum_2d_cells(       &
 428        ims2,ime2,jms2,jme2,    &
 429        its2,ite2,jts2,jte2,    &
 430        v2,                     &  ! input
 431        ims1,ime1,jms1,jme1,    &
 432        its1,ite1,jts1,jte1,    &
 433        v1)                        ! output
 434 implicit none
 435
 436 !*** purpose
 437 ! sum cell values in mesh2 to cell values of coarser mesh1
 438
 439 !*** arguments
 440 ! the dimensions are in cells, not nodes!
 441
 442 integer, intent(in)::its1,ite1,jts1,jte1,ims1,ime1,jms1,jme1
 443 real, intent(out)::v1(ims1:ime1,jms1:jme1)
 444 integer, intent(in)::its2,ite2,jts2,jte2,ims2,ime2,jms2,jme2
 445 real, intent(in)::v2(ims2:ime2,jms2:jme2)
 446
 447 !*** local
 448 integer:: i1,i2,j1,j2,ir,jr,isz1,isz2,jsz1,jsz2,ioff,joff,ibase,jbase
 449 real t
 450 character(len=128)msg
 451
 452 !*** executable
 453
 454 !check mesh dimensions and domain dimensions
 455 call check_mesh_2dim(its1,ite1,jts1,jte1,ims1,ime1,jms1,jme1)
 456 call check_mesh_2dim(its2,ite2,jts2,jte2,ims2,ime2,jms2,jme2)
 457
 458 ! compute mesh sizes
 459 isz1 = ite1-its1+1
 460 jsz1 = jte1-jts1+1
 461 isz2 = ite2-its2+1
 462 jsz2 = jte2-jts2+1
 463
 464 ! check mesh sizes
 465 if(isz1.le.0.or.jsz1.le.0.or.isz2.le.0.or.jsz2.le.0)then
 466     call message('all mesh sizes must be positive')
 467     goto 9
 468 endif
 469
 470 ! compute mesh ratios
 471 ir=isz2/isz1
 472 jr=jsz2/jsz1
 473
 474 if(isz2.ne.isz1*ir .or. jsz2.ne.jsz1*jr)then
 475     call message('input mesh size must be multiple of output mesh size')
 476     goto 9
 477 endif
 478
 479
 480 ! v1 = sum(v2)
 481 do j1=jts1,jte1
 482     jbase=jts2+jr*(j1-jts1)
 483     do i1=its1,ite1
 484        ibase=its2+ir*(i1-its1)
 485        t=0.
 486        do joff=0,jr-1
 487            j2=joff+jbase
 488            do ioff=0,ir-1
 489                i2=ioff+ibase
 490                t=t+v2(i2,j2)
 491            enddo
 492        enddo
 493        v1(i1,j1)=t
 494     enddo
 495 enddo
 496
 497 return
 498
 499 9 continue
 500 !$OMP CRITICAL(FIRE_UTIL_CRIT)
 501 write(msg,91)its2,ite2,jts2,jte2,ims2,ime2,jms2,jme2
 502 call message(msg)
 503 write(msg,91)its1,ite1,jts1,jte1,ims1,ime1,jms1,jme1
 504 call message(msg)
 505 write(msg,92)'input  mesh size:',isz2,jsz2
 506 call message(msg)
 507 91 format('dimensions: ',8i8)
 508 write(msg,92)'output mesh size:',isz1,jsz1
 509 call message(msg)
 510 92 format(a,2i8)
 511 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
 512 call crash('module_fr_spread_util:sum_mesh_cells: bad mesh sizes')
 513
 514 end subroutine sum_2d_cells
 515
 516
 517
 518 ! module_fr_fire_util%%interpolate_2d
 519 subroutine interpolate_2d(  &
 520     ims2,ime2,jms2,jme2, & ! array coarse grid
 521     its2,ite2,jts2,jte2, & ! dimensions coarse grid
 522     ims1,ime1,jms1,jme1, & ! array coarse grid
 523     its1,ite1,jts1,jte1, & ! dimensions fine grid
 524     ir,jr,               & ! refinement ratio
 525     rip2,rjp2,rip1,rjp1, & ! (rip2,rjp2) on grid 2 lines up with (rip1,rjp1) on grid 1
 526     v2, &                  ! in coarse grid
 527     v1  )                  ! out fine grid
 528 implicit none
 529
 530 !*** purpose
 531 ! interpolate nodal values in mesh2 to nodal values in mesh1
 532 ! interpolation runs over the mesh2 region its2:ite2,jts2:jte2
 533 ! only the part of mesh 1 in the region its1:ite1,jts1:jte1 is modified
 534
 535 !*** arguments
 536
 537 integer, intent(in)::its1,ite1,jts1,jte1,ims1,ime1,jms1,jme1
 538 integer, intent(in)::its2,ite2,jts2,jte2,ims2,ime2,jms2,jme2
 539 integer, intent(in)::ir,jr
 540 real,intent(in):: rjp1,rip1,rjp2,rip2
 541 real, intent(out)::v1(ims1:ime1,jms1:jme1)
 542 real, intent(in)::v2(ims2:ime2,jms2:jme2)
 543
 544 !*** local
 545 integer:: i1,i2,j1,j2,is,ie,js,je
 546 real:: tx,ty,rx,ry
 547 real:: rio,rjo
 548 intrinsic::ceiling,floor
 549
 550 !*** executable
 551
 552 !check mesh dimensions and domain dimensions
 553 call check_mesh_2dim(its1,ite1,jts1,jte1,ims1,ime1,jms1,jme1)
 554 call check_mesh_2dim(its2,ite2,jts2,jte2,ims2,ime2,jms2,jme2)
 555
 556 ! compute mesh ratios
 557 rx=1./ir
 558 ry=1./jr
 559
 560 do j2=jts2,jte2-1               ! loop over mesh 2 cells
 561     rjo=rjp1+jr*(j2-rjp2)       ! mesh 1 coordinate of the mesh 2 patch start
 562     js=max(jts1,ceiling(rjo))   ! lower bound of mesh 1 patch for this mesh 2 cell
 563     je=min(jte1,floor(rjo)+jr)  ! upper bound of mesh 1 patch for this mesh 2 cell
 564     do i2=its2,ite2-1
 565         rio=rip1+ir*(i2-rip2)
 566         is=max(its1,ceiling(rio))
 567         ie=min(ite1,floor(rio)+ir)
 568         do j1=js,je
 569             ty = (j1-rjo)*ry
 570             do i1=is,ie
 571                 ! in case mesh 1 node lies on the boundary of several mesh 2 cells
 572                 ! the result will be written multiple times with the same value
 573                 ! up to a rounding error
 574                 tx = (i1-rio)*rx
 575                 !print *,'coarse ',i2,j2,'to',i2+1,j2+1,' fine ',is,js,' to ',ie,je
 576                 v1(i1,j1)=                     &
 577                       (1-tx)*(1-ty)*v2(i2,j2)  &
 578                  +    (1-tx)*ty  *v2(i2,j2+1)  &
 579                  +      tx*(1-ty)*v2(i2+1,j2)  &
 580                  +        tx*ty  *v2(i2+1,j2+1)
 581                 !print *,'coarse ',i2,j2,' fine ',i1,j1, ' offset ',io,jo,' weights ',tx,ty, &
 582                 ! 'in ',v2(i2,j2),v2(i2,j2+1),v2(i2+1,j2),v2(i2+1,j2+1),' out ',v1(i1,j1)
 583                 !write(*,'(a,2i5,a,2f8.2,a,4f8.2,a,2i5,a,f8.2)') &
 584                 !'coarse ',i2,j2,' coord',rio,rjo,' val',v2(i2,j2),v2(i2,j2+1),v2(i2+1,j2),v2(i2+1,j2+1),&
 585                 !' fine ',i1,j1,' out ',v1(i1,j1)
 586            enddo
 587        enddo
 588     enddo
 589 enddo
 590
 591 end subroutine interpolate_2d
 592
 593 !
 594 !****************
 595 !
 596
 597 subroutine interpolate_2d_cells2cells(              &
 598       ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, &  ! in
 599       ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,v1  )  ! out
 600 implicit none
 601
 602 !*** purpose
 603 ! interpolate nodal values in mesh2 to nodal values in mesh1
 604 ! input mesh 2 is coarse output mesh 1 is fine
 605
 606 !*** arguments
 607
 608 integer, intent(in)::ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1
 609 real, intent(out)::v1(ims1:ime1,jms1:jme1)
 610 integer, intent(in)::ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2
 611 real, intent(in)::v2(ims2:ime2,jms2:jme2)
 612
 613 ! Example with mesh ratio=4. | = cell boundary,  x = cell center
 614 !
 615 !  mesh2   |-------x-------|-------x-------|
 616 !  mesh1   |-x-|-x-|-x-|-x-|-x-|-x-|-x-|-x-|
 617 !
 618
 619 !*** local
 620 integer:: ir,jr,isz1,isz2,jsz1,jsz2,ip,jp,ih,jh
 621 character(len=128)msg
 622
 623 !*** executable
 624
 625 !check mesh dimensions and domain dimensions
 626 call check_mesh_2dim(ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1)
 627 call check_mesh_2dim(ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2)
 628
 629 ! compute mesh sizes
 630 isz1 = ide1-ids1+1
 631 jsz1 = jde1-jds1+1
 632 isz2 = ide2-ids2+1
 633 jsz2 = jde2-jds2+1
 634
 635 ! check mesh sizes
 636 if(isz1.le.0.or.jsz1.le.0.or.isz2.le.0.or.jsz2.le.0)goto 9
 637 if(mod(isz1,isz2).ne.0.or.mod(jsz1,jsz2).ne.0)goto 9
 638
 639 ! compute mesh ratios
 640 ir=isz1/isz2
 641 jr=jsz1/jsz2
 642 !
 643 !  mesh2   |-------x-------|-------x-------|
 644 !  mesh1   |-x-|-x-|-x-|-x-|-x-|-x-|-x-|-x-|
 645
 646 !  mesh2   |-----x-----|-----x-----|  rx=3
 647 !  mesh1   |-x-|-x-|-x-|-x-|-x-|-x-|
 648 !  i2            1   1   1   2
 649 !  i1        1   2   3   4   5
 650 !  ioff          0   1   2   0
 651 !  tx            0  1/3 2/3
 652
 653 !  mesh2   |---x---|---x---| rx=2
 654 !  mesh1   |-x-|-x-|-x-|-x-|
 655 !  i2            1   1   2
 656 !  i1            2   3   4
 657 !  ioff          0   1   2
 658 !  tx           1/4 3/4
 659
 660
 661 ! offset of the last node in the 1st half of the cell
 662 ih=ir/2
 663 jh=jr/2
 664 ! 0 if coarse cell center coincides with fine, 1 if not
 665 ip=mod(ir+1,2)
 666 jp=mod(jr+1,2)
 667
 668 call interpolate_2d_w(ip,jp,ih,jh,ir,jr,              &
 669       ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, &  ! in
 670       ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,v1  )  ! out
 671
 672 return
 673
 674 9 continue
 675 !$OMP CRITICAL(FIRE_UTIL_CRIT)
 676 write(msg,91)ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2
 677 call message(msg)
 678 write(msg,91)ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1
 679 call message(msg)
 680 write(msg,92)'input  mesh size:',isz2,jsz2
 681 call message(msg)
 682 91 format('dimensions: ',8i8)
 683 write(msg,92)'output mesh size:',isz1,jsz1
 684 call message(msg)
 685 92 format(a,2i8)
 686 call crash("module_fr_fire_util:interpolate_2dmesh_cells: bad mesh sizes")
 687 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
 688 end subroutine interpolate_2d_cells2cells
 689
 690 !
 691 !****************
 692 !
 693
 694 subroutine interpolate_2d_cells2nodes(              &
 695       ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, &  ! in
 696       ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,v1  )  ! out
 697 implicit none
 698
 699 !*** purpose
 700 ! interpolate nodal values in mesh2 to nodal values in mesh1
 701 ! input mesh 2 is coarse output mesh 1 is fine
 702
 703 !*** arguments
 704
 705 integer, intent(in)::ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1
 706 real, intent(out)::v1(ims1:ime1,jms1:jme1)
 707 integer, intent(in)::ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2
 708 real, intent(in)::v2(ims2:ime2,jms2:jme2)
 709
 710 ! Example with mesh ratio=4. | = cell boundary,  x = cell center
 711 !
 712 !  mesh2   |-------x-------|-------x-------|
 713 !  mesh1   x-|-x-|-x-|-x-|-x-|-x-|-x-|-x-|-x
 714 !
 715
 716 !*** local
 717 integer:: ir,jr,isz1,isz2,jsz1,jsz2,ip,jp,ih,jh
 718 character(len=128)msg
 719
 720 !*** executable
 721
 722 !check mesh dimensions and domain dimensions
 723 call check_mesh_2dim(ids1,ide1+1,jds1,jde1+1,ims1,ime1,jms1,jme1)
 724 call check_mesh_2dim(ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2)
 725
 726 ! compute mesh sizes
 727 isz1 = ide1-ids1+1
 728 jsz1 = jde1-jds1+1
 729 isz2 = ide2-ids2+1
 730 jsz2 = jde2-jds2+1
 731
 732 ! check mesh sizes
 733 if(isz1.le.0.or.jsz1.le.0.or.isz2.le.0.or.jsz2.le.0)goto 9
 734 if(mod(isz1,isz2).ne.0.or.mod(jsz1,jsz2).ne.0)goto 9
 735
 736 ! compute mesh ratios
 737 ir=isz1/isz2
 738 jr=jsz1/jsz2
 739 !
 740 !  mesh2   |-------x-------|-------x-------|
 741 !  mesh1   x-|-x-|-x-|-x-|-x-|-x-|-x-|-x-|-x
 742
 743 !  mesh2   |-----x-----|-----x-----|  rx=3
 744 !  mesh1   x-|-x-|-x-|-x-|-x-|-x-|-x
 745
 746 !  mesh2   |---x---|---x---| rx=2
 747 !  mesh1   x-|-x-|-x-|-x-|-x
 748
 749 ! offset of the last node in the 1st half of the cell
 750 ih=(ir+1)/2
 751 jh=(jr+1)/2
 752 ! 0 if coarse cell center coincides with fine, 1 if not
 753 ip=mod(ir,2)
 754 jp=mod(jr,2)
 755
 756
 757 call interpolate_2d_w(ip,jp,ih,jh,ir,jr,              &
 758       ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, &  ! in
 759       ids1,ide1+1,jds1,jde1+1,ims1,ime1,jms1,jme1,v1  )  ! out
 760
 761
 762 return
 763 9 continue
 764 !$OMP CRITICAL(FIRE_UTIL_CRIT)
 765 write(msg,91)ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2
 766 call message(msg)
 767 write(msg,91)ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1
 768 call message(msg)
 769 write(msg,92)'input  mesh size:',isz2,jsz2
 770 call message(msg)
 771 91 format('dimensions: ',8i8)
 772 write(msg,92)'output mesh size:',isz1,jsz1
 773 call message(msg)
 774 92 format(a,2i8)
 775 call crash("module_fr_fire_util:interpolate_2d_cells2nodes: bad mesh sizes")
 776 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
 777 end subroutine interpolate_2d_cells2nodes
 778 !
 779 !****************
 780 !
 781
 782 subroutine interpolate_2d_w(ip,jp,ih,jh,ir,jr,             &
 783       ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, &  ! in
 784       ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,v1  )  ! out
 785 implicit none
 786 !*** EXCEPTION: THIS SUBROUTINE IS NEITHER CELL NOR NODE BASED.
 787
 788 integer, intent(in)::ip,jp,ih,jh,ir,jr
 789 integer, intent(in)::ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1
 790 real, intent(out)::v1(ims1:ime1,jms1:jme1)
 791 integer, intent(in)::ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2
 792 real, intent(in)::v2(ims2:ime2,jms2:jme2)
 793
 794 real:: tx,ty,rx,ry,half,xoff,yoff
 795 integer:: i1,i2,j1,j2,ioff,joff
 796 parameter(half=0.5)
 797
 798 rx=ir
 799 ry=jr
 800
 801 xoff = ip*half
 802 yoff = jp*half
 803
 804 ! the inside, ids1+ih:ide1-ih,jds1+jh:jde1-jh
 805 do j2=jds2,jde2-1     ! interpolate from nodes j2 and j2+1
 806     do i2=ids2,ide2-1
 807         do ioff=0,ir-ip
 808             do joff=0,jr-jp
 809                 ! compute fine mesh index
 810                 i1=ioff+(ih+ids1)+ir*(i2-ids2)
 811                 j1=joff+(jh+jds1)+jr*(j2-jds2)
 812                 ! weights
 813                 tx = (ioff+xoff)/rx
 814                 ty = (joff+yoff)/ry
 815                 ! interpolation
 816                 v1(i1,j1)=                     &
 817                       (1-tx)*(1-ty)*v2(i2,j2)  &
 818                  +    (1-tx)*ty  *v2(i2,j2+1)  &
 819                  +      tx*(1-ty)*v2(i2+1,j2)  &
 820                  +        tx*ty  *v2(i2+1,j2+1)
 821                 !write(*,'(3(a,2i5),a,2f7.4)')'coarse ',i2,j2,' fine ',i1,j1, &
 822                 ! ' offset ',ioff,joff,' weights ',tx,ty
 823                 !write(*,'(a,4f7.4,a,f7.4)')'in ',v2(i2,j2),v2(i2,j2+1),v2(i2+1,j2), &
 824                 !  v2(i2+1,j2+1),' out ',v1(i1,j1)
 825            enddo
 826        enddo
 827     enddo
 828 enddo
 829
 830 ! extend to the boundary strips from the nearest known
 831 do ioff=0,ih-1  ! top and bottom strips
 832     do j2=jds2,jde2-1
 833         do joff=0,jr-jp
 834            j1=joff+(jh+jds1)+jr*(j2-jds2)
 835            ! weights
 836            ty = (joff+yoff)/ry
 837            ! interpolation
 838            v1(ids1+ioff,j1)=(1-ty)*v2(ids2,j2)+ty*v2(ids2,j2+1)
 839            v1(ide1-ioff,j1)=(1-ty)*v2(ide2,j2)+ty*v2(ide2,j2+1)
 840        enddo
 841     enddo
 842 enddo
 843 do joff=0,jh-1  ! left and right strips
 844     do i2=ids2,ide2-1
 845         do ioff=0,ir-ip
 846            i1=ioff+(ih+ids1)+ir*(i2-ids2)
 847            ! weights
 848            tx = (ioff+xoff)/rx
 849            ! interpolation
 850            v1(i1,jds1+joff)=(1-tx)*v2(i2,jds2)+tx*v2(i2+1,jds2)
 851            v1(i1,jde1-joff)=(1-tx)*v2(i2,jde2)+tx*v2(i2+1,jde2)
 852        enddo
 853     enddo
 854 enddo
 855 ! extend to the 4 corner squares from the nearest known
 856 do ioff=0,ih-1
 857     do joff=0,jh-1
 858         v1(ids1+ioff,jds1+joff)=v2(ids2,jds2)
 859         v1(ide1-ioff,jds1+joff)=v2(ide2,jds2)
 860         v1(ids1+ioff,jde1-joff)=v2(ids2,jde2)
 861         v1(ide1-ioff,jde1-joff)=v2(ide2,jde2)
 862     enddo
 863 enddo
 864 end subroutine interpolate_2d_w
 865
 866 !
 867 !****************
 868 !
 869
 870 real function interp(ids,ide,jds,jde,ims,ime,jms,jme,x,y,v)
 871 implicit none
 872 !*** purpose
 873 ! general interpolation in a rectangular
 874
 875 !*** arguments
 876
 877 integer, intent(in)::ids,ide,jds,jde,ims,ime,jms,jme
 878 real, intent(in)::x,y,v(ims:ime,jms:jme)
 879 ! the mesh is cell based so the used dimension of v is ids:ide+1,jds:jde+1
 880
 881 !*** calls
 882 intrinsic floor,min,max
 883
 884 !*** local
 885 integer i,j
 886 real tx,ty
 887
 888 ! executable
 889
 890 ! indices of the lower left corner of the cell in the mesh that contains (x,y)
 891 i = floor(x)
 892 i=max(min(i,ide),ids)
 893 j = floor(y)
 894 j=max(min(j,jde),jds)
 895
 896 ! the leftover
 897 tx = x - real(i)
 898 ty = y - real(j)
 899
 900 ! interpolate the values
 901 interp = &
 902                     (1-tx)*(1-ty)*v(i,j)    &
 903                  +    tx*(1-ty)  *v(i+1,j)  &
 904                  +    (1-tx)*ty  *v(i,j+1)  &
 905                  +        tx*ty  *v(i+1,j+1)
 906
 907 !print *,'x,y=',x,y,'i1,i2=',i1,j1,'tx,ty=',tx,ty,' interp=',interp
 908 end function interp
 909
 910 subroutine meshdiffc_2d(ids, ide, jds,jde , &    ! mesh area used (in cells, end +1)
 911                    ims1,ime1,jms1,jme1, &       ! memory dimensiuons
 912                    dx,dy,               &       ! mesh spacing
 913                    lfn,                 &       ! input
 914                    diffCx,diffCy) ! output
 915 implicit none
 916
 917 !*** purpose
 918 ! central differences on a 2d mesh
 919
 920 !*** arguments
 921
 922 integer, intent(in)::ids,ide,jds,jde,ims1,ime1,jms1,jme1
 923 real, intent(in):: dx,dy
 924 real, intent(in), dimension(ims1:ime1,jms1:jme1):: lfn
 925 real, intent(out), dimension(ims1:ime1,jms1:jme1):: diffCx,diffCy
 926
 927 !*** local
 928 integer:: i,j
 929 real, dimension(ims1:ime1,jms1:jme1):: diffLx,diffRx,diffLy,diffRy
 930
 931 ! get one-sided differences; dumb but had that already...
 932 call meshdiff_2d(ids, ide, jds,jde , &    ! mesh area used (in cells, end +1)
 933                    ims1,ime1,jms1,jme1, &       ! dimensions of lfn
 934                    dx,dy,               &       ! mesh spacing
 935                    lfn,                 &       ! input
 936                    diffLx,diffRx,diffLy,diffRy) ! output
 937
 938 ! make into central
 939 do j=jds,jde+1
 940     do i=ids,ide+1
 941         diffCx(i,j)=0.5*(diffLx(i,j) + diffRx(i,j))
 942         diffCy(i,j)=0.5*(diffLy(i,j) + diffRy(i,j))
 943     enddo
 944 enddo
 945 end subroutine meshdiffc_2d
 946
 947 subroutine meshdiff_2d(ids, ide, jds,jde , &    ! mesh area used (in cells, end +1)
 948                    ims1,ime1,jms1,jme1, &       ! dimensions of lfn
 949                    dx,dy,               &       ! mesh spacing
 950                    lfn,                 &       ! input
 951                    diffLx,diffRx,diffLy,diffRy) ! output
 952 implicit none
 953
 954 !*** purpose
 955 ! one-sided differences on a 2d mesh
 956
 957 !*** arguments
 958
 959 integer, intent(in)::ids,ide,jds,jde,ims1,ime1,jms1,jme1
 960 real, intent(in):: dx,dy
 961 real, intent(in), dimension(ims1:ime1,jms1:jme1):: lfn
 962 real, intent(out), dimension(ims1:ime1,jms1:jme1):: diffLx,diffRx,diffLy,diffRy
 963
 964 !*** local
 965 integer:: i,j
 966 real:: tmpx,tmpy
 967
 968 !*** executable
 969
 970     call check_mesh_2dim(ids,ide+1,jds,jde+1,ims1,ime1,jms1,jme1)
 971
 972     ! the bulk of the work
 973     do j=jds,jde
 974         do i=ids,ide
 975             tmpx = (lfn(i+1,j)-lfn(i,j))/dx
 976             diffLx(i+1,j) = tmpx
 977             diffRx(i,j)   = tmpx
 978             tmpy = (lfn(i,j+1)-lfn(i,j))/dy
 979             diffLy(i,j+1) = tmpy
 980             diffRy(i,j)   = tmpy
 981         enddo
 982         ! missing values - put there the other one
 983         diffLx(ids,j)  = diffLx(ids+1,j)
 984         diffRx(ide+1,j)= diffRx(ide,j)
 985     enddo
 986     ! cleanup
 987     ! j=jde+1 from above loop
 988     do i=ids,ide
 989         tmpx = (lfn(i+1,j)-lfn(i,j))/dx
 990         diffLx(i+1,j) = tmpx
 991         diffRx(i,j)   = tmpx
 992     enddo
 993     ! i=ide+1 from above loop
 994     do j=jds,jde
 995         tmpy = (lfn(i,j+1)-lfn(i,j))/dy
 996         diffLy(i,j+1) = tmpy
 997         diffRy(i,j)   = tmpy
 998     enddo
 999     ! missing values - put there the other one
1000     ! j=jde+1 from above loop, j=jds:jde done before in main bulk loop
1001     diffLx(ids,j)   = diffLx(ids+1,j)
1002     diffRx(ide+1,j) = diffRx(ide,j)
1003     do i=ids,ide+1
1004         diffLy(i,jds)   = diffLy(i,jds+1)
1005         diffRy(i,jde+1) = diffRy(i,jde)
1006     enddo
1007
1008 end subroutine meshdiff_2d
1009
1010
1011
1012
1013 real pure function sum_2darray( its,ite,jts,jte,               &
1014                                 ims,ime,jms,jme,               &
1015                                 a)
1016 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme
1017 real, intent(in)::a(ims:ime,jms:jme)
1018 !*** local
1019 integer:: i,j
1020 real:: t
1021 t=0.
1022 do j=jts,jte
1023     do i=its,ite
1024         t=t+a(i,j)
1025     enddo
1026 enddo
1027 sum_2darray = t
1028 end function sum_2darray
1029
1030 real pure function max_2darray( its,ite,jts,jte,               &
1031                                 ims,ime,jms,jme,               &
1032                                 a)
1033 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme
1034 real, intent(in)::a(ims:ime,jms:jme)
1035 !*** local
1036 integer:: i,j
1037 real:: t
1038 t=0.
1039 do j=jts,jte
1040     do i=its,ite
1041         t=max(t,a(i,j))
1042     enddo
1043 enddo
1044 max_2darray = t
1045 end function max_2darray
1046
1047 subroutine print_2d_stats_vec(ips,ipe,jps,jpe, &
1048                          ims,ime,jms,jme, &
1049                          ax,ay,name)
1050 implicit none
1051 integer, intent(in)::ips,ipe,jps,jpe,ims,ime,jms,jme
1052 real, intent(in), dimension(ims:ime,jms:jme)::ax,ay
1053 character(len=*),intent(in)::name
1054 integer:: i,j
1055 real:: t
1056 real:: avg_a,max_a,min_a
1057 character(len=25)::id
1058 id=name
1059 call print_2d_stats(ips,ipe,jps,jpe, &
1060                          ims,ime,jms,jme, &
1061                          ax,id//'/x ')
1062 call print_2d_stats(ips,ipe,jps,jpe, &
1063                          ims,ime,jms,jme, &
1064                          ay,id//'/y ')
1065 avg_a=0
1066 max_a=-huge(max_a)
1067 min_a= huge(min_a)
1068 do j=jps,jpe
1069     do i=ips,ipe
1070         t=sqrt(ax(i,j)**2+ay(i,j)**2)
1071         max_a=max(max_a,t)
1072         min_a=min(min_a,t)
1073         avg_a=avg_a+t
1074     enddo
1075 enddo
1076 avg_a = avg_a/((ipe-ips+1)*(jpe-jps+1))
1077 call print_stat_line(id//'/sz',ips,ipe,jps,jpe,min_a,max_a,avg_a)
1078 end subroutine print_2d_stats_vec
1079
1080
1081 subroutine print_stat_line(name,ips,ipe,jps,jpe,min_a,max_a,avg_a)
1082 !*** encapsulate line with statistics
1083 implicit none
1084 !*** arguments
1085 integer, intent(in)::ips,ipe,jps,jpe
1086 character(len=*),intent(in)::name
1087 real,intent(in)::min_a,max_a,avg_a
1088 !*** local
1089 character(len=128)::msg
1090 character(len=24)::id
1091 if(fire_print_msg.eq.0)return
1092 id=name
1093 !$OMP CRITICAL(FIRE_UTIL_CRIT)
1094 write(msg,'(a,4i4,3g11.3)')id,ips,ipe,jps,jpe,min_a,max_a,avg_a
1095 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
1096 call message(msg)
1097 if(.not.avg_a.eq.avg_a)call crash('NaN detected')
1098 end subroutine print_stat_line
1099
1100
1101 subroutine print_3d_stats_by_slice(ips,ipe,kps,kpe,jps,jpe, &
1102                          ims,ime,kms,kme,jms,jme, &
1103                          a,name)
1104 implicit none
1105 integer, intent(in)::ips,ipe,jps,jpe,ims,ime,jms,jme,kms,kme,kps,kpe
1106 real, intent(in)::a(ims:ime,kms:kme,jms:jme)
1107 character(len=*),intent(in)::name
1108 integer::k
1109 character(len=128)::msg
1110 do k=kps,kpe
1111 ! print 3d stats for each horizontal slice separately
1112 !$OMP CRITICAL(SFIRE_UTIL_CRIT)
1113     write(msg,'(i2,1x,a)')k,name
1114 !$OMP END CRITICAL(SFIRE_UTIL_CRIT)
1115     call print_3d_stats(ips,ipe,k,k,jps,jpe, &
1116                          ims,ime,kms,kme,jms,jme, &
1117                          a,msg)
1118 enddo
1119 end subroutine print_3d_stats_by_slice
1120
1121
1122 subroutine print_3d_stats(ips,ipe,kps,kpe,jps,jpe, &
1123                          ims,ime,kms,kme,jms,jme, &
1124                          a,name)
1125 implicit none
1126 integer, intent(in)::ips,ipe,jps,jpe,ims,ime,jms,jme,kms,kme,kps,kpe
1127 real, intent(in)::a(ims:ime,kms:kme,jms:jme)
1128 character(len=*),intent(in)::name
1129 integer:: i,j,k
1130 real:: avg_a,max_a,min_a,t,aa,bb
1131 character(len=128)::msg
1132 ! if(fire_print_msg.eq.0)return
1133 ! check for nans in any case
1134 bb=0.
1135 do j=jps,jpe
1136   do k=kps,kpe
1137     do i=ips,ipe
1138        bb=bb+a(i,k,j)
1139     enddo
1140   enddo
1141 enddo
1142 if(bb.eq.bb.and.fire_print_msg.eq.0)return
1143 avg_a=0
1144 max_a=-huge(max_a)
1145 min_a= huge(min_a)
1146 t=huge(t)
1147 do j=jps,jpe
1148   do k=kps,kpe
1149     do i=ips,ipe
1150         aa=a(i,k,j)
1151         if(aa.ne.aa.or..not.aa.le.t.or..not.aa.ge.-t)goto 9
1152         max_a=max(max_a,aa)
1153         min_a=min(min_a,aa)
1154         avg_a=avg_a+aa
1155     enddo
1156   enddo
1157 enddo
1158 if(bb.ne.bb)goto 10 ! should never happen
1159 if(fire_print_msg.eq.0)return
1160 avg_a = avg_a/((ipe-ips+1)*(jpe-jps+1)*(kpe-kps+1))
1161 call print_stat_line(name,ips,ipe,jps,jpe,min_a,max_a,avg_a)
1162 return
1163 9 continue
1164 !$OMP CRITICAL(FIRE_UTIL_CRIT)
1165 write(msg,1)name,i,k,j,aa
1166 call message(msg)
1167 1 format(a30,'(',i6,',',i6,',',i6,') = ',g13.5)
1168 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
1169 call print_stat_line(name,ips,ipe,jps,jpe,aa,aa,aa)
1170 if(aa.ne.aa)goto 10
1171 msg='Invalid floating point number detected in '//name
1172 call crash(msg)
1173 10 msg='NaN detected in '//name
1174 call crash(msg)
1175 end subroutine print_3d_stats
1176
1177 subroutine print_2d_stats(ips,ipe,jps,jpe, &
1178                          ims,ime,jms,jme, &
1179                          a,name)
1180 implicit none
1181 integer, intent(in)::ips,ipe,jps,jpe,ims,ime,jms,jme
1182 real, intent(in)::a(ims:ime,jms:jme)
1183 character(len=*),intent(in)::name
1184 !!character(len=128)::msg
1185 !if(fire_print_msg.eq.0)return
1186 call print_3d_stats(ips,ipe,1,1,jps,jpe, &
1187                          ims,ime,1,1,jms,jme, &
1188                          a,name)
1189 !!write(msg,'(2a,z16)')name,' address =',loc(a)
1190 !!call message(msg)
1191 end subroutine print_2d_stats
1192
1193 real pure function avg_2darray( its,ite,jts,jte,               &
1194                                 ims,ime,jms,jme,               &
1195                                 a)
1196 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme
1197 real, intent(in)::a(ims:ime,jms:jme)
1198 !*** local
1199 !*** executable
1200 avg_2darray = sum_2darray( its,ite,jts,jte,               &
1201                            ims,ime,jms,jme,               &
1202                            a)/((ite-its+1)*(jte-jts+1))
1203 end function avg_2darray
1204
1205 real pure function avg_2darray_vec( its,ite,jts,jte,           &
1206                                 ims,ime,jms,jme,               &
1207                                 ax,ay)
1208 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme
1209 real, intent(in), dimension(ims:ime,jms:jme):: ax,ay
1210 !*** local
1211 integer:: i,j
1212 real:: t
1213 t=0.
1214 do j=jts,jte
1215     do i=its,ite
1216         t=t+sqrt(ax(i,j)**2+ay(i,j)**2)
1217     enddo
1218 enddo
1219 t = t/((ite-its+1)*(jte-jts+1))
1220 avg_2darray_vec = t
1221 end function avg_2darray_vec
1222
1223
1224 subroutine print_array(its,ite,jts,jte,           &
1225                          ims,ime,jms,jme,               &
1226                          a,name,id)
1227 ! debug
1228 !*** arguments
1229 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme,id
1230 real, intent(in), dimension(ims:ime,jms:jme):: a
1231 character(len=*),intent(in)::name
1232 !****
1233 integer i,j
1234 character(len=128)::msg
1235 !****
1236 !$OMP CRITICAL(FIRE_UTIL_CRIT)
1237 write(msg,*)name,' start ',id,' dim ',its,ite,jts,jte
1238 call message(msg)
1239 do j=jts,jte
1240     do i=its,ite
1241          write(msg,*)i,j,a(i,j)
1242          call message(msg)
1243     enddo
1244 enddo
1245 write(msg,*)name,' end ',id
1246 call message(msg)
1247 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
1248 end subroutine print_array
1249
1250 subroutine write_array_m(its,ite,jts,jte,           &
1251                          ims,ime,jms,jme,               &
1252                          a,name,id)
1253 ! debug
1254 !*** arguments
1255 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme,id
1256 real, intent(in), dimension(ims:ime,jms:jme):: a
1257 character(len=*),intent(in)::name
1258 !****
1259 call write_array_m3(its,ite,1,1,jts,jte,           &
1260                          ims,ime,1,1,jms,jme,               &
1261                          a,name,id)
1262 end subroutine write_array_m
1263
1264
1265 subroutine write_array_m3(its,ite,kts,kte,jts,jte,           &
1266                          ims,ime,kms,kme,jms,jme,               &
1267                          a,name,id)
1268
1269 implicit none
1270 ! debug
1271 !*** arguments
1272 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme,kts,kte,kms,kme,id
1273 real, intent(in), dimension(ims:ime,kms:kme,jms:jme):: a
1274 character(len=*),intent(in)::name
1275 !****
1276 integer i,j,k,iu,ilen,myproc,nprocs
1277 logical op
1278 character(len=128)::fname,msg
1279 !****
1280 if(fire_print_file.eq.0.or.id.le.0)return
1281 call check_mesh_2dim(its,ite,jts,jte,ims,ime,jms,jme)
1282 call wrf_get_nproc (nprocs)
1283 call wrf_get_myproc(myproc)
1284
1285 !$OMP CRITICAL(FIRE_UTIL_CRIT)
1286 if(nprocs.eq.1)then
1287     write(fname,3)name,'_',id,'.txt'
1288 else
1289     write(fname,4)name,'_',id,'.',myproc,'.txt'
1290 endif
1291
1292 iu=0
1293 do i=6,99
1294     inquire(unit=i,opened=op)
1295     if(.not.op.and.iu.le.0)iu=i
1296 enddo
1297 if(iu.gt.0)open(iu,file=trim(fname),form='formatted',status='unknown')
1298
1299 if(iu.le.0)call crash('write_array_m: cannot find available fortran unit')
1300
1301 write(iu,1)real(its)
1302 write(iu,1)real(ite)
1303 write(iu,1)real(jts)
1304 write(iu,1)real(jte)
1305 write(iu,1)real(kts)
1306 write(iu,1)real(kte)
1307 write(iu,1)(((a(i,k,j),i=its,ite),j=jts,jte),k=kts,kte)
1308 close(iu)
1309 write(msg,2)name,'(',its,':',ite,',',jts,':',jte,',', &
1310 kts,':',kte,') -> ',trim(fname)
1311 !$OMP END CRITICAL(FIRE_UTIL_CRIT)
1312 call message(msg)
1313 return
1314
1315 1 format(e20.12)
1316 2 format(2a,3(i5,a,i5,a),2a)
1317 3 format(a,a,i8.8,a)
1318 4 format(a,a,i8.8,a,i4.4,a)
1319
1320
1321 end subroutine write_array_m3
1322 !
1323 !***
1324 !
1325 subroutine read_array_2d_integer(filename,ia,its,ite,jts,jte,ims,ime,jms,jme)
1326 !*** arguments
1327 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme
1328 integer, intent(out), dimension(ims:ime,jms:jme):: ia
1329 character(len=*),intent(in)::filename
1330 !*** local
1331 real, allocatable, dimension(:,:)::a
1332 integer::i,j
1333 real::r
1334 character(len=256)msg
1335 !*** executable
1336 allocate(a(ims:ime,jms:jme))
1337 call read_array_2d_real(filename,a,its,ite,jts,jte,ims,ime,jms,jme)
1338 do j=jts,jte
1339     do i=its,ite
1340         ia(i,j)=a(i,j)
1341         r=ia(i,j)
1342         if(r.ne.a(i,j))then
1343            write(6,*)'Reading file ',trim(filename)
1344            write(6,*)'value at position ',i,j,' cannot be converted to integer'
1345            write(6,*)'read ',a(i,j),' converted as',ia(i,j),' converted as ',r
1346            msg=trim(filename)//' is not integer file'
1347            call crash(msg)
1348         endif
1349     enddo
1350 enddo
1351 end subroutine read_array_2d_integer
1352
1353
1354
1355 subroutine read_array_2d_real(filename,a,its,ite,jts,jte,ims,ime,jms,jme)
1356 #ifdef _OPENMP
1357 use OMP_LIB
1358 #endif
1359 implicit none
1360 !*** purpose: read a 2D array from a text file
1361 ! file format: line 1: array dimensions ni,nj
1362 !              following lines: one row of a at a time
1363 !              each row may be split between several lines
1364 !*** arguments
1365 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme
1366 real, intent(out), dimension(ims:ime,jms:jme):: a
1367 character(len=*),intent(in)::filename
1368 !*** local
1369 integer i,j,ni,nj,mi,mj,nprocs,myproc,mythread,iu
1370 logical op
1371 character(len=128)::fname,msg
1372 !*** executable
1373
1374 call wrf_get_nproc (nprocs)
1375 call wrf_get_myproc( myproc )
1376 mythread=0
1377 #ifdef _OPENMP
1378     mythread=omp_get_thread_num()
1379 #endif
1380 if(nprocs.ne.1.or.myproc.ne.0.or.mythread.ne.0) &
1381    call crash('read_array_2d: parallel execution not supported')
1382
1383 ! print line
1384 mi=ite-its+1
1385 mj=jte-jts+1
1386 write(msg,2)'reading array size ',mi,mj,' from file ',trim(filename)
1387 2 format(a,2i6,2a)
1388 call message(msg)
1389
1390 ! check array index overflow
1391 call check_mesh_2dim(its,ite,jts,jte,ims,ime,jms,jme)
1392
1393 ! find available unit
1394 iu=0
1395 do i=11,99
1396     inquire(unit=i,opened=op)
1397     if(.not.op.and.iu.le.0)iu=i
1398 enddo
1399 if(iu.le.0)call crash('read_array_2d: cannot find available fortran unit')
1400
1401 if(iu.gt.0)open(iu,file=filename,form='formatted',status='old',err=9)
1402 rewind(iu,err=9)
1403
1404 read(iu,*,err=10)ni,nj
1405 if(ni.ne.mi.or.nj.ne.mj)then
1406     write(msg,'(a,2i6,a,2i6)')'Array dimensions',ni,nj,' in the input file should be ',mi,mj
1407     call message(msg)
1408     goto 10
1409 endif
1410 do i=its,ite
1411    read(iu,*,err=10)(a(i,j),j=jts,jte)
1412 enddo
1413 close(iu,err=11)
1414 return
1415
1416 9  msg='Error opening file '//trim(filename)
1417 call crash(msg)
1418 10 msg='Error reading file '//trim(filename)
1419 call crash(msg)
1420 11 msg='Error closing file '//trim(filename)
1421 call crash(msg)
1422 end subroutine read_array_2d_real
1423 !
1424 !***
1425 !
1426
1427 ! general conditional expression
1428 pure integer function ifval(l,i,j)
1429 implicit none
1430 logical, intent(in)::l
1431 integer, intent(in)::i,j
1432 if(l)then
1433         ifval=i
1434 else
1435         ifval=j
1436 endif
1437 end function ifval
1438
1439 ! function to go beyond domain boundary if tile is next to it
1440 pure integer function snode(t,d,i)
1441 implicit none
1442 integer, intent(in)::t,d,i
1443 if(t.ne.d)then
1444     snode=t
1445 else
1446     snode=t+i
1447 endif
1448 end function snode
1449
1450 subroutine print_chsum( id, &
1451     ims,ime,kms,kme,jms,jme, &                ! memory dims
1452     ids,ide,kds,kde,jds,jde, &                ! domain dims
1453     ips,ipe,kps,kpe,jps,jpe, &                ! patch or tile dims
1454     istag,kstag,jstag,       &
1455     a,name)
1456
1457 #ifdef DM_PARALLEL
1458     USE module_dm , only : wrf_dm_bxor_integer
1459 #endif
1460
1461 integer, intent(in):: id, &
1462     ims,ime,kms,kme,jms,jme, &                ! memory dims
1463     ids,ide,kds,kde,jds,jde, &                ! domain dims
1464     ips,ipe,kps,kpe,jps,jpe, &                ! patch dims
1465     istag,kstag,jstag
1466 real, intent(in),dimension(ims:ime,kms:kme,jms:jme)::a
1467 character(len=*)::name
1468
1469 #ifdef _OPENMP
1470 !external, logical:: omp_in_parallel
1471 !external, integer:: omp_get_thread_num
1472 #endif
1473
1474 !*** local
1475 integer::lsum
1476 integer::i,j,k,n,ipe1,jpe1,kpe1,iel,thread,is,js,ks
1477 integer, save::psum,gsum
1478 real::rel
1479 equivalence(rel,iel)
1480 character(len=256)msg
1481
1482 if(fire_print_msg.le.0)return
1483
1484 ipe1=ifval(ipe.eq.ide.and.istag.ne.0,ipe+1,ipe)
1485 kpe1=ifval(kpe.eq.kde.and.kstag.ne.0,kpe+1,kpe)
1486 jpe1=ifval(jpe.eq.jde.and.jstag.ne.0,jpe+1,jpe)
1487 is=ifval(istag.ne.0,1,0)
1488 ks=ifval(kstag.ne.0,1,0)
1489 js=ifval(jstag.ne.0,1,0)
1490
1491 lsum=0
1492 do j=jps,jpe1
1493   do k=kps,kpe1-1
1494     do i=ips,ipe1
1495       rel=a(i,k,j)
1496       lsum=ieor(lsum,iel)
1497     enddo
1498   enddo
1499 enddo
1500
1501 ! get process sum over all threads
1502 thread=0
1503 #ifdef _OPENMP
1504 !thread=omp_get_thread_num()
1505 #endif
1506 if(thread.eq.0)psum=0
1507 !$OMP BARRIER
1508 !$OMP CRITICAL(CHSUM)
1509 psum=ieor(psum,lsum)
1510 !$OMP END CRITICAL(CHSUM)
1511 !$OMP BARRIER
1512
1513 ! get global sum over all processes
1514 if(thread.eq.0)then
1515 #ifdef DM_PARALLEL
1516     gsum = wrf_dm_bxor_integer ( psum )
1517 #else
1518     gsum = psum
1519 #endif
1520     write(msg,1)id,name,ids,ide+is,kds,kde+ks,jds,jde+js,gsum
1521 1   format(i6,1x,a10,' dims',6i5,' chsum ',z8.8)
1522     call message(msg)
1523 endif
1524
1525 end subroutine print_chsum
1526
1527
1528 real function fun_real(fun,  &
1529     ims,ime,kms,kme,jms,jme, &                ! memory dims
1530     ids,ide,kds,kde,jds,jde, &                ! domain dims
1531     ips,ipe,kps,kpe,jps,jpe, &                ! patch or tile dims
1532     istag,kstag,jstag,       &                ! staggering
1533     a,b)
1534
1535 #ifdef DM_PARALLEL
1536     USE module_dm , only : wrf_dm_sum_real , wrf_dm_max_real
1537 #endif
1538
1539 integer, intent(in)::  fun, &
1540     ims,ime,kms,kme,jms,jme, &                ! memory dims
1541     ids,ide,kds,kde,jds,jde, &                ! domain dims
1542     ips,ipe,kps,kpe,jps,jpe, &                ! patch dims
1543     istag,kstag,jstag                         ! staggering
1544 real, intent(in),dimension(ims:ime,kms:kme,jms:jme)::a,b
1545
1546 !*** local
1547 real::lsum,void
1548 integer::i,j,k,n,ipe1,jpe1,kpe1,iel,thread,is,js,ks
1549 real, save::psum,gsum
1550 real::rel
1551 logical:: dosum,domax
1552 character(len=256)msg
1553
1554 ipe1=ifval(ipe.eq.ide.and.istag.ne.0,ipe+1,ipe)
1555 kpe1=ifval(kpe.eq.kde.and.kstag.ne.0,kpe+1,kpe)
1556 jpe1=ifval(jpe.eq.jde.and.jstag.ne.0,jpe+1,jpe)
1557 is=ifval(istag.ne.0,1,0)
1558 ks=ifval(kstag.ne.0,1,0)
1559 js=ifval(jstag.ne.0,1,0)
1560
1561 if(fun.eq.REAL_SUM)then
1562   void=0.
1563   lsum=void
1564   do j=jps,jpe1
1565     do k=kps,kpe1
1566       do i=ips,ipe1
1567         lsum=lsum+a(i,k,j)
1568       enddo
1569     enddo
1570   enddo
1571 elseif(fun.eq.RNRM_SUM)then
1572   void=0.
1573   lsum=void
1574   do j=jps,jpe1
1575     do k=kps,kpe1
1576       do i=ips,ipe1
1577         lsum=lsum+sqrt(a(i,k,j)*a(i,k,j)+b(i,k,j)*b(i,k,j))
1578       enddo
1579     enddo
1580   enddo
1581 elseif(fun.eq.REAL_MAX)then
1582   void=-huge(lsum)
1583   lsum=void
1584   do j=jps,jpe1
1585     do k=kps,kpe1
1586       do i=ips,ipe1
1587         lsum=max(lsum,a(i,k,j))
1588       enddo
1589     enddo
1590   enddo
1591 elseif(fun.eq.RNRM_MAX)then
1592   void=0.
1593   lsum=void
1594   do j=jps,jpe1
1595     do k=kps,kpe1
1596       do i=ips,ipe1
1597         lsum=max(lsum,sqrt(a(i,k,j)*a(i,k,j)+b(i,k,j)*b(i,k,j)))
1598       enddo
1599     enddo
1600   enddo
1601 else
1602   call crash('fun_real: bad fun')
1603 endif
1604
1605 if(lsum.ne.lsum)call crash('fun_real: NaN detected')
1606
1607 dosum=fun.eq.REAL_SUM.or.fun.eq.RNRM_SUM
1608 domax=fun.eq.REAL_MAX.or.fun.eq.RNRM_MAX
1609
1610 ! get process sum over all threads
1611 !$OMP SINGLE
1612 ! only one thread should write to shared variable
1613 psum=void
1614 !$OMP END SINGLE
1615 !$OMP BARRIER
1616 ! now all threads know psum
1617
1618 !$OMP CRITICAL(RDSUM)
1619 ! each thread adds its own lsum
1620 if(dosum)psum=psum+lsum
1621 if(domax)psum=max(psum,lsum)
1622 !$OMP END CRITICAL(RDSUM)
1623
1624 ! wait till all theads are done
1625 !$OMP BARRIER
1626
1627 ! get global sum over all processes
1628 !$OMP SINGLE
1629 ! only one threads will do the mpi communication
1630 #ifdef DM_PARALLEL
1631     if(dosum) gsum = wrf_dm_sum_real ( psum )
1632     if(domax) gsum = wrf_dm_max_real ( psum )
1633 #else
1634     gsum = psum
1635 #endif
1636 if(gsum.ne.gsum)call crash('fun_real: NaN detected')
1637 !$OMP END SINGLE
1638
1639 !$OMP BARRIER
1640 ! now gsum is known to all threads
1641
1642 fun_real=gsum
1643
1644 end function fun_real
1645
1646 end module module_fr_fire_util