 |
forums.silverfrost.com
Welcome to the Silverfrost forums
|
| View previous topic :: View next topic |
| Author |
Message |
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
 Posted: Mon Jul 15, 2013 8:14 pm Post subject: Fortran 95 code problem |
 |
|
Hi all
The program I have fails with errors on the reals, grabbed a sample of some of the lines (not all)
MODULE MOD_PARAM
REAL(8), PARAMETER :: R_t=6370,pi=3.14159,G=6.67D-11
END MODULE
MODULE MOD_EXPAN
REAL(8) d_x,d_y,d_z,xc,yc,zc,Gp_as,d_rho
END MODULE
MODULE MOD_NEAR
REAL(8) rad,d_h(2,2),G_cp,rho_w,rho_c,h
END MODULE
USE MOD_EXPAN
USE MOD_NEAR
USE MOD_PARAM
REAL(8) del_xd,del_xi,cx,cy,del_xBg,R_d,R_i,del_xdet,del_x_pr
REAL(8) AB
REAL(8) x,y,z,del_g,Gpr,x1,y1,x2,y2,z1,z2,zm,E_pr
REAL(8) L_x,L_y,BS
REAL(8) BB,nu,al
REAL(8) l_i1,l_i2,l_j1,l_j2,dm
Compiling and linking file: FA2Boug_CGv1B.f95
C:\Temp\FA2Boug_CGv1B.F95(34) : error 62 - Invalid KIND specifier
C:\Temp\FA2Boug_CGv1B.F95(34) : comment 981 - Specifying the kind of the type REAL with a constant is non-portable - 'SELECTED_REAL_KIND(6,37)' would be better
C:\Temp\FA2Boug_CGv1B.F95(3Cool : error 62 - Invalid KIND specifier
Not sure if there is a way to post the whole code, seems to truncate it here.
Any ideas how to solve this issue?
Cheers
Lester |
|
| Back to top |
|
 |
simon
Joined: 05 Jul 2006
Posts: 299
|
| Posted: Mon Jul 15, 2013 8:51 pm Post subject: |
|
|
probably is not doing what you want it to; specifically it is not equivalent to . The 8 is not the number of bytes. If you want double precision, you could add the following line at the top of the program:
| Code: | INTEGER, PARAMETER :: dp=KIND(1.0_d0)
|
Then change all to . |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Mon Jul 15, 2013 9:12 pm Post subject: Re: |
|
|
| simon wrote: | probably is not doing what you want it to; specifically it is not equivalent to . The 8 is not the number of bytes. If you want double precision, you could add the following line at the top of the program:
| Code: | INTEGER, PARAMETER :: dp=KIND(1.0_d0)
|
Then change all to . |
Hello Simon.
Thanks for the idea. I tried that and it fixed most of the problem, but now get:
Compiling and linking file: FA2Boug_CGv1B.f95
C:\Temp\FA2Boug_CGv1B.F95(32) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(34) : error 1040 - MODULE cannot appear in a PROGRAM block
Compilation failed.
#32 INTEGER, PARAMETER :: dp=KIND(1.0_d0)
#33
#34 MODULE MOD_PARAM
#35 REAL(dp), PARAMETER :: R_t=6370,pi=3.14159,G=6.67D-11
#36 END MODULE
So the problem looks to beat lines 32-34. Is it an issue with the module paprameter statement which has integer and reals?
Cheers
Lester |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Mon Jul 15, 2013 10:59 pm Post subject: Re: |
|
|
| arctica wrote: | | simon wrote: | probably is not doing what you want it to; specifically it is not equivalent to . The 8 is not the number of bytes. If you want double precision, you could add the following line at the top of the program:
| Code: | INTEGER, PARAMETER :: dp=KIND(1.0_d0)
|
Then change all to . |
Hello Simon.
Thanks for the idea. I tried that and it fixed most of the problem, but now get:
Compiling and linking file: FA2Boug_CGv1B.f95
C:\Temp\FA2Boug_CGv1B.F95(32) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(34) : error 1040 - MODULE cannot appear in a PROGRAM block
Compilation failed.
#32 INTEGER, PARAMETER :: dp=KIND(1.0_d0)
#33
#34 MODULE MOD_PARAM
#35 REAL(dp), PARAMETER :: R_t=6370,pi=3.14159,G=6.67D-11
#36 END MODULE
So the problem looks to beat lines 32-34. Is it an issue with the module paprameter statement which has integer and reals?
Cheers
Lester |
Tried something else, based on the earlier error message:
!
! INTEGER, PARAMETER :: dp=KIND(1.0_d0)
INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(6,37)
MODULE MOD_PARAM
REAL(dp), PARAMETER :: R_t=6370,pi=3.14159,G=6.67D-11
END MODULE
MODULE MOD_EXPAN
REAL(dp) d_x,d_y,d_z,xc,yc,zc,Gp_as,d_rho
END MODULE
MODULE MOD_NEAR
REAL(dp) rad,d_h(2,2),G_cp,rho_w,rho_c,h
END MODULE
USE MOD_EXPAN
USE MOD_NEAR
USE MOD_PARAM
This time I got:
Compiling and linking file: FA2Boug_CGv1B.f95
C:\Temp\FA2Boug_CGv1B.F95(35) : error 1040 - MODULE cannot appear in a PROGRAM block
Compilation failed.
Not sure what changed and why - but still not compiling/building |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2615
Location: Sydney
|
| Posted: Tue Jul 16, 2013 12:57 am Post subject: |
|
|
@simon,
In real*8, the 8 is the number of bytes !!! This is my preferred and recomended way of declaring an 8 byte real.
@arctica,
You need to provide the definition of dp (correctly) to each of the modules. The following code would do that.
I have included some indication of the precision of the real types.
| Code: | MODULE MOD_Precision
INTEGER, PARAMETER :: sp=SELECTED_REAL_KIND(6,30)
INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(12,100)
INTEGER, PARAMETER :: qp=SELECTED_REAL_KIND(18,500)
END MODULE MOD_Precision
MODULE MOD_PARAM
USE MOD_Precision
REAL(dp), PARAMETER :: R_t=6370, G=6.67D-11 ! , pi=3.14159
REAL(dp) :: pi
REAL(qp) :: piq
REAL(sp) :: pis
END MODULE
MODULE MOD_EXPAN
USE MOD_Precision
REAL(dp) d_x,d_y,d_z,xc,yc,zc,Gp_as,d_rho
END MODULE
MODULE MOD_NEAR
USE MOD_Precision
REAL(dp) rad,d_h(2,2),G_cp,rho_w,rho_c,h
END MODULE
USE MOD_EXPAN
USE MOD_NEAR
USE MOD_PARAM
REAL(dp) del_xd,del_xi,cx,cy,del_xBg,R_d,R_i,del_xdet,del_x_pr
REAL(dp) AB
REAL(dp) x,y,z,del_g,Gpr,x1,y1,x2,y2,z1,z2,zm,E_pr
REAL(dp) L_x,L_y,BS
REAL(dp) BB,nu,al
REAL(dp) l_i1,l_i2,l_j1,l_j2,dm
!
pi = 4.0_dp * ATAN (1.0_dp)
piq = 4.0_qp * ATAN (1.0_qp)
pis = 4.0_sp * ATAN (1.0_sp)
write (*,*) pis, kind(pis), precision(pis)
write (*,*) pi, kind(pi), precision(pi)
write (*,*) piq, kind(piq), precision(piq)
write (*,'(f27.19)') pis
write (*,'(f27.19)') pi
write (*,'(f27.19)') piq
!
END
|
John |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2615
Location: Sydney
|
| Posted: Tue Jul 16, 2013 9:10 am Post subject: |
|
|
You did not give any examples of the types of errors being reported !
The definition and use of modules is somewhat restricted in it's placement in the code.
1) A module definition, is like a seperate routine. It can not occur inside the program or subroutine/function definition.
2) A module USE must be the first statement in a program or function. it must occur before any other declaration or computation statements.
The example I attached complies with this, as does the latest (limited) attachment you provided, where the module definitions are the first group of statements.
Make sure your USE statements are before any other declarations. This could be your problem.
The other possible issue is that all variable names must be unique. you should not declare the same variable name in different modules, or as a local variable in the routine. While this is permitted in some cases, it would be a nightmare to debug.
The introduction of MODULE to fortran has been a good change, as it provides a structure to define your variable database.
Include some error examples if I have not covered what is happening.
John |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 9:10 am Post subject: Re: |
|
|
| JohnCampbell wrote: | @simon,
In real*8, the 8 is the number of bytes !!! This is my preferred and recomended way of declaring an 8 byte real.
@arctica,
You need to provide the definition of dp (correctly) to each of the modules. The following code would do that.
I have included some indication of the precision of the real types.
| Code: | MODULE MOD_Precision
INTEGER, PARAMETER :: sp=SELECTED_REAL_KIND(6,30)
INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(12,100)
INTEGER, PARAMETER :: qp=SELECTED_REAL_KIND(18,500)
END MODULE MOD_Precision
MODULE MOD_PARAM
USE MOD_Precision
REAL(dp), PARAMETER :: R_t=6370, G=6.67D-11 ! , pi=3.14159
REAL(dp) :: pi
REAL(qp) :: piq
REAL(sp) :: pis
END MODULE
MODULE MOD_EXPAN
USE MOD_Precision
REAL(dp) d_x,d_y,d_z,xc,yc,zc,Gp_as,d_rho
END MODULE
MODULE MOD_NEAR
USE MOD_Precision
REAL(dp) rad,d_h(2,2),G_cp,rho_w,rho_c,h
END MODULE
USE MOD_EXPAN
USE MOD_NEAR
USE MOD_PARAM
REAL(dp) del_xd,del_xi,cx,cy,del_xBg,R_d,R_i,del_xdet,del_x_pr
REAL(dp) AB
REAL(dp) x,y,z,del_g,Gpr,x1,y1,x2,y2,z1,z2,zm,E_pr
REAL(dp) L_x,L_y,BS
REAL(dp) BB,nu,al
REAL(dp) l_i1,l_i2,l_j1,l_j2,dm
!
pi = 4.0_dp * ATAN (1.0_dp)
piq = 4.0_qp * ATAN (1.0_qp)
pis = 4.0_sp * ATAN (1.0_sp)
write (*,*) pis, kind(pis), precision(pis)
write (*,*) pi, kind(pi), precision(pi)
write (*,*) piq, kind(piq), precision(piq)
write (*,'(f27.19)') pis
write (*,'(f27.19)') pi
write (*,'(f27.19)') piq
!
END
|
John |
Still getting errors John. Is there anyway to upload the whole code? The "code" tags are not working here to do that for some reason.
Cheers
Lester |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2615
Location: Sydney
|
| Posted: Tue Jul 16, 2013 9:13 am Post subject: |
|
|
| how about listing a few errors. This would give me some idea of the type of problem you are having. |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 9:23 am Post subject: Re: |
|
|
| JohnCampbell wrote: | You did not give any examples of the types of errors being reported !
The definition and use of modules is somewhat restricted in it's placement in the code.
1) A module definition, is like a seperate routine. It can not occur inside the program or subroutine/function definition.
2) A module USE must be the first statement in a program or function. it must occur before any other declaration or computation statements.
The example I attached complies with this, as does the latest (limited) attachment you provided, where the module definitions are the first group of statements.
Make sure your USE statements are before any other declarations. This could be your problem.
The other possible issue is that all variable names must be unique. you should not declare the same variable name in different modules, or as a local variable in the routine. While this is permitted in some cases, it would be a nightmare to debug.
The introduction of MODULE to fortran has been a good change, as it provides a structure to define your variable database.
Include some error examples if I have not covered what is happening.
John |
I will have to retrieve the errors when I am home (at work just now), but I am wondering if the compiler setup may be different? I am working with the Personal Edition, just opened the code in the main Plato window and set it to build.
The errors I got this morning applying your code change, were all about types of variable, but will post the list of errors later on. Thanks for now
Cannot see a way to load the whole program which would be a big help, it's not that big only about 20k max.
Lester |
|
| Back to top |
|
|
davidb
Joined: 17 Jul 2009
Posts: 560
Location: UK
|
| Posted: Tue Jul 16, 2013 4:33 pm Post subject: |
|
|
The code for setting dp should be
| Code: |
integer, parameter :: dp = kind(1.0d0)
|
NOT kind(1.0_d0)
Your code above that uses selected_real_kind is nearly correct.
Just add :: to each real statement with a KIND specifier and it will compile.
e.g. use
real(dp) :: a
not
real(dp) a
_________________
Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 10:11 pm Post subject: Re: |
|
|
Errors from adding code from John:
[size=9]Compiling and linking file: FA2Boug_CGv1B.f95
C:\Temp\FA2Boug_CGv1B.F95(42) : error 62 - Invalid KIND specifier
C:\Temp\FA2Boug_CGv1B.F95(61) : error 636 - KIND parameter out of range, permitted KINDs are 1, 2, or 3
C:\Temp\FA2Boug_CGv1B.F95(61) : error 636 - KIND parameter out of range, permitted KINDs are 1, 2, or 3
C:\Temp\FA2Boug_CGv1B.F95(72) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(72) : error 775 - More than one main program encountered
C:\Temp\FA2Boug_CGv1B.F95(73) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(74) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(75) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(76) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(77) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(81) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(84) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(87) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(435) : error 381 - D_H is not an array
C:\Temp\FA2Boug_CGv1B.F95(436) : error 215 - Invalid expression on left hand side of assignment
C:\Temp\FA2Boug_CGv1B.F95(437) : error 215 - Invalid expression on left hand side of assignment
C:\Temp\FA2Boug_CGv1B.F95(438) : error 215 - Invalid expression on left hand side of assignment
C:\Temp\FA2Boug_CGv1B.F95(474) : error 176 - Only INTEGER constants are allowed as a KIND parameter
C:\Temp\FA2Boug_CGv1B.F95(491) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\FA2Boug_CGv1B.F95(533) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\FA2Boug_CGv1B.F95(535) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\FA2Boug_CGv1B.F95(537) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\FA2Boug_CGv1B.F95(539) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\FA2Boug_CGv1B.F95(541) : warning 179 - Comparing floating point quantities for equality may give misleading results
Compilation failed.
[/size] |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 10:23 pm Post subject: Fortran_code |
|
|
I am trying to add the whole program in small chunks
| Code: | !
! PROGRAM FA2Boug.f:
! Fortran 90 code to calculate Bouguer anomaly applying
! Bullard A, B and C corrections to a Free Air (FA)/slab-corrected
! Bouguer anomaly grid, provided an extra grid with elevation
! (topo/bathymetry)
! The program can calculate Bouguer anomaly in in either land or
! sea points
! Contributions to Bouguer anomaly are computed in different zones
! according the radial distance, R, to the calculation point:
! DISTANT ZONE: (R_d>R>R_i)
! INTERMEDIATE ZONE: (R_i>R>del_xi/2)
! INNER ZONE: (R<del_xi/2)
! Input files required by FA2Boug.f:
! *parameters.dat
! *topo_cart.xyz
! *gravi_cart.xyz
! Output files are:
! *bouguer.xyz
! *bouguer_slab.xyz
! Optionally a detailed topography file, *topo_cart_det.xyz,
! with a gridstep del_xdet(<del_xi) can be provided. In that
! case, a new calculation is performed in the
! DETAILED INTERMEDIATE ZONE (del_xi/2>R>del_xdet/2), and
! the INNER ZONE is redefined as (R<del_xdet/2)
! Also optionally, the program calculate the isostatic residual
! anomaly using an Airy compensation model for the crust
! Written by Javier Fullea Urchulutegui
! Instituto de Ciencias de la Tierra Jaume Almera (CSIC)
! Sept 2007
!
MODULE MOD_Precision
INTEGER, PARAMETER :: sp=SELECTED_REAL_KIND(6,30)
INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(12,100)
INTEGER, PARAMETER :: qp=SELECTED_REAL_KIND(18,500)
END MODULE MOD_Precision
MODULE MOD_PARAM
USE MOD_Precision
REAL(dp), PARAMETER :: R_t=6370,G=6.67D-11 ! ,pi=3.14159
REAL(dp) :: pi
REAL(qp) :: piq
REAL(sp) :: pis
END MODULE
MODULE MOD_EXPAN
USE MOD_Precision
REAL(dp) d_x,d_y,d_z,xc,yc,zc,Gp_as,d_rho
END MODULE
MODULE MOD_NEAR
USE MOD_Precision
REAL(dp) rad,d_h(2,2),G_cp,rho_w,rho_c,h
END MODULE
USE MOD_EXPAN
USE MOD_NEAR
USE MOD_PARAM
!
pi = 4.0_dp * ATAN (1.0_dp)
piq = 4.0_qp * ATAN (1.0_qp)
pis = 4.0_sp * ATAN (1.0_sp)
write (*,*) pis, kind(pis), precision(pis)
write (*,*) pi, kind(pi), precision(pi)
write (*,*) piq, kind(piq), precision(piq)
write (*,'(f27.19)') pis
write (*,'(f27.19)') pi
write (*,'(f27.19)') piq
!
END
REAL(dp) del_xd,del_xi,cx,cy,del_xBg,R_d,R_i,del_xdet,del_x_pr
REAL(dp) AB
REAL(dp) x,y,z,del_g,Gpr,x1,y1,x2,y2,z1,z2,zm,E_pr
REAL(dp) L_x,L_y,BS
REAL(dp) BB,nu,al
REAL(dp) l_i1,l_i2,l_j1,l_j2,dm
INTEGER N,M,k,l,IWIN_i,IWIN,rtio,rt,rt_i
INTEGER l_lst,l_fst,k_lst,k_fst,up_INT
INTEGER cont_x,cont_y,N_tot,N_nod,n_strg,slab_Boug
REAL(dp), ALLOCATABLE:: E(:,:),FA(:,:),E_det(:,:)
CHARACTER(10), ALLOCATABLE:: FA_strg(:,:),E_det_strg(:,:)
CHARACTER(10) COORD,land
REAL(dp) a(3,3),R_i_L,R_i_R,R_i_D,R_i_U
INTEGER det,N_det,M_det,IWIN_det
LOGICAL det_on
REAL(dp) d_rho_mho,t,ISOS,z_cref
!//////////////////INPUT///////////////////////////// |
Last edited by arctica on Tue Jul 16, 2013 10:37 pm; edited 1 time in total |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 10:26 pm Post subject: |
|
|
| Code: | ! Program inputs in parameters.dat:
! Perform a detailed topography calculation on the land areas
! present in the file topo_cart_det.xyz:
! det=1-->yes
! det=0-->no
! del_xd-->grid step Distant zone (km)
! del_xi-->grid step Intermediate zone (km)
! del_xBg-->grid step Bouguer output grid (km)
! rho_c-->crust reduction density (kg/m�)
! rho_w-->water reduction density (kg/m�)
! N-->number of nodes of elevation an FA input grids in x
! M-->number of nodes of elevation an FA input grids in y
! R_i-->Limit of the intermediate zone (km)
! R_d-->Limit of the distant zone (km)
! land='on'-->calculates Bouguer anomaly at land and sea points
! land='off'->calculates Bouguer anomaly only at sea points
! slab_Boug=1-->INPUT GRAVITY ANOMALY IS SLAB-CORRECTED BOUGUER
! slab_Boug=0-->INPUT GRAVITY ANOMALY IS FREE AIR
! If det=1, then, the input file must contain the following additional
! information:
! N_det-->number of nodes of topo_cart_det.xyz input grid in x
! M_det-->number of nodes of topo_cart_det.xyz input grid in y
! del_xdet-->grid step of the Detailed Intermediate zone (km)
OPEN(1, file='parameters.dat', status='OLD')
READ(1,'(I1)',ADVANCE='NO') det
IF (det==1) THEN
READ(1,*)del_xd,del_xi,del_xBg,rho_c,rho_w,N,M,R_d,R_i,land, &
slab_Boug,N_det,M_det,del_xdet
write(*,*)del_xd,del_xi,del_xBg,rho_c,rho_w,N,M,R_d,R_i,land, &
slab_Boug,N_det,M_det,del_xdet
ELSE
READ(1,*)del_xd,del_xi,del_xBg,rho_c,rho_w,N,M,R_d,R_i,land, &
slab_Boug
ENDIF
CLOSE(1)
!Open INPUT files
open(2,file='topo_cart.xyz', status='OLD')
open(5,file='gravi_cart.xyz',status='OLD')
!Open detailed topography, if necessary
IF (det==1) open(7,file='topo_cart_det.xyz', status='OLD')
!Open OUTPUT files
open(15,file='bouguer.xyz',status='UNKNOWN')
open(17,file='bouguer_slab.xyz',status='UNKNOWN')
!//////////////////INPUT/////////////////////////////
write(*,*)''
write(*,*)''
write(*,*)' ***FA2BOUG***'
IF (land=='on') THEN
write(*,*)'CALCULATES IN LAND AND SEA POINTS'
ELSE
write(*,*)'CALCULATES ONLY IN SEA POINTS'
ENDIF
write(*,*)''
IF (slab_Boug==1) THEN
write(*,*)'INPUT GRAVITY ANOMALY IS SLAB-CORRECTED BOUGUER'
ELSE
write(*,*)'INPUT GRAVITY ANOMALY IS FREE AIR'
ENDIF
IF (det==1) THEN
write(*,*)'DETAILED TOPOGRAPHY ON'
ELSE
write(*,*)'DETAILED TOPOGRAPHY OFF'
ENDIF
IWIN=nint(R_d/del_xi)
rtio=nint(del_xBg/del_xi)
rt=nint(del_xd/del_xi)
N_tot=N*M/rtio**2
!Dimensionate the elevation and FA matrices |
|
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 10:28 pm Post subject: |
|
|
| Code: | ALLOCATE(E(N,M),FA(N,M),FA_strg(N,M))
IF (det==1) ALLOCATE(E_det(N_det,M_det),E_det_strg(N_det,M_det))
! Read elevation and FA files
!
DO j=M,1,-1
DO i=1,N
read(2,*) E(i,j)
read(5,*) FA_strg(i,j)
ENDDO
ENDDO
REWIND (2)
!Read the detailed topography file (if present)
IF (det==1) THEN
call cpu_time(time1)
!Flag for NaN values
E_det=-1D6
DO j=M_det,1,-1
DO i=1,N_det
read(7,*) E_det_strg(i,j)
n_strg=ICHAR(E_det_strg(i,j))
!Skips NaN and store only the numerical values (i.e. points
!where the detailed topography is present in the file topo_cart_det.xyz)
IF ((n_strg<=47.OR.n_strg>57).AND.n_strg/=45) THEN
CYCLE
ELSE
READ(E_det_strg(i,j),*)E_det(i,j)
ENDIF
ENDDO
ENDDO
DEALLOCATE(E_det_strg)
call cpu_time(time2)
write(*,*)'DETAILED TOPO READ, ELAPSED TIME: ', time2-time1
IWIN_det=nint((del_xi-del_xdet)/(2*del_xdet))
ENDIF
R_tl1=R_i*1D3*sqrt(2.)
R_tl2=(R_d+1)*1D3*sqrt(2.)
!
! Loop for the FA anomaly grid
WRITE(*,*)'%%%%%%%%%%%%%%%%%%'
WRITE(*,*)'START CALCULATION'
WRITE(*,*)'%%%%%%%%%%%%%%%%%%'
cont_y=0
DO j=M,1,-rtio
cont_y=cont_y+1
cont_x=0
cy=(j-1)*del_xi
DO i=1,N,rtio
cont_x=cont_x+1
cx=(i-1)*del_xi
!Percentage of program execution
N_nod=(cont_y-1)*N/rtio+cont_x
IF (N_nod==NINT(N_tot*.1)) THEN
WRITE(*,*)'**10%**'
ELSEIF (N_nod==NINT(N_tot*.3)) THEN
WRITE(*,*)'****30%****'
ELSEIF (N_nod==NINT(N_tot*.5)) THEN
WRITE(*,*)'********50%********'
ELSEIF (N_nod==NINT(N_tot*.9)) THEN
WRITE(*,*)'***********80%************'
ENDIF
IF (i<IWIN+1.OR.i>N-IWIN.OR. &
j<IWIN+1.OR.j>M-IWIN) CYCLE
!Elevation of the calculation point
IF (det==1) THEN
det_on=.TRUE.
rt_i=nint(del_xi/del_xdet)
IWIN_i=nint(R_i/del_xi)*rt_i
rad=sqrt(2.)*del_xdet*1D3/2
i_det=NINT((cx-R_d+R_i)/del_xdet)+1
j_det=NINT((cy-R_d+R_i)/del_xdet)+1
IF (E_det(i_det,j_det)==-1D6.OR. &
E_det(i_det+IWIN_i,j_det)==-1D6.OR. &
E_det(i_det-IWIN_i,j_det)==-1D6.OR. &
E_det(i_det,j_det+IWIN_i)==-1D6.OR. &
E_det(i_det,j_det-IWIN_i)==-1D6.OR. &
E_det(i_det+IWIN_i,j_det+IWIN_i)==-1D6.OR. &
E_det(i_det+IWIN_i,j_det-IWIN_i)==-1D6.OR. &
E_det(i_det-IWIN_i,j_det+IWIN_i)==-1D6.OR. &
E_det(i_det-IWIN_i,j_det-IWIN_i)==-1D6.OR. &
h<0 ) det_on=.FALSE.
ENDIF |
|
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 10:30 pm Post subject: |
|
|
| Code: | !Elevation of the calculation point and parameters
!for the Intermediate zone
IF (det==1.AND.det_on) THEN
h=E_det(i_det,j_det)
ELSE
h=E(i,j)
IWIN_i=nint(R_i/del_xi)
rt_i=1
rad=sqrt(2.)*del_xi*1D3/2
ENDIF
!Skips land calculation if requiered
IF (h>0.AND.land=='off') CYCLE
!Skips NaN (FA no-data points)
n_strg=ICHAR(FA_strg(i,j))
IF ((n_strg<=47.OR.n_strg>57).AND.n_strg/=45) CYCLE
!Convert string to real
READ(FA_strg(i,j),*)FA(i,j)
IF (h<=0) THEN
rho=rho_c-rho_w
nu=-h/(R_t*1D3+h)
zm=0D0
ELSE
rho=rho_c
nu=h/(R_t*1D3+h)
zm=h
ENDIF
!Write simple Bouguer anomaly (i.e. only Bullard A correction)
BS=(2*pi*G*rho*h)*1D5
!Starts the value of Bouguer Anomaly
IF (slab_Boug==0) THEN
AB=FA(i,j)
WRITE(17,*)cx,cy,FA(i,j)-BS
ELSE
AB=FA(i,j)+BS
WRITE(17,*)cx,cy,AB
ENDIF
!Compute the curvature correction in mgal (Bullard B)
!using Whitman approximation (Whitman, Geophysics, 56 N12,
!pp 1980-1985,(1991))
al=R_d/R_t
BB=-2*pi*G*rho*h*(al/2-nu)*1D5
AB=AB+BB
!Loop for the external square (length=2*R_d) Distant zone
!
!Find the limits between the distant and intermediate zones
l_i1=INT((i-NINT((R_i+del_xd*.5)/del_xi)))
l_j1=INT((j-NINT((R_i+del_xd*.5)/del_xi)))
l_i2=INT((i+NINT((R_i+del_xd*.5)/del_xi)))
l_j2=INT((j+NINT((R_i+del_xd*.5)/del_xi)))
!The last and first nodes of the distant zone at each side of
!the calculation point (i,j)
l_lst=i-IWIN+INT((l_i1-(i-IWIN))/rt)*rt
l_fst=i+IWIN-INT(((i+IWIN)-l_i2)/rt)*rt
k_lst=j-IWIN+INT((l_j1-(j-IWIN))/rt)*rt
k_fst=j+IWIN-INT(((j+IWIN)-l_j2)/rt)*rt
!Limits R_i
R_i_L=((i-l_lst)*del_xi-del_xd*.5)*1D3
R_i_R=((l_fst-i)*del_xi-del_xd*.5)*1D3
R_i_D=((j-k_lst)*del_xi-del_xd*.5)*1D3
R_i_U=((k_fst-j)*del_xi-del_xd*.5)*1D3
DO k=j-IWIN,j+IWIN,rt
DO l=i-IWIN,i+IWIN,rt
!Select the DISTANT ZONE (R_d>R>R_i)
IF (l<=l_i1.OR.l>=l_i2.OR.k<=l_j1.OR.k>=l_j2) THEN
!Onshore/offshore prism
IF (E(l,k)>0) THEN
d_rho=rho_c
ELSE
d_rho=rho_c-rho_w
ENDIF
!Vertical limits of the prism
z1=zm
z2=ABS(E(l,k)-zm)
!Sides of the prism
d_x=del_xd*1D3
d_y=del_xd*1D3
d_z=z2-z1
!Cartesian coordinates of CM of the prism
xc=abs((l-i)*del_xi)*1D3
yc=abs((k-j)*del_xi)*1D3
zc=(z1+z2)*.5
CALL EXPAN
AB=AB+Gp_as
ENDIF!END of the IF for the DISTANT ZONE (R_d>R>R_i)
ENDDO
ENDDO
!
! |
|
|
| Back to top |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ
Search
Memberlist
Usergroups
Register
Profile
Log in to check your private messages
Log in
