Skip to content

Commit

Permalink
Merge pull request #643 from ketch/reformat_f90
Browse files Browse the repository at this point in the history
Improve formatting of f90 code in classic.
  • Loading branch information
rjleveque authored Jan 17, 2021
2 parents e094e91 + 1f61fc0 commit feab1b5
Show file tree
Hide file tree
Showing 14 changed files with 1,801 additions and 1,851 deletions.
34 changes: 16 additions & 18 deletions examples/shallow_sphere/qcor.f90
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
! =====================================================
subroutine qcor(ixy,i,m,aux,q,maxm,num_eqn,num_ghost,qc)
! =====================================================
! =====================================================
subroutine qcor(ixy,i,m,aux,q,maxm,num_eqn,num_ghost,qc)
! =====================================================

! # Compute the correction term to add to the solution for maintaining
! # the conservation on sphere.
! Compute the correction term to add to the solution for maintaining
! the conservation on sphere.

implicit double precision (a-h,o-z)

Expand All @@ -23,7 +23,7 @@ subroutine qcor(ixy,i,m,aux,q,maxm,num_eqn,num_ghost,qc)
dy = dxcom
endif

! # left interface
! left interface
etxl = aux(in+3,i)
etyl = aux(in+4,i)
etzl = aux(in+5,i)
Expand All @@ -32,7 +32,7 @@ subroutine qcor(ixy,i,m,aux,q,maxm,num_eqn,num_ghost,qc)
enyl = aux(in+1,i) * gammal
enzl = aux(in+2,i) * gammal

! # right interface
! right interface
etxr = aux(in+3,i+1)
etyr = aux(in+4,i+1)
etzr = aux(in+5,i+1)
Expand All @@ -41,23 +41,21 @@ subroutine qcor(ixy,i,m,aux,q,maxm,num_eqn,num_ghost,qc)
enyr = aux(in+1,i+1) * gammar
enzr = aux(in+2,i+1) * gammar

qc(1) = (enxr-enxl)*q(2,i) &
+(enyr-enyl)*q(3,i) &
+(enzr-enzl)*q(4,i)
qc(1) = (enxr-enxl)*q(2,i) + (enyr-enyl)*q(3,i) +(enzr-enzl)*q(4,i)

qc(2) = (enxr-enxl)*(q(2,i)**2.d0/q(1,i)+0.5d0*g*q(1,i)**2.d0) &
+(enyr-enyl)*(q(2,i)*q(3,i)/q(1,i)) &
+(enzr-enzl)*(q(2,i)*q(4,i)/q(1,i))
+(enyr-enyl)*(q(2,i)*q(3,i)/q(1,i)) &
+(enzr-enzl)*(q(2,i)*q(4,i)/q(1,i))

qc(3) = (enxr-enxl)*(q(2,i)*q(3,i)/q(1,i)) &
+(enyr-enyl)*(q(3,i)**2.d0/q(1,i)+0.5d0*g*q(1,i)**2.d0) &
+(enzr-enzl)*(q(3,i)*q(4,i)/q(1,i))
+(enyr-enyl)*(q(3,i)**2.d0/q(1,i)+0.5d0*g*q(1,i)**2.d0) &
+(enzr-enzl)*(q(3,i)*q(4,i)/q(1,i))

qc(4) = (enxr-enxl)*(q(2,i)*q(4,i)/q(1,i)) &
+(enyr-enyl)*(q(3,i)*q(4,i)/q(1,i)) &
+(enzr-enzl)*(q(4,i)**2.d0/q(1,i)+0.5d0*g*q(1,i)**2.d0)
+(enyr-enyl)*(q(3,i)*q(4,i)/q(1,i)) &
+(enzr-enzl)*(q(4,i)**2.d0/q(1,i)+0.5d0*g*q(1,i)**2.d0)

! # project qc to tangent plane:
! project qc to tangent plane:
erx = aux(14,i)
ery = aux(15,i)
erz = aux(16,i)
Expand All @@ -69,4 +67,4 @@ subroutine qcor(ixy,i,m,aux,q,maxm,num_eqn,num_ghost,qc)
qc(4) = qc(4) - qcn*erz

return
end subroutine qcor
end subroutine qcor
57 changes: 29 additions & 28 deletions examples/shallow_sphere/qinit.f90
Original file line number Diff line number Diff line change
@@ -1,17 +1,18 @@
! =====================================================
subroutine qinit(maxmx,maxmy,num_eqn,num_ghost,mx,my,xlower,ylower, &
dx,dy,q,num_aux,aux,Rsphere)
! =====================================================
! =====================================================
subroutine qinit(maxmx,maxmy,num_eqn,num_ghost,mx,my,xlower,ylower, &
dx,dy,q,num_aux,aux,Rsphere)
! =====================================================

! # Set initial conditions for q.
! Set initial conditions for q.

! # -------4-Rossby-Haurwitz wave-----------------------
! -------4-Rossby-Haurwitz wave-----------------------

