Skip to content

Commit

Permalink
final fix in the merged version of geos/outPV
Browse files Browse the repository at this point in the history 
- also remove the deprecated outPV3 from the sphinx doc
  • Loading branch information
@tgastine
tgastine committed Feb 8, 2017
1 parent bb8c2d3 commit 9147b74
Show file tree
Hide file tree
Showing 2 changed files with 42 additions and 45 deletions.
10 changes: 3 additions & 7 deletions doc/sphinx/apiFortran/IOAdd.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,15 +55,11 @@ IO: RMS force balance, torsional oscillations, misc

.. f:automodule:: radial_spectra
``Egeos.f90``
-------------

.. f:automodule:: egeos_mod
``outGeos.f90``
---------------

``outPV3.f90``
--------------
.. f:automodule:: geos_mod
.. f:automodule:: outpv3
``chebInt.f90``
---------------
Expand Down
77 changes: 39 additions & 38 deletions src/outGeos.f90
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,8 @@
module geos_mod
!
! This module is used to calculate the outputs when either l_par
! or l_PV are set to .true.
!

use precision_mod
use parallel_mod
Expand Down Expand Up @@ -26,16 +30,16 @@ module geos_mod
private

real(cp) :: timeOld
real(cp), allocatable :: PlmS(:,:,:), PlmZ(:,:,:)
real(cp), allocatable :: dPlmS(:,:,:), dPlmZ(:,:,:)
real(cp), allocatable :: OsinTS(:,:)
real(cp), allocatable :: PlmS(:,:,:), PlmZ(:,:,:), PlmSPV(:,:,:)
real(cp), allocatable :: dPlmS(:,:,:), dPlmZ(:,:,:), dPlmSPV(:,:,:)
real(cp), allocatable :: OsinTS(:,:), OsinTSPV(:,:)

complex(cp), allocatable :: wS_global(:,:), dwS_global(:,:), ddwS_global(:,:)
complex(cp), allocatable :: zS_global(:,:), dzS_global(:,:)

type(costf_odd_t), allocatable :: chebt_Z(:)
integer, allocatable :: nZmaxS(:), nZC(:), nZC_Sloc(:), nZ2(:,:)
real(cp), allocatable :: zZ(:,:), rZ(:,:)
real(cp), allocatable :: zZ(:,:), rZ(:,:), rZPV(:,:)
real(cp), allocatable :: VorOld(:,:,:)
real(cp), parameter :: eps = 10.0_cp*epsilon(one)

Expand Down Expand Up @@ -71,22 +75,20 @@ subroutine initialize_geos_mod(l_geos, l_PV)
end if
nS_on_last_rank = nS_per_rank + nS_remaining

allocate( OsinTS(nZmaxA/2+1,nSstart:nSstop) )
allocate( PlmS(lm_max,nZmaxA/2+1,nSstart:nSstop) )
allocate( dPlmS(lm_max,nZmaxA/2+1,nSstart:nSstop) )
allocate( zZ(nZmaxA,nSstart:nSstop) )
allocate( rZ(nZmaxA,nSstart:nSstop) )
bytes_allocated = bytes_allocated+ ((nZmaxA/2+1)*(nSstop-nSstart+1)* &
& (1+2*lm_max))*SIZEOF_DEF_REAL + &
& 2*nZmaxA*(nSstop-nSstart+1)*SIZEOF_DEF_REAL

!-- The following global arrays are required in getDVptr
allocate( wS_global(lm_max,n_r_max), dwS_global(lm_max,n_r_max) )
allocate( ddwS_global(lm_max,n_r_max), zS_global(lm_max,n_r_max) )
allocate( dzS_global(lm_max,n_r_max) )
bytes_allocated = bytes_allocated+ 5*lm_max*n_r_max*SIZEOF_DEF_COMPLEX

if ( l_PV ) then
allocate( OsinTSPV(nZmaxA/2+1,nSstart:nSstop) )
allocate( PlmSPV(lm_max,nZmaxA/2+1,nSstart:nSstop) )
allocate( dPlmSPV(lm_max,nZmaxA/2+1,nSstart:nSstop) )
allocate( rZPV(nZmaxA,nSstart:nSstop) )
bytes_allocated = bytes_allocated+ ((nZmaxA/2+1)*(nSstop-nSstart+1)* &
& (1+2*lm_max))*SIZEOF_DEF_REAL + &
& nZmaxA*(nSstop-nSstart+1)*SIZEOF_DEF_REAL
allocate( PlmZ(l_max+1,nZmaxA/2+1,nSstart:nSstop) )
allocate( dPlmZ(l_max+1,nZmaxA/2+1,nSstart:nSstop) )
bytes_allocated = bytes_allocated + 2*(l_max+1)*(nZmaxA/2+1)* &
Expand All @@ -107,6 +109,15 @@ subroutine initialize_geos_mod(l_geos, l_PV)
end if

