program main

use test, only: test_value
use io, only: stdout
use trig, only: d2r

use okada92

implicit none

double precision :: disp(3), strain(3,3), sta_coord(3), evdp, dip, moment(4), lambda, shear_modulus
double precision :: slip_vec(3), wid, len
! double precision :: ux,uy,uz,xin,yin,stdp,dipin,rakin,slip,vp,vs,dens,area
double precision :: uxx, uxy, uxz, uyx, uyy, uyz, uzx, uzy, uzz

#ifdef UNIT_TEST
  lambda = 40.0d0
  shear_modulus = 40.0d0
  call halfspace_vars(lambda,shear_modulus)
  call test_value(a,0.66666666666666663d0,'halfspace_vars(): a')
  call test_value(CA1,0.16666666666666669d0,'halfspace_vars(): CA1')
  call test_value(CA2,0.33333333333333331d0,'halfspace_vars(): CA2')
  call test_value(CB,0.50000000000000011d0,'halfspace_vars(): CB')
  call test_value(CC,0.33333333333333337d0,'halfspace_vars(): CC')
  write(stdout,*) 'subroutine halfspace_vars() passed unit test'
  write(stdout,*)

  dip = 70.0d0
  call dip_vars(dip)
  call test_value(sd,0.93969262078590832d0,'dip_vars(): sd')
  call test_value(cd,0.34202014332566882d0,'dip_vars(): cd')
  call test_value(s2d,0.64278760968653947d0,'dip_vars(): s2d')
  call test_value(c2d,-0.76604444311897790d0,'dip_vars(): c2d')
  call test_value(sdsd,0.88302222155948884d0,'dip_vars(): sdsd')
  call test_value(cdcd,0.11697777844051105d0,'dip_vars(): cdcd')
  call test_value(cdsd,0.32139380484326974d0,'dip_vars(): cdsd')
  write(stdout,*) 'subroutine dip_vars() passed unit test'
  write(stdout,*)

  dip = 70.0d0
  call dip_vars(dip)
  sta_coord(1) = 2.0d0
  sta_coord(2) = 3.0d0
  sta_coord(3) = 0.0d0
  evdp = 4.0d0
  call pt_src_vars(sta_coord(1),sta_coord(2),sta_coord(3),evdp)
  call test_value(p,4.7848309131206399d0,'pt_src_vars(): p')
  call test_value(q,1.4509972890550495d0,'pt_src_vars(): q')
  call test_value(s,4.9925406015355298d0,'pt_src_vars(): s')
  call test_value(t,0.27301710938922441d0,'pt_src_vars(): t')
  call test_value(xx,4.0d0,'pt_src_vars(): xx')
  call test_value(xy,6.0d0,'pt_src_vars(): xy')
  call test_value(yy,9.0d0,'pt_src_vars(): yy')
  call test_value(dd,16.0d0,'pt_src_vars(): dd')
  call test_value(xd,8.0d0,'pt_src_vars(): xd')
  call test_value(xc,8.0d0,'pt_src_vars(): xc')
  call test_value(xq,2.9019945781100991d0,'pt_src_vars(): xq')
  call test_value(yq,4.3529918671651489d0,'pt_src_vars(): yq')
  call test_value(dq,5.8039891562201982d0,'pt_src_vars(): dq')
  call test_value(pq,6.9427766835248459d0,'pt_src_vars(): pq')
  call test_value(xs,9.9850812030710596d0,'pt_src_vars(): xs')
  call test_value(xt,0.54603421877844882d0,'pt_src_vars(): xt')
  call test_value(ys,14.977621804606589d0,'pt_src_vars(): ys')
  call test_value(yt,0.81905132816767323d0,'pt_src_vars(): yt')
  call test_value(qq,2.1053931328451032d0,'pt_src_vars(): qq')
  call test_value(R,5.3851648071345037d0,'pt_src_vars(): R')
  call test_value(R2,28.999999999999996d0,'pt_src_vars(): R2')
  call test_value(R3,156.16977940690060d0,'pt_src_vars(): R3')
  call test_value(R4,840.99999999999989d0,'pt_src_vars(): R4')
  call test_value(R5,4528.9236028001169d0,'pt_src_vars(): R5')
  call test_value(R7,131338.78448120339d0,'pt_src_vars(): R7')
  call test_value(Rd,9.3851648071345046d0,'pt_src_vars(): Rd')
  call test_value(A3,0.58620689655172409d0,'pt_src_vars(): A3')
  call test_value(A5,0.31034482758620685d0,'pt_src_vars(): A5')
  call test_value(A7,3.4482758620689502d-2,'pt_src_vars(): A7')
  call test_value(B3,6.8965517241379226d-2,'pt_src_vars(): B3')
  call test_value(B5,-0.55172413793103470d0,'pt_src_vars(): B5')
  call test_value(B7,-1.1724137931034484d0,'pt_src_vars(): B7')
  call test_value(C3,-0.65517241379310365d0,'pt_src_vars(): C3')
  call test_value(C5,-1.7586206896551726d0,'pt_src_vars(): C5')
  call test_value(C7,-2.8620689655172420d0,'pt_src_vars(): C7')
  call test_value(I1,4.4511857695225003d-003,'pt_src_vars(): I1')
  call test_value(I2,1.4062132997373915d-003,'pt_src_vars(): I2')
  call test_value(I3,1.1400361746955842d-002,'pt_src_vars(): I3')
  call test_value(I4,-6.4425920722921085d-003,'pt_src_vars(): I4')
  call test_value(I5,1.5490988410148277d-002,'pt_src_vars(): I5')
  call test_value(J1,-2.3037701706801377d-003,'pt_src_vars(): J1')
  call test_value(J2,8.3831231772032186d-004,'pt_src_vars(): J2')
  call test_value(J3,2.9153389890690740d-003,'pt_src_vars(): J3')
  call test_value(J4,-1.6706841541556938d-003,'pt_src_vars(): J4')
  call test_value(K1,-1.7024420377114214d-003,'pt_src_vars(): K1')
  call test_value(K2,1.3075030688791331d-004,'pt_src_vars(): K2')
  call test_value(K3,-5.4300227400023550d-003,'pt_src_vars(): K3')
  call test_value(U2,0.18917678161950324d0,'pt_src_vars(): U2')
  call test_value(V2,1.4014492135054364d0,'pt_src_vars(): V2')
  call test_value(W2,1.1288694024054116d0,'pt_src_vars(): W2')
  call test_value(U3,1.3427079288808756d0,'pt_src_vars(): U3')
  call test_value(V3,5.0611389600960157d0,'pt_src_vars(): V3')
  call test_value(W3,1.6847280722065445d0,'pt_src_vars(): W3')
  write(stdout,*) 'subroutine pt_src_vars() passed unit test'
  write(stdout,*)

  dip = 70.0d0
  call dip_vars(dip)
  sta_coord(1) = 2.0d0
  sta_coord(2) = 3.0d0
  sta_coord(3) = 0.0d0
  evdp = 4.0d0
  len = 3.0d0
  wid = 2.0d0
  call rect_src_coords(sta_coord(1),sta_coord(2),sta_coord(3),evdp,wid,len)
  call test_value(p,4.7848309131206399d0,'rect_src_coords(): p')
  call test_value(q,1.4509972890550495d0,'rect_src_coords(): q')
  call test_value(eta_vec(1),3.7848309131206399d0,'rect_src_coords(): eta_vec(1)')
  call test_value(eta_vec(2),5.7848309131206399d0,'rect_src_coords(): eta_vec(2)')
  call test_value(ksi_vec(1),0.5000000000000000d0,'rect_src_coords(): ksi_vec(1)')
  call test_value(ksi_vec(2),3.5000000000000000d0,'rect_src_coords(): ksi_vec(2)')
  write(stdout,*) 'subroutine rect_src_coords() passed unit test'
  call rect_src_vars(ksi_vec(1),eta_vec(1))
  call test_value(ybar,2.6579798566743307d0,'rect_src_vars(): ybar')
  call test_value(dbar,3.0603073792140916d0,'rect_src_vars(): cbar')
  call test_value(cbar,3.0603073792140916d0,'rect_src_vars(): dbar')
  call test_value(R,4.0841569722231199d0,'rect_src_vars(): R')
  call test_value(R2,16.680338173758720d0,'rect_src_vars(): R2')
  call test_value(R3,68.125119451396145d0,'rect_src_vars(): R3')
  ! call test_value(R4,840.99999999999989d0,'rect_src_vars(): R4')
  call test_value(R5,1136.3500305769958d0,'rect_src_vars(): R5')
  ! call test_value(R7,131338.78448120339d0,'rect_src_vars(): R7')
  call test_value(Rd,7.1444643514372110d0,'rect_src_vars(): Rd')
  call test_value(qq,2.1053931328451032d0,'rect_src_vars(): qq')
  call test_value(TH,0.30910022531517090d0,'rect_src_vars(): TH')
  call test_value(Rk,4.5841569722231199d0,'rect_src_vars(): Rk')
  call test_value(X11,5.3411908204318771d-002,'rect_src_vars(): X11')
  call test_value(X32,6.0549197882804537d-003,'rect_src_vars(): X32')
  call test_value(X53,1.3937426671487758d-003,'rect_src_vars(): X53')
  call test_value(logRk,1.5226062223303141d0,'rect_src_vars(): logRk')
  call test_value(Re,7.8689878853437598d0,'rect_src_vars(): Re')
  call test_value(Y11,3.1115637101260175d-002,'rect_src_vars(): Y11')
  call test_value(Y32,2.8335909219300561d-003,'rect_src_vars(): Y32')
  call test_value(Y53,5.6987957651285224d-004,'rect_src_vars(): Y53')
  call test_value(logRe,2.0629294500096305d0,'rect_src_vars(): logRe')
  call test_value(Y0,3.0407239370777661d-002,'rect_src_vars(): Y0')
  call test_value(I1,-24.464842377362366d0,'rect_src_vars(): I1')
  call test_value(I2,1.2599616181546565d0,'rect_src_vars(): I2')
  call test_value(I3,-0.75170987457349314d0,'rect_src_vars(): I3')
  call test_value(I4,26.009469280193130d0,'rect_src_vars(): I4')
  call test_value(Z32,1.2387401944116128d-002,'rect_src_vars(): Z53')
  call test_value(Z53,2.1980039131027692d-003,'rect_src_vars(): Z32')
  call test_value(Z0,1.1837900965840436d-002,'rect_src_vars(): Z0')
  call test_value(D11,3.4271088684950102d-002,'rect_src_vars(): D11')
  call test_value(K1,7.3562247838898460d-003,'rect_src_vars(): K1')
  call test_value(K2,0.11862078136010840d0,'rect_src_vars(): K2')
  call test_value(K3,-0.13432880842342962d0,'rect_src_vars(): K3')
  call test_value(K4,-1.5915028157424500d-003,'rect_src_vars(): K4')
  call test_value(J1,-3.6668885560592841d-002,'rect_src_vars(): J1')
  call test_value(J2,6.3749960046040766d-003,'rect_src_vars(): J2')
  call test_value(J3,3.9930048199627510d-003,'rect_src_vars(): J3')
  call test_value(J4,-1.3985998433743778d-002,'rect_src_vars(): J4')
  call test_value(J6,-1.1485677136015674d-002,'rect_src_vars(): J6')
  call test_value(E2,0.17347006835163784d0,'rect_src_vars(): E2')
  call test_value(F2,4.5587540534369766d-002,'rect_src_vars(): F2')
  call test_value(G2,7.7029412273117825d-002,'rect_src_vars(): G2')
  call test_value(H2,2.8818645950405460d-002,'rect_src_vars(): H2')
  call test_value(P2,8.8840473466618580d-003,'rect_src_vars(): P2')
  call test_value(Q2,1.9608281698289438d-003,'rect_src_vars(): Q2')
  call test_value(E3,0.14892464696612687d0,'rect_src_vars(): E3')
  call test_value(F3,3.9258435961432761d-002,'rect_src_vars(): F3')
  call test_value(G3,6.3422754458309291d-002,'rect_src_vars(): G3')
  call test_value(H3,2.4055253239807024d-002,'rect_src_vars(): H3')
  call test_value(P3,1.2387401944116128d-002,'rect_src_vars(): P3')
  call test_value(Q3,1.7300629665110766d-002,'rect_src_vars(): Q3')
  write(stdout,*) 'subroutine rect_src_vars() passed unit test'
  write(stdout,*)

  ! Half-space
  lambda = 40.0d0
  shear_modulus = 40.0d0
  call halfspace_vars(lambda,shear_modulus)
  ! Fault dip
  dip = 70.0d0
  call dip_vars(dip)
  ! Geometry
  sta_coord(1) = -12.0d0
  sta_coord(2) = 6.2d0
  sta_coord(3) = 13.0d0
  evdp = 7.0d0
  wid = 3.2d0
  len = 6.5d0
  call rect_src_coords(sta_coord(1),sta_coord(2),-sta_coord(3),evdp,wid,len)
  call rect_src_vars(ksi_vec(2),eta_vec(1))
  disp = uA_ss_rect()
  call test_value(disp(1), 0.72547327384003113d0,'uA_ss_rect()(1)')
  call test_value(disp(2),-1.5927062399294267d-2,'uA_ss_rect()(2)')
  call test_value(disp(3), 0.61666031833564772d0,'uA_ss_rect()(3)')
  disp = uB_ss_rect()
  call test_value(disp(1),-11.993002387482605d0,'uB_ss_rect()(1)')
  call test_value(disp(2),0.11463965563116107d0,'uB_ss_rect()(2)')
  call test_value(disp(3),-0.99745107154736778d0,'uB_ss_rect()(3)')
  disp = uC_ss_rect()
  call test_value(disp(1),-1.4461853668621687d-003,'uC_ss_rect()(1)')
  call test_value(disp(2),5.0207492407645738d-003,'uC_ss_rect()(2)')
  call test_value(disp(3),-2.0079399191054680d-002,'uC_ss_rect()(3)')
  call rect_src_coords(sta_coord(1),sta_coord(2),-sta_coord(3),evdp,wid,len)
  call rect_src_vars(ksi_vec(1),eta_vec(2))
  disp = uA_ds_rect()
  call test_value(disp(1),-1.2424808371844004d-2,'uA_ds_rect()(1)')
  call test_value(disp(2), 0.72190403451379292d0,'uA_ds_rect()(2)')
  call test_value(disp(3), 0.41256779800532334d0,'uA_ds_rect()(3)')
  disp = uB_ds_rect()
  call test_value(disp(1),-0.24780594882338475d0,'uB_ds_rect()(1)')
  call test_value(disp(2),-1.3701174902198523d0,'uB_ds_rect()(2)')
  call test_value(disp(3),-3.8188429813460623d0,'uB_ds_rect()(3)')
  disp = uC_ds_rect()
  call test_value(disp(1), 5.1793376586471750d-003,'uC_ds_rect()(1)')
  call test_value(disp(2), 8.6684162903643731d-003,'uC_ds_rect()(2)')
  call test_value(disp(3),-7.2807977939472246d-002,'uC_ds_rect()(3)')
  write(stdout,*) 'displacement functions for rectangular source passed unit tests'
  write(stdout,*)

  ! Half-space
  lambda = 40.0d0
  shear_modulus = 40.0d0
  call halfspace_vars(lambda,shear_modulus)
  ! Fault dip
  dip = 70.0d0
  call dip_vars(dip)
  ! Geometry
  sta_coord(1) = -12.0d0
  sta_coord(2) = 6.2d0
  sta_coord(3) = 13.0d0
  evdp = 7.0d0
  wid = 3.2d0
  len = 6.5d0
  call rect_src_coords(sta_coord(1),sta_coord(2),-sta_coord(3),evdp,wid,len)
  call rect_src_vars(ksi_vec(1),eta_vec(1))
  ! X-derivatives of components of displacement
  disp = duAdx_ss_rect()
  call test_value(disp(1), 3.4306426504990750d-004,'duAdx_ss_rect(): (1)')
  call test_value(disp(2),-3.4508486032947076d-004,'duAdx_ss_rect(): (2)')
  call test_value(disp(3),-2.3607160020503188d-003,'duAdx_ss_rect(): (3)')
  disp = duAdx_ds_rect()
  call test_value(disp(1),-3.4508486032947076d-004,'duAdx_ds_rect(): (1)')
  call test_value(disp(2), 9.0562367901465455d-004,'duAdx_ds_rect(): (2)')
  call test_value(disp(3), 6.7897701068703959d-003,'duAdx_ds_rect(): (3)')
  ! disp = duAdx_ts_rect()
  ! write(0,*) 'duAdx_ts_rect()',disp
  disp = duBdx_ss_rect()
  call test_value(disp(1), 6.9599894428799694d-004,'duBdx_ss_rect(): (1)')
  call test_value(disp(2), 1.9194218761293762d-003,'duBdx_ss_rect(): (2)')
  call test_value(disp(3), 2.0942277171944572d-003,'duBdx_ss_rect(): (3)')
  disp = duBdx_ds_rect()
  call test_value(disp(1), 2.7417691868664293d-003,'duBdx_ds_rect(): (1)')
  call test_value(disp(2),-5.1586589356422364d-003,'duBdx_ds_rect(): (2)')
  call test_value(disp(3),-6.7076617173327598d-004,'duBdx_ds_rect(): (3)')
  ! disp = duBdx_ts_rect()
  ! write(0,*) 'duBdx_ts_rect()',disp
  disp = duCdx_ss_rect()
  call test_value(disp(1), 4.2126939364670769d-005,'duCdx_ss_rect(): (1)')
  call test_value(disp(2), 1.1209090419535319d-004,'duCdx_ss_rect(): (2)')
  call test_value(disp(3),-3.4475338947556833d-004,'duCdx_ss_rect(): (3)')
  disp = duCdx_ds_rect()
  call test_value(disp(1), 1.6966631277932049d-004,'duCdx_ds_rect(): (1)')
  call test_value(disp(2),-1.4986921869860253d-004,'duCdx_ds_rect(): (2)')
  call test_value(disp(3), 1.1331156239589844d-003,'duCdx_ds_rect(): (3)')
  ! disp = duCdx_ts_rect()
  ! write(0,*) 'duCdx_ts_rect()',disp
  ! Y-derivatives of components of displacement
  disp = duAdy_ss_rect()
  call test_value(disp(1),2.8890966553217295d-002,'duAdy_ss_rect(): (1)')
  call test_value(disp(2),1.2845371264737372d-002,'duAdy_ss_rect(): (2)')
  call test_value(disp(3),2.7617711810937568d-003,'duAdy_ss_rect(): (3)')
  disp = duAdy_ds_rect()
  call test_value(disp(1), 1.2845371264737372d-002,'duAdy_ds_rect(): (1)')
  call test_value(disp(2), 6.0051670894042106d-002,'duAdy_ds_rect(): (2)')
  call test_value(disp(3), 6.8785272909704029d-003,'duAdy_ds_rect(): (3)')
  ! disp = duAdy_ts_rect()
  disp = duBdy_ss_rect()
  call test_value(disp(1),-5.4892739626334362d-002,'duBdy_ss_rect(): (1)')
  call test_value(disp(2),-2.7969238876413290d-002,'duBdy_ss_rect(): (2)')
  call test_value(disp(3),-3.1020145458715599d-003,'duBdy_ss_rect(): (3)')
  disp = duBdy_ds_rect()
  call test_value(disp(1),-3.8965904427687062d-002,'duBdy_ds_rect(): (1)')
  call test_value(disp(2),-0.14704676790658636d0,'duBdy_ds_rect(): (2)')
  call test_value(disp(3),-9.1052106567444321d-003,'duBdy_ds_rect(): (3)')
  ! disp = duBdy_ts_rect()
  disp = duCdy_ss_rect()
  call test_value(disp(1),3.6229956319952023d-004,'duCdy_ss_rect(): (1)')
  call test_value(disp(2),2.7645676601602590d-004,'duCdy_ss_rect(): (2)')
  call test_value(disp(3),5.2251433363535108d-004,'duCdy_ss_rect(): (3)')
  disp = duCdy_ds_rect()
  call test_value(disp(1),-1.1160386819908798d-003,'duCdy_ds_rect(): (1)')
  call test_value(disp(2),-6.3814177065406037d-004,'duCdy_ds_rect(): (2)')
  call test_value(disp(3), 2.7583611996049840d-003,'duCdy_ds_rect(): (3)')
  ! disp = duCdy_ts_rect()
  ! Z-derivatives of components of displacement
  disp = duAdz_ss_rect()
  call test_value(disp(1), 8.5466041996219398d-003,'duAdz_ss_rect(): (1)')
  call test_value(disp(2), 4.2102278432794181d-003,'duAdz_ss_rect(): (2)')
  call test_value(disp(3),-6.2203199241805318d-003,'duAdz_ss_rect(): (3)')
  disp = duAdz_ds_rect()
  call test_value(disp(1), 4.2102278432794181d-003,'duAdz_ds_rect(): (1)')
  call test_value(disp(2), 1.7605900194629188d-002,'duAdz_ds_rect(): (2)')
  call test_value(disp(3),-1.2506584807178457d-002,'duAdz_ds_rect(): (3)')
  ! disp = duAdy_ts_rect()
  disp = duBdz_ss_rect()
  call test_value(disp(1),-1.7349770792071152d-002,'duBdz_ss_rect(): (1)')
  call test_value(disp(2),-1.0130288362773653d-002,'duBdz_ss_rect(): (2)')
  call test_value(disp(3), 1.0421179348291247d-002,'duBdz_ss_rect(): (3)')
  disp = duBdz_ds_rect()
  call test_value(disp(1),-9.6899784763487366d-003,'duBdz_ds_rect(): (1)')
  call test_value(disp(2),-4.0686184682977758d-002,'duBdz_ds_rect(): (2)')
  call test_value(disp(3),-2.2296367147699551d-003,'duBdz_ds_rect(): (3)')
  ! disp = duBdy_ts_rect()
  disp = duCdz_ss_rect()
  call test_value(disp(1),-3.6434767894171044d-005,'duCdz_ss_rect(): (1)')
  call test_value(disp(2), 2.3748100977727307d-004,'duCdz_ss_rect(): (2)')
  call test_value(disp(3),-7.0849983200360817d-004,'duCdz_ss_rect(): (3)')
  disp = duCdz_ds_rect()
  call test_value(disp(1),-1.5620330749121738d-004,'duCdz_ds_rect(): (1)')
  call test_value(disp(2), 5.7029847969910998d-004,'duCdz_ds_rect(): (2)')
  call test_value(disp(3),-5.4083583036719013d-003,'duCdz_ds_rect(): (3)')
  ! disp = duCdz_ts_rect()
  write(stdout,*) 'displacement derivative functions for rectangular source passed unit tests'
  write(stdout,*)

