Merge branch 'MFJansen-dev/gfdl' into dev/gfdl

src/user/Neverland_initialization.F90

@@ -15,6 +15,8 @@ module Neverland_initialization
 use MOM_verticalGrid, only : verticalGrid_type
 use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type
 
+use random_numbers_mod, only: initializeRandomNumberStream, getRandomNumbers, randomNumberStream
+
 implicit none ; private
 
 #include <MOM_memory.h>
@@ -122,6 +124,11 @@ subroutine Neverland_initialize_thickness(h, G, GV, US, param_file, eqn_of_state
                             ! usually negative because it is positive upward.
   real, dimension(SZK_(G)) :: h_profile ! Vector of initial thickness profile (Z)
   real :: e_interface ! Current interface position (m)
+  real :: x,y,r1,r2 ! x,y and radial coordinates for computation of initial pert.
+  real :: pert_amp ! Amplitude of perturbations measured in Angstrom_H
+  real :: h_noise ! Amplitude of noise to scale h by
+  real :: noise ! Noise
+  type(randomNumberStream) :: rns ! Random numbers for stochastic tidal parameterization
   character(len=40)  :: mdl = "Neverland_initialize_thickness" ! This subroutine's name.
   integer :: i, j, k, k1, is, ie, js, je, nz, itt
 
@@ -131,6 +138,11 @@ subroutine Neverland_initialize_thickness(h, G, GV, US, param_file, eqn_of_state
   call get_param(param_file, mdl, "INIT_THICKNESS_PROFILE", h_profile, &
                  "Profile of initial layer thicknesses.", units="m", scale=US%m_to_Z, &
                  fail_if_missing=.true.)
+  call get_param(param_file, mdl, "NL_THICKNESS_PERT_AMP", pert_amp, &
+                 "Amplitude of finite scale perturbations as fraction of depth.", &
+                 units="nondim", default=0.)
+  call get_param(param_file, mdl, "NL_THICKNESS_NOISE_AMP", h_noise, &
+                 "Amplitude of noise to scale layer by.", units="nondim", default=0.)
 
   ! e0 is the notional position of interfaces
   e0(1) = 0. ! The surface
@@ -140,10 +152,24 @@ subroutine Neverland_initialize_thickness(h, G, GV, US, param_file, eqn_of_state
 
   do j=js,je ; do i=is,ie
     e_interface = -G%bathyT(i,j)
-    do k=nz,1,-1
-      h(i,j,k) = max( GV%Angstrom_H, GV%Z_to_H * (e0(k) - e_interface) )
-      e_interface = max( e0(k), e_interface - GV%H_to_Z * h(i,j,k) )
+    do k=nz,2,-1
+      h(i,j,k) = GV%Z_to_H * (e0(k) - e_interface) ! Nominal thickness
+      x=(G%geoLonT(i,j)-G%west_lon)/G%len_lon
+      y=(G%geoLatT(i,j)-G%south_lat)/G%len_lat
+      r1=sqrt((x-0.7)**2+(y-0.2)**2)
+      r2=sqrt((x-0.3)**2+(y-0.25)**2)
+      h(i,j,k) = h(i,j,k) + pert_amp*(e0(k) - e0(nz+1))*GV%Z_to_H*(spike(r1,0.15)-spike(r2,0.15)) ! Prescribed perturbation
+      if (h_noise /= 0.) then
+        rns = initializeRandomNumberStream( int( 4096*(x + (y+1.)) ) )
+        call getRandomNumbers(rns, noise) ! x will be in (0,1)
+        noise = h_noise * 2. * ( noise - 0.5 ) ! range -h_noise to h_noise
+        h(i,j,k) = ( 1. + noise ) * h(i,j,k)
+      endif
+      h(i,j,k) = max( GV%Angstrom_H, h(i,j,k) ) ! Limit to non-negative
+      e_interface = e_interface + GV%H_to_Z * h(i,j,k) ! Actual position of upper interface
     enddo
+    h(i,j,1) = GV%Z_to_H * (e0(1) - e_interface) ! Nominal thickness
+    h(i,j,1) = max( GV%Angstrom_H, h(i,j,1) ) ! Limit to non-negative
   enddo ; enddo
 
 end subroutine Neverland_initialize_thickness