309 subroutine print_ex_3d_arr(var_names,file_name_i,z,x,y,draw)
313 character(len=*),
intent(in) :: var_names(:)
314 character(len=*),
intent(in) :: file_name_i
315 real(dp),
intent(in) :: z(1:,1:,1:)
316 real(dp),
intent(in),
optional :: x(1:,1:,1:)
317 real(dp),
intent(in),
optional :: y(1:,1:,1:)
318 logical,
intent(in),
optional :: draw
322 integer :: iplt, ipntx, ipnty, nplt, npntx, npnty
323 real(dp),
allocatable :: x_fin(:,:,:)
324 real(dp),
allocatable :: y_fin(:,:,:)
326 character(len=max_str_ln) :: file_name
327 character(len=max_str_ln) :: write_format
334 if (
size(x,1).ne.
size(z,1) .or.
size(x,2).ne.
size(z,2) .or. &
335 &(
size(x,3).ne.1 .and.
size(x,3).ne.
size(z,3)))
then
336 write(*,*)
'shape(x)', shape(x)
337 write(*,*)
'shape(z)', shape(z)
338 call writo(
'In print_ex, the size of x and z has to be the &
339 &same... Skipping plot',persistent=.true.,warning=.true.)
344 if (
size(y,1).ne.
size(z,1) .or.
size(y,2).ne.
size(z,2) .or. &
345 &(
size(y,3).ne.1 .and.
size(y,3).ne.
size(z,3)))
then
346 write(*,*)
'shape(y)', shape(y)
347 write(*,*)
'shape(z)', shape(z)
348 call writo(
'In print_ex, the size of y and z has to be the &
349 &same... Skipping plot',persistent=.true.,warning=.true.)
355 allocate(x_fin(npntx,npnty,nplt))
356 if (
present (x))
then
357 if (
size(x,3).eq.1)
then
359 x_fin(:,:,iplt) = x(:,:,1)
366 x_fin(ipntx,:,:) = ipntx
370 allocate(y_fin(npntx,npnty,nplt))
371 if (
present (y))
then
372 if (
size(y,3).eq.1)
then
374 y_fin(:,:,iplt) = y(:,:,1)
381 y_fin(:,ipnty,:) = ipnty
386 if (trim(file_name_i).eq.
'')
then
387 file_name =
'temp_data_print_ex_3D_'//trim(
i2str(
rank))//
'_'//&
388 trim(
i2str(temp_id(1)))
389 temp_id(1) = temp_id(1)+1
391 file_name = trim(file_name_i)
400 write(file_i,
'(A)',iostat=istat) &
402 write_format=
"("//trim(
i2str(nplt*3))//
"(ES23.16,' '))"
405 write(file_i,fmt=trim(write_format),iostat=istat) &
406 &(x_fin(ipntx,ipnty,iplt), iplt = 1,nplt), &
407 &(y_fin(ipntx,ipnty,iplt), iplt = 1,nplt), &
408 &(z(ipntx,ipnty,iplt), iplt = 1,nplt)
410 write(file_i,
'(A)',iostat=istat)
''
412 write(file_i,
'(A)',iostat=istat)
''
421 if (
present(draw))
then
422 if (.not.draw)
return
424 call draw_ex(var_names,file_name,nplt,2,.true.)
426 if (trim(file_name_i).eq.
'')
then
427 call use_execute_command_line(
'rm '//
data_dir//
'/'//&
428 &trim(file_name)//
'.dat')
433 subroutine plot_hdf5_arr(var_names,file_name,vars,tot_dim,loc_offset,&
434 &X,Y,Z,col_id,col,sym_type,cont_plot,descr)
444 character(len=*),
intent(in) :: var_names(:)
445 character(len=*),
intent(in) :: file_name
446 real(dp),
intent(in),
target :: vars(:,:,:,:)
447 integer,
intent(in),
optional :: tot_dim(4)
448 integer,
intent(in),
optional :: loc_offset(4)
449 real(dp),
intent(in),
target,
optional :: X(:,:,:,:)
450 real(dp),
intent(in),
target,
optional :: Y(:,:,:,:)
451 real(dp),
intent(in),
target,
optional :: Z(:,:,:,:)
452 integer,
intent(in),
optional :: col_id
453 integer,
intent(in),
optional :: col
454 integer,
intent(in),
optional :: sym_type
455 logical,
intent(in),
optional :: cont_plot
456 character(len=*),
intent(in),
optional :: descr
461 integer :: col_id_loc
465 integer :: sym_type_loc
466 integer :: sym_pol, sym_tor
468 integer,
allocatable :: tot_sym_pol(:), tot_sym_tor(:)
469 integer :: tot_dim_loc(4)
470 integer :: loc_offset_loc(4)
471 integer :: tot_dim_3D(3)
472 integer :: loc_dim_3D(3)
473 integer :: loc_offset_3D(3)
481 logical :: col_mask(4)
483 logical :: cont_plot_loc
484 logical :: first_vc, last_vc
485 real(dp),
allocatable :: sym_ang(:,:,:)
486 real(dp) :: tol_sym = 1.e-8_dp
487 real(dp),
pointer :: var_3D(:,:,:) => null()
488 real(dp),
pointer :: X_3D(:,:,:) => null()
489 real(dp),
pointer :: Y_3D(:,:,:) => null()
490 real(dp),
pointer :: Z_3D(:,:,:) => null()
491 character(len=max_str_ln),
allocatable :: grd_names(:)
492 character(len=max_str_ln),
allocatable :: att_names(:)
496 if (
present(col_id)) col_id_loc = col_id
498 if (
present(col)) col_loc = col
501 n_plot =
size(vars,col_id_loc)
504 tot_dim_loc = shape(vars)
505 if (
present(tot_dim)) tot_dim_loc = tot_dim
507 if (
present(loc_offset)) loc_offset_loc = loc_offset
510 if (tot_dim_loc(col_id_loc).ne.n_plot)
then
512 call writo(
'In plot_HDF5, all the processes need to have the full &
513 &range in the dimension given by col_id',persistent=.true.,&
517 if (n_plot.eq.1 .and. col_loc.ne.1)
then
519 call writo(
'In plot_HDF5, if single plot, the collection type &
520 &needs to be one',persistent=.true.,warning=.true.)
526 col_mask(col_id_loc) = .true.
527 tot_dim_3d = pack(tot_dim_loc,.not.col_mask)
528 loc_dim_3d = pack(shape(vars),.not.col_mask)
529 loc_offset_3d = pack(loc_offset_loc,.not.col_mask)
532 if (tot_dim_3d(1).eq.loc_dim_3d(1) .and. &
533 &tot_dim_3d(2).eq.loc_dim_3d(2) .and. &
534 &tot_dim_3d(3).eq.loc_dim_3d(3))
then
541 cont_plot_loc = .false.
542 if (
present(cont_plot)) cont_plot_loc = cont_plot
550 if (
present(sym_type))
then
551 sym_type_loc = sym_type
552 else if (minval(tot_dim_3d).eq.1)
then
554 allocate(sym_ang(loc_dim_3d(1),loc_dim_3d(2),loc_dim_3d(3)))
561 call assign_pointers(id)
563 sym_ang = atan2(y_3d,x_3d)
564 if (maxval(sym_ang)-minval(sym_ang).lt.tol_sym .and. &
565 &(maxval(x_3d).ge.0._dp .neqv. &
566 &minval(x_3d).lt.0._dp).and.&
567 &(maxval(y_3d).ge.0._dp .neqv. &
568 &minval(y_3d).lt.0._dp)) &
571 sym_ang = atan2(sqrt(z_3d**2),sqrt(x_3d**2+y_3d**2))
572 if (maxval(sym_ang)-minval(sym_ang).lt.tol_sym) &
578 tot_sym_pol = sym_pol
579 tot_sym_tor = sym_tor
581 istat =
get_ser_var([sym_pol],tot_sym_pol,scatter=.true.)
583 istat =
get_ser_var([sym_tor],tot_sym_tor,scatter=.true.)
588 if (sum(tot_sym_pol).eq.
n_procs*n_plot)
then
590 else if (sum(tot_sym_tor).eq.
n_procs*n_plot)
then
598 if (col_loc.eq.4)
then
599 if (sym_type_loc.eq.1)
then
600 if (n_plot.ne.3)
then
602 call writo(
'For vector field plots, need 3 dimensions',&
607 if (n_plot.ne.2)
then
609 call writo(
'For symmetric vector field plots, need 2 &
610 &dimensions',warning=.true.)
617 if (.not.cont_plot_loc)
then
619 allocate(grd_names(n_plot))
620 allocate(att_names(n_plot))
621 if (col_loc.eq.1)
then
622 if (n_plot.eq.1)
then
623 if (
size(var_names).eq.n_plot)
then
624 att_names = var_names
625 else if (
size(var_names).gt.n_plot)
then
626 att_names = var_names(1:n_plot)
627 call writo(
'Too many variable names provided',&
628 &persistent=.true.,warning=.true.)
630 att_names(1:
size(var_names)) = var_names
631 do id =
size(var_names)+1,n_plot
632 att_names(id) =
'unnamed variable '//trim(
i2str(id))
634 call writo(
'Not enough variable names provided',&
635 &persistent=.true.,warning=.true.)
637 grd_names =
'default_grid_name'
639 if (
size(var_names).eq.n_plot)
then
640 grd_names = var_names
641 else if (
size(var_names).gt.n_plot)
then
642 grd_names = var_names(1:n_plot)
643 call writo(
'Too many variable names provided',&
644 &persistent=.true.,warning=.true.)
646 grd_names(1:
size(var_names)) = var_names
647 do id =
size(var_names)+1,n_plot
648 grd_names(id) =
'unnamed variable '//trim(
i2str(id))
650 call writo(
'Not enough variable names provided',&
651 &persistent=.true.,warning=.true.)
653 att_names =
'default_att_name'
656 att_names = var_names(1)
657 if (
size(var_names).gt.1)
call writo(
'For collections, only &
658 &the first variable name is used',persistent=.true.,&
660 grd_names =
'default_grid_name'
665 istat =
open_hdf5_file(file_info,file_name,sym_type=sym_type_loc,&
666 &descr=descr,ind_plot=ind_plot,cont_plot=cont_plot_loc)
670 if (.not.cont_plot_loc)
allocate(grids(n_plot))
673 if (sym_type_loc.eq.1)
then
683 first_vc = (id.eq.1 .or. col_loc.ne.4)
684 last_vc = (id.eq.n_plot .or. col_loc.ne.4)
687 if (.not.cont_plot_loc .and. last_vc)
then
688 if (sym_type_loc.eq.1)
then
696 call assign_pointers(id)
702 select case (sym_type_loc)
705 &
'X_'//trim(
i2str(id)),x_3d,tot_dim_3d,loc_dim_3d,&
706 &loc_offset_3d,ind_plot=ind_plot,&
707 &cont_plot=cont_plot_loc)
710 &
'Y_'//trim(
i2str(id)),y_3d,tot_dim_3d,loc_dim_3d,&
711 &loc_offset_3d,ind_plot=ind_plot,&
712 &cont_plot=cont_plot_loc)
715 &
'Z_'//trim(
i2str(id)),z_3d,tot_dim_3d,loc_dim_3d,&
716 &loc_offset_3d,ind_plot=ind_plot,&
717 &cont_plot=cont_plot_loc)
721 &
'R_'//trim(
i2str(id)),sqrt(x_3d**2+y_3d**2),&
722 &tot_dim_3d,loc_dim_3d,loc_offset_3d,&
723 &ind_plot=ind_plot,cont_plot=cont_plot_loc)
726 &
'Z_'//trim(
i2str(id)),z_3d,tot_dim_3d,loc_dim_3d,&
727 &loc_offset_3d,ind_plot=ind_plot,&
728 &cont_plot=cont_plot_loc)
732 &
'X_'//trim(
i2str(id)),x_3d,tot_dim_3d,loc_dim_3d,&
733 &loc_offset_3d,ind_plot=ind_plot,&
734 &cont_plot=cont_plot_loc)
737 &
'Y_'//trim(
i2str(id)),y_3d,tot_dim_3d,loc_dim_3d,&
738 &loc_offset_3d,ind_plot=ind_plot,&
739 &cont_plot=cont_plot_loc)
743 call writo(
'symmetry type '//&
744 &trim(
i2str(sym_type_loc))//
' not recognized',&
745 &persistent=.true.,warning=.true.)
752 if (.not.cont_plot_loc .and. first_vc)
then
753 if (sym_type_loc.eq.1)
then
763 if (col_loc.eq.4)
then
769 &trim(
i2str(id)),var_3d,tot_dim_3d,loc_dim_3d,&
770 &loc_offset_3d,ind_plot=ind_plot,cont_plot=cont_plot_loc)
774 if (.not.cont_plot_loc)
then
775 if (col_loc.eq.4)
then
778 &tot_dim_3d,reset=.true.,ind_plot=ind_plot,&
779 &cont_plot=cont_plot_loc)
781 &att_names(1),1,2,reset=.true.,ind_plot=ind_plot)
785 &att_names(id),1,1,reset=.true.,ind_plot=ind_plot)
792 if (.not.cont_plot_loc .and. last_vc)
then
793 if (col_loc.eq.2)
then
795 &grid_time=id*1._dp,grid_top=top,grid_geom=geom,&
796 &grid_atts=att,reset=.true.,ind_plot=ind_plot)
800 &grid_top=top,grid_geom=geom,grid_atts=att,&
801 &reset=.true.,ind_plot=ind_plot)
809 if (.not.cont_plot_loc)
then
810 if (col_loc.eq.1 .or. col_loc.eq.4)
then
813 if (col_loc.eq.4 .and. id.lt.n_plot) cycle
821 &grid_grids=grids,reset=.true.,ind_plot=ind_plot)
833 &cont_plot=cont_plot_loc)
838 nullify(x_3d,y_3d,z_3d)
840 if (.not.cont_plot_loc)
then
849 subroutine assign_pointers(id)
859 if (col_id_loc.eq.1)
then
860 if (
size(x,1).eq.1) id_loc = 1
861 x_3d => x(id_loc,:,:,:)
862 else if (col_id_loc.eq.2)
then
863 if (
size(x,2).eq.1) id_loc = 1
864 x_3d => x(:,id_loc,:,:)
865 else if (col_id_loc.eq.3)
then
866 if (
size(x,3).eq.1) id_loc = 1
867 x_3d => x(:,:,id_loc,:)
868 else if (col_id_loc.eq.4)
then
869 if (
size(x,4).eq.1) id_loc = 1
870 x_3d => x(:,:,:,id_loc)
876 allocate(x_3d(loc_dim_3d(1),loc_dim_3d(2),loc_dim_3d(3)))
877 do jd = 1,loc_dim_3d(1)
878 x_3d(jd,:,:) = loc_offset_3d(1) + jd - 1
884 if (col_id_loc.eq.1)
then
885 if (
size(y,1).eq.1) id_loc = 1
886 y_3d => y(id_loc,:,:,:)
887 else if (col_id_loc.eq.2)
then
888 if (
size(y,2).eq.1) id_loc = 1
889 y_3d => y(:,id_loc,:,:)
890 else if (col_id_loc.eq.3)
then
891 if (
size(y,3).eq.1) id_loc = 1
892 y_3d => y(:,:,id_loc,:)
893 else if (col_id_loc.eq.4)
then
894 if (
size(y,4).eq.1) id_loc = 1
895 y_3d => y(:,:,:,id_loc)
901 allocate(y_3d(loc_dim_3d(1),loc_dim_3d(2),loc_dim_3d(3)))
902 do jd = 1,loc_dim_3d(2)
903 y_3d(:,jd,:) = loc_offset_3d(2) + jd - 1
909 if (col_id_loc.eq.1)
then
910 if (
size(z,1).eq.1) id_loc = 1
911 z_3d => z(id_loc,:,:,:)
912 else if (col_id_loc.eq.2)
then
913 if (
size(z,2).eq.1) id_loc = 1
914 z_3d => z(:,id_loc,:,:)
915 else if (col_id_loc.eq.3)
then
916 if (
size(z,3).eq.1) id_loc = 1
917 z_3d => z(:,:,id_loc,:)
918 else if (col_id_loc.eq.4)
then
919 if (
size(z,4).eq.1) id_loc = 1
920 z_3d => z(:,:,:,id_loc)
926 allocate(z_3d(loc_dim_3d(1),loc_dim_3d(2),loc_dim_3d(3)))
927 do jd = 1,loc_dim_3d(3)
928 z_3d(:,:,jd) = loc_offset_3d(3) + jd - 1
932 if (col_id_loc.eq.1)
then
933 var_3d => vars(id,:,:,:)
934 else if (col_id_loc.eq.2)
then
935 var_3d => vars(:,id,:,:)
936 else if (col_id_loc.eq.3)
then
937 var_3d => vars(:,:,id,:)
938 else if (col_id_loc.eq.4)
then
939 var_3d => vars(:,:,:,id)
944 end subroutine assign_pointers
947 subroutine plot_hdf5_ind(var_name,file_name,var,tot_dim,loc_offset,&
948 &X,Y,Z,cont_plot,descr)
951 character(len=*),
intent(in) :: var_name
952 character(len=*),
intent(in) :: file_name
953 real(dp),
intent(in) :: var(:,:,:)
954 integer,
intent(in),
optional :: tot_dim(3)
955 integer,
intent(in),
optional :: loc_offset(3)
956 real(dp),
intent(in),
target,
optional :: X(:,:,:)
957 real(dp),
intent(in),
target,
optional :: Y(:,:,:)
958 real(dp),
intent(in),
target,
optional :: Z(:,:,:)
959 logical,
intent(in),
optional :: cont_plot
960 character(len=*),
intent(in),
optional :: descr
963 integer :: tot_dim_loc(4), loc_offset_loc(4)
966 tot_dim_loc = [shape(var),1]
967 if (
present(tot_dim)) tot_dim_loc = [tot_dim,1]
969 if (
present(loc_offset)) loc_offset_loc = [loc_offset,0]
974 call plot_hdf5_arr([var_name],file_name,&
975 &reshape(var,[
size(var,1),
size(var,2),
size(var,3),1]),&
976 &tot_dim_loc,loc_offset_loc,&
977 &x=reshape(x,[
size(x,1),
size(x,2),
size(x,3),1]),&
978 &y=reshape(y,[
size(y,1),
size(y,2),
size(y,3),1]),&
979 &z=reshape(z,[
size(z,1),
size(z,2),
size(z,3),1]),&
980 &col=1,cont_plot=cont_plot,descr=descr)
982 call plot_hdf5_arr([var_name],file_name,&
983 &reshape(var,[
size(var,1),
size(var,2),
size(var,3),1]),&
984 &tot_dim_loc,loc_offset_loc,&
985 &x=reshape(x,[
size(x,1),
size(x,2),
size(x,3),1]),&
986 &y=reshape(y,[
size(y,1),
size(y,2),
size(y,3),1]),&
987 &col=1,cont_plot=cont_plot,descr=descr)
991 call plot_hdf5_arr([var_name],file_name,&
992 &reshape(var,[
size(var,1),
size(var,2),
size(var,3),1]),&
993 &tot_dim_loc,loc_offset_loc,&
994 &x=reshape(x,[
size(x,1),
size(x,2),
size(x,3),1]),&
995 &z=reshape(z,[
size(z,1),
size(z,2),
size(z,3),1]),&
996 &col=1,cont_plot=cont_plot,descr=descr)
998 call plot_hdf5_arr([var_name],file_name,&
999 &reshape(var,[
size(var,1),
size(var,2),
size(var,3),1]),&
1000 &tot_dim_loc,loc_offset_loc,&
1001 &x=reshape(x,[
size(x,1),
size(x,2),
size(x,3),1]),&
1002 &col=1,cont_plot=cont_plot,descr=descr)
1006 if (
present(y))
then
1007 if (
present(z))
then
1008 call plot_hdf5_arr([var_name],file_name,&
1009 &reshape(var,[
size(var,1),
size(var,2),
size(var,3),1]),&
1010 &tot_dim_loc,loc_offset_loc,&
1011 &y=reshape(y,[
size(y,1),
size(y,2),
size(y,3),1]),&
1012 &z=reshape(z,[
size(z,1),
size(z,2),
size(z,3),1]),&
1013 &col=1,cont_plot=cont_plot,descr=descr)
1015 call plot_hdf5_arr([var_name],file_name,&
1016 &reshape(var,[
size(var,1),
size(var,2),
size(var,3),1]),&
1017 &tot_dim_loc,loc_offset_loc,&
1018 &y=reshape(y,[
size(y,1),
size(y,2),
size(y,3),1]),&
1019 &col=1,cont_plot=cont_plot,descr=descr)
1022 if (
present(z))
then
1023 call plot_hdf5_arr([var_name],file_name,&
1024 &reshape(var,[
size(var,1),
size(var,2),
size(var,3),1]),&
1025 &tot_dim_loc,loc_offset_loc,&
1026 &col=1,z=reshape(z,[
size(z,1),
size(z,2),
size(z,3),1]),&
1027 &cont_plot=cont_plot,descr=descr)
1029 call plot_hdf5_arr([var_name],file_name,&
1030 &reshape(var,[
size(var,1),
size(var,2),
size(var,3),1]),&
1031 &tot_dim_loc,loc_offset_loc,&
1032 &col=1,cont_plot=cont_plot,descr=descr)
1076 subroutine draw_ex(var_names,draw_name,nplt,draw_dim,plot_on_screen,&
1077 &ex_plot_style,data_name,draw_ops,extra_ops,is_animated,ranges,delay,&
1083 character(len=*),
intent(in) :: var_names(:)
1084 character(len=*),
intent(in) :: draw_name
1085 integer,
intent(in) :: nplt
1086 integer,
intent(in) :: draw_dim
1087 logical,
intent(in) :: plot_on_screen
1089 character(len=*),
intent(in),
optional :: data_name
1090 character(len=*),
intent(in),
optional :: draw_ops(:)
1091 character(len=*),
intent(in),
optional :: extra_ops
1092 logical,
intent(in),
optional :: is_animated
1093 real(
dp),
intent(in),
optional :: ranges(:,:)
1094 integer,
intent(in),
optional :: delay
1095 logical,
intent(in),
optional :: persistent
1098 character(len=4) :: ex_ext
1099 character(len=5*max_str_ln) :: cmdmsg
1100 character(len=max_str_ln) :: script_name
1101 character(len=max_str_ln) :: ex_prog_name
1102 character(len=max_str_ln) :: data_name_loc
1103 character(len=max_str_ln),
allocatable :: var_names_loc(:)
1104 logical :: is_animated_loc
1105 logical :: run_shell_necessary
1106 logical :: persistent_loc
1107 real(
dp),
allocatable :: ranges_loc(:,:)
1108 integer :: delay_loc
1109 integer :: n_draw_ops
1114 integer :: ex_plot_style_loc
1120 data_name_loc = draw_name
1121 if (
present(data_name)) data_name_loc = data_name
1122 allocate(var_names_loc(nplt))
1123 if (
size(var_names).eq.nplt)
then
1124 var_names_loc = var_names
1127 var_names_loc(iplt) = trim(var_names(1))//
' ('//&
1128 &trim(
i2str(iplt))//
' / '//trim(
i2str(nplt))//
')'
1131 is_animated_loc = .false.
1132 if (
present(is_animated)) is_animated_loc = is_animated
1133 ex_plot_style_loc = ex_plot_style_global
1135 persistent_loc = .false.
1136 if (
present(persistent)) persistent_loc = persistent
1139 run_shell_necessary = .true.
1142 if (is_animated_loc)
then
1143 if (draw_dim.ne.1)
then
1144 call writo(
'Animations are only possible in 2D',warning=.true.,&
1149 if (
present(delay))
then
1152 delay_loc = max(250/nplt,10)
1158 if (plot_on_screen .and. &
1159 &(ex_plot_style_loc.ne.1 .and. is_animated_loc))
then
1160 script_name = trim(
script_dir)//
'/'//
'temp_script_draw_ex_'//&
1162 temp_id(2) = temp_id(2)+1
1164 script_name = trim(
script_dir)//
'/'//trim(draw_name)
1166 select case (ex_plot_style_loc)
1168 script_name = trim(script_name)//
'.gnu'
1170 ex_prog_name =
'gnuplot'
1172 script_name = trim(script_name)//
'.py'
1173 if (draw_dim.eq.1)
then
1178 ex_prog_name =
'python'
1182 open(
nextunit(cmd_i),file=trim(script_name),iostat=istat)
1186 if (
present(draw_ops))
then
1187 n_draw_ops =
size(draw_ops)
1192 select case (ex_plot_style_loc)
1194 if (is_animated_loc)
then
1195 call draw_ex_animated_gnuplot()
1197 call draw_ex_gnuplot()
1200 select case (draw_dim)
1203 if (is_animated_loc)
then
1204 call draw_ex_animated_bokeh()
1206 call draw_ex_bokeh()
1209 call draw_ex_mayavi()
1219 if (run_shell_necessary) &
1220 &
call use_execute_command_line(trim(ex_prog_name)//
' "'//&
1221 &trim(script_name)//
'"'//err_output_str,exitstat=istat,&
1222 &cmdstat=cmdstat,cmdmsg=cmdmsg)
1224 if (istat.ne.0)
then
1225 call writo(
'Failed to plot '//trim(draw_name)//
'.'//&
1226 &trim(ex_ext),persistent=.true.,warning=.true.)
1228 if (cmdstat.ne.0)
then
1229 call writo(
'Failed to plot '//trim(draw_name)//
'.'//&
1230 &trim(ex_ext),persistent=.true.,warning=.true.)
1232 call writo(
'System message: "'//trim(cmdmsg)//
'"',&
1234 if (.not.plot_on_screen)
call writo(&
1235 &
'Try running "'//trim(ex_prog_name)//
' "'//&
1236 &trim(script_name)//
'"'//
'" manually',persistent=.true.)
1239 if (plot_on_screen .and. run_shell_necessary .and. &
1240 &.not.persistent_loc)
then
1241 call use_execute_command_line(
'rm '//&
1242 &trim(script_name),exitstat=istat,&
1243 &cmdstat=cmdstat,cmdmsg=cmdmsg)
1246 call writo(
'Plot in output file "'//&
1247 &trim(
plot_dir)//
'/'//trim(draw_name)//&
1248 &
'.'//trim(ex_ext)//
'"',persistent=.true.)
1254 subroutine draw_ex_gnuplot
1256 write(cmd_i,
"(A)",iostat=istat)
'set grid'
1257 write(cmd_i,
"(A)",iostat=istat)
'set border 4095 front linetype -1 &
1259 if (plot_on_screen)
then
1260 if (persistent_loc)
then
1261 write(cmd_i,
"(A)",iostat=istat)
'set terminal wxt persist'
1263 write(cmd_i,
"(A)",iostat=istat)
'set terminal wxt'
1266 write(cmd_i,
"(A)",iostat=istat)
'set terminal pdf size '//&
1268 write(cmd_i,
"(A)",iostat=istat)
'set output "'//&
1269 &trim(
plot_dir)//
'/'//trim(draw_name)//
'.pdf"'
1273 if (nplt.gt.10)
write(cmd_i,
"(A)",iostat=istat)
'set nokey;'
1276 if (
present(extra_ops))
write(cmd_i,
"(A)",iostat=istat) &
1280 if (n_draw_ops.eq.0)
then
1281 do iplt = 1,min(nplt,
size(line_clrs))
1282 write(cmd_i,
"(A)",iostat=istat)
'set style line '//&
1283 &trim(
i2str(iplt))//
' lc rgb '//&
1284 &line_clrs(iplt)//
' '//trim(line_style)
1289 select case (draw_dim)
1291 write(cmd_i,
"(A)",iostat=istat)
'plot \'
1293 write(cmd_i,
"(A)",iostat=istat,advance=
"no") &
1294 &
' "'//trim(
data_dir)//
'/'//trim(data_name_loc)//&
1295 &
'.dat" using '//trim(
i2str(iplt))//
':'//&
1296 &trim(
i2str(nplt+iplt))//
' title "'//&
1297 &trim(var_names_loc(iplt))//
'" '//&
1298 &trim(loc_draw_op())
1299 if (iplt.eq.nplt)
then
1300 write(cmd_i,
"(A)",iostat=istat)
''
1302 write(cmd_i,
"(A)",iostat=istat)
', \'
1306 write(cmd_i,
"(A)",iostat=istat)
'splot \'
1308 write(cmd_i,
"(A)",iostat=istat,advance=
"no") &
1309 &
' "'//trim(
data_dir)//
'/'//trim(data_name_loc)//&
1310 &
'.dat" using '//trim(
i2str(iplt))//
':'//&
1311 &trim(
i2str(nplt+iplt))//
':'//&
1312 &trim(
i2str(2*nplt+iplt))//
' title "'//&
1313 &trim(var_names_loc(iplt))//&
1314 &
'" '//trim(loc_draw_op())
1315 if (iplt.eq.nplt)
then
1316 write(cmd_i,
"(A)",iostat=istat)
''
1318 write(cmd_i,
"(A)",iostat=istat)
', \'
1322 write(cmd_i,
"(A)",iostat=istat)
'splot \'
1324 write(cmd_i,
"(A)",iostat=istat,advance=
"no") &
1325 &
' "'//trim(
data_dir)//
'/'//trim(data_name_loc)//&
1326 &
'.dat" using ('//trim(
i2str(iplt))//
'):'//&
1327 &trim(
i2str(iplt))//
':'//trim(
i2str(nplt+iplt))//&
1328 &
' title "'//trim(var_names_loc(iplt))//
'" '//&
1329 &trim(loc_draw_op())
1330 if (iplt.eq.nplt)
then
1331 write(cmd_i,
"(A)",iostat=istat)
''
1333 write(cmd_i,
"(A)",iostat=istat)
', \'
1337 call writo(
'No draw_dim associated with '//&
1338 &trim(
i2str(draw_dim)),persistent=.true.)
1344 write(cmd_i,
"(A)",iostat=istat)
''
1345 if (plot_on_screen)
write(cmd_i,
"(A)",iostat=istat)
'pause -1'
1346 end subroutine draw_ex_gnuplot
1352 subroutine draw_ex_bokeh
1354 write(cmd_i,
"(A)",iostat=istat)
'from numpy import genfromtxt'
1355 write(cmd_i,
"(A)",iostat=istat)
'from bokeh.plotting import &
1356 &figure, output_file, save, show'
1357 write(cmd_i,
"(A)",iostat=istat)
'from bokeh.layouts import &
1359 write(cmd_i,
"(A)",iostat=istat)
'from bokeh.models.widgets import &
1361 write(cmd_i,
"(A)",iostat=istat)
'from bokeh.models import CustomJS'
1362 write(cmd_i,
"(A)",iostat=istat)
''
1363 write(cmd_i,
"(A)",iostat=istat)
'TOOLS="hover,pan,wheel_zoom,&
1364 &box_zoom,reset,save"'
1365 write(cmd_i,
"(A)",iostat=istat)
''
1366 write(cmd_i,
"(A)",iostat=istat)
'output_file("'//trim(
plot_dir)//&
1367 &
'/'//trim(draw_name)//
'.html", title="'//trim(draw_name)//&
1369 write(cmd_i,
"(A)",iostat=istat)
''
1370 write(cmd_i,
"(A)",iostat=istat)
'data = genfromtxt("'//&
1371 &trim(
data_dir)//
'/'//trim(data_name_loc)//
'.dat")'
1372 write(cmd_i,
"(A)",iostat=istat)
''
1373 write(cmd_i,
"(A)",iostat=istat)
'p = figure(toolbar_location=&
1374 &"above",tools=TOOLS,active_drag="box_zoom",active_scroll=&
1375 &"wheel_zoom",plot_width=600, plot_height=600)'
1376 write(cmd_i,
"(A)",iostat=istat)
''
1379 if (
present(extra_ops))
write(cmd_i,
"(A)",iostat=istat) &
1385 write(cmd_i,
"(A)",iostat=istat)
'l'//trim(
i2str(iplt))//&
1386 &
' = p.line(data[:,'//trim(
i2str(iplt-1))//
'],data[:,'//&
1387 &trim(
i2str(nplt+iplt-1))//
'],'//trim(loc_draw_op(1))//
')'
1389 write(cmd_i,
"(A)",iostat=istat)
'c'//trim(
i2str(iplt))//&
1390 &
' = p.circle(data[:,'//&
1391 &trim(
i2str(iplt-1))//
'],data[:,'//&
1392 &trim(
i2str(nplt+iplt-1))//
'],'//trim(loc_draw_op(2))//
')'
1394 write(cmd_i,
"(A)",iostat=istat)
''
1397 write(cmd_i,
"(A)",iostat=istat)
'labels=list()'
1399 write(cmd_i,
"(A)",iostat=istat)
'labels.append("'//&
1400 &trim(var_names_loc(iplt))//
'")'
1402 write(cmd_i,
"(A)",iostat=istat)
'active=list()'
1404 write(cmd_i,
"(A)",iostat=istat)
'active.append('//&
1405 &trim(
i2str(iplt-1))//
')'
1407 write(cmd_i,
"(A)",iostat=istat)
'checkbox = CheckboxGroup(&
1408 &labels=labels, active=active)'
1409 write(cmd_i,
"(A)",iostat=istat)
''
1412 write(cmd_i,
"(A)",iostat=istat)
'args=dict()'
1414 write(cmd_i,
"(A)",iostat=istat)
'args["l'//&
1415 &trim(
i2str(iplt))//
'"] = l'//trim(
i2str(iplt))
1416 write(cmd_i,
"(A)",iostat=istat)
'args["c'//&
1417 &trim(
i2str(iplt))//
'"] = c'//trim(
i2str(iplt))
1420 write(cmd_i,
"(A)",iostat=istat)
'code="""'
1421 write(cmd_i,
"(A)",iostat=istat)
'//console.log(cb_obj.active);'
1423 write(cmd_i,
"(A)",iostat=istat)
'l'//trim(
i2str(iplt))//&
1424 &
'.visible = false;'
1425 write(cmd_i,
"(A)",iostat=istat)
'c'//trim(
i2str(iplt))//&
1426 &
'.visible = false;'
1428 write(cmd_i,
"(A)",iostat=istat)
'for (i in cb_obj.active) {'
1429 write(cmd_i,
"(A)",iostat=istat) &
1430 &
' //console.log(cb_obj.active[i]);'
1432 write(cmd_i,
"(A)",iostat=istat)
' if (cb_obj.active[i] == '&
1433 &//trim(
i2str(iplt-1))//
') {'
1434 write(cmd_i,
"(A)",iostat=istat)
' l'//trim(
i2str(iplt))//&
1436 write(cmd_i,
"(A)",iostat=istat)
' c'//trim(
i2str(iplt))//&
1438 write(cmd_i,
"(A)",iostat=istat)
' }'
1440 write(cmd_i,
"(A)",iostat=istat)
'}'
1441 write(cmd_i,
"(A)",iostat=istat)
'"""'
1443 write(cmd_i,
"(A)",iostat=istat)
'checkbox.callback = CustomJS(&
1444 &args=args, code=code)'
1447 write(cmd_i,
"(A)",iostat=istat) &
1448 &
'layout = row(p, widgetbox(checkbox))'
1449 write(cmd_i,
"(A)",iostat=istat)
''
1452 write(cmd_i,
"(A)",iostat=istat)
''
1453 if (plot_on_screen)
then
1454 write(cmd_i,
"(A)",iostat=istat)
'show(layout)'
1456 write(cmd_i,
"(A)",iostat=istat)
'save(layout)'
1458 end subroutine draw_ex_bokeh
1461 subroutine draw_ex_mayavi
1463 write(cmd_i,
"(A)",iostat=istat)
'from numpy import genfromtxt, &
1464 &array, size, zeros, amin, amax'
1465 write(cmd_i,
"(A)",iostat=istat)
'from mayavi.mlab import outline, &
1466 &mesh, points3d, savefig, colorbar, axes, savefig, show'
1467 write(cmd_i,
"(A)",iostat=istat)
''
1468 write(cmd_i,
"(A)",iostat=istat)
'nplt = '//trim(
i2str(nplt))
1471 write(cmd_i,
"(A)",iostat=istat)
'npnty = -1'
1472 write(cmd_i,
"(A)",iostat=istat)
'with open("'//trim(
data_dir)//&
1473 &
'/'//trim(data_name_loc)//
'.dat") as f:'
1474 write(cmd_i,
"(A)",iostat=istat)
' for line in f:'
1475 write(cmd_i,
"(A)",iostat=istat)
' if (npnty < 0):'
1476 write(cmd_i,
"(A)",iostat=istat)
' if &
1477 &(line.strip()[0] == "#"):'
1478 write(cmd_i,
"(A)",iostat=istat)
' npnty = 0'
1479 write(cmd_i,
"(A)",iostat=istat)
' else:'
1480 write(cmd_i,
"(A)",iostat=istat)
' if (not line.strip()):'
1481 write(cmd_i,
"(A)",iostat=istat)
' break'
1482 write(cmd_i,
"(A)",iostat=istat)
' else:'
1483 write(cmd_i,
"(A)",iostat=istat)
' npnty += 1'
1486 write(cmd_i,
"(A)",iostat=istat)
'data = genfromtxt("'//&
1487 &trim(
data_dir)//
'/'//trim(data_name_loc)//
'.dat")'
1488 write(cmd_i,
"(A)",iostat=istat)
'dims = &
1489 &array([size(data,0)/npnty,npnty,nplt])'
1490 write(cmd_i,
"(A)",iostat=istat)
''
1491 write(cmd_i,
"(A)",iostat=istat)
'X = zeros(dims)'
1492 write(cmd_i,
"(A)",iostat=istat)
'Y = zeros(dims)'
1493 write(cmd_i,
"(A)",iostat=istat)
'Z = zeros(dims)'
1494 write(cmd_i,
"(A)",iostat=istat)
'for i in range(0,dims[0]):'
1495 write(cmd_i,
"(A)",iostat=istat)
' for j in range(0,dims[2]):'
1496 write(cmd_i,
"(A)",iostat=istat)
' X[i,:,j] = &
1497 &data[npnty*i:npnty*(i+1),j]'
1498 write(cmd_i,
"(A)",iostat=istat)
' Y[i,:,j] = &
1499 &data[npnty*i:npnty*(i+1),nplt+j]'
1500 write(cmd_i,
"(A)",iostat=istat)
' Z[i,:,j] = &
1501 &data[npnty*i:npnty*(i+1),nplt*2+j]'
1502 write(cmd_i,
"(A)",iostat=istat)
''
1503 write(cmd_i,
"(A)",iostat=istat)
'minZ = amin(Z)'
1504 write(cmd_i,
"(A)",iostat=istat)
'maxZ = amax(Z)'
1505 write(cmd_i,
"(A)",iostat=istat)
''
1508 if (
present(extra_ops))
write(cmd_i,
"(A)",iostat=istat) &
1512 select case (draw_dim)
1514 write(cmd_i,
"(A)",iostat=istat)
'for j in range(0,dims[2]):'
1515 write(cmd_i,
"(A)",iostat=istat)
' &
1516 &mesh(X[:,:,j],Y[:,:,j],Z[:,:,j],vmin=minZ,vmax=maxZ)'
1518 write(cmd_i,
"(A)",iostat=istat)
'for j in range(0,dims[2]):'
1519 write(cmd_i,
"(A)",iostat=istat)
' &
1520 &points3d(X[:,:,j],Y[:,:,j],Z[:,:,j],mode="sphere",&
1521 &scale_factor=0.05,vmin=minZ,vmax=maxZ)'
1523 write(cmd_i,
"(A)",iostat=istat)
'outline()'
1524 write(cmd_i,
"(A)",iostat=istat)
'colorbar()'
1525 write(cmd_i,
"(A)",iostat=istat)
'axes()'
1528 if (plot_on_screen)
then
1529 write(cmd_i,
"(A)",iostat=istat)
''
1530 write(cmd_i,
"(A)",iostat=istat)
'show()'
1532 write(cmd_i,
"(A)",iostat=istat)
''
1533 write(cmd_i,
"(A)",iostat=istat)
'#show()'
1534 write(cmd_i,
"(A)",iostat=istat)
'savefig("'//trim(
plot_dir)//&
1535 &
'/'//trim(draw_name)//
'.png"'//
',magnification=4)'
1537 end subroutine draw_ex_mayavi
1540 subroutine draw_ex_animated_gnuplot
1542 write(cmd_i,
"(A)",iostat=istat)
'set grid'
1543 write(cmd_i,
"(A)",iostat=istat)
'set border 4095 front linetype -1 &
1545 write(cmd_i,
"(A)",iostat=istat)
'set terminal gif animate &
1546 &delay '//trim(
i2str(delay_loc))//
' size '//&
1549 write(cmd_i,
"(A)",iostat=istat)
'set output "'//trim(
plot_dir)//&
1550 &
'/'//trim(draw_name)//
'.gif"'
1553 allocate(ranges_loc(2,2))
1554 if (
present(ranges))
then
1555 if (
size(ranges,1).eq.2 .and.
size(ranges,2).eq.2)
then
1558 call writo(
'invalid ranges given &
1559 &to draw_ex_animated',persistent=.true.,&
1564 ranges_loc(:,1) = huge(1._dp)
1565 ranges_loc(:,2) = -huge(1._dp)
1567 call get_ranges(ranges_loc)
1571 write(cmd_i,
"(A)",iostat=istat)
'set xrange ['//&
1572 &trim(
r2str(ranges_loc(1,1)))//
':'//&
1573 &trim(
r2str(ranges_loc(1,2)))//
'];'
1574 write(cmd_i,
"(A)",iostat=istat)
'set yrange ['//&
1575 &trim(
r2str(ranges_loc(2,1)))//
':'//&
1576 &trim(
r2str(ranges_loc(2,2)))//
'];'
1579 if (
present(extra_ops))
write(cmd_i,
"(A)",iostat=istat) &
1583 if (n_draw_ops.eq.0)
then
1584 do iplt = 1,min(nplt,
size(line_clrs))
1585 write(cmd_i,
"(A)",iostat=istat)
'set style line '//&
1586 &trim(
i2str(iplt))//
' lc rgb '//&
1587 &line_clrs(iplt)//
' '//trim(line_style)
1593 write(cmd_i,
"(A)",iostat=istat)
'plot "'//trim(
data_dir)//
'/'//&
1594 &trim(data_name_loc)//
'.dat" using '//trim(
i2str(iplt))//&
1595 &
':'//trim(
i2str(nplt+iplt))//
' title "'//&
1596 &trim(var_names_loc(iplt))//
'" '//&
1597 &trim(loc_draw_op())
1601 write(cmd_i,
"(A)",iostat=istat)
''
1602 end subroutine draw_ex_animated_gnuplot
1605 subroutine draw_ex_animated_bokeh
1607 write(cmd_i,
"(A)",iostat=istat)
'# Note: it is necessary to first &
1608 &run "bokeh serve" before calling python'
1609 write(cmd_i,
"(A)",iostat=istat)
''
1610 write(cmd_i,
"(A)",iostat=istat)
'from numpy import genfromtxt'
1611 write(cmd_i,
"(A)",iostat=istat)
'from bokeh.plotting import &
1612 &figure, output_file, save, show, curdoc'
1613 write(cmd_i,
"(A)",iostat=istat)
'from bokeh.client import &
1615 write(cmd_i,
"(A)",iostat=istat)
'from bokeh.driving import repeat'
1616 write(cmd_i,
"(A)",iostat=istat)
''
1617 write(cmd_i,
"(A)",iostat=istat)
'TOOLS="crosshair,pan,wheel_zoom,&
1618 &box_zoom,reset,save"'
1619 write(cmd_i,
"(A)",iostat=istat)
''
1620 write(cmd_i,
"(A)",iostat=istat)
'output_file("'//trim(
plot_dir)//&
1621 &
'/'//trim(draw_name)//
'.html", title="'//trim(draw_name)//&
1623 write(cmd_i,
"(A)",iostat=istat)
''
1624 write(cmd_i,
"(A)",iostat=istat)
'data = genfromtxt("'//&
1625 &trim(
data_dir)//
'/'//trim(data_name_loc)//
'.dat")'
1626 write(cmd_i,
"(A)",iostat=istat)
''
1627 write(cmd_i,
"(A)",iostat=istat)
'p = figure(toolbar_location=&
1628 &"above",'//
'tools=TOOLS,active_drag="pan",active_scroll=&
1629 &"wheel_zoom",'//
'plot_width=600, plot_height=600)'
1630 write(cmd_i,
"(A)",iostat=istat)
''
1633 if (
present(extra_ops))
write(cmd_i,
"(A)",iostat=istat) &
1637 write(cmd_i,
"(A)",iostat=istat)
'line = p.line(data[:,0],data[:,'//&
1638 &trim(
i2str(nplt))//
'],'//trim(loc_draw_op(1))//
',legend="'//&
1639 &trim(var_names_loc(1))//
'")'
1641 write(cmd_i,
"(A)",iostat=istat)
'circle = p.circle(data[:,'//&
1642 &trim(
i2str(0))//
'],data[:,'//trim(
i2str(nplt))//
'],'//&
1643 &trim(loc_draw_op(2))//
')'
1646 write(cmd_i,
"(A)",iostat=istat)
''
1647 write(cmd_i,
"(A)",iostat=istat)
'session = push_session(curdoc())'
1648 write(cmd_i,
"(A)",iostat=istat)
'session.show(p)'
1649 write(cmd_i,
"(A)",iostat=istat)
''
1650 write(cmd_i,
"(A)",iostat=istat)
'seq = range(0,'//&
1651 &trim(
i2str(nplt))//
')'
1652 write(cmd_i,
"(A)",iostat=istat)
'@repeat(seq)'
1653 write(cmd_i,
"(A)",iostat=istat)
'def update(i):'
1654 write(cmd_i,
"(A)",iostat=istat)
' line.data_source.data["x"] = &
1656 write(cmd_i,
"(A)",iostat=istat)
' line.data_source.data["y"] = &
1657 &data[:,'//trim(
i2str(nplt))//
'+i]'
1658 write(cmd_i,
"(A)",iostat=istat)
' circle.data_source.data["x"] = &
1660 write(cmd_i,
"(A)",iostat=istat)
' circle.data_source.data["y"] = &
1662 &trim(
i2str(nplt))//
'+i]'
1663 write(cmd_i,
"(A)",iostat=istat)
''
1664 write(cmd_i,
"(A)",iostat=istat)
'curdoc().add_periodic_callback(&
1665 &update, '//trim(
i2str(delay_loc*10))//
')'
1666 write(cmd_i,
"(A)",iostat=istat)
'session.loop_until_closed()'
1667 end subroutine draw_ex_animated_bokeh
1670 subroutine get_ranges(ranges)
1672 real(
dp),
intent(inout) :: ranges(:,:)
1676 real(
dp),
allocatable :: loc_data(:)
1677 character(len=1) :: loc_data_char
1682 &
'.dat',status=
'old',iostat=istat)
1684 if (istat.eq.0)
then
1686 allocate(loc_data(2*nplt))
1690 do while (istat.eq.0)
1691 read(data_i,*,iostat=istat) loc_data_char
1692 if (istat.eq.0)
then
1693 if (loc_data_char.ne.
'#')
then
1694 backspace(unit=data_i)
1695 read(data_i,*,iostat=istat) loc_data
1697 &min(ranges(1,1),minval(loc_data(1:nplt)))
1699 &max(ranges(1,2),maxval(loc_data(1:nplt)))
1702 &minval(loc_data(nplt+1:2*nplt)))
1705 &maxval(loc_data(nplt+1:2*nplt)))
1714 close(data_i,iostat=istat)
1716 end subroutine get_ranges
1722 function loc_draw_op(Bokeh_style)
1724 character(len=max_str_ln) :: loc_draw_op
1725 integer,
intent(in),
optional :: bokeh_style
1727 select case (ex_plot_style_loc)
1729 if (n_draw_ops.gt.0)
then
1730 loc_draw_op = trim(draw_ops(mod(iplt-1,n_draw_ops)+1))
1732 loc_draw_op =
'with linespoints linestyle '//&
1733 &trim(
i2str(mod(iplt-1,
size(line_clrs))+1))
1736 if (n_draw_ops.gt.0)
then
1737 loc_draw_op = trim(draw_ops(mod(iplt-1,n_draw_ops)+1))
1739 select case (bokeh_style)
1741 loc_draw_op =
'line_color='//&
1742 &line_clrs(mod(iplt-1,
size(line_clrs))+1)//&
1745 loc_draw_op =
'color='//&
1746 &line_clrs(mod(iplt-1,
size(line_clrs))+1)//&
1753 end function loc_draw_op
1767 real(dp),
intent(in) :: a(:,:,:)
1768 real(dp),
intent(in) :: b(:,:,:)
1769 character(len=*),
intent(in) :: file_name
1770 integer,
intent(in),
optional :: tot_dim(3)
1771 integer,
intent(in),
optional :: loc_offset(3)
1772 character(len=*),
intent(in),
optional :: descr
1773 logical,
intent(in),
optional :: output_message
1776 real(dp),
allocatable :: plot_var(:,:,:,:)
1777 integer :: tot_dim_loc(4)
1778 integer :: loc_offset_loc(4)
1779 character(len=max_str_ln) :: var_names(5)
1780 logical :: output_message_loc
1784 real(dp),
allocatable :: tot_lim(:)
1785 real(dp),
allocatable :: tot_err(:)
1790 tot_dim_loc = [shape(a),5]
1791 if (
present(tot_dim)) tot_dim_loc = [tot_dim,5]
1793 if (
present(loc_offset)) loc_offset_loc = [loc_offset,5]
1796 if (
size(a,1).ne.
size(b,1) .or.
size(a,2).ne.
size(b,2) .or. &
1797 &
size(a,3).ne.
size(b,3))
then
1798 call writo(
'in plot_diff_HDF5, A and B need to have the correct &
1799 &size',persistent=.true.,warning=.true.)
1804 output_message_loc = .false.
1805 if (
present(output_message)) output_message_loc = output_message
1808 if (tot_dim_loc(1).eq.
size(a,1) .and. tot_dim_loc(2).eq.
size(a,2) &
1809 &.and. tot_dim_loc(3).eq.
size(a,3)) ind_plot = .true.
1812 allocate(plot_var(
size(a,1),
size(a,2),
size(a,3),5))
1813 plot_var(:,:,:,1) = a
1814 plot_var(:,:,:,2) = b
1815 plot_var(:,:,:,3) = diff(a,b,shape(a),rel=.true.)
1816 plot_var(:,:,:,4) = log10(abs(diff(a,b,shape(a),rel=.true.)))
1817 plot_var(:,:,:,5) = diff(a,b,shape(a),rel=.false.)
1822 var_names(3) =
'rel v1 - v2'
1823 var_names(4) =
'log abs rel v1 - v2'
1824 var_names(5) =
'abs v1 - v2'
1827 call plot_hdf5_arr(var_names,file_name,plot_var,tot_dim=tot_dim_loc,&
1828 &loc_offset=loc_offset_loc,col=1,descr=descr)
1831 if (output_message_loc)
then
1832 call writo(
'Information about errors:',persistent=.true.)
1835 call stats(tot_dim_loc(1:3),plot_var(:,:,:,3),lim_lo,lim_hi,err_av)
1836 call writo(trim(
r2strt(lim_lo))//
' < rel. err. < '//&
1837 &trim(
r2strt(lim_hi))//
', average value: '//&
1838 &trim(
r2strt(err_av)),persistent=.true.)
1840 call stats(tot_dim_loc(1:3),plot_var(:,:,:,4),lim_lo,lim_hi,err_av)
1841 call writo(trim(
r2strt(lim_lo))//
' < log(abs(rel. err.)) < '//&
1842 &trim(
r2strt(lim_hi))//
', average value: '//&
1843 &trim(
r2strt(err_av)),persistent=.true.)
1845 call stats(tot_dim_loc(1:3),plot_var(:,:,:,5),lim_lo,lim_hi,err_av)
1846 call writo(trim(
r2strt(lim_lo))//
' < abs. err. < '//&
1847 &trim(
r2strt(lim_hi))//
', average value: '//&
1848 &trim(
r2strt(err_av)),persistent=.true.)
1854 function diff(A,B,dims,rel)
result(C)
1858 real(
dp),
intent(in) :: a(:,:,:)
1859 real(
dp),
intent(in) :: b(:,:,:)
1860 integer,
intent(in) :: dims(3)
1861 logical,
intent(in) :: rel
1862 real(
dp) :: c(dims(1),dims(2),dims(3))
1866 c = 2*(a-b)/max(
tol_zero,(abs(a)+abs(b)))
1874 subroutine stats(tot_dims,var,lim_lo,lim_hi,err_av)
1878 integer,
intent(in) :: tot_dims(3)
1879 real(dp),
intent(in) :: var(:,:,:)
1880 real(dp),
intent(inout) :: lim_lo
1881 real(dp),
intent(inout) :: lim_hi
1882 real(dp),
intent(inout) :: err_av
1888 lim_lo = minval(var)
1889 lim_hi = maxval(var)
1893 if (.not.ind_plot)
then
1896 lim_lo = minval(tot_lim)
1899 lim_hi = maxval(tot_lim)
1902 sum_err = sum(tot_err)
1906 err_av = sum_err/product(tot_dims)