nrutil.F

Go to the documentation of this file.

        MODULE nrutil
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS 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