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