<p><b>mpetersen@lanl.gov</b> 2012-10-16 15:13:19 -0600 (Tue, 16 Oct 2012)</p><p>branch commit, basin_pbc_mrp: add Ilicak2_overflow case, and flags to choose between test cases for topography and initial conditions. Tested with the generation of 120km global mesh, matches previous ocean.nc using unix diff.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/basin_pbc_mrp/src/basin.F
===================================================================
--- branches/ocean_projects/basin_pbc_mrp/src/basin.F 2012-10-15 22:03:46 UTC (rev 2214)
+++ branches/ocean_projects/basin_pbc_mrp/src/basin.F 2012-10-16 21:13:19 UTC (rev 2215)
@@ -21,19 +21,13 @@
! Step 1: Set the number of Vertical levels
! Step 2: Set if the grid is on a sphere or not, and it's radius
! Step 3: Specify some Parameters
-! Step 4: Check the Initial conditions routine get_init_conditions
-! Step 5: Check the depth routine define_kmt
+! Step 4: Check get_init_conditions routine for initial T&S, thickness, etc.
+! Step 5: Check define_kmt routine for bottomDepth and kmt (maxLevelCell) variables
+! Step 6: Check get_dz routine for hZLevel variable
!
-!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
-integer, parameter :: nx = 50
-! This needs to be changed for correct periodic boundaries
-! Lx is the TOTAL domain width, and needs to be exact for correct periodic
-! boundaries in x.
-real, parameter :: Lx = 3200.0e3 ! 40x80km=3200km
-
! original grid variables
integer :: time, nCells, nEdges, nVertices
integer :: maxEdges, maxEdges2, TWO, vertexDegree, nVertLevels
@@ -65,21 +59,20 @@
real, dimension(40) :: dz
integer :: nMonths = 1
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
! Step 1: Set the number of Vertical levels
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
integer, parameter :: nVertLevelsMOD = 40
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! basin-mod
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
-! Need to uncomment the options appropriate for the input grid file. If it's on
-! a sphere, specify the flag "on_a_sphere" and "sphere_radius". Otherwise set
-! them equal to NO and 0.0 respectively
+! Step 2: Set if the grid is on a sphere or not, and it's radius
!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! Step 2: Set if the grid is on a sphere or not, and it's radius
character (len=16) :: on_a_sphere = 'YES '
real*8, parameter :: sphere_radius = 6.37122e6
!real*8, parameter :: sphere_radius = 1.0
@@ -87,32 +80,50 @@
!character (len=16) :: on_a_sphere = 'NO '
!real*8, parameter :: sphere_radius = 0.0
-logical, parameter :: real_bathymetry=.true.
+! zLevel thickness options:
+! 'POP_40_zLevel', 'equally_spaced', 'isopycnal'
+character (len=32) :: zLevel_thickness = 'POP_40_zLevel'
+
+! bottom topography options:
+! 'realistic_ETOPO', 'flat_bottom', 'Ilicak2_overflow'
+character (len=32) :: bottom_topography = 'realistic_ETOPO'
+
+! initial temperature and salinity options:
+! 'realistic_WOCE', 'constant', 'Ilicak2_overflow'
+character (len=32) :: initial_conditions = 'realistic_WOCE'
+
logical, parameter :: eliminate_inland_seas=.true.
logical, parameter :: l_woce = .true.
-logical, parameter :: write_OpenDX_flag = .true.
+logical, parameter :: write_OpenDX_flag = .false.
logical, parameter :: monthly_forcing = .false.
! Set the number of top layers that are not allowed to have land, usually three.
integer, parameter :: top_layers_without_land = 3
-
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
! Step 3: Specify some Parameters
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real (kind=8), parameter :: &
- h_total_max = 2000.0, &
- u_max = 0.0, &
+ h_total_max = 2000.0, & ! total layer thickness, for equally spaced case
+ u_max = 0.0, & ! max initial velocity
u_src_max = 0.1, & ! max wind stress, N/m2
- beta = 1.4e-11, &
- f0 = -1.1e-4, &
- omega = 7.29212e-5
+ f0 = -1.1e-4, & ! Coriolis parameter
+ beta = 1.4e-11, &
+ omega = 7.29212e-5 ! rotation rate of earth
+! This needs to be changed for correct periodic boundaries
+! Lx is the TOTAL domain width, and needs to be exact for correct periodic
+! boundaries in x.
+ real (kind=8), parameter :: Lx = 3200.0e3 ! 40x80km=3200km
+
+
real (kind=8) :: ymid, ytmp, ymax, xmid, xloc, yloc, pert, ymin, distance, r, c1(3), c2(3)
real (kind=8) :: latmid, lattmp, latmax, latmin
integer :: cell1, cell2
-
! new grid variables
real, dimension(nVertLevelsMOD) :: hZLevel, refBottomDepth
integer :: nCellsNew, nEdgesNew, nVerticesNew
@@ -430,7 +441,7 @@
end subroutine write_graph
-!Step 4: Check the Initial conditions routine get_init_conditions
+! Step 4: Check get_init_conditions routine for initial T&S, thickness, etc.
subroutine get_init_conditions
implicit none
real :: halfwidth, dtr, pi, p(3), q(3), xin, yin, zin, ulon, ulat, stress, n1, n2, distance, r, temp_t, temp_s
@@ -448,18 +459,53 @@
uNew = 0
vNew = 0
-if(.not.real_bathymetry) then
+if (initial_conditions.eq.'uniform_TS') then
- write(6,*) ' not using real bathymetry'
+ temperatureNew = 10.0
+ salinityNew = 34.0
+ tracer1New = 1.0
+ uNew = 0.0
+ u_srcNew = 0.0
+ ! change thickness of bottom level
+ do iCell=1,nCellsNew
+ do k=1, nVertLevelsMod
+ hNew(1,k,iCell) = hZLevel(k)
+ enddo
+ enddo
+
+elseif (initial_conditions.eq.'Ilicak2_overflow') then
+
+ write(6,*) ' setting temperature'
+ do i = 1,nCellsNew
+ if(yCellNew(i) < 20000) then
+ temperatureNew(1,:,i) = 10.0
+ else
+ temperatureNew(1,:,i) = 20.0
+ endif
+ enddo
+
+ salinityNew(1,:,:) = 35.0
+ Tracer1New(1,:,:) = 1.0
+ uNew = 0.0
+ u_srcNew = 0.0
+
+ do iCell=1,nCellsNew
+ do k=1,maxLevelCellNew(iCell)
+ hNew(1,k,iCell) = hZLevel(k)
+ enddo
+ do k=maxLevelCellNew(iCell)+1, nVertLevelsMod
+ hNew(1,k,iCell) = 0.0
+ enddo
+ enddo
+
+elseif (initial_conditions.eq.'isopycnal') then
+
fEdgeNew(:) = 0.0
fVertexNew(:) = 0.0
bottomDepthNew(:) = 0.0
uNew(:,:,:) = 0.0
- vNew(:,:,:) = 0.0
-!!!!!!!!!!!! mrp need to define hZLevel in get_dz, below.
-
! basin-mod
! setting for three levels - Set h values for isopycnal system
write(6,*) ' setting three levels for isopycnal system'
@@ -472,15 +518,8 @@
hNew(1,1,:) = 3250.0
bottomDepthNew(:) = -( hNew(1,1,:) )
endif
- hZLevel = hNew(1,:,1)
-!!!!!!!!!!! remove these lines mrp
- refBottomDepth(1) = hZLevel(1)
- do k = 2,nVertLevels
- refBottomDepth(k) = refBottomDepth(k-1) + hZLevel(k)
- end do
-
! basin-mod
! Noise is meant to make the flow unstable at some point
! Not needed for all simulations
@@ -621,33 +660,10 @@
endif
write(6,*) ' u_srcNew ', minval(u_srcNew), maxval(u_srcNew)
- ! basin-mod
- ! set coriolis parameter for grid
- write(6,*) ' setting Coriolis parameter'
- if(on_a_sphere.eq.'YES ') then
- do i = 1,nVerticesNew
- fVertexNew(i) = 2.0 * omega * sin(latVertexNew(i))
- enddo
- do i = 1,nEdgesNew
- fEdgeNew(i) = 2.0 * omega * sin(latEdgeNew(i))
- enddo
- else
- do i = 1,nVerticesNew
- fVertexNew(i) = f0 + (yVertexNew(i) - ymid) * beta
- enddo
+elseif (initial_conditions.eq.'realistic_WOCE') then
- do i = 1,nEdgesNew
- fEdgeNew(i) = f0 + (yEdgeNew(i) - ymid) * beta
- enddo
- endif
- write(6,*) ' done not real bathymetry'
-endif ! if(.not.real_bathymetry) then
-
-
-if(real_bathymetry) then
-
u_srcNew = 0.0
windStressMonthlyNew = 0.0
do iEdge=1,nEdgesNew
@@ -890,8 +906,39 @@
! salinityRestoreMonthlyNew(iMonth,:) = salinityRestoreNew(:)
!enddo
-endif ! real_bathymetry
+else
+ print *, ' Incorrect choice of initial_conditions: ',initial_conditions
+ stop
+
+endif ! initial_conditions
+
+
+if (trim(initial_conditions).ne.'realistic_WOCE') then
+
+ ! basin-mod
+ ! set coriolis parameter for grid
+ write(6,*) ' setting Coriolis parameter'
+ if(on_a_sphere.eq.'YES ') then
+ do i = 1,nVerticesNew
+ fVertexNew(i) = 2.0 * omega * sin(latVertexNew(i))
+ enddo
+
+ do i = 1,nEdgesNew
+ fEdgeNew(i) = 2.0 * omega * sin(latEdgeNew(i))
+ enddo
+ else
+ do i = 1,nVerticesNew
+ fVertexNew(i) = f0 + (yVertexNew(i) - ymid) * beta
+ enddo
+
+ do i = 1,nEdgesNew
+ fEdgeNew(i) = f0 + (yEdgeNew(i) - ymid) * beta
+ enddo
+ endif
+
+endif
+
write(6,*) ' done get_init_conditions'
end subroutine get_init_conditions
@@ -1346,14 +1393,14 @@
end subroutine write_grid
-! Step 5: Check the depth routine define_kmt
+! Step 5: Check define_kmt routine for bottomDepth and kmt (maxLevelCell) variables
subroutine define_kmt
implicit none
real (kind=4), allocatable, dimension(:) :: x,y, work_kmt
real (kind=4), allocatable, dimension(:,:) :: ztopo
integer :: nx, ny, inx, iny, ix, iy, kmt_neighbor_max
real :: pi, dtr, zdata, rlon, rlat, r, ymin, ymax, xmin, xmax
-real :: latmin, latmax, lonmin, lonmax
+real :: latmin, latmax, lonmin, lonmax, ridgeDepth
logical :: flag, kmt_flag
pi = 4.0*atan(1.0)
dtr = pi / 180.0
@@ -1361,7 +1408,7 @@
allocate(kmt(nCells))
kmt = 0
-if(.not.real_bathymetry) then
+if (trim(bottom_topography).ne.'realistic_ETOPO') then
kmt = nVertLevelsMOD
if(on_a_sphere.eq.'YES ') then
write(6,*) 'Working on a sphere'
@@ -1399,9 +1446,12 @@
where(kmt.eq.0) work_kmt=1.0
write(6,*) 'number of cells culled ',sum(work_kmt)
deallocate(work_kmt)
+
endif
-if(real_bathymetry) then
+
+if (bottom_topography.eq.'realistic_ETOPO') then
+
nx = 10800
ny = 5400
allocate(x(nx))
@@ -1470,6 +1520,45 @@
deallocate(x)
deallocate(y)
deallocate(ztopo)
+
+elseif (bottom_topography.eq.'Ilicak2_overflow') then
+
+ kmt = nVertLevelsMOD
+
+ ridgeDepth = 500.0
+ do iCell = 1,nCells
+ ! From Mehmet Ilicak:
+ ! depth=2000
+ ! val1 = 500 is top of ridge
+ ! h(i,j) = val1 + 0.5*(depth-val1) * (1.0+TANH((lon(i,j)-40000.0)/7000.0))
+ bottomDepth(iCell) = ridgeDepth + 0.5*(h_total_max-ridgeDepth) * (1.0+tanh((yCell(iCell)-40000.0)/7000.0))
+
+ if (bottomDepth(iCell).lt.0.0.or. &
+ bottomDepth(iCell).gt.refBottomDepth(nVertLevelsMOD)) then
+ print *, 'error: bottomDepth cannot be less than zero or greater than refBottomDepth(nVertLevels)'
+ print *, 'iCell, bottomDepth(iCell):', iCell, bottomDepth(iCell)
+ exit
+ end if
+
+ do k=1,nVertLevelsMOD
+ if (bottomDepth(iCell).le.refBottomDepth(k)) then
+ kmt(iCell) = k
+ exit
+ endif
+ end do
+
+ enddo
+
+
+elseif (bottom_topography.eq.'flat_bottom') then
+
+ kmt = nVertLevelsMOD
+
+else
+
+ print *, ' Incorrect choice of bottom_topography: ',bottom_topography
+ stop
+
endif
! Eliminate isolated ocean cells, and make these isolated deep cells
@@ -1921,20 +2010,27 @@
end subroutine map_connectivity
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! Step 6: Check get_dz routine for hZLevel variable
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine get_dz
integer k
-if(.not.real_bathymetry) then
+if (zLevel_thickness.eq.'isopycnal') then
-!!!!!!!!!! mrp add cases here, and a step to change this.
+ ! hZLevel not used for isopycnal setting, just set to zero.
+ hZLevel = 0.0
- write(6,*) ' equally spaced layers'
+elseif (zLevel_thickness.eq.'equally_spaced') then
+
+ write(6,*) ' equally spaced zLevels'
do i = 1, nVertLevelsMOD
hZLevel(i) = h_total_max / nVertLevelsMOD
end do
-else
+elseif(zLevel_thickness.eq.'POP_40_zLevel') then
dz( 1) = 1001.244 ! 5.006218 10.01244
dz( 2) = 1011.258 ! 15.06873 20.12502
@@ -1981,6 +2077,11 @@
hZLevel = dz
+else
+
+ print *, ' Incorrect choice of zLevel_thickness: ',zLevel_thickness
+ stop
+
endif
refBottomDepth(1) = hZLevel(1)
</font>
</pre>