#endif

!----
! Point source tests
!----
! Point source example from Okada (1985)
lambda = 40.0d0            ! lambda = shear_modulus implies
shear_modulus = 40.0d0     !     a Poisson's ratio of 0.25
dip = 70.0d0
sta_coord(1) = 2.0d0
sta_coord(2) = 3.0d0
sta_coord(3) = 0.0d0
evdp = 4.0d0

! Strike-slip
moment(1) = 40.0d0
moment(2) =  0.0d0
moment(3) =  0.0d0
call o92_pt_disp(disp,sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(disp(1),-9.4474186064630631d-4,'o92_pt_disp(): disp(1) (strike-slip)')
call test_value(disp(2),-1.0229885427654646d-3,'o92_pt_disp(): disp(2) (strike-slip)')
call test_value(disp(3),-7.4201246532591128d-4,'o92_pt_disp(): disp(3) (strike-slip)')
write(stdout,*) 'subroutine o92_pt_disp() passed strike-slip unit test'
call o92_pt_partials(uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz,&
                     sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(uxx,-2.2862094355883986d-004,'o92_pt_partials(): uxx (strike-slip)')
call test_value(uxy,-1.4246398245262621d-004,'o92_pt_partials(): uxy (strike-slip)')
call test_value(uxz, 6.2590905699873601d-005,'o92_pt_partials(): uxz (strike-slip)')
call test_value(uyx,-2.0511017314976633d-004,'o92_pt_partials(): uyx (strike-slip)')
call test_value(uyy,-3.0067825210386353d-004,'o92_pt_partials(): uyy (strike-slip)')
call test_value(uyz, 1.6933001353308135d-004,'o92_pt_partials(): uyz (strike-slip)')
call test_value(uzx,-6.2590905699873859d-005,'o92_pt_partials(): uzx (strike-slip)')
call test_value(uzy,-1.6933001353308129d-004,'o92_pt_partials(): uzy (strike-slip)')
call test_value(uzz, 1.7643306522090115d-004,'o92_pt_partials(): uzz (strike-slip)')
write(stdout,*) 'subroutine o92_pt_partials() passed strike-slip unit test'
call o92_pt_strain(strain,sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(strain(1,1),-2.2862094355883986d-004,'o92_pt_strain(): strain(1,1) (strike-slip)')
call test_value(strain(1,2),-1.7378707780119627d-004,'o92_pt_strain(): strain(1,2) (strike-slip)')
call test_value(strain(1,3),-1.2874900798265365d-019,'o92_pt_strain(): strain(1,3) (strike-slip)')
call test_value(strain(2,1),-1.7378707780119627d-004,'o92_pt_strain(): strain(2,1) (strike-slip)')
call test_value(strain(2,2),-3.0067825210386353d-004,'o92_pt_strain(): strain(2,2) (strike-slip)')
call test_value(strain(2,3), 2.7105054312137611d-020,'o92_pt_strain(): strain(2,3) (strike-slip)')
call test_value(strain(3,1),-1.2874900798265365d-019,'o92_pt_strain(): strain(3,1) (strike-slip)')
call test_value(strain(3,2), 2.7105054312137611d-020,'o92_pt_strain(): strain(3,2) (strike-slip)')
call test_value(strain(3,3), 1.7643306522090115d-004,'o92_pt_strain(): strain(3,3) (strike-slip)')
write(stdout,*) 'subroutine o92_pt_strain() passed strike-slip unit test'

! Dip-slip
moment(1) =  0.0d0
moment(2) = 40.0d0
moment(3) =  0.0d0
call o92_pt_disp(disp,sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(disp(1),-1.1723215364075840d-3,'o92_pt_disp(): disp(1) (dip-slip)')
call test_value(disp(2),-2.0819985410726201d-3,'o92_pt_disp(): disp(2) (dip-slip)')
call test_value(disp(3),-2.5315947108865925d-3,'o92_pt_disp(): disp(3) (dip-slip)')
write(stdout,*) 'subroutine o92_pt_disp() passed dip-slip unit test'
call o92_pt_partials(uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz,&
                     sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(uxx,-1.5259408689661085d-004,'o92_pt_partials(): uxx (dip-slip)')
call test_value(uxy,-3.5441800891339358d-004,'o92_pt_partials(): uxy (dip-slip)')
call test_value(uxz,-8.7070529207727333d-004,'o92_pt_partials(): uxz (dip-slip)')
call test_value(uyx, 6.9826591999645748d-004,'o92_pt_partials(): uyx (dip-slip)')
call test_value(uyy,-1.1537526800074350d-003,'o92_pt_partials(): uyy (dip-slip)')
call test_value(uyz, 6.3453607666584869d-004,'o92_pt_partials(): uyz (dip-slip)')
call test_value(uzx, 8.7070529207727322d-004,'o92_pt_partials(): uzx (dip-slip)')
call test_value(uzy,-6.3453607666584847d-004,'o92_pt_partials(): uzy (dip-slip)')
call test_value(uzz, 4.3544892230134886d-004,'o92_pt_partials(): uzz (dip-slip)')
write(stdout,*) 'subroutine o92_pt_partials() passed dip-slip unit test'
call o92_pt_strain(strain,sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(strain(1,1),-1.5259408689661085d-004,'o92_pt_strain(): strain(1,1) (dip-slip)')
call test_value(strain(1,2), 1.7192395554153195d-004,'o92_pt_strain(): strain(1,2) (dip-slip)')
call test_value(strain(1,3),-5.4210108624275222d-020,'o92_pt_strain(): strain(1,3) (dip-slip)')
call test_value(strain(2,1), 1.7192395554153195d-004,'o92_pt_strain(): strain(2,1) (dip-slip)')
call test_value(strain(2,2),-1.1537526800074350d-003,'o92_pt_strain(): strain(2,2) (dip-slip)')
call test_value(strain(2,3), 1.0842021724855044d-019,'o92_pt_strain(): strain(2,3) (dip-slip)')
call test_value(strain(3,1),-5.4210108624275222d-020,'o92_pt_strain(): strain(3,1) (dip-slip)')
call test_value(strain(3,2), 1.0842021724855044d-019,'o92_pt_strain(): strain(3,2) (dip-slip)')
call test_value(strain(3,3), 4.3544892230134886d-004,'o92_pt_strain(): strain(3,3) (dip-slip)')
write(stdout,*) 'subroutine o92_pt_strain() passed dip-slip unit test'

! Tensile-slip
moment(1) =  0.0d0
moment(2) =  0.0d0
moment(3) = 40.0d0
call o92_pt_disp(disp,sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(disp(1),-3.5716333691232548d-4,'o92_pt_disp(): disp(1) (tensile-slip)')
call test_value(disp(2), 3.5310854913572190d-4,'o92_pt_disp(): disp(2) (tensile-slip)')
call test_value(disp(3),-2.0068126969335952d-4,'o92_pt_disp(): disp(3) (tensile-slip)')
write(stdout,*) 'subroutine o92_pt_disp() passed tensile-slip unit test'
call o92_pt_partials(uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz,&
                     sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(uxx,-1.3597209844017221d-004,'o92_pt_partials(): uxx (tensile-slip)')
call test_value(uxy, 5.0725453276876352d-004,'o92_pt_partials(): uxy (tensile-slip)')
call test_value(uxz,-7.5405344488748576d-005,'o92_pt_partials(): uxz (tensile-slip)')
call test_value(uyx,-6.7733492861951009d-005,'o92_pt_partials(): uyx (tensile-slip)')
call test_value(uyy, 6.8111287703614769d-004,'o92_pt_partials(): uyy (tensile-slip)')
call test_value(uyz,-8.1037165493104131d-004,'o92_pt_partials(): uyz (tensile-slip)')
call test_value(uzx, 7.5405344488748590d-005,'o92_pt_partials(): uzx (tensile-slip)')
call test_value(uzy, 8.1037165493104120d-004,'o92_pt_partials(): uzy (tensile-slip)')
call test_value(uzz,-1.8171359286532499d-004,'o92_pt_partials(): uzz (tensile-slip)')
write(stdout,*) 'subroutine o92_pt_partials() passed tensile-slip unit test'
call o92_pt_strain(strain,sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(strain(1,1),-1.3597209844017221d-004,'o92_pt_strain(): strain(1,1) (tensile-slip)')
call test_value(strain(1,2), 2.1976051995340625d-004,'o92_pt_strain(): strain(1,2) (tensile-slip)')
call test_value(strain(1,3), 6.7762635780344027d-021,'o92_pt_strain(): strain(1,3) (tensile-slip)')
call test_value(strain(2,1), 2.1976051995340625d-004,'o92_pt_strain(): strain(2,1) (tensile-slip)')
call test_value(strain(2,2), 6.8111287703614769d-004,'o92_pt_strain(): strain(2,2) (tensile-slip)')
call test_value(strain(2,3),-5.4210108624275222d-020,'o92_pt_strain(): strain(2,3) (tensile-slip)')
call test_value(strain(3,1), 6.7762635780344027d-021,'o92_pt_strain(): strain(3,1) (tensile-slip)')
call test_value(strain(3,2),-5.4210108624275222d-020,'o92_pt_strain(): strain(3,2) (tensile-slip)')
call test_value(strain(3,3),-1.8171359286532499d-004,'o92_pt_strain(): strain(3,3) (tensile-slip)')
write(stdout,*) 'subroutine o92_pt_strain() passed tensile-slip unit test'
write(stdout,*) 'point source example from Okada (1985) passed unit test'
write(stdout,*)

! Point source, station depth greater than zero, oblique slip
lambda = 40.0d0
shear_modulus = 40.0d0
sta_coord(1) = -12.0d0
sta_coord(2) = 6.2d0
sta_coord(3) = 13.0d0
evdp = 7.0d0
moment(1) = -27.0d0
moment(2) = 13.0d0
moment(3) = 0.0d0

! Displacement
call o92_pt_disp(disp,sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(disp(1), 1.9245175992850267d-004,'o92_pt_disp(): disp(1) (z>0; oblique-slip)')
call test_value(disp(2),-1.9098044959056914d-004,'o92_pt_disp(): disp(2) (z>0; oblique-slip)')
call test_value(disp(3), 8.1576627770328566d-005,'o92_pt_disp(): disp(3) (z>0; oblique-slip)')

! Partial derivatives of displacement
call o92_pt_partials(uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz,&
                     sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(uxx, 2.0874574675642794d-005,'o92_pt_partials(): uxx (z>0; oblique-slip)')
call test_value(uxy, 2.3620493327193696d-006,'o92_pt_partials(): uxy (z>0; oblique-slip)')
call test_value(uxz, 2.2899507819320573d-005,'o92_pt_partials(): uxz (z>0; oblique-slip)')
call test_value(uyx,-2.1572559953257881d-005,'o92_pt_partials(): uyx (z>0; oblique-slip)')
call test_value(uyy,-6.8032602399881741d-006,'o92_pt_partials(): uyy (z>0; oblique-slip)')
call test_value(uyz,-2.1249535452408235d-005,'o92_pt_partials(): uyz (z>0; oblique-slip)')
call test_value(uzx, 1.2032014067111290d-005,'o92_pt_partials(): uzx (z>0; oblique-slip)')
call test_value(uzy,-6.7523163638339241d-006,'o92_pt_partials(): uzy (z>0; oblique-slip)')
call test_value(uzz,-1.7420298917986620d-006,'o92_pt_partials(): uzz (z>0; oblique-slip)')

! Strain tensor
call o92_pt_strain(strain,sta_coord,evdp,dip,moment,lambda,shear_modulus)
call test_value(strain(1,1), 2.0874574675642794d-5,'o92_pt_strain(): e11 (z>0; oblique-slip)')
call test_value(strain(1,2),-9.6052553102692557d-6,'o92_pt_strain(): e12 (z>0; oblique-slip)')
call test_value(strain(1,3), 1.7465760943215932d-5,'o92_pt_strain(): e13 (z>0; oblique-slip)')
call test_value(strain(2,1),-9.6052553102692557d-6,'o92_pt_strain(): e21 (z>0; oblique-slip)')
call test_value(strain(2,2),-6.8032602399881741d-6,'o92_pt_strain(): e22 (z>0; oblique-slip)')
call test_value(strain(2,3),-1.4000925908121080d-5,'o92_pt_strain(): e23 (z>0; oblique-slip)')
call test_value(strain(3,1), 1.7465760943215932d-5,'o92_pt_strain(): e31 (z>0; oblique-slip)')
call test_value(strain(3,2),-1.4000925908121080d-5,'o92_pt_strain(): e32 (z>0; oblique-slip)')
call test_value(strain(3,3),-1.7420298917986620d-6,'o92_pt_strain(): e33 (z>0; oblique-slip)')
write(stdout,*) 'point source example with depth>0, oblique slip passed unit test'
write(stdout,*)

!----
! Rectangular source
!----
! Example from Okada (1985)
lambda = 40.0d0            ! lambda = shear_modulus implies
shear_modulus = 40.0d0     !     a Poisson's ratio of 0.25
dip = 70.0d0
sta_coord(1) = 2.0d0
sta_coord(2) = 3.0d0
sta_coord(3) = 0.0d0
evdp = 4.0d0
len = 3.0d0
wid = 2.0d0

! Original Okada derivation has fault origin at the bottom left corner, but I put the fault
! center at the coordinates. To compare with Okada (1985), shift fault.
sta_coord(1) = sta_coord(1)-len/2.0d0
sta_coord(2) = sta_coord(2)-wid/2.0d0*dcos(dip*d2r)
evdp = evdp-wid/2.0d0*dsin(dip*d2r)

! Strike-slip
slip_vec(1) = 1.0d0
slip_vec(2) = 0.0d0
slip_vec(3) = 0.0d0
call o92_rect_disp(disp,sta_coord,evdp,dip,slip_vec,wid,len,lambda,shear_modulus)
call test_value(disp(1),-8.6891650042561549d-003,'o92_rect_disp(): disp(1) (strike-slip)')
call test_value(disp(2),-4.2975821897415613d-003,'o92_rect_disp(): disp(2) (strike-slip)')
call test_value(disp(3),-2.7474058276388967d-003,'o92_rect_disp(): disp(3) (strike-slip)')
write(stdout,*) 'subroutine o92_rect_disp() passed strike-slip unit test'

! Dip-slip
slip_vec(1) = 0.0d0
slip_vec(2) = 1.0d0
slip_vec(3) = 0.0d0
call o92_rect_disp(disp,sta_coord,evdp,dip,slip_vec,wid,len,lambda,shear_modulus)
call test_value(disp(1),-4.6823487628354045d-003,'o92_rect_disp(): disp(1) (dip-slip)')
call test_value(disp(2),-3.5267267968718020d-002,'o92_rect_disp(): disp(2) (dip-slip)')
call test_value(disp(3),-3.5638557673268470d-002,'o92_rect_disp(): disp(3) (dip-slip)')
write(stdout,*) 'subroutine o92_rect_disp() passed dip-slip unit test'

! Tensile-slip
slip_vec(1) = 0.0d0
slip_vec(2) = 0.0d0
slip_vec(3) = 1.0d0
call o92_rect_disp(disp,sta_coord,evdp,dip,slip_vec,wid,len,lambda,shear_modulus)
call test_value(disp(1),-2.6599600964428191d-004,'o92_rect_disp(): disp(1) (tensile-slip)')
call test_value(disp(2), 1.0564074876983659d-002,'o92_rect_disp(): disp(2) (tensile-slip)')
call test_value(disp(3), 3.2141931142207676d-003,'o92_rect_disp(): disp(3) (tensile-slip)')
write(stdout,*) 'subroutine o92_rect_disp() passed tensile-slip unit test'
write(stdout,*) 'rectangular source example from Okada (1985) passed unit test'
write(stdout,*)

! Example with rectangular source, station depth greater than zero
lambda = 40.0d0
shear_modulus = 40.0d0
dip = 70.0d0
sta_coord(1) = -12.0d0
sta_coord(2) = 6.2d0
sta_coord(3) = 13.0d0
evdp = 7.0d0
wid = 3.2d0
len = 6.5d0

! Strike-slip displacement
slip_vec(1) = 1.0d0
slip_vec(2) = 0.0d0
slip_vec(3) = 0.0d0
call o92_rect_disp(disp,sta_coord,evdp,dip,slip_vec,wid,len,lambda,shear_modulus)
call test_value(disp(1),-6.6858712769214453d-003,'o92_rect_disp(): disp(1) (z>0; strike-slip)')
call test_value(disp(2), 5.6833739558216096d-003,'o92_rect_disp(): disp(2) (z>0; strike-slip)')
call test_value(disp(3),-3.8337220569204383d-003,'o92_rect_disp(): disp(3) (z>0; strike-slip)')
write(stdout,*)

! Dip-slip displacement
slip_vec(1) = 0.0d0
slip_vec(2) = 1.0d0
slip_vec(3) = 0.0d0
call o92_rect_disp(disp,sta_coord,evdp,dip,slip_vec,wid,len,lambda,shear_modulus)
call test_value(disp(1),-1.4891278409154987d-003,'o92_rect_disp(): disp(1) (z>0; dip-slip)')
call test_value(disp(2),-6.6628749834942423d-004,'o92_rect_disp(): disp(2) (z>0; dip-slip)')
call test_value(disp(3),-2.6808905339339140d-003,'o92_rect_disp(): disp(3) (z>0; dip-slip)')
write(stdout,*)

! Strike-slip partial derivatives of displacement
slip_vec(1) = 1.0d0
slip_vec(2) = 0.0d0
slip_vec(3) = 0.0d0
call o92_rect_partials(uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz, &
                       sta_coord,evdp,dip,slip_vec,wid,len,lambda,shear_modulus)
call test_value(uxx,-7.4224535411782154d-004,'o92_rect_partials(): uxx (z>0; strike-slip)')
call test_value(uyx, 7.0705631705157462d-004,'o92_rect_partials(): uyx (z>0; strike-slip)')
call test_value(uzx,-5.7940045023416157d-004,'o92_rect_partials(): uzx (z>0; strike-slip)')
call test_value(uxy, 6.4347351147793649d-005,'o92_rect_partials(): uxy (z>0; strike-slip)')
call test_value(uyy, 1.5409553602864544d-004,'o92_rect_partials(): uyy (z>0; strike-slip)')
call test_value(uzy, 2.2074703722434528d-004,'o92_rect_partials(): uzy (z>0; strike-slip)')
call test_value(uxz,-6.3176025002264938d-004,'o92_rect_partials(): uxz (z>0; strike-slip)')
call test_value(uyz, 5.5337876385493990d-004,'o92_rect_partials(): uyz (z>0; strike-slip)')
call test_value(uzz, 1.0685389278657297d-004,'o92_rect_partials(): uzz (z>0; strike-slip)')

slip_vec(1) = 0.0d0
slip_vec(2) = 1.0d0
slip_vec(3) = 0.0d0
call o92_rect_partials(uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz, &
                       sta_coord,evdp,dip,slip_vec,wid,len,lambda,shear_modulus)
call test_value(uxx,-2.7933167653439758d-004,'o92_rect_partials(): uxx (z>0; dip-slip)')
call test_value(uyx, 9.6965877919922243d-005,'o92_rect_partials(): uyx (z>0; dip-slip)')
call test_value(uzx,-4.7400064111388521d-004,'o92_rect_partials(): uzx (z>0; dip-slip)')
call test_value(uxy, 1.8756328884862554d-004,'o92_rect_partials(): uxy (z>0; dip-slip)')
call test_value(uyy,-6.2549348110294185d-005,'o92_rect_partials(): uyy (z>0; dip-slip)')
call test_value(uzy,-2.0932522925848662d-005,'o92_rect_partials(): uzy (z>0; dip-slip)')
call test_value(uxz, 2.4695902664444695d-004,'o92_rect_partials(): uxz (z>0; dip-slip)')
call test_value(uyz,-3.3607226779516899d-004,'o92_rect_partials(): uyz (z>0; dip-slip)')
call test_value(uzz, 1.6913995963324557d-004,'o92_rect_partials(): uzz (z>0; dip-slip)')

! slip_vec(1) = 0.0d0
! slip_vec(2) = 0.0d0
! slip_vec(3) = 1.0d0
! call o92_rect_partials(uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz, &
!                        sta_coord,evdp,dip,slip_vec,wid,len,lambda,shear_modulus)
! call test_value(uxx,-7.4224535411782154d-004,'o92_rect_partials(): uxx (z>0; tensile-slip)')
! call test_value(uyx, 7.0705631705157462d-004,'o92_rect_partials(): uyx (z>0; tensile-slip)')
! call test_value(uzx,-5.7940045023416157d-004,'o92_rect_partials(): uzx (z>0; tensile-slip)')
! call test_value(uxy, 6.4347351147793649d-005,'o92_rect_partials(): uxy (z>0; tensile-slip)')
! call test_value(uyy, 1.5409553602864544d-004,'o92_rect_partials(): uyy (z>0; tensile-slip)')
! call test_value(uzy, 2.2074703722434528d-004,'o92_rect_partials(): uzy (z>0; tensile-slip)')
! call test_value(uxz,-6.3176025002264938d-004,'o92_rect_partials(): uxz (z>0; tensile-slip)')
! call test_value(uyz, 5.5337876385493990d-004,'o92_rect_partials(): uyz (z>0; tensile-slip)')
! call test_value(uzz, 1.0685389278657297d-004,'o92_rect_partials(): uzz (z>0; tensile-slip)')

slip_vec(1) = 3.0d0
slip_vec(2) = -8.0d0
slip_vec(3) = 0.0d0
call o92_rect_strain(strain,sta_coord,evdp,dip,slip_vec,wid,len,lambda,shear_modulus)
call test_value(strain(1,1), 7.9173499217153817d-006,'o92_rect_strain(): e11 (z>0; oblique-slip)')
call test_value(strain(1,2), 1.8988835224861660d-005,'o92_rect_strain(): e12 (z>0; oblique-slip)')
call test_value(strain(1,3),-9.0857459250746377d-004,'o92_rect_strain(): e13 (z>0; oblique-slip)')
call test_value(strain(2,1), 1.8988835224861660d-005,'o92_rect_strain(): e21 (z>0; oblique-slip)')
call test_value(strain(2,2), 9.6268139296828981d-004,'o92_rect_strain(): e22 (z>0; oblique-slip)')
call test_value(strain(2,3), 2.5892078645029983d-003,'o92_rect_strain(): e23 (z>0; oblique-slip)')
call test_value(strain(3,1),-9.0857459250746377d-004,'o92_rect_strain(): e31 (z>0; oblique-slip)')
call test_value(strain(3,2), 2.5892078645029983d-003,'o92_rect_strain(): e32 (z>0; oblique-slip)')
call test_value(strain(3,3),-1.0325579987062454d-003,'o92_rect_strain(): e33 (z>0; oblique-slip)')

!----
! Compare with original code results
!----
! sta_coord(1) = 1111.949d0
! sta_coord(2) = -1111.949d0
! sta_coord(3) = 0.0d0
! evdp = 10000.0d0
! dip = 90.0d0
! wid = 4000.0d0
! len = 6000.0d0
! lambda = 4.0d10
! shear_modulus = 4.0d10
! slip_vec(1) = -1.0d0
! slip_vec(2) = 0.0d0
! slip_vec(3) = 0.0d0
! xin = sta_coord(1)
! yin = sta_coord(2)
! stdp = sta_coord(3)
! dipin = dip
! rakin = datan2(moment(2),moment(1))/d2r
! area = 1.0d0
! slip = dsqrt(moment(1)*moment(1)+moment(2)*moment(2))/shear_modulus/area
! dens = 1.0d0
! vp = dsqrt((lambda+2.0d0*shear_modulus)/dens)
! vs = dsqrt(shear_modulus/dens)
! call o92pt(ux,uy,uz,xin,yin,stdp,evdp,dipin,rakin,area,slip,vp,vs,dens)
! write(6,*) ux,uy,uz
! call o92ptstn(strain,xin,yin,stdp,evdp,dipin,rakin,area,slip,vp,vs,dens)
! write(6,*) strain(1,:)
! write(6,*) strain(2,:)
! write(6,*) strain(3,:)

! rakin = datan2(slip_vec(2),slip_vec(1))/d2r
! slip = dsqrt(slip_vec(1)*slip_vec(1)+slip_vec(2)*slip_vec(2))
! write(0,*) 'sta_coord:',sta_coord,'    x,y,z:',xin,yin,stdp
! write(0,*) 'evdp:',evdp
! write(0,*) 'dip:',dip,dipin
! write(0,*) 'slip_vec:',slip_vec,'    slip,rakin:',slip,rakin
! write(0,*) 'wid:',wid
! write(0,*) 'len:',len
! write(0,*) 'shear:',shear_modulus
! write(0,*) 'vp,vs,dens:',vp,vs,dens
! call o92_rect_disp(disp,sta_coord,evdp,dip,slip_vec,wid,len,lambda,shear_modulus)
! write(0,*) disp
! call o92rect(ux,uy,uz,xin,yin,stdp,evdp,dipin,rakin,wid,len,slip,vp,vs,dens)
! write(0,*) ux,uy,uz
! call o92rectstn(strain,xin,yin,stdp,evdp,dipin,rakin,wid,len,slip,vp,vs,dens)

write(stdout,*) 'okada92_module unit test passed'
end


      SUBROUTINE SELECT(L,EI,DE,EMAIN,REL,VR,RNOT,DX,JRI,ZZ)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
!        Arguments
!
      INTEGER            JRI, L
      DOUBLE PRECISION   DE, DX, EI, RNOT, ZZ
      DOUBLE PRECISION   VR(*)
      LOGICAL            REL
      CHARACTER*4        EMAIN
!
!     ..................................................................
!
!        automatic energy search:
!           finds energy parameter by looking where value and slope of
!           wavefunction changes sign  (E-top:    u  changes sign)
!                                      (E-bottom: u' changes sign)
!
!     ..................................................................
!
!        Local Scalars
!
      INTEGER            NODE
      DOUBLE PRECISION   DDNOLD, DNOLD, DU, DUDN, DUPOLD, DUTST, DUUP
      DOUBLE PRECISION   E1, E2, EIDN, EIUP, FL, U, UDN, UPOLD, UTST
      DOUBLE PRECISION   UUP
      CHARACTER*67       ERRMSG
!
!        External Subroutines
!
      EXTERNAL           OUTERR, OUTWIN
!
!        Intrinsic Functions
!
      INTRINSIC          MAX
!
      E1 = -100.0D+0
      E2 = -100.0D+0
      if (de.lt.1.d-8) goto 600
      FL = L
      EI = EI*0.5D+0
      DE = DE*0.5D+0
      CALL DIRACOUT(REL,VR,RNOT,DX,JRI,EI,L,1,U,DU,NODE,ZZ)
      EIUP = EI
      EIDN = EI
      UPOLD = U
      DUPOLD = DU
      DNOLD = U
      DDNOLD = DU
   10 CONTINUE
      EIUP = EIUP + DE
      CALL DIRACOUT(REL,VR,RNOT,DX,JRI,EIUP,L,1,UUP,DUUP,NODE,ZZ)
      UTST = UPOLD*UUP
      DUTST = DUPOLD*DUUP
!      WRITE(6,99)L,EIUP,UTST,DUTST
!   99 FORMAT(' L,EIUP,UTST,DUTST',I2,F8.4,2E12.3)
      UPOLD = UUP
      DUPOLD = DUUP
      IF (UTST .LT. 0.0D+0) THEN
         E2 = EIUP
         IF (E1 .GT. -30.0D+0) GOTO 30
      ENDIF
      IF (DUTST .LT. 0.0D+0) THEN
         E1 = EIUP
         IF (E2 .GT. -30.0D+0) GOTO 30
      ENDIF
      IF ((E1 .LT. -30.0D+0) .AND. (EIUP .LT. 2.5D+0)) THEN
   20    CONTINUE
         EIDN = EIDN - DE
         CALL DIRACOUT(REL,VR,RNOT,DX,JRI,EIDN,L,1,UDN,DUDN,NODE,ZZ)
         UTST = DNOLD*UDN
         DUTST = DDNOLD*DUDN
!         WRITE(6,98)L,EIDN,UTST,DUTST
!   98    FORMAT(' L,EIDN,UTST,DUTST',I2,F8.4,2E12.3)
         DNOLD = UDN
         DDNOLD = DUDN
         IF (UTST .LT. 0.0D+0) THEN
            E2 = EIDN
            IF (E1 .GT. -30.0D+0) GOTO 30
         ENDIF
         IF (DUTST .LT. 0.) THEN
            E1 = EIDN
            IF (E2 .GT. -30.0D+0) GOTO 30
         ENDIF
         IF (E2 .LT. -30.0D+0) THEN
            GOTO 10
         ELSEIF (EIDN .GT. (EI - 3.0D+0)) THEN
            GOTO 20
         ENDIF
      ELSEIF (EIUP .LT. 0.5D+0) THEN
         GOTO 10
      ENDIF
   30 CONTINUE
      EI = EI*2.0D+0
      E1 = E1*2.0D+0
      E2 = E2*2.0D+0
      IF ((E1 .LT. -60.0D+0) .AND. (E2 .LT. -60.0D+0)) THEN
         IF (EMAIN .EQ. 'STOP') THEN
!
!        Error: no energy limits found
!
            GOTO 900
         ELSE
            EI = 2.0D+0
!           WRITE (6,6000)  L, EI, E1, E2
!           WRITE (21,6000) L, EI, E1, E2
         ENDIF
      ELSEIF (E2 .LT. -60.0D+0) THEN
         IF (EMAIN .EQ. 'STOP') THEN
!
!        Error: no energy limits found
!
            GOTO 900
         ELSE
            EI = MAX(E1,EI)
!           WRITE (6,6000)  L, EI, E1, E2
!           WRITE (21,6000) L, EI, E1, E2
         ENDIF
      ELSEIF (E1 .LT. -60.0D+0) THEN
!
!        Error: no energy limits found
!
         GOTO 900
      ELSE
         EI = (E1+E2)*0.5D+0
!        WRITE (6,6000)  L, EI, E1, E2
!        WRITE (21,6000) L, EI, E1, E2
      ENDIF
!
      RETURN
!
  600 continue
      WRITE (21,6000) L, EI, E1, E2
!     WRITE (6,6000)  L, EI, E1, E2
      return
!        Error messages
!
  900 WRITE (ERRMSG,9000) L
      CALL OUTERR('SELECT',ERRMSG)
      WRITE (ERRMSG,9010) E1, E2
      CALL OUTERR('SELECT',ERRMSG)
      STOP 'SELECT - Error'
!
 6000 FORMAT (10X,'E(',I2,')=',F10.4,3X,'E(BOTTOM)=',F9.3,3X,'E(TOP)=', &
              F9.3)
 9000 FORMAT ('no energy limits found for L=',I2)
 9010 FORMAT ('E-bottom ',F10.5,3X,'E-top ',F10.5)
!
!        End of 'SELECT'
!
      END


text:
!
!  Written by Leandro Martínez, 2009-2011.
!  Copyright (c) 2009-2011, Leandro Martínez, Jose Mario Martinez,
!  Ernesto G. Birgin.
!

!
! Function that returns a real random number between 0. and 1.
!

double precision function rnd()

  call random_number(rnd)

  return
end function rnd

!
! Subroutine that initializes the random number generator given a seed
!

subroutine init_random_number(iseed)
  integer :: size
  integer :: i, iseed
  integer, allocatable :: seed(:)
  call random_seed(size=size)
  allocate(seed(size))
  do i = 1, size
    seed(i) = i*iseed
  end do
  call random_seed(put=seed)
  deallocate(seed)
  return
end subroutine init_random_number

!
! Subroutine that uses the date to create a random seed
!

subroutine seed_from_time(seed)

  implicit none
  integer :: seed, value(8)
  character(len=10) :: b(3)
  call date_and_time( b(1), b(2), b(3), value )
  seed = value(1)+value(2)+value(3)+value(4)+value(5)+value(6)+value(7)+value(8)
  seed = seed + value(1)+value(2)+value(3)+value(4)+value(5)/100+value(6)*100+value(7)/10+value(8)*10

end subroutine seed_from_time

max_stars_repo_name: susilehtola/packmol


text:
#include "eiscor.h"
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! d_1Darray_check
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This routine checks a one dimensional double array for INFs and NANs.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! INPUT VARIABLES:
!
!  N               INTEGER
!                    positive integer, dimension of A
!
!  A               REAL(8) array of dimension (N)
!                    real array to be checked
!
! OUTPUT VARIABLES:
!
!  FLAG            LOGICAL
!                    .TRUE. implies A contains valid numbers
!                    .FALSE. implies N < 1 or A contains at least one invalid
!                            number
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine d_1Darray_check(N,A,FLAG)

  implicit none

  ! input variables
  integer, intent(in) :: N
  real(8), intent(in) :: A(N)
  logical, intent(inout) :: FLAG

  ! compute variables
  integer :: ii

  ! initialize FLAG
  FLAG = .FALSE.

  ! check N
  if (N < 1) then
    return
  end if

  ! check array
  do ii = 1,N

    ! check for INF or NAN
    call d_scalar_check(A(ii),FLAG)
    if (.NOT.FLAG) then
      return
    end if

  end do

end subroutine d_1Darray_check

max_stars_repo_name: trcameron/FPML

max_issues_repo_name: trcameron/FPML

max_forks_repo_name: trcameron/FPML


text:
*
       SUBROUTINE GMRESREVCOM(N, B, X, RESTRT, WORK, LDW, WORK2,
     $                  LDW2, ITER, RESID, INFO, NDX1, NDX2, SCLR1,
     $                  SCLR2, IJOB)
*
*  -- Iterative template routine --
*     Univ. of Tennessee and Oak Ridge National Laboratory
*     October 1, 1993
*     Details of this algorithm are described in "Templates for the
*     Solution of Linear Systems: Building Blocks for Iterative
*     Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
*     EiITERkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
*     1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
*
*     .. Scalar Arguments ..
      INTEGER            N, RESTRT, LDW, LDW2, ITER, INFO
      DOUBLE PRECISION   RESID
      INTEGER            NDX1, NDX2
      DOUBLE PRECISION   SCLR1, SCLR2
      INTEGER            IJOB
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   B( * ), X( * ), WORK( LDW,* ), WORK2( LDW2,* )
*     ..
*
*  Purpose
*  =======
*
*  GMRES solves the linear system Ax = b using the
*  Generalized Minimal Residual iterative method with preconditioning.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER.
*          On entry, the dimension of the matrix.
*          Unchanged on exit.
*
*  B       (input) DOUBLE PRECISION array, dimension N.
*          On entry, right hand side vector B.
*          Unchanged on exit.
*
*  X       (input/output) DOUBLE PRECISION array, dimension N.
*          On input, the initial guess; on exit, the iterated solution.
*
*  RESTRT  (input) INTEGER
*          Restart parameter, <= N. This parameter controls the amount
*          of memory required for matrix WORK2.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (LDW,5).
*          Note that if the initial guess is the zero vector, then
*          storing the initial residual is not necessary.
*
*  LDW     (input) INTEGER
*          The leading dimension of the array WORK. LDW >= max(1,N).
*
*  WORK2   (workspace) DOUBLE PRECISION array, dimension (LDW2,2*RESTRT+2).
*          This workspace is used for constructing and storing the
*          upper Hessenberg matrix. The two extra columns are used to
*          store the Givens rotation matrices.
*
*  LDW2    (input) INTEGER
*          The leading dimension of the array WORK2.
*          LDW2 >= max(1,RESTRT).
*
*  ITER    (input/output) INTEGER
*          On input, the maximum iterations to be performed.
*          On output, actual number of iterations performed.
*
*  RESID   (input/output) DOUBLE PRECISION
*          On input, the allowable error tolerance.
*          On ouput, the norm of the residual vector if solution
*          approximated to tolerance, otherwise reset to input
*          tolerance.
*
*  INFO    (output) INTEGER
*          =  0:  successful exit
*          =  1:  maximum number of iterations performed;
*                 convergence not achieved.
*            -5: Erroneous NDX1/NDX2 in INIT call.
*            -6: Erroneous RLBL.
*
*  NDX1    (input/output) INTEGER.
*  NDX2    On entry in INIT call contain indices required by interface
*          level for stopping test.
*          All other times, used as output, to indicate indices into
*          WORK[] for the MATVEC, PSOLVE done by the interface level.
*
*  SCLR1   (output) DOUBLE PRECISION.
*  SCLR2   Used to pass the scalars used in MATVEC. Scalars are reqd because
*          original routines use dgemv.
*
*  IJOB    (input/output) INTEGER.
*          Used to communicate job code between the two levels.
*
*  ============================================================
*
*     .. Parameters ..
      DOUBLE PRECISION    ZERO, ONE
      PARAMETER         ( ZERO = 0.0D+0, ONE = 1.0D+0 )
      INTEGER             OFSET
      PARAMETER         ( OFSET = 1000 )
*     ..
*     .. Local Scalars ..
      INTEGER             I, MAXIT, AV, GIV, H, R, S, V, W, Y,
     $                    NEED1, NEED2
      DOUBLE PRECISION    BNRM2, RNORM, TOL, APPROXRES, DDOT, DNRM2
*
*     indicates where to resume from. Only valid when IJOB = 2!
      INTEGER RLBL
*
*     saving all.
      SAVE
*
*     ..
*     .. External Routines ..
      EXTERNAL            DAXPY, DCOPY, DDOT, DNRM2
*     ..
*     .. Executable Statements ..
*
* Entry point, so test IJOB
      IF (IJOB .eq. 1) THEN
         GOTO 1
      ELSEIF (IJOB .eq. 2) THEN
*        here we do resumption handling
         IF (RLBL .eq. 2) GOTO 2
         IF (RLBL .eq. 3) GOTO 3
         IF (RLBL .eq. 4) GOTO 4
         IF (RLBL .eq. 5) GOTO 5
         IF (RLBL .eq. 6) GOTO 6
         IF (RLBL .eq. 7) GOTO 7
*        if neither of these, then error
         INFO = -6
         GOTO 200
      ENDIF
*
* init.
*****************
 1    CONTINUE
*****************
*
      INFO = 0
      MAXIT = ITER
      TOL   = RESID
*
*     Alias workspace columns.
*
      R   = 1
      S   = 2
      W   = 3
      Y   = 4
      AV  = 5
      V   = 6
*
      H   = 1
      GIV = H + RESTRT
*
*     Check if caller will need indexing info.
*
      IF( NDX1.NE.-1 ) THEN
         IF( NDX1.EQ.1 ) THEN
            NEED1 = ((R - 1) * LDW) + 1
         ELSEIF( NDX1.EQ.2 ) THEN
            NEED1 = ((S - 1) * LDW) + 1
         ELSEIF( NDX1.EQ.3 ) THEN
            NEED1 = ((W - 1) * LDW) + 1
         ELSEIF( NDX1.EQ.4 ) THEN
            NEED1 = ((Y - 1) * LDW) + 1
         ELSEIF( NDX1.EQ.5 ) THEN
            NEED1 = ((AV - 1) * LDW) + 1
         ELSEIF( NDX1.EQ.6 ) THEN
            NEED1 = ((V - 1) * LDW) + 1
         ELSEIF( ( NDX1.GT.V*OFSET ) .AND.
     $           ( NDX1.LE.V*OFSET+RESTRT ) ) THEN
            NEED1 = ((NDX1-V*OFSET - 1) * LDW) + 1
         ELSEIF( ( NDX1.GT.GIV*OFSET ) .AND.
     $           ( NDX1.LE.GIV*OFSET+RESTRT ) ) THEN
            NEED1 = ((NDX1-GIV*OFSET - 1) * LDW) + 1
         ELSE
*           report error
            INFO = -5
            GO TO 100
         ENDIF
      ELSE
         NEED1 = NDX1
      ENDIF
*
      IF( NDX2.NE.-1 ) THEN
         IF( NDX2.EQ.1 ) THEN
            NEED2 = ((R - 1) * LDW) + 1
         ELSEIF( NDX2.EQ.2 ) THEN
            NEED2 = ((S - 1) * LDW) + 1
         ELSEIF( NDX2.EQ.3 ) THEN
            NEED2 = ((W - 1) * LDW) + 1
         ELSEIF( NDX2.EQ.4 ) THEN
            NEED2 = ((Y - 1) * LDW) + 1
         ELSEIF( NDX2.EQ.5 ) THEN
            NEED2 = ((AV - 1) * LDW) + 1
         ELSEIF( NDX2.EQ.6 ) THEN
            NEED2 = ((V - 1) * LDW) + 1
         ELSEIF( ( NDX2.GT.V*OFSET ) .AND.
     $           ( NDX2.LE.V*OFSET+RESTRT ) ) THEN
            NEED2 = ((NDX2-V*OFSET - 1) * LDW) + 1
         ELSEIF( ( NDX2.GT.GIV*OFSET ) .AND.
     $           ( NDX2.LE.GIV*OFSET+RESTRT ) ) THEN
            NEED2 = ((NDX2-GIV*OFSET - 1) * LDW) + 1
         ELSE
*           report error
            INFO = -5
            GO TO 100
         ENDIF
      ELSE
         NEED2 = NDX2
      ENDIF
*
*     Set initial residual.
*
      CALL DCOPY( N, B, 1, WORK(1,R), 1 )
      IF ( DNRM2( N, X, 1 ).NE.ZERO ) THEN
*********CALL MATVEC( -ONE, X, ONE, WORK(1,R) )
*        Note: using X directly
         SCLR1 = -ONE
         SCLR2 = ONE
         NDX1 = -1
         NDX2 = ((R - 1) * LDW) + 1
*
*        Prepare for resumption & return
         RLBL = 2
         IJOB = 1
         RETURN
*
*****************
 2       CONTINUE
*****************
*
         IF ( DNRM2( N, WORK(1,R), 1 ).LT.TOL ) GO TO 200
      ENDIF
      BNRM2 = DNRM2( N, B, 1 )
      IF ( BNRM2.EQ.ZERO ) BNRM2 = ONE
*
      ITER = 0
   10 CONTINUE
*
         ITER = ITER + 1
*
*        Construct the first column of V, and initialize S to the
*        elementary vector E1 scaled by RNORM.
*
*********CALL PSOLVE( WORK( 1,V ), WORK( 1,R ) )
*
         NDX1 = ((V - 1) * LDW) + 1
         NDX2 = ((R - 1) * LDW) + 1
*        Prepare for return & return
         RLBL = 3
         IJOB = 2
         RETURN
*
*****************
 3       CONTINUE
*****************
*
         RNORM = DNRM2( N, WORK( 1,V ), 1 )
         CALL DSCAL( N, RNORM, WORK( 1,V ), 1 )
         CALL ELEMVEC( 1, N, RNORM, WORK( 1,S ) )
*
         DO 50 I = 1, RESTRT
************CALL MATVEC( ONE, WORK( 1,V+I-1 ), ZERO, WORK( 1,AV ) )
*
         NDX1 = ((V+I-1 - 1) * LDW) + 1
         NDX2 = ((AV    - 1) * LDW) + 1
*        Prepare for return & return
         SCLR1 = ONE
         SCLR2 = ZERO
         RLBL = 4
         IJOB = 3
         RETURN
*
*****************
 4       CONTINUE
*****************
*
*********CALL PSOLVE( WORK( 1,W ), WORK( 1,AV ) )
*
         NDX1 = ((W  - 1) * LDW) + 1
         NDX2 = ((AV - 1) * LDW) + 1
*        Prepare for return & return
         RLBL = 5
         IJOB = 2
         RETURN
*
*****************
 5       CONTINUE
*****************
*
*           Construct I-th column of H so that it is orthnormal to
*           the previous I-1 columns.
*
            CALL ORTHOH( I, N, WORK2( 1,I+H-1 ), WORK( 1,V ), LDW,
     $                   WORK( 1,W ) )
*
            IF ( I.GT.1 )
*
*              Apply Givens rotations to the I-th column of H. This
*              effectively reduces the Hessenberg matrix to upper
*              triangular form during the RESTRT iterations.
*
     $         CALL APPLYGIVENS( I, WORK2( 1,I+H-1 ), WORK2( 1,GIV ),
     $                           LDW2 )
*
*           Approxiyate residual norm. Check tolerance. If okay, compute
*           final approximation vector X and quit.
*
            RESID = APPROXRES( I, WORK2( 1,I+H-1 ), WORK( 1,S ),
     $                         WORK2( 1,GIV ), LDW2 ) / BNRM2
            IF ( RESID.LE.TOL ) THEN
               CALL UPDATE(I, N, X, WORK2( 1,H ), LDW2,
     $                     WORK(1,Y), WORK(1,S), WORK( 1,V ), LDW)
               GO TO 200
            ENDIF
   50    CONTINUE
*
*        Compute current solution vector X.
*
         CALL UPDATE(RESTRT, N, X, WORK2( 1,H ), LDW2,
     $               WORK(1,Y), WORK( 1,S ), WORK( 1,V ), LDW )
*
*        Compute residual vector R, find norm,
*        then check for tolerance.
*
         CALL DCOPY( N, B, 1, WORK( 1,R ), 1 )
*********CALL MATVEC( -ONE, X, ONE, WORK( 1,R ) )
*
         NDX1 = -1
         NDX2 = ((R - 1) * LDW) + 1
*        Prepare for return & return
         SCLR1 = -ONE
         SCLR2 = ONE
         RLBL = 6
         IJOB = 1
         RETURN
*
*****************
 6       CONTINUE
*****************
*
         WORK( I+1,S ) = DNRM2( N, WORK( 1,R ), 1 )
*
*********RESID = WORK( I+1,S ) / BNRM2
*********IF ( RESID.LE.TOL  ) GO TO 200
*
         NDX1 = NEED1
         NDX2 = NEED2
*        Prepare for resumption & return
         RLBL = 7
         IJOB = 4
         RETURN
*
*****************
 7       CONTINUE
*****************
         IF( INFO.EQ.1 ) GO TO 200
*
         IF ( ITER.EQ.MAXIT ) THEN
            INFO = 1
            GO TO 100
         ENDIF
*
         GO TO 10
*
  100 CONTINUE
*
*     Iteration fails.
*
      RLBL = -1
      IJOB = -1
      RETURN
*
  200 CONTINUE
*
*     Iteration successful; return.
*
      INFO = 0
      RLBL = -1
      IJOB = -1

      RETURN
*
*     End of GMRESREVCOM
*
      END
*
*     =========================================================
      SUBROUTINE ORTHOH( I, N, H, V, LDV, W )
*
      INTEGER            I, N, LDV
      DOUBLE PRECISION   H( * ), W( * ), V( LDV,* )
*
*     Construct the I-th column of the upper Hessenberg matrix H
*     using the Gram-Schmidt process on V and W.
*
      INTEGER            K
      DOUBLE PRECISION   DDOT, DNRM2, ONE
      PARAMETER        ( ONE = 1.0D+0 )
      EXTERNAL           DAXPY, DCOPY, DDOT, DNRM2, DSCAL
*
      DO 10 K = 1, I
         H( K ) = DDOT( N, W, 1, V( 1,K ), 1 )
         CALL DAXPY( N, -H( K ), V( 1,K ), 1, W, 1 )
   10 CONTINUE
      H( I+1 ) = DNRM2( N, W, 1 )
      CALL DCOPY( N, W, 1, V( 1,I+1 ), 1 )
      CALL DSCAL( N, ONE / H( I+1 ), V( 1,I+1 ), 1 )
*
      RETURN
*
      END
*     =========================================================
      SUBROUTINE APPLYGIVENS( I, H, GIVENS, LDG )
*
      INTEGER            I, LDG
      DOUBLE PRECISION   H( * ), GIVENS( LDG,* )
*
*     This routine applies a sequence of I-1 Givens rotations to
*     the I-th column of H. The Givens parameters are stored, so that
*     the first I-2 Givens rotation matrices are known. The I-1st
*     Givens rotation is computed using BLAS 1 routine DROTG. Each
*     rotation is applied to the 2x1 vector [H( J ), H( J+1 )]',
*     which results in H( J+1 ) = 0.
*
      INTEGER            J
*      DOUBLE PRECISION   TEMP
      EXTERNAL           DROTG
*
*     .. Executable Statements ..
*
*     Construct I-1st rotation matrix.
*
*     CALL DROTG( H( I ), H( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
      CALL GETGIV( H( I ), H( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
*
*     Apply 1,...,I-1st rotation matrices to the I-th column of H.
*
      DO 10 J = 1, I-1
         CALL ROTVEC( H( J ), H( J+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
*        TEMP     =  GIVENS( J,1 ) * H( J ) + GIVENS( J,2 ) * H( J+1 )
*        H( J+1 ) = -GIVENS( J,2 ) * H( J ) + GIVENS( J,1 ) * H( J+1 )
*        H( J ) = TEMP
 10   CONTINUE
*
      RETURN
*
      END
*
*     ===============================================================
      DOUBLE PRECISION FUNCTION APPROXRES( I, H, S, GIVENS, LDG )
*
      INTEGER            I, LDG
      DOUBLE PRECISION   H( * ), S( * ), GIVENS( LDG,* )
*
*     This function allows the user to approximate the residual
*     using an updating scheme involving Givens rotations. The
*     rotation matrix is formed using [H( I ),H( I+1 )]' with the
*     intent of zeroing H( I+1 ), but here is applied to the 2x1
*     vector [S(I), S(I+1)]'.
*
      INTRINSIC          DABS
      EXTERNAL           DROTG
*
*     .. Executable Statements ..
*
*     CALL DROTG( H( I ), H( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
      CALL GETGIV( H( I ), H( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
      CALL ROTVEC( S( I ), S( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
*
      APPROXRES = DABS( S( I+1 ) )
*
      RETURN
*
      END
*     ===============================================================
      SUBROUTINE UPDATE( I, N, X, H, LDH, Y, S, V, LDV )
*
      INTEGER            N, I, J, LDH, LDV
      DOUBLE PRECISION   X( * ), Y( * ), S( * ), H( LDH,* ), V( LDV,* )
      EXTERNAL           DAXPY, DCOPY, DTRSV
*
*     Solve H*y = s for upper triangualar H.
*
      CALL DCOPY( I, S, 1, Y, 1 )
      CALL DTRSV( 'UPPER', 'NOTRANS', 'NONUNIT', I, H, LDH, Y, 1 )
*
*     Compute current solution vector X.
*
      DO 10 J = 1, I
         CALL DAXPY( N, Y( J ), V( 1,J ), 1, X, 1 )
   10 CONTINUE
*
      RETURN
*
      END
*
*     ===============================================================
      SUBROUTINE GETGIV( A, B, C, S )
*
      DOUBLE PRECISION   A, B, C, S, TEMP, ZERO, ONE
      PARAMETER        ( ZERO = 0.0D+0, ONE = 1.0D+0 )
*
      IF ( B.EQ.ZERO ) THEN
         C = ONE
         S = ZERO
      ELSE IF ( ABS( B ).GT.ABS( A ) ) THEN
         TEMP = -A / B
         S = ONE / SQRT( ONE + TEMP**2 )
         C = TEMP * S
      ELSE
         TEMP = -B / A
         C = ONE / SQRT( ONE + TEMP**2 )
         S = TEMP * C
      ENDIF
*
      RETURN
*
      END
*
*     ================================================================
      SUBROUTINE ROTVEC( X, Y, C, S )
*
      DOUBLE PRECISION   X, Y, C, S, TEMP

*
      TEMP = C * X - S * Y
      Y    = S * X + C * Y
      X    = TEMP
*
      RETURN
*
      END
*
*     ===============================================================
      SUBROUTINE ELEMVEC( I, N, ALPHA, E )
*
*     Construct the I-th elementary vector E, scaled by ALPHA.
*
      INTEGER            I, J, N
      DOUBLE PRECISION   ALPHA, E( * )
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      PARAMETER        ( ZERO = 0.0D+0 )
*
      DO 10 J = 1, N
         E( J ) = ZERO
   10 CONTINUE
      E( I ) = ALPHA
*
      RETURN
*
      END

ext: f

max_stars_repo_name: ljktest/siconos

max_issues_repo_name: ljktest/siconos

max_forks_repo_name: ljktest/siconos


! { dg-do run }
  real :: a(30), m
  real, allocatable :: c(:)
  integer :: e(30), n, ia(1)
  integer, allocatable :: g(:)
  logical :: l(30)
  allocate (c (30))
  allocate (g (30))
  a = 7.0
  c = 7.0
  e = 7
  g = 7
  m = huge(m)
  n = huge(n)
  a(7) = 6.0
  c(7) = 6.0
  e(7) = 6
  g(7) = 6
  ia = minloc (a)
  if (ia(1).ne.7) call abort
  ia = minloc (a(::2))
  if (ia(1).ne.4) call abort
  if (any (minloc (a).ne.(/ 7 /))) call abort
  if (any (minloc (a(::2)).ne.(/ 4 /))) call abort
  ia = minloc (c)
  if (ia(1).ne.7) call abort
  ia = minloc (c(::2))
  if (ia(1).ne.4) call abort
  if (any (minloc (c).ne.(/ 7 /))) call abort
  if (any (minloc (c(::2)).ne.(/ 4 /))) call abort
  ia = minloc (e)
  if (ia(1).ne.7) call abort
  ia = minloc (e(::2))
  if (ia(1).ne.4) call abort
  if (any (minloc (e).ne.(/ 7 /))) call abort
  if (any (minloc (e(::2)).ne.(/ 4 /))) call abort
  ia = minloc (g)
  if (ia(1).ne.7) call abort
  ia = minloc (g(::2))
  if (ia(1).ne.4) call abort
  if (any (minloc (g).ne.(/ 7 /))) call abort
  if (any (minloc (g(::2)).ne.(/ 4 /))) call abort
  l = .true.
  ia = minloc (a, mask = l)
  if (ia(1).ne.7) call abort
  ia = minloc (a(::2), mask = l(::2))
  if (ia(1).ne.4) call abort
  if (any (minloc (a, mask = l).ne.(/ 7 /))) call abort
  if (any (minloc (a(::2), mask = l(::2)).ne.(/ 4 /))) call abort
  ia = minloc (c, mask = l)
  if (ia(1).ne.7) call abort
  ia = minloc (c(::2), mask = l(::2))
  if (ia(1).ne.4) call abort
  if (any (minloc (c, mask = l).ne.(/ 7 /))) call abort
  if (any (minloc (c(::2), mask = l(::2)).ne.(/ 4 /))) call abort
  ia = minloc (e, mask = l)
  if (ia(1).ne.7) call abort
  ia = minloc (e(::2), mask = l(::2))
  if (ia(1).ne.4) call abort
  if (any (minloc (e, mask = l).ne.(/ 7 /))) call abort
  if (any (minloc (e(::2), mask = l(::2)).ne.(/ 4 /))) call abort
  ia = minloc (g, mask = l)
  if (ia(1).ne.7) call abort
  ia = minloc (g(::2), mask = l(::2))
  if (ia(1).ne.4) call abort
  if (any (minloc (g, mask = l).ne.(/ 7 /))) call abort
  if (any (minloc (g(::2), mask = l(::2)).ne.(/ 4 /))) call abort
  l = .false.
  ia = minloc (a, mask = l)
  if (ia(1).ne.0) call abort
  ia = minloc (a(::2), mask = l(::2))
  if (ia(1).ne.0) call abort
  if (any (minloc (a, mask = l).ne.(/ 0 /))) call abort
  if (any (minloc (a(::2), mask = l(::2)).ne.(/ 0 /))) call abort
  ia = minloc (c, mask = l)
  if (ia(1).ne.0) call abort
  ia = minloc (c(::2), mask = l(::2))
  if (ia(1).ne.0) call abort
  if (any (minloc (c, mask = l).ne.(/ 0 /))) call abort
  if (any (minloc (c(::2), mask = l(::2)).ne.(/ 0 /))) call abort
  ia = minloc (e, mask = l)
  if (ia(1).ne.0) call abort
  ia = minloc (e(::2), mask = l(::2))
  if (ia(1).ne.0) call abort
  if (any (minloc (e, mask = l).ne.(/ 0 /))) call abort
  if (any (minloc (e(::2), mask = l(::2)).ne.(/ 0 /))) call abort
  ia = minloc (g, mask = l)
  if (ia(1).ne.0) call abort
  ia = minloc (g(::2), mask = l(::2))
  if (ia(1).ne.0) call abort
  if (any (minloc (g, mask = l).ne.(/ 0 /))) call abort
  if (any (minloc (g(::2), mask = l(::2)).ne.(/ 0 /))) call abort
  a = 7.0
  c = 7.0
end


program dumper

   use ieee_arithmetic
   implicit none

   character(len=999) :: fname
   integer, external :: iargc
   character(len=999) :: buf
   character(len=16) dataset, dataver
   character(len=64) datamod
   character(len=132) comment
   character(len=8) axtz,pmap,datum
   character(len=12) daten
   character(len=4) utmhem
   character(len=8) clab
   character(len=8), allocatable, dimension(:) :: clabs

   character(len=16), parameter :: prgname = 'calmet_dump'

   integer:: idh0, isec0, idh1, isec1
   integer :: novar
   character(len=8), dimension(100) :: ovar

   integer :: idbg, jdbg,kdbg
   integer :: nclab, ncheader
   logical :: lnoisy

   integer :: ncom
   integer ::ibyr,ibmo,ibdy,ibhr,ibsec,ieyr,iemo,iedy,iehr,iesec,irlg,irtype,nx,ny,nz,&
           nssta,nusta,npsta,nowsta,nlu,iwat1,iwat2,iutmzn
   real :: dgrid, xorigr,yorigr,feast,fnorth,rnlat0,relon0,xlat1,xlat2
   logical :: lcalgrd
   real, allocatable, dimension(:) :: zfacem,xssta,yssta,xusta,yusta,xpsta,ypsta
   real, allocatable, dimension(:,:) :: z0,elev,xlai,ustar,zi,el,wstar,rmm,tempk,rho,qsw
   integer, allocatable, dimension(:,:) :: ilandu,nears,ipgt,irh,ipcode
   real, allocatable, dimension(:,:,:) :: uu,vv,ww,tt

   integer,dimension(99) :: idum
   integer :: n,i,k,it,ios

   ! get file name

   lnoisy=.true.
   idbg=1
   jdbg=1
   kdbg=1
   idh0= 190000000
   idh1= 209936624
   isec0=0
   isec1=3600
   novar = -99
   !n = iargc()
   n = command_argument_count()
   if (n == 1) then
      !call getarg(1, fname)
      call get_command_argument(1, fname)
   elseif (n == 3) then
      !call getarg(1, fname)
      call get_command_argument(1, fname)
      read(fname, *, iostat=ios) idbg
      if (ios /= 0) then
         print*,'arg1 parse error: ', trim(fname)
         call usage
      endif
      !call getarg(2, fname)
      call get_command_argument(2, fname)
      read(fname, *, iostat=ios) jdbg
      if (ios /= 0) then
         print*,'arg2 parse error: ', trim(fname)
         call usage
      endif
      !call getarg(3, fname)
      call get_command_argument(3, fname)
   else
      call usage
   endif
   if (trim(fname) == '-') then
      ! read input from stdin
      call userin
   endif
   if (lnoisy) print*,trim(fname)
   open(11,file=fname,form='unformatted', status='old')

   ! scan entire file first, read for clabs
   call scn

   !print*,nclab,ncheader
   !print*,clabs(1:20)

   allocate(z0(nx,ny),ilandu(nx,ny),elev(nx,ny),xlai(nx,ny),nears(nx,ny))
   allocate(uu(nx,ny,nz),vv(nx,ny,nz),ww(nx,ny,nz),tt(nx,ny,nz))
   allocate(ipgt(nx,ny),ustar(nx,ny),zi(nx,ny),el(nx,ny),wstar(nx,ny),rmm(nx,ny), &
      tempk(nx,ny),rho(nx,ny),qsw(nx,ny),irh(nx,ny),ipcode(nx,ny))
   rewind(11)

   call dmp

contains

   subroutine scn
      integer :: ic, ios, i
      ! header, metadata part
      read(11)dataset, dataver, datamod
      read(11)ncom
      if (lnoisy) print*,ncom
      if (ncom <= 10) then
         do i=1,ncom
            read(11) comment
            if (lnoisy) print*,comment
         enddo
      else
         do i=1,ncom
            read(11) comment
            if (lnoisy) then
               if (i <= 6 .or. i >=(ncom-2)) then
                  print*,comment
               elseif (i==7) then
                       print*
                       print*,'   ...   '
                       print*
               endif
            endif
         enddo
      endif

      ! header, scalar part
      read(11)ibyr,ibmo,ibdy,ibhr,ibsec,ieyr,iemo,iedy,iehr,iesec,axtz,&
              irlg,irtype,nx,ny,nz,dgrid,xorigr,yorigr,&
              nssta,nusta,npsta,nowsta,nlu,iwat1,iwat2,lcalgrd,pmap,datum,daten,&
              feast,fnorth,utmhem,iutmzn,rnlat0,relon0,xlat1,xlat2

      if (lnoisy) then
         print*
         print*,'i,jdbg  :', idbg,jdbg
         print*,'begin   :', ibyr,ibmo,ibdy,ibhr,ibsec
         print*,'end     :', ieyr,iemo,iedy,iehr,iesec
         print*,'irlg    :', irlg
         print*,'irtype  :', irlg
         print*,'nx/ny/nz:', nx,ny,nz
         print*,'xor/yor/dx:', xorigr,yorigr,dgrid
         print*
      endif

      ! i dont know why i cannot preduct # of records...
      !allocate(clabs(nx*ny*irlg*4+100))
      !allocate(clabs(99999))
      allocate(clabs(999999))
      !print*,nx*ny*irlg*4+100

      nclab = 0
      do while(.true.)
         read(11,iostat=ios)clab,idum(1:5)
         if (ios/=0) exit
         !print*,clab,idum(1:5)

         nclab = nclab + 1
         clabs(nclab) = clab
      enddo
      !print*,nclab

      do ic=1,nclab
         if (clabs(ic)(1:5)=='U-LEV') then
            ncheader = ic - 1
            exit
         endif
      enddo

   end subroutine

   subroutine dmp
      integer :: ic,ityp, ilab
      character(len=80) :: fmtf, fmti, fmta
      read(11)dataset, dataver, datamod
      read(11)ncom
      do i=1,ncom
         read(11) comment
      enddo

      ! header, scalar part
      read(11)ibyr,ibmo,ibdy,ibhr,ibsec,ieyr,iemo,iedy,iehr,iesec,axtz,&
              irlg,irtype,nx,ny,nz,dgrid,xorigr,yorigr,&
              nssta,nusta,npsta,nowsta,nlu,iwat1,iwat2,lcalgrd,pmap,datum,daten,&
              feast,fnorth,utmhem,iutmzn,rnlat0,relon0,xlat1,xlat2

      fmtf = "(a8,',',2(i9,',',i4,','),3(i4,','),e14.7)"
      fmti = "(a8,',',2(i9,',',i4,','),3(i4,','),i14)"
      fmta = "(a8,',',2(a9,',',a4,','),3(a4,','),a14)"
      ! data
      !print*,ncheader
      if (.not.lnoisy) write(*,fmta) 'variable', 'bdh', 'bsec', 'edh', 'esec','i','j','k', 'value'
      do ic=1,ncheader
         clab=clabs(ic)
         !print*,ic,clab
         if (clab == 'ZFACE') then
            read(11) clab, idum(1:4), zfacem
            if(lnoisy) print*, clab,idum(1:4),zfacem
         elseif (clab == 'XSSTA') then
            read(11)clab,idum(1:4),xssta
            if(lnoisy) print*, clab,idum(1:4)
         elseif (clab == 'YSSTA') then
            read(11)clab,idum(1:4),yssta
            if(lnoisy) print*, clab,idum(1:4)
         elseif (clab == 'XUSTA') then
            read(11)clab,idum(1:4),xusta
            if(lnoisy) print*, clab,idum(1:4)
         elseif (clab == 'YUSTA') then
            read(11)clab,idum(1:4),yusta
            if(lnoisy) print*, clab,idum(1:4)
         elseif (clab == 'XPSTA') then
            read(11)clab,idum(1:4),xpsta
            if(lnoisy) print*, clab,idum(1:4)
         elseif (clab == 'YPSTA') then
            read(11)clab,idum(1:4),ypsta
            if(lnoisy) print*, clab,idum(1:4)
         elseif (clab == 'Z0') then
            read(11)clab,idum(1:4),z0
            if(lnoisy) print*,clab, idum(1:4),z0(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               write(*,fmtf) clab, 0, 0, 0, 0,idbg,jdbg,0, z0(idbg,jdbg)
            endif
         elseif (clab == 'ILANDU') then
            read(11)clab,idum(1:4),ilandu
            if(lnoisy) print*,clab, idum(1:4),ilandu(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               write(*,fmti) clab, 0, 0, 0, 0,idbg,jdbg,0, ilandu(idbg,jdbg)
            endif
         elseif (clab == 'ELEV') then
            read(11)clab,idum(1:4),elev
            if(lnoisy) print*,clab, idum(1:4),elev(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               write(*,fmtf) clab, 0, 0, 0, 0,idbg,jdbg,0, elev(idbg,jdbg)
            endif
         elseif (clab == 'XLAI') then
            read(11)clab,idum(1:4),xlai
            if(lnoisy) print*,clab, idum(1:4),xlai(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               write(*,fmtf) clab, 0, 0, 0, 0,idbg,jdbg,0, xlai(idbg,jdbg)
            endif
         elseif (clab == 'NEARS') then
            read(11)clab,idum(1:4),nears
            if(lnoisy) print*,clab, idum(1:4),nears(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               write(*,fmti) clab, 0, 0, 0, 0,idbg,jdbg,0, nears(idbg,jdbg)
            endif
         else
            print*, 'unknown clab (header):', clab
            stop
         endif

      enddo
      do ic=ncheader+1, nclab
         clab=clabs(ic)
         !print*,ic,clab
         if (clab(1:5) == 'U-LEV') then
            read(clab(6:8),*)k
            read(11)clab,idum(1:4),uu(:,:,k)
            if(lnoisy) print*, clab,idum(1:4),uu(idbg,jdbg,k)
            if (novar == 0 .or. any(clab == ovar)) then
               write(clab(6:8), '("   ")')
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,uu(idbg,jdbg,k)
               endif
            endif
         elseif (clab(1:5) == 'V-LEV') then
            read(clab(6:8),*)k
            read(11)clab,idum(1:4),vv(:,:,k)
            if(lnoisy) print*, clab,idum(1:4),vv(idbg,jdbg,k)
            if (novar == 0 .or. any(clab == ovar)) then
               write(clab(6:8), '("   ")')
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,vv(idbg,jdbg,k)
               endif
            endif
         elseif (clab(1:5) == 'WFACE') then
            read(clab(6:8),*)k
            read(11)clab,idum(1:4),ww(:,:,k)
            if(lnoisy) print*, clab,idum(1:4),ww(idbg,jdbg,k)
            if (novar == 0 .or. any(clab == ovar)) then
               write(clab(6:8), '("   ")')
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,ww(idbg,jdbg,k)
               endif
            endif
         elseif (clab(1:5) == 'T-LEV') then
            read(clab(6:8),*)k
            read(11)clab,idum(1:4),tt(:,:,k)
            if(lnoisy) print*, clab,idum(1:4),tt(idbg,jdbg,k)
            if (novar == 0 .or. any(clab == ovar)) then
               write(clab(6:8), '("   ")')
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,tt(idbg,jdbg,k)
               endif
            endif
         elseif (clab(1:5) == 'IPGT') then
            read(11)clab,idum(1:4),ipgt
            k = 1
            if(lnoisy) print*, clab,idum(1:4),ipgt(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmti) clab, idum(1:4), idbg,jdbg,k,ipgt(idbg,jdbg)
               endif
            endif
         elseif (clab(1:5) == 'USTAR') then
            read(11)clab,idum(1:4),ustar
            k = 1
            if(lnoisy) print*, clab,idum(1:4),ustar(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,ustar(idbg,jdbg)
               endif
            endif
         elseif (clab(1:5) == 'ZI') then
            read(11)clab,idum(1:4),zi
            k = 1
            if(lnoisy) print*, clab,idum(1:4),zi(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,zi(idbg,jdbg)
               endif
            endif
         elseif (clab(1:5) == 'EL') then
            read(11)clab,idum(1:4),el
            k = 1
            if(lnoisy) print*, clab,idum(1:4),el(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,el(idbg,jdbg)
               endif
            endif
         elseif (clab(1:5) == 'WSTAR') then
            read(11)clab,idum(1:4),wstar
            k = 1
            if(lnoisy) print*, clab,idum(1:4),wstar(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,wstar(idbg,jdbg)
               endif
            endif
         elseif (clab(1:5) == 'RMM') then
            read(11)clab,idum(1:4),rmm
            k = 1
            if(lnoisy) print*, clab,idum(1:4),rmm(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,rmm(idbg,jdbg)
               endif
            endif
         elseif (clab(1:5) == 'TEMPK') then
            read(11)clab,idum(1:4),tempk
            k = 1
            if(lnoisy) print*, clab,idum(1:4),tempk(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,tempk(idbg,jdbg)
               endif
            endif
         elseif (clab(1:5) == 'RHO') then
            read(11)clab,idum(1:4),rho
            k = 1
            if(lnoisy) print*, clab,idum(1:4),rho(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,rho(idbg,jdbg)
               endif
            endif
         elseif (clab(1:5) == 'QSW') then
            read(11)clab,idum(1:4),qsw
            k = 1
            if(lnoisy) print*, clab,idum(1:4),qsw(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmtf) clab, idum(1:4), idbg,jdbg,k,qsw(idbg,jdbg)
               endif
            endif
         elseif (clab(1:5) == 'IRH') then
            read(11)clab,idum(1:4),irh
            k = 1
            if(lnoisy) print*, clab,idum(1:4),irh(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmti) clab, idum(1:4), idbg,jdbg,k,irh(idbg,jdbg)
               endif
            endif
         elseif (clab(1:6) == 'IPCODE') then
            read(11)clab,idum(1:4),ipcode
            k = 1
            if(lnoisy) print*, clab,idum(1:4),ipcode(idbg,jdbg)
            if (novar == 0 .or. any(clab == ovar)) then
               if (k==kdbg .and. &
                  (idum(1) > idh0 .or. (idum(1) == idh0 .and. idum(2) >= isec0)) .and. &
                  (idum(1) < idh1 .or. (idum(1) == idh1 .and. idum(2) <= isec1)) ) then
                  write(*,fmti) clab, idum(1:4), idbg,jdbg,k,ipcode(idbg,jdbg)
               endif
            endif
         else
            print*, 'unknown clab (data):', ic, clab
            stop
         endif

      enddo

   end subroutine

   subroutine userin
      character(len=999) s
      lnoisy = .false.
      write(0,'(a)', advance='no')'fname: '
      read(*, '(a)') fname
      write(0, *) trim(fname)

      write(0,'(a)', advance='no')'i,j: '
      s = ''
      read(*, '(a)', iostat=ios) s
      if (ios /= 0 .or. s == '') then
         write(0, *) 'not specified, using default'
      else
         read(s, *, iostat=ios) idbg, jdbg
         if (ios /= 0 ) then
            write(0,*) 'cannot parse ijdx: ', trim(s)
            stop
         endif
      endif
      write(0, *) idbg, jdbg

      write(0,'(a)', advance='no') 'btime (yyyyjjjhh ssss): '
      s = ''
      read(*, '(a)', iostat=ios) s
      if (ios /= 0 .or. s == '') then
         write(0, *) 'not specified, using default'
      else
         read(s, *, iostat=ios) idh0, isec0
         if (ios /= 0 ) then
            write(0,*) 'cannot parse btime: ', trim(s)
            stop
         endif
         ! update default for endtime
         idh1 = idh0
         isec1 = isec0
      endif
      write(0,*) idh0, isec0

      write(0,'(a)', advance='no') 'etime (yyyyjjjhh ssss): '
      s = ''
      read(*, '(a)', iostat=ios) s
      if (ios /= 0 .or. s == '') then
         write(0, *) 'not specified, using default'
      else
         read(s, *, iostat=ios) idh1,isec1
         if (ios /= 0 ) then
            write(0,*) 'cannot parse etime: ', trim(s)
            stop
         endif
      endif
      write(0,*) idh1, isec1

      novar = 0
      write(0,'(a)', advance='no') 'variables (blank to all): '
      do while(.true.)
        read(*, '(a)', iostat=ios) s
        if (ios/=0) then
           s= ''
        endif
        if (s=='') then
           if (novar==0) then
              write(0, *) 'not specified, output all var'
           else
              write(0, *) 'finished, n=', novar
           endif
           exit
        endif
        novar=novar+1
        ovar(novar) = trim(s)
        write(0, '(a)', advance='no') ' more variable (blank to end): '
      enddo

      ! filename
      ! i,j cell to dump
      ! time to dump
      ! variable to dump

   end subroutine

   subroutine usage()
      print '(99a)', 'usage:  ', trim(prgname) , ' [ idbg jdbg ] _m3d_file_'
      print '(99a)', '        ', trim(prgname) , ' -'
      print*
      print*, 'first method dumps values of all variables to screan, more of debugging output'
      print*, 'second method prompt further input from STDIN, and create csv format output'
      stop
   end subroutine

end program

! vim: et sw=3