 |
forums.silverfrost.com
Welcome to the Silverfrost forums
|
| View previous topic :: View next topic |
| Author |
Message |
arctica
Joined: 10 Sep 2006
Posts: 148
Location: United Kingdom
|
 Posted: Tue Jul 16, 2013 10:31 pm Post subject: |
 |
|
| Code: |
!Loop for the internal square (length=2*R_i)
!INTERMEDIATE ZONE (R_i>R>del_xi/2)
IF (det==1.AND.det_on) THEN
I_in=i_det
J_in=j_det
ELSE
I_in=i
J_in=j
ENDIF
DO k=J_in-IWIN_i,J_in+IWIN_i,rt_i
DO l=I_in-IWIN_i,I_in+IWIN_i,rt_i
!Onshore/offshore prism
IF (det==1.AND.det_on) THEN
E_pr=E_det(l,k)
del_x_pr=del_xdet
ELSE
E_pr=E(l,k)
del_x_pr=del_xi
ENDIF
IF (E_pr==-1D6) E_pr=zm
IF (E_pr>0) THEN
d_rho=rho_c
ELSE
d_rho=rho_c-rho_w
ENDIF
!Coordinates defining the Flat-Toped-Prism (in m)
x1=((l-I_in)*del_x_pr-.5*del_xi)*1D3
y1=((k-J_in)*del_x_pr-.5*del_xi)*1D3
x2=((l-I_in)*del_x_pr+.5*del_xi)*1D3
y2=((k-J_in)*del_x_pr+.5*del_xi)*1D3
! write(*,*)'x1,x2',x1,x2
!Extend the prisms in the limit of the intermediate zone
!in order to fill the space between the two grids
IF (l==I_in-IWIN_i) THEN
IF(ABS(x1)>R_i_L.AND.ABS(x2)>R_i_L) CYCLE
x1=-R_i_L
ENDIF
IF (l==I_in+IWIN_i) THEN
IF(ABS(x2)>R_i_R.AND.ABS(x1)>R_i_R) CYCLE
x2=R_i_R
ENDIF
IF (k==J_in-IWIN_i) THEN
IF(ABS(y1)>R_i_D.AND.ABS(y1)>R_i_D) CYCLE
y1=-R_i_D
ENDIF
IF (k==J_in+IWIN_i) THEN
IF(ABS(y2)>R_i_U.AND.ABS(y1)>R_i_U) CYCLE
y2=R_i_U
ENDIF
z1=zm
z2=ABS(E_pr-zm)
IF (det==1.AND.det_on.AND.k==j_det.AND.l==i_det) CYCLE
CALL atrac_prisma(x1,x2,y1,y2,z1,z2,d_rho,Gpr)
AB=AB+Gpr
!
ENDDO
ENDDO
!
!Calculate using the detailed topography
IF (det==1.AND.det_on) THEN
DO jj=j_det-IWIN_det,j_det+IWIN_det
DO ii=i_det-IWIN_det,i_det+IWIN_det
d_rho=rho_c
!Coordinates defining the Flat-Toped-Prism (in m)
x1=((ii-i_det)-.5)*del_xdet*1D3
y1=((jj-j_det)-.5)*del_xdet*1D3
x2=((ii-i_det)+.5)*del_xdet*1D3
y2=((jj-j_det)+.5)*del_xdet*1D3
z1=zm
z2=ABS(E_det(ii,jj)-zm)
!Extend the prisms in the limit in order to fill the space
!between the two grids
IF (ii==i_det-IWIN_det) x1=x1-(del_xi*.5-del_xdet*(IWIN_det+.5)) &
*1D3
IF (ii==i_det+IWIN_det) x2=x2+(del_xi*.5-del_xdet*(IWIN_det+.5)) &
*1D3
IF (jj==j_det-IWIN_det) y1=y1-(del_xi*.5-del_xdet*(IWIN_det+.5)) &
*1D3
IF (jj==j_det+IWIN_det) y2=y2+(del_xi*.5-del_xdet*(IWIN_det+.5)) &
*1D3
CALL atrac_prisma(x1,x2,y1,y2,z1,z2,d_rho,Gpr)
AB=AB+Gpr
ENDDO
ENDDO
ENDIF
! |
|
|
| Back to top |
|
 |
