The index of the loop that outputs the model file

[MistBornPow]

Dear Dr. Fang,

Thank you for your reply. I have additional questions regarding issue #22. There are two segments of the code responsible for writing the model output file that I find confusing. The first is the segment in main.f90:

open(64,file=outmodel) do k=1,nz-1 do j=1,ny-2 do i=1,nx-2 write(64,'(5f10.5)') gozd+(j-1)dvzd,goxd-(i-1)dvxd,depz(k),vsf(i+1,j+1,k) enddo enddo enddo close(64)

In this loop, the indices for x and y axis start from 1 end at nx-2, ny-2, respectively. The written latitude and longitude start from "gozd" and "goxd", while the model grid starts from "vsf(2,2,k)". Since the indices for latitude, longitude and the model "vsf" seems inconsistent, I'm not quite sure about this. Does this mean the model grid (2,2,:) corresponds to "gozd" and "goxd", while the model grid (1,1,:) corresponds to "gozd-dvzd" and "goxd+dvxd"? I would appreciate your clarification on this matter.

Second is the segment in CalSurfG.f90, subroutine "synthetic":

if(kmaxRc.gt.0) then iwave=2 igr=0 call caldespersion(nx,ny,nz,vels,pvRc, & iwave,igr,kmaxRc,tRc,depz,minthk)

open(62,file='velmap2dRc.dat') do k = 1,kmaxRc do j=1,ny-2 do i=1,nx-2 write(62,'(5f8.4)') gozd+(j-1)dvzd,goxd-(i-1)dvxd,tRc(k),pvRc((j+1)*nx+i+1,k) enddo enddo enddo close(62) endif

Here, I would like to assume that my previous understanding is correct. Since the first model grid to be written is at (2,2,k) but the dimension of "pvRc" is (nxny, kmaxRc), I believe the first written value of "pvRc" should be "pvRc(nx+2, k)". However, my interpretation differs from the code which starts from "pvRc(2 nx+2, k)". I'm not sure where my understanding goes wrong when trying to comprehend the code. Or, are there dummy values in "pvRc"? Could you give me some hints? It will be very helpful. I've attached a handwritten note to further illustrate my understanding. Could you please confirm if my interpretation is correct?

Best regards, Wei-Li Chen

[HongjianFang]

Dear Wei-Li, Sorry again for the confusion (I did not realize someone would care to read the code, so I just tested it and left it that way). Anyway, you are right that 'gozd' and 'goxd' are the actual starting points for the inverted model, there is another boundary layer with one dvxd/dvzd outside. For the second segment, I thought I had deleted those as they are only for testing: comparing dispersion data. You can ignore those as they are of little use. May I ask why you are interested in reading the code so carefully? I spent lots of hours reading the original FMM codes, in Fortran. Nowadays, there are many readable codes in Python, such as PyKonal written by Malcolm White, which can be easily translated for your own use, with more ease, I think. Let me know if you have any other questions.

Cheers, Hongjian Fang

[MistBornPow]

Dear Dr. Fang,

Thank you for your reply.

There are several reasons why I am studying the code in detail. Firstly, I am currently working on outputting the sensitivity kernel from the code. In doing so, I have realized that this is closely related to the model grid. Previously, I was able to generate the dispersion map output from the final Vs model and successfully convert the data structure to make it compatible with the program fm2dss in FMST package. However, I have encountered some obstacles when trying to ensure that the coordinates are correct. Although I simply copied the original code from DSurfTomo, slight modifications are still needed to meet my requirements. Thus I have to understand the code.

Secondly, I'm currently studying for master degree in National Central University. I use DSurfTomo and other programs or software as a tool in my research. I expect to understand the tools I use so that I can be sure of what I'm really doing. Furthermore, I prefer not to just "copy and paste" from others' theses. I hope to truly learn during my master's program and acquire the ability to write the thesis in my own words. That’s why I feel the need to study the code in detail.

For now, I haven't found a package that can complete the entire inversion routine for surface wave tomography like DSurfTomo. If I'm going to use Direct Inversion scheme for my research but in Python, I will need to combine several packages and write scripts to make it work. However, to do this effectively, I still need to understand the original code of DSurfTomo to integrate these different packages, but I currently lack the ability to do so.

I apologize for taking your time to read through these lengthy messages, and I appreciate your valuable feedback.

Best regards, Wei-Li Chen

[MistBornPow]

I forgot to mention an important detail. The problem I encountered while working on outputting the sensitivity kernel is that I noticed the array structure of "sen_vsRc" is similar with "pvRc". Both compress the dimensions of the x and y axes in the first dimension, and I don't need the sensitivity kernel of the model boundary. Therefore, I planned to follow the code segment that generates the output of "pvRc," but I ran into the problems I mentioned earlier.

Wei-Li Chen