MODULE nrutil
!
!git $Id$
!svn $Id: nrutil.F 1151 2023-02-09 03:08:53Z arango $
!================================================== Hernan G. Arango ===
! Copyright (c) 2002-2023 The ROMS/TOMS Group !
! Licensed under a MIT/X style license !
! See License_ROMS.md !
!=======================================================================
! !
! Numerical Recipies Utility. !
! !
! Adapted from Numerical Recepies. !
! !
! Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, !
! 1996: Numerical Recipes in Fortran 90, The Art of Parallel !
! Scientific Computing, 2nd Edition, Cambridge Univ. Press. !
! !
!=======================================================================
!
USE mod_kinds
implicit none
PUBLIC
integer(i8b), parameter :: NPAR_ARTH = 16
integer(i8b), parameter :: NPAR2_ARTH = 8
INTERFACE array_copy
MODULE PROCEDURE array_copy_r, array_copy_d, array_copy_i
END INTERFACE
INTERFACE arth
MODULE PROCEDURE arth_r, arth_d, arth_i
END INTERFACE
INTERFACE reallocate
MODULE PROCEDURE reallocate_rv, reallocate_rm, &
& reallocate_iv, reallocate_im, &
& reallocate_hv
END INTERFACE
INTERFACE gasdev
SUBROUTINE gasdev_s (harvest)
USE mod_kinds
real(r8), intent(out) :: harvest
END SUBROUTINE gasdev_s
!
SUBROUTINE gasdev_v (harvest)
USE mod_kinds
real(r8), dimension(:), intent(out) :: harvest
END SUBROUTINE gasdev_v
END INTERFACE
INTERFACE ran1
SUBROUTINE ran1_s (harvest)
USE mod_kinds
real(r8), intent(out) :: harvest
END SUBROUTINE ran1_s
!
SUBROUTINE ran1_v (harvest)
USE mod_kinds
real(r8), dimension(:), intent(out) :: harvest
END SUBROUTINE ran1_v
END INTERFACE
CONTAINS
SUBROUTINE array_copy_r (src, dest, n_copied, n_not_copied)
!
!=======================================================================
! !
! Copy single precision array where size of source not known in !
! advance. !
! !
!=======================================================================
!
! Imported variable declarations.
!
real(r4), intent(in) :: src(:)
real(r4), intent(out) :: dest(:)
integer(i8b), intent(out) :: n_copied, n_not_copied
!
!-----------------------------------------------------------------------
! Copy single precision array.
!-----------------------------------------------------------------------
!
n_copied=MIN(SIZE(src), SIZE(dest))
n_not_copied=SIZE(src)-n_copied
dest(1:n_copied)=src(1:n_copied)
RETURN
END SUBROUTINE array_copy_r
SUBROUTINE array_copy_d (src, dest, n_copied, n_not_copied)
!
!=======================================================================
! !
! Copy double precision array where size of source not known in !
! advance. !
! !
!=======================================================================
!
! Imported variable declarations.
!
real(dp), intent(in) :: src(:)
real(dp), intent(out) :: dest(:)
integer(i8b), intent(out) :: n_copied, n_not_copied
!
!-----------------------------------------------------------------------
! Copy double precision array.
!-----------------------------------------------------------------------
!
n_copied=MIN(SIZE(src), SIZE(dest))
n_not_copied=size(src)-n_copied
dest(1:n_copied)=src(1:n_copied)
RETURN
END SUBROUTINE array_copy_d
SUBROUTINE array_copy_i (src, dest, n_copied, n_not_copied)
!
!=======================================================================
! !
! Copy integer array where size of source not known in advance. !
! !
!=======================================================================
!
! Imported variable declarations.
!
integer(i8b), intent(in) :: src(:)
integer(i8b), intent(out) :: dest(:)
integer(i8b), intent(out) :: n_copied, n_not_copied
!
!-----------------------------------------------------------------------
! Copy integer array.
!-----------------------------------------------------------------------
!
n_copied=min(size(src),size(dest))
n_not_copied=size(src)-n_copied
dest(1:n_copied)=src(1:n_copied)
RETURN
END SUBROUTINE array_copy_i
FUNCTION arth_r (first, increment, n)
!
!=======================================================================
! !
! Array function returning an arithmetic progression, single !
! precision. !
! !
!=======================================================================
!
! Imported variable declarations.
!
integer(i8b), intent(in) :: n
real(r4), intent(in) :: first, increment
real(r4), dimension(n) :: arth_r
!
! Local variable declarations.
!
integer(i8b) :: k, k2
real(r4) :: temp
!
!----------------------------------------------------------------------
! Set arithmetic progression.
!----------------------------------------------------------------------
!
IF (n.gt.0) arth_r(1)=first
IF (n.le.NPAR_ARTH) THEN
DO k=2,n
arth_r(k)=arth_r(k-1)+increment
END DO
ELSE
DO k=2,NPAR2_ARTH
arth_r(k)=arth_r(k-1)+increment
END DO
temp=increment*NPAR2_ARTH
k=NPAR2_ARTH
DO
IF (k.ge.n) EXIT
k2=k+k
arth_r(k+1:MIN(k2,n))=temp+arth_r(1:MIN(k,n-k))
temp=temp+temp
k=k2
END DO
END IF
RETURN
END FUNCTION arth_r
FUNCTION arth_d (first, increment, n)
!
!=======================================================================
! !
! Array function returning an arithmetic progression, double !
! precision. !
! !
!=======================================================================
!
! Imported variable declarations.
!
integer(i8b), intent(in) :: n
real(dp), intent(in) :: first, increment
real(dp), dimension(n) :: arth_d
!
! Local variable declarations.
!
integer(i8b) :: k, k2
real(dp) :: temp
!
!----------------------------------------------------------------------
! Set arithmetic progression.
!----------------------------------------------------------------------
!
IF (n.gt.0) arth_d(1)=first
IF (n.le.NPAR_ARTH) THEN
DO k=2,n
arth_d(k)=arth_d(k-1)+increment
END DO
ELSE
DO k=2,NPAR2_ARTH
arth_d(k)=arth_d(k-1)+increment
END DO
temp=increment*NPAR2_ARTH
k=NPAR2_ARTH
DO
IF (k.ge.n) EXIT
k2=k+k
arth_d(k+1:MIN(k2,n))=temp+arth_d(1:MIN(k,n-k))
temp=temp+temp
k=k2
END DO
END IF
RETURN
END FUNCTION arth_d
FUNCTION arth_i (first, increment, n)
!
!=======================================================================
! !
! Integer array function returning an arithmetic progression. !
! !
!=======================================================================
!
! Imported variable declarations.
!
integer(i8b), intent(in) :: first, increment, n
integer(i8b), dimension(n) :: arth_i
!
! Local variable declarations.
!
integer(i8b) :: k, k2, temp
!
!----------------------------------------------------------------------
! Set arithmetic progression.
!----------------------------------------------------------------------
!
IF (n.gt.0) arth_i(1)=first
IF (n.le.NPAR_ARTH) THEN
DO k=2,n
arth_i(k)=arth_i(k-1)+increment
END DO
ELSE
DO k=2,NPAR2_ARTH
arth_i(k)=arth_i(k-1)+increment
END DO
temp=increment*NPAR2_ARTH
k=NPAR2_ARTH
DO
IF (k.ge.n) EXIT
k2=k+k
arth_i(k+1:MIN(k2,n))=temp+arth_i(1:MIN(k,n-k))
temp=temp+temp
k=k2
END DO
END IF
RETURN
END FUNCTION arth_i
SUBROUTINE nrerror (string)
!
!=======================================================================
! !
! Report an error message and the die. !
! !
!=======================================================================
!
USE mod_iounits, ONLY : stdout
!
! Imported variable declarations.
!
character(len=*), intent(in) :: string
!
!-----------------------------------------------------------------------
! Report error message to standard output and terminate execution.
!-----------------------------------------------------------------------
!
WRITE (stdout,10) string, 'program terminated by NRERROR'
10 FORMAT (/,1x,a,/20x,a)
STOP
END SUBROUTINE nrerror
FUNCTION reallocate_rv (p, n)
!
!=======================================================================
! !
! Reallocate a pointer of a single precision vector to a new size, !
! preserving its previous content. !
! !
!=======================================================================
!
! Imported variable declarations.
!
real(r4), pointer :: p(:)
real(r4), pointer :: reallocate_rv(:)
integer(i8b), intent(in) :: n
!
! Local variable declarations.
!
integer(i8b) :: nold, ierr
!
!-----------------------------------------------------------------------
! Reallocate pointer for a single precision vector.
!-----------------------------------------------------------------------
!
ALLOCATE (reallocate_rv(n), STAT=ierr)
IF (ierr.ne.0) &
& CALL nrerror ('REALLOCATE_RV: error while allocating memory')
IF (.not.ASSOCIATED(p)) RETURN
nold=SIZE(p)
reallocate_rv(1:MIN(nold,n))=p(1:MIN(nold,n))
DEALLOCATE (p)
RETURN
END FUNCTION reallocate_rv
FUNCTION reallocate_iv (p,n)
!
!=======================================================================
! !
! Reallocate a pointer of a integer vector to a new size, preserving !
! its previous content. !
! !
!=======================================================================
!
! Imported variable declarations.
!
integer(i8b), pointer :: p(:)
integer(i8b), pointer :: reallocate_iv(:)
integer(i8b), intent(in) :: n
!
! Local variable declarations.
!
integer(i8b) :: nold, ierr
!
!-----------------------------------------------------------------------
! Reallocate pointer for a integer vector.
!-----------------------------------------------------------------------
!
ALLOCATE (reallocate_iv(n), STAT=ierr)
IF (ierr.ne.0) &
& CALL nrerror ('REALLOCATE_IV: error while allocating memory')
IF (.not.ASSOCIATED(p)) RETURN
nold=SIZE(p)
reallocate_iv(1:MIN(nold,n))=p(1:MIN(nold,n))
DEALLOCATE (p)
RETURN
END FUNCTION reallocate_iv
FUNCTION reallocate_hv (p, n)
!
!=======================================================================
! !
! Reallocate a pointer of a character vector to a new size, !
! preserving its previous content. !
! !
!=======================================================================
!
! Imported variable declarations.
!
character (len=1), pointer :: p(:)
character (len=1), pointer :: reallocate_hv(:)
integer(i8b), intent(in) :: n
!
! Local variable declarations.
!
integer(i8b) :: nold, ierr
!
!-----------------------------------------------------------------------
! Reallocate pointer for a integer vector.
!-----------------------------------------------------------------------
!
ALLOCATE (reallocate_hv(n),stat=ierr)
IF (ierr.ne.0) &
& CALL nrerror ('REALLOCATE_HV: error while allocating memory')
IF (.not.ASSOCIATED(p)) RETURN
nold=SIZE(p)
reallocate_hv(1:MIN(nold,n))=p(1:MIN(nold,n))
DEALLOCATE (p)
RETURN
END FUNCTION reallocate_hv
FUNCTION reallocate_rm (p, n, m)
!
!=======================================================================
! !
! Reallocate a pointer of a single precision matrix to a new size, !
! preserving its previous content. !
! !
!=======================================================================
!
! Imported variable declarations.
!
real(r4), pointer :: p(:,:)
real(r4), pointer :: reallocate_rm(:,:)
integer(i8b), intent(in) :: n, m
!
! Local variable declarations.
!
integer(i8b) :: nold, mold, ierr
!
!-----------------------------------------------------------------------
! Reallocate pointer for a single precision matrix.
!-----------------------------------------------------------------------
!
ALLOCATE (reallocate_rm(n,m), STAT=ierr)
IF (ierr.ne.0) &
& CALL nrerror ('REALLOCATE_RM: error while allocating memory')
IF (.not.ASSOCIATED(p)) RETURN
nold=SIZE(p,1)
mold=SIZE(p,2)
reallocate_rm(1:MIN(nold,n),1:MIN(mold,m))= &
& p(1:MIN(nold,n),1:MIN(mold,m))
DEALLOCATE (p)
RETURN
END FUNCTION reallocate_rm
FUNCTION reallocate_im (p, n, m)
!
!=======================================================================
! !
! Reallocate a pointer of a integer matrix to a new size, preserving !
! its previous content. !
! !
!=======================================================================
!
! Imported variable declarations.
!
integer(i8b), pointer :: p(:,:)
integer(i8b), pointer :: reallocate_im(:,:)
integer(i8b), intent(in) :: n, m
!
! Local variable declarations.
!
integer(i8b) :: nold, mold, ierr
!
!-----------------------------------------------------------------------
! Reallocate pointer for a integer matrix.
!-----------------------------------------------------------------------
!
ALLOCATE (reallocate_im(n,m), STAT=ierr)
IF (ierr.ne.0) &
& CALL nrerror ('REALLOCATE_IM: error while allocating memory')
IF (.not.ASSOCIATED(p)) RETURN
nold=SIZE(p,1)
mold=SIZE(p,2)
reallocate_im(1:MIN(nold,n),1:MIN(mold,m))= &
& p(1:MIN(nold,n),1:MIN(mold,m))
DEALLOCATE (p)
RETURN
END FUNCTION reallocate_im
END MODULE nrutil