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