[cig-commits] r22276 - seismo/3D/GRD_CMT3D/cmt3d

[liuqy at geodynamics.org]

Author: liuqy
Date: 2013-06-14 15:01:39 -0700 (Fri, 14 Jun 2013)
New Revision: 22276

Modified:
   seismo/3D/GRD_CMT3D/cmt3d/INVERSION.PAR
   seismo/3D/GRD_CMT3D/cmt3d/Makefile
   seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub.f90
Log:
Update the way optional parameters are read in. Now the code should be able to handle old style parameter files (without global_coord or read_weight) exactly the same as before (in local coord).



Modified: seismo/3D/GRD_CMT3D/cmt3d/INVERSION.PAR
===================================================================
--- seismo/3D/GRD_CMT3D/cmt3d/INVERSION.PAR	2013-06-14 21:52:26 UTC (rev 22275)
+++ seismo/3D/GRD_CMT3D/cmt3d/INVERSION.PAR	2013-06-14 22:01:39 UTC (rev 22276)
@@ -1,7 +1,7 @@
 CMTSOLUTION
 CMTSOLUTION_NEW
-6   11  0.  0.                  -- npar[, utm_zone, utm_center_x,y]
-0.5 0.5  1.0e22             -- dX/Y(km),ddepth(km),dmoment(dyne.cm)
+6   .false.                  -- npar[, global_coord]
+0.5 0.5  1.0e22             -- dlocation,ddepth,dmoment(dyne.cm)
 flexwin.out
 .true. .false.                   -- weigh_data_files[, read_weight_from_file]
 2 2 1  0.5  1.15 0.55 0.78       -- weights of data comp, az(exp), dist(exp)

Modified: seismo/3D/GRD_CMT3D/cmt3d/Makefile
===================================================================
--- seismo/3D/GRD_CMT3D/cmt3d/Makefile	2013-06-14 21:52:26 UTC (rev 22275)
+++ seismo/3D/GRD_CMT3D/cmt3d/Makefile	2013-06-14 22:01:39 UTC (rev 22276)
@@ -5,7 +5,7 @@
 
 LIB = /opt/sac/lib/libsacio.a
 
-CMT_F90_SRC = get_cmt utm_geo \
+CMT_F90_SRC = get_cmt  \
     cmt3d_constants cmt3d_variables cmt3d_sub4 cmt3d_sub3 cmt3d_sub2 cmt3d_sub
 CMT_F77_SRC = cmt3d_utils distaz
 CMT_SRC = cmt3d_flexwin

Modified: seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub.f90
===================================================================
--- seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub.f90	2013-06-14 21:52:26 UTC (rev 22275)
+++ seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub.f90	2013-06-14 22:01:39 UTC (rev 22276)
@@ -15,26 +15,28 @@
     
     character(len=*), intent(in) :: par_file
     integer ios
+    character(len=250) :: line
 
     open(IOPAR,file=par_file,status='old',iostat=ios)    
     if (ios /= 0) stop 'Error opening inversion parameter file'
 
     read(IOPAR,'(a)') cmt_file
     read(IOPAR,'(a)') new_cmt_file
-    read(IOPAR,'(i2)',iostat=ios,advance='no') npar
+    read(IOPAR,'(a)') line
+    read(line,*, iostat=ios) npar
     if (ios /= 0) stop 'Error reading npar'
-    global_coord=.false.
-
-    read(IOPAR,*,iostat=ios) global_coord
+    read(line(3:),*,iostat=ios) global_coord
+    ! if no global_coord specified
+    if (ios /= 0)  global_coord=.false.
     
-    if (.not. global_coord) then
+    if (.not. global_coord) then ! local coordinates (DEP, LON, LAT)
        par_name = (/'Mrr','Mtt','Mpp','Mrt','Mrp', 'Mtp','dep','lon','lat', &
             'ctm','hdr'/) 
        SCALE_PAR =  &
             (/ SCALE_MOMENT, SCALE_MOMENT, SCALE_MOMENT, SCALE_MOMENT, &
             SCALE_MOMENT, SCALE_MOMENT, SCALE_DEPTH, SCALE_DELTA, SCALE_DELTA, &
             SCALE_CTIME, SCALE_HDUR /)
-    else 
+    else  ! global coordinates (X, Y, Z)
        par_name = (/'Mxx','Myy','Mzz','Mxy','Mxz', 'Myz','xxx','yyy','zzz', &
             'ctm','hdr'/) 
        if (npar == 7) stop 'depth only inversion is not allowed in global coordinates'
@@ -45,18 +47,26 @@
     endif
      
     read(IOPAR,*) ddelta,ddepth,dmoment
-    dcmt_par = (/dble(dmoment),dble(dmoment),dble(dmoment), &
-                 dble(dmoment),dble(dmoment),dble(dmoment),&
-                 dble(ddepth),dble(ddelta),dble(ddelta), &
-                 1.0d0, 1.0d0/) / SCALE_PAR
+    if (.not. global_coord) then
+       dcmt_par = (/dble(dmoment),dble(dmoment),dble(dmoment), &
+            dble(dmoment),dble(dmoment),dble(dmoment),&
+            dble(ddepth),dble(ddelta),dble(ddelta), &
+            1.0d0, 1.0d0/) / SCALE_PAR
+    else ! global code have different dcmt_par
+       dcmt_par = (/dble(dmoment),dble(dmoment),dble(dmoment), &
+            dble(dmoment),dble(dmoment),dble(dmoment),&
+            dble(ddepth),dble(ddepth),dble(ddepth), &
+            1.0d0, 1.0d0/) / SCALE_PAR
+    endif
 
     read(IOPAR,'(a)') flexwin_out_file
 
 ! now allow the possibility of reading in weights from input file
-    read(IOPAR,'(l7)',iostat=ios, advance='no') weigh_data_files
+    read(IOPAR, '(a)') line
+    read(line,*,iostat=ios) weigh_data_files
     if (ios /= 0) stop 'Error reading weigh_data_files'
 
-    read(IOPAR,'(l7)',iostat=ios) read_weight
+    read(line(8:),'(l7)',iostat=ios) read_weight
     if (ios /= 0) read_weight=.false.
     
 ! here we assume that Pnl wave window is the first window out 
@@ -162,8 +172,9 @@
        endif
        do j = 1, nwins(i)
           if (read_weight) then
-            ! LQY: tjunk: time shift for data/syn within the window file
-            ! maybe used in the future to replace local corr subroutine
+            ! LQY: tjunk: time shift for data/syn within the window file is not 
+            ! used at this stage, 
+            ! it may be used in the future to replace local corr subroutine
             read(IOWIN,*,iostat=ios) tstart, tend, tjunk, data_weights(nwin_total+j)
             if (ios /= 0) stop 'ts, te, tshift, weight are expected!'
           else