wrfv2_fire/phys/module_fr_sfire_util.F

   1 !
   2 !*** Jan Mandel August-October 2007
   3 !*** email: jmandel@ucar.edu or Jan.Mandel@gmail.com or Jan.Mandel@cudenver.edu
   4 !
   5 ! This module contains general purpose utilities and WRF wrappers because I want the
   6 ! model to be able to run standalone. No physics here.
   7 ! Some are dependent on WRF indexing scheme. Some violate WRF conventions but these
   8 ! are not called from the WRF dependent code. Some are not called at all.
   9 !
  10
  11 module module_fr_sfire_util
  12
  13 ! various method selection parameters
  14 ! 1. add the parameter and its static default here
  15 ! 2. add copy from config_flags in module_fr_sfire_driver%%set_flags
  16 ! 3. (optional) add a line in Registry.EM to define the variable and set default value
  17
  18
  19 integer,save::          &
  20  fire_print_msg=1,      &  ! print SFIRE progress
  21  fire_print_file=1,     &  ! write many files by write_array_m; compile with DEBUG_OUT, do not run in parallel
  22  fuel_left_method=1,    &  ! 1=simple, 2=exact in linear case
  23  fuel_left_irl=2,       &  ! refinement for fuel calculation, must be even
  24  fuel_left_jrl=2,       &
  25  boundary_guard=-1,     &  ! crash if fire gets this many cells to domain boundary, -1=off
  26  fire_grows_only=1,     &  ! fire can spread out only (level set functions may not increase)
  27  fire_upwinding=3,      &  ! upwind normal spread: 1=standard, 2=godunov, 3=eno, 4=sethian
  28  fire_upwind_split=0,   &  ! 1=upwind advection separately from normal direction spread
  29  fire_test_steps=0,     &  ! 0=no fire, 1=normal, >1 = do specified number of steps and terminate (testing only)
  30  fire_topo_from_atm=1      ! 0 = expect ZSF set correctly on entry, 1 = populate by interploating from atmosphere
  31
  32
  33
  34
  35 real, save::            &
  36  fire_atm_feedback=1. , &  ! 1 = normal, 0. = one way coupling atmosphere -> fire only
  37  fire_back_weight=0.5,  &  ! RK parameter, 0 = Euler method, 0.5 = Heun, 1 = fake backward Euler
  38  fire_viscosity=0.4,    &  ! artificial viscosity
  39  fire_lfn_ext_up=1.        ! 0.=extend level set function at boundary by reflection, 1.=always up
  40
  41 contains
  42
  43 !
  44 !****************
  45 !
  46 subroutine crash(s)
  47 use module_wrf_error
  48 implicit none
  49 character(len=*), intent(in)::s
  50 character(len=128)msg
  51 msg='SFIRE:'//s
  52 call wrf_error_fatal(msg)
  53 end subroutine crash
  54
  55 !
  56 !****************
  57 !
  58
  59 subroutine message(s)
  60 use module_wrf_error
  61 !$ use OMP_LIB
  62 implicit none
  63 character(len=*), intent(in)::s
  64 character(len=128)msg
  65 character(len=1)thread
  66 integer m
  67 m=-1
  68 !$     m=omp_get_thread_num()
  69 if (m<0) then
  70     msg='SFIRE:'//s
  71 else
  72     write(thread,'(i1,a)')m
  73     msg='SFIRE'//thread//s
  74 endif
  75 if(fire_print_msg.gt.0)then
  76 call wrf_message(msg)
  77 endif
  78 end subroutine message
  79
  80 !
  81 !****************
  82 !
  83
  84 subroutine set_ideal_coord( dxf,dyf, &
  85                 ifds,ifde,jfds,jfde,  &
  86                 ifms,ifme,jfms,jfme,  &
  87                 ifts,ifte,jfts,jfte,  &
  88                 fxlong,fxlat          &
  89             )
  90 implicit none
  91 ! arguments
  92 real, intent(in)::dxf,dyf
  93 integer, intent(in):: &
  94                 ifds,ifde,jfds,jfde,  &
  95                 ifms,ifme,jfms,jfme,  &
  96                 ifts,ifte,jfts,jfte
  97 real, intent(out),dimension(ifms:ifme,jfms:jfme)::fxlong,fxlat
  98 ! local
  99 integer::i,j
 100                 ! set fake  coordinates, in m
 101                 do j=jfts,jfte
 102                     do i=ifts,ifte
 103                         ! uniform mesh, lower left domain corner is (0,0)
 104                         fxlong(i,j)=(i-ifds+0.5)*dxf
 105                         fxlat (i,j)=(j-jfds+0.5)*dyf
 106                     enddo
 107                 enddo
 108 end subroutine set_ideal_coord
 109
 110 !
 111 !****************
 112 !
 113
 114
 115 subroutine continue_at_boundary(ix,iy,bias, & ! do x direction or y direction
 116     ims,ime,jms,jme, &                ! memory dims
 117     ids,ide,jds,jde, &                ! domain dims
 118     ips,ipe,jps,jpe, &                ! patch dims
 119     its,ite,jts,jte, &                ! tile dims
 120     lfn)                              ! array
 121 implicit none
 122 !*** description
 123 ! extend array by one beyond the domain by linear continuation
 124 !*** arguments
 125 integer, intent(in)::ix,iy              ! not 0 = do x or y (1 or 2) direction
 126 real,intent(in)::bias                   ! 0=none, 1.=max
 127 integer, intent(in)::ims,ime,jms,jme, &                ! memory dims
 128     ids,ide,jds,jde, &                ! domain dims
 129     ips,ipe,jps,jpe, &                ! patch dims
 130     its,ite,jts,jte                 ! tile dims
 131 real,intent(inout),dimension(ims:ime,jms:jme)::lfn
 132 !*** local
 133 integer i,j,itso,jtso,iteo,jteo
 134 character(len=128)::msg
 135 integer::its1,ite1,jts1,jte1
 136 integer,parameter::halo=1           ! width of halo region to update
 137 !*** executable
 138 ! check if there is space for the extension
 139 call check_mesh_2dim(its-1,ite+1,jts-1,jte+1,ims,ime,jms,jme)
 140 ! for dislay only
 141 itso=its
 142 jtso=jts
 143 iteo=ite
 144 jteo=jte
 145 ! go halo width beyond if at patch boundary but not at domain boudnary
 146 ! assume we have halo need to compute the value we do not have
 147 ! the next thread that would conveniently computer the extended values at patch corners
 148 ! besides halo may not transfer values outside of the domain
 149 !
 150 its1=its
 151 jts1=jts
 152 ite1=ite
 153 jte1=jte
 154 if(its.eq.ips.and..not.its.eq.ids)its1=its-halo
 155 if(jts.eq.jps.and..not.jts.eq.jds)jts1=jts-halo
 156 if(ite.eq.ipe.and..not.ite.eq.ide)ite1=ite+halo
 157 if(jte.eq.jpe.and..not.jte.eq.jde)jte1=jte+halo
 158 write(msg,'(a,2i5,a,f5.2)')'continue_at_boundary: directions',ix,iy,' bias ',bias
 159 call message(msg)
 160 if(ix.ne.0)then
 161     if(its.eq.ids)then
 162         do j=jts1,jte1
 163             lfn(ids-1,j)=EX(lfn(ids,j),lfn(ids+1,j))
 164         enddo
 165         itso=ids-1
 166     endif
 167     if(ite.eq.ide)then
 168         do j=jts1,jte1
 169             lfn(ide+1,j)=EX(lfn(ide,j),lfn(ide-1,j))
 170         enddo
 171         iteo=ide+1
 172     endif
 173     write(msg,'(8(a,i5))')'continue_at_boundary: x:',its,':',ite,',',jts,':',jte,' ->',itso,':',iteo,',',jts1,':',jte1
 174     call message(msg)
 175 endif
 176 if(iy.ne.0)then
 177     if(jts.eq.jds)then
 178         do i=its1,ite1
 179             lfn(i,jds-1)=EX(lfn(i,jds),lfn(i,jds+1))
 180         enddo
 181         jtso=jds-1
 182     endif
 183     if(jte.eq.jde)then
 184         do i=its1,ite1
 185             lfn(i,jde+1)=EX(lfn(i,jde),lfn(i,jde-1))
 186         enddo
 187         jteo=jde+1
 188     endif
 189     write(msg,'(8(a,i5))')'continue_at_boundary: y:',its,':',ite,',',jts,':',jte,' ->',its1,':',ite1,',',jtso,':',jteo
 190     call message(msg)
 191 endif
 192 ! corners of the domain
 193 if(ix.ne.0.and.iy.ne.0)then
 194     if(its.eq.ids.and.jts.eq.jds)lfn(ids-1,jds-1)=EX(lfn(ids,jds),lfn(ids+1,jds+1))
 195     if(its.eq.ids.and.jte.eq.jde)lfn(ids-1,jde+1)=EX(lfn(ids,jde),lfn(ids+1,jde-1))
 196     if(ite.eq.ide.and.jts.eq.jds)lfn(ide+1,jds-1)=EX(lfn(ide,jds),lfn(ide-1,jds+1))
 197     if(ite.eq.ide.and.jte.eq.jde)lfn(ide+1,jde+1)=EX(lfn(ide,jde),lfn(ide-1,jde-1))
 198 endif
 199 return
 200 contains
 201 real function EX(a,b)
 202 !*** statement function
 203 real a,b
 204 EX=(1.-bias)*(2.*a-b)+bias*max(2.*a-b,a,b)   ! extrapolation, max quarded
 205 end function EX
 206 end subroutine continue_at_boundary
 207
 208 !
 209 !*****************************
 210 !
 211
 212 subroutine check_mesh_2dim(ids,ide,jds,jde,ims,ime,jms,jme)
 213 implicit none
 214 integer, intent(in)::ids,ide,jds,jde,ims,ime,jms,jme
 215 character(len=128)msg
 216 if(ids<ims.or.ide>ime.or.jds<jms.or.jde>jme)then
 217     write(msg,*)'mesh dimensions:  ',ids,ide,jds,jde
 218     call message(msg)
 219     write(msg,*)'memory dimensions:',ims,ime,jms,jme
 220     call message(msg)
 221     call crash('check_mesh_2dim: memory dimensions too small')
 222 endif
 223 end subroutine check_mesh_2dim
 224
 225 !
 226 !****************
 227 !
 228
 229 subroutine check_mesh_3dim(ids,ide,kds,kde,jds,jde,ims,ime,kms,kme,jms,jme)
 230 integer, intent(in)::ids,ide,jds,jde,ims,ime,jms,jme,kds,kde,kms,kme
 231 if(ids<ims.or.ide>ime.or.jds<jms.or.jde>jme.or.kds<kms.or.kde>kme) then
 232     call crash('memory dimensions too small')
 233 endif
 234 end subroutine check_mesh_3dim
 235
 236 !
 237 !****************
 238 !
 239
 240 subroutine sum_2d_cells(       &
 241        ims2,ime2,jms2,jme2,    &
 242        its2,ite2,jts2,jte2,    &
 243        v2,                     &  ! input
 244        ims1,ime1,jms1,jme1,    &
 245        its1,ite1,jts1,jte1,    &
 246        v1)                        ! output
 247 implicit none
 248
 249 !*** purpose
 250 ! sum cell values in mesh2 to cell values of coarser mesh1
 251
 252 !*** arguments
 253 ! the dimensions are in cells, not nodes!
 254
 255 integer, intent(in)::its1,ite1,jts1,jte1,ims1,ime1,jms1,jme1
 256 real, intent(out)::v1(ims1:ime1,jms1:jme1)
 257 integer, intent(in)::its2,ite2,jts2,jte2,ims2,ime2,jms2,jme2
 258 real, intent(in)::v2(ims2:ime2,jms2:jme2)
 259
 260 !*** local
 261 integer:: i1,i2,j1,j2,ir,jr,isz1,isz2,jsz1,jsz2,ioff,joff,ibase,jbase
 262 real t
 263 character(len=128)msg
 264
 265 !*** executable
 266
 267 !check mesh dimensions and domain dimensions
 268 call check_mesh_2dim(its1,ite1,jts1,jte1,ims1,ime1,jms1,jme1)
 269 call check_mesh_2dim(its2,ite2,jts2,jte2,ims2,ime2,jms2,jme2)
 270
 271 ! compute mesh sizes
 272 isz1 = ite1-its1+1
 273 jsz1 = jte1-jts1+1
 274 isz2 = ite2-its2+1
 275 jsz2 = jte2-jts2+1
 276
 277 ! check mesh sizes
 278 if(isz1.le.0.or.jsz1.le.0.or.isz2.le.0.or.jsz2.le.0)then
 279     call message('all mesh sizes must be positive')
 280     goto 9
 281 endif
 282
 283 ! compute mesh ratios
 284 ir=isz2/isz1
 285 jr=jsz2/jsz1
 286
 287 if(isz2.ne.isz1*ir .or. jsz2.ne.jsz1*jr)then
 288     call message('input mesh size must be multiple of output mesh size')
 289     goto 9
 290 endif
 291
 292
 293 ! v1 = sum(v2)
 294 do j1=jts1,jte1
 295     jbase=jts2+jr*(j1-jts1)
 296     do i1=its1,ite1
 297        ibase=its2+ir*(i1-its1)
 298        t=0.
 299        do joff=0,jr-1
 300            j2=joff+jbase
 301            do ioff=0,ir-1
 302                i2=ioff+ibase
 303                t=t+v2(i2,j2)
 304            enddo
 305        enddo
 306        v1(i1,j1)=t
 307     enddo
 308 enddo
 309
 310 return
 311
 312 9 continue
 313 write(msg,91)its2,ite2,jts2,jte2,ims2,ime2,jms2,jme2
 314 call message(msg)
 315 write(msg,91)its1,ite1,jts1,jte1,ims1,ime1,jms1,jme1
 316 call message(msg)
 317 write(msg,92)'input  mesh size:',isz2,jsz2
 318 call message(msg)
 319 91 format('dimensions: ',8i8)
 320 write(msg,92)'output mesh size:',isz1,jsz1
 321 call message(msg)
 322 92 format(a,2i8)
 323 call crash('module_fr_spread_util:sum_mesh_cells: bad mesh sizes')
 324
 325 end subroutine sum_2d_cells
 326
 327
 328
 329 ! module_fr_sfire_util%%interpolate_2d
 330 subroutine interpolate_2d(  &
 331     ims2,ime2,jms2,jme2, & ! array coarse grid
 332     its2,ite2,jts2,jte2, & ! dimensions coarse grid
 333     ims1,ime1,jms1,jme1, & ! array coarse grid
 334     its1,ite1,jts1,jte1, & ! dimensions fine grid
 335     ir,jr,               & ! refinement ratio
 336     rip2,rjp2,rip1,rjp1, & ! (rip2,rjp2) on grid 2 lines up with (rip1,rjp1) on grid 1
 337     v2, &                  ! in coarse grid
 338     v1  )                  ! out fine grid
 339 implicit none
 340
 341 !*** purpose
 342 ! interpolate nodal values in mesh2 to nodal values in mesh1
 343 ! input mesh 2 is coarse output mesh 1 is fine
 344 ! interpolation runs over mesh2 region its2:ite2,jts2:jte2
 345 ! only the part of mesh 1 region its1:ite1,jts1:jte1 is modified
 346
 347 !*** arguments
 348
 349 integer, intent(in)::its1,ite1,jts1,jte1,ims1,ime1,jms1,jme1
 350 integer, intent(in)::its2,ite2,jts2,jte2,ims2,ime2,jms2,jme2
 351 integer, intent(in)::ir,jr
 352 real,intent(in):: rjp1,rip1,rjp2,rip2
 353 real, intent(out)::v1(ims1:ime1,jms1:jme1)
 354 real, intent(in)::v2(ims2:ime2,jms2:jme2)
 355
 356 !*** local
 357 integer:: i1,i2,j1,j2,is,ie,js,je
 358 real:: tx,ty,rx,ry
 359 real:: rio,rjo
 360 intrinsic::ceiling,floor
 361
 362 !*** executable
 363
 364 !check mesh dimensions and domain dimensions
 365 call check_mesh_2dim(its1,ite1,jts1,jte1,ims1,ime1,jms1,jme1)
 366 call check_mesh_2dim(its2,ite2,jts2,jte2,ims2,ime2,jms2,jme2)
 367
 368 ! compute mesh ratios
 369 rx=1./ir
 370 ry=1./jr
 371
 372 do j2=jts2,jte2-1       ! loop over coarse mesh cells
 373     rjo=rjp1+jr*(j2-rjp2)  ! fine grid coordinate of the coarse grid patch start
 374     js=max(jts1,ceiling(rjo))     ! lower bound of fine grid patch for this cell
 375     je=min(jte1,floor(rjo)+jr)  ! upper bound of fine grid patch for this cell
 376     do i2=its2,ite2-1
 377         rio=rip1+ir*(i2-rip2)
 378         is=max(its1,ceiling(rio))
 379         ie=min(ite1,floor(rio)+ir)
 380         do j1=js,je
 381             ty = (j1-rjo)*ry
 382             do i1=is,ie
 383                 ! in case fine grid lies on the boundary of several cells
 384                 ! the result will be written multiple times with the same value
 385                 ! up to a rounding error
 386                 tx = (i1-rio)*rx
 387                 !print *,'coarse ',i2,j2,'to',i2+1,j2+1,' fine ',is,js,' to ',ie,je
 388                 v1(i1,j1)=                     &
 389                       (1-tx)*(1-ty)*v2(i2,j2)  &
 390                  +    (1-tx)*ty  *v2(i2,j2+1)  &
 391                  +      tx*(1-ty)*v2(i2+1,j2)  &
 392                  +        tx*ty  *v2(i2+1,j2+1)
 393                 !print *,'coarse ',i2,j2,' fine ',i1,j1, ' offset ',io,jo,' weights ',tx,ty, &
 394                 ! 'in ',v2(i2,j2),v2(i2,j2+1),v2(i2+1,j2),v2(i2+1,j2+1),' out ',v1(i1,j1)
 395                 !write(*,'(a,2i5,a,2f8.2,a,4f8.2,a,2i5,a,f8.2)') &
 396                 !'coarse ',i2,j2,' coord',rio,rjo,' val',v2(i2,j2),v2(i2,j2+1),v2(i2+1,j2),v2(i2+1,j2+1),&
 397                 !' fine ',i1,j1,' out ',v1(i1,j1)
 398            enddo
 399        enddo
 400     enddo
 401 enddo
 402
 403 end subroutine interpolate_2d
 404
 405 !
 406 !****************
 407 !
 408
 409 subroutine interpolate_2d_cells2cells(              &
 410       ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, &  ! in
 411       ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,v1  )  ! out
 412 implicit none
 413
 414 !*** purpose
 415 ! interpolate nodal values in mesh2 to nodal values in mesh1
 416 ! input mesh 2 is coarse output mesh 1 is fine
 417
 418 !*** arguments
 419
 420 integer, intent(in)::ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1
 421 real, intent(out)::v1(ims1:ime1,jms1:jme1)
 422 integer, intent(in)::ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2
 423 real, intent(in)::v2(ims2:ime2,jms2:jme2)
 424
 425 ! Example with mesh ratio=4. | = cell boundary,  x = cell center
 426 !
 427 !  mesh2   |-------x-------|-------x-------|
 428 !  mesh1   |-x-|-x-|-x-|-x-|-x-|-x-|-x-|-x-|
 429 !
 430
 431 !*** local
 432 integer:: ir,jr,isz1,isz2,jsz1,jsz2,ip,jp,ih,jh
 433 character(len=128)msg
 434
 435 !*** executable
 436
 437 !check mesh dimensions and domain dimensions
 438 call check_mesh_2dim(ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1)
 439 call check_mesh_2dim(ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2)
 440
 441 ! compute mesh sizes
 442 isz1 = ide1-ids1+1
 443 jsz1 = jde1-jds1+1
 444 isz2 = ide2-ids2+1
 445 jsz2 = jde2-jds2+1
 446
 447 ! check mesh sizes
 448 if(isz1.le.0.or.jsz1.le.0.or.isz2.le.0.or.jsz2.le.0)goto 9
 449 if(mod(isz1,isz2).ne.0.or.mod(jsz1,jsz2).ne.0)goto 9
 450
 451 ! compute mesh ratios
 452 ir=isz1/isz2
 453 jr=jsz1/jsz2
 454 !
 455 !  mesh2   |-------x-------|-------x-------|
 456 !  mesh1   |-x-|-x-|-x-|-x-|-x-|-x-|-x-|-x-|
 457
 458 !  mesh2   |-----x-----|-----x-----|  rx=3
 459 !  mesh1   |-x-|-x-|-x-|-x-|-x-|-x-|
 460 !  i2            1   1   1   2
 461 !  i1        1   2   3   4   5
 462 !  ioff          0   1   2   0
 463 !  tx            0  1/3 2/3
 464
 465 !  mesh2   |---x---|---x---| rx=2
 466 !  mesh1   |-x-|-x-|-x-|-x-|
 467 !  i2            1   1   2
 468 !  i1            2   3   4
 469 !  ioff          0   1   2
 470 !  tx           1/4 3/4
 471
 472
 473 ! offset of the last node in the 1st half of the cell
 474 ih=ir/2
 475 jh=jr/2
 476 ! 0 if coarse cell center coincides with fine, 1 if not
 477 ip=mod(ir+1,2)
 478 jp=mod(jr+1,2)
 479
 480 call interpolate_2d_w(ip,jp,ih,jh,ir,jr,              &
 481       ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, &  ! in
 482       ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,v1  )  ! out
 483
 484 return
 485
 486 9 continue
 487 write(msg,91)ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2
 488 call message(msg)
 489 write(msg,91)ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1
 490 call message(msg)
 491 write(msg,92)'input  mesh size:',isz2,jsz2
 492 call message(msg)
 493 91 format('dimensions: ',8i8)
 494 write(msg,92)'output mesh size:',isz1,jsz1
 495 call message(msg)
 496 92 format(a,2i8)
 497 call crash("module_fr_sfire_util:interpolate_2dmesh_cells: bad mesh sizes")
 498 end subroutine interpolate_2d_cells2cells
 499
 500 !
 501 !****************
 502 !
 503
 504 subroutine interpolate_2d_cells2nodes(              &
 505       ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, &  ! in
 506       ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,v1  )  ! out
 507 implicit none
 508
 509 !*** purpose
 510 ! interpolate nodal values in mesh2 to nodal values in mesh1
 511 ! input mesh 2 is coarse output mesh 1 is fine
 512
 513 !*** arguments
 514
 515 integer, intent(in)::ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1
 516 real, intent(out)::v1(ims1:ime1,jms1:jme1)
 517 integer, intent(in)::ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2
 518 real, intent(in)::v2(ims2:ime2,jms2:jme2)
 519
 520 ! Example with mesh ratio=4. | = cell boundary,  x = cell center
 521 !
 522 !  mesh2   |-------x-------|-------x-------|
 523 !  mesh1   x-|-x-|-x-|-x-|-x-|-x-|-x-|-x-|-x
 524 !
 525
 526 !*** local
 527 integer:: ir,jr,isz1,isz2,jsz1,jsz2,ip,jp,ih,jh
 528 character(len=128)msg
 529
 530 !*** executable
 531
 532 !check mesh dimensions and domain dimensions
 533 call check_mesh_2dim(ids1,ide1+1,jds1,jde1+1,ims1,ime1,jms1,jme1)
 534 call check_mesh_2dim(ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2)
 535
 536 ! compute mesh sizes
 537 isz1 = ide1-ids1+1
 538 jsz1 = jde1-jds1+1
 539 isz2 = ide2-ids2+1
 540 jsz2 = jde2-jds2+1
 541
 542 ! check mesh sizes
 543 if(isz1.le.0.or.jsz1.le.0.or.isz2.le.0.or.jsz2.le.0)goto 9
 544 if(mod(isz1,isz2).ne.0.or.mod(jsz1,jsz2).ne.0)goto 9
 545
 546 ! compute mesh ratios
 547 ir=isz1/isz2
 548 jr=jsz1/jsz2
 549 !
 550 !  mesh2   |-------x-------|-------x-------|
 551 !  mesh1   x-|-x-|-x-|-x-|-x-|-x-|-x-|-x-|-x
 552
 553 !  mesh2   |-----x-----|-----x-----|  rx=3
 554 !  mesh1   x-|-x-|-x-|-x-|-x-|-x-|-x
 555
 556 !  mesh2   |---x---|---x---| rx=2
 557 !  mesh1   x-|-x-|-x-|-x-|-x
 558
 559 ! offset of the last node in the 1st half of the cell
 560 ih=(ir+1)/2
 561 jh=(jr+1)/2
 562 ! 0 if coarse cell center coincides with fine, 1 if not
 563 ip=mod(ir,2)
 564 jp=mod(jr,2)
 565
 566
 567 call interpolate_2d_w(ip,jp,ih,jh,ir,jr,              &
 568       ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, &  ! in
 569       ids1,ide1+1,jds1,jde1+1,ims1,ime1,jms1,jme1,v1  )  ! out
 570
 571
 572 return
 573 9 continue
 574 write(msg,91)ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2
 575 call message(msg)
 576 write(msg,91)ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1
 577 call message(msg)
 578 write(msg,92)'input  mesh size:',isz2,jsz2
 579 call message(msg)
 580 91 format('dimensions: ',8i8)
 581 write(msg,92)'output mesh size:',isz1,jsz1
 582 call message(msg)
 583 92 format(a,2i8)
 584 call crash("module_fr_sfire_util:interpolate_2d_cells2nodes: bad mesh sizes")
 585 end subroutine interpolate_2d_cells2nodes
 586 !
 587 !****************
 588 !
 589
 590 subroutine interpolate_2d_w(ip,jp,ih,jh,ir,jr,             &
 591       ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, &  ! in
 592       ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,v1  )  ! out
 593 implicit none
 594 !*** EXCEPTION: THIS SUBROUTINE IS NEITHER CELL NOR NODE BASED.
 595
 596 integer, intent(in)::ip,jp,ih,jh,ir,jr
 597 integer, intent(in)::ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1
 598 real, intent(out)::v1(ims1:ime1,jms1:jme1)
 599 integer, intent(in)::ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2
 600 real, intent(in)::v2(ims2:ime2,jms2:jme2)
 601
 602 real:: tx,ty,rx,ry,half,xoff,yoff
 603 integer:: i1,i2,j1,j2,ioff,joff
 604 parameter(half=0.5)
 605
 606 rx=ir
 607 ry=jr
 608
 609 xoff = ip*half
 610 yoff = jp*half
 611
 612 ! the inside, ids1+ih:ide1-ih,jds1+jh:jde1-jh
 613 do j2=jds2,jde2-1     ! interpolate from nodes j2 and j2+1
 614     do i2=ids2,ide2-1
 615         do ioff=0,ir-ip
 616             do joff=0,jr-jp
 617                 ! compute fine mesh index
 618                 i1=ioff+(ih+ids1)+ir*(i2-ids2)
 619                 j1=joff+(jh+jds1)+jr*(j2-jds2)
 620                 ! weights
 621                 tx = (ioff+xoff)/rx
 622                 ty = (joff+yoff)/ry
 623                 ! interpolation
 624                 v1(i1,j1)=                     &
 625                       (1-tx)*(1-ty)*v2(i2,j2)  &
 626                  +    (1-tx)*ty  *v2(i2,j2+1)  &
 627                  +      tx*(1-ty)*v2(i2+1,j2)  &
 628                  +        tx*ty  *v2(i2+1,j2+1)
 629                 !write(*,'(3(a,2i5),a,2f7.4)')'coarse ',i2,j2,' fine ',i1,j1, &
 630                 ! ' offset ',ioff,joff,' weights ',tx,ty
 631                 !write(*,'(a,4f7.4,a,f7.4)')'in ',v2(i2,j2),v2(i2,j2+1),v2(i2+1,j2), &
 632                 !  v2(i2+1,j2+1),' out ',v1(i1,j1)
 633            enddo
 634        enddo
 635     enddo
 636 enddo
 637
 638 ! extend to the boundary strips from the nearest known
 639 do ioff=0,ih-1  ! top and bottom strips
 640     do j2=jds2,jde2-1
 641         do joff=0,jr-jp
 642            j1=joff+(jh+jds1)+jr*(j2-jds2)
 643            ! weights
 644            ty = (joff+yoff)/ry
 645            ! interpolation
 646            v1(ids1+ioff,j1)=(1-ty)*v2(ids2,j2)+ty*v2(ids2,j2+1)
 647            v1(ide1-ioff,j1)=(1-ty)*v2(ide2,j2)+ty*v2(ide2,j2+1)
 648        enddo
 649     enddo
 650 enddo
 651 do joff=0,jh-1  ! left and right strips
 652     do i2=ids2,ide2-1
 653         do ioff=0,ir-ip
 654            i1=ioff+(ih+ids1)+ir*(i2-ids2)
 655            ! weights
 656            tx = (ioff+xoff)/rx
 657            ! interpolation
 658            v1(i1,jds1+joff)=(1-tx)*v2(i2,jds2)+tx*v2(i2+1,jds2)
 659            v1(i1,jde1-joff)=(1-tx)*v2(i2,jde2)+tx*v2(i2+1,jde2)
 660        enddo
 661     enddo
 662 enddo
 663 ! extend to the 4 corner squares from the nearest known
 664 do ioff=0,ih-1
 665     do joff=0,jh-1
 666         v1(ids1+ioff,jds1+joff)=v2(ids2,jds2)
 667         v1(ide1-ioff,jds1+joff)=v2(ide2,jds2)
 668         v1(ids1+ioff,jde1-joff)=v2(ids2,jde2)
 669         v1(ide1-ioff,jde1-joff)=v2(ide2,jde2)
 670     enddo
 671 enddo
 672 end subroutine interpolate_2d_w
 673
 674 !
 675 !****************
 676 !
 677
 678 real function interp(ids,ide,jds,jde,ims,ime,jms,jme,x,y,v)
 679 implicit none
 680 !*** purpose
 681 ! general interpolation in a rectangular
 682
 683 !*** arguments
 684
 685 integer, intent(in)::ids,ide,jds,jde,ims,ime,jms,jme
 686 real, intent(in)::x,y,v(ims:ime,jms:jme)
 687 ! the mesh is cell based so the used dimension of v is ids:ide+1,jds:jde+1
 688
 689 !*** calls
 690 intrinsic floor,min,max
 691
 692 !*** local
 693 integer i,j
 694 real tx,ty
 695
 696 ! executable
 697
 698 ! indices of the lower left corner of the cell in the mesh that contains (x,y)
 699 i = floor(x)
 700 i=max(min(i,ide),ids)
 701 j = floor(y)
 702 j=max(min(j,jde),jds)
 703
 704 ! the leftover
 705 tx = x - real(i)
 706 ty = y - real(j)
 707
 708 ! interpolate the values
 709 interp = &
 710                     (1-tx)*(1-ty)*v(i,j)    &
 711                  +    tx*(1-ty)  *v(i+1,j)  &
 712                  +    (1-tx)*ty  *v(i,j+1)  &
 713                  +        tx*ty  *v(i+1,j+1)
 714
 715 !print *,'x,y=',x,y,'i1,i2=',i1,j1,'tx,ty=',tx,ty,' interp=',interp
 716 end function interp
 717
 718 subroutine meshdiffc_2d(ids, ide, jds,jde , &    ! mesh area used (in cells, end +1)
 719                    ims1,ime1,jms1,jme1, &       ! memory dimensiuons
 720                    dx,dy,               &       ! mesh spacing
 721                    lfn,                 &       ! input
 722                    diffCx,diffCy) ! output
 723 implicit none
 724
 725 !*** purpose
 726 ! central differences on a 2d mesh
 727
 728 !*** arguments
 729
 730 integer, intent(in)::ids,ide,jds,jde,ims1,ime1,jms1,jme1
 731 real, intent(in):: dx,dy
 732 real, intent(in), dimension(ims1:ime1,jms1:jme1):: lfn
 733 real, intent(out), dimension(ims1:ime1,jms1:jme1):: diffCx,diffCy
 734
 735 !*** local
 736 integer:: i,j
 737 real, dimension(ims1:ime1,jms1:jme1):: diffLx,diffRx,diffLy,diffRy
 738
 739 ! get one-sided differences; dumb but had that already...
 740 call meshdiff_2d(ids, ide, jds,jde , &    ! mesh area used (in cells, end +1)
 741                    ims1,ime1,jms1,jme1, &       ! dimensions of lfn
 742                    dx,dy,               &       ! mesh spacing
 743                    lfn,                 &       ! input
 744                    diffLx,diffRx,diffLy,diffRy) ! output
 745
 746 ! make into central
 747 do j=jds,jde+1
 748     do i=ids,ide+1
 749         diffCx(i,j)=0.5*(diffLx(i,j) + diffRx(i,j))
 750         diffCy(i,j)=0.5*(diffLy(i,j) + diffRy(i,j))
 751     enddo
 752 enddo
 753 end subroutine meshdiffc_2d
 754
 755 subroutine meshdiff_2d(ids, ide, jds,jde , &    ! mesh area used (in cells, end +1)
 756                    ims1,ime1,jms1,jme1, &       ! dimensions of lfn
 757                    dx,dy,               &       ! mesh spacing
 758                    lfn,                 &       ! input
 759                    diffLx,diffRx,diffLy,diffRy) ! output
 760 implicit none
 761
 762 !*** purpose
 763 ! one-sided differences on a 2d mesh
 764
 765 !*** arguments
 766
 767 integer, intent(in)::ids,ide,jds,jde,ims1,ime1,jms1,jme1
 768 real, intent(in):: dx,dy
 769 real, intent(in), dimension(ims1:ime1,jms1:jme1):: lfn
 770 real, intent(out), dimension(ims1:ime1,jms1:jme1):: diffLx,diffRx,diffLy,diffRy
 771
 772 !*** local
 773 integer:: i,j
 774 real:: tmpx,tmpy
 775
 776 !*** executable
 777
 778     call check_mesh_2dim(ids,ide+1,jds,jde+1,ims1,ime1,jms1,jme1)
 779
 780     ! the bulk of the work
 781     do j=jds,jde
 782         do i=ids,ide
 783             tmpx = (lfn(i+1,j)-lfn(i,j))/dx
 784             diffLx(i+1,j) = tmpx
 785             diffRx(i,j)   = tmpx
 786             tmpy = (lfn(i,j+1)-lfn(i,j))/dy
 787             diffLy(i,j+1) = tmpy
 788             diffRy(i,j)   = tmpy
 789         enddo
 790         ! missing values - put there the other one
 791         diffLx(ids,j)  = diffLx(ids+1,j)
 792         diffRx(ide+1,j)= diffRx(ide,j)
 793     enddo
 794     ! cleanup
 795     ! j=jde+1 from above loop
 796     do i=ids,ide
 797         tmpx = (lfn(i+1,j)-lfn(i,j))/dx
 798         diffLx(i+1,j) = tmpx
 799         diffRx(i,j)   = tmpx
 800     enddo
 801     ! i=ide+1 from above loop
 802     do j=jds,jde
 803         tmpy = (lfn(i,j+1)-lfn(i,j))/dy
 804         diffLy(i,j+1) = tmpy
 805         diffRy(i,j)   = tmpy
 806     enddo
 807     ! missing values - put there the other one
 808     ! j=jde+1 from above loop, j=jds:jde done before in main bulk loop
 809     diffLx(ids,j)   = diffLx(ids+1,j)
 810     diffRx(ide+1,j) = diffRx(ide,j)
 811     do i=ids,ide+1
 812         diffLy(i,jds)   = diffLy(i,jds+1)
 813         diffRy(i,jde+1) = diffRy(i,jde)
 814     enddo
 815
 816 end subroutine meshdiff_2d
 817
 818
 819
 820
 821 real pure function sum_2darray( its,ite,jts,jte,               &
 822                                 ims,ime,jms,jme,               &
 823                                 a)
 824 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme
 825 real, intent(in)::a(ims:ime,jms:jme)
 826 !*** local
 827 integer:: i,j
 828 real:: t
 829 t=0.
 830 do j=jts,jte
 831     do i=its,ite
 832         t=t+a(i,j)
 833     enddo
 834 enddo
 835 sum_2darray = t
 836 end function sum_2darray
 837
 838 real pure function max_2darray( its,ite,jts,jte,               &
 839                                 ims,ime,jms,jme,               &
 840                                 a)
 841 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme
 842 real, intent(in)::a(ims:ime,jms:jme)
 843 !*** local
 844 integer:: i,j
 845 real:: t
 846 t=0.
 847 do j=jts,jte
 848     do i=its,ite
 849         t=max(t,a(i,j))
 850     enddo
 851 enddo
 852 max_2darray = t
 853 end function max_2darray
 854
 855 subroutine print_2d_stats_vec(ips,ipe,jps,jpe, &
 856                          ims,ime,jms,jme, &
 857                          ax,ay,name)
 858 implicit none
 859 integer, intent(in)::ips,ipe,jps,jpe,ims,ime,jms,jme
 860 real, intent(in), dimension(ims:ime,jms:jme)::ax,ay
 861 character(len=*),intent(in)::name
 862 integer:: i,j
 863 real:: t
 864 real:: avg_a,max_a,min_a
 865 character(len=25)id
 866 id=name
 867 call print_2d_stats(ips,ipe,jps,jpe, &
 868                          ims,ime,jms,jme, &
 869                          ax,id//'/x ')
 870 call print_2d_stats(ips,ipe,jps,jpe, &
 871                          ims,ime,jms,jme, &
 872                          ay,id//'/y ')
 873 avg_a=0
 874 max_a=-huge(max_a)
 875 min_a= huge(min_a)
 876 do j=jps,jpe
 877     do i=ips,ipe
 878         t=sqrt(ax(i,j)**2+ay(i,j)**2)
 879         max_a=max(max_a,t)
 880         min_a=min(min_a,t)
 881         avg_a=avg_a+t
 882     enddo
 883 enddo
 884 avg_a = avg_a/((ipe-ips+1)*(jpe-jps+1))
 885 call print_stat_line(id//'/sz',ips,ipe,jps,jpe,min_a,max_a,avg_a)
 886 end subroutine print_2d_stats_vec
 887
 888
 889 subroutine print_stat_line(name,ips,ipe,jps,jpe,min_a,max_a,avg_a)
 890 !*** encapsulate line with statistics
 891 implicit none
 892 !*** arguments
 893 integer, intent(in)::ips,ipe,jps,jpe
 894 character(len=*),intent(in)::name
 895 real,intent(in)::min_a,max_a,avg_a
 896 !*** local
 897 character(len=128)msg
 898 character(len=24)id
 899 if(fire_print_msg.eq.0)return
 900 id=name
 901 write(msg,'(a,4i4,3g11.3)')id,ips,ipe,jps,jpe,min_a,max_a,avg_a
 902 call message(msg)
 903 if(.not.avg_a.eq.avg_a)call crash('NaN detected')
 904 end subroutine print_stat_line
 905
 906
 907 subroutine print_3d_stats(ips,ipe,kps,kpe,jps,jpe, &
 908                          ims,ime,kms,kme,jms,jme, &
 909                          a,name)
 910 implicit none
 911 integer, intent(in)::ips,ipe,jps,jpe,ims,ime,jms,jme,kms,kme,kps,kpe
 912 real, intent(in)::a(ims:ime,kms:kme,jms:jme)
 913 character(len=*),intent(in)::name
 914 integer:: i,j,k
 915 real:: avg_a,max_a,min_a,t,aa
 916 character(len=128)msg
 917 if(fire_print_msg.eq.0)return
 918 avg_a=0
 919 max_a=-huge(max_a)
 920 min_a= huge(min_a)
 921 t=huge(t)
 922 do j=jps,jpe
 923   do k=kps,kpe
 924     do i=ips,ipe
 925         aa=a(i,k,j)
 926         if(aa.ne.aa)then
 927             write(msg,1)name,i,k,j,'NaN'
 928             call crash(msg)
 929         endif
 930         if(.not.aa.le.t)then
 931             write(msg,1)name,i,k,j,'Inf'
 932             call crash(msg)
 933         endif
 934         if(.not.aa.ge.-t)then
 935             write(msg,1)name,i,k,j,'-Inf'
 936             call crash(msg)
 937         endif
 938         max_a=max(max_a,aa)
 939         min_a=min(min_a,aa)
 940         avg_a=avg_a+aa
 941     enddo
 942   enddo
 943 enddo
 944 1 format(a,'(',i6,',',i6,',',i6,') = ',a)
 945 avg_a = avg_a/((ipe-ips+1)*(jpe-jps+1)*(kpe-kps+1))
 946 call print_stat_line(name,ips,ipe,jps,jpe,min_a,max_a,avg_a)
 947 end subroutine print_3d_stats
 948
 949 subroutine print_2d_stats(ips,ipe,jps,jpe, &
 950                          ims,ime,jms,jme, &
 951                          a,name)
 952 implicit none
 953 integer, intent(in)::ips,ipe,jps,jpe,ims,ime,jms,jme
 954 real, intent(in)::a(ims:ime,jms:jme)
 955 character(len=*),intent(in)::name
 956 !!character(len=128)msg
 957 if(fire_print_msg.eq.0)return
 958 call print_3d_stats(ips,ipe,1,1,jps,jpe, &
 959                          ims,ime,1,1,jms,jme, &
 960                          a,name)
 961 !!write(msg,'(2a,z16)')name,' address =',loc(a)
 962 !!call message(msg)
 963 end subroutine print_2d_stats
 964
 965 real pure function avg_2darray( its,ite,jts,jte,               &
 966                                 ims,ime,jms,jme,               &
 967                                 a)
 968 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme
 969 real, intent(in)::a(ims:ime,jms:jme)
 970 !*** local
 971 !*** executable
 972 avg_2darray = sum_2darray( its,ite,jts,jte,               &
 973                            ims,ime,jms,jme,               &
 974                            a)/((ite-its+1)*(jte-jts+1))
 975 end function avg_2darray
 976
 977 real pure function avg_2darray_vec( its,ite,jts,jte,           &
 978                                 ims,ime,jms,jme,               &
 979                                 ax,ay)
 980 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme
 981 real, intent(in), dimension(ims:ime,jms:jme):: ax,ay
 982 !*** local
 983 integer:: i,j
 984 real:: t
 985 t=0.
 986 do j=jts,jte
 987     do i=its,ite
 988         t=t+sqrt(ax(i,j)**2+ay(i,j)**2)
 989     enddo
 990 enddo
 991 t = t/((ite-its+1)*(jte-jts+1))
 992 avg_2darray_vec = t
 993 end function avg_2darray_vec
 994
 995
 996 subroutine print_array(its,ite,jts,jte,           &
 997                          ims,ime,jms,jme,               &
 998                          a,name,id)
 999 ! debug
1000 !*** arguments
1001 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme,id
1002 real, intent(in), dimension(ims:ime,jms:jme):: a
1003 character(len=*),intent(in)::name
1004 !****
1005 integer i,j
1006 character(len=128)::msg
1007 !****
1008 write(msg,*)name,' start ',id,' dim ',its,ite,jts,jte
1009 call message(msg)
1010 do j=jts,jte
1011     do i=its,ite
1012          write(msg,*)i,j,a(i,j)
1013          call message(msg)
1014     enddo
1015 enddo
1016 write(msg,*)name,' end ',id
1017 call message(msg)
1018 end subroutine print_array
1019
1020 subroutine write_array_m(its,ite,jts,jte,           &
1021                          ims,ime,jms,jme,               &
1022                          a,name,id)
1023 ! debug
1024 !*** arguments
1025 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme,id
1026 real, intent(in), dimension(ims:ime,jms:jme):: a
1027 character(len=*),intent(in)::name
1028 !****
1029 call write_array_m3(its,ite,1,1,jts,jte,           &
1030                          ims,ime,1,1,jms,jme,               &
1031                          a,name,id)
1032 end subroutine write_array_m
1033
1034
1035 subroutine write_array_m3(its,ite,kts,kte,jts,jte,           &
1036                          ims,ime,kms,kme,jms,jme,               &
1037                          a,name,id)
1038 use module_dm
1039
1040 implicit none
1041 ! debug
1042 !*** arguments
1043 integer, intent(in)::its,ite,jts,jte,ims,ime,jms,jme,kts,kte,kms,kme,id
1044 real, intent(in), dimension(ims:ime,kms:kme,jms:jme):: a
1045 character(len=*),intent(in)::name
1046 !****
1047 integer i,j,k,iu,ilen,myproc,nprocs
1048 logical op
1049 character(len=128)::fname,msg
1050 !****
1051 if(fire_print_file.eq.0.or.id.le.0)return
1052 call check_mesh_2dim(its,ite,jts,jte,ims,ime,jms,jme)
1053 call wrf_get_nproc (nprocs)
1054 call wrf_get_myproc( myproc )
1055
1056 if(nprocs.eq.1)then
1057     write(fname,3)name,'_',id,'.txt'
1058 else
1059     write(fname,4)name,'_',id,'.',myproc,'.txt'
1060 endif
1061 do ilen=len(fname),2,-1
1062     if(fname(ilen:ilen).ne.' ')goto 19
1063 enddo
1064 19 continue
1065
1066 !$OMP CRITICAL
1067
1068 iu=0
1069 do i=6,99
1070     inquire(unit=i,opened=op)
1071     if(.not.op.and.iu.le.0)iu=i
1072 enddo
1073 if(iu.gt.0)open(iu,file=fname(1:ilen),form='formatted',status='unknown')
1074
1075 !$OMP END CRITICAL
1076 if(iu.le.0)call crash('write_array_m: cannot find available fortran unit')
1077
1078 write(iu,1)real(its)
1079 write(iu,1)real(ite)
1080 write(iu,1)real(jts)
1081 write(iu,1)real(jte)
1082 write(iu,1)real(kts)
1083 write(iu,1)real(kte)
1084 write(iu,1)(((a(i,k,j),i=its,ite),j=jts,jte),k=kts,kte)
1085 close(iu)
1086 write(msg,2)name,'(',its,':',ite,',',jts,':',jte,',', &
1087 kts,':',kte,') -> ',fname(1:ilen)
1088 call message(msg)
1089 return
1090
1091 1 format(e20.12)
1092 2 format(2a,3(i5,a,i5,a),2a)
1093 3 format(a,a,i5.5,a)
1094 4 format(a,a,i5.5,a,i4.4,a)
1095
1096 9 write(msg,'(3a)')'Cannot open file ',fname(1:ilen),' for writing'
1097 call crash(msg)
1098
1099 end subroutine write_array_m3
1100
1101 ! function to go beyond domain boundary if tile is next to it
1102 pure integer function snode(t,d,i)
1103 implicit none
1104 integer, intent(in)::t,d,i
1105 if(t.ne.d)then
1106     snode=t
1107 else
1108     snode=t+i
1109 endif
1110 end function snode
1111
1112 end module module_fr_sfire_util