From 14376099fccc36b9525d0e2a40779fa0a4fbd4d2 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 12 Dec 2019 16:02:35 +0000 Subject: [PATCH 1/5] Correct spelling in HI doxumentation --- src/framework/MOM_hor_index.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/framework/MOM_hor_index.F90 b/src/framework/MOM_hor_index.F90 index 2fda7bd68d..8f7ce1db8e 100644 --- a/src/framework/MOM_hor_index.F90 +++ b/src/framework/MOM_hor_index.F90 @@ -113,7 +113,7 @@ end subroutine HIT_assign !> \namespace mom_hor_index !! -!! The hor_index_type provides the decalarations and loop ranges for almost all data with horizontal extent. +!! The hor_index_type provides the declarations and loop ranges for almost all data with horizontal extent. !! !! Declarations and loop ranges should always be coded with the symmetric memory model in mind. !! The non-symmetric memory mode will then also work, albeit with a different (less efficient) communication pattern. From a2dfe21cbadbc999f9293705256495c3182a4fcb Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 12 Dec 2019 16:08:28 +0000 Subject: [PATCH 2/5] Adds global extents to hor_index type - niglobal,njglobal are available from the Domain type and usually accessed via a ocean_grid_type function. Modules that do not need the Domain but do use the hor_index type are obliged to use the ocean_grid_type to get the global extent, even though it's just integer information most closely related to the indices. - This commit adds niglobal,njglobal to the hor_index type so that modules that use HI, and only need indices, do not need to the full ocean_grid_type. --- src/framework/MOM_hor_index.F90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/framework/MOM_hor_index.F90 b/src/framework/MOM_hor_index.F90 index 8f7ce1db8e..db52afcdd8 100644 --- a/src/framework/MOM_hor_index.F90 +++ b/src/framework/MOM_hor_index.F90 @@ -3,7 +3,7 @@ module MOM_hor_index ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_domains, only : MOM_domain_type, get_domain_extent +use MOM_domains, only : MOM_domain_type, get_domain_extent, get_global_shape use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL use MOM_file_parser, only : get_param, log_param, log_version, param_file_type @@ -46,6 +46,9 @@ module MOM_hor_index integer :: idg_offset !< The offset between the corresponding global and local i-indices. integer :: jdg_offset !< The offset between the corresponding global and local j-indices. logical :: symmetric !< True if symmetric memory is used. + + integer :: niglobal !< The global number of h-cells in the i-direction + integer :: njglobal !< The global number of h-cells in the j-direction end type hor_index_type !> Copy the contents of one horizontal index type into another @@ -71,6 +74,7 @@ subroutine hor_index_init(Domain, HI, param_file, local_indexing, index_offset) HI%isg, HI%ieg, HI%jsg, HI%jeg, & HI%idg_offset, HI%jdg_offset, HI%symmetric, & local_indexing=local_indexing) + call get_global_shape(Domain, HI%niglobal, HI%njglobal) ! Read all relevant parameters and write them to the model log. call log_version(param_file, "MOM_hor_index", version, & @@ -108,6 +112,7 @@ subroutine HIT_assign(HI1, HI2) HI1%idg_offset = HI2%idg_offset ; HI1%jdg_offset = HI2%jdg_offset HI1%symmetric = HI2%symmetric + HI1%niglobal = HI2%niglobal ; HI1%njglobal = HI2%njglobal end subroutine HIT_assign From 0409fb239b034fa7baf35ca5b5bc6fc096097bde Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 20 Jan 2020 18:07:17 +0000 Subject: [PATCH 3/5] Adds MOM_random module - MOM_random provides a wrapper to the Mersenne twister pseudo-random number generator in FMS. In particular it allows for reproducible 2d fields of numbers using the MOM6 domain decomposition. - Commit includes unit tests that check specific values as well as the statistics (the latter are not unit tests but there's no harm in doing them). --- src/core/MOM_unit_tests.F90 | 3 + src/framework/MOM_random.F90 | 458 +++++++++++++++++++++++++++++++++++ 2 files changed, 461 insertions(+) create mode 100644 src/framework/MOM_random.F90 diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90 index ff5a93a62c..e01348f087 100644 --- a/src/core/MOM_unit_tests.F90 +++ b/src/core/MOM_unit_tests.F90 @@ -9,6 +9,7 @@ module MOM_unit_tests use MOM_remapping, only : remapping_unit_tests use MOM_neutral_diffusion, only : neutral_diffusion_unit_tests use MOM_diag_vkernels, only : diag_vkernels_unit_tests +use MOM_random, only : random_unit_tests implicit none ; private @@ -35,6 +36,8 @@ subroutine unit_tests(verbosity) "MOM_unit_tests: neutralDiffusionUnitTests FAILED") if (diag_vkernels_unit_tests(verbose)) call MOM_error(FATAL, & "MOM_unit_tests: diag_vkernels_unit_tests FAILED") + if (random_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: random_unit_tests FAILED") endif end subroutine unit_tests diff --git a/src/framework/MOM_random.F90 b/src/framework/MOM_random.F90 new file mode 100644 index 0000000000..6e254abed2 --- /dev/null +++ b/src/framework/MOM_random.F90 @@ -0,0 +1,458 @@ +!> Provides gridded random number capability +module MOM_random + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_hor_index, only : hor_index_type +use MOM_time_manager, only : time_type, set_date, get_date + +use MersenneTwister_mod, only : randomNumberSequence ! Random number class from FMS +use MersenneTwister_mod, only : new_RandomNumberSequence ! Constructor/initializer +use MersenneTwister_mod, only : getRandomReal ! Generates a random number +use MersenneTwister_mod, only : getRandomPositiveInt ! Generates a random positive integer + +implicit none ; private + +public :: random_0d_constructor +public :: random_01 +public :: random_norm +public :: random_2d_constructor +public :: random_2d_01 +public :: random_2d_norm +public :: random_unit_tests + +#include + +!> Container for pseudo-random number generators +type, public :: PRNG ; private + + !> Scalar random number generator for whole model + type(randomNumberSequence) :: stream0d + + !> Random number generator for each cell on horizontal grid + type(randomNumberSequence), dimension(:,:), allocatable :: stream2d + +end type PRNG + +contains + +!> Returns a random number between 0 and 1 +real function random_01(CS) + type(PRNG), intent(inout) :: CS !< Container for pseudo-random number generators + + random_01 = getRandomReal(CS%stream0d) + +end function random_01 + +!> Returns an approximately normally distributed random number with mean 0 and variance 1 +real function random_norm(CS) + type(PRNG), intent(inout) :: CS !< Container for pseudo-random number generators + ! Local variables + integer :: i + + random_norm = getRandomReal(CS%stream0d) - 0.5 + do i = 1,11 + random_norm = random_norm + ( getRandomReal(CS%stream0d) - 0.5 ) + enddo + +end function random_norm + +!> Generates random numbers between 0 and 1 for each cell of the model grid +subroutine random_2d_01(CS, HI, rand) + type(PRNG), intent(inout) :: CS !< Container for pseudo-random number generators + type(hor_index_type), intent(in) :: HI !< Horizontal index structure + real, dimension(SZI_(HI),SZJ_(HI)), intent(out) :: rand !< Random numbers between 0 and 1 + ! Local variables + integer :: i,j + + do j = HI%jsd,HI%jed + do i = HI%isd,HI%ied + rand(i,j) = getRandomReal( CS%stream2d(i,j) ) + enddo + enddo + +end subroutine random_2d_01 + +!> Returns an approximately normally distributed random number with mean 0 and variance 1 +!! for each cell of the model grid +subroutine random_2d_norm(CS, HI, rand) + type(PRNG), intent(inout) :: CS !< Container for pseudo-random number generators + type(hor_index_type), intent(in) :: HI !< Horizontal index structure + real, dimension(SZI_(HI),SZJ_(HI)), intent(out) :: rand !< Random numbers between 0 and 1 + ! Local variables + integer :: i,j,n + + do j = HI%jsd,HI%jed + do i = HI%isd,HI%ied + rand(i,j) = getRandomReal( CS%stream2d(i,j) ) - 0.5 + enddo + do n = 1,11 + do i = HI%isd,HI%ied + rand(i,j) = rand(i,j) + ( getRandomReal( CS%stream2d(i,j) ) - 0.5 ) + enddo + enddo + enddo + +end subroutine random_2d_norm + +!> Constructor for scalar PRNG. Can be used to reset the sequence. +subroutine random_0d_constructor(CS, Time, seed) + type(PRNG), intent(inout) :: CS !< Container for pseudo-random number generators + type(time_type), intent(in) :: Time !< Current model time + integer, intent(in) :: seed !< Seed for PRNG + ! Local variables + integer :: tseed + + tseed = seed_from_time(Time) + tseed = ieor(tseed, seed) + CS%stream0d = new_RandomNumberSequence(tseed) + +end subroutine random_0d_constructor + +!> Constructor for gridded PRNG. Can be used to reset the sequence. +subroutine random_2d_constructor(CS, HI, Time, seed) + type(PRNG), intent(inout) :: CS !< Container for pseudo-random number generators + type(hor_index_type), intent(in) :: HI !< Horizontal index structure + type(time_type), intent(in) :: Time !< Current model time + integer, intent(in) :: seed !< Seed for PRNG + ! Local variables + integer :: i,j,sseed,tseed + + if (.not. allocated(CS%stream2d)) allocate( CS%stream2d(HI%isd:HI%ied,HI%jsd:HI%jed) ) + + tseed = seed_from_time(Time) + tseed = ieor(tseed*9007, seed) + do j = HI%jsd,HI%jed + do i = HI%isd,HI%ied + sseed = seed_from_index(HI, i, j) + sseed = ieor(tseed, sseed*7993) + CS%stream2d(i,j) = new_RandomNumberSequence(sseed) + enddo + enddo + +end subroutine random_2d_constructor + +!> Return a seed derived as hash of values in Time +integer function seed_from_time(Time) + type(time_type), intent(in) :: Time !< Current model time + ! Local variables + integer :: yr,mo,dy,hr,mn,sc,s1,s2 + + call get_date(Time,yr,mo,dy,hr,mn,sc) + s1 = sc + 61*(mn + 61*hr) + 379 ! Range 379 .. 89620 + ! Fun fact: 2147483647 is the eighth Mersenne prime. + ! This is not the reason for using 2147483647+1 here. + s2 = mod(dy + 32*(mo + 13*yr), 2147483648) ! Range 0 .. 2147483647 + seed_from_time = ieor(s1*4111, s2) + +end function seed_from_time + +!> Create seed from position index +integer function seed_from_index(HI, i, j) + type(hor_index_type), intent(in) :: HI !< Horizontal index structure + integer, intent(in) :: i !< i-index (of h-cell) + integer, intent(in) :: j !< j-index (of h-cell) + ! Local variables + integer :: ig, jg, ni, nj, ij + + ni = HI%niglobal + nj = HI%njglobal + ! Periodicity is assumed here but does not break non-periodic models + ig = mod(HI%idg_offset + i - 1 + ni, ni)+1 + jg = max(HI%jdg_offset + j, 0) + if (jg>nj) then ! Tri-polar hard-coded until we put needed info in HI **TODO** + jg = 2*nj+1-jg + ig = ni+1-ig + endif + seed_from_index = ig + ni*(jg-1) + +end function seed_from_index + +!> Destructor for PRNG +subroutine random_destruct(CS) + type(PRNG), pointer :: CS !< Container for pseudo-random number generators + + if (allocated(CS%stream2d)) deallocate(CS%stream2d) + !deallocate(CS) +end subroutine random_destruct + +!> Runs some statistical tests on the PRNG +logical function random_unit_tests(verbose) + logical :: verbose !< True if results should be written to stdout + ! Local variables + type(PRNG) :: test_rng ! Generator + type(time_type) :: Time ! Model time + real :: r1, r2, r3 ! Some random numbers and re-used work variables + real :: mean, var, ar1, std ! Some statistics + integer :: stdunit ! For messages + integer, parameter :: n_samples = 800 + integer :: i, j, ni, nj + ! Fake being on a decomposed domain + type(hor_index_type), pointer :: HI => null() !< Not the real HI + real, dimension(:,:), allocatable :: r2d ! Random numbers + + ! Fake a decomposed domain + ni = 6 + nj = 9 + allocate(HI) + HI%isd = 0 + HI%ied = ni+1 + HI%jsd = 0 + HI%jed = nj+1 + HI%niglobal = ni + HI%njglobal = nj + HI%idg_offset = 0 + HI%jdg_offset = 0 + + random_unit_tests = .false. + stdunit = 6 + write(stdunit,'(1x,a)') '==== MOM_random: random_unit_tests =======================' + + if (verbose) write(stdunit,'(1x,"random: ",a)') '-- Time-based seeds ---------------------' + ! Check time-based seed generation + Time = set_date(1903, 11, 21, 13, 47, 29) + i = seed_from_time(Time) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, i==212584341, 'time seed 1903/11/21 13:47:29', ivalue=i) + Time = set_date(1903, 11, 22, 13, 47, 29) + i = seed_from_time(Time) + random_unit_tests = random_unit_tests .or.& + test_fn(verbose, i==212584342, 'time seed 1903/11/22 13:47:29', ivalue=i) + Time = set_date(1903, 11, 21, 13, 47, 30) + i = seed_from_time(Time) + random_unit_tests = random_unit_tests .or.& + test_fn(verbose, i==212596634, 'time seed 1903/11/21 13:47:30', ivalue=i) + + if (verbose) write(stdunit,'(1x,"random: ",a)') '-- PRNG tests ---------------------------' + ! Generate a random number, r1 + call random_0d_constructor(test_rng, Time, 1) + r1 = random_01(test_rng) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r1-4.75310122e-2)<1.e-9, 'first call', r1) + + ! Check that we get a different number, r2, on a second call + r2 = random_01(test_rng) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r2-2.71289742e-1)<1.e-9, 'consecutive test', r2) + + ! Check that we can reproduce r1 by resetting the seed + call random_0d_constructor(test_rng, Time, 1) + r2 = random_01(test_rng) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r2-r1)==0., 'reproduce test', r2) + + ! Check that we get a different number, r2, with a different seed but same date + call random_0d_constructor(test_rng, Time, 2) + r2 = random_01(test_rng) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r2-7.15508473e-1)<1.e-9, 'different seed test', r2) + + ! Check that we get a different number, r2, for a different date but same seed + Time = set_date(1903, 11, 21, 13, 0, 29) + call random_0d_constructor(test_rng, Time, 1) + r2 = random_01(test_rng) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r2-9.56667163e-1)<1.e-9, 'different date test', r2) + + if (verbose) write(stdunit,'(1x,"random: ",a)') '-- index-based seeds --------------------' + ! Check index-based seed + i = seed_from_index(HI,1,1) + random_unit_tests = random_unit_tests .or. test_fn(verbose, i==1, 'seed from index (1,1)', ivalue=i) + j = seed_from_index(HI,ni+1,1) + random_unit_tests = random_unit_tests .or. test_fn(verbose, j==i, 'seed from index (n+1,1)', ivalue=j) + i = seed_from_index(HI,ni,1) + random_unit_tests = random_unit_tests .or. test_fn(verbose, i==6, 'seed from index (n,1)', ivalue=i) + j = seed_from_index(HI,0,1) + random_unit_tests = random_unit_tests .or. test_fn(verbose, j==i, 'seed from index (0,1)', ivalue=j) + i = seed_from_index(HI,1,nj) + random_unit_tests = random_unit_tests .or. test_fn(verbose, i==49, 'seed from index (1,n)', ivalue=i) + j = seed_from_index(HI,ni,nj+1) + random_unit_tests = random_unit_tests .or. test_fn(verbose, j==i, 'seed from index (n,n+1)', ivalue=j) + i = seed_from_index(HI,ni,nj) + random_unit_tests = random_unit_tests .or. test_fn(verbose, i==54, 'seed from index (n,n)', ivalue=i) + j = seed_from_index(HI,1,nj+1) + random_unit_tests = random_unit_tests .or. test_fn(verbose, j==i, 'seed from index (1,n+1)', ivalue=j) + + if (.not.random_unit_tests) write(stdunit,'(1x,a)') 'Passed unit tests' + ! The rest of these are not unit tests but statistical tests and as such + ! could fail for different sample sizes but happen to pass here. + + ! Check statistics of large samples for uniform generator + mean = 0. ; var = 0. ; ar1 = 0. ; r2 = 0. + do i = 1, n_samples + r1 = random_01(test_rng) - 0.5 + mean = mean + r1 + var = var + r1**2 + ar1 = ar1 + r1*r2 + r2 = r1 ! Keep copy of last value + enddo + mean = mean / real(n_samples) ! Expected mean is 0 + var = var / real(n_samples) ! Expected variance is 1/12 + ar1 = ar1 / real(n_samples-1) ! Autocovariance + std = sqrt(var) ! Expected std is sqrt(1/12) + r2 = mean*sqrt(real(12*n_samples)) ! Normalized error in mean + r3 = std*sqrt(12.) ! Normalized standard deviation + r1 = ( ar1 * sqrt(real(n_samples-1)) ) / var + if (verbose) then + write(stdunit,'(1x,"random: ",a)') '-- Uniform -0.5 .. 0.5 generator --------' + write(stdunit,'(1x,"random: ",a,f12.9)') 'mean =',mean,'std =',std,'AR1 =',ar1 + write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. mean =',r2, & + 'norm. std =',r3,'norm. AR1 =',r1 + endif + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r2)<2., & + 'n>>1, mean within 2 sigma [uniform]', r2) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r3-1.)<1./sqrt(real(n_samples)), & + 'n>>1, std ~ 1/sqrt(12) [uniform]', r3-1.) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r1)<2., & + 'n>>1, AR1 < std/sqrt(n) [uniform]', r1) + + ! Check statistics of large samples for normal generator + mean = 0. ; var = 0. ; ar1 = 0. ; r2 = 0. + do i = 1, n_samples + r1 = random_norm(test_rng) + mean = mean + r1 + var = var + r1**2 + ar1 = ar1 + r1*r2 + r2 = r1 ! Keep copy of last value for AR calculation + enddo + mean = mean / real(n_samples) + var = var / real(n_samples) + ar1 = ar1 / real(n_samples) + std = sqrt(var) + r3 = 1./sqrt(real(n_samples)) ! Standard error of mean + r2 = mean*sqrt(real(n_samples)) ! Normalized error in mean + r3 = std ! Normalized standard deviation + r1 = ( ar1 * sqrt(real(n_samples-1)) ) / var + if (verbose) then + write(stdunit,'(1x,"random: ",a)') '-- Normal distribution generator --------' + write(stdunit,'(1x,"random: ",a,f12.9)') 'mean =',mean,'std =',std,'AR1 =',ar1 + write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. error in mean =',r2, & + 'norm. standard deviation =',r3,'norm. AR1 =',r1 + endif + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r2)<2., & + 'n>>1, mean within 2 sigma [norm]', r2) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r3-1.)<1./sqrt(real(n_samples)), & + 'n>>1, std ~ 1 [norm]', r3-1.) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r1)<2., & + 'n>>1, AR1 < std/sqrt(n) [norm]', r1) + + if (verbose) write(stdunit,'(1x,"random: ",a)') '-- 2d PRNG ------------------------------' + ! Check 2d random number generator 0..1 + allocate( r2d(HI%isd:HI%ied,HI%jsd:HI%jed) ) + call random_2d_constructor(test_rng, HI, Time, 123) + r2d(:,:) = -999. ! Use -9. to detect unset values + call random_2d_01(test_rng, HI, r2d) + if (any(abs(r2d(:,:)+999.)<=0.)) random_unit_tests=.true. + r1 = minval(r2d) + r2 = maxval(r2d) + random_unit_tests = random_unit_tests .or. test_fn(verbose, r1>=0., '2d all set', r1) + random_unit_tests = random_unit_tests .or. test_fn(verbose, r2<=1., '2d all valid', r2) + mean = sum( r2d(1:ni,1:nj) - 0.5 )/real(ni*nj) + var = sum( (r2d(1:ni,1:nj) - 0.5 - mean)**2 )/real(ni*nj) + std = sqrt(var) + r3 = 1./sqrt(real(12*ni*nj)) ! Standard error of mean + r2 = mean*sqrt(real(12*ni*nj)) ! Normalized error in mean + r3 = std*sqrt(12.) ! Normalized standard deviation + if (verbose) then + write(stdunit,'(1x,"random: ",a)') '2D uniform 0..1 generator' + write(stdunit,'(1x,"random: ",a,f12.9)') 'mean =',mean,'std =',std + write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. error in mean =',r2 + write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. standard deviation =',r3 + endif + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r2)<2., & + '2d, mean within 2 sigma [uniform]', r2) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r3-1.)<1./sqrt(real(ni*nj)), & + '2d, std ~ 1/sqrt(12) [uniform]', r3-1.) + if (verbose) then + write(stdunit,'(1x,"random:")') + write(stdunit,'(1x,"random:",8f8.5)') r2d + write(stdunit,'(1x,"random:")') + endif + + ! Check 2d normal random number generator + call random_2d_norm(test_rng, HI, r2d) + mean = sum( r2d(1:ni,1:nj) )/real(ni*nj) + var = sum( r2d(1:ni,1:nj)**2 )/real(ni*nj) + std = sqrt(var) + r3 = 1./sqrt(real(ni*nj)) ! Standard error of mean + r2 = mean*sqrt(real(ni*nj)) ! Normalized error in mean + r3 = std ! Normalized standard deviation + if (verbose) then + write(stdunit,'(1x,"random: ",a)') '2D normal generator' + write(stdunit,'(1x,"random: ",a,f12.9)') 'mean =',mean,'std =',std + write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. error in mean =',r2 + write(stdunit,'(1x,"random: ",a,f12.9)') 'norm. standard deviation =',r3 + endif + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r2)<2., & + '2d, mean within 2 sigma [norm]', r2) + random_unit_tests = random_unit_tests .or. & + test_fn(verbose, abs(r3-1.)<1./sqrt(real(ni*nj)), & + '2d, std ~ 1/sqrt(12) [norm]', r3-1.) + + ! Clean up + deallocate(r2d) + deallocate(HI) + + if (.not.random_unit_tests) write(stdunit,'(1x,a)') 'Passed statistical tests' + +end function random_unit_tests + +!> Convenience function for reporting result of test +logical function test_fn(verbose, good, label, rvalue, ivalue) + logical, intent(in) :: verbose !< Verbosity + logical, intent(in) :: good !< True if pass, false otherwise + character(len=*), intent(in) :: label !< Label for messages + real, intent(in) :: rvalue !< Result of calculation + integer, intent(in) :: ivalue !< Result of calculation + optional :: rvalue, ivalue + + if (present(ivalue)) then + if (.not. good) then + write(0,'(1x,a,i10,1x,a,a)') 'random: result =',ivalue,label,' <------- FAIL!' + elseif (verbose) then + write(6,'(1x,a,i10,1x,a)') 'random: result =',ivalue,label + endif + else + if (.not. good) then + write(0,'(1x,a,1pe15.8,1x,a,a)') 'random: result =',rvalue,label,' <------- FAIL!' + elseif (verbose) then + write(6,'(1x,a,1pe15.8,1x,a)') 'random: result =',rvalue,label + endif + endif + test_fn = .not. good + +end function test_fn + +end module MOM_random + +!> \namespace mom_random +!! +!! Provides MOM6 wrappers to the FMS implementation of the Mersenne twister. +!! +!! Example usage: +!! \code +!! type(PRNG) :: rng +!! real :: rn +!! call random_0d_constructor(rng, Time, seed) ! Call this each time-step +!! rn = random_01(rng) +!! rn = random_norm(rng) +!! +!! type(PRNG) :: rng +!! real, dimension(:,:) :: rn2d +!! call random_2d_constructor(rng, HI, Time, seed) ! Call this each time-step +!! call random_2d_01(rng, HI, rn2d) +!! call random_2d_norm(rng, HI, rn2d) +!! +!! Note: reproducibility across restarts is implemented by using time-derived +!! seeds to pass to the Mersenne twister. It is therefore important that any +!! PRNG type be re-initialized each time-step. +!! \endcode From 4ce3ffc87e267bd86d51f58eb1ec8025f476a261 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Sat, 25 Jan 2020 18:37:39 +0000 Subject: [PATCH 4/5] Replaced intent(out) with intent(inout) in EOS - For arrays using intent(inout) avoids potential copies and also allows aliasing of arguments which can reduce the overhead of multiple spare arrays in the calling code. --- src/equation_of_state/MOM_EOS.F90 | 201 +++++++++++------------ src/equation_of_state/MOM_EOS_Wright.F90 | 110 ++++++------- 2 files changed, 155 insertions(+), 156 deletions(-) diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 5d3d33534b..3d0cb9abc4 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -136,14 +136,14 @@ module MOM_EOS !> Calls the appropriate subroutine to calculate density of sea water for scalar inputs. !! If rho_ref is present, the anomaly with respect to rho_ref is returned. subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, scale) - real, intent(in) :: T !< Potential temperature referenced to the surface [degC] - real, intent(in) :: S !< Salinity [ppt] + real, intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, intent(in) :: S !< Salinity [ppt] real, intent(in) :: pressure !< Pressure [Pa] - real, intent(out) :: rho !< Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3] - type(EOS_type), pointer :: EOS !< Equation of state structure + real, intent(out) :: rho !< Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3] + type(EOS_type), pointer :: EOS !< Equation of state structure real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. - real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density - !! from kg m-3 to the desired units [R m3 kg-1] + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density + !! from kg m-3 to the desired units [R m3 kg-1] if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_scalar called with an unassociated EOS_type EOS.") @@ -174,16 +174,16 @@ end subroutine calculate_density_scalar !> Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs. !! If rho_ref is present, the anomaly with respect to rho_ref is returned. subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_ref, scale) - real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] - real, dimension(:), intent(in) :: S !< Salinity [ppt] - real, dimension(:), intent(in) :: pressure !< Pressure [Pa] - real, dimension(:), intent(out) :: rho !< Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3] - integer, intent(in) :: start !< Start index for computation - integer, intent(in) :: npts !< Number of point to compute - type(EOS_type), pointer :: EOS !< Equation of state structure - real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. - real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density - !! from kg m-3 to the desired units [R m3 kg-1] + real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(:), intent(in) :: S !< Salinity [ppt] + real, dimension(:), intent(in) :: pressure !< Pressure [Pa] + real, dimension(:), intent(inout) :: rho !< Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3] + integer, intent(in) :: start !< Start index for computation + integer, intent(in) :: npts !< Number of point to compute + type(EOS_type), pointer :: EOS !< Equation of state structure + real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density + !! from kg m-3 to the desired units [R m3 kg-1] integer :: j if (.not.associated(EOS)) call MOM_error(FATAL, & @@ -215,14 +215,14 @@ end subroutine calculate_density_array !> Calls the appropriate subroutine to calculate specific volume of sea water !! for scalar inputs. subroutine calculate_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, scale) - real, intent(in) :: T !< Potential temperature referenced to the surface [degC] - real, intent(in) :: S !< Salinity [ppt] + real, intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, intent(in) :: S !< Salinity [ppt] real, intent(in) :: pressure !< Pressure [Pa] real, intent(out) :: specvol !< In situ? specific volume [m3 kg-1] or [R-1 ~> m3 kg-1] type(EOS_type), pointer :: EOS !< Equation of state structure real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. - real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific volume - !! from m3 kg-1 to the desired units [kg m-3 R-1] + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific + !! volume from m3 kg-1 to the desired units [kg m-3 R-1] real :: rho @@ -261,17 +261,16 @@ end subroutine calculate_spec_vol_scalar !> Calls the appropriate subroutine to calculate the specific volume of sea water !! for 1-D array inputs. subroutine calculate_spec_vol_array(T, S, pressure, specvol, start, npts, EOS, spv_ref, scale) - real, dimension(:), intent(in) :: T !< potential temperature relative to the surface - !! [degC]. - real, dimension(:), intent(in) :: S !< salinity [ppt]. - real, dimension(:), intent(in) :: pressure !< pressure [Pa]. - real, dimension(:), intent(out) :: specvol !< in situ specific volume [kg m-3] or [R-1 ~> m3 kg-1]. - integer, intent(in) :: start !< the starting point in the arrays. - integer, intent(in) :: npts !< the number of values to calculate. - type(EOS_type), pointer :: EOS !< Equation of state structure - real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. - real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific volume - !! from m3 kg-1 to the desired units [kg m-3 R-1] + real, dimension(:), intent(in) :: T !< potential temperature relative to the surface [degC]. + real, dimension(:), intent(in) :: S !< salinity [ppt]. + real, dimension(:), intent(in) :: pressure !< pressure [Pa]. + real, dimension(:), intent(inout) :: specvol !< in situ specific volume [kg m-3] or [R-1 ~> m3 kg-1]. + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + type(EOS_type), pointer :: EOS !< Equation of state structure + real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific volume + !! from m3 kg-1 to the desired units [kg m-3 R-1] real, dimension(size(specvol)) :: rho integer :: j @@ -336,13 +335,13 @@ end subroutine calculate_TFreeze_scalar !> Calls the appropriate subroutine to calculate the freezing point for a 1-D array. subroutine calculate_TFreeze_array(S, pressure, T_fr, start, npts, EOS) - real, dimension(:), intent(in) :: S !< Salinity [ppt] - real, dimension(:), intent(in) :: pressure !< Pressure [Pa] - real, dimension(:), intent(out) :: T_fr !< Freezing point potential temperature referenced - !! to the surface [degC] - integer, intent(in) :: start !< Starting index within the array - integer, intent(in) :: npts !< The number of values to calculate - type(EOS_type), pointer :: EOS !< Equation of state structure + real, dimension(:), intent(in) :: S !< Salinity [ppt] + real, dimension(:), intent(in) :: pressure !< Pressure [Pa] + real, dimension(:), intent(inout) :: T_fr !< Freezing point potential temperature referenced + !! to the surface [degC] + integer, intent(in) :: start !< Starting index within the array + integer, intent(in) :: npts !< The number of values to calculate + type(EOS_type), pointer :: EOS !< Equation of state structure if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_TFreeze_scalar called with an unassociated EOS_type EOS.") @@ -364,18 +363,18 @@ end subroutine calculate_TFreeze_array !> Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs. subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, start, npts, EOS, scale) - real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] - real, dimension(:), intent(in) :: S !< Salinity [ppt] - real, dimension(:), intent(in) :: pressure !< Pressure [Pa] - real, dimension(:), intent(out) :: drho_dT !< The partial derivative of density with potential - !! temperature [kg m-3 degC-1] or [R degC-1 ~> kg m-3 degC-1]. - real, dimension(:), intent(out) :: drho_dS !< The partial derivative of density with salinity, - !! in [kg m-3 ppt-1] or [R degC-1 ~> kg m-3 ppt-1]. - integer, intent(in) :: start !< Starting index within the array - integer, intent(in) :: npts !< The number of values to calculate - type(EOS_type), pointer :: EOS !< Equation of state structure - real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density - !! from kg m-3 to the desired units [R m3 kg-1] + real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(:), intent(in) :: S !< Salinity [ppt] + real, dimension(:), intent(in) :: pressure !< Pressure [Pa] + real, dimension(:), intent(inout) :: drho_dT !< The partial derivative of density with potential + !! temperature [kg m-3 degC-1] or [R degC-1 ~> kg m-3 degC-1]. + real, dimension(:), intent(inout) :: drho_dS !< The partial derivative of density with salinity, + !! in [kg m-3 ppt-1] or [R degC-1 ~> kg m-3 ppt-1]. + integer, intent(in) :: start !< Starting index within the array + integer, intent(in) :: npts !< The number of values to calculate + type(EOS_type), pointer :: EOS !< Equation of state structure + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density + !! from kg m-3 to the desired units [R m3 kg-1] integer :: j if (.not.associated(EOS)) call MOM_error(FATAL, & @@ -448,16 +447,16 @@ subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drh real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] real, dimension(:), intent(in) :: S !< Salinity [ppt] real, dimension(:), intent(in) :: pressure !< Pressure [Pa] - real, dimension(:), intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S - !! [kg m-3 ppt-2] or [R ppt-2 ~> kg m-3 ppt-2] - real, dimension(:), intent(out) :: drho_dS_dT !< Partial derivative of beta with respect to T - !! [kg m-3 ppt-1 degC-1] or [R ppt-1 degC-1 ~> kg m-3 ppt-1 degC-1] - real, dimension(:), intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T - !! [kg m-3 degC-2] or [R degC-2 ~> kg m-3 degC-2] - real, dimension(:), intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure - !! [kg m-3 ppt-1 Pa-1] or [R ppt-1 Pa-1 ~> kg m-3 ppt-1 Pa-1] - real, dimension(:), intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure - !! [kg m-3 degC-1 Pa-1] or [R degC-1 Pa-1 ~> kg m-3 degC-1 Pa-1] + real, dimension(:), intent(inout) :: drho_dS_dS !< Partial derivative of beta with respect to S + !! [kg m-3 ppt-2] or [R ppt-2 ~> kg m-3 ppt-2] + real, dimension(:), intent(inout) :: drho_dS_dT !< Partial derivative of beta with respect to T + !! [kg m-3 ppt-1 degC-1] or [R ppt-1 degC-1 ~> kg m-3 ppt-1 degC-1] + real, dimension(:), intent(inout) :: drho_dT_dT !< Partial derivative of alpha with respect to T + !! [kg m-3 degC-2] or [R degC-2 ~> kg m-3 degC-2] + real, dimension(:), intent(inout) :: drho_dS_dP !< Partial derivative of beta with respect to pressure + !! [kg m-3 ppt-1 Pa-1] or [R ppt-1 Pa-1 ~> kg m-3 ppt-1 Pa-1] + real, dimension(:), intent(inout) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure + !! [kg m-3 degC-1 Pa-1] or [R degC-1 Pa-1 ~> kg m-3 degC-1 Pa-1] integer, intent(in) :: start !< Starting index within the array integer, intent(in) :: npts !< The number of values to calculate type(EOS_type), pointer :: EOS !< Equation of state structure @@ -546,10 +545,10 @@ subroutine calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] real, dimension(:), intent(in) :: S !< Salinity [ppt] real, dimension(:), intent(in) :: pressure !< Pressure [Pa] - real, dimension(:), intent(out) :: dSV_dT !< The partial derivative of specific volume with potential - !! temperature [m3 kg-1 degC-1] or [R-1 degC-1 ~> m3 kg-1 degC-1] - real, dimension(:), intent(out) :: dSV_dS !< The partial derivative of specific volume with salinity - !! [m3 kg-1 ppt-1] or [R-1 ppt-1 ~> m3 kg-1 ppt-1]. + real, dimension(:), intent(inout) :: dSV_dT !< The partial derivative of specific volume with potential + !! temperature [m3 kg-1 degC-1] or [R-1 degC-1 ~> m3 kg-1 degC-1] + real, dimension(:), intent(inout) :: dSV_dS !< The partial derivative of specific volume with salinity + !! [m3 kg-1 ppt-1] or [R-1 ppt-1 ~> m3 kg-1 ppt-1]. integer, intent(in) :: start !< Starting index within the array integer, intent(in) :: npts !< The number of values to calculate type(EOS_type), pointer :: EOS !< Equation of state structure @@ -603,9 +602,9 @@ subroutine calculate_compress_array(T, S, pressure, rho, drho_dp, start, npts, E real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] real, dimension(:), intent(in) :: S !< Salinity [PSU] real, dimension(:), intent(in) :: pressure !< Pressure [Pa] - real, dimension(:), intent(out) :: rho !< In situ density [kg m-3]. - real, dimension(:), intent(out) :: drho_dp !< The partial derivative of density with pressure - !! (also the inverse of the square of sound speed) [s2 m-2]. + real, dimension(:), intent(inout) :: rho !< In situ density [kg m-3]. + real, dimension(:), intent(inout) :: drho_dp !< The partial derivative of density with pressure + !! (also the inverse of the square of sound speed) [s2 m-2]. integer, intent(in) :: start !< Starting index within the array integer, intent(in) :: npts !< The number of values to calculate type(EOS_type), pointer :: EOS !< Equation of state structure @@ -677,18 +676,18 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & !! alpha_ref, but this reduces the effects of roundoff. type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(out) :: dza !< The change in the geopotential anomaly across + intent(inout) :: dza !< The change in the geopotential anomaly across !! the layer [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(out) :: intp_dza !< The integral in pressure through the layer of the + optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of the !! geopotential anomaly relative to the anomaly at the bottom of the !! layer [Pa m2 s-2]. real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & - optional, intent(out) :: intx_dza !< The integral in x of the difference between the + optional, intent(inout) :: intx_dza !< The integral in x of the difference between the !! geopotential anomaly at the top and bottom of the layer divided by !! the x grid spacing [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & - optional, intent(out) :: inty_dza !< The integral in y of the difference between the + optional, intent(inout) :: inty_dza !< The integral in y of the difference between the !! geopotential anomaly at the top and bottom of the layer divided by !! the y grid spacing [m2 s-2]. integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza. @@ -747,17 +746,17 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, & real, intent(in) :: G_e !< The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2]. type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa]. + intent(inout) :: dpa !< The change in the pressure anomaly across the layer [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of + optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of !! the pressure anomaly relative to the anomaly at the !! top of the layer [Pa Z ~> Pa m]. real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), & - optional, intent(out) :: intx_dpa !< The integral in x of the difference between the + optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the !! pressure anomaly at the top and bottom of the layer !! divided by the x grid spacing [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & - optional, intent(out) :: inty_dpa !< The integral in y of the difference between the + optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the !! pressure anomaly at the top and bottom of the layer !! divided by the y grid spacing [Pa]. real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & @@ -1000,18 +999,18 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, real, intent(in) :: G_e !< The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2]. type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - intent(out) :: dpa !< The change in the pressure anomaly + intent(inout) :: dpa !< The change in the pressure anomaly !! across the layer [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - optional, intent(out) :: intz_dpa !< The integral through the thickness of the + optional, intent(inout) :: intz_dpa !< The integral through the thickness of the !! layer of the pressure anomaly relative to the !! anomaly at the top of the layer [Pa Z ~> Pa m]. real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), & - optional, intent(out) :: intx_dpa !< The integral in x of the difference between + optional, intent(inout) :: intx_dpa !< The integral in x of the difference between !! the pressure anomaly at the top and bottom of the !! layer divided by the x grid spacing [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & - optional, intent(out) :: inty_dpa !< The integral in y of the difference between + optional, intent(inout) :: inty_dpa !< The integral in y of the difference between !! the pressure anomaly at the top and bottom of the !! layer divided by the y grid spacing [Pa]. real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & @@ -1193,17 +1192,17 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]. type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa]. + intent(inout) :: dpa !< The change in the pressure anomaly across the layer [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of + optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of !! the pressure anomaly relative to the anomaly at the !! top of the layer [Pa Z]. real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), & - optional, intent(out) :: intx_dpa !< The integral in x of the difference between the + optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the !! pressure anomaly at the top and bottom of the layer !! divided by the x grid spacing [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & - optional, intent(out) :: inty_dpa !< The integral in y of the difference between the + optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the !! pressure anomaly at the top and bottom of the layer !! divided by the y grid spacing [Pa]. logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to @@ -1608,17 +1607,17 @@ subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, & real, intent(in) :: G_e !< The Earth's gravitational acceleration [m s-2] type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa]. + intent(inout) :: dpa !< The change in the pressure anomaly across the layer [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of + optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of !! the pressure anomaly relative to the anomaly at the !! top of the layer [Pa Z ~> Pa m]. real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), & - optional, intent(out) :: intx_dpa !< The integral in x of the difference between the + optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the !! pressure anomaly at the top and bottom of the layer !! divided by the x grid spacing [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & - optional, intent(out) :: inty_dpa !< The integral in y of the difference between the + optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the !! pressure anomaly at the top and bottom of the layer !! divided by the y grid spacing [Pa]. @@ -1916,10 +1915,10 @@ end subroutine compute_integral_quadratic subroutine evaluate_shape_bilinear( xi, eta, phi, dphidxi, dphideta ) real, intent(in) :: xi !< The x position to evaluate real, intent(in) :: eta !< The z position to evaluate - real, dimension(4), intent(out) :: phi !< The weights of the four corners at this point - real, dimension(4), intent(out) :: dphidxi !< The x-gradient of the weights of the four + real, dimension(4), intent(inout) :: phi !< The weights of the four corners at this point + real, dimension(4), intent(inout) :: dphidxi !< The x-gradient of the weights of the four !! corners at this point - real, dimension(4), intent(out) :: dphideta !< The z-gradient of the weights of the four + real, dimension(4), intent(inout) :: dphideta !< The z-gradient of the weights of the four !! corners at this point ! The shape functions within the parent element are defined as shown here: @@ -1957,11 +1956,11 @@ subroutine evaluate_shape_quadratic ( xi, eta, phi, dphidxi, dphideta ) ! Arguments real, intent(in) :: xi !< The x position to evaluate real, intent(in) :: eta !< The z position to evaluate - real, dimension(9), intent(out) :: phi !< The weights of the 9 bilinear quadrature points + real, dimension(9), intent(inout) :: phi !< The weights of the 9 bilinear quadrature points !! at this point - real, dimension(9), intent(out) :: dphidxi !< The x-gradient of the weights of the 9 bilinear + real, dimension(9), intent(inout) :: dphidxi !< The x-gradient of the weights of the 9 bilinear !! quadrature points corners at this point - real, dimension(9), intent(out) :: dphideta !< The z-gradient of the weights of the 9 bilinear + real, dimension(9), intent(inout) :: dphideta !< The z-gradient of the weights of the 9 bilinear !! quadrature points corners at this point ! The quadratic shape functions within the parent element are defined as shown here: @@ -2037,18 +2036,18 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & !! alters the effects of roundoff, and answers do change. type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(out) :: dza !< The change in the geopotential anomaly + intent(inout) :: dza !< The change in the geopotential anomaly !! across the layer [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(out) :: intp_dza !< The integral in pressure through the + optional, intent(inout) :: intp_dza !< The integral in pressure through the !! layer of the geopotential anomaly relative to the anomaly !! at the bottom of the layer [Pa m2 s-2]. real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & - optional, intent(out) :: intx_dza !< The integral in x of the difference + optional, intent(inout) :: intx_dza !< The integral in x of the difference !! between the geopotential anomaly at the top and bottom of !! the layer divided by the x grid spacing [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & - optional, intent(out) :: inty_dza !< The integral in y of the difference + optional, intent(inout) :: inty_dza !< The integral in y of the difference !! between the geopotential anomaly at the top and bottom of !! the layer divided by the y grid spacing [m2 s-2]. integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza. @@ -2235,18 +2234,18 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, intent(in) :: bathyP !< The pressure at the bathymetry [Pa] type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(out) :: dza !< The change in the geopotential anomaly + intent(inout) :: dza !< The change in the geopotential anomaly !! across the layer [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(out) :: intp_dza !< The integral in pressure through the + optional, intent(inout) :: intp_dza !< The integral in pressure through the !! layer of the geopotential anomaly relative to the anomaly !! at the bottom of the layer [Pa m2 s-2]. real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & - optional, intent(out) :: intx_dza !< The integral in x of the difference + optional, intent(inout) :: intx_dza !< The integral in x of the difference !! between the geopotential anomaly at the top and bottom of !! the layer divided by the x grid spacing [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & - optional, intent(out) :: inty_dza !< The integral in y of the difference + optional, intent(inout) :: inty_dza !< The integral in y of the difference !! between the geopotential anomaly at the top and bottom of !! the layer divided by the y grid spacing [m2 s-2]. logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index 899f32b27d..8d29e08f92 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -108,13 +108,13 @@ end subroutine calculate_density_scalar_wright !! (T [degC]), and pressure [Pa]. It uses the expression from !! Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. subroutine calculate_density_array_wright(T, S, pressure, rho, start, npts, rho_ref) - real, dimension(:), intent(in) :: T !< potential temperature relative to the surface [degC]. - real, dimension(:), intent(in) :: S !< salinity [PSU]. - real, dimension(:), intent(in) :: pressure !< pressure [Pa]. - real, dimension(:), intent(out) :: rho !< in situ density [kg m-3]. - integer, intent(in) :: start !< the starting point in the arrays. - integer, intent(in) :: npts !< the number of values to calculate. - real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. + real, dimension(:), intent(in) :: T !< potential temperature relative to the surface [degC]. + real, dimension(:), intent(in) :: S !< salinity [PSU]. + real, dimension(:), intent(in) :: pressure !< pressure [Pa]. + real, dimension(:), intent(inout) :: rho !< in situ density [kg m-3]. + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. ! Original coded by R. Hallberg, 7/00, anomaly coded in 3/18. ! Local variables @@ -169,14 +169,14 @@ end subroutine calculate_spec_vol_scalar_wright !! Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. !! If spv_ref is present, specvol is an anomaly from spv_ref. subroutine calculate_spec_vol_array_wright(T, S, pressure, specvol, start, npts, spv_ref) - real, dimension(:), intent(in) :: T !< potential temperature relative to the + real, dimension(:), intent(in) :: T !< potential temperature relative to the !! surface [degC]. - real, dimension(:), intent(in) :: S !< salinity [PSU]. - real, dimension(:), intent(in) :: pressure !< pressure [Pa]. - real, dimension(:), intent(out) :: specvol !< in situ specific volume [m3 kg-1]. - integer, intent(in) :: start !< the starting point in the arrays. - integer, intent(in) :: npts !< the number of values to calculate. - real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. + real, dimension(:), intent(in) :: S !< salinity [PSU]. + real, dimension(:), intent(in) :: pressure !< pressure [Pa]. + real, dimension(:), intent(inout) :: specvol !< in situ specific volume [m3 kg-1]. + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. ! Local variables real :: al0, p0, lambda @@ -197,16 +197,16 @@ end subroutine calculate_spec_vol_array_wright !> For a given thermodynamic state, return the thermal/haline expansion coefficients subroutine calculate_density_derivs_array_wright(T, S, pressure, drho_dT, drho_dS, start, npts) - real, intent(in), dimension(:) :: T !< Potential temperature relative to the - !! surface [degC]. - real, intent(in), dimension(:) :: S !< Salinity [PSU]. - real, intent(in), dimension(:) :: pressure !< pressure [Pa]. - real, intent(out), dimension(:) :: drho_dT !< The partial derivative of density with potential - !! temperature [kg m-3 degC-1]. - real, intent(out), dimension(:) :: drho_dS !< The partial derivative of density with salinity, - !! in [kg m-3 PSU-1]. - integer, intent(in) :: start !< The starting point in the arrays. - integer, intent(in) :: npts !< The number of values to calculate. + real, intent(in), dimension(:) :: T !< Potential temperature relative to the + !! surface [degC]. + real, intent(in), dimension(:) :: S !< Salinity [PSU]. + real, intent(in), dimension(:) :: pressure !< pressure [Pa]. + real, intent(inout), dimension(:) :: drho_dT !< The partial derivative of density with potential + !! temperature [kg m-3 degC-1]. + real, intent(inout), dimension(:) :: drho_dS !< The partial derivative of density with salinity, + !! in [kg m-3 PSU-1]. + integer, intent(in) :: start !< The starting point in the arrays. + integer, intent(in) :: npts !< The number of values to calculate. ! Local variables real :: al0, p0, lambda, I_denom2 @@ -259,15 +259,15 @@ subroutine calculate_density_second_derivs_array_wright(T, S, P, drho_ds_ds, drh real, dimension(:), intent(in ) :: T !< Potential temperature referenced to 0 dbar [degC] real, dimension(:), intent(in ) :: S !< Salinity [PSU] real, dimension(:), intent(in ) :: P !< Pressure [Pa] - real, dimension(:), intent( out) :: drho_ds_ds !< Partial derivative of beta with respect + real, dimension(:), intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect !! to S [kg m-3 PSU-2] - real, dimension(:), intent( out) :: drho_ds_dt !< Partial derivative of beta with respcct + real, dimension(:), intent(inout) :: drho_ds_dt !< Partial derivative of beta with respcct !! to T [kg m-3 PSU-1 degC-1] - real, dimension(:), intent( out) :: drho_dt_dt !< Partial derivative of alpha with respect + real, dimension(:), intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect !! to T [kg m-3 degC-2] - real, dimension(:), intent( out) :: drho_ds_dp !< Partial derivative of beta with respect + real, dimension(:), intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect !! to pressure [kg m-3 PSU-1 Pa-1] - real, dimension(:), intent( out) :: drho_dt_dp !< Partial derivative of alpha with respect + real, dimension(:), intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect !! to pressure [kg m-3 degC-1 Pa-1] integer, intent(in ) :: start !< Starting index in T,S,P integer, intent(in ) :: npts !< Number of points to loop over @@ -340,15 +340,15 @@ end subroutine calculate_density_second_derivs_scalar_wright !> For a given thermodynamic state, return the partial derivatives of specific volume !! with temperature and salinity subroutine calculate_specvol_derivs_wright(T, S, pressure, dSV_dT, dSV_dS, start, npts) - real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface [degC]. - real, intent(in), dimension(:) :: S !< Salinity [PSU]. - real, intent(in), dimension(:) :: pressure !< pressure [Pa]. - real, intent(out), dimension(:) :: dSV_dT !< The partial derivative of specific volume with - !! potential temperature [m3 kg-1 degC-1]. - real, intent(out), dimension(:) :: dSV_dS !< The partial derivative of specific volume with - !! salinity [m3 kg-1 / Pa]. - integer, intent(in) :: start !< The starting point in the arrays. - integer, intent(in) :: npts !< The number of values to calculate. + real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface [degC]. + real, intent(in), dimension(:) :: S !< Salinity [PSU]. + real, intent(in), dimension(:) :: pressure !< pressure [Pa]. + real, intent(inout), dimension(:) :: dSV_dT !< The partial derivative of specific volume with + !! potential temperature [m3 kg-1 degC-1]. + real, intent(inout), dimension(:) :: dSV_dS !< The partial derivative of specific volume with + !! salinity [m3 kg-1 / Pa]. + integer, intent(in) :: start !< The starting point in the arrays. + integer, intent(in) :: npts !< The number of values to calculate. ! Local variables real :: al0, p0, lambda, I_denom @@ -377,15 +377,15 @@ end subroutine calculate_specvol_derivs_wright !! from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. !! Coded by R. Hallberg, 1/01 subroutine calculate_compress_wright(T, S, pressure, rho, drho_dp, start, npts) - real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface [degC]. - real, intent(in), dimension(:) :: S !< Salinity [PSU]. - real, intent(in), dimension(:) :: pressure !< pressure [Pa]. - real, intent(out), dimension(:) :: rho !< In situ density [kg m-3]. - real, intent(out), dimension(:) :: drho_dp !< The partial derivative of density with pressure - !! (also the inverse of the square of sound speed) - !! [s2 m-2]. - integer, intent(in) :: start !< The starting point in the arrays. - integer, intent(in) :: npts !< The number of values to calculate. + real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface [degC]. + real, intent(in), dimension(:) :: S !< Salinity [PSU]. + real, intent(in), dimension(:) :: pressure !< pressure [Pa]. + real, intent(inout), dimension(:) :: rho !< In situ density [kg m-3]. + real, intent(inout), dimension(:) :: drho_dp !< The partial derivative of density with pressure + !! (also the inverse of the square of sound speed) + !! [s2 m-2]. + integer, intent(in) :: start !< The starting point in the arrays. + integer, intent(in) :: npts !< The number of values to calculate. ! Coded by R. Hallberg, 1/01 ! Local variables @@ -428,18 +428,18 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, !! equation of state. real, intent(in) :: G_e !< The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2]. real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - intent(out) :: dpa !< The change in the pressure anomaly across the + intent(inout) :: dpa !< The change in the pressure anomaly across the !! layer [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer + optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer !! of the pressure anomaly relative to the anomaly !! at the top of the layer [Pa Z ~> Pa m]. real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), & - optional, intent(out) :: intx_dpa !< The integral in x of the difference between the + optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the !! pressure anomaly at the top and bottom of the !! layer divided by the x grid spacing [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & - optional, intent(out) :: inty_dpa !< The integral in y of the difference between the + optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the !! pressure anomaly at the top and bottom of the !! layer divided by the y grid spacing [Pa]. real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & @@ -629,18 +629,18 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & !! mathematically identical with different values of spv_ref, but this reduces the !! effects of roundoff. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(out) :: dza !< The change in the geopotential anomaly across + intent(inout) :: dza !< The change in the geopotential anomaly across !! the layer [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(out) :: intp_dza !< The integral in pressure through the layer of + optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of !! the geopotential anomaly relative to the anomaly !! at the bottom of the layer [Pa m2 s-2]. real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & - optional, intent(out) :: intx_dza !< The integral in x of the difference between the + optional, intent(inout) :: intx_dza !< The integral in x of the difference between the !! geopotential anomaly at the top and bottom of !! the layer divided by the x grid spacing [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & - optional, intent(out) :: inty_dza !< The integral in y of the difference between the + optional, intent(inout) :: inty_dza !< The integral in y of the difference between the !! geopotential anomaly at the top and bottom of !! the layer divided by the y grid spacing [m2 s-2]. integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate From 4d245d6b7448182c38c737a2e2f558f3965e692a Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 5 Feb 2020 13:59:46 -0500 Subject: [PATCH 5/5] Segment index added for several OBC segment flags In setup_[uv]_point_obc, a number of flags for each segment were assigned without index (OBC%segment%flag), which was causing all segments to be assigned the same value. This could problems with, say, `[uv]_value_needed` flags, which have different values depending on if they are E-W or N-S oriented. This patch adds the index to the assignments. --- src/core/MOM_open_boundary.F90 | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 822ca6486f..2fdfad4b59 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -904,8 +904,8 @@ subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y) OBC%segment(l_seg)%open = .true. OBC%Flather_u_BCs_exist_globally = .true. OBC%open_u_BCs_exist_globally = .true. - OBC%segment%z_values_needed = .true. - OBC%segment%u_values_needed = .true. + OBC%segment(l_seg)%z_values_needed = .true. + OBC%segment(l_seg)%u_values_needed = .true. elseif (trim(action_str(a_loop)) == 'ORLANSKI') then OBC%segment(l_seg)%radiation = .true. OBC%segment(l_seg)%open = .true. @@ -933,14 +933,14 @@ subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y) elseif (trim(action_str(a_loop)) == 'NUDGED') then OBC%segment(l_seg)%nudged = .true. OBC%nudged_u_BCs_exist_globally = .true. - OBC%segment%u_values_needed = .true. + OBC%segment(l_seg)%u_values_needed = .true. elseif (trim(action_str(a_loop)) == 'NUDGED_TAN') then OBC%segment(l_seg)%nudged_tan = .true. OBC%nudged_u_BCs_exist_globally = .true. - OBC%segment%v_values_needed = .true. + OBC%segment(l_seg)%v_values_needed = .true. elseif (trim(action_str(a_loop)) == 'NUDGED_GRAD') then OBC%segment(l_seg)%nudged_grad = .true. - OBC%segment%g_values_needed = .true. + OBC%segment(l_seg)%g_values_needed = .true. elseif (trim(action_str(a_loop)) == 'GRADIENT') then OBC%segment(l_seg)%gradient = .true. OBC%segment(l_seg)%open = .true. @@ -948,13 +948,13 @@ subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y) elseif (trim(action_str(a_loop)) == 'SIMPLE') then OBC%segment(l_seg)%specified = .true. OBC%specified_u_BCs_exist_globally = .true. ! This avoids deallocation - OBC%segment%u_values_needed = .true. + OBC%segment(l_seg)%u_values_needed = .true. elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then OBC%segment(l_seg)%specified_tan = .true. - OBC%segment%v_values_needed = .true. + OBC%segment(l_seg)%v_values_needed = .true. elseif (trim(action_str(a_loop)) == 'SIMPLE_GRAD') then OBC%segment(l_seg)%specified_grad = .true. - OBC%segment%g_values_needed = .true. + OBC%segment(l_seg)%g_values_needed = .true. else call MOM_error(FATAL, "MOM_open_boundary.F90, setup_u_point_obc: "//& "String '"//trim(action_str(a_loop))//"' not understood.") @@ -1045,8 +1045,8 @@ subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x) OBC%segment(l_seg)%open = .true. OBC%Flather_v_BCs_exist_globally = .true. OBC%open_v_BCs_exist_globally = .true. - OBC%segment%z_values_needed = .true. - OBC%segment%v_values_needed = .true. + OBC%segment(l_seg)%z_values_needed = .true. + OBC%segment(l_seg)%v_values_needed = .true. elseif (trim(action_str(a_loop)) == 'ORLANSKI') then OBC%segment(l_seg)%radiation = .true. OBC%segment(l_seg)%open = .true. @@ -1074,14 +1074,14 @@ subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x) elseif (trim(action_str(a_loop)) == 'NUDGED') then OBC%segment(l_seg)%nudged = .true. OBC%nudged_v_BCs_exist_globally = .true. - OBC%segment%v_values_needed = .true. + OBC%segment(l_seg)%v_values_needed = .true. elseif (trim(action_str(a_loop)) == 'NUDGED_TAN') then OBC%segment(l_seg)%nudged_tan = .true. OBC%nudged_v_BCs_exist_globally = .true. - OBC%segment%u_values_needed = .true. + OBC%segment(l_seg)%u_values_needed = .true. elseif (trim(action_str(a_loop)) == 'NUDGED_GRAD') then OBC%segment(l_seg)%nudged_grad = .true. - OBC%segment%g_values_needed = .true. + OBC%segment(l_seg)%g_values_needed = .true. elseif (trim(action_str(a_loop)) == 'GRADIENT') then OBC%segment(l_seg)%gradient = .true. OBC%segment(l_seg)%open = .true. @@ -1089,13 +1089,13 @@ subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x) elseif (trim(action_str(a_loop)) == 'SIMPLE') then OBC%segment(l_seg)%specified = .true. OBC%specified_v_BCs_exist_globally = .true. ! This avoids deallocation - OBC%segment%v_values_needed = .true. + OBC%segment(l_seg)%v_values_needed = .true. elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then OBC%segment(l_seg)%specified_tan = .true. - OBC%segment%u_values_needed = .true. + OBC%segment(l_seg)%u_values_needed = .true. elseif (trim(action_str(a_loop)) == 'SIMPLE_GRAD') then OBC%segment(l_seg)%specified_grad = .true. - OBC%segment%g_values_needed = .true. + OBC%segment(l_seg)%g_values_needed = .true. else call MOM_error(FATAL, "MOM_open_boundary.F90, setup_v_point_obc: "//& "String '"//trim(action_str(a_loop))//"' not understood.")