implicit double precision (a-h,o-z)
dimension q(num_eqn, 1-num_ghost:maxmx+num_ghost, 1-num_ghost:maxmy+num_ghost)
dimension aux(num_aux, 1-num_ghost:maxmx+num_ghost, 1-num_ghost:maxmy+num_ghost)
double precision :: Uin(3),Uout(3)
double precision :: K

!f2py integer intent(in) maxmx
!f2py integer intent(in) maxmy
!f2py integer optional,intent(in) num_eqn
Expand All @@ -29,7 +30,7 @@ subroutine qinit(maxmx,maxmy,num_eqn,num_ghost,mx,my,xlower,ylower, &

pi = 4.d0*datan(1.d0)

! # Parameters
! Parameters
a = 6.37122d6
K = 7.848d-6
Omega = 7.292d-5
Expand All @@ -38,12 +39,12 @@ subroutine qinit(maxmx,maxmy,num_eqn,num_ghost,mx,my,xlower,ylower, &
h0 = 8.d3
R = 4.d0

do 20 i=1,mx
do i=1,mx
xc = xlower + (i-0.5d0)*dx
do 20 j=1,my
do j=1,my
yc = ylower + (j-0.5d0)*dy
call mapc2p(xc,yc,xp,yp,zp,Rsphere)
! # compute longitude theta from positive x axis:
! compute longitude theta from positive x axis:
rad = dmax1(dsqrt(xp**2 + yp**2),1.d-6)

if(xp > 0.d0 .AND. yp > 0.d0) then
Expand All @@ -56,7 +57,7 @@ subroutine qinit(maxmx,maxmy,num_eqn,num_ghost,mx,my,xlower,ylower, &
theta = -dasin(-yp/rad)
endif

! # compute phi, at north pole: pi/2 at south pool: -pi/2
! compute phi, at north pole: pi/2 at south pool: -pi/2
if (zp > 0.d0) then
phi = dasin(zp/Rsphere)
else
Expand All @@ -67,41 +68,41 @@ subroutine qinit(maxmx,maxmy,num_eqn,num_ghost,mx,my,xlower,ylower, &
yp = phi

bigA = 0.5d0*K*(2.d0*Omega+K)*dcos(yp)**2.d0 + &
0.25d0*K*K*dcos(yp)**(2.d0*R)*( &
(1.d0*R+1.d0)*dcos(yp)**2.d0 &
+(2.d0*R*R - 1.d0*R - 2.d0) &
- 2.d0*R*R*(dcos(yp))**(-2.d0))
0.25d0*K*K*dcos(yp)**(2.d0*R)*( &
(1.d0*R+1.d0)*dcos(yp)**2.d0 &
+(2.d0*R*R - 1.d0*R - 2.d0) &
- 2.d0*R*R*(dcos(yp))**(-2.d0))
bigB = (2.d0*(Omega+K)*K)/((1.d0*R+1.d0)*(1.d0*R+2.d0)) &
*dcos(yp)**R*( (1.d0*R*R + 2.d0*R + 2.d0) &
- (1.d0*R+1.d0)**(2)*dcos(yp)**2 )
*dcos(yp)**R*( (1.d0*R*R + 2.d0*R + 2.d0) &
- (1.d0*R+1.d0)**(2)*dcos(yp)**2 )
bigC = 0.25d0*K*K*dcos(yp)**(2*R)*( (1.d0*R + 1.d0)* &
dcos(yp)**2 - (1.d0*R + 2.d0))
dcos(yp)**2 - (1.d0*R + 2.d0))

! # longitude (angular) velocity component
! longitude (angular) velocity component
Uin(1) = (K*dcos(yp)+K*dcos(yp)**(R-1.)*( R*dsin(yp)**2. &
- dcos(yp)**2. )*dcos(R*xp))*t0
- dcos(yp)**2. )*dcos(R*xp))*t0

! # latitude (angular) velocity component
! latitude (angular) velocity component
Uin(2) = (-K*R*dcos(yp)**(R-1.)*dsin(yp)*dsin(R*xp))*t0

! # radial velocity component
! radial velocity component
Uin(3) = 0.d0


! # calculate velocity vetor in cartesian coordinates
! calculate velocity vetor in cartesian coordinates
Uout(1) = (-dsin(xp)*Uin(1)-dsin(yp)*dcos(xp)*Uin(2))
Uout(2) = (dcos(xp)*Uin(1)-dsin(yp)*dsin(xp)*Uin(2))
Uout(3) = dcos(yp)*Uin(2)

! # set the clawpack initial values:
! set the clawpack initial values:
q(1,i,j) = h0/a + (a/G)*( bigA + bigB*dcos(R*xp) &
+ bigC*dcos(2.d0*R*xp))
+ bigC*dcos(2.d0*R*xp))
q(2,i,j) = q(1,i,j)*Uout(1)
q(3,i,j) = q(1,i,j)*Uout(2)
q(4,i,j) = q(1,i,j)*Uout(3)


20 END DO
end do
end do

return
end subroutine qinit
end subroutine qinit
Loading

0 comments on commit feab1b5

Please sign in to comment.