arctica
Joined: 10 Sep 2006
Posts: 148
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 10:33 pm Post subject: |
|
|
| Code: |
!Compute Bullard C correction in the NEAR ZONE (R<del_xi/2 or del_xdet/2)
DO ih=1,3
DO ik=1,3
IF (det==1.AND.det_on) THEN
a(ih,ik)=E_det(i_det+ih-2,j_det+ik-2)
ELSE
a(ih,ik)=E(i+ih-2,j+ik-2)
ENDIF
ENDDO
ENDDO
d_h(1,1)=SUM(a(1:2,1:2))/4-h
d_h(1,2)=SUM(a(1:2,2:3))/4-h
d_h(2,1)=SUM(a(2:3,1:2))/4-h
d_h(2,2)=SUM(a(2:3,2:3))/4-h
IF (det==1.AND.det_on) THEN
rad=sqrt(2.)*del_xdet*1D3/2
ELSE
rad=sqrt(2.)*del_xi*1D3/2
ENDIF
CALL NEAR
AB=AB+G_cp
!
!Write Bouguer anomaly value
write(15,*)cx,cy,AB
ENDDO
ENDDO
!End of loop for the FA anomaly grid
!
L_x=cont_x*rtio*del_xi
L_y=cont_y*rtio*del_xi
OPEN(1, file='lat_ex.dat', status='UNKNOWN')
WRITE(1,*) L_x,L_y,rtio*del_xi
CLOSE(1)
write(*,*)'*****************************'
write(*,*)'BOUGUER ANOMALY CALCULATED'
write(*,*)'*****************************'
DEALLOCATE (E,FA,FA_strg)
IF (det==1) DEALLOCATE (E_det)
CLOSE(15)
CLOSE(17)
STOP
END
INTEGER FUNCTION up_INT(n)
REAL(dp) :: n
IF (INT(n)-n==0) THEN
up_INT=INT(n)
ELSE
up_INT=INT(n)+1
ENDIF
END FUNCTION
SUBROUTINE EXPAN
!Calculates the vertical atraction of a rigth rectangular
!prism considering an spherical harmonic expansion (Mc Millan, 1958)
!to order 1/r**6, being r the distance between the Mass Centre (MC)
!of the prism and the calculation point
!The result is expresed in mgal=1e-5m/s²
USE MOD_EXPAN
USE MOD_PARAM
REAL(dp) ax,ay,az,ra,Gp,rh
IF (zc==0) THEN
Gp_as=0D0
RETURN
ENDIF
! Distance to MC of the prism (m)
ra=sqrt(xc**2+yc**2+zc**2)
! Parámeters of the power expansion
ax=(2*d_x**2-d_y**2-d_z**2)*xc**2
ay=(2*d_y**2-d_x**2-d_z**2)*yc**2
az=(2*d_z**2-d_y**2-d_x**2)*zc**2 |
|
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 148
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 10:35 pm Post subject: |
|
|
| Code: |
Gp=(G*d_rho*d_x*d_y*d_z*zc)/(ra**3)
Gp_as=(Gp+(ax+ay+az*(1-2*ra**2/(5*zc**2)))*5*G*d_rho* &
d_x*d_y*d_z*zc/(24*ra**7))*1D5
END SUBROUTINE
SUBROUTINE atrac_prisma(x1,x2,y1,y2,z1,z2,d_r,g_z)
USE MOD_PARAM
!Calculates the vertical atraction, g_z, of rigth rectangular
!prism of density d_r, defined by the coordinates of its six
!vertex x1,x2,y1,y2,z1,z2 (Nagy et al., J Geodesy 2000, 74)
!All input units must be referred to SI (m,kg/m³,m/s²).
!Output is in mgal (1 mgal=1e-5 m/s²)
REAL(dp) x1,x2,y1,y2,z1,z2,d_r,x(2),y(2),z(2),v,g_z
INTEGER i,j,k,sgn
x(1)=x1
x(2)=x2
y(1)=y1
y(2)=y2
z(1)=z1
z(2)=z2
v=0
do i=1,2
do j=1,2
do k=1,2
sgn=(-1)**(i+j+k+1)
r=sqrt(x(i)**2+y(j)**2+z(k)**2)
if ((x(i)==0).AND.(y(j)==0).AND.(z(k)==0)) then
g_z=0
elseif ((x(i)==0).AND.(y(j)==0)) then
g_z=-z(k)*atan(x(i)*y(j)/(z(k)*r))*sgn
elseif (((x(i)==0).AND.(z(k)==0))) then
g_z=y(j)*log(abs(y(j)))*sgn
elseif ((y(j)==0).AND.(z(k)==0)) then
g_z=x(i)*log(abs(x(i)))*sgn
elseif ((z(k)==0)) then
g_z=(x(i)*log(y(j)+r)+y(j)*log(x(i)+r))*sgn
else
g_z=x(i)*log(y(j)+r)+y(j)*log(x(i)+r)-z(k)*atan(x(i)*y(j)/(z(k)* &
r))
g_z=g_z*sgn
endif
v=v+g_z
enddo
enddo
enddo
g_z=v*d_r*G*1D5
END SUBROUTINE
SUBROUTINE NEAR
!Calculates the vertical atraction of four
!quadrants of a conic prism which slope continuously from
!the vertex of the near zone to the calculation point.
!If the calculation point is near the coast line,i.e., E and h
!have different signs, the quarter of the conic prism is splitted in
!three parts (see Fullea et al., Comp & Geosc)
!The result is expresed in mgal=1e-5m/s² |
|
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 148
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 10:36 pm Post subject: |
|
|
Final part of the program code - hope it si easy to follow
| Code: |
USE MOD_PARAM
USE MOD_NEAR
REAL(dp) root,g_cp1,g_cp2,g_c,d_r,gdm
G_cp=0
DO i=1,2
DO j=1,2
root=SQRT(rad**2+d_h(i,j)**2)
gdm= pi*G*rad*(root-rad)/(2*root)
IF (h>0.AND.d_h(i,j)+h>0) THEN
d_r=rho_c
G_cp=G_cp+gdm*d_r
ELSEIF(h<0.AND.d_h(i,j)+h<0) THEN
d_r=rho_c-rho_w
IF (d_h(i,j)>0) THEN
G_cp=G_cp-gdm*d_r
ELSE
G_cp=G_cp+gdm*d_r
ENDIF
ELSE
g_cp1=-(h/d_h(i,j))*gdm
g_cp2=(pi*G/(2*d_h(i,j)))*(-rad**2*(d_h(i,j)+h)/root+d_h(i,j)* &
sqrt(rad**2+h**2)+h*root)
g_c=pi*G*(rad*(d_h(i,j)+h)-d_h(i,j)*sqrt(rad**2+h**2) &
-h*root)/(2*d_h(i,j))
IF (h>0.AND.d_h(i,j)+h<0) THEN
G_cp=G_cp+g_cp1*rho_c+g_cp2*(rho_c-rho_w)+g_c*rho_c
ELSEIF (h<0.AND.d_h(i,j)+h>0) THEN
G_cp=G_cp-g_cp1*(rho_c-rho_w)+g_cp2*rho_c-g_c*(rho_c-rho_w)
ENDIF
ENDIF
ENDDO
ENDDO
!Convert to mgals
G_cp=G_cp*1D5
END SUBROUTINE
|
|
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 148
Location: United Kingdom
|
| Posted: Tue Jul 16, 2013 10:40 pm Post subject: |
|
|
Hi all
Apologies for all the code postings, but there is no other way of showing the code and to try and understand why I get the errors on compilation as shown previously.
Copied and pasted together in the sequence they loaded and you get the whole thing. The code is public domain and published, and on closer inspection of the original looks to be a conversion of Fortran77 as I had to correct the continuation syntax and comments.
@John: this includes your code suggestions placed ahead of the main declarations as you indicated.
Hopefully this will make things a little clearer
Cheers
Lester |
|
| Back to top |
|
|
simon
Joined: 05 Jul 2006
Posts: 308
|
| Posted: Wed Jul 17, 2013 12:58 pm Post subject: |
|
|
Your definition of qp is not accepted by the compiler because the numbers are out of range - try printing qp and you should find that it is likely to be a negative number, indicating an error (see the FTN95 help pages on SELECTED_REAL_KIND).
The other error messages are either a result of the invalid value for qp, or a result of the fact that the structure of the program seems to break down about 1/2 of the way through the first posting - there's an END statement at line 70 that appears to be unmatched. At line 52 - is this meant to be the start of a new module or of the PROGRAM? What is the structure after line 70? |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 148
Location: United Kingdom
|
| Posted: Wed Jul 17, 2013 2:15 pm Post subject: Re: |
|
|
| simon wrote: |
Your definition of qp is not accepted by the compiler because the numbers are out of range - try printing qp and you should find that it is likely to be a negative number, indicating an error (see the FTN95 help pages on SELECTED_REAL_KIND).
The other error messages are either a result of the invalid value for qp, or a result of the fact that the structure of the program seems to break down about 1/2 of the way through the first posting - there's an END statement at line 70 that appears to be unmatched. At line 52 - is this meant to be the start of a new module or of the PROGRAM? What is the structure after line 70? |
The structure from line 70 on is defining the variables and inputs that will be read in from the input parameter file, and I guess dimensioning arrays?
The END statement was referenced in John's code, so I just transferred that over, can it be deleted?
The code is public domain Fortran and is a geophysical program that is used to compute the full Bouguer anomaly (including local terrain correction). The reference if you need to check is:
Fullea et al 2008. FA2BOUG -A Fortran 90 code to compute Bouguer gravity anomalies from gridded free-air anomalies: Application to the Atlantic-Mediterranean transition zone. Computers and Geosciences, v34, 1665-1681.
I have tried to compile the code using various compilers, FTN95 and gFortran, and each gives a different response in terms of the errors. GFortran for example, takes exception to the ICHAR statements, whereas FTN95 does not. I am not an expert in Fortran, and why I am trying to see if the code can be compiled, clearly it must be possible because other people in papers published over the last few years have used it. Tried writing to the author, but to date have not had a response.
If needed, I can e-mail both the paper and the original code as published.
Cheers
Lester |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2623
Location: Sydney
|
| Posted: Thu Jul 18, 2013 1:52 am Post subject: |
|
|
Lester,
You should not have placed the code I wrote before your declarations.
There are only a few changes,
1) remove my PI code
2) include USE MOD_Precision in routine up_INT(n)
3) change references to file units 1,2 & 5 to 21,22 & 25, as units 1:8 can be reserved file unit numbers with some compilers.
I see no problem with the definition of sp, dp or qp you now use ?
I copied your code then did my updates and used FC to identify all changes. The following is this list.
| Code: |
Comparing files arctica_ver0.f90 and ARCTICA_VER1.F90
***** arctica_ver0.f90
58: USE MOD_PARAM
59: !
60: pi = 4.0_dp * ATAN (1.0_dp)
61: piq = 4.0_qp * ATAN (1.0_qp)
62: pis = 4.0_sp * ATAN (1.0_sp)
63: write (*,*) pis, kind(pis), precision(pis)
64: write (*,*) pi, kind(pi), precision(pi)
65: write (*,*) piq, kind(piq), precision(piq)
66: write (*,'(f27.19)') pis
67: write (*,'(f27.19)') pi
68: write (*,'(f27.19)') piq
69: !
70: END
71:
72: REAL(dp) del_xd,del_xi,cx,cy,del_xBg,R_d,R_i,del_xdet,del_x_pr
***** ARCTICA_VER1.F90
58: USE MOD_PARAM
59: !
60: REAL(dp) del_xd,del_xi,cx,cy,del_xBg,R_d,R_i,del_xdet,del_x_pr
*****
***** arctica_ver0.f90
90: !
91:
92: !Last edited by arctica on Tue Jul 16, 2013 10:37 pm; edited 1 time in total
***** ARCTICA_VER1.F90
78: !
79: call write_pi
80: !
81: !Last edited by arctica on Tue Jul 16, 2013 10:37 pm; edited 1 time in total
*****
***** arctica_ver0.f90
119:
120:
121: OPEN(1, file='parameters.dat', status='OLD')
122: READ(1,'(I1)',ADVANCE='NO') det
123: IF (det==1) THEN
124: READ(1,*)del_xd,del_xi,del_xBg,rho_c,rho_w,N,M,R_d,R_i,land, &
125: slab_Boug,N_det,M_det,del_xdet
***** ARCTICA_VER1.F90
108:
109: OPEN(21, file='parameters.dat', status='OLD')
110: READ(21,'(I1)',ADVANCE='NO') det
111: IF (det==1) THEN
112: READ(21,*)del_xd,del_xi,del_xBg,rho_c,rho_w,N,M,R_d,R_i,land, &
113: slab_Boug,N_det,M_det,del_xdet
*****
***** arctica_ver0.f90
128: ELSE
129: READ(1,*)del_xd,del_xi,del_xBg,rho_c,rho_w,N,M,R_d,R_i,land, &
130: slab_Boug
***** ARCTICA_VER1.F90
116: ELSE
117: READ(21,*)del_xd,del_xi,del_xBg,rho_c,rho_w,N,M,R_d,R_i,land, &
118: slab_Boug
*****
***** arctica_ver0.f90
131: ENDIF
132: CLOSE(1)
133: !Open INPUT files
134: open(2,file='topo_cart.xyz', status='OLD')
135: open(5,file='gravi_cart.xyz',status='OLD')
136: !Open detailed topography, if necessary
137: IF (det==1) open(7,file='topo_cart_det.xyz', status='OLD')
138: !Open OUTPUT files
***** ARCTICA_VER1.F90
119: ENDIF
120: CLOSE(21)
121: !Open INPUT files
122: open(22,file='topo_cart.xyz', status='OLD')
123: open(25,file='gravi_cart.xyz',status='OLD')
124: !Open detailed topography, if necessary
125: IF (det==1) open(27,file='topo_cart_det.xyz', status='OLD')
126: !Open OUTPUT files
*****
|
|
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2623
Location: Sydney
|
| Posted: Thu Jul 18, 2013 1:54 am Post subject: |
|
|
the rest of the comparison :
| Code: |
***** arctica_ver0.f90
178: DO i=1,N
179: read(2,*) E(i,j)
180: read(5,*) FA_strg(i,j)
181: ENDDO
***** ARCTICA_VER1.F90
166: DO i=1,N
167: read(22,*) E(i,j)
168: read(25,*) FA_strg(i,j)
169: ENDDO
*****
***** arctica_ver0.f90
182: ENDDO
183: REWIND (2)
184: !Read the detailed topography file (if present)
***** ARCTICA_VER1.F90
170: ENDDO
171: REWIND (22)
172: !Read the detailed topography file (if present)
*****
***** arctica_ver0.f90
190: DO i=1,N_det
191: read(7,*) E_det_strg(i,j)
192: n_strg=ICHAR(E_det_strg(i,j))
***** ARCTICA_VER1.F90
178: DO i=1,N_det
179: read(27,*) E_det_strg(i,j)
180: n_strg=ICHAR(E_det_strg(i,j))
*****
***** arctica_ver0.f90
472:
473: OPEN(1, file='lat_ex.dat', status='UNKNOWN')
474: WRITE(1,*) L_x,L_y,rtio*del_xi
475: CLOSE(1)
476: write(*,*)'*****************************'
***** ARCTICA_VER1.F90
460:
461: OPEN(21, file='lat_ex.dat', status='UNKNOWN')
462: WRITE(21,*) L_x,L_y,rtio*del_xi
463: CLOSE(21)
464: write(*,*)'*****************************'
*****
***** arctica_ver0.f90
488: INTEGER FUNCTION up_INT(n)
489: REAL(dp) :: n
***** ARCTICA_VER1.F90
476: INTEGER FUNCTION up_INT(n)
477: USE MOD_Precision
478: REAL(dp) :: n
*****
***** arctica_ver0.f90
***** ARCTICA_VER1.F90
612:
613: subroutine write_pi
614: ! example of kind options
615: USE MOD_Precision
616: real(sp) :: pis
617: real(dp) :: pi
618: real(qp) :: piq
619: !
620: pi = 4.0_dp * ATAN (1.0_dp)
621: piq = 4.0_qp * ATAN (1.0_qp)
622: pis = 4.0_sp * ATAN (1.0_sp)
623: write (*,*) pis, kind(pis), sp, precision(pis)
624: write (*,*) pi, kind(pi), dp, precision(pi)
625: write (*,*) piq, kind(piq), qp, precision(piq)
626: write (*,'(f27.19)') pis
627: write (*,'(f27.19)') pi
628: write (*,'(f27.19)') piq
629: !
630: END
*****
|
This size limitation, or no capability of including, even small, attachments is annoying for us all.
David's comment about changing:
REAL(dp) AB
to
REAL(dp) :: AB
is probably a change you could make generally through the code, as it is the standard.
ICHAR should receive a single character, not a character string. FTN95 allows this as an extension and returns the value of the first character only.
John
Last edited by JohnCampbell on Thu Jul 18, 2013 2:17 am; edited 1 time in total |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2623
Location: Sydney
|
| Posted: Thu Jul 18, 2013 2:04 am Post subject: |
|
|
You might be interested, I once wrote the following program to test the valid values for defining KIND
| Code: |
module kind_precision
integer, parameter :: int_1byte = selected_int_kind (2) ! 1
integer, parameter :: int_2byte = selected_int_kind (4) ! 2
integer, parameter :: int_4byte = selected_int_kind (9) ! 3
integer, parameter :: int_8byte = selected_int_kind (18) ! 4
integer, parameter :: real_float = selected_real_kind (6,37) ! 1
integer, parameter :: real_double = selected_real_kind (15,307) ! 2
integer, parameter :: real_long = selected_real_kind (18,4931) ! 3
end module kind_precision
use kind_precision
!
integer ( kind=int_4byte ) :: i,k, last_k
character kind_name*10
real ( kind=real_double ) :: x
real ( kind=real_long ) :: y
!
! confirm integer precision
write (*,*) ' '
write (*,*) 'Test of integer KIND'
last_k = -99
do i = 30,0,-1
k = selected_int_kind (i)
if (k == last_k) cycle
kind_name = ' undefined'
if ( k == int_1byte ) kind_name = '1 byte'
if ( k == int_2byte ) kind_name = '2 byte'
if ( k == int_4byte ) kind_name = '4 byte'
if ( k == int_8byte ) kind_name = '8 byte'
write (*,*) kind_name, ' precision =', i, ' kind =',k
last_k = k
end do
!
! confirm real precision
write (*,*) ' '
write (*,*) 'Test of Real precision KIND'
last_k = -99
do i = 30,0,-1
k = selected_real_kind (i,1)
if (k == last_k) cycle
kind_name = ' undefined'
if ( k == real_float ) kind_name = '4 byte'
if ( k == real_double ) kind_name = '8 byte'
if ( k == real_long ) kind_name = '10 byte'
write (*,*) kind_name, ' precision =', i, ' kind =',k
last_k = k
end do
!
! confirm real exponent
write (*,*) ' '
write (*,*) 'Test of Real exponent KIND'
last_k = -99
do i = 5000,0,-1
k = selected_real_kind (1,i)
if (k == last_k) cycle
kind_name = ' undefined'
if ( k == real_float ) kind_name = '4 byte'
if ( k == real_double ) kind_name = '8 byte'
if ( k == real_long ) kind_name = '10 byte'
write (*,*) kind_name, ' exponent =', i, ' kind =',k
last_k = k
end do
!
! test real constants
!
x = 1.0e300_real_double ; write (*,*) 'x=',x
y = 1.0e3000_real_long ; write (*,*) 'y=',y
x = 1.0d300 ; write (*,*) 'x=',x
! x = 1.0e300 ; write (*,*) 'x=',x
! y = 1.0e3000 ; write (*,*) 'y=',y
!
end
|
There are so many ways of defining precision, which I think is totally over the top. As I have said, my preference is REAL*8 as it clearly states what is required and you don't have to go searching for what dp is or how it was defined.
John |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 148
Location: United Kingdom
|
| Posted: Thu Jul 18, 2013 1:27 pm Post subject: |
|
|
Thanks for the pointers John, will try the code edits out later.
With the MODULE statement, is it optional to leave off the module name from the END MODULE statement, i.e. is FTN95 flexible on that? Only MOD_Precision has the name at both ends, others do not.
Cheers
Lester |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2623
Location: Sydney
|
| Posted: Thu Jul 18, 2013 1:41 pm Post subject: |
|
|
Lester,
The module name is optional for the END MODULE statement, but I would recommend its use.
I find the program source is easier to navigate where end statements are named, but that is my preference.
John |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 148
Location: United Kingdom
|
| Posted: Thu Jul 18, 2013 6:56 pm Post subject: |
|
|
Hi John,
I added you changes and the errors are greatly reduced, just warnings (especially with pi).
598 SUBROUTINE write_pi
599 ! Example of kind options
600 USE MOD_Precision
601 real(sp) :: pis
602 real(dp) :: pi
603 real(qp) :: piq *************
604
605 pi = 4.0_dp * ATAN (1.0_dp)
606 piq = 4.0_qp * ATAN (1.0_qp) ************
607 pis = 4.0_sp * ATAN (1.0_sp)
608 write (*,*) pis, kind(pis), precision(pis)
609 write (*,*) pi, kind(pi), precision(pi)
610 write (*,*) piq, kind(piq), precision(piq)
611 write (*,'(f27.19)') pis
612 write (*,'(f27.19)') pi
613 write (*,'(f27.19)') piq
614!
615 END
Compiling and linking file: FA2Boug_CGv1D.f95
C:\Temp\Fullea\FA2Boug_CGv1D.F95(63) : warning 197 - Variable X has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(63) : warning 197 - Variable Y has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(63) : warning 197 - Variable Z has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(63) : warning 197 - Variable DEL_G has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(66) : warning 197 - Variable DM has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(68) : warning 197 - Variable UP_INT has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(72) : warning 197 - Variable COORD has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(76) : warning 197 - Variable D_RHO_MHO has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(76) : warning 197 - Variable T has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(76) : warning 197 - Variable ISOS has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(76) : warning 197 - Variable Z_CREF has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(190) : warning 242 - Variable R_TL1 has been given a value but never used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(191) : warning 242 - Variable R_TL2 has been given a value but never used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(267) : warning 298 - Variable PI has been used without being given an initial value *************
C:\Temp\Fullea\FA2Boug_CGv1D.F95(237) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\Fullea\FA2Boug_CGv1D.F95(349) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\Fullea\FA2Boug_CGv1D.F95(469) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\Fullea\FA2Boug_CGv1D.F95(484) : warning 197 - Variable RH has been declared but not used
C:\Temp\Fullea\FA2Boug_CGv1D.F95(485) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\Fullea\FA2Boug_CGv1D.F95(527) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\Fullea\FA2Boug_CGv1D.F95(529) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\Fullea\FA2Boug_CGv1D.F95(531) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\Fullea\FA2Boug_CGv1D.F95(533) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\Fullea\FA2Boug_CGv1D.F95(535) : warning 179 - Comparing floating point quantities for equality may give misleading results
C:\Temp\Fullea\FA2Boug_CGv1D.F95(566) : warning 298 - Variable PI has been used without being given an initial value *************
C:\Temp\Fullea\FA2Boug_CGv1D.F95(603) : error 62 - Invalid KIND specifier
C:\Temp\Fullea\FA2Boug_CGv1D.F95(606) : error 636 - KIND parameter out of range, permitted KINDs are 1, 2, or 3
C:\Temp\Fullea\FA2Boug_CGv1D.F95(606) : error 636 - KIND parameter out of range, permitted KINDs are 1, 2, or 3
Compilation failed.
It does not seem to like the real(qp) setting.
It is almost there!
Cheers
Lester |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 148
Location: United Kingdom
|
| Posted: Thu Jul 18, 2013 9:13 pm Post subject: |
|
|
Hi all,
I made a small change to the code:
SUBROUTINE write_pi
! Example of kind options
USE MOD_Precision
real(sp) :: pis
real(dp) :: pi
! real(qp) :: piq
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
Just commented out the real(qp) and its value piq, and the code then compiled and ran as expected. So I am not sure what it did not like about the entries and the precision - @ JohnC
Would be good to understand what the issue is, as I am going through the code putting commentary so I can get to grips with the flow.
Cheers
Lester |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2623
Location: Sydney
|
| Posted: Fri Jul 19, 2013 1:38 am Post subject: |
|
|
Lester,
to overcome the warning about PI, change write_pi to:
| Code: |
subroutine write_pi
! example of kind options
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), sp, precision(pis)
write (*,*) pi, kind(pi), dp, precision(pi)
write (*,*) piq, kind(piq), qp, precision(piq)
write (*,'(f27.19)') pis
write (*,'(f27.19)') pi
write (*,'(f27.19)') piq
!
END |
The only warning that is important is #298, which implies that PI, that is defined in MOD_PARAM is not being used at line 267 and 566 ?
You should make sure you either call write_pi, or include it's definition at the start.
I really do not understand the comments about qp. The following is a valid definition for real*10, which should be in MOD_Precision
INTEGER, PARAMETER :: qp=SELECTED_REAL_KIND(18,500)
John |
|
| 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