if ( l_geos ) then
allocate( OsinTS(nZmaxA/2+1,nSstart:nSstop) )
allocate( PlmS(lm_max,nZmaxA/2+1,nSstart:nSstop) )
allocate( dPlmS(lm_max,nZmaxA/2+1,nSstart:nSstop) )
allocate( zZ(nZmaxA,nSstart:nSstop) )
allocate( rZ(nZmaxA,nSstart:nSstop) )
bytes_allocated = bytes_allocated+ ((nZmaxA/2+1)*(nSstop-nSstart+1)* &
& (1+2*lm_max))*SIZEOF_DEF_REAL + &
& 2*nZmaxA*(nSstop-nSstart+1)*SIZEOF_DEF_REAL

allocate( nZmaxS(nSstart:nSstop) )
allocate( chebt_Z(nSstart:nSstop) )
bytes_allocated = bytes_allocated+(nSstop-nSstart+1)*SIZEOF_INTEGER
Expand All @@ -126,15 +137,18 @@ subroutine finalize_geos_mod(l_geos, l_PV)
logical, intent(in) :: l_geos
logical, intent(in) :: l_PV

deallocate( OsinTS, PlmS, dPlmS, zZ, rZ )
deallocate( wS_global, dwS_global, ddwS_global, zS_global, dzS_global)

if ( l_geos ) then
deallocate( OsinTS, PlmS, dPlmS, zZ, rZ )
deallocate( chebt_Z, nZmaxS )
if ( rank == 0 .and. (.not. l_save_out) ) close(n_geos_file)
end if

if ( l_PV ) deallocate( PlmZ, dPlmZ, VorOld, nZC_Sloc, nZC, nZ2 )
if ( l_PV ) then
deallocate( PlmZ, dPlmZ, VorOld, nZC_Sloc, nZC, nZ2 )
deallocate( OsinTSPV, PlmSPV, dPlmSPV, rZPV )
end if

end subroutine finalize_geos_mod
!----------------------------------------------------------------------------
Expand Down Expand Up @@ -171,7 +185,6 @@ subroutine getEgeos(time,nGeosSets,w,dw,ddw,z,dz,Geos,dpFlow,dzFlow)
real(cp) :: zNorm ! Norm z interval
integer :: nNorm ! No. of grid points for norm interval
real(cp) :: zMin,zMax,help ! integration boundarie, help variable
logical :: lAS ! .true. if axisymmetric (m=0) functions
real(cp) :: sZ(nSmax),dsZ ! cylindrical radius s and s-step
integer :: nPhi,nI
real(cp) :: phiNorm
Expand All @@ -195,8 +208,6 @@ subroutine getEgeos(time,nGeosSets,w,dw,ddw,z,dz,Geos,dpFlow,dzFlow)
real(cp) :: dpEk_s(nSstart:nSstop),dzEk_s(nSstart:nSstop)
real(cp) :: thetaZ

logical :: lStopRun

!-- Correlation (new Oct. 4 2007)
logical :: lCorrel
real(cp) :: VzS,VzN,VorS,VorN,surf,delz
Expand Down Expand Up @@ -227,17 +238,14 @@ subroutine getEgeos(time,nGeosSets,w,dw,ddw,z,dz,Geos,dpFlow,dzFlow)
call costf_arrays(w,dw,ddw,z,dz)

lCorrel=.true. ! Calculate Vz and Vorz north/south correlation

phiNorm=two*pi/n_phi_max
lStopRun=.false.
lDeriv=.true.
kindCalc=1

!---- Get resolution in s and z for z-integral:
zNorm=one ! This is r_CMB-r_ICB
nNorm=int(zDens*n_r_max) ! Covered with nNorm points !
dsZ =r_CMB/real(nSmax,cp) ! Step in s controlled by nSmax
lAS=.false.

