Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_fr_sfire_util.F @ 1544

Last change on this file since 1544 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 41.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format