<p><b>dwj07@fsu.edu</b> 2012-12-04 10:57:26 -0700 (Tue, 04 Dec 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Updating temporal convergence test case to work with newer versions of basin.<br>
<br>
Changing how the setup scripts work.<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/BASIN-namlist.basin.template
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/BASIN-namlist.basin.template (rev 0)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/BASIN-namlist.basin.template 2012-12-04 17:57:26 UTC (rev 2339)
@@ -0,0 +1,24 @@
+&basin
+ nVertLevelsMOD = *VERTLEVS
+ on_a_sphere = 'NO'
+ sphere_radius = 0.0
+ expand_from_unit_sphere = .false.
+ zLevel_thickness = 'equally_spaced'
+ bottom_topography = 'flat_bottom'
+ initial_conditions = 'none'
+ eliminate_inland_seas=.false.
+ load_woce_IC = .false.
+ write_OpenDX_flag = .false.
+ monthly_forcing = .false.
+ check_mesh = .true.
+ cut_domain_from_sphere = .false.
+ solid_boundary_in_y = .true.
+ solid_boundary_in_x = .false.
+ top_layers_without_land = 3
+ h_total_max = 1000.0
+ u_src_max = 0.0
+ f0 = 0.0
+ beta = 0.0
+ omega = 0.0
+ Lx = 160.0e3
+/
Modified: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/Makefile.end
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/Makefile.end 2012-12-03 21:05:55 UTC (rev 2338)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/Makefile.end 2012-12-04 17:57:26 UTC (rev 2339)
@@ -15,18 +15,20 @@
.SUFFIXES: .F .o
-
OBJS = basin.o \
utilities.o \
module_read_netcdf.o \
module_read_topo.o \
module_read_TS.o \
+ module_read_monthly.o \
module_cullLoops.o \
module_write_netcdf.o
+
all: map
-basin.o: utilities.o module_write_netcdf.o module_read_netcdf.o module_read_topo.o module_read_TS.o module_cullLoops.o
+#basin.o: utilities.o module_write_netcdf.o module_read_netcdf.o module_read_topo.o module_read_TS.o module_cullLoops.o
+basin.o: utilities.o module_write_netcdf.o module_read_netcdf.o module_read_topo.o module_read_TS.o module_read_monthly.o module_cullLoops.o
map: $(OBJS)
$(SFC) $(LDFLAGS) -o $@ $(OBJS) $(LIBS)
Added: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/Makefile.front
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/Makefile.front (rev 0)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/Makefile.front 2012-12-04 17:57:26 UTC (rev 2339)
@@ -0,0 +1,10 @@
+FILE_OFFSET = -DOFFSET64BIT
+gfortran:
+ ( $(MAKE) all \
+ "SFC = gfortran" \
+ "SCC = gcc" \
+ "FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8 -fconvert=big-endian -ffree-form" \
+ "CFLAGS = -O3 -m64" \
+ "LDFLAGS = -O3 -m64" \
+ "CORE = $(CORE)" \
+ "CPPFLAGS = $(MODEL_FORMULATION) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
Deleted: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/basin-template.F
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/basin-template.F 2012-12-03 21:05:55 UTC (rev 2338)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/basin-template.F 2012-12-04 17:57:26 UTC (rev 2339)
@@ -1,1633 +0,0 @@
-program map_to_basin
-
-use read_netcdf
-use read_topo
-use read_TS
-use write_netcdf
-use utilities
-use cullLoops
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! Program: basin.F
-!
-! This program is meant to add land to grids, as well as initial conditions.
-!
-! This program is used to take a specific mesh, and remove Cells from it
-! It can be used to change a planar grid into a Channel or a basin grid, or to
-! Change a spherical grid into a Limited area spherical grid.
-!
-! How to use:
-! 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
-!
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-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 = 160.0e3
-
-! original grid variables
-integer :: time, nCells, nEdges, nVertices
-integer :: maxEdges, maxEdges2, TWO, vertexDegree, nVertLevels
-integer, allocatable, dimension(:) :: indexToCellID, indexToEdgeID, indexToVertexID
-real, allocatable, dimension(:) :: xCell, yCell, zCell, latCell, lonCell, meshDensity
-real, allocatable, dimension(:) :: xEdge, yEdge, zEdge, latEdge, lonEdge
-real, allocatable, dimension(:) :: xVertex, yVertex, zVertex, latVertex, lonVertex
-integer, allocatable, dimension(:) :: nEdgesOnCell, nEdgesOnEdge
-integer, allocatable, dimension(:,:) :: cellsOnCell, edgesOnCell, verticesOnCell
-integer, allocatable, dimension(:,:) :: cellsOnEdge, verticesOnEdge, edgesOnEdge
-integer, allocatable, dimension(:,:) :: cellsOnVertex, edgesOnVertex
-real, allocatable, dimension(:) :: areaCell, areaTriangle, dcEdge, dvEdge, angleEdge
-real, allocatable, dimension(:,:) :: kiteAreasOnVertex, weightsOnEdge
-
-real, allocatable, dimension(:) :: fEdge, fVertex, h_s, work1
-real, allocatable, dimension(:,:) :: u_src
-real, allocatable, dimension(:,:,:) :: u, v, h
-real, allocatable, dimension(:,:,:) :: rho
-
-integer nlon, nlat, ndepth
-real(kind=4), allocatable, dimension(:) :: t_lon, t_lat, depth_t
-real(kind=4), allocatable, dimension(:,:) :: mTEMP, mSALT
-real(kind=4), allocatable, dimension(:,:,:) :: TEMP, SALT
-real(kind=4), allocatable, dimension(:,:) :: TAUX, TAUY
-
-real, dimension(40) :: dz
-
-! Step 1: Set the number of Vertical levels
-integer, parameter :: nVertLevelsMOD = *VERTLEVS
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! 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
-!character (len=16) :: on_a_sphere = 'YES '
-!real*8, parameter :: sphere_radius = 6.37122e6
-!real*8, parameter :: sphere_radius = 1.0
-
-character (len=16) :: on_a_sphere = 'NO '
-real*8, parameter :: sphere_radius = 0.0
-
-logical, parameter :: real_bathymetry=.false.
-logical, parameter :: eliminate_inland_seas=.false.
-logical, parameter :: l_woce = .false.
-
-
-! Step 3: Specify some Parameters
- real (kind=8), parameter :: &
- h_total_max = 1000.0, &
- u_max = 0.0, &
- u_src_max = 0.0, & ! max wind stress, N/m2
- beta = 0.0, &
- f0 = 1.2e-4, &
- omega = 0.0
-
- 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
-integer :: nCellsNew, nEdgesNew, nVerticesNew
-integer :: maxEdgesNew, maxEdges2New, TWONew, vertexDegreeNew, nVertLevelsNew
-integer, allocatable, dimension(:) :: indexToCellIDNew, indexToEdgeIDNew, indexToVertexIDNew
-real, allocatable, dimension(:) :: xCellNew, yCellNew, zCellNew, latCellNew, lonCellNew, meshDensityNew
-real, allocatable, dimension(:) :: xEdgeNew, yEdgeNew, zEdgeNew, latEdgeNew, lonEdgeNew
-real, allocatable, dimension(:) :: xVertexNew, yVertexNew, zVertexNew, latVertexNew, lonVertexNew
-integer, allocatable, dimension(:) :: nEdgesOnCellNew, nEdgesOnEdgeNew, flipVerticesOnEdgeOrdering
-integer, allocatable, dimension(:,:) :: cellsOnCellNew, edgesOnCellNew, verticesOnCellNew
-integer, allocatable, dimension(:,:) :: cellsOnEdgeNew, verticesOnEdgeNew, edgesOnEdgeNew
-integer, allocatable, dimension(:,:) :: cellsOnVertexNew, edgesOnVertexNew
-integer, allocatable, dimension(:,:) :: boundaryEdgeNew, boundaryVertexNew
-real, allocatable, dimension(:) :: areaCellNew, areaTriangleNew, dcEdgeNew, dvEdgeNew, angleEdgeNew
-real, allocatable, dimension(:,:) :: kiteAreasOnVertexNew, weightsOnEdgeNew, normalsNew
-
-real, allocatable, dimension(:) :: fEdgeNew, fVertexNew, h_sNew, hZLevel
-real, allocatable, dimension(:,:) :: u_srcNew
-real, allocatable, dimension(:,:,:) :: uNew, vNew, hNew
-real, allocatable, dimension(:,:,:) :: rhoNew, temperatureNew, salinityNew, tracer1New
-real, allocatable, dimension(:) :: temperatureRestoreNew, salinityRestoreNew
-
-! mapping variables
-integer, allocatable, dimension(:) :: kmt, maxLevelCellNew
-integer, allocatable, dimension(:) :: cellMap, edgeMap, vertexMap
-real, allocatable, dimension(:) :: depthCell
-
-! work variables
-integer :: i,j,jNew,k,jEdge,jEdgeNew,iVertex1New,iVertex2New,iCell1New,iCell2New
-integer :: iCell, iCell1, iCell2, iCell3, iEdge, iVertex, iVertex1, iVertex2
-integer :: iCellNew, iEdgeNew, iVertexNew, ndata, jCell1, jCell2, jCell, iter
-real :: xin, yin, zin, ulon, ulat, ux, uy, uz, rlon, rlat, temp_t, temp_s
-
-! get to work
-write(6,*) ' starting'
-write(6,*)
-
-! get depth profile for later
-!call get_dz
-
-! get grid
-write(6,*) ' calling read_grid'
-write(6,*)
-call read_grid
-write(6,*) ' xCell 1: ',minval(xCell), maxval(xCell)
-
-! copy dimensions
-write(6,*) ' copy dimensions'
-write(6,*)
-call copy_dimensions
-write(6,*) ' xCell 1: ',minval(xCell), maxval(xCell)
-
-! define the kmt array
-write(6,*) ' calling define_kmt'
-write(6,*)
-call define_kmt
-
-! define the mapping between original and new cells, edges and vertices
-write(6,*) ' calling define_mapping'
-write(6,*)
-call define_mapping
-
-! copy the vector arrays form the original to new arrays
-write(6,*) ' calling map_vectors'
-write(6,*)
-call map_vectors
-
-! define the new connectivity variables
-write(6,*) ' calling map_connectivity'
-write(6,*)
-call map_connectivity
-
-! check the mesh
-call error_checking
-
-!write(6,*) ' getting woce t and s '
-!call read_TS_init(nlon, nlat, ndepth)
-!write(6,*) ' TS INIT ', nlon, nlat, ndepth
-!allocate(t_lon(nlon), t_lat(nlat), depth_t(ndepth), TEMP(nlon,nlat,ndepth), SALT(nlon,nlat,ndepth))
-!allocate(TAUX(nlon,nlat), TAUY(nlon,nlat))
-!allocate(mTEMP(nlat,ndepth), mSALT(nlat,ndepth))
-!call read_TS_fields(t_lon, t_lat, depth_t, TEMP, SALT, TAUX, TAUY)
-!call read_TS_finalize()
-!do k=1,ndepth
-! ndata = 0; temp_t=0; temp_s=0
-! do j=1,nlat
-! do i=1,nlon
-! if(TEMP(i,j,k).gt.-10.0) then
-! ndata = ndata + 1
-! temp_t = temp_t + TEMP(i,j,k)
-! temp_s = temp_s + SALT(i,j,k)
-! endif
-! enddo
-! enddo
-! mTEMP(:,k) = temp_t / float(ndata)
-! mSALT(:,k) = temp_s / float(ndata)
-! write(6,*) ndata,mTemp(1,k),mSalt(1,k)
-!enddo
-
-! generate initial conditions
-call get_init_conditions
-
-! dump new grid to netCDF
-write(6,*) ' calling write_grid'
-write(6,*)
-call write_grid
-
-! dump graph for partioning
-write(6,*) ' call write_graph'
-write(6,*)
-call write_graph
-
-! write OpenDx file
-!write(6,*) ' calling write_OpenDX'
-!write(6,*)
- call write_OpenDX( on_a_sphere, &
- nCellsNew, &
- nVerticesNew, &
- nEdgesNew, &
- vertexDegreeNew, &
- maxEdgesNew, &
- xCellNew, &
- yCellNew, &
- zCellNew, &
- xVertexNew, &
- yVertexNew, &
- zVertexNew, &
- xEdgeNew, &
- yEdgeNew, &
- zEdgeNew, &
- nEdgesOnCellNew, &
- verticesOnCellNew, &
- verticesOnEdgeNew, &
- cellsOnVertexNew, &
- edgesOnCellNew, &
- areaCellNew, &
- maxLevelCellNew, &
- depthCell, &
- temperatureNew(1,1,:), &
- kiteAreasOnVertexNew )
-
-
-!do iCell=1,nCellsNew
- !ulon = 1.0; ulat = 0.0
- !xin = xCellNew(iCell); yin = yCellNew(iCell); zin = zCellNew(iCell)
- !call transform_from_lonlat_to_xyz(xin, yin, zin, ulon, ulat, ux, uy, uz)
- !if(abs(ux).lt.1.0e-10) ux=0.0
- !if(abs(uy).lt.1.0e-10) uy=0.0
- !if(abs(uz).lt.1.0e-10) uz=0.0
- !write(20,10) ux, uy, uz
- !10 format(3e25.10)
-!enddo
-
-write(6,*) ' finished'
-
-contains
-
-subroutine write_graph
-implicit none
-integer :: m,itmp(maxEdgesNew),k
-
- m=nEdgesNew
- do i=1,nCellsNew
- do j=1,nEdgesOnCellNew(i)
- if(cellsOnCellNew(j,i).eq.0) m=m-1
- enddo
- enddo
-
- open(42,file='graph.info',form='formatted')
- write(42,*) nCellsNew, m
- do i=1,nCellsNew
- itmp = 0; k = 0;
- do j=1,nEdgesOnCellNew(i)
- if(cellsOnCellNew(j,i).gt.0) then
- k=k+1; itmp(k)=cellsOnCellNew(j,i)
- endif
- enddo
- write(42,'(1x,12i8)',advance='no') (itmp(m),m=1,k)
- write(42,'(1x)')
- end do
- close(42)
-end subroutine write_graph
-
-
-!Step 4: Check the Initial conditions routine get_init_conditions
-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
-real :: dotProd
-real :: surfMaxTemp, surfMinTemp
-real :: y_0, x_0, x_1, x_2, x_3, width, cff1
-real :: yTrans, vertTempChange, maxTemp
-real, dimension(:), allocatable :: vertCoord
-real, dimension(:,:,:), allocatable :: tempHolder
-integer :: iTracer, ix, iy, ndata, i, j, k, ixt, iyt, ncull, jcount, iNoData, kdata(nVertLevelsMod)
-logical :: flag_lat
-
-pi = 4.0*atan(1.0)
-dtr = pi/180.0
-
-hNew = 100.0
-temperatureNew = 1.0
-salinityNew = 1.0
-tracer1New = 1.0
-uNew = 0
-vNew = 0
-
-allocate(vertCoord(nVertLevelsMOD))
-allocate(tempHolder(1, nVertLevelsMOD, nCellsNew))
-
-if(.not.real_bathymetry) then
-
- write(6,*) ' not using real bathymetry'
-
- fEdgeNew(:) = 0.0
- fVertexNew(:) = 0.0
- h_sNew(:) = 0.0
- uNew(:,:,:) = 0.0
- vNew(:,:,:) = 0.0
-
- ! basin-mod
- ! setting for three levels - Set h values for isopycnal system
- write(6,*) ' setting three levels for isopycnal system'
- vertCoord = 0
- do i = 1, nVertLevelsMOD
- hNew(1,i,:) = h_total_max / nVertLevelsMOD
- hZLevel(i) = h_total_max / nVertLevelsMOD
- vertCoord(i) = i * h_total_max / nVertLevelsMOD
- end do
-
- h_sNew(:) = -h_total_max
-
- ! basin-mod
- !Specify Density values for isopycnal levels
- write(6,*) ' setting density'
- rhoNew(1,:,:) = 1000.0
-! if(nVertLevelsMOD .eq. 3) then
-! rhoNew(1,2,:) = 1011.0
-! rhoNew(1,3,:) = 1012.0
-! endif
-
- ! basin-mod
- ! set temperature for isopycnal levels
- write(6,*) ' setting temperature'
- vertTempChange = 2.5
- surfMaxTemp = 13.1
- surfMinTemp = 10.1
- do k = 1, nVertLevelsMOD
- temperatureNew(1,k,:) = (surfMaxTemp - surfMinTemp) * ((-vertCoord(k)+h_total_max)/h_total_max)+ surfMinTemp
- enddo
-
- tempHolder = temperatureNew
-
- y_0 = 200.0e3
- x_0 = 0.0e3
- x_1 = 160.0e3
- x_2 = 110.0e3
- x_3 = 130.0e3
- width = 40.0e3
- do i = 1, nCellsNew
- cff1 = width * sin (6.0 * 3.141592 * (xCellNew(i) - x_0)/(x_1 - x_0))
-
- if( yCellNew(i) < y_0 - cff1 ) then
- do k = 1, nVertLevelsMOD
- temperatureNew(1,k,i) = temperatureNew(1,k,i) - 1.2
- end do
- else if( yCellNew(i) .ge. y_0 - cff1 .and. yCellNew(i) .le. y_0 - cff1+width) then
- do k = 1, nVertLevelsMOD
- temperatureNew(1,k,i) = tempHolder(1,k,i) - 1.2*(1.0 -( yCellNew(i) - (y_0 - cff1)) / (1.0 * width))
- end do
- endif
- enddo
-
- do i = 1, nCellsNew
- cff1 = 0.5 * width * sin(1.0 * 3.141592 * (xCellNew(i) - x_2)/(x_3 - x_2))
-
- if( yCellNew(i) .ge. y_0 - cff1-0.5*width .and. yCellNew(i) .le. y_0 - cff1+0.5*width .and. xCellNew(i) .ge. x_2 .and. xCellNew(i) .le. x_3) then
- do k = 1, nVertLevelsMOD
- temperatureNew(1,k,i) = temperatureNew(1,k,i) + 0.3 * (1.0 - ( (yCellNew(i)-(y_0-cff1))/(0.5*width)))
- end do
- endif
- end do
-! do i = 1,nCellsNew
-! yTrans = 250000.0 + 25000.0 * sin(2.0 * xCellNew(i)/15000.0)
-
-! if(yCellNew(i) > yTrans + 25000.0) then
-! maxTemp = 13.0
-! else if (yCellNew(i) < yTrans - 25000.0) then
-! maxTemp = 11.8
-! else
-
-! maxTemp = surfMinTemp + (surfMaxTemp - surfMinTemp) * (yCellNew(i) - (yTrans - 25000.0))/50000.0
-! print *,'maxTemp', maxTemp, (yCellNew(i) - (yTrans - 25000.0))/50000.0
-! endif
-
-
-! do k = 1,nVertLevelsMOD
-! temperatureNew(1,k,i) = maxTemp - ((k-1)*vertTempChange/nVertLevelsMOD)
-! end do
-! enddo
-
- ! basin-mod
- ! set salinity for levels
- salinityNew(1,:,:) = 35.0
-
- ! Updating density with linear EOS
- do i = 1,nCellsNew
- rhoNew(1,:,i) = 1000.0*(1.0 - 2.5e-4*temperatureNew(1,1,i) + 7.6e-4*salinityNew(1,1,i))
- enddo
-
- ! basin-mod
- ! set forcing for isopycnal levels
- write(6,*) 'setting u_src - wind forcing'
- u_srcNew = 0.0
- write(6,*) ' u_srcNew ', minval(u_srcNew), maxval(u_srcNew)
-
- ! basin-mod
- ! set coriolis parameter for grid
- write(6,*) ' setting Coriolis parameter'
- ymid = (maxval(yVertexNew(:)) - minval(yVertexNew(:)))/2.0
- do i = 1,nVerticesNew
- fVertexNew(i) = f0 + (yVertexNew(i) - ymid) * beta
- enddo
-
- ymid = (maxval(yEdgeNew(:)) - minval(yEdgeNew(:)))/2.0
- do i = 1,nEdgesNew
- fEdgeNew(i) = f0 + (yEdgeNew(i) - ymid) * beta
- enddo
-
-
- write(6,*) ' done not real bathymetry'
-endif ! if(.not.real_bathymetry) then
-
-
-if(real_bathymetry) then
-
-u_srcNew = 0.0
-do iEdge=1,nEdgesNew
- xin = xEdgeNew(iEdge)
- yin = yEdgeNew(iEdge)
- zin = zEdgeNew(iEdge)
- rlon = lonEdgeNew(iEdge)/dtr
- rlat = latEdgeNew(iEdge)/dtr
- ix = nint(rlon/0.1 - 0.05) + nlon + 1
- ix = mod(ix,nlon)+1
- iy = nlat
- do jcount=1,nlat
- if(t_lat(jcount).gt.rlat) then
- iy = jcount
- exit
- endif
- enddo
- !stress = -0.1*cos(3.0*latEdgeNew(iEdge))
- !ulon = stress
- !ulat = 0.0
- ulon = TAUX(ix,iy)
- ulat = TAUY(ix,iy)
- write(6,*) rlon, t_lon(ix), rlat, t_lat(iy)
- call transform_from_lonlat_to_xyz(xin,yin,zin,ulon,ulat,ux,uy,uz)
- if(boundaryEdgeNew(1,iEdge).eq.1) then
- u_srcNew(1,iEdge) = 0.0
- else
- iCell1 = cellsOnEdgeNew(1,iEdge)
- iCell2 = cellsOnEdgeNew(2,iEdge)
- p(1) = xCellNew(iCell1); p(2) = yCellNew(iCell1); p(3) = zCellNew(iCell1)
- q(1) = xCellNew(iCell2); q(2) = yCellNew(iCell2); q(3) = zCellNew(iCell2)
- q = q - p
- call unit_vector_in_3space(q)
- u_srcNew(1,iEdge) = ux*q(1) + uy*q(2) + uz*q(3)
- endif
-enddo
-
-
-!set tracers at a first guess
-temperatureNew = -99.0
-salinityNew = -99.0
-do iCell=1,nCellsNew
-do k = 1,maxLevelCellNew(iCell)
- temperatureNew(1,k,iCell) = 20.0 - 10.0*k/nVertLevelsMod
- salinityNew(1,k,iCell) = 34.0 ! salinity
-enddo
-enddo
-
-! update T and S field with WOCE data
-if(l_woce) then
-iNoData = 0
-do iCell=1,nCellsNew
- hNew(1,:,iCell) = dz(:)
- hZLevel = dz
- if(mod(iCell,100).eq.0) write(6,*) 'l_woce t and s',iCell
- rlon = lonCellNew(iCell)/dtr
- rlat = latCellNew(iCell)/dtr
- ix = nint(rlon/0.1 - 0.05) + nlon + 1
- ix = mod(ix,nlon)+1
- iy = nlat
- do j=1,nlat
- if(t_lat(j).gt.rlat) then
- iy = j
- exit
- endif
- enddo
- do k=1,maxLevelCellNew(iCell)
- ndata = 0; temp_t = 0; temp_s = 0; kdata(:) = 0
- do i=-15,15
- ixt = ix + 8*i
- if(ixt.lt.1) then
- ixt = ixt + nlon
- elseif(ixt.gt.nlon) then
- ixt = ixt - nlon
- endif
- do j=-15,15
- iyt = iy + 8*j
- flag_lat = .true.
- if(iyt.lt.1.or.iyt.gt.nlat) then
- iyt = 1
- flag_lat = .false.
- endif
- if(TEMP(ixt,iyt,k).gt.-10.0.and.flag_lat) then
- ndata = ndata + 1
- temp_t = temp_t + TEMP(ixt,iyt,k)
- temp_s = temp_s + SALT(ixt,iyt,k)
- endif
- enddo
- enddo
-
- if(ndata.gt.0) then
- temperatureNew(1,k,iCell) = temp_t / float(ndata)
- salinityNEW(1,k,iCell) = temp_s / float(ndata)
- kdata(k) = 1
- else
- if(k.eq.1) iNoData = iNoData + 1
- if(k.ge.3) then
- if(kdata(k-1).eq.1) maxLevelCellNew(iCell) = k-1
- endif
- endif
-
- enddo
-
-enddo
-
-! do a couple of smoothing passes
-do iter=1,5
-do iCell=1,nCellsNew
-do k=1,maxLevelCellNew(iCell)
- ndata=1
- temp_t = temperatureNew(1,k,iCell)
- temp_s = salinityNew(1,k,iCell)
- do j=1,nEdgesOnCellNew(iCell)
- jCell = cellsOnCellNew(j,iCell)
- if(jCell.gt.0) then
- if(maxLevelCellNew(jCell).ge.k) then
- temp_t = temp_t + temperatureNew(1,k,jCell)
- temp_s = temp_s + salinityNew(1,k,jCell)
- ndata = ndata + 1
- endif
- endif
- enddo
- temperatureNew(1,k,iCell) = temp_t / ndata
- salinityNew(1,k,iCell) = temp_s / ndata
-enddo
-enddo
-write(6,*) maxval(temperatureNew(1,1,:)),maxval(salinityNew(1,1,:))
-enddo
-
-write(6,*) iNoData, nCellsNew
-
-temperatureRestoreNew(:) = temperatureNew(1,1,:)
-salinityRestoreNew(:) = salinityNew(1,1,:)
-
-endif
-
-endif
-
-write(6,*) ' done get_init_conditions'
-
-deallocate(tempHolder)
-deallocate(vertCoord)
-
-end subroutine get_init_conditions
-
-
-subroutine error_checking
-real :: p(3), q(3), r(3), angle, s(3), t(3), dot, mindot, maxdot, b(vertexDegree)
-real :: work(nCellsNew)
-
-
-! write
-write(6,*)
-write(6,*) ' error checking '
-write(6,*)
-
-! check to see if every edge is normal to associated cells
-mindot = 2
-maxdot = -2
-do iEdge=1,nEdgesNew
- if(boundaryEdgeNew(1,iEdge).eq.1) cycle
- iCell1 = cellsOnEdgeNew(1,iEdge)
- iCell2 = cellsOnEdgeNew(2,iEdge)
- p(1)=xCellNew(iCell1); p(2)=yCellNew(iCell1); p(3)=zCellNew(iCell1)
- q(1)=xCellNew(iCell2); q(2)=yCellNew(iCell2); q(3)=zCellNew(iCell2)
- r(1)=xEdgeNew(iEdge); r(2)=yEdgeNew(iEdge); r(3)=zEdgeNew(iEdge)
- call unit_vector_in_3space(p)
- call unit_vector_in_3space(q)
- call unit_vector_in_3space(r)
- t = q - p
- s = r - p
- call unit_vector_in_3space(t)
- call unit_vector_in_3space(s)
- dot = s(1)*t(1)+s(2)*t(2)+s(3)*t(3)
- if(dot.lt.mindot) mindot=dot
- if(dot.gt.maxdot) maxdot=dot
-enddo
-write(6,10) 'alignment of edges and cells (should be ones)', mindot, maxdot
-10 format(a60,5x,2e15.5)
-
-! check to see if every segments connecting cells and vertices are orothogonal'
-mindot = 2
-maxdot = -2
-do iEdge=1,nEdgesNew
- if(boundaryEdgeNew(1,iEdge).eq.1) cycle
- iCell1 = cellsOnEdgeNew(1,iEdge)
- iCell2 = cellsOnEdgeNew(2,iEdge)
- iVertex1 = verticesOnEdgeNew(1,iEdge)
- iVertex2 = verticesOnEdgeNew(2,iEdge)
- p(1)=xCellNew(iCell1); p(2)=yCellNew(iCell1); p(3)=zCellNew(iCell1)
- q(1)=xCellNew(iCell2); q(2)=yCellNew(iCell2); q(3)=zCellNew(iCell2)
- r(1)=xVertexNew(iVertex1); r(2)=yVertexNew(iVertex1); r(3)=zVertexNew(iVertex1)
- s(1)=xVertexNew(iVertex2); s(2)=yVertexNew(iVertex2); s(3)=zVertexNew(iVertex2)
- call unit_vector_in_3space(p)
- call unit_vector_in_3space(q)
- call unit_vector_in_3space(r)
- call unit_vector_in_3space(s)
- t = q - p
- s = s - r
- call unit_vector_in_3space(t)
- call unit_vector_in_3space(s)
- dot = s(1)*t(1)+s(2)*t(2)+s(3)*t(3)
- if(dot.lt.mindot) mindot=dot
- if(dot.gt.maxdot) maxdot=dot
-enddo
-write(6,10) 'orthogonality of cell and vertex edges (should be zeros)', mindot, maxdot
-
-! check that the kiteareas sum to the areatriangle
-mindot = 2
-maxdot = -2
-do iVertex=1,nVerticesNew
- b = 0
- do i=1,vertexDegree
- b(i) = kiteAreasOnVertexNew(i,iVertex)
- enddo
- angle = sum(b)
- if(angle - areaTriangleNew(iVertex).lt.mindot) mindot = angle - areaTriangleNew(iVertex)
- if(angle - areaTriangleNew(iVertex).gt.maxdot) maxdot = angle - areaTriangleNew(iVertex)
-enddo
-write(6,10) ' error in sum of kites and triangles (should be zeroes)', mindot, maxdot
-
-! check that the kiteareas sum to the areaCell
-mindot = 2
-maxdot = -2
-work = 0
-do iVertex=1,nVerticesNew
- iCell1 = cellsOnVertexNew(1,iVertex)
- iCell2 = cellsOnVertexNew(2,iVertex)
- iCell3 = cellsOnVertexNew(3,iVertex)
- if(iCell1.ne.0) work(iCell1) = work(iCell1) + kiteAreasOnVertexNew(1,iVertex)
- if(iCell2.ne.0) work(iCell2) = work(iCell2) + kiteAreasOnVertexNew(2,iVertex)
- if(iCell3.ne.0) work(iCell3) = work(iCell3) + kiteAreasOnVertexNew(3,iVertex)
-enddo
-mindot = minval(areaCellNew - work)
-maxdot = maxval(areaCellNew - work)
-write(6,10) ' error in sum of kites and cells (should be zeroes)', mindot, maxdot
-
-!check for connectivity inverses for cells/edges
-do iCell=1,nCellsNew
- do i=1,nEdgesOnCellNew(iCell)
- iEdge=edgesOnCellNew(i,iCell)
- if(iEdge.le.0) stop ' iEdge le 0'
- iCell1 = cellsOnEdgeNew(1,iEdge)
- iCell2 = cellsOnEdgeNew(2,iEdge)
- if(iCell1.ne.iCell.and.iCell2.ne.iCell) stop ' cells/edges inverse failed'
- enddo
-enddo
-write(6,*) ' cellsOnEdge and edgesOnCell are duals for every cell/edge combination'
-
-!check for connectivity inverses for cells/vertices
-do iCell=1,nCellsNew
- do i=1,nEdgesOnCellNew(iCell)
- iVertex = verticesOnCellNew(i,iCell)
- if(iVertex.le.0) stop ' iVertex le 0'
- iCell1 = cellsOnVertexNew(1,iVertex)
- iCell2 = cellsOnVertexNew(2,iVertex)
- iCell3 = cellsOnVertexNew(3,iVertex)
- if(iCell1.ne.iCell.and.iCell2.ne.iCell.and.iCell3.ne.iCell) stop ' cells/vertices inverse failed'
- enddo
-enddo
-write(6,*) ' cellsOnVertex and verticesOnCell are duals for every cell/vertex combination'
-
-!check edgesOnEdge
-do iEdge=1,nEdgesNew
- iCell1 = cellsOnEdgeNew(1,iEdge)
- iCell2 = cellsOnEdgeNew(2,iEdge)
- if(nEdgesOnEdgeNew(iEdge).eq.0) then
- if(boundaryEdgeNew(1,iEdge).ne.1) stop ' stopping boundaryEdgeNew'
- endif
- do i=1,nEdgesOnEdgeNew(iEdge)
- jEdge = edgesOnEdgeNew(i,iEdge)
- jCell1 = cellsOnEdgeNew(1,jEdge)
- jCell2 = cellsOnEdgeNew(2,jEdge)
- if(jCell1.ne.iCell1.and.jCell1.ne.iCell2) then
- if(jCell2.ne.iCell1.and.jCell2.ne.iCell2) then
- write(6,*) 'error in edgesOnEdge'
- write(6,*) iCell1, iCell2, jCell1, jCell2
- stop
- endif
- endif
- enddo
-enddo
-write(6,*) ' edgesOnEdge is consistent with cellsOnEdge'
-
-end subroutine error_checking
-
-
-subroutine copy_dimensions
-
-maxEdgesNew = maxEdges
-maxEdges2New = maxEdges2
-TWONew = TWO
-vertexDegreeNew = vertexDegree
-nVertLevelsNew = nVertLevelsMod
-
-write(6,*)
-write(6,*) ' new dimensions '
-write(6,*) ' maxEdgesNew : ', maxEdgesNew
-write(6,*) ' maxEdges2New : ', maxEdges2New
-write(6,*) ' TWONew : ', TWONew
-write(6,*) ' vertexDegreeNew : ', vertexDegreeNew
-write(6,*) ' nVertLevelsNew : ', nVertLevelsNew
-
-end subroutine copy_dimensions
-
-
-
-subroutine read_grid
-implicit none
-
-call read_netcdf_init(nCells, nEdges, nVertices, maxEdges,maxEdges2,&
- nVertLevels,TWO,vertexDegree)
-
-write(6,*) ' init from grid '
-write(6,*) 'nCells :', nCells
-write(6,*) 'nEdges :', nEdges
-write(6,*) 'nVertices :', nVertices
-write(6,*) 'maxEdges :', maxEdges
-write(6,*) 'maxEdges2 :', maxEdges2
-write(6,*) 'nVertLevels :', nVertLevels
-write(6,*) 'vertexDegree :', vertexDegree
-write(6,*) 'TWO :', TWO
-
-allocate(xCell(nCells))
-allocate(yCell(nCells))
-allocate(zCell(nCells))
-allocate(latCell(nCells))
-allocate(lonCell(nCells))
-allocate(meshDensity(nCells))
-allocate(xEdge(nEdges))
-allocate(yEdge(nEdges))
-allocate(zEdge(nEdges))
-allocate(latEdge(nEdges))
-allocate(lonEdge(nEdges))
-allocate(xVertex(nVertices))
-allocate(yVertex(nVertices))
-allocate(zVertex(nVertices))
-allocate(latVertex(nVertices))
-allocate(lonVertex(nVertices))
-allocate(dcEdge(nEdges))
-allocate(dvEdge(nEdges))
-
-allocate(indexToCellID(nCells))
-allocate(indexToEdgeID(nEdges))
-allocate(indexToVertexID(nVertices))
-
-allocate(cellsOnEdge(TWO,nEdges))
-allocate(nEdgesOnCell(nCells))
-allocate(nEdgesOnEdge(nEdges))
-allocate(edgesOnCell(maxEdges,nCells))
-allocate(edgesOnEdge(maxEdges2,nEdges))
-allocate(weightsOnEdge(maxEdges2,nEdges))
-
-allocate(angleEdge(nEdges))
-allocate(areaCell(nCells))
-allocate(areaTriangle(nVertices))
-allocate(cellsOnCell(maxEdges,nCells))
-allocate(verticesOnCell(maxEdges,nCells))
-allocate(verticesOnEdge(TWO,nEdges))
-allocate(edgesOnVertex(vertexDegree,nVertices))
-allocate(cellsOnVertex(vertexDegree,nVertices))
-allocate(kiteAreasOnVertex(vertexDegree,nVertices))
-
-allocate(fEdge(nEdges))
-allocate(fVertex(nVertices))
-allocate(h_s(nCells))
-allocate(work1(nCells))
-allocate(u_src(nVertLevels,nEdges))
-allocate(u(1,nVertLevels,nEdges))
-allocate(v(1,nVertLevels,nEdges))
-allocate(h(1,nVertLevels,nCells))
-allocate(rho(1,nVertLevels,nCells))
-
-xCell=0; yCell=0; zCell=0; latCell=0; lonCell=0; meshDensity=1.0
-xEdge=0; yEdge=0; zEdge=0; latEdge=0; lonEdge=0
-xVertex=0; yVertex=0; zVertex=0; latVertex=0; lonVertex=0
-
-indexToCellID=0; indexToEdgeID=0; indexToVertexID=0
-cellsOnEdge=0; nEdgesOnCell=0; edgesOnCell=0
-edgesOnEdge=0; weightsOnEdge=0
-angleEdge=0; areaCell=0; areaTriangle=0
-cellsOnCell=0; verticesOnCell=0; verticesOnEdge=0
-edgesOnVertex=0; cellsOnVertex=0; kiteAreasOnVertex=0
-
-fEdge=0; fVertex=0; h_s=0; u_src=0; work1=0
-u=0; v=0; h=0; rho=0
-
-
-call read_netcdf_fields( &
- time, &
- latCell, &
- lonCell, &
- meshDensity, &
- xCell, &
- yCell, &
- zCell, &
- indexToCellID, &
- latEdge, &
- lonEdge, &
- xEdge, &
- yEdge, &
- zEdge, &
- indexToEdgeID, &
- latVertex, &
- lonVertex, &
- xVertex, &
- yVertex, &
- zVertex, &
- indexToVertexID, &
- cellsOnEdge, &
- nEdgesOnCell, &
- nEdgesOnEdge, &
- edgesOnCell, &
- edgesOnEdge, &
- weightsOnEdge, &
- dvEdge, &
- dcEdge, &
- angleEdge, &
- areaCell, &
- areaTriangle, &
- cellsOnCell, &
- verticesOnCell, &
- verticesOnEdge, &
- edgesOnVertex, &
- cellsOnVertex, &
- kiteAreasOnVertex, &
- fEdge, &
- fVertex, &
- h_s, &
- u, &
- v, &
- h &
- )
-
-write(6,*) ' values from read grid, min/max'
-write(6,*) ' latCell : ', minval(latCell), maxval(latCell)
-write(6,*) ' lonCell : ', minval(lonCell), maxval(lonCell)
-write(6,*) ' meshDensity : ', minval(meshDensity),maxval(meshDensity)
-write(6,*) ' xCell : ', minval(xCell), maxval(xCell)
-write(6,*) ' yCell : ', minval(yCell), maxval(yCell)
-write(6,*) ' zCell : ', minval(zCell), maxval(zCell)
-write(6,*) ' indexToCellID : ', minval(indexToCellID), maxval(indexToCellID)
-write(6,*) ' latEdge : ', minval(latEdge), maxval(latEdge)
-write(6,*) ' lonEdge : ', minval(lonEdge), maxval(lonEdge)
-write(6,*) ' xEdge : ', minval(xEdge), maxval(xEdge)
-write(6,*) ' yEdge : ', minval(yEdge), maxval(yEdge)
-write(6,*) ' zEdge : ', minval(zEdge), maxval(zEdge)
-write(6,*) ' indexToEdgeID : ', minval(indexToEdgeID), maxval(indexToEdgeID)
-write(6,*) ' latVertex : ', minval(latVertex), maxval(latVertex)
-write(6,*) ' lonVertex : ', minval(lonVertex), maxval(lonVertex)
-write(6,*) ' xVertex : ', minval(xVertex), maxval(xVertex)
-write(6,*) ' yVertex : ', minval(yVertex), maxval(yVertex)
-write(6,*) ' zVertex : ', minval(zVertex), maxval(zVertex)
-write(6,*) ' indexToVertexID : ', minval(indexToVertexID), maxval(indexToVertexID)
-write(6,*) ' cellsOnEdge : ', minval(cellsOnEdge), maxval(cellsOnEdge)
-write(6,*) ' nEdgesOnCell : ', minval(nEdgesOnCell), maxval(nEdgesOnCell)
-write(6,*) ' nEdgesOnEdge : ', minval(nEdgesOnEdge), maxval(nEdgesOnEdge)
-write(6,*) ' edgesOnCell : ', minval(edgesOnCell), maxval(edgesOnCell)
-write(6,*) ' edgesOnEdge : ', minval(edgesOnEdge), maxval(edgesOnEdge)
-write(6,*) ' weightsOnEdge : ', minval(weightsOnEdge), maxval(weightsOnEdge)
-write(6,*) ' dvEdge : ', minval(dvEdge), maxval(dvEdge)
-write(6,*) ' dcEdge : ', minval(dcEdge), maxval(dcEdge)
-write(6,*) ' angleEdge : ', minval(angleEdge), maxval(angleEdge)
-write(6,*) ' areaCell : ', minval(areaCell), maxval(areaCell)
-write(6,*) ' areaTriangle : ', minval(areaTriangle), maxval(areaTriangle)
-write(6,*) ' cellsOnCell : ', minval(cellsOnCell), maxval(cellsOnCell)
-write(6,*) ' verticesOnCell : ', minval(verticesOnCell), maxval(verticesOnCell)
-write(6,*) ' verticesOnEdge : ', minval(verticesOnEdge), maxval(verticesOnEdge)
-write(6,*) ' edgesOnVertex : ', minval(edgesOnVertex), maxval(edgesOnVertex)
-write(6,*) ' cellsOnVertex : ', minval(cellsOnVertex), maxval(cellsOnVertex)
-write(6,*) ' kiteAreasOnVertex : ', minval(kiteAreasOnVertex), maxval(kiteAreasOnVertex)
-write(6,*) ' fEdge : ', minval(fEdge), maxval(fEdge)
-write(6,*) ' fVertex : ', minval(fVertex), maxval(fVertex)
-write(6,*) ' h_s : ', minval(h_s), maxval(h_s)
-write(6,*) ' u : ', minval(u), maxval(u)
-write(6,*) ' v : ', minval(v), maxval(v)
-write(6,*) ' h : ', minval(h), maxval(h)
-
-end subroutine read_grid
-
-
-subroutine write_grid
-implicit none
-
-if(on_a_sphere.eq.'YES ') then
- xCellNew = xCellNew * sphere_radius
- yCellNew = yCellNew * sphere_radius
- zCellNew = zCellNew * sphere_radius
- xEdgeNew = xEdgeNew * sphere_radius
- yEdgeNew = yEdgeNew * sphere_radius
- zEdgeNew = zEdgeNew * sphere_radius
- xVertexNew = xVertexNew * sphere_radius
- yVertexNew = yVertexNew * sphere_radius
- zVertexNew = zVertexNew * sphere_radius
- dcEdgeNew = dcEdgeNew * sphere_radius
- dvEdgeNew = dvEdgeNew * sphere_radius
- areaCellNew = areaCellNew * (sphere_radius)**2
- areaTriangleNew = areaTriangleNew * (sphere_radius)**2
- kiteAreasOnVertexNew = kiteAreasOnVertexNew * (sphere_radius)**2
-endif
-
-call write_netcdf_init( &
- nCellsNew, &
- nEdgesNew, &
- nVerticesNew, &
- maxEdgesNew, &
- nVertLevelsNew, &
- vertexDegreeNew, &
- sphere_radius, &
- on_a_sphere &
- )
-
-call write_netcdf_fields( &
- 1, &
- latCellNew, &
- lonCellNew, &
- meshDensityNew, &
- xCellNew, &
- yCellNew, &
- zCellNew, &
- indexToCellIDNew, &
- latEdgeNew, &
- lonEdgeNew, &
- xEdgeNew, &
- yEdgeNew, &
- zEdgeNew, &
- indexToEdgeIDNew, &
- latVertexNew, &
- lonVertexNew, &
- xVertexNew, &
- yVertexNew, &
- zVertexNew, &
- indexToVertexIDNew, &
- maxLevelCellNew, &
- cellsOnEdgeNew, &
- nEdgesOnCellNew, &
- nEdgesOnEdgeNew, &
- edgesOnCellNew, &
- edgesOnEdgeNew, &
- weightsOnEdgeNew, &
- dvEdgeNew, &
- dcEdgeNew, &
- angleEdgeNew, &
- areaCellNew, &
- areaTriangleNew, &
- cellsOnCellNew, &
- verticesOnCellNew, &
- verticesOnEdgeNew, &
- edgesOnVertexNew, &
- cellsOnVertexNew, &
- kiteAreasOnVertexNew, &
- fEdgeNew, &
- fVertexNew, &
- h_sNew, &
- boundaryEdgeNew, &
- boundaryVertexNew, &
- u_srcNew, &
- uNew, &
- vNew, &
- hNew, &
- rhoNew, &
- temperatureNew, &
- salinityNew, &
- tracer1New, &
- temperatureRestoreNew, &
- salinityRestoreNew, &
- hZLevel &
- )
-
-call write_netcdf_finalize
-
-if(on_a_sphere.eq.'YES ') then
- xCellNew = xCellNew / sphere_radius
- yCellNew = yCellNew / sphere_radius
- zCellNew = zCellNew / sphere_radius
- xEdgeNew = xEdgeNew / sphere_radius
- yEdgeNew = yEdgeNew / sphere_radius
- zEdgeNew = zEdgeNew / sphere_radius
- xVertexNew = xVertexNew / sphere_radius
- yVertexNew = yVertexNew / sphere_radius
- zVertexNew = zVertexNew / sphere_radius
- dcEdgeNew = dcEdgeNew / sphere_radius
- dvEdgeNew = dvEdgeNew / sphere_radius
- areaCellNew = areaCellNew / (sphere_radius)**2
- areaTriangleNew = areaTriangleNew / (sphere_radius)**2
- kiteAreasOnVertexNew = kiteAreasOnVertexNew / (sphere_radius)**2
-endif
-
-end subroutine write_grid
-
-! Step 5: Check the depth routine define_kmt
-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
-logical :: flag, kmt_flag
-pi = 4.0*atan(1.0)
-dtr = pi / 180.0
-
-allocate(kmt(nCells))
-kmt = 0
-
-if(.not.real_bathymetry) then
- kmt = nVertLevelsMOD
- if(on_a_sphere.eq.'YES ') then
- write(6,*) 'Working on a sphere'
- latmin = -30*dtr
- latmax = +30*dtr
- lonmin = +10*dtr
- lonmax = +70*dtr
- write(6,*) ' lat min ', latmin
- write(6,*) ' lat max ', latmax
- where(latCell.lt.latmin) kmt = 0
- where(latCell.gt.latmax) kmt = 0
- where(lonCell.lt.lonmin) kmt = 0
- where(lonCell.gt.lonmax) kmt = 0
- else
- ! solid boundary in y
- ymin = minval(yCell)
- write(6,*) ' minimum yCell ', ymin
- ymax = maxval(yCell)
- write(6,*) ' maximum yCell ', ymax
- where(yCell.lt.1.001*ymin) kmt = 0
- where(yCell.gt.0.999*ymax) kmt = 0
-
-! xmin = minval(xCell)
-! xmax = maxval(xCell)
-! where(xCell.lt.1.001*xmin) kmt = 0
-! where(xCell.gt.0.999*xmax) kmt = 0
-
- ! ! solid boundary in x
- ! xmin = minval(xCell)
- ! write(6,*) ' minimum xCell ', xmin
- ! xmax = maxval(xCell)
- ! write(6,*) ' maximum xCell ', xmax
- ! where(xCell.lt.xmin+dc/1.5) kmt = 0
- ! where(xCell.gt.xmax-dc/1.5) kmt = 0
- endif
-
-
- allocate(work_kmt(nCells))
- work_kmt = 0.0
- 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
- nx = 10800
- ny = 5400
- allocate(x(nx))
- allocate(y(ny))
- allocate(ztopo(nx,ny))
- x = 0.0
- y = 0.0
- ztopo = 0.0
- write(6,*) ' ztopo ', minval(ztopo), maxval(ztopo)
- call read_topo_init( inx, iny)
- if(inx.ne.nx) stop ' nx topo'
- if(iny.ne.ny) stop ' ny topo'
- call read_topo_fields(x,y,ztopo)
- call read_topo_finalize()
- write(6,*) minval(x), maxval(x), x(1)
- write(6,*) minval(y), maxval(y), y(1)
- write(6,*) minval(ztopo), maxval(ztopo)
-
- do iCell=1,nCells
- rlon = lonCell(iCell) / dtr
- rlat = latCell(iCell) / dtr
- ix = nint((rlon+180)*30) + nx
- ix = mod(ix,nx)+1
- iy = nint((rlat+90 )*30)
- ix = max(1,ix); ix = min(nx,ix)
- iy = max(1,iy); iy = min(ny,iy)
-
- zdata = ztopo(ix,iy)
-
- if(zdata.lt.0.0) then
- zdata = -zdata
- r = 0
- kmt_flag=.false.
- do k=1,nVertLevelsMod
- if(.not.kmt_flag) then
- r = r + dz(k)
- if(r.gt.zdata) then
- kmt(iCell) = k
- kmt_flag = .true.
- endif
- endif
- enddo
- if(kmt(iCell).eq.0) kmt(iCell)=nVertLevelsMod
- ! write(6,*) kmt(iCell)
- endif
-
- ! if(zdata.lt.0.0) kmt(iCell) = nVertLevelsMod
-
- enddo
-
- deallocate(x)
- deallocate(y)
- deallocate(ztopo)
-endif
-
-! Eliminate isolated ocean cells, and make these isolated deep cells
-! flush with the deepest neighbor.
-do iCell=1,nCells
- kmt_neighbor_max = 0
- do j=1,nEdgesOnCell(iCell)
- iCell1 = cellsOnCell(j,iCell)
- kmt_neighbor_max = max(kmt_neighbor_max,kmt(iCell1))
- enddo
- kmt(iCell) = min(kmt(iCell),kmt_neighbor_max)
-enddo
-
-if(eliminate_inland_seas) then
-call eliminateLoops(nCells,nEdges,nVertices,maxEdges,vertexDegree, &
- nEdgesOnCell, cellsOnCell, verticesOnEdge, cellsOnVertex, edgesOnCell, lonCell, latCell, &
- xCell, yCell, zCell, xEdge, yEdge, zEdge, xVertex, yVertex, zVertex, &
- KMT)
-endif
-
-if(real_bathymetry) then
- where(kmt.eq.1) kmt=3
- where(kmt.eq.2) kmt=3
-endif
-
-end subroutine define_kmt
-
-
-
-subroutine define_mapping
-implicit none
-
-allocate(cellMap(nCells))
-allocate(edgeMap(nEdges))
-allocate(vertexMap(nVertices))
-cellMap = 0; edgeMap = 0; vertexMap = 0
-
-j=1
-do i=1,nCells
-if(kmt(i).ne.0) then
- cellMap(i) = j
- j=j+1
-endif
-write(10,*) i, cellMap(i)
-enddo
-
-j=1
-do i=1,nEdges
-iCell1 = cellsOnEdge(1,i)
-iCell2 = cellsOnEdge(2,i)
-if(kmt(iCell1).ne.0.or.kmt(iCell2).ne.0) then
- edgeMap(i)=j
- j=j+1
-endif
-write(11,*) i,edgeMap(i)
-enddo
-
-j=1
-do i=1,nVertices
-iCell1 = cellsOnVertex(1,i)
-iCell2 = cellsOnVertex(2,i)
-iCell3 = cellsOnVertex(3,i)
-if(kmt(iCell1).ne.0.or.kmt(iCell2).ne.0.or.kmt(iCell3).ne.0) then
- vertexMap(i)=j
- j=j+1
-endif
-write(12,*) i,vertexMap(i)
-enddo
-
-nCellsNew = 0
-do i=1,nCells
-if(cellMap(i).ne.0) nCellsNew = nCellsNew + 1
-enddo
-
-nEdgesNew = 0
-do i=1,nEdges
-if(edgeMap(i).ne.0) nEdgesNew = nEdgesNew + 1
-enddo
-
-nVerticesNew = 0
-do i=1,nVertices
-if(vertexMap(i).ne.0) nVerticesNew = nVerticesNew + 1
-enddo
-
-write(6,*) ' mesh mapping found '
-write(6,*) nCells, nCellsNew
-write(6,*) nEdges, nEdgesNew
-write(6,*) nVertices, nVerticesNew
-
-allocate(indexToCellIDNew(nCellsNew))
-allocate(indexToEdgeIDNew(nEdgesNew))
-allocate(indexToVertexIDNew(nVerticesNew))
-indextoCellIDNew = 0; indexToEdgeIDNew = 0; indexToVertexIDNew = 0
-
-do i=1,nCellsNew
-indexToCellIDNew(i)=i
-enddo
-
-do i=1,nEdgesNew
-indexToEdgeIDNew(i)=i
-enddo
-
-do i=1,nVerticesNew
-indexToVertexIDNew(i)=i
-enddo
-
-end subroutine define_mapping
-
-
-subroutine map_vectors
-implicit none
-
-allocate(xCellNew(nCellsNew))
-allocate(yCellNew(nCellsNew))
-allocate(zCellNew(nCellsNew))
-allocate(normalsNew(3,nEdgesNew))
-allocate(latCellNew(nCellsNew))
-allocate(lonCellNew(nCellsNew))
-allocate(meshDensityNew(nCellsNew))
-allocate(xEdgeNew(nEdgesNew))
-allocate(yEdgeNew(nEdgesNew))
-allocate(zEdgeNew(nEdgesNew))
-allocate(latEdgeNew(nEdgesNew))
-allocate(lonEdgeNew(nEdgesNew))
-allocate(xVertexNew(nVerticesNew))
-allocate(yVertexNew(nVerticesNew))
-allocate(zVertexNew(nVerticesNew))
-allocate(latVertexNew(nVerticesNew))
-allocate(lonVertexNew(nVerticesNew))
-allocate(dcEdgeNew(nEdgesNew))
-allocate(dvEdgeNew(nEdgesNew))
-allocate(angleEdgeNew(nEdgesNew))
-allocate(areaCellNew(nCellsNew))
-allocate(areaTriangleNew(nVerticesNew))
-allocate(maxLevelCellNew(nCellsNew))
-allocate(depthCell(nCellsNew))
-
-allocate(fEdgeNew(nEdgesNew))
-allocate(fVertexNew(nVerticesNew))
-allocate(h_sNew(nCellsNew))
-allocate(u_srcNew(nVertLevelsNew,nEdgesNew))
-allocate(uNew(1,nVertLevelsNew,nEdgesNew))
-allocate(vNew(1,nVertLevelsNew,nEdgesNew))
-allocate(hNew(1,nVertLevelsNew,nCellsNew))
-allocate(hZLevel(nVertLevelsNew))
-allocate(rhoNew(1,nVertLevelsNew,nCellsNew))
-allocate(temperatureNew(1,nVertLevelsNew,nCellsNew))
-allocate(salinityNew(1,nVertLevelsNew,nCellsNew))
-allocate(tracer1New(1,nVertLevelsNew,nCellsNew))
-
-allocate(temperatureRestoreNew(nCellsNew))
-allocate(salinityRestoreNew(nCellsNew))
-
-xCellNew=0; yCellNew=0; zCellNew=0; latCellNew=0; lonCellNew=0; meshDensityNew=1.0
-xEdgeNew=0; yEdgeNew=0; zEdgeNew=0; latEdgeNew=0; lonEdgeNew=0
-xVertexNew=0; yVertexNew=0; zVertexNew=0; latVertexNew=0; lonVertexNew=0
-
-fEdgeNew=0; fVertexNew=0; h_sNew=0; u_srcNew=0
-uNew=0; vNew=0; hNew=0; rhoNew=0
-temperatureNew=0; salinityNew=0; tracer1New=0;
-
-temperatureRestoreNew = 0.0
-salinityRestoreNew = 0.0
-
-
-do i=1,nCells
-jNew = cellMap(i)
-if(jNew.ne.0) then
- xCellNew(jNew)=xCell(i)
- yCellNew(jNew)=yCell(i)
- zCellNew(jNew)=zCell(i)
- latCellNew(jNew)=latCell(i)
- lonCellNew(jNew)=lonCell(i)
- meshDensityNew(jNew)=meshDensity(i)
- areaCellNew(jNew)=areaCell(i)
- maxLevelCellNew(jNew) = kmt(i)
- depthCell(jNew) = kmt(i)
-endif
-enddo
-
-do i=1,nEdges
-jNew = edgeMap(i)
-if(jNew.ne.0) then
- xEdgeNew(jNew)=xEdge(i)
- yEdgeNew(jNew)=yEdge(i)
- zEdgeNew(jNew)=zEdge(i)
- latEdgeNew(jNew)=latEdge(i)
- lonEdgeNew(jNew)=lonEdge(i)
- dcEdgeNew(jNew) = dcEdge(i)
- dvEdgeNew(jNew) = dvEdge(i)
- fEdgeNew(jNew) = fEdge(i)
- angleEdgeNew(jNew) = angleEdge(i)
-endif
-enddo
-
-do i=1,nVertices
-jNew = vertexMap(i)
-if(jNew.ne.0) then
- xVertexNew(jNew)=xVertex(i)
- yVertexNew(jNew)=yVertex(i)
- zVertexNew(jNew)=zVertex(i)
- latVertexNew(jNew)=latVertex(i)
- lonVertexNew(jNew)=lonVertex(i)
- fVertexNew(jNew)=fVertex(i)
- areaTriangleNew(jNew)=areaTriangle(i)
-endif
-enddo
-
-deallocate(xCell)
-deallocate(yCell)
-deallocate(zCell)
-deallocate(latCell)
-deallocate(lonCell)
-deallocate(meshDensity)
-deallocate(xEdge)
-deallocate(yEdge)
-deallocate(zEdge)
-deallocate(latEdge)
-deallocate(lonEdge)
-deallocate(xVertex)
-deallocate(yVertex)
-deallocate(zVertex)
-deallocate(latVertex)
-deallocate(lonVertex)
-deallocate(dcEdge)
-deallocate(dvEdge)
-
-end subroutine map_vectors
-
-
-
-subroutine map_connectivity
-implicit none
-
-allocate(cellsOnEdgeNew(TWONew,nEdgesNew))
-allocate(boundaryEdgeNew(nVertLevelsNew,nEdgesNew))
-allocate(flipVerticesOnEdgeOrdering(nEdgesNew))
-cellsOnEdgeNew(:,:) = 0
-boundaryEdgeNew(:,:) = 0
-flipVerticesOnEdgeOrdering(:) = 0
-do iEdge=1,nEdges
-if(edgeMap(iEdge).eq.0) cycle
-iEdgeNew = edgeMap(iEdge)
-iCell1 = cellsOnEdge(1,iEdge)
-iCell2 = cellsOnEdge(2,iEdge)
-iCell1New = cellMap(iCell1)
-iCell2New = cellMap(iCell2)
-cellsOnEdgeNew(1,iEdgeNew) = iCell1New
-cellsOnEdgeNew(2,iEdgeNew) = iCell2New
-if(iCell1New.eq.0.or.iCell2New.eq.0) boundaryEdgeNew(:,iEdgeNew) = 1
-if(iCell1New.eq.0.and.iCell2New.eq.0) stop "cellsOnEdge"
-if(iCell1New.eq.0) then
- cellsOnEdgeNew(1,iEdgeNew) = iCell2New
- cellsOnEdgeNew(2,iEdgeNew) = iCell1New
- flipVerticesOnEdgeOrdering(iEdgeNew) = 1
-endif
-enddo
-deallocate(cellsOnEdge)
-
-allocate(verticesOnEdgeNew(TWONew,nEdgesNew))
-allocate(boundaryVertexNew(nVertLevelsNew,nVerticesNew))
-verticesOnEdgeNew(:,:) = 0
-boundaryVertexNew(:,:) = 0
-do iEdge=1,nEdges
-if(edgeMap(iEdge).eq.0) cycle
-iEdgeNew = edgeMap(iEdge)
-iVertex1 = VerticesOnEdge(1,iEdge)
-iVertex2 = VerticesOnEdge(2,iEdge)
-iVertex1New = vertexMap(iVertex1)
-iVertex2New = vertexMap(iVertex2)
-if(iVertex1New.eq.0.or.iVertex2New.eq.0) stop "verticesOnEdge"
-if(flipVerticesOnEdgeOrdering(iEdgeNew).eq.0) then
- verticesOnEdgeNew(1,iEdgeNew) = iVertex1New
- verticesOnEdgeNew(2,iEdgeNew) = iVertex2New
-else
- verticesOnEdgeNew(1,iEdgeNew) = iVertex2New
- verticesOnEdgeNew(2,iEdgeNew) = iVertex1New
-endif
-if(boundaryEdgeNew(1,iEdgeNew).eq.1) then
- boundaryVertexNew(:,iVertex1New)=1
- boundaryVertexNew(:,iVertex2New)=1
-endif
-enddo
-deallocate(verticesOnEdge)
-
-allocate(nEdgesOnEdgeNew(nEdgesNew))
-allocate(edgesOnEdgeNew(maxEdges2,nEdgesNew))
-allocate(weightsOnEdgeNew(maxEdges2,nEdgesNew))
-nEdgesOnEdgeNew(:) = 0
-edgesOnEdgeNew(:,:) = 0
-weightsOnEdgeNew(:,:) = 0.0
-do iEdge=1,nEdges
-if(edgeMap(iEdge).eq.0) cycle
-iEdgeNew = edgeMap(iEdge)
-if(boundaryEdgeNew(1,iEdgeNew).eq.1) then
- nEdgesOnEdgeNew(iEdgeNew) = 0
- edgesOnEdgeNew(:,iEdgeNew) = 0
- weightsOnEdgeNew(:,iEdgeNew) = 0.0
-else
- nEdgesOnEdgeNew(iEdgeNew) = nEdgesOnEdge(iEdge)
- do i=1,nEdgesOnEdgeNew(iEdgeNew)
- jEdge = edgesOnEdge(i,iEdge)
- jEdgeNew = edgeMap(jEdge)
- if(jEdgeNew.eq.0) stop "jEdgeNew"
- edgesOnEdgeNew(i,iEdgeNew)=jEdgeNew
- weightsOnEdgeNew(i,iEdgeNew) = weightsOnEdge(i,iEdge)
- enddo
-endif
-enddo
-deallocate(nEdgesOnEdge)
-deallocate(edgesOnEdge)
-deallocate(weightsOnEdge)
-
-allocate(cellsOnCellNew(maxEdges,nCellsNew))
-allocate(nEdgesOnCellNew(nCellsNew))
-cellsOnCellNew = 0
-nEdgesOnCellNew = 0
-do iCell=1,nCells
-if(cellMap(iCell).eq.0) cycle
-iCellNew = cellMap(iCell)
-nEdgesOnCellNew(iCellNew)=nEdgesOnCell(iCell)
-do i=1,nEdgesOnCellNew(iCellNew)
-j = cellsOnCell(i,iCell)
-jNew = cellMap(j)
-cellsOnCellNew(i,iCellNew) = jNew
-enddo
-enddo
-deallocate(cellsOnCell)
-deallocate(nEdgesOnCell)
-
-allocate(edgesOnCellNew(maxEdgesNew,nCellsNew))
-edgesOnCellNew(:,:) = 0
-do iCell=1,nCells
-if(cellMap(iCell).eq.0) cycle
-iCellNew = cellMap(iCell)
-do i=1,nEdgesOnCellNew(iCellNew)
-j = edgesOnCell(i,iCell)
-jNew = edgeMap(j)
-if(jNew.eq.0) stop "edgesOnCell"
-edgesOnCellNew(i,iCellNew) = jNew
-enddo
-enddo
-deallocate(edgesOnCell)
-
-allocate(verticesOnCellNew(maxEdgesNew,nCellsNew))
-verticesOnCellNew(:,:)=0
-do iCell=1,nCells
-if(cellMap(iCell).eq.0) cycle
-iCellNew = cellMap(iCell)
-do i=1,nEdgesOnCellNew(iCellNew)
-j=verticesOnCell(i,iCell)
-jNew = vertexMap(j)
-if(jNew.eq.0) stop "verticesOnCell"
-verticesOnCellNew(i,iCellNew) = jNew
-enddo
-enddo
-deallocate(verticesOnCell)
-
-allocate(cellsOnVertexNew(vertexDegreeNew,nVerticesNew))
-allocate(kiteAreasOnVertexNew(vertexDegreeNew,nVerticesNew))
-cellsOnVertexNew = 0
-kiteAreasOnVertexNew = 0
-do iVertex=1,nVertices
-if(vertexMap(iVertex).eq.0) cycle
-iVertexNew = vertexMap(iVertex)
-do i=1,vertexDegree
-j=cellsOnVertex(i,iVertex)
-jNew=cellMap(j)
-if(jNew.eq.0) then
- kiteAreasOnVertexNew(i,iVertexNew)=0
-else
- kiteAreasOnVertexNew(i,iVertexNew)=kiteAreasOnVertex(i,iVertex)
-endif
-cellsOnVertexNew(i,iVertexNew)=jNew
-enddo
-enddo
-deallocate(cellsOnVertex)
-deallocate(kiteAreasOnVertex)
-
-areaTriangleNew = 0
-do iVertex=1,nVerticesNew
-do i=1,vertexDegree
-areaTriangleNew(iVertex) = areaTriangleNew(iVertex) + kiteAreasOnVertexNew(i,iVertex)
-enddo
-enddo
-
-allocate(edgesOnVertexNew(vertexDegreeNew, nVerticesNew))
-edgesOnVertexNew = 0
-do iVertex=1,nVertices
-if(vertexMap(iVertex).eq.0) cycle
-iVertexNew = vertexMap(iVertex)
-do i=1,vertexDegree
-j=edgesOnVertex(i,iVertex)
-jNew=edgeMap(j)
-edgesOnVertexNew(i,iVertexNew)=jNew
-enddo
-enddo
-deallocate(edgesOnVertex)
-
-! find normals
-normalsNew = 0.0
-do iEdge=1,nEdgesNew
-cell1 = cellsOnEdgeNew(1,iEdge)
-cell2 = cellsOnEdgeNew(2,iEdge)
-if(cell1.eq.0.or.cell2.eq.0) cycle
-c1(1) = xCellNew(cell1); c1(2) = yCellNew(cell1); c1(3) = zCellNew(cell1)
-c2(1) = xCellNew(cell2); c2(2) = yCellNew(cell2); c2(3) = zCellNew(cell2)
-distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
-
-if(on_a_sphere.eq.'YES ') then
- normalsNew(1,iEdge) = c2(1) - c1(1)
- normalsNew(2,iEdge) = c2(2) - c1(2)
- normalsNew(3,iEdge) = c2(3) - c1(3)
- distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
- normalsNew(:,iEdge) = normalsNew(:,iEdge) / distance
-else
- if(distance.gt.0.5*Lx) then
- write(6,*) ' periodic edge ', iEdge, distance
- write(6,10) ' c1 ', c1(:)
- write(6,10) ' c2 ', c2(:)
- r = c2(1) - c1(1)
- if(r.gt.0) c2(1) = c2(1) - Lx
- if(r.lt.0) c2(1) = c2(1) + Lx
- distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
- write(6,*) ' periodic edge fix ', iEdge, r, distance
- endif
- normalsNew(1,iEdge) = c2(1) - c1(1)
- normalsNew(2,iEdge) = c2(2) - c1(2)
- normalsNew(3,iEdge) = c2(3) - c1(3)
- distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
- normalsNew(:,iEdge) = normalsNew(:,iEdge) / distance
-endif
-enddo
-10 format(a20,3e15.5)
-
-end subroutine map_connectivity
-
-
-subroutine get_dz
-integer k
-
- dz( 1) = 1001.244 ! 5.006218 10.01244
- dz( 2) = 1011.258 ! 15.06873 20.12502
- dz( 3) = 1031.682 ! 25.28342 30.44183
- dz( 4) = 1063.330 ! 35.75848 41.07513
- dz( 5) = 1107.512 ! 46.61269 52.15025
- dz( 6) = 1166.145 ! 57.98098 63.81171
- dz( 7) = 1241.928 ! 70.02135 76.23099
- dz( 8) = 1338.612 ! 82.92405 89.61711
- dz( 9) = 1461.401 ! 96.92412 104.2311
- dz(10) = 1617.561 ! 112.3189 120.4067
- dz(11) = 1817.368 ! 129.4936 138.5804
- dz(12) = 2075.558 ! 148.9582 159.3360
- dz(13) = 2413.680 ! 171.4044 183.4728
- dz(14) = 2863.821 ! 197.7919 212.1110
- dz(15) = 3474.644 ! 229.4842 246.8575
- dz(16) = 4320.857 ! 268.4617 290.0660
- dz(17) = 5516.812 ! 317.6501 345.2342
- dz(18) = 7230.458 ! 381.3865 417.5388
- dz(19) = 9674.901 ! 465.9133 514.2878
- dz(20) = 13003.92 ! 579.3074 644.3270
- dz(21) = 17004.89 ! 729.3514 814.3759
- dz(22) = 20799.33 ! 918.3725 1022.369
- dz(23) = 23356.94 ! 1139.154 1255.939
- dz(24) = 24527.19 ! 1378.574 1501.210
- dz(25) = 24898.04 ! 1625.701 1750.191
- dz(26) = 24983.22 ! 1875.107 2000.023
- dz(27) = 24997.87 ! 2125.012 2250.002
- dz(28) = 24999.79 ! 2375.000 2500.000
- dz(29) = 24999.98 ! 2625.000 2749.999
- dz(30) = 25000.00 ! 2874.999 2999.999
- dz(31) = 25000.00 ! 3124.999 3249.999
- dz(32) = 25000.00 ! 3374.999 3499.999
- dz(33) = 25000.00 ! 3624.999 3749.999
- dz(34) = 25000.00 ! 3874.999 3999.999
- dz(35) = 25000.00 ! 4124.999 4249.999
- dz(36) = 25000.00 ! 4374.999 4499.999
- dz(37) = 25000.00 ! 4624.999 4749.999
- dz(38) = 25000.00 ! 4874.999 4999.999
- dz(39) = 25000.00 ! 5124.999 5249.999
- dz(40) = 25000.00 ! 5374.999 5499.999
-
- dz = dz / 100.0
-
- write(6,*)
- do k=1,40
- write(6,*) k,dz(k)
- enddo
- write(6,*)
-
-end subroutine get_dz
-end program map_to_basin
Added: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/get_init_conds.F
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/get_init_conds.F (rev 0)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/get_init_conds.F 2012-12-04 17:57:26 UTC (rev 2339)
@@ -0,0 +1,123 @@
+implicit none
+real :: halfwidth, dtr, pi, p(3), q(3), xin, yin, zin, ulon, ulat, stress, n1, n2, distance, r, temp_t, temp_s
+real :: dotProd
+real :: surfMaxTemp, surfMinTemp
+real :: y_0, x_0, x_1, x_2, x_3, width, cff1
+real :: yTrans, vertTempChange, maxTemp
+real, dimension(:), allocatable :: vertCoord
+real, dimension(:,:,:), allocatable :: tempHolder
+integer :: iTracer, ix, iy, ndata, i, j, k, ixt, iyt, ncull, jcount, iNoData, kdata(nVertLevelsMod)
+logical :: flag_lat
+
+pi = 4.0*atan(1.0)
+dtr = pi/180.0
+
+hNew = 100.0
+temperatureNew = 1.0
+salinityNew = 1.0
+tracer1New = 1.0
+uNew = 0
+
+allocate(vertCoord(nVertLevelsMOD))
+allocate(tempHolder(1, nVertLevelsMOD, nCellsNew))
+
+fEdgeNew(:) = 0.0
+fVertexNew(:) = 0.0
+bottomDepthNew(:) = 0.0
+uNew(:,:,:) = 0.0
+
+! basin-mod
+! setting for three levels - Set h values for isopycnal system
+write(6,*) ' setting three levels for isopycnal system'
+vertCoord = 0
+do i = 1, nVertLevelsMOD
+ hNew(1,i,:) = h_total_max / nVertLevelsMOD
+ hZLevel(i) = h_total_max / nVertLevelsMOD
+ vertCoord(i) = i * h_total_max / nVertLevelsMOD
+end do
+
+bottomDepthNew(:) = -h_total_max
+
+! basin-mod
+!Specify Density values for isopycnal levels
+write(6,*) ' setting density'
+rhoNew(1,:,:) = 1000.0
+if(nVertLevelsMOD .eq. 3) then
+ rhoNew(1,2,:) = 1011.0
+ rhoNew(1,3,:) = 1012.0
+endif
+
+! basin-mod
+! set temperature for isopycnal levels
+write(6,*) ' setting temperature'
+vertTempChange = 2.5
+surfMaxTemp = 13.1
+surfMinTemp = 10.1
+do k = 1, nVertLevelsMOD
+ temperatureNew(1,k,:) = (surfMaxTemp - surfMinTemp) * ((-vertCoord(k)+h_total_max)/h_total_max)+ surfMinTemp
+enddo
+
+tempHolder = temperatureNew
+
+y_0 = 200.0e3
+x_0 = 0.0e3
+x_1 = 160.0e3
+x_2 = 110.0e3
+x_3 = 130.0e3
+width = 40.0e3
+do i = 1, nCellsNew
+ cff1 = width * sin (6.0 * 3.141592 * (xCellNew(i) - x_0)/(x_1 - x_0))
+
+ if( yCellNew(i) < y_0 - cff1 ) then
+ do k = 1, nVertLevelsMOD
+ temperatureNew(1,k,i) = temperatureNew(1,k,i) - 1.2
+ end do
+ else if( yCellNew(i) .ge. y_0 - cff1 .and. yCellNew(i) .le. y_0 - cff1+width) then
+ do k = 1, nVertLevelsMOD
+ temperatureNew(1,k,i) = tempHolder(1,k,i) - 1.2*(1.0 -( yCellNew(i) - (y_0 - cff1)) / (1.0 * width))
+ end do
+ endif
+enddo
+
+do i = 1, nCellsNew
+ cff1 = 0.5 * width * sin(1.0 * 3.141592 * (xCellNew(i) - x_2)/(x_3 - x_2))
+
+ if( yCellNew(i) .ge. y_0 - cff1-0.5*width .and. yCellNew(i) .le. y_0 - cff1+0.5*width .and. xCellNew(i) .ge. x_2 .and. xCellNew(i) .le. x_3) then
+ do k = 1, nVertLevelsMOD
+ temperatureNew(1,k,i) = temperatureNew(1,k,i) + 0.3 * (1.0 - ( (yCellNew(i)-(y_0-cff1))/(0.5*width)))
+ end do
+ endif
+end do
+
+! basin-mod
+! set salinity for levels
+salinityNew(1,:,:) = 35.0
+
+! Updating density with linear EOS
+do i = 1,nCellsNew
+ rhoNew(1,:,i) = 1000.0*(1.0 - 2.5e-4*temperatureNew(1,1,i) + 7.6e-4*salinityNew(1,1,i))
+enddo
+
+! basin-mod
+! set forcing for isopycnal levels
+write(6,*) 'setting u_src - wind forcing'
+u_srcNew = 0.0
+write(6,*) ' u_srcNew ', minval(u_srcNew), maxval(u_srcNew)
+
+! basin-mod
+! set coriolis parameter for grid
+write(6,*) ' setting Coriolis parameter'
+ymid = (maxval(yVertexNew(:)) - minval(yVertexNew(:)))/2.0
+do i = 1,nVerticesNew
+ fVertexNew(i) = f0 + (yVertexNew(i) - ymid) * beta
+enddo
+
+ymid = (maxval(yEdgeNew(:)) - minval(yEdgeNew(:)))/2.0
+do i = 1,nEdgesNew
+ fEdgeNew(i) = f0 + (yEdgeNew(i) - ymid) * beta
+enddo
+
+write(6,*) ' done get_init_conditions'
+
+deallocate(tempHolder)
+deallocate(vertCoord)
Deleted: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_cullLoops.F
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_cullLoops.F 2012-12-03 21:05:55 UTC (rev 2338)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_cullLoops.F 2012-12-04 17:57:26 UTC (rev 2339)
@@ -1,293 +0,0 @@
-module cullLoops
-
- public :: eliminateLoops
-
- contains
-
- subroutine eliminateLoops(nCells,nEdges,nVertices,maxEdges,vertexDegree, &
- nEdgesOnCell, cellsOnCell, verticesOnEdge, cellsOnVertex, edgesOnCell, lonCell, latCell, &
- xCell, yCell, zCell, xEdge, yEdge, zEdge, xVertex, yVertex, zVertex, &
- KMT)
-
- implicit none
-
- ! intent (in)
- integer :: nCells, nEdges, nVertices, maxEdges, vertexDegree
- integer :: nEdgesOnCell(nCells), cellsOnCell(maxEdges,nCells), verticesOnEdge(2,nEdges)
- integer :: cellsOnVertex(vertexDegree,nVertices), edgesOnCell(maxEdges,nCells)
- real :: lonCell(nCells), latCell(nCells)
- real :: xCell(nCells), yCell(nCells), zCell(nCells)
- real :: xEdge(nEdges), yEdge(nEdges), zEdge(nEdges)
- real :: xVertex(nVertices), yVertex(nVertices), zVertex(nVertices)
- integer :: edgeList(nEdges), iCellMask(nCells)
-
- ! intent(inout)
- integer, intent(inout) :: KMT(ncells)
-
- ! local workspace
- integer :: iCell, jCell, oCell, lCell, iEdge, i, kCell, iSharedEdge, iStartEdge, iSave, iSweep
- integer :: iEdgeCounter, nEdgesInLoop(nCells), iCellAhead, LeftTurns, RightTurns
- logical :: connected, atBoundary, moveSouth, moveEast, atGrenwich
- real :: lat, rlat, rlon, rCenter(3), s(3), t(3), q(3), rCross, mylon, mylat, pi
-
- pi = 4.0*atan(1.0)
-
- ! we loop over all cells and count the number of edges in the loop containing iCell
- ! there is no coupling between iCell, so this can be inside an openMP directive
- iCellMask(:) = 0
- moveSouth = .true.
- do iCell=1,nCells
-
- write(6,*) 'working on : ',iCell, KMT(iCell)
-
- ! skip over land cells
- if(KMT(iCell).eq.0) then
- write(6,*) ' skipping : ', iCell
- cycle
- endif
-
- ! the working cell will be jCell, so set jCell=iCell to start
- jCell = iCell
- ! write(6,*) 'setting jCell: ', jCell
-
- atBoundary=.false. ! are we at a boundary?
- lCell = -1 ! when at a boundary, what is the index of the land cell or our right?
- oCell = -1 ! when at a boundary, what is the index of the ocean cell to our left?
-
- do while (.not.atBoundary)
-
- ! check to see if any edges of jCell are along the boundary
- do i=1,nEdgesOnCell(jCell)
- kCell = cellsOnCell(i,jCell)
- if(KMT(kCell).eq.0) then
- lCell = kCell ! this is a land cell
- oCell = jCell ! this is an ocean cell
- atBoundary = .true.
- write(6,*) ' found boundary : ',lCell,oCell
- endif
- enddo
-
- ! choose the next cell to be the one with the most southern latitude
- ! this jCell will only be used if atBoundary=.false., thus jCell must be an ocean cell
- if(moveSouth) then
- rlat = 10000.0
- mylat = latCell(jCell)
- do i=1,nEdgesOnCell(jCell)
- kCell = cellsOnCell(i,jCell)
- if(latCell(kCell).lt.rlat) then
- rlat = latCell(kCell)
- iSave = kCell
- endif
- enddo
- jCell = iSave
- endif
-
- ! if(moveSouth.and..not.atBoundary) write(6,*) ' pushing on to the south : ', jCell
-
- enddo ! .not.atBoundary
-
- ! OK, we hit a boundary ..... trace out the full loop in the CCW direction
- write(6,*) ' OK we hit a boundary, let us trace out this loop '
- write(6,*) ' ocean cell ', oCell, KMT(oCell)
- write(6,*) ' land cell ', lCell, KMT(lCell)
-
- ! start the counter at 1
- iEdgeCounter = 1
- edgeList(:) = 0
-
- ! find the shared edge where we are starting and save the starting edge index
- iSharedEdge = sharedEdge(oCell,lCell,nCells,maxEdges,nEdgesOnCell,edgesOnCell)
- iStartEdge = iSharedEdge
- edgeList(iEdgeCounter) = iSharedEdge
-
- connected = .false.
- LeftTurns = 0; RightTurns = 0
- do while (.not.connected)
-
- call moveAhead(xCell,yCell,zCell,xVertex,yVertex,zVertex, &
- oCell,lCell,iSharedEdge,nCells,nEdges,nVertices,maxEdges,vertexDegree, &
- verticesOnEdge,cellsOnVertex,iCellAhead)
-
- ! if the cell ahead is ocean, then the boundary is shared between lCell and iCellAhead
- ! if the cell ahead is land, then the boundary is shared between oCell and iCellAhead
- if(KMT(iCellAhead).gt.0) then
- oCell = iCellAhead
- RightTurns = RightTurns + 1
- ! write(6,*) ' the cell ahead is ocean, will turn right ', lCell, oCell
- else
- lCell = iCellAhead
- LeftTurns = LeftTurns + 1
- ! write(6,*) ' the cell ahead is land, will turn left ', lCell, oCell
- endif
- iSharedEdge = sharedEdge(oCell,lCell,nCells,maxEdges,nEdgesOnCell,edgesOnCell)
-
- ! check to see if we are where we started
- if(iSharedEdge.eq.iStartEdge) then
- connected=.true.
- write(6,*) ' we are back to the where we started '
- else
- iEdgeCounter=iEdgeCounter+1
- edgeList(iEdgeCounter) = iSharedEdge
- endif
-
- enddo ! .not.connected
-
- ! OK, we now have a loop .... but did we circle an inland see or a land mass?
- rCenter(:) = 0.0
- do iEdge=1,iEdgeCounter
- rCenter(1) = rCenter(1) + xEdge(edgeList(iEdge))/iEdgeCounter
- rCenter(2) = rCenter(2) + yEdge(edgeList(iEdge))/iEdgeCounter
- rCenter(3) = rCenter(3) + zEdge(edgeList(iEdge))/iEdgeCounter
- enddo
- rCenter(:) = rCenter(:) / sqrt ( rCenter(1)**2 + rCenter(2)**2 + rCenter(3)**2 )
-
- rCross = 0.0
- do iEdge=1,iEdgeCounter-1
- t(1) = xEdge(edgeList(iEdge+1)) - xEdge(edgeList(iEdge))
- t(2) = yEdge(edgeList(iEdge+1)) - yEdge(edgeList(iEdge))
- t(3) = zEdge(edgeList(iEdge+1)) - zEdge(edgeList(iEdge))
- s(1) = rCenter(1) - xEdge(edgeList(iEdge))
- s(2) = rCenter(2) - yEdge(edgeList(iEdge))
- s(3) = rCenter(3) - zEdge(edgeList(iEdge))
- t(:) = t(:) / sqrt( t(1)**2 + t(2)**2 + t(3)**2 )
- s(:) = s(:) / sqrt( s(1)**2 + s(2)**2 + s(3)**2 )
- call cross_product_in_R3(t,s,q)
- rCross = rCross + q(1)*rCenter(1) + q(2)*rCenter(2) + q(3)*rCenter(3)
- enddo
-
- write(6,*)
- write(6,*) ' edges and cross ', iEdgeCounter, rCross, LeftTurns, RightTurns
- write(6,*)
-
- if(LeftTurns-RightTurns.gt.0.and.rCross.gt.0.0) then
- iCellMask(iCell) = 1
- write(50,11) iCell, lonCell(iCell), latCell(iCell)
- 11 format(i8,2f12.4)
- endif
-
- enddo ! iCell
-
- ! cull all inland seas
- do iSweep=1,nCells/10
- write(6,*) iSweep
- do iCell=1,nCells
- if(iCellMask(iCell).eq.1) then
- do i=1,nEdgesOnCell(iCell)
- kCell=cellsOnCell(i,iCell)
- if(KMT(kCell).gt.0) iCellMask(kCell)=1
- enddo
- endif
- enddo
- enddo
-
- write(6,*) ' total cells culled ', sum(iCellMask)
- where(iCellMask(:).eq.1) KMT(:)=0
-
- end subroutine eliminateLoops
-
-
- subroutine moveAhead(xCell,yCell,zCell,xVertex,yVertex,zVertex, &
- oCell,lCell,iSharedEdge,nCells,nEdges,nVertices,maxEdges,vertexDegree, &
- verticesOnEdge,cellsOnVertex,iCellAhead)
- implicit none
- integer, intent(in) :: oCell, lCell,iSharedEdge, nCells, nEdges, nVertices, maxEdges, vertexDegree
- integer, intent(in) :: verticesOnEdge(2,nEdges), cellsOnVertex(vertexDegree,nVertices)
- real, intent(in), dimension(nCells) :: xCell, yCell, zCell
- real, intent(in), dimension(nVertices) :: xVertex, yVertex, zVertex
- integer, intent(out) :: iCellAhead
- integer :: iVertex1,iVertex2, i, kCell
- real :: v1(3), v2(3), ocean(3), land(3), d1, d2, cross1(3), cross2(3)
-
- ! solution assumes a CCW ordering of cellsOnVertex
- ! the cell ahead will be connected to the vertex that lists cellsOnVertex with lCell following oCell
-
- ! the vertex moving in the CCW direction has to be one of the two vertices connected to iSharedEdge
- iVertex1 = verticesOnEdge(1,iSharedEdge)
- iVertex2 = verticesOnEdge(2,iSharedEdge)
- !write(6,*) ' iVertex1, iVertex2 ', iVertex1, iVertex2
- !write(6,*) cellsOnVertex(:,iVertex1)
- !write(6,*) cellsOnVertex(:,iVertex2)
-
- v1(1)=xVertex(iVertex1)
- v1(2)=yVertex(iVertex1)
- v1(3)=zVertex(iVertex1)
-
- v2(1)=xVertex(iVertex2)
- v2(2)=yVertex(iVertex2)
- v2(3)=zVertex(iVertex2)
-
- ocean(1) = xCell(oCell)
- ocean(2) = yCell(oCell)
- ocean(3) = zCell(oCell)
-
- land(1) = xCell(lCell)
- land(2) = yCell(lCell)
- land(3) = zCell(lCell)
-
- ocean = land - ocean
- v1 = v1 - ocean
- v2 = v2 - ocean
-
- ocean(:) = ocean(:) / sqrt( ocean(1)**2 + ocean(2)**2 + ocean(3)**2)
- v1(:) = v1(:) / sqrt( v1(1)**2 + v1(2)**2 + v1(3)**2)
- v2(:) = v2(:) / sqrt( v2(1)**2 + v2(2)**2 + v2(3)**2)
-
- call cross_product_in_R3(ocean,v1,cross1)
- call cross_product_in_R3(ocean,v2,cross2)
-
- d1 = (ocean(1)+land(1))*cross1(1) + (ocean(2)+land(2))*cross1(2) + (ocean(3)+land(3))*cross1(3)
- d2 = (ocean(1)+land(1))*cross2(1) + (ocean(2)+land(2))*cross2(2) + (ocean(3)+land(3))*cross2(3)
-
- if(d1.gt.0.0) then
- do i=1,vertexDegree
- kCell=cellsOnVertex(i,iVertex1)
- if(kCell.ne.oCell.and.kCell.ne.lCell) iCellAhead=kCell
- enddo
- endif
-
- if(d2.gt.0.0) then
- do i=1,vertexDegree
- kCell=cellsOnVertex(i,iVertex2)
- if(kCell.ne.oCell.and.kCell.ne.lCell) iCellAhead=kCell
- enddo
- endif
-
- end subroutine moveAhead
-
-
-
- function sharedEdge(oCell,lCell,nCells,maxEdges,nEdgesOnCell,edgesOnCell)
- implicit none
- integer, intent(in) :: oCell, lCell, nCells, maxEdges, nEdgesOnCell(nCells), edgesOnCell(maxEdges,nCells)
- integer :: sharedEdge
- integer :: i,j,iEdge,jEdge
-
- sharedEdge=-1
- do i=1,nEdgesOnCell(oCell)
- iEdge = edgesOnCell(i,oCell)
- do j=1,nEdgesOnCell(lCell)
- jEdge = edgesOnCell(j,lCell)
- if(iEdge.eq.jEdge) then
- sharedEdge = jEdge
- exit
- endif
- enddo
- enddo
-
- if(SharedEdge.eq.-1) then
- write(6,*) ' problem with SharedEdge ',oCell,lCell
- stop
- endif
-
- end function sharedEdge
-
- subroutine cross_product_in_R3(p_1,p_2,p_out)
- real , intent(in) :: p_1 (3), p_2 (3)
- real , intent(out) :: p_out (3)
- p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
- p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
- p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
- end subroutine cross_product_in_R3
-
-
-end module cullLoops
Deleted: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_read_TS.F
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_read_TS.F 2012-12-03 21:05:55 UTC (rev 2338)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_read_TS.F 2012-12-04 17:57:26 UTC (rev 2339)
@@ -1,165 +0,0 @@
-module read_TS
-
- integer :: rd_ncid, rd_ncids, rd_ncidu
- integer :: rdDimIDt_lon
- integer :: rdDimIDt_lat
- integer :: rdDimIDdepth_t
- integer :: rdVarIDt_lon
- integer :: rdVarIDt_lat
- integer :: rdVarIDdepth_t
- integer :: rdVarIDTEMP
- integer :: rdVarIDSALT
- integer :: rdVarIDTAUX
- integer :: rdVarIDTAUY
-
- integer :: rdLocalt_lon
- integer :: rdLocalt_lat
- integer :: rdLocaldepth_t
-
- contains
-
- subroutine read_TS_init(nx, ny, nz)
-
- implicit none
-
- include 'netcdf.inc'
-
- integer, intent(out) :: nx, ny, nz
-
- integer :: nferr, nferrs, nferru
-
- nferr = nf_open('TS/woce_t_ann.3600x2431x42interp.r4.nc', NF_SHARE, rd_ncid)
- write(6,*) ' nferr ', nferr, rd_ncid
-
- !
- ! Get IDs for variable dimensions
- !
- nferr = nf_inq_dimid(rd_ncid, 't_lon', rdDimIDt_lon)
- write(6,*) ' nferr ', nferr, rdDimIDt_lon
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDt_lon, rdLocalt_lon)
- write(6,*) ' nferr ', nferr, rdLocalt_lon
- nferr = nf_inq_dimid(rd_ncid, 't_lat', rdDimIDt_lat)
- write(6,*) ' nferr ', nferr, rdDimIDt_lat
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDt_lat, rdLocalt_lat)
- write(6,*) ' nferr ', nferr, rdLocalt_lat
- nferr = nf_inq_dimid(rd_ncid, 'depth_t', rdDimIDdepth_t)
- write(6,*) ' nferr ', nferr, rdDimIDdepth_t
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDdepth_t, rdLocaldepth_t)
- write(6,*) ' nferr ', nferr, rdLocaldepth_t
-
- nx = rdLocalt_lon
- ny = rdLocalt_lat
- nz = rdLocaldepth_t
-
- write(6,*) nx, ny, nz
-
- !
- ! Get IDs for variables
- !
- nferr = nf_inq_varid(rd_ncid, 't_lon', rdVarIDt_lon)
- write(6,*) ' nferr ', nferr, rdVarIDt_lon
- nferr = nf_inq_varid(rd_ncid, 't_lat', rdVarIDt_lat)
- write(6,*) ' nferr ', nferr, rdVarIDt_lat
- nferr = nf_inq_varid(rd_ncid, 'depth_t', rdVarIDdepth_t)
- write(6,*) ' nferr ', nferr, rdVarIDdepth_t
- nferr = nf_inq_varid(rd_ncid, 'TEMP', rdVarIDTEMP)
- write(6,*) ' nferr ', nferr, rdVarIDTEMP
-
- nferrs = nf_open('TS/woce_s_ann.3600x2431x42interp.r4.nc', NF_SHARE, rd_ncids)
- nferrs = nf_inq_varid(rd_ncids, 'SALT', rdVarIDSALT)
- write(6,*) ' nferrs ', nferrs, rdVarIDSALT
-
- nferru = nf_open('TS/ws.old_ncep_1958-2000avg.interp3600x2431.nc', NF_SHARE, rd_ncidu)
- nferru = nf_inq_varid(rd_ncidu, 'TAUX', rdVarIDTAUX)
- nferru = nf_inq_varid(rd_ncidu, 'TAUY', rdVarIDTAUY)
- write(6,*) ' nferru ', nferru, rdVarIDTAUX, rdVarIDTAUY
-
- end subroutine read_TS_init
-
- subroutine read_TS_fields(t_lon, t_lat, depth_t, TEMP, SALT, TAUX, TAUY)
-
- implicit none
-
- include 'netcdf.inc'
-
- real (kind=4), dimension(:), intent(out) :: t_lon, t_lat, depth_t
- real (kind=4), dimension(:,:,:), intent(out) :: TEMP, SALT
- real (kind=4), dimension(:,:), intent(out) :: TAUX, TAUY
-
- integer, dimension(1) :: start1, count1
- integer, dimension(2) :: start2, count2
- integer, dimension(3) :: start3, count3
- integer, dimension(4) :: start4, count4
-
- integer :: nferr, nferrs, nferru
-
- start1(1) = 1
- count1(1) = rdLocalt_lon
- nferr = nf_get_vara_real(rd_ncid, rdVarIDt_lon, start1, count1, t_lon)
- write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDt_lon
-
- start1(1) = 1
- count1(1) = rdLocalt_lat
- nferr = nf_get_vara_real(rd_ncid, rdVarIDt_lat, start1, count1, t_lat)
- write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDt_lat
-
- start1(1) = 1
- count1(1) = rdLocaldepth_t
- nferr = nf_get_vara_real(rd_ncid, rdVarIDdepth_t, start1, count1, depth_t)
- write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDdepth_t
-
- start3(1) = 1
- start3(2) = 1
- start3(3) = 1
- count3(1) = rdLocalt_lon
- count3(2) = rdLocalt_lat
- count3(3) = rdLocaldepth_t
- nferr = nf_get_vara_real(rd_ncid, rdVarIDTEMP, start3, count3, TEMP)
- write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDTEMP
- write(6,*) ' temperature' , minval(TEMP), maxval(TEMP)
-
- start3(1) = 1
- start3(2) = 1
- start3(3) = 1
- count3(1) = rdLocalt_lon
- count3(2) = rdLocalt_lat
- count3(3) = rdLocaldepth_t
- nferrs = nf_get_vara_real(rd_ncids, rdVarIDSALT, start3, count3, SALT)
- write(6,*) ' nferrs ', nferrs, rd_ncids, rdVarIDSALT
- write(6,*) ' salinity' , minval(SALT), maxval(SALT)
-
- start2(1) = 1
- start2(2) = 1
- count2(1) = rdLocalt_lon
- count2(2) = rdLocalt_lat
- nferru = nf_get_vara_real(rd_ncidu, rdVarIDTAUX, start2, count2, TAUX)
- nferru = nf_get_vara_real(rd_ncidu, rdVarIDTAUY, start2, count2, TAUY)
- write(6,*) ' nferru ', nferru, rd_ncidu, rdVarIDTAUX, rdVarIDTAUY
- write(6,*) ' TAUX' , minval(TAUX), maxval(TAUX)
- write(6,*) ' TAUY' , minval(TAUY), maxval(TAUY)
-
-
- end subroutine read_TS_fields
-
-
- subroutine read_TS_finalize()
-
- implicit none
-
- include 'netcdf.inc'
-
- integer :: nferr, nferrs, nferru
-
- nferr = nf_close(rd_ncid)
- write(6,*) ' nferr ', nferr
-
- nferrs = nf_close(rd_ncids)
- write(6,*) ' nferrs ', nferrs
-
- nferru = nf_close(rd_ncidu)
- write(6,*) ' nferru ', nferru
-
-
- end subroutine read_TS_finalize
-
-end module read_TS
Deleted: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_read_netcdf.F
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_read_netcdf.F 2012-12-03 21:05:55 UTC (rev 2338)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_read_netcdf.F 2012-12-04 17:57:26 UTC (rev 2339)
@@ -1,523 +0,0 @@
-module read_netcdf
-
- integer :: rd_ncid
- integer :: rdDimIDTime
- integer :: rdDimIDnCells
- integer :: rdDimIDnEdges
- integer :: rdDimIDnVertices
- integer :: rdDimIDmaxEdges
- integer :: rdDimIDmaxEdges2
- integer :: rdDimIDnVertLevels
- integer :: rdDimIDTWO
- integer :: rdDimIDvertexDegree
- integer :: rdVarIDlatCell
- integer :: rdVarIDlonCell
- integer :: rdVarIDmeshDensity
- integer :: rdVarIDxCell
- integer :: rdVarIDyCell
- integer :: rdVarIDzCell
- integer :: rdVarIDindexToCellID
- integer :: rdVarIDlatEdge
- integer :: rdVarIDlonEdge
- integer :: rdVarIDxEdge
- integer :: rdVarIDyEdge
- integer :: rdVarIDzEdge
- integer :: rdVarIDindexToEdgeID
- integer :: rdVarIDlatVertex
- integer :: rdVarIDlonVertex
- integer :: rdVarIDxVertex
- integer :: rdVarIDyVertex
- integer :: rdVarIDzVertex
- integer :: rdVarIDindexToVertexID
- integer :: rdVarIDcellsOnEdge
- integer :: rdVarIDnEdgesOnCell
- integer :: rdVarIDnEdgesOnEdge
- integer :: rdVarIDedgesOnCell
- integer :: rdVarIDedgesOnEdge
- integer :: rdVarIDweightsOnEdge
- integer :: rdVarIDdvEdge
- integer :: rdVarIDdcEdge
- integer :: rdVarIDangleEdge
- integer :: rdVarIDareaCell
- integer :: rdVarIDareaTriangle
- integer :: rdVarIDcellsOnCell
- integer :: rdVarIDverticesOnCell
- integer :: rdVarIDverticesOnEdge
- integer :: rdVarIDedgesOnVertex
- integer :: rdVarIDcellsOnVertex
- integer :: rdVarIDkiteAreasOnVertex
- integer :: rdVarIDfEdge
- integer :: rdVarIDfVertex
- integer :: rdVarIDh_s
- integer :: rdVarIDu
- integer :: rdVarIDv
- integer :: rdVarIDh
-
- integer :: rdLocalnCells
- integer :: rdLocalnEdges
- integer :: rdLocalnVertices
- integer :: rdLocalmaxEdges
- integer :: rdLocalmaxEdges2
- integer :: rdLocalnVertLevels
- integer :: rdLocalTWO
- integer :: rdLocalvertexDegree
-
- contains
-
- subroutine read_netcdf_init( &
- nCells, &
- nEdges, &
- nVertices, &
- maxEdges, &
- maxEdges2, &
- nVertLevels, &
- TWO, &
- vertexDegree &
- )
-
- implicit none
-
- include 'netcdf.inc'
-
- integer, intent(out) :: nCells
- integer, intent(out) :: nEdges
- integer, intent(out) :: nVertices
- integer, intent(out) :: maxEdges
- integer, intent(out) :: maxEdges2
- integer, intent(out) :: nVertLevels
- integer, intent(out) :: TWO
- integer, intent(out) :: vertexDegree
-
- integer :: nferr
-
-
- nferr = nf_open('grid.nc', NF_SHARE, rd_ncid)
-
- !
- ! Get IDs for variable dimensions
- !
- nferr = nf_inq_unlimdim(rd_ncid, rdDimIDTime)
- nferr = nf_inq_dimid(rd_ncid, 'nCells', rdDimIDnCells)
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDnCells, rdLocalnCells)
- nferr = nf_inq_dimid(rd_ncid, 'nEdges', rdDimIDnEdges)
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDnEdges, rdLocalnEdges)
- nferr = nf_inq_dimid(rd_ncid, 'nVertices', rdDimIDnVertices)
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDnVertices, rdLocalnVertices)
- nferr = nf_inq_dimid(rd_ncid, 'maxEdges', rdDimIDmaxEdges)
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDmaxEdges, rdLocalmaxEdges)
- nferr = nf_inq_dimid(rd_ncid, 'maxEdges2', rdDimIDmaxEdges2)
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDmaxEdges2, rdLocalmaxEdges2)
- nferr = nf_inq_dimid(rd_ncid, 'nVertLevels', rdDimIDnVertLevels)
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDnVertLevels, rdLocalnVertLevels)
- nferr = nf_inq_dimid(rd_ncid, 'vertexDegree', rdDimIDvertexDegree)
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDvertexDegree, rdLocalvertexDegree)
- nferr = nf_inq_dimid(rd_ncid, 'TWO', rdDimIDTWO)
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDTWO, rdLocalTWO)
-
-
- nCells = rdLocalnCells
- nEdges = rdLocalnEdges
- nVertices = rdLocalnVertices
- maxEdges = rdLocalmaxEdges
- maxEdges2 = rdLocalmaxEdges2
- nVertLevels = rdLocalnVertLevels
- vertexDegree = rdLocalvertexDegree
- TWO = rdLocalTWO
-
- !
- ! Get IDs for variables
- !
- nferr = nf_inq_varid(rd_ncid, 'latCell', rdVarIDlatCell)
- nferr = nf_inq_varid(rd_ncid, 'lonCell', rdVarIDlonCell)
- nferr = nf_inq_varid(rd_ncid, 'meshDensity', rdVarIDmeshDensity)
- nferr = nf_inq_varid(rd_ncid, 'xCell', rdVarIDxCell)
- nferr = nf_inq_varid(rd_ncid, 'yCell', rdVarIDyCell)
- nferr = nf_inq_varid(rd_ncid, 'zCell', rdVarIDzCell)
- nferr = nf_inq_varid(rd_ncid, 'indexToCellID', rdVarIDindexToCellID)
- nferr = nf_inq_varid(rd_ncid, 'latEdge', rdVarIDlatEdge)
- nferr = nf_inq_varid(rd_ncid, 'lonEdge', rdVarIDlonEdge)
- nferr = nf_inq_varid(rd_ncid, 'xEdge', rdVarIDxEdge)
- nferr = nf_inq_varid(rd_ncid, 'yEdge', rdVarIDyEdge)
- nferr = nf_inq_varid(rd_ncid, 'zEdge', rdVarIDzEdge)
- nferr = nf_inq_varid(rd_ncid, 'indexToEdgeID', rdVarIDindexToEdgeID)
- nferr = nf_inq_varid(rd_ncid, 'latVertex', rdVarIDlatVertex)
- nferr = nf_inq_varid(rd_ncid, 'lonVertex', rdVarIDlonVertex)
- nferr = nf_inq_varid(rd_ncid, 'xVertex', rdVarIDxVertex)
- nferr = nf_inq_varid(rd_ncid, 'yVertex', rdVarIDyVertex)
- nferr = nf_inq_varid(rd_ncid, 'zVertex', rdVarIDzVertex)
- nferr = nf_inq_varid(rd_ncid, 'indexToVertexID', rdVarIDindexToVertexID)
- nferr = nf_inq_varid(rd_ncid, 'cellsOnEdge', rdVarIDcellsOnEdge)
- nferr = nf_inq_varid(rd_ncid, 'nEdgesOnCell', rdVarIDnEdgesOnCell)
- nferr = nf_inq_varid(rd_ncid, 'nEdgesOnEdge', rdVarIDnEdgesOnEdge)
- nferr = nf_inq_varid(rd_ncid, 'edgesOnCell', rdVarIDedgesOnCell)
- nferr = nf_inq_varid(rd_ncid, 'edgesOnEdge', rdVarIDedgesOnEdge)
- nferr = nf_inq_varid(rd_ncid, 'weightsOnEdge', rdVarIDweightsOnEdge)
- nferr = nf_inq_varid(rd_ncid, 'dvEdge', rdVarIDdvEdge)
- nferr = nf_inq_varid(rd_ncid, 'dcEdge', rdVarIDdcEdge)
- nferr = nf_inq_varid(rd_ncid, 'angleEdge', rdVarIDangleEdge)
- nferr = nf_inq_varid(rd_ncid, 'areaCell', rdVarIDareaCell)
- nferr = nf_inq_varid(rd_ncid, 'areaTriangle', rdVarIDareaTriangle)
- nferr = nf_inq_varid(rd_ncid, 'cellsOnCell', rdVarIDcellsOnCell)
- nferr = nf_inq_varid(rd_ncid, 'verticesOnCell', rdVarIDverticesOnCell)
- nferr = nf_inq_varid(rd_ncid, 'verticesOnEdge', rdVarIDverticesOnEdge)
- nferr = nf_inq_varid(rd_ncid, 'edgesOnVertex', rdVarIDedgesOnVertex)
- nferr = nf_inq_varid(rd_ncid, 'cellsOnVertex', rdVarIDcellsOnVertex)
- nferr = nf_inq_varid(rd_ncid, 'kiteAreasOnVertex', rdVarIDkiteAreasOnVertex)
- nferr = nf_inq_varid(rd_ncid, 'fEdge', rdVarIDfEdge)
- nferr = nf_inq_varid(rd_ncid, 'fVertex', rdVarIDfVertex)
- nferr = nf_inq_varid(rd_ncid, 'h_s', rdVarIDh_s)
- nferr = nf_inq_varid(rd_ncid, 'u', rdVarIDu)
- nferr = nf_inq_varid(rd_ncid, 'v', rdVarIDv)
- nferr = nf_inq_varid(rd_ncid, 'h', rdVarIDh)
-
- end subroutine read_netcdf_init
-
-
- subroutine read_netcdf_fields( &
- time, &
- latCell, &
- lonCell, &
- meshDensity, &
- xCell, &
- yCell, &
- zCell, &
- indexToCellID, &
- latEdge, &
- lonEdge, &
- xEdge, &
- yEdge, &
- zEdge, &
- indexToEdgeID, &
- latVertex, &
- lonVertex, &
- xVertex, &
- yVertex, &
- zVertex, &
- indexToVertexID, &
- cellsOnEdge, &
- nEdgesOnCell, &
- nEdgesOnEdge, &
- edgesOnCell, &
- edgesOnEdge, &
- weightsOnEdge, &
- dvEdge, &
- dcEdge, &
- angleEdge, &
- areaCell, &
- areaTriangle, &
- cellsOnCell, &
- verticesOnCell, &
- verticesOnEdge, &
- edgesOnVertex, &
- cellsOnVertex, &
- kiteAreasOnVertex, &
- fEdge, &
- fVertex, &
- h_s, &
- u, &
- v, &
- h &
- )
-
- implicit none
-
- include 'netcdf.inc'
-
- integer, intent(in) :: time
- real (kind=8), dimension(:), intent(out) :: latCell
- real (kind=8), dimension(:), intent(out) :: lonCell
- real (kind=8), dimension(:), intent(out) :: meshDensity
- real (kind=8), dimension(:), intent(out) :: xCell
- real (kind=8), dimension(:), intent(out) :: yCell
- real (kind=8), dimension(:), intent(out) :: zCell
- integer, dimension(:), intent(out) :: indexToCellID
- real (kind=8), dimension(:), intent(out) :: latEdge
- real (kind=8), dimension(:), intent(out) :: lonEdge
- real (kind=8), dimension(:), intent(out) :: xEdge
- real (kind=8), dimension(:), intent(out) :: yEdge
- real (kind=8), dimension(:), intent(out) :: zEdge
- integer, dimension(:), intent(out) :: indexToEdgeID
- real (kind=8), dimension(:), intent(out) :: latVertex
- real (kind=8), dimension(:), intent(out) :: lonVertex
- real (kind=8), dimension(:), intent(out) :: xVertex
- real (kind=8), dimension(:), intent(out) :: yVertex
- real (kind=8), dimension(:), intent(out) :: zVertex
- integer, dimension(:), intent(out) :: indexToVertexID
- integer, dimension(:,:), intent(out) :: cellsOnEdge
- integer, dimension(:), intent(out) :: nEdgesOnCell
- integer, dimension(:), intent(out) :: nEdgesOnEdge
- integer, dimension(:,:), intent(out) :: edgesOnCell
- integer, dimension(:,:), intent(out) :: edgesOnEdge
- real (kind=8), dimension(:,:), intent(out) :: weightsOnEdge
- real (kind=8), dimension(:), intent(out) :: dvEdge
- real (kind=8), dimension(:), intent(out) :: dcEdge
- real (kind=8), dimension(:), intent(out) :: angleEdge
- real (kind=8), dimension(:), intent(out) :: areaCell
- real (kind=8), dimension(:), intent(out) :: areaTriangle
- integer, dimension(:,:), intent(out) :: cellsOnCell
- integer, dimension(:,:), intent(out) :: verticesOnCell
- integer, dimension(:,:), intent(out) :: verticesOnEdge
- integer, dimension(:,:), intent(out) :: edgesOnVertex
- integer, dimension(:,:), intent(out) :: cellsOnVertex
- real (kind=8), dimension(:,:), intent(out) :: kiteAreasOnVertex
- real (kind=8), dimension(:), intent(out) :: fEdge
- real (kind=8), dimension(:), intent(out) :: fVertex
- real (kind=8), dimension(:), intent(out) :: h_s
- real (kind=8), dimension(:,:,:), intent(out) :: u
- real (kind=8), dimension(:,:,:), intent(out) :: v
- real (kind=8), dimension(:,:,:), intent(out) :: h
-
- logical :: meshDensityPresent
-
- integer :: nferr
- integer, dimension(1) :: start1, count1
- integer, dimension(2) :: start2, count2
- integer, dimension(3) :: start3, count3
- integer, dimension(4) :: start4, count4
-
- meshDensityPresent = .false.
-
- start1(1) = 1
-
- start2(1) = 1
- start2(2) = 1
-
- start3(1) = 1
- start3(2) = 1
- start3(3) = 1
-
- start1(1) = 1
- count1( 1) = rdLocalnCells
- count1( 1) = rdLocalnCells
- nferr = nf_get_vara_double(rd_ncid, rdVarIDlatCell, start1, count1, latCell)
-
- start1(1) = 1
- count1( 1) = rdLocalnCells
- count1( 1) = rdLocalnCells
- nferr = nf_get_vara_double(rd_ncid, rdVarIDlonCell, start1, count1, lonCell)
-
- start1(1) = 1
- count1( 1) = rdLocalnCells
- count1( 1) = rdLocalnCells
- nferr = nf_inq_varid(rd_ncid, 'meshDensity', rdVarIDmeshDensity)
- if(nferr.eq.0) then
- nferr = nf_get_vara_double(rd_ncid, rdVarIDmeshDensity, start1, count1, meshDensity)
- else
- meshDensity=1.0
- write(6,*) ' mesh density not present ', nferr, rdVarIDmeshDensity
- endif
-
- start1(1) = 1
- count1( 1) = rdLocalnCells
- count1( 1) = rdLocalnCells
- nferr = nf_get_vara_double(rd_ncid, rdVarIDxCell, start1, count1, xCell)
-
- start1(1) = 1
- count1( 1) = rdLocalnCells
- count1( 1) = rdLocalnCells
- nferr = nf_get_vara_double(rd_ncid, rdVarIDyCell, start1, count1, yCell)
-
- start1(1) = 1
- count1( 1) = rdLocalnCells
- count1( 1) = rdLocalnCells
- nferr = nf_get_vara_double(rd_ncid, rdVarIDzCell, start1, count1, zCell)
-
- start1(1) = 1
- count1( 1) = rdLocalnCells
- count1( 1) = rdLocalnCells
- nferr = nf_get_vara_int(rd_ncid, rdVarIDindexToCellID, start1, count1, indexToCellID)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_double(rd_ncid, rdVarIDlatEdge, start1, count1, latEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_double(rd_ncid, rdVarIDlonEdge, start1, count1, lonEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_double(rd_ncid, rdVarIDxEdge, start1, count1, xEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_double(rd_ncid, rdVarIDyEdge, start1, count1, yEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_double(rd_ncid, rdVarIDzEdge, start1, count1, zEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_int(rd_ncid, rdVarIDindexToEdgeID, start1, count1, indexToEdgeID)
-
- start1(1) = 1
- count1( 1) = rdLocalnVertices
- count1( 1) = rdLocalnVertices
- nferr = nf_get_vara_double(rd_ncid, rdVarIDlatVertex, start1, count1, latVertex)
-
- start1(1) = 1
- count1( 1) = rdLocalnVertices
- count1( 1) = rdLocalnVertices
- nferr = nf_get_vara_double(rd_ncid, rdVarIDlonVertex, start1, count1, lonVertex)
-
- start1(1) = 1
- count1( 1) = rdLocalnVertices
- count1( 1) = rdLocalnVertices
- nferr = nf_get_vara_double(rd_ncid, rdVarIDxVertex, start1, count1, xVertex)
-
- start1(1) = 1
- count1( 1) = rdLocalnVertices
- count1( 1) = rdLocalnVertices
- nferr = nf_get_vara_double(rd_ncid, rdVarIDyVertex, start1, count1, yVertex)
-
- start1(1) = 1
- count1( 1) = rdLocalnVertices
- count1( 1) = rdLocalnVertices
- nferr = nf_get_vara_double(rd_ncid, rdVarIDzVertex, start1, count1, zVertex)
-
- start1(1) = 1
- count1( 1) = rdLocalnVertices
- count1( 1) = rdLocalnVertices
- nferr = nf_get_vara_int(rd_ncid, rdVarIDindexToVertexID, start1, count1, indexToVertexID)
-
- start2(2) = 1
- count2( 1) = rdLocalTWO
- count2( 2) = rdLocalnEdges
- nferr = nf_get_vara_int(rd_ncid, rdVarIDcellsOnEdge, start2, count2, cellsOnEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnCells
- count1( 1) = rdLocalnCells
- nferr = nf_get_vara_int(rd_ncid, rdVarIDnEdgesOnCell, start1, count1, nEdgesOnCell)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_int(rd_ncid, rdVarIDnEdgesOnEdge, start1, count1, nEdgesOnEdge)
-
- start2(2) = 1
- count2( 1) = rdLocalmaxEdges
- count2( 2) = rdLocalnCells
- nferr = nf_get_vara_int(rd_ncid, rdVarIDedgesOnCell, start2, count2, edgesOnCell)
-
- start2(2) = 1
- count2( 1) = rdLocalmaxEdges2
- count2( 2) = rdLocalnEdges
- nferr = nf_get_vara_int(rd_ncid, rdVarIDedgesOnEdge, start2, count2, edgesOnEdge)
-
- start2(2) = 1
- count2( 1) = rdLocalmaxEdges2
- count2( 2) = rdLocalnEdges
- nferr = nf_get_vara_double(rd_ncid, rdVarIDweightsOnEdge, start2, count2, weightsOnEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_double(rd_ncid, rdVarIDdvEdge, start1, count1, dvEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_double(rd_ncid, rdVarIDdcEdge, start1, count1, dcEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_double(rd_ncid, rdVarIDangleEdge, start1, count1, angleEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnCells
- count1( 1) = rdLocalnCells
- nferr = nf_get_vara_double(rd_ncid, rdVarIDareaCell, start1, count1, areaCell)
-
- start1(1) = 1
- count1( 1) = rdLocalnVertices
- count1( 1) = rdLocalnVertices
- nferr = nf_get_vara_double(rd_ncid, rdVarIDareaTriangle, start1, count1, areaTriangle)
-
- start2(2) = 1
- count2( 1) = rdLocalmaxEdges
- count2( 2) = rdLocalnCells
- nferr = nf_get_vara_int(rd_ncid, rdVarIDcellsOnCell, start2, count2, cellsOnCell)
-
- start2(2) = 1
- count2( 1) = rdLocalmaxEdges
- count2( 2) = rdLocalnCells
- nferr = nf_get_vara_int(rd_ncid, rdVarIDverticesOnCell, start2, count2, verticesOnCell)
-
- start2(2) = 1
- count2( 1) = rdLocalTWO
- count2( 2) = rdLocalnEdges
- nferr = nf_get_vara_int(rd_ncid, rdVarIDverticesOnEdge, start2, count2, verticesOnEdge)
-
- start2(2) = 1
- count2( 1) = rdLocalvertexDegree
- count2( 2) = rdLocalnVertices
- nferr = nf_get_vara_int(rd_ncid, rdVarIDedgesOnVertex, start2, count2, edgesOnVertex)
-
- start2(2) = 1
- count2( 1) = rdLocalvertexDegree
- count2( 2) = rdLocalnVertices
- nferr = nf_get_vara_int(rd_ncid, rdVarIDcellsOnVertex, start2, count2, cellsOnVertex)
-
- start2(2) = 1
- count2( 1) = rdLocalvertexDegree
- count2( 2) = rdLocalnVertices
- nferr = nf_get_vara_double(rd_ncid, rdVarIDkiteAreasOnVertex, start2, count2, kiteAreasOnVertex)
-
- start1(1) = 1
- count1( 1) = rdLocalnEdges
- count1( 1) = rdLocalnEdges
- nferr = nf_get_vara_double(rd_ncid, rdVarIDfEdge, start1, count1, fEdge)
-
- start1(1) = 1
- count1( 1) = rdLocalnVertices
- count1( 1) = rdLocalnVertices
- nferr = nf_get_vara_double(rd_ncid, rdVarIDfVertex, start1, count1, fVertex)
-
- start1(1) = 1
- count1( 1) = rdLocalnCells
- count1( 1) = rdLocalnCells
- nferr = nf_get_vara_double(rd_ncid, rdVarIDh_s, start1, count1, h_s)
-
- start3(3) = time
- count3( 1) = rdLocalnVertLevels
- count3( 2) = rdLocalnEdges
- count3( 3) = 1
- nferr = nf_get_vara_double(rd_ncid, rdVarIDu, start3, count3, u)
-
- start3(3) = time
- count3( 1) = rdLocalnVertLevels
- count3( 2) = rdLocalnEdges
- count3( 3) = 1
- nferr = nf_get_vara_double(rd_ncid, rdVarIDv, start3, count3, v)
-
- start3(3) = time
- count3( 1) = rdLocalnVertLevels
- count3( 2) = rdLocalnCells
- count3( 3) = 1
- nferr = nf_get_vara_double(rd_ncid, rdVarIDh, start3, count3, h)
-
- end subroutine read_netcdf_fields
-
-
- subroutine read_netcdf_finalize()
-
- implicit none
-
- include 'netcdf.inc'
-
- integer :: nferr
-
- nferr = nf_close(rd_ncid)
-
- end subroutine read_netcdf_finalize
-
-end module read_netcdf
Deleted: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_read_topo.F
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_read_topo.F 2012-12-03 21:05:55 UTC (rev 2338)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_read_topo.F 2012-12-04 17:57:26 UTC (rev 2339)
@@ -1,109 +0,0 @@
-module read_topo
-
- integer :: rd_ncid
- integer :: rdDimIDnx
- integer :: rdDimIDny
- integer :: rdVarIDz
- integer :: rdVarIDx
- integer :: rdVarIDy
-
- integer :: rdLocalnx
- integer :: rdLocalny
-
- contains
-
- subroutine read_topo_init( nx, ny)
-
- implicit none
-
- include 'netcdf.inc'
-
- integer, intent(out) :: nx, ny
-
- integer :: nferr
-
-
- nferr = nf_open('topo/ETOPO2v2c_f4.nc', NF_SHARE, rd_ncid)
- write(6,*) ' nferr ', nferr, rd_ncid
-
- !
- ! Get IDs for variable dimensions
- !
- nferr = nf_inq_dimid(rd_ncid, 'x', rdDimIDnx)
- write(6,*) ' nferr ', nferr, rdDimIDnx
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDnx, rdLocalnx)
- write(6,*) ' nferr ', nferr, rdLocalnx
- nferr = nf_inq_dimid(rd_ncid, 'y', rdDimIDny)
- write(6,*) ' nferr ', nferr, rdDimIDny
- nferr = nf_inq_dimlen(rd_ncid, rdDimIDny, rdLocalny)
- write(6,*) ' nferr ', nferr, rdLocalny
-
- nx = rdLocalnx
- ny = rdLocalny
-
- write(6,*) nx, ny
-
- !
- ! Get IDs for variables
- !
- nferr = nf_inq_varid(rd_ncid, 'x', rdVarIDx)
- write(6,*) ' nferr ', nferr, rdVarIDx
- nferr = nf_inq_varid(rd_ncid, 'y', rdVarIDy)
- write(6,*) ' nferr ', nferr, rdVarIDy
- nferr = nf_inq_varid(rd_ncid, 'z', rdVarIDz)
- write(6,*) ' nferr ', nferr, rdVarIDz
-
- end subroutine read_topo_init
-
-
- subroutine read_topo_fields(x,y,z)
-
- implicit none
-
- include 'netcdf.inc'
-
- real (kind=4), dimension(:), intent(out) :: x,y
- real (kind=4), dimension(:,:), intent(out) :: z
-
- integer, dimension(1) :: start1, count1
- integer, dimension(2) :: start2, count2
- integer, dimension(3) :: start3, count3
- integer, dimension(4) :: start4, count4
-
- integer :: nferr
-
- start1(1) = 1
- count1(1) = rdLocalnx
- nferr = nf_get_vara_real(rd_ncid, rdVarIDx, start1, count1, x)
- write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDx
-
- start1(1) = 1
- count1(1) = rdLocalny
- nferr = nf_get_vara_real(rd_ncid, rdVarIDy, start1, count1, y)
- write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDy
-
- start2(1) = 1
- start2(2) = 1
- count2(1) = rdLocalnx
- count2(2) = rdLocalny
- nferr = nf_get_vara_real(rd_ncid, rdVarIDz, start2, count2, z)
- write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDz, rdLocalnx
-
- end subroutine read_topo_fields
-
-
- subroutine read_topo_finalize()
-
- implicit none
-
- include 'netcdf.inc'
-
- integer :: nferr
-
- nferr = nf_close(rd_ncid)
- write(6,*) ' nferr ', nferr
-
-
- end subroutine read_topo_finalize
-
-end module read_topo
Deleted: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_write_netcdf.F
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_write_netcdf.F 2012-12-03 21:05:55 UTC (rev 2338)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/module_write_netcdf.F 2012-12-04 17:57:26 UTC (rev 2339)
@@ -1,666 +0,0 @@
-module write_netcdf
-
- integer :: wr_ncid
- integer :: wrDimIDTime
- integer :: wrDimIDnCells
- integer :: wrDimIDnEdges
- integer :: wrDimIDnVertices
- integer :: wrDimIDmaxEdges
- integer :: wrDimIDmaxEdges2
- integer :: wrDimIDTWO
- integer :: wrDimIDvertexDegree
- integer :: wrDimIDnVertLevels
- integer :: wrVarIDlatCell
- integer :: wrVarIDlonCell
- integer :: wrVarIDmeshDensity
- integer :: wrVarIDxCell
- integer :: wrVarIDyCell
- integer :: wrVarIDzCell
- integer :: wrVarIDindexToCellID
- integer :: wrVarIDlatEdge
- integer :: wrVarIDlonEdge
- integer :: wrVarIDxEdge
- integer :: wrVarIDyEdge
- integer :: wrVarIDzEdge
- integer :: wrVarIDindexToEdgeID
- integer :: wrVarIDlatVertex
- integer :: wrVarIDlonVertex
- integer :: wrVarIDxVertex
- integer :: wrVarIDyVertex
- integer :: wrVarIDzVertex
- integer :: wrVarIDindexToVertexID
- integer :: wrVarIDmaxLevelCell
- integer :: wrVarIDcellsOnEdge
- integer :: wrVarIDnEdgesOnCell
- integer :: wrVarIDnEdgesOnEdge
- integer :: wrVarIDedgesOnCell
- integer :: wrVarIDedgesOnEdge
- integer :: wrVarIDweightsOnEdge
- integer :: wrVarIDdvEdge
- integer :: wrVarIDdcEdge
- integer :: wrVarIDangleEdge
- integer :: wrVarIDareaCell
- integer :: wrVarIDareaTriangle
- integer :: wrVarIDcellsOnCell
- integer :: wrVarIDverticesOnCell
- integer :: wrVarIDverticesOnEdge
- integer :: wrVarIDedgesOnVertex
- integer :: wrVarIDcellsOnVertex
- integer :: wrVarIDkiteAreasOnVertex
- integer :: wrVarIDfEdge
- integer :: wrVarIDfVertex
- integer :: wrVarIDh_s
- integer :: wrVarIDu
- integer :: wrVarIDboundaryEdge
- integer :: wrVarIDboundaryVertex
- integer :: wrVarIDu_src
- integer :: wrVarIDv
- integer :: wrVarIDh
- integer :: wrVarIDrho
- integer :: wrVarIDtemperature
- integer :: wrVarIDsalinity
- integer :: wrVarIDtracer1
- integer :: wrVarIDtemperatureRestore
- integer :: wrVarIDsalinityRestore
- integer :: wrVarIDhZLevel
-
- integer :: wrLocalnCells
- integer :: wrLocalnEdges
- integer :: wrLocalnVertices
- integer :: wrLocalmaxEdges
- integer :: wrLocalnVertLevels
- integer :: wrLocalvertexDegree
-
- contains
-
- subroutine write_netcdf_init( &
- nCells, &
- nEdges, &
- nVertices, &
- maxEdges, &
- nVertLevels, &
- vertexDegree, &
- sphere_radius, &
- on_a_sphere &
- )
-
- implicit none
-
- include 'netcdf.inc'
-
- integer, intent(in) :: nCells
- integer, intent(in) :: nEdges
- integer, intent(in) :: nVertices
- integer, intent(in) :: maxEdges
- integer, intent(in) :: nVertLevels
- integer, intent(in) :: vertexDegree
- character (len=16) :: on_a_sphere
- real*8 :: sphere_radius
-
-
- integer :: nferr
- integer, dimension(10) :: dimlist
-
-
- wrLocalnCells = nCells
- wrLocalnEdges = nEdges
- wrLocalnVertices = nVertices
- wrLocalmaxEdges = maxEdges
- wrLocalnVertLevels = nVertLevels
- wrLocalvertexDegree = vertexDegree
-
- nferr = nf_create('ocean.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), wr_ncid)
-
- !
- ! Define dimensions
- !
- nferr = nf_def_dim(wr_ncid, 'nCells', nCells, wrDimIDnCells)
- nferr = nf_def_dim(wr_ncid, 'nEdges', nEdges, wrDimIDnEdges)
- nferr = nf_def_dim(wr_ncid, 'nVertices', nVertices, wrDimIDnVertices)
- nferr = nf_def_dim(wr_ncid, 'maxEdges', maxEdges, wrDimIDmaxEdges)
- nferr = nf_def_dim(wr_ncid, 'maxEdges2', 2*maxEdges, wrDimIDmaxEdges2)
- nferr = nf_def_dim(wr_ncid, 'TWO', 2, wrDimIDTWO)
- nferr = nf_def_dim(wr_ncid, 'vertexDegree', vertexDegree, wrDimIDvertexDegree)
- nferr = nf_def_dim(wr_ncid, 'nVertLevels', nVertLevels, wrDimIDnVertLevels)
- nferr = nf_def_dim(wr_ncid, 'Time', NF_UNLIMITED, wrDimIDTime)
-
- !
- ! Define variables
- !
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'latCell', NF_DOUBLE, 1, dimlist, wrVarIDlatCell)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'lonCell', NF_DOUBLE, 1, dimlist, wrVarIDlonCell)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'meshDensity', NF_DOUBLE, 1, dimlist, wrVarIDmeshDensity)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'xCell', NF_DOUBLE, 1, dimlist, wrVarIDxCell)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'yCell', NF_DOUBLE, 1, dimlist, wrVarIDyCell)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'zCell', NF_DOUBLE, 1, dimlist, wrVarIDzCell)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'indexToCellID', NF_INT, 1, dimlist, wrVarIDindexToCellID)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'latEdge', NF_DOUBLE, 1, dimlist, wrVarIDlatEdge)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'lonEdge', NF_DOUBLE, 1, dimlist, wrVarIDlonEdge)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'xEdge', NF_DOUBLE, 1, dimlist, wrVarIDxEdge)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'yEdge', NF_DOUBLE, 1, dimlist, wrVarIDyEdge)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'zEdge', NF_DOUBLE, 1, dimlist, wrVarIDzEdge)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'indexToEdgeID', NF_INT, 1, dimlist, wrVarIDindexToEdgeID)
- dimlist( 1) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'latVertex', NF_DOUBLE, 1, dimlist, wrVarIDlatVertex)
- dimlist( 1) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'lonVertex', NF_DOUBLE, 1, dimlist, wrVarIDlonVertex)
- dimlist( 1) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'xVertex', NF_DOUBLE, 1, dimlist, wrVarIDxVertex)
- dimlist( 1) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'yVertex', NF_DOUBLE, 1, dimlist, wrVarIDyVertex)
- dimlist( 1) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'zVertex', NF_DOUBLE, 1, dimlist, wrVarIDzVertex)
- dimlist( 1) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'indexToVertexID', NF_INT, 1, dimlist, wrVarIDindexToVertexID)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'maxLevelCell', NF_INT, 1, dimlist, wrVarIDmaxLevelCell)
- dimlist( 1) = wrDimIDTWO
- dimlist( 2) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'cellsOnEdge', NF_INT, 2, dimlist, wrVarIDcellsOnEdge)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'nEdgesOnCell', NF_INT, 1, dimlist, wrVarIDnEdgesOnCell)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'nEdgesOnEdge', NF_INT, 1, dimlist, wrVarIDnEdgesOnEdge)
- dimlist( 1) = wrDimIDmaxEdges
- dimlist( 2) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'edgesOnCell', NF_INT, 2, dimlist, wrVarIDedgesOnCell)
- dimlist( 1) = wrDimIDmaxEdges2
- dimlist( 2) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'edgesOnEdge', NF_INT, 2, dimlist, wrVarIDedgesOnEdge)
- dimlist( 1) = wrDimIDmaxEdges2
- dimlist( 2) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'weightsOnEdge', NF_DOUBLE, 2, dimlist, wrVarIDweightsOnEdge)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'dvEdge', NF_DOUBLE, 1, dimlist, wrVarIDdvEdge)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'dcEdge', NF_DOUBLE, 1, dimlist, wrVarIDdcEdge)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'angleEdge', NF_DOUBLE, 1, dimlist, wrVarIDangleEdge)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'areaCell', NF_DOUBLE, 1, dimlist, wrVarIDareaCell)
- dimlist( 1) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'areaTriangle', NF_DOUBLE, 1, dimlist, wrVarIDareaTriangle)
- dimlist( 1) = wrDimIDmaxEdges
- dimlist( 2) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'cellsOnCell', NF_INT, 2, dimlist, wrVarIDcellsOnCell)
- dimlist( 1) = wrDimIDmaxEdges
- dimlist( 2) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'verticesOnCell', NF_INT, 2, dimlist, wrVarIDverticesOnCell)
- dimlist( 1) = wrDimIDTWO
- dimlist( 2) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'verticesOnEdge', NF_INT, 2, dimlist, wrVarIDverticesOnEdge)
- dimlist( 1) = wrDimIDvertexDegree
- dimlist( 2) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'edgesOnVertex', NF_INT, 2, dimlist, wrVarIDedgesOnVertex)
- dimlist( 1) = wrDimIDvertexDegree
- dimlist( 2) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'cellsOnVertex', NF_INT, 2, dimlist, wrVarIDcellsOnVertex)
- dimlist( 1) = wrDimIDvertexDegree
- dimlist( 2) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'kiteAreasOnVertex', NF_DOUBLE, 2, dimlist, wrVarIDkiteAreasOnVertex)
- dimlist( 1) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'fEdge', NF_DOUBLE, 1, dimlist, wrVarIDfEdge)
- dimlist( 1) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'fVertex', NF_DOUBLE, 1, dimlist, wrVarIDfVertex)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'h_s', NF_DOUBLE, 1, dimlist, wrVarIDh_s)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'temperatureRestore', NF_DOUBLE, 1, dimlist, wrVarIDtemperatureRestore)
- dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'salinityRestore', NF_DOUBLE, 1, dimlist, wrVarIDsalinityRestore)
- dimlist( 1) = wrDimIDnVertLevels
- nferr = nf_def_var(wr_ncid, 'hZLevel', NF_DOUBLE, 1, dimlist, wrVarIDhZLevel)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnEdges
- dimlist( 3) = wrDimIDTime
- nferr = nf_def_var(wr_ncid, 'u', NF_DOUBLE, 3, dimlist, wrVarIDu)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'boundaryEdge', NF_INT, 2, dimlist, wrVarIDboundaryEdge)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnVertices
- nferr = nf_def_var(wr_ncid, 'boundaryVertex', NF_INT, 2, dimlist, wrVarIDboundaryVertex)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnEdges
- nferr = nf_def_var(wr_ncid, 'u_src', NF_DOUBLE, 2, dimlist, wrVarIDu_src)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnEdges
- dimlist( 3) = wrDimIDTime
- nferr = nf_def_var(wr_ncid, 'v', NF_DOUBLE, 3, dimlist, wrVarIDv)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnCells
- dimlist( 3) = wrDimIDTime
- nferr = nf_def_var(wr_ncid, 'h', NF_DOUBLE, 3, dimlist, wrVarIDh)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnCells
- dimlist( 3) = wrDimIDTime
- nferr = nf_def_var(wr_ncid, 'rho', NF_DOUBLE, 3, dimlist, wrVarIDrho)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnCells
- dimlist( 3) = wrDimIDTime
- nferr = nf_def_var(wr_ncid, 'temperature', NF_DOUBLE, 3, dimlist, wrVarIDtemperature)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnCells
- dimlist( 3) = wrDimIDTime
- nferr = nf_def_var(wr_ncid, 'salinity', NF_DOUBLE, 3, dimlist, wrVarIDsalinity)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnCells
- dimlist( 3) = wrDimIDTime
- ! If you do not want tracer1 in your input file, simply comment out these two lines (one of two)
- nferr = nf_def_var(wr_ncid, 'tracer1', NF_DOUBLE, 3, dimlist, wrVarIDtracer1)
-
-
- nferr = nf_put_att_text(wr_ncid, NF_GLOBAL, 'on_a_sphere', 16, on_a_sphere)
- nferr = nf_put_att_double(wr_ncid, NF_GLOBAL, 'sphere_radius', NF_DOUBLE, 1, sphere_radius)
-
- nferr = nf_enddef(wr_ncid)
-
- end subroutine write_netcdf_init
-
-
- subroutine write_netcdf_fields( &
- time, &
- latCell, &
- lonCell, &
- meshDensity, &
- xCell, &
- yCell, &
- zCell, &
- indexToCellID, &
- latEdge, &
- lonEdge, &
- xEdge, &
- yEdge, &
- zEdge, &
- indexToEdgeID, &
- latVertex, &
- lonVertex, &
- xVertex, &
- yVertex, &
- zVertex, &
- indexToVertexID, &
- maxLevelCell, &
- cellsOnEdge, &
- nEdgesOnCell, &
- nEdgesOnEdge, &
- edgesOnCell, &
- edgesOnEdge, &
- weightsOnEdge, &
- dvEdge, &
- dcEdge, &
- angleEdge, &
- areaCell, &
- areaTriangle, &
- cellsOnCell, &
- verticesOnCell, &
- verticesOnEdge, &
- edgesOnVertex, &
- cellsOnVertex, &
- kiteAreasOnVertex, &
- fEdge, &
- fVertex, &
- h_s, &
- boundaryEdge, &
- boundaryVertex, &
- u_src, &
- u, &
- v, &
- h, &
- rho, &
- temperature, &
- salinity, &
- tracer1, &
- temperatureRestore, &
- salinityRestore, &
- hZLevel &
- )
-
- implicit none
-
- include 'netcdf.inc'
-
- integer, intent(in) :: time
- real (kind=8), dimension(:), intent(in) :: latCell
- real (kind=8), dimension(:), intent(in) :: lonCell
- real (kind=8), dimension(:), intent(in) :: meshDensity
- real (kind=8), dimension(:), intent(in) :: xCell
- real (kind=8), dimension(:), intent(in) :: yCell
- real (kind=8), dimension(:), intent(in) :: zCell
- integer, dimension(:), intent(in) :: indexToCellID
- real (kind=8), dimension(:), intent(in) :: latEdge
- real (kind=8), dimension(:), intent(in) :: lonEdge
- real (kind=8), dimension(:), intent(in) :: xEdge
- real (kind=8), dimension(:), intent(in) :: yEdge
- real (kind=8), dimension(:), intent(in) :: zEdge
- integer, dimension(:), intent(in) :: indexToEdgeID
- real (kind=8), dimension(:), intent(in) :: latVertex
- real (kind=8), dimension(:), intent(in) :: lonVertex
- real (kind=8), dimension(:), intent(in) :: xVertex
- real (kind=8), dimension(:), intent(in) :: yVertex
- real (kind=8), dimension(:), intent(in) :: zVertex
- integer, dimension(:), intent(in) :: indexToVertexID
- integer, dimension(:), intent(in) :: maxLevelCell
- integer, dimension(:,:), intent(in) :: cellsOnEdge
- integer, dimension(:), intent(in) :: nEdgesOnCell
- integer, dimension(:), intent(in) :: nEdgesOnEdge
- integer, dimension(:,:), intent(in) :: edgesOnCell
- integer, dimension(:,:), intent(in) :: edgesOnEdge
- real (kind=8), dimension(:,:), intent(in) :: weightsOnEdge
- real (kind=8), dimension(:), intent(in) :: dvEdge
- real (kind=8), dimension(:), intent(in) :: dcEdge
- real (kind=8), dimension(:), intent(in) :: angleEdge
- real (kind=8), dimension(:), intent(in) :: areaCell
- real (kind=8), dimension(:), intent(in) :: areaTriangle
- integer, dimension(:,:), intent(in) :: cellsOnCell
- integer, dimension(:,:), intent(in) :: verticesOnCell
- integer, dimension(:,:), intent(in) :: verticesOnEdge
- integer, dimension(:,:), intent(in) :: edgesOnVertex
- integer, dimension(:,:), intent(in) :: cellsOnVertex
- real (kind=8), dimension(:,:), intent(in) :: kiteAreasOnVertex
- real (kind=8), dimension(:), intent(in) :: fEdge
- real (kind=8), dimension(:), intent(in) :: fVertex
- real (kind=8), dimension(:), intent(in) :: h_s
- integer, dimension(:,:), intent(in) :: boundaryEdge
- integer, dimension(:,:), intent(in) :: boundaryVertex
- real (kind=8), dimension(:,:), intent(in) :: u_src
- real (kind=8), dimension(:,:,:), intent(in) :: u
- real (kind=8), dimension(:,:,:), intent(in) :: v
- real (kind=8), dimension(:,:,:), intent(in) :: h
- real (kind=8), dimension(:,:,:), intent(in) :: rho
- real (kind=8), dimension(:,:,:), intent(in) :: temperature
- real (kind=8), dimension(:,:,:), intent(in) :: salinity
- real (kind=8), dimension(:,:,:), intent(in) :: tracer1
- real (kind=8), dimension(:), intent(in) :: temperatureRestore
- real (kind=8), dimension(:), intent(in) :: salinityRestore
- real (kind=8), dimension(:), intent(in) :: hZLevel
-
-
- integer :: nferr
- integer, dimension(1) :: start1, count1
- integer, dimension(2) :: start2, count2
- integer, dimension(3) :: start3, count3
- integer, dimension(4) :: start4, count4
-
- start1(1) = 1
-
- start2(1) = 1
- start2(2) = 1
-
- start3(1) = 1
- start3(2) = 1
- start3(3) = 1
-
- start4(1) = 1
- start4(2) = 1
- start4(3) = 1
- start4(4) = 1
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDlatCell, start1, count1, latCell)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDlonCell, start1, count1, lonCell)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDmeshDensity, start1, count1, meshDensity)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDxCell, start1, count1, xCell)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDyCell, start1, count1, yCell)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDzCell, start1, count1, zCell)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToCellID, start1, count1, indexToCellID)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDlatEdge, start1, count1, latEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDlonEdge, start1, count1, lonEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDxEdge, start1, count1, xEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDyEdge, start1, count1, yEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDzEdge, start1, count1, zEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToEdgeID, start1, count1, indexToEdgeID)
-
- start1(1) = 1
- count1( 1) = wrLocalnVertices
- nferr = nf_put_vara_double(wr_ncid, wrVarIDlatVertex, start1, count1, latVertex)
-
- start1(1) = 1
- count1( 1) = wrLocalnVertices
- nferr = nf_put_vara_double(wr_ncid, wrVarIDlonVertex, start1, count1, lonVertex)
-
- start1(1) = 1
- count1( 1) = wrLocalnVertices
- nferr = nf_put_vara_double(wr_ncid, wrVarIDxVertex, start1, count1, xVertex)
-
- start1(1) = 1
- count1( 1) = wrLocalnVertices
- nferr = nf_put_vara_double(wr_ncid, wrVarIDyVertex, start1, count1, yVertex)
-
- start1(1) = 1
- count1( 1) = wrLocalnVertices
- nferr = nf_put_vara_double(wr_ncid, wrVarIDzVertex, start1, count1, zVertex)
-
- start1(1) = 1
- count1( 1) = wrLocalnVertices
- nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToVertexID, start1, count1, indexToVertexID)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_int(wr_ncid, wrVarIDmaxLevelCell, start1, count1, maxLevelCell)
-
- start2(2) = 1
- count2( 1) = 2
- count2( 2) = wrLocalnEdges
- nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnEdge, start2, count2, cellsOnEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_int(wr_ncid, wrVarIDnEdgesOnCell, start1, count1, nEdgesOnCell)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_int(wr_ncid, wrVarIDnEdgesOnEdge, start1, count1, nEdgesOnEdge)
-
- start2(2) = 1
- count2( 1) = wrLocalmaxEdges
- count2( 2) = wrLocalnCells
- nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnCell, start2, count2, edgesOnCell)
-
- start2(2) = 1
- count2( 1) = 2*wrLocalmaxEdges
- count2( 2) = wrLocalnEdges
- nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnEdge, start2, count2, edgesOnEdge)
-
- start2(2) = 1
- count2( 1) = 2*wrLocalmaxEdges
- count2( 2) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDweightsOnEdge, start2, count2, weightsOnEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDdvEdge, start1, count1, dvEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDdcEdge, start1, count1, dcEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDangleEdge, start1, count1, angleEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDareaCell, start1, count1, areaCell)
-
- start1(1) = 1
- count1( 1) = wrLocalnVertices
- nferr = nf_put_vara_double(wr_ncid, wrVarIDareaTriangle, start1, count1, areaTriangle)
-
- start2(2) = 1
- count2( 1) = wrLocalmaxEdges
- count2( 2) = wrLocalnCells
- nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnCell, start2, count2, cellsOnCell)
-
- start2(2) = 1
- count2( 1) = wrLocalmaxEdges
- count2( 2) = wrLocalnCells
- nferr = nf_put_vara_int(wr_ncid, wrVarIDverticesOnCell, start2, count2, verticesOnCell)
-
- start2(2) = 1
- count2( 1) = 2
- count2( 2) = wrLocalnEdges
- nferr = nf_put_vara_int(wr_ncid, wrVarIDverticesOnEdge, start2, count2, verticesOnEdge)
-
- start2(2) = 1
- count2( 1) = wrLocalvertexDegree
- count2( 2) = wrLocalnVertices
- nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnVertex, start2, count2, edgesOnVertex)
-
- start2(2) = 1
- count2( 1) = wrLocalvertexDegree
- count2( 2) = wrLocalnVertices
- nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnVertex, start2, count2, cellsOnVertex)
-
- start2(2) = 1
- count2( 1) = wrLocalvertexDegree
- count2( 2) = wrLocalnVertices
- nferr = nf_put_vara_double(wr_ncid, wrVarIDkiteAreasOnVertex, start2, count2, kiteAreasOnVertex)
-
- start1(1) = 1
- count1( 1) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDfEdge, start1, count1, fEdge)
-
- start1(1) = 1
- count1( 1) = wrLocalnVertices
- nferr = nf_put_vara_double(wr_ncid, wrVarIDfVertex, start1, count1, fVertex)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDh_s, start1, count1, h_s)
-
- start2(2) = 1
- count2( 1) = wrLocalnVertLevels
- count2( 2) = wrLocalnEdges
- nferr = nf_put_vara_int(wr_ncid, wrVarIDboundaryEdge, start2, count2, boundaryEdge)
-
- start2(2) = 1
- count2( 1) = wrLocalnVertLevels
- count2( 2) = wrLocalnVertices
- nferr = nf_put_vara_int(wr_ncid, wrVarIDboundaryVertex, start2, count2, boundaryVertex)
-
- start2(2) = 1
- count2( 1) = wrLocalnVertLevels
- count2( 2) = wrLocalnEdges
- nferr = nf_put_vara_double(wr_ncid, wrVarIDu_src, start2, count2, u_src)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDtemperatureRestore, start1, count1, temperatureRestore)
-
- start1(1) = 1
- count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDsalinityRestore, start1, count1, salinityRestore)
-
- start1(1) = 1
- count1( 1) = wrLocalnVertLevels
- nferr = nf_put_vara_double(wr_ncid, wrVarIDhZLevel, start1, count1, hZLevel)
-
- start3(3) = time
- count3( 1) = wrLocalnVertLevels
- count3( 2) = wrLocalnEdges
- count3( 3) = 1
- nferr = nf_put_vara_double(wr_ncid, wrVarIDu, start3, count3, u)
-
- start3(3) = time
- count3( 1) = wrLocalnVertLevels
- count3( 2) = wrLocalnEdges
- count3( 3) = 1
- nferr = nf_put_vara_double(wr_ncid, wrVarIDv, start3, count3, v)
-
- start3(3) = time
- count3( 1) = wrLocalnVertLevels
- count3( 2) = wrLocalnCells
- count3( 3) = 1
- nferr = nf_put_vara_double(wr_ncid, wrVarIDh, start3, count3, h)
-
- start3(3) = time
- count3( 1) = wrLocalnVertLevels
- count3( 2) = wrLocalnCells
- count3( 3) = 1
- nferr = nf_put_vara_double(wr_ncid, wrVarIDrho, start3, count3, rho)
-
- start3(3) = time
- count3( 1) = wrLocalnVertLevels
- count3( 2) = wrLocalnCells
- count3( 3) = 1
- nferr = nf_put_vara_double(wr_ncid, wrVarIDtemperature, start3, count3, temperature)
-
- start3(3) = time
- count3( 1) = wrLocalnVertLevels
- count3( 2) = wrLocalnCells
- count3( 3) = 1
- nferr = nf_put_vara_double(wr_ncid, wrVarIDsalinity, start3, count3, salinity)
-
- start3(3) = time
- count3( 1) = wrLocalnVertLevels
- count3( 2) = wrLocalnCells
- count3( 3) = 1
- ! If you do not want tracer1 in your input file, simply comment out these two lines (two of two)
- nferr = nf_put_vara_double(wr_ncid, wrVarIDtracer1, start3, count3, tracer1)
-
- end subroutine write_netcdf_fields
-
-
- subroutine write_netcdf_finalize()
-
- implicit none
-
- include 'netcdf.inc'
-
- integer :: nferr
-
- nferr = nf_close(wr_ncid)
-
- end subroutine write_netcdf_finalize
-
-end module write_netcdf
Deleted: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/utilities.F
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/utilities.F 2012-12-03 21:05:55 UTC (rev 2338)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/basin_src/utilities.F 2012-12-04 17:57:26 UTC (rev 2339)
@@ -1,776 +0,0 @@
-module utilities
-
-contains
-
-subroutine write_OpenDX( on_a_sphere, &
- nCells, &
- nVertices, &
- nEdges, &
- vertexDegree, &
- maxEdges, &
- xCell, &
- yCell, &
- zCell, &
- xVertex, &
- yVertex, &
- zVertex, &
- xEdge, &
- yEdge, &
- zEdge, &
- nEdgesOnCell, &
- verticesOnCell, &
- verticesOnEdge, &
- cellsOnVertex, &
- edgesOnCell, &
- areaCell, &
- maxLevelCell, &
- depthCell, &
- SST, &
- kiteAreasOnVertex )
-
- implicit none
-
- character (len=16), intent(in) :: on_a_sphere
- integer, intent(in) :: nCells, nVertices, vertexDegree, nEdges, maxEdges
- real (kind=8), dimension(nCells), intent(inout) :: xCell
- real (kind=8), dimension(nCells), intent(inout) :: yCell
- real (kind=8), dimension(nCells), intent(inout) :: zCell
- real (kind=8), dimension(nVertices), intent(inout) :: xVertex
- real (kind=8), dimension(nVertices), intent(inout) :: yVertex
- real (kind=8), dimension(nVertices), intent(inout) :: zVertex
- real (kind=8), dimension(nEdges), intent(inout) :: xEdge
- real (kind=8), dimension(nEdges), intent(inout) :: yEdge
- real (kind=8), dimension(nEdges), intent(inout) :: zEdge
- integer, dimension(nCells), intent(in) :: nEdgesOnCell
- integer, dimension(maxEdges,nCells), intent(in) :: verticesOnCell
- integer, dimension(maxEdges,nCells), intent(in) :: edgesOnCell
- integer, dimension(2,nEdges), intent(in) :: verticesOnEdge
- integer, dimension(vertexDegree, nVertices), intent(in) :: cellsOnVertex
- integer, dimension(nCells), intent(in) :: maxLevelCell
- real (kind=8), dimension(nCells), intent(in) :: areaCell
- real (kind=8), dimension(nCells), intent(in) :: depthCell, SST
- real (kind=8), dimension(vertexDegree,nVertices), intent(in) :: kiteAreasOnVertex
-
- character(len=80) :: a, b, c, d, e, f
- integer :: i, j, k, nVerticesTotal, iEdge, iLoop, iFace, Vert(4), Edge(4), iVertex, i1, i2, jp1
- integer :: nKitesTotal, iCell, iEdge1, iEdge2, iVertex11, iVertex12, iVertex21, iVertex22, ksave
- real (kind=8) :: x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xscale, work(nCells), work1(nCells)
- real (kind=8) :: xv, yv, zv, xc, yc, zc, dist
- logical (kind=8) :: eflag
-
- if(on_a_sphere.eq.'NO ') then
- write(6,*) ' write_dx, not on a sphere '
- endif
-
- xscale = 1.00
- xCell = xCell*xscale
- yCell = yCell*xscale
- zCell = zCell*xscale
- xVertex = xVertex*xscale
- yVertex = yVertex*xscale
- zVertex = zVertex*xscale
- xEdge = xEdge*xscale
- yEdge = yEdge*xscale
- zEdge = zEdge*xscale
-
- write(6,*) 'xCell', minval(xCell), maxval(xCell)
- write(6,*) ' nCells', nCells
- write(6,*) ' nEdges', nEdges
- write(6,*) ' nVertices', nVertices
- write(6,*) ' nEdgesOnCell',minval(nEdgesOnCell), maxval(nEdgesOnCell)
-
- open(unit=1,file='dx/vector.dx',form='formatted',status='unknown')
-
- a = trim('object "positions list" class array type float rank 1 shape 3 items')
- b = trim('ascii data file vector.position.data')
- write(1,10) a, nCells
- write(1,10) b
- write(1,*)
-
- a = trim('object 0 class array type float rank 1 shape 3 items')
- b = trim('ascii data file vector.data')
- c = trim('attribute "dep" string "positions"')
- write(1,10) a, nCells
- write(1,10) b
- write(1,10) c
- write(1,*)
-
- a = trim('object "vector" class field')
- b = trim('component "positions" "positions list"')
- c = trim('component "data" 0')
- write(1,10) a
- write(1,10) b
- write(1,10) c
-
- close(1)
-
- open(unit=14,file='dx/vector.position.data',form='formatted',status='unknown')
- do i=1,nCells
- write(14,22) xCell(i), yCell(i), zCell(i)
- enddo
- close(14)
-
-
-
- nVerticesTotal = 0
- do i=1,nCells
- nVerticesTotal = nVerticesTotal + nEdgesOnCell(i)
- enddo
- write(6,*) 'total number of vertices', nVerticesTotal
-
- open(unit=1,file='dx/ocean.dx',form='formatted',status='unknown')
-
- a = trim('object "positions list" class array type float rank 1 shape 3 items')
- b = trim('ascii data file ocean.position.data')
- write(1,10) a, nVerticesTotal
- write(1,10) b
- write(1,*)
- 10 format(a70,i10)
-
- a = trim('object "edge list" class array type int rank 0 items')
- b = trim('ascii data file ocean.edge.data')
- c = trim('attribute "ref" string "positions"')
- write(1,10) a, nVerticesTotal
- write(1,10) b
- write(1,10) c
- write(1,*)
-
- a = trim('object "loops list" class array type int rank 0 items')
- b = trim('ascii data file ocean.loop.data')
- c = trim('attribute "ref" string "edges"')
- write(1,10) a, nCells
- write(1,10) b
- write(1,10) c
- write(1,*)
-
- a = trim('object "face list" class array type int rank 0 items')
- b = trim('ascii data file ocean.face.data')
- c = trim('attribute "ref" string "loops"')
- write(1,10) a, nCells
- write(1,10) b
- write(1,10) c
- write(1,*)
-
- a = trim('object 0 class array type float rank 0 items')
- b = trim('data file ocean.area.data')
- c = trim('attribute "dep" string "faces"')
- write(1,10) a, nCells
- write(1,10) b
- write(1,10) c
- write(1,*)
-
- a = trim('object "area" class field')
- b = trim('component "positions" "positions list"')
- c = trim('component "edges" "edge list"')
- d = trim('component "loops" "loops list"')
- e = trim('component "faces" "face list"')
- f = trim('component "data" 0')
- write(1,10) a
- write(1,10) b
- write(1,10) c
- write(1,10) d
- write(1,10) e
- write(1,10) f
-
- close(1)
-
- work1 = depthCell
- work = SST
-
- open(unit= 9,file='dx/ocean.depth.data',form='formatted',status='unknown')
- open(unit=10,file='dx/ocean.area.data',form='formatted',status='unknown')
- open(unit=11,file='dx/ocean.face.data',form='formatted',status='unknown')
- open(unit=12,file='dx/ocean.loop.data',form='formatted',status='unknown')
- open(unit=13,file='dx/ocean.edge.data',form='formatted',status='unknown')
- open(unit=14,file='dx/ocean.position.data',form='formatted',status='unknown')
-
- iLoop = 0
- iEdge = 0
- do i=1,nCells
- write(9,20) work1(i)
- write(10,20) work(i)
- write(11,21) i-1
- write(12,21) iLoop
- iLoop = iLoop + nEdgesOnCell(i)
-
- eflag = .false.
- do j=1,nEdgesOnCell(i)
- k = verticesOnCell(j,i)
- xv = xVertex(k); yv = yVertex(k); zv = zVertex(k)
- xc = xCell(i); yc = yCell(i); zc = zCell(i)
- dist = sqrt( (xc-xv)**2 + (yc-yv)**2 + (zc-zv)**2 )
- if(dist.gt.5.0e5.and.on_a_sphere.eq.'NO ') then
- eflag = .true.
- endif
- enddo
-
- if(eflag) then
-
- do j=1,nEdgesOnCell(i)
- write(13,21) iEdge
- iEdge = iEdge + 1
- k = verticesOnCell(j,i)
- xv = xVertex(k); yv = yVertex(k); zv = zVertex(k)
- xc = xCell(i); yc = yCell(i); zc = zCell(i)
- dist = sqrt( (xc-xv)**2 + (yc-yv)**2 + (zc-zv)**2 )
- if(dist.gt.5.0e5) then
- write(14,22) xc, yc, zc
- else
- write(14,22) xv, yv, zv
- endif
- enddo
-
- else
-
- do j=1,nEdgesOnCell(i)
- write(13,21) iEdge
- iEdge = iEdge + 1
- k = verticesOnCell(j,i)
- if(k.le.0) write(6,*) ' vert1 ',k, verticesOnCell(:,i)
- write(14,22) xVertex(k), yVertex(k), zVertex(k)
- write(15,23) j,i,k,xVertex(k), yVertex(k), zVertex(k)
- enddo
- endif
- enddo
-
- 20 format(e20.10)
- 21 format(i20)
- 22 format(3e20.10)
- 23 format(3i8, 3e20.10)
-
- close(9)
- close(10)
- close(11)
- close(12)
- close(13)
- close(14)
-
- ! nVerticesTotal = 0
- ! nKitesTotal = 0
- ! do i=1,nCells
- ! nKitesTotal = nKitesTotal + nEdgesOnCell(i)
- ! enddo
- ! nVerticesTotal = nKitesTotal*4
- ! write(6,*) nKitesTotal, nVerticesTotal
-
- ! open(unit=1,file='dx/kite.dx',form='formatted',status='unknown')
-
- ! a = trim('object "positions list" class array type float rank 1 shape 3 items')
- ! b = trim('ascii data file kite.position.data')
- ! write(1,10) a, nVerticesTotal
- ! write(1,10) b
- ! write(1,*)
-
- ! a = trim('object "edge list" class array type int rank 0 items')
- ! b = trim('ascii data file kite.edge.data')
- ! c = trim('attribute "ref" string "positions"')
- ! write(1,10) a, nVerticesTotal
- ! write(1,10) b
- ! write(1,10) c
- ! write(1,*)
-
- ! a = trim('object "loops list" class array type int rank 0 items')
- ! b = trim('ascii data file kite.loop.data')
- ! c = trim('attribute "ref" string "edges"')
- ! write(1,10) a, nKitesTotal
- ! write(1,10) b
- ! write(1,10) c
- ! write(1,*)
-
- ! a = trim('object "face list" class array type int rank 0 items')
- ! b = trim('ascii data file kite.face.data')
- ! c = trim('attribute "ref" string "loops"')
- ! write(1,10) a, nKitesTotal
- ! write(1,10) b
- ! write(1,10) c
- ! write(1,*)
-
- ! a = trim('object 0 class array type float rank 0 items')
- ! b = trim('data file kite.area.data')
- ! c = trim('attribute "dep" string "faces"')
- ! write(1,10) a, nKitesTotal
- ! write(1,10) b
- ! write(1,10) c
- ! write(1,*)
-
- ! a = trim('object "area" class field')
- ! b = trim('component "positions" "positions list"')
- ! c = trim('component "edges" "edge list"')
- ! d = trim('component "loops" "loops list"')
- ! e = trim('component "faces" "face list"')
- ! f = trim('component "data" 0')
- ! write(1,10) a
- ! write(1,10) b
- ! write(1,10) c
- ! write(1,10) d
- ! write(1,10) e
- ! write(1,10) f
-
- ! close(1)
-
- ! open(unit=10,file='dx/kite.area.data',form='formatted',status='unknown')
- ! open(unit=11,file='dx/kite.face.data',form='formatted',status='unknown')
- ! open(unit=12,file='dx/kite.loop.data',form='formatted',status='unknown')
- ! open(unit=13,file='dx/kite.edge.data',form='formatted',status='unknown')
- ! open(unit=14,file='dx/kite.position.data',form='formatted',status='unknown')
-
- ! iLoop = 0
- ! iEdge = 0
- ! iFace = 0
-
- ! do iCell=1,nCells
- ! do j=1,nEdgesOnCell(iCell)
- ! iEdge1 = edgesOnCell(j,iCell)
- ! jp1 = j+1
- ! if(j.eq.nEdgesOnCell(iCell)) jp1=1
- ! iEdge2 = edgesOnCell(jp1,iCell)
-
- ! iVertex11 = verticesOnEdge(1,iEdge1)
- ! iVertex21 = verticesOnEdge(2,iEdge1)
- ! iVertex12 = verticesOnEdge(1,iEdge2)
- ! ivertex22 = verticesOnEdge(2,iEdge2)
-
- ! if(iVertex11.eq.iVertex12.or.iVertex11.eq.iVertex22) then
- ! iVertex = iVertex11
- ! elseif(iVertex21.eq.iVertex12.or.iVertex21.eq.iVertex22) then
- ! iVertex = iVertex21
- ! else
- ! write(6,*) iVertex11, iVertex21, iVertex12, iVertex22
- ! stop
- ! endif
-
- ! ksave = 0
- ! do k=1,vertexDegree
- ! if(cellsOnVertex(k,iVertex).eq.iCell) ksave=k
- ! enddo
- ! if(ksave.eq.0) then
- ! write(6,*) ' can not find iCell'
- ! write(6,*) cellsOnVertex(:,iVertex)
- ! write(6,*) iCell
- ! write(6,*) iEdge1, iEdge2
- ! write(6,*) iVertex11, iVertex21, iVertex21, iVertex22
- ! write(6,*) iVertex
- ! stop
- ! endif
-
- ! write(11,21) iFace
- ! write(12,21) iLoop
- ! iFace = iFace + 1
- ! iLoop = iLoop + 4
- ! do k=1,4
- ! write(13,21) iEdge
- ! iEdge = iEdge + 1
- ! enddo
- !
- ! x1 = xCell(iCell) ; y1 = yCell(iCell) ; z1 = zCell(iCell)
- ! x2 = xEdge(iEdge1) ; y2 = yEdge(iEdge1) ; z2 = zEdge(iEdge1)
- ! x3 = xVertex(iVertex); y3 = yVertex(iVertex); z3 = zVertex(iVertex)
- ! x4 = xEdge(iEdge2) ; y4 = yEdge(iEdge2) ; z4 = zEdge(iEdge2)
- !
- ! write(14,22) x1, y1, z1
- ! write(14,22) x2, y2, z2
- ! write(14,22) x3, y3, z3
- ! write(14,22) x4, y4, z4
- ! write(10,22) kiteAreasOnVertex(ksave,iVertex)
-
- ! enddo
- ! enddo
-
-end subroutine write_OpenDX
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! SUBROUTINE CONVERT_LX
-!
-! Convert (lat,lon) to an (x, y, z) location on a sphere with specified radius.
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-subroutine convert_lx(x, y, z, radius, lat, lon)
-
- implicit none
-
- real, intent(in) :: radius
- real, intent(in) :: lat, lon
- real, intent(out) :: x, y, z
-
- z = radius * sin(lat)
- x = radius * cos(lon) * cos(lat)
- y = radius * sin(lon) * cos(lat)
-
-end subroutine convert_lx
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! SUBROUTINE CONVERT_XL
-!
-! Convert (x, y, z) to a (lat, lon) location on a sphere with
-! radius sqrt(x^2 + y^2 + z^2).
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-subroutine convert_xl(x, y, z, lat,lon)
-
- implicit none
-
- real, intent(in) :: x, y, z
- real, intent(out) :: lat, lon
-
- real :: dl, clat, pii, rtod
- real :: eps
- parameter (eps=1.e-10)
-
- pii = 2.*asin(1.0)
- rtod=180./pii
- dl = sqrt(x*x + y*y + z*z)
-
- lat = asin(z/dl)
-
-! check for being close to either pole
-
- if (abs(x) > eps) then
-
- if (abs(y) > eps) then
-
- lon = atan(abs(y/x))
-
- if ((x <= 0.) .and. (y >= 0.)) then
- lon = pii-lon
- else if ((x <= 0.) .and. (y < 0.)) then
- lon = lon+pii
- else if ((x >= 0.) .and. (y <= 0.)) then
- lon = 2*pii-lon
- end if
-
- else ! we're either on longitude 0 or 180
-
- if (x > 0) then
- lon = 0.
- else
- lon = pii
- end if
-
- end if
-
- else if (abs(y) > eps) then
-
- if (y > 0) then
- lon = pii/2.
- else
- lon = 3.*pii/2.
- end if
-
- else ! we are at a pole
-
- lon = 0.
-
- end if
-
-end subroutine convert_xl
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-subroutine transform_from_lonlat_to_xyz(xin, yin, zin, ulon, ulat, ux, uy, uz)
-!
-! transform vector measured in latitude/longitude space to a vector measured in x,y,z
-!
-! INTENT(IN)
-! xin = x position
-! yin = y position
-! zin = z position
-! ulon = east component of vector
-! ulat = north component of vector
-!
-! INTENT(OUT)
-! ux = x component of vector
-! uy = y component of vector
-! uz = z component of vector
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-implicit none
-real, intent(in) :: xin, yin, zin, ulon, ulat
-real, intent(out) :: ux, uy, uz
-real :: h(3,3), p(3), q(3), g(3), X1(3,3), X2(3,3), trans_X2_to_X1(3,3), r
-integer :: i,j,k
-logical :: l_Pole
-real, parameter :: epsvt = 1.0e-10
-
-!-----------------------------------------------------------------------
-! define the e1, e2, and e3 directions
-!-----------------------------------------------------------------------
- X1(1,1) = 1.0; X1(1,2) = 0.0; X1(1,3) = 0.0
- X1(2,1) = 0.0; X1(2,2) = 1.0; X1(2,3) = 0.0
- X1(3,1) = 0.0; X1(3,2) = 0.0; X1(3,3) = 1.0
-
-!-----------------------------------------------------------------------
-! find the vectors (measured in X1) that point in the local
-! east (h(1,:)), north (h(2,:)), and vertical (h(3,:)) direction
-!-----------------------------------------------------------------------
- h(3,1) = xin; h(3,2) = yin; h(3,3) = zin
- call unit_vector_in_3space(h(3,:))
-
-!-----------------------------------------------------------------------
-! g(:) is a work array and holds the vector pointing to the North Pole.
-! measured in X1
-!-----------------------------------------------------------------------
- g(:) = X1(3,:)
-
-!-----------------------------------------------------------------------
-! determine if the local vertical hits a pole
-!-----------------------------------------------------------------------
- l_Pole = .false.
- r = g(1)*h(3,1) + g(2)*h(3,2) + g(3)*h(3,3)
- r = abs(r) + epsvt
- if(r.gt.1.0) then
- l_Pole = .true.
- h(3,:) = h(3,:) + epsvt
- call unit_vector_in_3space(h(3,:))
- endif
-
-!-----------------------------------------------------------------------
-! find the vector that is perpendicular to the local vertical vector
-! and points in the direction of of the North pole, this defines the local
-! north direction. measured in X1
-!-----------------------------------------------------------------------
- call vector_on_tangent_plane ( h(3,:), g(:), h(2,:) )
-
-!-----------------------------------------------------------------------
-! take the cross product of the local North direction and the local vertical
-! to find the local east vector. still in X1
-!-----------------------------------------------------------------------
- call cross_product_in_3space ( h(2,:), h(3,:), h(1,:) )
-
-!-----------------------------------------------------------------------
-! put these 3 vectors into a matrix X2
-!-----------------------------------------------------------------------
- X2(1,:) = h(1,:) ! local east (measured in X1)
- X2(2,:) = h(2,:) ! local north (measured in X1)
- X2(3,:) = h(3,:) ! local vertical (measured in X1)
-
-!-----------------------------------------------------------------------
-! compute the transformation matrix
-!-----------------------------------------------------------------------
- trans_X2_to_X1(:,:) = matmul(X1,transpose(X2))
-
-!-----------------------------------------------------------------------
-! transform (ulon, ulat) into (x,y,z)
-!-----------------------------------------------------------------------
- p(1) = ulon; p(2) = ulat; p(3) = 0
- g(:) = matmul(trans_X2_to_X1(:, :), p(:))
- ux = g(1); uy = g(2); uz = g(3)
-
-end subroutine transform_from_lonlat_to_xyz
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-subroutine transform_from_xyz_to_lonlat(xin, yin, zin, ux, uy, uz, ulon, ulat)
-!
-! transform vector measured in x,y,z space to a vector measured in latitude/longitude space
-!
-! INTENT(IN)
-! xin = x position
-! yin = y position
-! zin = z position
-! ux = x component of vector
-! uy = y component of vector
-! uz = z component of vector
-!
-! INTENT(OUT)
-! ulon = east component of vector
-! ulat = north component of vector
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-implicit none
-real, intent(in) :: xin, yin, zin, ux, uy, uz
-real, intent(out) :: ulon, ulat
-real :: h(3,3), p(3), q(3), g(3), X1(3,3), X2(3,3), trans_X1_to_X2(3,3), r
-integer :: i,j,k
-logical :: l_Pole
-real, parameter :: epsvt = 1.0e-10
-
-!-----------------------------------------------------------------------
-! define the e1, e2, and e3 directions
-!-----------------------------------------------------------------------
- X1(1,1) = 1.0; X1(1,2) = 0.0; X1(1,3) = 0.0
- X1(2,1) = 0.0; X1(2,2) = 1.0; X1(2,3) = 0.0
- X1(3,1) = 0.0; X1(3,2) = 0.0; X1(3,3) = 1.0
-
-!-----------------------------------------------------------------------
-! find the vectors (measured in X1) that point in the local
-! east (h(1,:)), north (h(2,:)), and vertical (h(3,:)) direction
-!-----------------------------------------------------------------------
- h(3,1) = xin; h(3,2) = yin; h(3,3) = zin
- call unit_vector_in_3space(h(3,:))
-
-!-----------------------------------------------------------------------
-! g(:) is a work array and holds the vector pointing to the North Pole.
-! measured in X1
-!-----------------------------------------------------------------------
- g(:) = X1(3,:)
-
-!-----------------------------------------------------------------------
-! determine if the local vertical hits a pole
-!-----------------------------------------------------------------------
- l_Pole = .false.
- r = g(1)*h(3,1) + g(2)*h(3,2) + g(3)*h(3,3)
- r = abs(r) + epsvt
- if(r.gt.1.0) then
- l_Pole = .true.
- h(3,:) = h(3,:) + epsvt
- call unit_vector_in_3space(h(3,:))
- endif
-
-!-----------------------------------------------------------------------
-! find the vector that is perpendicular to the local vertical vector
-! and points in the direction of of the North pole, this defines the local
-! north direction. measured in X1
-!-----------------------------------------------------------------------
- call vector_on_tangent_plane ( h(3,:), g(:), h(2,:) )
-
-!-----------------------------------------------------------------------
-! take the cross product of the local North direction and the local vertical
-! to find the local east vector. still in X1
-!-----------------------------------------------------------------------
- call cross_product_in_3space ( h(2,:), h(3,:), h(1,:) )
-
-!-----------------------------------------------------------------------
-! put these 3 vectors into a matrix X2
-!-----------------------------------------------------------------------
- X2(1,:) = h(1,:) ! local east (measured in X1)
- X2(2,:) = h(2,:) ! local north (measured in X1)
- X2(3,:) = h(3,:) ! local vertical (measured in X1)
-
-!-----------------------------------------------------------------------
-! compute the transformation matrix
-!-----------------------------------------------------------------------
- trans_X1_to_X2(:,:) = matmul(X2,transpose(X1))
-
-!-----------------------------------------------------------------------
-! transform (ulon, ulat) into (x,y,z)
-!-----------------------------------------------------------------------
- p(1) = ux; p(2) = uy; p(3) = uz
- g(:) = matmul(trans_X1_to_X2(:, :), p(:))
- ulon = g(1); ulat= g(2);
-
-end subroutine transform_from_xyz_to_lonlat
-
-!======================================================================
-! BEGINNING OF UNIT_VECTOR_IN_3SPACE
-!======================================================================
- subroutine unit_vector_in_3space (p_1)
-
-!-----------------------------------------------------------------------
-! PURPOSE : normalize p_1 to unit length and overwrite p_1
-!-----------------------------------------------------------------------
-
-!-----------------------------------------------------------------------
-! intent(inout)
-!-----------------------------------------------------------------------
- real , intent(inout) :: &
- p_1 (:)
-
-!-----------------------------------------------------------------------
-! local
-!-----------------------------------------------------------------------
- real :: length
-
- length = SQRT (p_1(1)**2 + p_1(2)**2 + p_1(3)**2 )
- length = 1.0/length
- p_1(1) = p_1(1)*length
- p_1(2) = p_1(2)*length
- p_1(3) = p_1(3)*length
-
- end subroutine unit_vector_in_3space
-!======================================================================
-! END OF UNIT_VECTOR_IN_3SPACE
-!======================================================================
-
-!======================================================================
-! BEGINNING OF CROSS_PRODUCT_IN_3SPACE
-!======================================================================
- subroutine cross_product_in_3space(p_1,p_2,p_out)
-
-!-----------------------------------------------------------------------
-! PURPOSE: compute p_1 cross p_2 and place in p_out
-!-----------------------------------------------------------------------
-
-!-----------------------------------------------------------------------
-! intent(in)
-!-----------------------------------------------------------------------
- real , intent(in) :: &
- p_1 (:), &
- p_2 (:)
-
-!-----------------------------------------------------------------------
-! intent(out)
-!-----------------------------------------------------------------------
- real , intent(out) :: &
- p_out (:)
-
- p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
- p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
- p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
-
- end subroutine cross_product_in_3space
-!======================================================================
-! END OF CROSS_PRODUCT_IN_3SPACE
-!======================================================================
-
-!======================================================================
-! BEGINNING OF VECTOR_ON_TANGENT_PLANE
-!======================================================================
- subroutine vector_on_tangent_plane(p_1, p_2, p_out)
-
-!-----------------------------------------------------------------------
-! PURPOSE : given two points measured in (x,y,z) and lying on
-! the unit sphere, find the vector (p_out) that lies on the plane
-! perpendicular to the p_1 vector and points in the direction of
-! the projection of p_2 onto the tangent plane.
-!
-! NOTE : p_1 and p_2 are assumed to be of unit length
-! NOTE : p_out is normalized to unit length
-!-----------------------------------------------------------------------
-
-!-----------------------------------------------------------------------
-! intent(in)
-!-----------------------------------------------------------------------
- real , intent(in) :: &
- p_1 (:), &
- p_2 (:)
-
-!-----------------------------------------------------------------------
-! intent(out)
-!-----------------------------------------------------------------------
- real , intent(out) :: &
- p_out (:)
-
-!-----------------------------------------------------------------------
-! local
-!-----------------------------------------------------------------------
- real :: &
- work (3), t1(3), t2(3)
-
-! work (1) = - p_1(2) * ( -p_1(2) * p_2(1) + p_1(1) * p_2(2) ) &
-! + p_1(3) * ( p_1(3) * p_2(1) - p_1(1) * p_2(3) )
-
-! work (2) = + p_1(1) * ( -p_1(2) * p_2(1) + p_1(1) * p_2(2) ) &
-! - p_1(3) * ( -p_1(3) * p_2(2) + p_1(2) * p_2(3) )
-
-! work (3) = - p_1(1) * ( p_1(3) * p_2(1) - p_1(1) * p_2(3) ) &
-! + p_1(2) * ( -p_1(3) * p_2(2) + p_1(2) * p_2(3) )
-
-
- t1(:) = p_2(:) - p_1(:)
- t2(:) = p_1
-
- call unit_vector_in_3space (t1)
- call unit_vector_in_3space (t2)
-
- call cross_product_in_3space(t1(:), t2(:), work(:))
- call unit_vector_in_3space (work)
- call cross_product_in_3space(t2(:),work(:),p_out(:))
- call unit_vector_in_3space (p_out)
-
- end subroutine vector_on_tangent_plane
-!======================================================================
-! END OF VECTOR_ON_TANGENT_PLANE
-!======================================================================
-
-end module utilities
Modified: branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/makeMeshes.sh
===================================================================
--- branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/makeMeshes.sh 2012-12-03 21:05:55 UTC (rev 2338)
+++ branches/ocean_projects/ocean_test_cases_staging/ocean/temporal_convergence/makeMeshes.sh 2012-12-04 17:57:26 UTC (rev 2339)
@@ -130,35 +130,37 @@
#################################################################
## Generate full meshes, with initial conditions, using basin. ##
#################################################################
+echo " Checking out basin"
+svn co https://svn-mpas-model.cgd.ucar.edu/branches/ocean_projects/basin/src basin_checkout 1> /dev/null 2> /dev/null
+cp basin_src/* basin_checkout/.
+echo " Building Basin"
+sed -ne '/subroutine get_init_conditions/ {p; r basin_src/get_init_conds.F' -e ':a; n; /end subroutine get_init_conditions/ {p; b}; ba}; p' basin_checkout/basin.F > temp.F
+mv temp.F basin_checkout/basin.F
+
+if [ -a map ]; then
+ rm map
+fi
+
+cd basin_checkout
+if [ -a Makefile.front ]; then
+ cat Makefile.front > Makefile
+ cat Makefile.end >> Makefile
+else
+ cp Makefile.bak Makefile
+fi
+make clean > /dev/null
+make > /dev/null
+
+cd ../
+cp basin_checkout/map .
+
for VERTLEV in $VERTLEVS
do
if [ -a grid.nc ]; then
unlink grid.nc
fi
- if [ -a map ]; then
- rm map
- fi
-
- echo " Bulding basin for ${VERTLEV} levels"
-
- ## Build Basin
- cd basin_src
- cat basin-template.F | sed "s/*VERTLEVS/${VERTLEV}/g" > basin.F
-
- if [ -a Makefile.front ]; then
- cat Makefile.front > Makefile
- cat Makefile.end >> Makefile
- else
- cp Makefile.bak Makefile
- fi
-
- make clean > /dev/null
- make > /dev/null
- cd ../
- cp basin_src/map .
-
## Call basin, for each perfect hex mesh.
for SPACING in $SPACINGS
do
@@ -181,7 +183,9 @@
ln -s ${TCNAME}_${NAME}.grid.nc grid.nc
mkdir -p dx
+ sed "s/*VERTLEVS/${VERTLEVS}/g" BASIN-namelist.basin.template > namelist.basin
./map > /dev/null
+ rm namelist.basin
unlink grid.nc
@@ -275,17 +279,14 @@
| sed "s/config_stats_interval .*/config_stats_interval = ${STATS}/g" \
> ${BASE_DIR}/namelist.input
done
-
- rm map
done
+rm map
rm -rf dx
rm MPAS-namelist.input.template
rm ${TCNAME}*
rm fort.*
-cd basin_src
-make clean > /dev/null
-rm Makefile
+rm -rf basin_checkout
cd ${CUR_DIR}
</font>
</pre>