da_scale_length.f90
References to this file elsewhere.
1 program da_scale_length
2
3 !-----------------------------------------------------------------------
4 ! Purpose: TBD
5 !-----------------------------------------------------------------------
6
7 implicit none
8
9 integer, parameter :: plot_style = 1, &
10 input_unit = 10, &
11 input_scnd = 12, &
12 nmlst_uint = 9, &
13 nt = 62
14
15 integer, parameter :: plot_switch = 0
16
17 real, parameter :: xwb=0.00, xwe=1.00, ywb=0.00, ywe=1.00, &
18 xlb=0.50, ylb=0.97
19
20 integer :: iy, jx, kz, ier, num, n
21
22 real, dimension(:,:,:,:), allocatable :: var
23 real, dimension(:,:,:), allocatable :: tmp
24
25 character(len=7) :: varname
26 real :: cut_dist_km, resolution_km
27 namelist /control_param/ varname, resolution_km, cut_dist_km
28 integer :: dgrid
29
30 call opngks
31
32 call da_setup_color_table
33
34 ! set background white
35
36 call gscr(1, 0, 1.00, 1.00, 1.00)
37
38 open(unit=nmlst_uint, file='namelist.input', status='old')
39 read(unit=nmlst_uint, nml=control_param, iostat=ier)
40 print*,' Doing job for : ',varname
41 print*,' Cut off Grids : ',dgrid
42 if (ier /= 0) stop 'Wrong namelist.input.'
43 close(unit=nmlst_uint)
44 dgrid = cut_dist_km/resolution_km
45
46 read(input_unit, iostat=ier) iy, jx, kz
47 read(input_scnd, iostat=ier) iy, jx, kz
48 if (ier /= 0) stop 'Wrong input file.'
49
50 allocate(var(1:iy, 1:jx, 1:kz, 1:nt))
51 allocate(tmp(1:iy, 1:jx, 1:kz))
52 ier = 0
53 num = 0
54
55
56 do
57 num=num+1
58
59 write(unit=*, fmt='(a, 4i6)') &
60 'iy, jx, kz, num=', iy, jx, kz, num
61
62 if (varname(1:3) == 'PSI' .or. varname(1:3) == 'psi') then
63 read(input_unit) var(:,:,:,num)
64 read(input_unit) tmp(:,:,:)
65 read(input_unit) tmp(:,:,:)
66 read(input_unit) tmp(:,:,:)
67 else if (varname(1:3) == 'CHI' .or. varname(1:3) == 'chi') then
68 read(input_unit) tmp(:,:,:)
69 read(input_unit) var(:,:,:,num)
70 read(input_unit) tmp(:,:,:)
71 read(input_unit) tmp(:,:,:)
72 else if (varname(1:3) == 'P_U' .or. varname(1:3) == 'p_u') then
73 read(input_unit) tmp(:,:,:)
74 read(input_unit) tmp(:,:,:)
75 read(input_unit) var(:,:,:,num)
76 read(input_unit) tmp(:,:,:)
77 else if (varname(1:3) == 'Q_M' .or. varname(1:3) == 'q_m') then
78 read(input_unit) tmp(:,:,:)
79 read(input_unit) tmp(:,:,:)
80 read(input_unit) tmp(:,:,:)
81 read(input_unit) var(:,:,:,num)
82 else
83 do n = 1,4
84 read(input_unit) tmp(:,:,:)
85 end do
86 end if
87
88 if (varname(1:3) == 'RHM' .or. varname(1:3) == 'rhm') then
89 read(input_scnd) tmp(:,:,:)
90 read(input_scnd) tmp(:,:,:)
91 read(input_scnd) tmp(:,:,:)
92 read(input_scnd) var(:,:,:,num)
93 else
94 do n =1,4
95 read(input_scnd) tmp(:,:,:)
96 end do
97 end if
98 do n = 1,6
99 read(input_unit) tmp(:,:,:)
100 read(input_scnd) tmp(:,:,:)
101 end do
102 read(input_unit, iostat=ier) iy, jx, kz
103 read(input_scnd, iostat=ier) iy, jx, kz
104
105 if (ier /= 0) exit
106 end do
107
108 close(input_unit)
109
110 deallocate(tmp)
111
112 call da_process_single_variable(var,varname,dgrid,num,iy,jx,kz,nt)
113
114 call clsgks
115
116 contains
117
118 subroutine da_process_single_variable(var,varname,dgrid,num,iy,jx,kz,nt)
119
120 implicit none
121
122 integer, intent(in) :: dgrid, num,iy,jx,kz,nt
123 real, dimension(iy,jx,kz,nt), intent(in) :: var
124 character(len=7), intent(in) :: varname
125
126 real, dimension(iy,jx,nt) :: plt
127
128 integer :: i, j, k, n, output_unit
129
130 character(len=12) :: flnm
131
132 output_unit = 20
133
134 write(flnm(1:12), fmt='(a,a)') 'sl_print.', varname(5:7)
135
136 open(unit=output_unit, &
137 file=flnm, &
138 form='formatted', &
139 action='write', &
140 access='sequential', &
141 ! position='rewind', &
142 status='replace')
143
144 write(unit=output_unit, fmt='(2a/a)') &
145 varname(1:3), ' scale length', &
146 ' Lvl m scale-length'
147
148 do k=1,kz
149 do n=1,num
150 do j=1,jx
151 do i=1,iy
152 plt(i,j,n)=var(i,j,k,n)
153 end do
154 end do
155 end do
156
157 call da_make_scale_length(plt,num,k,varname, dgrid, &
158 iy,jx,nt,output_unit)
159
160 call plot_it(plt,num,k,varname,iy,jx,nt)
161 end do
162
163 close(unit=output_unit, status='keep')
164
165 end subroutine da_process_single_variable
166
167 subroutine da_make_scale_length(plt,num,k,varname,dgrid, nx,ny,nt,output_unit)
168
169 implicit none
170
171 integer, intent(in) :: dgrid, num,k,nx,ny,nt, &
172 output_unit
173 real, dimension(nx,ny,nt), intent(inout) :: plt
174 character(len=7), intent(in) :: varname
175
176 real :: radius, value, tinv, rmax, rmin, valmin,valmax
177
178 integer :: i, j, ib, jb, ie, je, m, n, nn, nl, nnn
179
180 ! Steps of fitting Gaussian Distribution:
181
182 ! B(r) = B(0) exp(-r**2/(8*s**2) (1)
183 ! Log at both side of (1):
184 ! ln[B(0)/B(r)] = r**2/(8*s**2) (2)
185 ! {8*ln[B(0)/B(r)]}**0.5 = r/s = m * r
186
187 ! Let:
188 ! y(r) = {8*ln[B(0)/B(r)]}**0.5
189 ! m = sum[r * y(r)]/sum[r*r]
190
191
192
193 real(kind=8), dimension(:), allocatable :: yr
194 real(kind=8), dimension(:), allocatable :: nr, r, bb
195
196 real(kind=8) :: criteria, ml, sl, cl, a, b, c, d, e
197
198 call da_zero_mean(plt, nx, ny, nt, num)
199
200 nn=nx+ny
201 allocate( yr(0:nn))
202
203 allocate(nr(0:nn))
204 allocate(r(0:nn))
205 allocate(bb(0:nn))
206
207 yr = 0.0
208
209 nr = 0.0
210 r = 0.0
211 bb = 0.0
212
213 ib = 4
214 jb = 6
215 ie = nx - 4
216 je= ny - 6
217 ! start of the time loop
218
219 j_loop: do j=jb,je
220 i_loop: do i=ib,ie
221 n_loop: do n=jb, min((j+dgrid),je)
222 m_loop: do m=ib, min((i+dgrid),ie)
223 radius=sqrt(real((m-i)*(m-i)+(n-j)*(n-j)))
224 if (radius == 0) then
225 nl = 0
226 else
227 nl = int(radius + 0.5)
228 end if
229 t_loop: do nn=1,num
230 value = plt(m,n,nn)*plt(i,j,nn)
231 bb(nl)=bb(nl) + value
232 nr(nl)=nr(nl) + 1.0
233 r(nl)= r(nl) + radius
234 end do t_loop
235 end do m_loop
236 end do n_loop
237 end do i_loop
238 end do j_loop
239 nn=nx+ny
240
241 write(unit=*, fmt='(2a/)') varname(1:3), &
242 ' dist points radius y(r)'
243 r(0)=0.0
244 nl = 0
245 bb(0)=bb(0)/nr(0)
246 do n=1,nn
247 if (nr(n) < 1.0) exit
248 if (bb(n) <= 0.0 ) then
249 ! print*,n,' ---> -ve corr bb= ',bb(n),' nr= ',nr(n)
250 ! if ( n == 1) cycle
251 exit
252 end if
253 radius=r(n)/nr(n)
254
255 nl=nl+1
256 bb(nl)=bb(n)/nr(n)
257 r(nl)=radius
258 if (bb(nl) < bb(0)) then
259 yr(nl)=sqrt(8.0*(log(bb(0)/bb(nl))))
260 else
261 yr(nl)=0.0
262 end if
263
264 ! r(n)=radius
265 ! nl=nl+1
266
267 write(unit=*, fmt='(i4,f12.0,f20.4,e20.8)') &
268 nl, nr(n), radius, yr(nl)
269
270 ! if (nl > 1 .and. yr(n) > 5.0) exit
271 if (nl > 1 .and. yr(nl) > 3.0) exit
272 end do
273
274 ! a=nr(0)
275 ! b=0.0
276 ! c=0.0
277 d=0.0
278 e=0.0
279 if ( nl == 0 ) stop ' Did not get any point with +ve corr'
280 do n=1,nl
281 ! a=a+nr(n)
282 ! b=b+nr(n)*yr(n)
283 ! c=c+nr(n)* r(n)
284 d=d+nr(n)* r(n)*yr(n)
285 e=e+nr(n)* r(n)* r(n)
286 end do
287
288 ! ml=(a*d-c*b)/(a*e-c*c)
289 ! cl=(b*e-d*c)/(a*e-c*c)
290 print*,' ml num= ',d,' ml denom= ',e
291 ml=d/e
292 print*,'inverse of scale length = ',ml
293 ! cl=0.0
294
295 sl=1.0/ml
296
297 write(unit=*, fmt='(/2a,i4,3e30.8/)') varname(1:3), &
298 ' scale-length at mode:', k, ml, sl
299
300 write(unit=output_unit, fmt='(i4,2e20.8)') k, ml, sl
301
302 if (nl > 1) &
303 call plot_sl(yr,r,nl,nn,ml,cl,k,varname)
304
305 end subroutine da_make_scale_length
306
307 subroutine da_zero_mean(a, nx, ny, nt, nm)
308
309 implicit none
310
311 integer, intent(in) :: nx, ny, nt, nm
312 real, dimension(nx,ny,nt), intent(inout) :: a
313
314 real :: sum
315
316 integer :: i, j, n
317
318 do n=1,nm
319 sum = 0.0
320
321 do j=1,ny
322 do i=1,nx
323 sum=sum+a(i,j,n)
324 end do
325 end do
326
327 sum = sum/real(nx*ny)
328
329 do j=1,ny
330 do i=1,nx
331 a(i,j,n)=a(i,j,n)-sum
332 end do
333 end do
334 end do
335
336 end subroutine da_zero_mean
337
338 end program da_scale_length
339
340 subroutine da_plot_it(plt,num,k,varname,nx,ny,nt,plot_switch)
341
342 implicit none
343
344 integer, intent(in) :: num,k,nx,ny,nt, plot_switch
345 real, dimension(nx,ny,nt), intent(inout) :: plt
346 character(len=7), intent(in) :: varname
347
348 real :: radius, value, xpb, xpe, ypb, ype
349
350 integer :: i, j, ib, jb, ie, je, m, n, mm, nn
351
352 real, dimension(nx,ny,nt) :: pltsqr
353
354 ! JRB hack. I don't know what I'm doing
355 real :: xwb, xwe, ywb, ywe,xlb,ylb
356 integer :: plot_style
357 real :: red
358 ! JRB
359
360 character(len=20) :: pltlab
361
362 ib=4
363 jb=6
364
365 ie=nx-4
366 je=ny-6
367
368 write(pltlab(1:20),fmt='(2a,i5)') &
369 varname(1:3), ' FOR MODE : ', k
370
371 write(unit=*, fmt='(a)') pltlab
372
373 call set(xwb,xwe,ywb,ywe,xwb,xwe,ywb,ywe,plot_style)
374
375 call gsplci(red)
376 call gspmci(red)
377 call gstxci(red)
378
379 call pwritx(xlb,ylb,pltlab,20,1,0,0)
380
381 xpb=0.0
382 xpe=sqrt(real((nx-1)*(nx-1)+(ny-1)*(ny-1)))
383
384 if (plot_switch > 0) then
385 call da_point_plot(plt,num,nx,ny,nt,ib,jb,ie,je, &
386 xpb,xpe,ypb,ype)
387 else
388 call da_line_plot(plt,num,nx,ny,nt,ib,jb,ie,je, &
389 xpb,xpe,ypb,ype)
390 end if
391
392 call gsplci(red)
393 call gspmci(red)
394 call gstxci(red)
395
396 call line(xpb,0.0,xpe,0.0)
397
398 call frame
399
400 end subroutine da_plot_it
401
402 subroutine da_point_plot(plt,num,nx,ny,nt,ib,jb,ie,je, &
403 xpb,xpe,ypb,ype)
404
405 implicit none
406
407 integer, intent(in) :: num,nx,ny,nt, &
408 ib,jb,ie,je
409 real, dimension(nx,ny,nt), intent(in) :: plt
410 real, intent(in) :: xpb,xpe
411 real, intent(out) :: ypb,ype
412
413
414 real :: radius, value
415
416 integer :: i, j, m, n, mm, nn
417
418 real, dimension(nx,ny,nt) :: pltsqr
419
420 character(len=1), parameter :: symbol='.'
421
422 ! JRB hack
423 real :: xfb,xfe,yfb,yfe,blue
424 integer :: plot_style
425 ! JRB
426
427 do n=1,num
428 do j=1,ny
429 do i=1,nx
430 pltsqr(i,j,n)=plt(i,j,n)*plt(i,j,n)
431 end do
432 end do
433 end do
434
435 ype=maxval(pltsqr(ib:ie, jb:je, 1:num))*1.25
436 ypb=-ype
437
438 if (abs(ype - ypb) < 1.0e-5) then
439 call frame
440 stop 'ype - ypb is too small.'
441 end if
442
443 call set(xfb,xfe,yfb,yfe,xpb,xpe,ypb,ype,plot_style)
444
445 call line(xpb,ypb,xpe,ypb)
446
447 call line(xpb,ypb,xpb,ype)
448
449 i=int(xpe) + 1
450 j=2
451
452 m=int(ype-ypb) + 1
453 n=m/10
454
455 ! call perim(i,j,m,m)
456
457 value = ypb+0.02*(ype-ypb)
458
459 do m=2,i,2
460 radius=real(m)
461 call line(radius,ypb,radius,value)
462 end do
463
464 call gsplci(blue)
465 call gspmci(blue)
466 call gstxci(blue)
467
468 do j=jb,je,2
469 do i=ib,ie,2
470 do n=j,je,2
471 do m=i,ie,2
472 radius=sqrt(real((m-i)*(m-i)+(n-j)*(n-j)))
473
474 do nn=1,num
475 value=plt(m,n,nn)*plt(i,j,nn)
476
477 call pwritx(radius,value,symbol,1,1,0,0)
478 end do
479 end do
480 end do
481 end do
482 end do
483
484 end subroutine da_point_plot
485
486 subroutine da_line_plot(plt,num,nx,ny,nt,ib,jb,ie,je, &
487 xpb,xpe,ypb,ype)
488
489 implicit none
490
491 integer, intent(in) :: num,nx,ny,nt, &
492 ib,jb,ie,je
493 real, dimension(nx,ny,nt), intent(in) :: plt
494 real, intent(in) :: xpb,xpe
495 real, intent(out) :: ypb,ype
496
497 real, dimension(nx+ny) :: avg, sum
498
499 real :: radius, value
500
501 integer :: i, j, m, n, mm, nn
502
503 ! JRB hack
504 real :: xfb,xfe,yfb,yfe,blue
505 integer :: plot_style
506 ! JRB
507
508
509 sum = 0.0
510 avg = 0.0
511
512
513 do j=jb,je
514 do i=ib,ie
515 do n=j,je
516 do m=i,ie
517 radius=sqrt(real((m-i)*(m-i)+(n-j)*(n-j)))
518 mm=int(radius+0.5) + 1
519
520 do nn=1,num
521 value=plt(m,n,nn)*plt(i,j,nn)
522
523 avg(mm)=avg(mm)+value
524 sum(mm)=sum(mm)+1.0
525 end do
526 end do
527 end do
528 end do
529 end do
530
531 n = 0
532
533 do i=1,nx+ny
534 if (sum(i) < 0.5) exit
535
536 avg(i)=avg(i)/sum(i)
537 n=i
538 end do
539
540 ypb=minval(avg)*1.25
541 ype=maxval(avg)*1.25
542
543 call set(xfb,xfe,yfb,yfe,xpb,xpe,ypb,ype,plot_style)
544
545 call line(xpb,ypb,xpe,ypb)
546
547 call line(xpb,ypb,xpb,ype)
548
549 i=int(xpe) + 1
550 j=2
551
552 m=int(ype-ypb) + 1
553 n=m/10
554
555 value = ypb+0.02*(ype-ypb)
556
557 do m=2,i,2
558 radius=real(m)
559 call line(radius,ypb,radius,value)
560 end do
561
562 call gsplci(blue)
563 call gspmci(blue)
564 call gstxci(blue)
565
566 do i=2,nx+ny
567 if (sum(i) < 0.5) exit
568
569 call line(real(i-2), avg(i-1), real(i-1), avg(i))
570 end do
571
572 end subroutine da_line_plot
573
574 subroutine da_plot_sl(yr,r,nm,nn,slnt,cnst,k,varname)
575
576 implicit none
577
578 integer, intent(in) :: nm,nn,k
579 real(kind=8), dimension(0:nn), intent(in) :: yr,r
580 real(kind=8), intent(in) :: slnt,cnst
581 character(len=7), intent(in) :: varname
582
583 real :: x,y,xpb,xpe,ypb,ype
584
585 integer :: i
586
587 character(len=1), parameter :: symbol='.'
588 character(len=9) :: label
589
590 ! JRB hack
591 real :: xfb,xfe,yfb,yfe,red,xwb,xwe,ywb,ywe,xlb,ylb,blue
592 integer :: plot_style
593 ! JRB
594
595 call set(xwb,xwe,ywb,ywe,xwb,xwe,ywb,ywe,plot_style)
596
597 call gsplci(red)
598 call gspmci(red)
599 call gstxci(red)
600
601 write(label(1:9), fmt='(2a, i3)') varname(1:3), ' M=', k
602
603 call pwritx(xlb,ylb,label,9,1,0,0)
604
605 xpb=0.0
606 xpe=r(nm)*1.05
607
608 write(unit=*, fmt='(a, i5, f18.8)') &
609 'nm,r(nm)=', nm,r(nm)
610
611 ype=maxval(yr(0:nm))*1.05
612 ypb=0.0
613
614
615 write(unit=*, fmt='(a, 2f18.8)') &
616 'xpe,ype=',xpe,ype
617
618 call set(xfb,xfe,yfb,yfe,xpb,xpe,ypb,ype,plot_style)
619
620 call line(xpb,ypb,xpe,ypb)
621
622 call line(xpb,ypb,xpb,ype)
623
624 y = ypb+0.02*(ype-ypb)
625
626 do i=2,nm,2
627 x=real(i)
628 call line(x,ypb,x,y)
629 end do
630
631 x = xpb+0.02*(xpe-xpb)
632
633 y=0.0
634
635 do
636 y=y+1.0
637 if (y > ype) exit
638 call line(xpb,y,x,y)
639 end do
640
641 call gsplci(blue)
642 call gspmci(blue)
643 call gstxci(blue)
644
645 do i=1,nm
646 x= r(i)
647 y=yr(i)
648 call pwritx(x,y,symbol,1,1,0,0)
649
650 write(unit=*, fmt='(a,i3,2(f8.4,f18.8))') &
651 'i,x,y,r,yr=',i,x,y,r(i),yr(i)
652 end do
653
654 xpb=0.0
655 ypb=cnst
656
657 do i=1,nm
658 x=real(i)
659 y=slnt*x+cnst
660
661 if (y > ype) exit
662
663 call line(xpb,ypb,x,y)
664
665 xpb=x
666 ypb=y
667 end do
668
669 call frame
670
671 end subroutine da_plot_sl