!-- Global array
do nS=1,nSmax
Expand Down Expand Up @@ -296,11 +304,6 @@ subroutine getEgeos(time,nGeosSets,w,dw,ddw,z,dz,Geos,dpFlow,dzFlow)
! all together nZmaxS(nS) from
! south to north including equator
end if
if ( 2*nZmax > nZmaxA ) then ! TC case most critical
write(*,*) '! nZmaxA too small in getEgeos!'
write(*,*) '! Should be at least:',2*nZmax
lStopRun=.true.
end if
do nZ=1,nZmax
rZ(nZ,nS) =sqrt(zZ(nZ,nS)**2+sZ(nS)**2)
thetaZ =atan2(sZ(nS),zZ(nZ,nS))
Expand Down Expand Up @@ -727,8 +730,8 @@ subroutine outPV(time,l_stop_time,nPVsets,w,dw,ddw,z,dz,omega_IC,omega_MA)
real(cp) :: fac

!--- define Grid
integer :: nS,nSI
real(cp) :: sZ(nSmax),dsZ ! cylindrical radius s and s-step
integer :: nS
real(cp) :: sZ(nSmax),dsZ ! cylindrical radius s and s-step

integer :: nZ,nZmaxNS
integer, save :: nZS
Expand Down Expand Up @@ -796,10 +799,8 @@ subroutine outPV(time,l_stop_time,nPVsets,w,dw,ddw,z,dz,omega_IC,omega_MA)

!-- Global arrays
dsZ=r_CMB/real(nSmax,kind=cp) ! Step in s controlled by nSmax
nSI=0 ! Inner core position
do nS=1,nSmax
sZ(nS)=(nS-half)*dsZ
if ( sZ(nS) < r_ICB .and. nS > nSI ) nSI=nS
end do
zstep=2*r_CMB/real(nZmaxA-1,kind=cp)
do nZ=1,nZmaxA
Expand All @@ -818,14 +819,14 @@ subroutine outPV(time,l_stop_time,nPVsets,w,dw,ddw,z,dz,omega_IC,omega_MA)
nZC_Sloc(nS)=nZC_Sloc(nS)+1 ! Counts all z within shell
nZ2(nZ,nS)=nZC_Sloc(nS) ! No of point within shell
if ( zZ(nZ) > 0 ) then ! Onl north hemisphere
rZ(nZC_Sloc(nS),nS)=rZS
rZPV(nZC_Sloc(nS),nS)=rZS
thetaZ=atan2(sZ(nS),zZ(nZ))
OsinTS(nZC_Sloc(nS),nS)=one/sin(thetaZ)
OsinTSPV(nZC_Sloc(nS),nS)=one/sin(thetaZ)
call plm_theta(thetaZ,l_max,0,minc,PlmZ(:,nZC_Sloc(nS),nS),&
& dPlmZ(:,nZC_Sloc(nS),nS),l_max+1,2)
call plm_theta(thetaZ,l_max,m_max,minc, &
& PlmS(:,nZC_Sloc(nS),nS), &
& dPlmS(:,nZC_Sloc(nS),nS),lm_max,2)
call plm_theta(thetaZ,l_max,m_max,minc, &
& PlmSPV(:,nZC_Sloc(nS),nS), &
& dPlmSPV(:,nZC_Sloc(nS),nS),lm_max,2)
end if
else
nZ2(nZ,nS)=-1 ! No z found within shell !
Expand All @@ -837,7 +838,7 @@ subroutine outPV(time,l_stop_time,nPVsets,w,dw,ddw,z,dz,omega_IC,omega_MA)
nZmaxNS=nZC_Sloc(nS) ! all z points within shell
if ( l_stop_time ) then
call getPAStr(VpAS,dzVpLMr,nZmaxNS,nZmaxA,l_max,r_ICB,r_CMB,n_r_max, &
& rZ(:,nS),dPlmZ(:,:,nS),OsinTS(:,nS))
& rZPV(:,nS),dPlmZ(:,:,nS),OsinTSPV(:,nS))

!-- Copy to array with all z-points
do nZ=1,nZmaxA
Expand All @@ -857,8 +858,8 @@ subroutine outPV(time,l_stop_time,nPVsets,w,dw,ddw,z,dz,omega_IC,omega_MA)

!-- Get all three components in the shell
call getDVptr(wS_global,dwS_global,ddwS_global,zS_global,dzS_global, &
& r_ICB,r_CMB,rZ(:,nS),nZmaxNS,nZmaxA,PlmS(:,:,nS), &
& dPlmS(:,:,nS),OsinTS(:,nS),kindCalc,VsS,VzS,VpS,VorS)
& r_ICB,r_CMB,rZPV(:,nS),nZmaxNS,nZmaxA,PlmSPV(:,:,nS), &
& dPlmSPV(:,:,nS),OsinTSPV(:,nS),kindCalc,VsS,VzS,VpS,VorS)

if ( l_stop_time ) then
nC=0
Expand Down

0 comments on commit 9147b74

Please sign in to comment.