<p><b>mpetersen@lanl.gov</b> 2012-10-23 07:31:10 -0600 (Tue, 23 Oct 2012)</p><p>branch commit: Merge basin_pbc_mrp to basin. This adds the ability to produce the bottomDepth variable that is required for partial bottom cells. In addition, basin now reads a namelist, and namelist.basin can point to a standard namelist in the namelists directory.<br>
</p><hr noshade><pre><font color="gray">Index: branches/ocean_projects/basin
===================================================================
--- branches/ocean_projects/basin 2012-10-23 04:36:40 UTC (rev 2258)
+++ branches/ocean_projects/basin 2012-10-23 13:31:10 UTC (rev 2259)
Property changes on: branches/ocean_projects/basin
___________________________________________________________________
Added: svn:mergeinfo
## -0,0 +1 ##
+/branches/ocean_projects/basin_pbc_mrp:2208-2258
\ No newline at end of property
Modified: branches/ocean_projects/basin/README
===================================================================
--- branches/ocean_projects/basin/README 2012-10-23 04:36:40 UTC (rev 2258)
+++ branches/ocean_projects/basin/README 2012-10-23 13:31:10 UTC (rev 2259)
@@ -3,11 +3,17 @@
This program reads in a MPAS grid file (grid.nc) and outputs a modified grid file
(ocean.nc) and it's associated graph.info files for partitions from 2 up to 1024 in powers of 2.
-The purpose of this code is to remove grid cells from any valid MPAS grid.
+The purpose of this code is to remove grid cells from any valid MPAS grid, and to
+add initial condition variables like h, u, u_src, forcing, and tracers.
Please see source file src/basin.F to define the specifics of the output grid file.
-After a grid.nc file has been placed in this directory, simply run the script runit.
+The required files are:
+ grid.nc is an mpas grid, either spherical or Cartesian
+ namelist.basin may point to a standard case in the namelists directory.
+After grid.nc and namelist.basin files has been placed in this directory,
+simply run the script runit.
+
This script will compile basin, run basin (producing an ocean.nc file) and use metis
to partition the graph.info file.
Copied: branches/ocean_projects/basin/namelist.basin (from rev 2258, branches/ocean_projects/basin_pbc_mrp/namelist.basin)
===================================================================
--- branches/ocean_projects/basin/namelist.basin (rev 0)
+++ branches/ocean_projects/basin/namelist.basin 2012-10-23 13:31:10 UTC (rev 2259)
@@ -0,0 +1 @@
+link namelists/namelist.global_realistic
\ No newline at end of file
Deleted: branches/ocean_projects/basin/namelists/namelist.Ilicak2_overflow
===================================================================
--- branches/ocean_projects/basin_pbc_mrp/namelists/namelist.Ilicak2_overflow 2012-10-23 04:36:40 UTC (rev 2258)
+++ branches/ocean_projects/basin/namelists/namelist.Ilicak2_overflow 2012-10-23 13:31:10 UTC (rev 2259)
@@ -1,22 +0,0 @@
-&basin
- nVertLevelsMOD = 10
- on_a_sphere = 'NO'
- sphere_radius = 0.0
- zLevel_thickness = 'equally_spaced'
- bottom_topography = 'Ilicak2_overflow'
- initial_conditions = 'Ilicak2_overflow'
- eliminate_inland_seas=.false.
- load_woce_IC = .false.
- write_OpenDX_flag = .false.
- monthly_forcing = .false.
- cut_domain_from_sphere = .false.
- solid_boundary_in_y = .true.
- solid_boundary_in_x = .false.
- top_layers_without_land = 3
- h_total_max = 2000.0
- u_src_max = 0.0
- f0 = 0.0
- beta = 0.0
- omega = 0.0
- Lx = 3.0e3
-/
Copied: branches/ocean_projects/basin/namelists/namelist.Ilicak2_overflow (from rev 2258, branches/ocean_projects/basin_pbc_mrp/namelists/namelist.Ilicak2_overflow)
===================================================================
--- branches/ocean_projects/basin/namelists/namelist.Ilicak2_overflow (rev 0)
+++ branches/ocean_projects/basin/namelists/namelist.Ilicak2_overflow 2012-10-23 13:31:10 UTC (rev 2259)
@@ -0,0 +1,22 @@
+&basin
+ nVertLevelsMOD = 10
+ on_a_sphere = 'NO'
+ sphere_radius = 0.0
+ zLevel_thickness = 'equally_spaced'
+ bottom_topography = 'Ilicak2_overflow'
+ initial_conditions = 'Ilicak2_overflow'
+ eliminate_inland_seas=.false.
+ load_woce_IC = .false.
+ write_OpenDX_flag = .false.
+ monthly_forcing = .false.
+ cut_domain_from_sphere = .false.
+ solid_boundary_in_y = .true.
+ solid_boundary_in_x = .false.
+ top_layers_without_land = 3
+ h_total_max = 2000.0
+ u_src_max = 0.0
+ f0 = 0.0
+ beta = 0.0
+ omega = 0.0
+ Lx = 3.0e3
+/
Deleted: branches/ocean_projects/basin/namelists/namelist.Ilicak2_overflow_sigma
===================================================================
--- branches/ocean_projects/basin_pbc_mrp/namelists/namelist.Ilicak2_overflow_sigma 2012-10-23 04:36:40 UTC (rev 2258)
+++ branches/ocean_projects/basin/namelists/namelist.Ilicak2_overflow_sigma 2012-10-23 13:31:10 UTC (rev 2259)
@@ -1,22 +0,0 @@
-&basin
- nVertLevelsMOD = 10
- on_a_sphere = 'NO'
- sphere_radius = 0.0
- zLevel_thickness = 'equally_spaced'
- bottom_topography = 'Ilicak2_overflow_sigma'
- initial_conditions = 'Ilicak2_overflow_sigma'
- eliminate_inland_seas=.false.
- load_woce_IC = .false.
- write_OpenDX_flag = .false.
- monthly_forcing = .false.
- cut_domain_from_sphere = .false.
- solid_boundary_in_y = .true.
- solid_boundary_in_x = .false.
- top_layers_without_land = 3
- h_total_max = 2000.0
- u_src_max = 0.0
- f0 = 0.0
- beta = 0.0
- omega = 0.0
- Lx = 3.0e3
-/
Copied: branches/ocean_projects/basin/namelists/namelist.Ilicak2_overflow_sigma (from rev 2258, branches/ocean_projects/basin_pbc_mrp/namelists/namelist.Ilicak2_overflow_sigma)
===================================================================
--- branches/ocean_projects/basin/namelists/namelist.Ilicak2_overflow_sigma (rev 0)
+++ branches/ocean_projects/basin/namelists/namelist.Ilicak2_overflow_sigma 2012-10-23 13:31:10 UTC (rev 2259)
@@ -0,0 +1,22 @@
+&basin
+ nVertLevelsMOD = 10
+ on_a_sphere = 'NO'
+ sphere_radius = 0.0
+ zLevel_thickness = 'equally_spaced'
+ bottom_topography = 'Ilicak2_overflow_sigma'
+ initial_conditions = 'Ilicak2_overflow_sigma'
+ eliminate_inland_seas=.false.
+ load_woce_IC = .false.
+ write_OpenDX_flag = .false.
+ monthly_forcing = .false.
+ cut_domain_from_sphere = .false.
+ solid_boundary_in_y = .true.
+ solid_boundary_in_x = .false.
+ top_layers_without_land = 3
+ h_total_max = 2000.0
+ u_src_max = 0.0
+ f0 = 0.0
+ beta = 0.0
+ omega = 0.0
+ Lx = 3.0e3
+/
Deleted: branches/ocean_projects/basin/namelists/namelist.global_realistic
===================================================================
--- branches/ocean_projects/basin_pbc_mrp/namelists/namelist.global_realistic 2012-10-23 04:36:40 UTC (rev 2258)
+++ branches/ocean_projects/basin/namelists/namelist.global_realistic 2012-10-23 13:31:10 UTC (rev 2259)
@@ -1,25 +0,0 @@
-&basin
- nVertLevelsMOD = 40
- on_a_sphere = 'YES'
- sphere_radius = 6.37122e6
- zLevel_thickness = 'POP_40_zLevel'
- bottom_topography = 'realistic_ETOPO'
- initial_conditions = 'realistic_WOCE'
- eliminate_inland_seas=.true.
- load_woce_IC = .true.
- write_OpenDX_flag = .false.
- monthly_forcing = .true.
- cut_domain_from_sphere = .false.
- solid_boundary_in_y = .false.
- solid_boundary_in_x = .false.
- top_layers_without_land = 3
-
- ! These variables are not needed for realistic global topography:
-
- ! h_total_max = 2000.0
- ! u_src_max = 0.1
- ! f0 = -1.1e-4
- ! beta = 1.4e-11
- ! omega = 7.29212e-5
- ! Lx = 3200.0e3
-/
Copied: branches/ocean_projects/basin/namelists/namelist.global_realistic (from rev 2258, branches/ocean_projects/basin_pbc_mrp/namelists/namelist.global_realistic)
===================================================================
--- branches/ocean_projects/basin/namelists/namelist.global_realistic (rev 0)
+++ branches/ocean_projects/basin/namelists/namelist.global_realistic 2012-10-23 13:31:10 UTC (rev 2259)
@@ -0,0 +1,25 @@
+&basin
+ nVertLevelsMOD = 40
+ on_a_sphere = 'YES'
+ sphere_radius = 6.37122e6
+ zLevel_thickness = 'POP_40_zLevel'
+ bottom_topography = 'realistic_ETOPO'
+ initial_conditions = 'realistic_WOCE'
+ eliminate_inland_seas=.true.
+ load_woce_IC = .true.
+ write_OpenDX_flag = .false.
+ monthly_forcing = .true.
+ cut_domain_from_sphere = .false.
+ solid_boundary_in_y = .false.
+ solid_boundary_in_x = .false.
+ top_layers_without_land = 3
+
+ ! These variables are not needed for realistic global topography:
+
+ ! h_total_max = 2000.0
+ ! u_src_max = 0.1
+ ! f0 = -1.1e-4
+ ! beta = 1.4e-11
+ ! omega = 7.29212e-5
+ ! Lx = 3200.0e3
+/
Modified: branches/ocean_projects/basin/src/basin.F
===================================================================
--- branches/ocean_projects/basin/src/basin.F 2012-10-23 04:36:40 UTC (rev 2258)
+++ branches/ocean_projects/basin/src/basin.F 2012-10-23 13:31:10 UTC (rev 2259)
@@ -18,22 +18,15 @@
! 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
+! Step 1: Link namelist.basin to the correct namelist file.
+! Step 2: Change parameters and flags in namelist file as needed.
+! Step 3: Check get_init_conditions routine for initial T&S, thickness, etc.
+! Step 4: Check define_kmt routine for bottomDepth and kmt (maxLevelCell) variables.
+! Step 5: Check get_dz routine for hZLevel variable.
!
-!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
-integer, parameter :: nx = 50
-! This needs to be changed for correct periodic boundaries
-! Lx is the TOTAL domain width, and needs to be exact for correct periodic
-! boundaries in x.
-real, parameter :: Lx = 3200.0e3 ! 40x80km=3200km
-
! original grid variables
integer :: time, nCells, nEdges, nVertices
integer :: maxEdges, maxEdges2, TWO, vertexDegree, nVertLevels
@@ -48,7 +41,7 @@
real, allocatable, dimension(:) :: areaCell, areaTriangle, dcEdge, dvEdge, angleEdge
real, allocatable, dimension(:,:) :: kiteAreasOnVertex, weightsOnEdge
-real, allocatable, dimension(:) :: fEdge, fVertex, h_s, work1
+real, allocatable, dimension(:) :: fEdge, fVertex, bottomDepth, work1
real, allocatable, dimension(:,:) :: u_src
real, allocatable, dimension(:,:,:) :: u, v, h
real, allocatable, dimension(:,:,:) :: rho
@@ -63,50 +56,14 @@
real(kind=4), allocatable, dimension(:,:,:) :: TAUX_MONTHLY, TAUY_MONTHLY
real, dimension(40) :: dz
-integer, parameter :: nMonths = 12
+integer :: nMonths = 1
-! Step 1: Set the number of Vertical levels
-integer, parameter :: nVertLevelsMOD = 40
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! basin-mod
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-! Need to uncomment the options appropriate for the input grid file. If it's on
-! a sphere, specify the flag "on_a_sphere" and "sphere_radius". Otherwise set
-! them equal to NO and 0.0 respectively
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-! Step 2: Set if the grid is on a sphere or not, and it's radius
-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=.true.
-logical, parameter :: eliminate_inland_seas=.true.
-logical, parameter :: l_woce = .true.
-
-
-! Step 3: Specify some Parameters
- real (kind=8), parameter :: &
- h_total_max = 2000.0, &
- u_max = 0.0, &
- u_src_max = 0.1, & ! max wind stress, N/m2
- beta = 1.4e-11, &
- f0 = -1.1e-4, &
- omega = 7.29212e-5
-
real (kind=8) :: ymid, ytmp, ymax, xmid, xloc, yloc, pert, ymin, distance, r, c1(3), c2(3)
real (kind=8) :: latmid, lattmp, latmax, latmin
integer :: cell1, cell2
-
! new grid variables
+real, allocatable, dimension(:) :: hZLevel, refBottomDepth
integer :: nCellsNew, nEdgesNew, nVerticesNew
integer :: maxEdgesNew, maxEdges2New, TWONew, vertexDegreeNew, nVertLevelsNew
integer, allocatable, dimension(:) :: indexToCellIDNew, indexToEdgeIDNew, indexToVertexIDNew
@@ -121,10 +78,10 @@
real, allocatable, dimension(:) :: areaCellNew, areaTriangleNew, dcEdgeNew, dvEdgeNew, angleEdgeNew
real, allocatable, dimension(:,:) :: kiteAreasOnVertexNew, weightsOnEdgeNew, normalsNew
-real, allocatable, dimension(:) :: fEdgeNew, fVertexNew, h_sNew, hZLevel, referenceBottomDepth
+real, allocatable, dimension(:) :: fEdgeNew, fVertexNew, bottomDepthNew
real, allocatable, dimension(:,:) :: u_srcNew
real, allocatable, dimension(:,:) :: windStressMonthlyNew
-real, allocatable, dimension(:,:,:) :: uNew, vNew, hNew
+real, allocatable, dimension(:,:,:) :: uNew, hNew
real, allocatable, dimension(:,:,:) :: rhoNew, temperatureNew, salinityNew, tracer1New
real, allocatable, dimension(:) :: temperatureRestoreNew, salinityRestoreNew
real, allocatable, dimension(:,:) :: temperatureRestoreMonthlyNew, salinityRestoreMonthlyNew
@@ -132,7 +89,6 @@
! 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
@@ -142,15 +98,90 @@
integer :: iMonth
character(len=80) :: fileNameT, fileNameS, fileNameU
-fileNameT = 'TS/woce_t_ann.3600x2431x42interp.r4.nc'
-fileNameS = 'TS/woce_s_ann.3600x2431x42interp.r4.nc'
-fileNameU = 'TS/ws.old_ncep_1958-2000avg.interp3600x2431.nc'
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! Namelist variables
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! Variables in namelist file
+character (len=32) :: on_a_sphere, zLevel_thickness,bottom_topography, initial_conditions
+logical :: eliminate_inland_seas, load_woce_IC, write_OpenDX_flag, monthly_forcing, &
+ cut_domain_from_sphere, solid_boundary_in_y, solid_boundary_in_x
+integer :: nVertLevelsMOD, top_layers_without_land
+real (kind=8) :: sphere_radius, h_total_max, u_src_max, f0, beta, omega, Lx
+
+! specify namelist
+namelist /basin/ nVertLevelsMOD, on_a_sphere, sphere_radius, &
+ zLevel_thickness, bottom_topography, initial_conditions, &
+ eliminate_inland_seas, load_woce_IC, write_OpenDX_flag, monthly_forcing, &
+ cut_domain_from_sphere, solid_boundary_in_y, solid_boundary_in_x, &
+ top_layers_without_land, h_total_max, u_src_max, f0, beta, omega, Lx
+
+! Default namelist values. Default set for realistic global IC.
+nVertLevelsMOD = 40
+on_a_sphere = 'YES'
+sphere_radius = 6.37122e6
+
+! zLevel thickness options:
+! 'POP_40_zLevel', 'equally_spaced', 'zero'
+zLevel_thickness = 'POP_40_zLevel'
+
+! bottom topography options:
+! 'realistic_ETOPO', 'flat_bottom', 'Ilicak2_overflow'
+bottom_topography = 'realistic_ETOPO'
+
+! initial temperature and salinity options:
+! 'realistic_WOCE', 'constant', 'Ilicak2_overflow', 'Ilicak2_overflow_sigma'
+initial_conditions = 'realistic_WOCE'
+
+eliminate_inland_seas=.true.
+load_woce_IC = .true.
+write_OpenDX_flag = .false.
+monthly_forcing = .true.
+cut_domain_from_sphere = .false.
+solid_boundary_in_y = .false.
+solid_boundary_in_x = .false.
+
+! Set the number of top layers that are not allowed to have land, usually three.
+top_layers_without_land = 3
+
+h_total_max = 2000.0 ! total layer thickness, for equally spaced case
+u_src_max = 0.1 ! max wind stress, N/m2
+f0 = -1.1e-4 ! Coriolis parameter
+beta = 1.4e-11
+omega = 7.29212e-5 ! rotation rate of earth
+
+! This needs to be changed for correct periodic boundaries
+! Lx is the TOTAL domain width, and needs to be exact for correct periodic
+! boundaries in x.
+Lx = 3200.0e3 ! 40x80km=3200km
+
+! Read in namelist
+ open(20,file='namelist.basin',status='old')
+ read(20,basin)
+ close(20)
+
+allocate (hZLevel(nVertLevelsMOD), refBottomDepth(nVertLevelsMOD))
+
+
+if(monthly_forcing) then
+ nMonths = 12
+else
+ nMonths = 1
+end if
+
+fileNameT = 'TS/annual/woce_t_ann.3600x2431x42interp.r4.nc'
+fileNameS = 'TS/annual/woce_s_ann.3600x2431x42interp.r4.nc'
+fileNameU = 'TS/annual/ws.old_ncep_1958-2000avg.interp3600x2431.nc'
+
! get to work
write(6,*) ' starting'
write(6,*)
! get depth profile for later
+write(6,*) ' calling get_dz'
call get_dz
! get grid
@@ -188,129 +219,18 @@
! check the mesh
call error_checking
-write(6,*) ' getting woce t and s '
+if (load_woce_IC) then
+ write(6,*) ' getting woce t and s '
-call read_TS_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-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()
+ call read_TS_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ 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()
-allocate(SST_MONTHLY(nlon,nlat,12), SSS_MONTHLY(nlon,nlat,12))
-allocate(TAUX_MONTHLY(nlon,nlat,12), TAUY_MONTHLY(nlon,nlat,12))
-SST_MONTHLY=0; SSS_MONTHLY=0; TAUX_MONTHLY=0; TAUY_MONTHLY=0
-
-iMonth=1
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.01.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly01.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.01.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=2
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.02.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly02.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.02.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=3
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.03.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly03.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.03.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=4
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.04.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly04.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.04.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=5
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.05.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly05.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.05.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=6
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.06.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly06.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.06.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=7
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.07.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly07.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.07.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=8
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.08.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly08.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.08.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=9
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.09.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly09.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.09.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=10
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.10.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly10.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.10.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=11
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.11.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly11.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.11.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=12
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.12.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly12.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.12.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-do k=1,ndepth
+ do k=1,ndepth
ndata = 0; temp_t=0; temp_s=0
do j=1,nlat
do i=1,nlon
@@ -324,8 +244,123 @@
mTEMP(:,k) = temp_t / float(ndata)
mSALT(:,k) = temp_s / float(ndata)
write(6,*) ndata,mTemp(1,k),mSalt(1,k)
-enddo
+ enddo
+endif
+
+if(monthly_forcing) then
+ allocate(SST_MONTHLY(nlon,nlat,nMonths), SSS_MONTHLY(nlon,nlat,nMonths))
+ allocate(TAUX_MONTHLY(nlon,nlat,nMonths), TAUY_MONTHLY(nlon,nlat,nMonths))
+ SST_MONTHLY=0; SSS_MONTHLY=0; TAUX_MONTHLY=0; TAUY_MONTHLY=0
+ iMonth=1
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.01.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly01.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.01.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=2
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.02.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly02.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.02.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=3
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.03.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly03.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.03.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=4
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.04.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly04.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.04.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=5
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.05.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly05.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.05.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=6
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.06.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly06.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.06.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=7
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.07.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly07.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.07.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=8
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.08.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly08.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.08.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=9
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.09.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly09.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.09.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=10
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.10.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly10.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.10.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=11
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.11.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly11.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.11.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=12
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.12.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly12.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.12.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+end if
+
! generate initial conditions
call get_init_conditions
@@ -340,9 +375,10 @@
call write_graph
! write OpenDx file
-write(6,*) ' calling write_OpenDX'
-write(6,*)
-call write_OpenDX( on_a_sphere, &
+if (write_OpenDX_flag) then
+ write(6,*) ' calling write_OpenDX'
+ write(6,*)
+ call write_OpenDX( on_a_sphere, &
nCellsNew, &
nVerticesNew, &
nEdgesNew, &
@@ -365,11 +401,11 @@
areaCellNew, &
maxLevelCellNew, &
meshDensityNew, &
- depthCell, &
+ bottomDepthNew, &
temperatureNew(1,1,:), &
kiteAreasOnVertexNew )
+endif
-
!do iCell=1,nCellsNew
!ulon = 1.0; ulat = 0.0
!xin = xCellNew(iCell); yin = yCellNew(iCell); zin = zCellNew(iCell)
@@ -412,7 +448,11 @@
end subroutine write_graph
-!Step 4: Check the Initial conditions routine get_init_conditions
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! Step 3: Check get_init_conditions routine for initial T&S, thickness, etc.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine get_init_conditions
implicit none
real :: halfwidth, dtr, pi, p(3), q(3), xin, yin, zin, ulon, ulat, stress, n1, n2, distance, r, temp_t, temp_s
@@ -423,22 +463,78 @@
pi = 4.0*atan(1.0)
dtr = pi/180.0
+! defaults
hNew = 100.0
temperatureNew = 1.0
salinityNew = 1.0
tracer1New = 1.0
uNew = 0
-vNew = 0
+u_srcNew = 0
+rhoNew = 1025.0
-if(.not.real_bathymetry) then
+if (initial_conditions.eq.'uniform_TS') then
- write(6,*) ' not using real bathymetry'
+ temperatureNew = 10.0
+ salinityNew = 34.0
+ tracer1New = 1.0
+ uNew = 0.0
+ u_srcNew = 0.0
+ ! change thickness of bottom level
+ do iCell=1,nCellsNew
+ do k=1, nVertLevelsMod
+ hNew(1,k,iCell) = hZLevel(k)
+ enddo
+ enddo
+
+elseif (initial_conditions.eq.'Ilicak2_overflow') then
+
+ do i = 1,nCellsNew
+ if(yCellNew(i) < 20000) then
+ temperatureNew(1,:,i) = 10.0
+ else
+ temperatureNew(1,:,i) = 20.0
+ endif
+ enddo
+
+ salinityNew(1,:,:) = 35.0
+ Tracer1New(1,:,:) = 1.0
+ uNew = 0.0
+ u_srcNew = 0.0
+
+ do iCell=1,nCellsNew
+ do k=1,nVertLevelsMod
+ hNew(1,k,iCell) = hZLevel(k)
+ enddo
+ enddo
+
+elseif (initial_conditions.eq.'Ilicak2_overflow_sigma') then
+
+ do i = 1,nCellsNew
+ if(yCellNew(i) < 20000) then
+ temperatureNew(1,:,i) = 10.0
+ else
+ temperatureNew(1,:,i) = 20.0
+ endif
+ enddo
+
+ salinityNew(1,:,:) = 35.0
+ Tracer1New(1,:,:) = 1.0
+ uNew = 0.0
+ u_srcNew = 0.0
+
+ do iCell=1,nCellsNew
+ do k=1,nVertLevelsMod
+ hNew(1,k,iCell) = bottomDepthNew(iCell) / nVertLevelsMOD
+ enddo
+ enddo
+
+elseif (initial_conditions.eq.'isopycnal') then
+
fEdgeNew(:) = 0.0
fVertexNew(:) = 0.0
- h_sNew(:) = 0.0
+ bottomDepthNew(:) = 0.0
uNew(:,:,:) = 0.0
- vNew(:,:,:) = 0.0
! basin-mod
! setting for three levels - Set h values for isopycnal system
@@ -447,17 +543,12 @@
hNew(1,1,:) = 500.0
hNew(1,2,:) = 1250.0
hNew(1,3,:) = 3250.0
- h_sNew(:) = -( hNew(1,1,:) + hNew(1,2,:) + hNew(1,3,:) )
+ bottomDepthNew(:) = -( hNew(1,1,:) + hNew(1,2,:) + hNew(1,3,:) )
else
hNew(1,1,:) = 3250.0
- h_sNew(:) = -( hNew(1,1,:) )
+ bottomDepthNew(:) = -( hNew(1,1,:) )
endif
- hZLevel = hNew(1,:,1)
- referenceBottomDepth(1) = hZLevel(1)
- do k = 2,nVertLevels
- referenceBottomDepth(k) = referenceBottomDepth(k-1) + hZLevel(k)
- end do
! basin-mod
! Noise is meant to make the flow unstable at some point
@@ -500,7 +591,7 @@
! basin-mod
! set salinity for isopycnal levels
write(6,*) ' setting salinity'
- if(on_a_sphere.eq.'YES ') then
+ if(on_a_sphere.eq.'YES') then
latmin = minval(latCellNew)
latmax = maxval(latCellNew)
r = 10.0*dtr
@@ -544,7 +635,7 @@
! set forcing for isopycnal levels
write(6,*) 'setting u_src - wind forcing'
u_srcNew = 0.0
- if(on_a_sphere.eq.'YES ') then
+ if(on_a_sphere.eq.'YES') then
latmin = -60*dtr
latmax = -10*dtr
latmid = -35*dtr
@@ -599,33 +690,10 @@
endif
write(6,*) ' u_srcNew ', minval(u_srcNew), maxval(u_srcNew)
- ! basin-mod
- ! set coriolis parameter for grid
- write(6,*) ' setting Coriolis parameter'
- if(on_a_sphere.eq.'YES ') then
- do i = 1,nVerticesNew
- fVertexNew(i) = 2.0 * omega * sin(latVertexNew(i))
- enddo
- do i = 1,nEdgesNew
- fEdgeNew(i) = 2.0 * omega * sin(latEdgeNew(i))
- enddo
- else
- do i = 1,nVerticesNew
- fVertexNew(i) = f0 + (yVertexNew(i) - ymid) * beta
- enddo
+elseif (initial_conditions.eq.'realistic_WOCE') then
- do i = 1,nEdgesNew
- fEdgeNew(i) = f0 + (yEdgeNew(i) - ymid) * beta
- enddo
- endif
- write(6,*) ' done not real bathymetry'
-endif ! if(.not.real_bathymetry) then
-
-
-if(real_bathymetry) then
-
u_srcNew = 0.0
windStressMonthlyNew = 0.0
do iEdge=1,nEdgesNew
@@ -659,28 +727,32 @@
u_srcNew(1,iEdge) = ux*q(1) + uy*q(2) + uz*q(3)
endif
- do iMonth=1,12
- ulon = TAUX_MONTHLY(ix,iy,iMonth)
- ulat = TAUY_MONTHLY(ix,iy,iMonth)
-! if(abs(ulon).gt.1.0.or.abs(ulat).gt.1.0) then
-! ulon=0.0
-! ulat=0.0
-! endif
- call transform_from_lonlat_to_xyz(xin,yin,zin,ulon,ulat,ux,uy,uz)
- if(boundaryEdgeNew(1,iEdge).eq.1) then
- windStressMonthlyNew(iMonth,iEdge) = 0.0
+ if(monthly_forcing) then
+ do iMonth=1,nMonths
+ ulon = TAUX_MONTHLY(ix,iy,iMonth)
+ ulat = TAUY_MONTHLY(ix,iy,iMonth)
+ ! if(abs(ulon).gt.1.0.or.abs(ulat).gt.1.0) then
+ ! ulon=0.0
+ ! ulat=0.0
+ ! endif
+ call transform_from_lonlat_to_xyz(xin,yin,zin,ulon,ulat,ux,uy,uz)
+ if(boundaryEdgeNew(1,iEdge).eq.1) then
+ windStressMonthlyNew(iMonth,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)
+ ! repeat
+ windStressMonthlyNew(iMonth,iEdge) = ux*q(1) + uy*q(2) + uz*q(3)
+ ! windStressMonthlyNew(iMonth,iEdge) = u_srcNew(1,iEdge)
+ endif
+ enddo
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)
- ! repeat
- windStressMonthlyNew(iMonth,iEdge) = ux*q(1) + uy*q(2) + uz*q(3)
- ! windStressMonthlyNew(iMonth,iEdge) = u_srcNew(1,iEdge)
- endif
- enddo
+ windStressMonthlyNew(:,:) = 0.0
+ end if
enddo
@@ -695,12 +767,11 @@
enddo
! update T and S field with WOCE data
-if(l_woce) then
+if(load_woce_IC) 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
+ if(mod(iCell,100).eq.0) write(6,*) 'load_woce_IC t and s',iCell
rlon = lonCellNew(iCell)/dtr
rlat = latCellNew(iCell)/dtr
ix = nint(rlon/0.1 - 0.05) + nlon + 1
@@ -780,88 +851,124 @@
temperatureRestoreNew(:) = temperatureNew(1,1,:)
salinityRestoreNew(:) = salinityNew(1,1,:)
-do iMonth=1,12
-iNoData = 0
-do iCell=1,nCellsNew
- if(mod(iCell,100).eq.0) write(6,*) 'l_woce t and s RESTORE',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 ! j
- k=1
- ndata = 0; temp_t = 0; temp_s = 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
+if(monthly_forcing) then
+ do iMonth=1,nMonths
+ iNoData = 0
+ do iCell=1,nCellsNew
+ if(mod(iCell,100).eq.0) write(6,*) 'load_woce_IC t and s RESTORE',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 ! j
+ k=1
+ ndata = 0; temp_t = 0; temp_s = 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(SST_MONTHLY(ixt,iyt,iMonth).gt.-10.0.and.flag_lat) then
+ ndata = ndata + 1
+ temp_t = temp_t + SST_MONTHLY(ixt,iyt,iMonth)
+ temp_s = temp_s + SSS_MONTHLY(ixt,iyt,iMonth)
+ endif
+ enddo !j
+ enddo !i
+
+ if(ndata.gt.0) then
+ temperatureRestoreMonthlyNew(iMonth,iCell) = temp_t / float(ndata)
+ salinityRestoreMonthlyNew(iMonth,iCell) = temp_s / float(ndata)
+ else
+ temperatureRestoreMonthlyNew(iMonth,iCell) = temperatureNew(1,1,iCell)
+ salinityRestoreMonthlyNew(iMonth,iCell) = salinityNew(1,1,iCell)
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(SST_MONTHLY(ixt,iyt,iMonth).gt.-10.0.and.flag_lat) then
- ndata = ndata + 1
- temp_t = temp_t + SST_MONTHLY(ixt,iyt,iMonth)
- temp_s = temp_s + SSS_MONTHLY(ixt,iyt,iMonth)
- endif
- enddo !j
- enddo !i
-
- if(ndata.gt.0) then
- temperatureRestoreMonthlyNew(iMonth,iCell) = temp_t / float(ndata)
- salinityRestoreMonthlyNew(iMonth,iCell) = temp_s / float(ndata)
- else
- temperatureRestoreMonthlyNew(iMonth,iCell) = temperatureNew(1,1,iCell)
- salinityRestoreMonthlyNew(iMonth,iCell) = salinityNew(1,1,iCell)
- endif
-
+
+ enddo ! iCell
+ enddo ! iMonth
+
+ ! do a couple of smoothing passes
+ do iter=1,5
+ do iCell=1,nCellsNew
+ k=1
+ ndata=1
+ temp_t = temperatureRestoreMonthlyNew(iMonth,iCell)
+ temp_s = salinityRestoreMonthlyNew(iMonth,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 + temperatureRestoreMonthlyNew(iMonth,iCell)
+ temp_s = temp_s + salinityRestoreMonthlyNew(iMonth,iCell)
+ ndata = ndata + 1
+ endif
+ endif
+ enddo ! j
+ temperatureRestoreMonthlyNew(iMonth,iCell) = temp_t / ndata
+ salinityRestoreMonthlyNew(iMonth,iCell) = temp_s / ndata
enddo ! iCell
-enddo ! iMonth
+ enddo ! iter
+else
+ temperatureRestoreMonthlyNew(:,:) = 0.0
+ salinityRestoreMonthlyNew(:,:) = 0.0
+end if
-! do a couple of smoothing passes
-do iter=1,5
-do iCell=1,nCellsNew
- k=1
- ndata=1
- temp_t = temperatureRestoreMonthlyNew(iMonth,iCell)
- temp_s = salinityRestoreMonthlyNew(iMonth,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 + temperatureRestoreMonthlyNew(iMonth,iCell)
- temp_s = temp_s + salinityRestoreMonthlyNew(iMonth,iCell)
- ndata = ndata + 1
- endif
- endif
- enddo ! j
- temperatureRestoreMonthlyNew(iMonth,iCell) = temp_t / ndata
- salinityRestoreMonthlyNew(iMonth,iCell) = temp_s / ndata
-enddo ! iCell
-enddo ! iter
+endif ! load_woce_IC
-endif ! l_woce
-
!repeat
!do iMonth=1,12
! temperatureRestoreMonthlyNew(iMonth,:) = temperatureRestoreNew(:)
! salinityRestoreMonthlyNew(iMonth,:) = salinityRestoreNew(:)
!enddo
-endif ! real_bathymetry
+else
+ print *, ' Incorrect choice of initial_conditions: ',initial_conditions
+ stop
+
+endif ! initial_conditions
+
+
+if (trim(initial_conditions).ne.'realistic_WOCE') then
+
+ ! basin-mod
+ ! set coriolis parameter for grid
+ write(6,*) ' setting Coriolis parameter'
+ if(on_a_sphere.eq.'YES') then
+ do i = 1,nVerticesNew
+ fVertexNew(i) = 2.0 * omega * sin(latVertexNew(i))
+ enddo
+
+ do i = 1,nEdgesNew
+ fEdgeNew(i) = 2.0 * omega * sin(latEdgeNew(i))
+ enddo
+ else
+ do i = 1,nVerticesNew
+ fVertexNew(i) = f0 + (yVertexNew(i) - ymid) * beta
+ enddo
+
+ do i = 1,nEdgesNew
+ fEdgeNew(i) = f0 + (yEdgeNew(i) - ymid) * beta
+ enddo
+ endif
+
+endif
+
write(6,*) ' done get_init_conditions'
end subroutine get_init_conditions
@@ -1086,7 +1193,7 @@
allocate(fEdge(nEdges))
allocate(fVertex(nVertices))
-allocate(h_s(nCells))
+allocate(bottomDepth(nCells))
allocate(work1(nCells))
allocate(u_src(nVertLevels,nEdges))
allocate(u(1,nVertLevels,nEdges))
@@ -1105,7 +1212,7 @@
cellsOnCell=0; verticesOnCell=0; verticesOnEdge=0
edgesOnVertex=0; cellsOnVertex=0; kiteAreasOnVertex=0
-fEdge=0; fVertex=0; h_s=0; u_src=0; work1=0
+fEdge=0; fVertex=0; bottomDepth=0; u_src=0; work1=0
u=0; v=0; h=0; rho=0
@@ -1149,7 +1256,7 @@
kiteAreasOnVertex, &
fEdge, &
fVertex, &
- h_s, &
+ bottomDepth, &
u, &
v, &
h &
@@ -1194,7 +1301,7 @@
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,*) ' bottomDepth : ', minval(bottomDepth), maxval(bottomDepth)
write(6,*) ' u : ', minval(u), maxval(u)
write(6,*) ' v : ', minval(v), maxval(v)
write(6,*) ' h : ', minval(h), maxval(h)
@@ -1205,7 +1312,7 @@
subroutine write_grid
implicit none
-if(on_a_sphere.eq.'YES ') then
+if(on_a_sphere.eq.'YES') then
xCellNew = xCellNew * sphere_radius
yCellNew = yCellNew * sphere_radius
zCellNew = zCellNew * sphere_radius
@@ -1275,13 +1382,12 @@
kiteAreasOnVertexNew, &
fEdgeNew, &
fVertexNew, &
- h_sNew, &
+ bottomDepthNew, &
boundaryEdgeNew, &
boundaryVertexNew, &
u_srcNew, &
windStressMonthlyNew, &
uNew, &
- vNew, &
hNew, &
rhoNew, &
temperatureNew, &
@@ -1292,12 +1398,12 @@
temperatureRestoreMonthlyNew, &
salinityRestoreMonthlyNew, &
hZLevel, &
- referenceBottomDepth &
+ refBottomDepth &
)
call write_netcdf_finalize
-if(on_a_sphere.eq.'YES ') then
+if(on_a_sphere.eq.'YES') then
xCellNew = xCellNew / sphere_radius
yCellNew = yCellNew / sphere_radius
zCellNew = zCellNew / sphere_radius
@@ -1316,62 +1422,28 @@
end subroutine write_grid
-! Step 5: Check the depth routine define_kmt
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! Step 4: Check define_kmt routine for bottomDepth and kmt (maxLevelCell) variables
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine define_kmt
implicit none
real (kind=4), allocatable, dimension(:) :: x,y, work_kmt
real (kind=4), allocatable, dimension(:,:) :: ztopo
integer :: nx, ny, inx, iny, ix, iy, kmt_neighbor_max
real :: pi, dtr, zdata, rlon, rlat, r, ymin, ymax, xmin, xmax
-real :: latmin, latmax, lonmin, lonmax
+real :: latmin, latmax, lonmin, lonmax, ridgeDepth, maxdc
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
+if (bottom_topography.eq.'realistic_ETOPO') then
- ! ! 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))
@@ -1391,8 +1463,14 @@
write(6,*) minval(ztopo), maxval(ztopo)
do iCell=1,nCells
+
+ ! Convert from radians to degrees
rlon = lonCell(iCell) / dtr
rlat = latCell(iCell) / dtr
+
+ ! Find nearest coordinate in topo file.
+ ! This is 1/30th degree topo data, so multiply degrees by 30 and
+ ! round to get index.
ix = nint((rlon+180)*30) + nx
ix = mod(ix,nx)+1
iy = nint((rlat+90 )*30)
@@ -1401,8 +1479,10 @@
zdata = ztopo(ix,iy)
+ ! zdata is less than zero for ocean points.
if(zdata.lt.0.0) then
zdata = -zdata
+ bottomDepth(iCell) = zdata
r = 0
kmt_flag=.false.
do k=1,nVertLevelsMod
@@ -1414,8 +1494,15 @@
endif
endif
enddo
- if(kmt(iCell).eq.0) kmt(iCell)=nVertLevelsMod
+
+ ! zdata is deeper than deepest cell
+ if (kmt(iCell).eq.0) then
+ kmt(iCell)=nVertLevelsMod
+ bottomDepth(iCell) = refBottomDepth(nVertLevelsMod)
+ endif
+
! write(6,*) kmt(iCell)
+
endif
! if(zdata.lt.0.0) kmt(iCell) = nVertLevelsMod
@@ -1425,8 +1512,108 @@
deallocate(x)
deallocate(y)
deallocate(ztopo)
+
+elseif (bottom_topography.eq.'Ilicak2_overflow') then
+
+ kmt = nVertLevelsMOD
+
+ ridgeDepth = 500.0
+ do iCell = 1,nCells
+ ! From Mehmet Ilicak:
+ ! depth=2000
+ ! val1 = 500 is top of ridge
+ ! h(i,j) = val1 + 0.5*(depth-val1) * (1.0+TANH((lon(i,j)-40000.0)/7000.0))
+ bottomDepth(iCell) = ridgeDepth + 0.5*(h_total_max-ridgeDepth) * (1.0+tanh((yCell(iCell)-40000.0)/7000.0))
+
+ if (bottomDepth(iCell).lt.0.0.or. &
+ bottomDepth(iCell).gt.refBottomDepth(nVertLevelsMOD)) then
+ print *, 'error: bottomDepth cannot be less than zero or greater than refBottomDepth(nVertLevels)'
+ print *, 'iCell, bottomDepth(iCell):', iCell, bottomDepth(iCell)
+ exit
+ end if
+
+ do k=1,nVertLevelsMOD
+ if (bottomDepth(iCell).le.refBottomDepth(k)) then
+ kmt(iCell) = k
+ exit
+ endif
+ end do
+
+ enddo
+
+elseif (bottom_topography.eq.'Ilicak2_overflow_sigma') then
+
+ ridgeDepth = 500.0
+ do iCell = 1,nCells
+ ! From Mehmet Ilicak:
+ ! depth=2000
+ ! val1 = 500 is top of ridge
+ ! h(i,j) = val1 + 0.5*(depth-val1) * (1.0+TANH((lon(i,j)-40000.0)/7000.0))
+ bottomDepth(iCell) = ridgeDepth + 0.5*(h_total_max-ridgeDepth) * (1.0+tanh((yCell(iCell)-40000.0)/7000.0))
+
+ if (bottomDepth(iCell).lt.0.0.or. &
+ bottomDepth(iCell).gt.refBottomDepth(nVertLevelsMOD)) then
+ print *, 'error: bottomDepth cannot be less than zero or greater than refBottomDepth(nVertLevels)'
+ print *, 'iCell, bottomDepth(iCell):', iCell, bottomDepth(iCell)
+ exit
+ end if
+
+ enddo
+
+ ! for sigma coordinates, set kmt to the max level.
+ kmt = nVertLevelsMOD
+
+elseif (bottom_topography.eq.'flat_bottom') then
+
+ kmt = nVertLevelsMOD
+
+else
+
+ print *, ' Incorrect choice of bottom_topography: ',bottom_topography
+ stop
+
endif
+if (cut_domain_from_sphere) then
+ 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
+endif
+
+if (solid_boundary_in_y) then
+ 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
+endif
+
+if (solid_boundary_in_x) then
+ maxdc = maxval(dcEdge)
+ xmin = minval(xCell)
+ write(6,*) ' minimum xCell ', xmin
+ xmax = maxval(xCell)
+ write(6,*) ' maximum xCell ', xmax
+ where(xCell.lt.xmin+maxdc/1.5) kmt = 0
+ where(xCell.gt.xmax-maxdc/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)
+
+
! Eliminate isolated ocean cells, and make these isolated deep cells
! flush with the deepest neighbor.
do iCell=1,nCells
@@ -1435,7 +1622,10 @@
iCell1 = cellsOnCell(j,iCell)
kmt_neighbor_max = max(kmt_neighbor_max,kmt(iCell1))
enddo
- kmt(iCell) = min(kmt(iCell),kmt_neighbor_max)
+ if (kmt(iCell).gt.kmt_neighbor_max) then
+ kmt(iCell) = kmt_neighbor_max
+ bottomDepth(iCell) = refBottomDepth(kmt(iCell))
+ endif
enddo
if(eliminate_inland_seas) then
@@ -1445,10 +1635,12 @@
KMT)
endif
-if(real_bathymetry) then
- where(kmt.eq.1) kmt=3
- where(kmt.eq.2) kmt=3
-endif
+! do not allow land or PBCs in top layers
+k = top_layers_without_land
+where(kmt.gt.0.and.kmt.le.k)
+ bottomDepth = refBottomDepth(k)
+ kmt=k
+endwhere
end subroutine define_kmt
@@ -1561,18 +1753,14 @@
allocate(areaCellNew(nCellsNew))
allocate(areaTriangleNew(nVerticesNew))
allocate(maxLevelCellNew(nCellsNew))
-allocate(depthCell(nCellsNew))
allocate(fEdgeNew(nEdgesNew))
allocate(fVertexNew(nVerticesNew))
-allocate(h_sNew(nCellsNew))
+allocate(bottomDepthNew(nCellsNew))
allocate(u_srcNew(nVertLevelsNew,nEdgesNew))
-allocate(windStressMonthlyNew(12,nEdgesNew))
+allocate(windStressMonthlyNew(nMonths,nEdgesNew))
allocate(uNew(1,nVertLevelsNew,nEdgesNew))
-allocate(vNew(1,nVertLevelsNew,nEdgesNew))
allocate(hNew(1,nVertLevelsNew,nCellsNew))
-allocate(hZLevel(nVertLevelsNew))
-allocate(referenceBottomDepth(nVertLevelsNew))
allocate(rhoNew(1,nVertLevelsNew,nCellsNew))
allocate(temperatureNew(1,nVertLevelsNew,nCellsNew))
allocate(salinityNew(1,nVertLevelsNew,nCellsNew))
@@ -1581,16 +1769,16 @@
allocate(temperatureRestoreNew(nCellsNew))
allocate(salinityRestoreNew(nCellsNew))
-allocate(temperatureRestoreMonthlyNew(12,nCellsNew))
-allocate(salinityRestoreMonthlyNew(12,nCellsNew))
+allocate(temperatureRestoreMonthlyNew(nMonths,nCellsNew))
+allocate(salinityRestoreMonthlyNew(nMonths,nCellsNew))
xCellNew=0; yCellNew=0; zCellNew=0; latCellNew=0; lonCellNew=0; meshDensityNew=1.0; meshSpacingNew=0.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; windStressMonthlyNew = 0
-uNew=0; vNew=0; hNew=0; rhoNew=0
+fEdgeNew=0; fVertexNew=0; bottomDepthNew=0; u_srcNew=0; windStressMonthlyNew = 0
+uNew=0; hNew=0; rhoNew=0
temperatureNew=0; salinityNew=0; tracer1New=0;
temperatureRestoreNew = 0.0
@@ -1611,7 +1799,7 @@
meshDensityNew(jNew)=meshDensity(i)
areaCellNew(jNew)=areaCell(i)
maxLevelCellNew(jNew) = kmt(i)
- depthCell(jNew) = kmt(i)
+ bottomDepthNew(jNew) = bottomDepth(i)
endif
enddo
@@ -1661,6 +1849,7 @@
deallocate(lonVertex)
deallocate(dcEdge)
deallocate(dvEdge)
+!deallocate(bottomDepth)
end subroutine map_vectors
@@ -1846,7 +2035,7 @@
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
+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)
@@ -1874,10 +2063,28 @@
end subroutine map_connectivity
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! Step 5: Check get_dz routine for hZLevel variable
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine get_dz
integer k
+if (zLevel_thickness.eq.'zero') then
+
+ ! hZLevel not used for isopycnal setting, just set to zero.
+ hZLevel = 0.0
+
+elseif (zLevel_thickness.eq.'equally_spaced') then
+
+ write(6,*) ' equally spaced zLevels'
+ do i = 1, nVertLevelsMOD
+ hZLevel(i) = h_total_max / nVertLevelsMOD
+ end do
+
+elseif(zLevel_thickness.eq.'POP_40_zLevel') then
+
dz( 1) = 1001.244 ! 5.006218 10.01244
dz( 2) = 1011.258 ! 15.06873 20.12502
dz( 3) = 1031.682 ! 25.28342 30.44183
@@ -1921,11 +2128,26 @@
dz = dz / 100.0
- write(6,*)
- do k=1,40
- write(6,*) k,dz(k)
+ hZLevel = dz
+
+else
+
+ print *, ' Incorrect choice of zLevel_thickness: ',zLevel_thickness
+ stop
+
+endif
+
+ refBottomDepth(1) = hZLevel(1)
+ do k = 2,nVertLevelsMod
+ refBottomDepth(k) = refBottomDepth(k-1) + hZLevel(k)
+ end do
+
+ write(6,*) ' k hZLevel refBottomDepth'
+ do k=1,nVertLevelsMod
+ write(6,'(i5,2f10.2)') k,hZLevel(k), refBottomDepth(k)
enddo
write(6,*)
end subroutine get_dz
+
end program map_to_basin
Modified: branches/ocean_projects/basin/src/module_cullLoops.F
===================================================================
--- branches/ocean_projects/basin/src/module_cullLoops.F 2012-10-23 04:36:40 UTC (rev 2258)
+++ branches/ocean_projects/basin/src/module_cullLoops.F 2012-10-23 13:31:10 UTC (rev 2259)
@@ -1,13 +1,13 @@
module cullLoops
- public :: eliminateLoops
+ public :: eliminateLoops
- contains
+ 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)
+ 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
@@ -30,274 +30,55 @@
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)
+ integer, dimension(:), pointer :: cellStack
+ integer, dimension(:), pointer :: oceanMask
+ integer :: iCellStart, nStack, addedCells
+ real :: latStart, lonStart
- ! 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,*) 'Culling inland seas.....'
- write(6,*) 'working on : ',iCell, KMT(iCell)
-
- ! skip over land cells
- if(KMT(iCell).eq.0) then
- write(6,*) ' skipping : ', iCell
- cycle
- endif
+ allocate(cellStack(nCells/2))
+ allocate(oceanMask(nCells))
- ! the working cell will be jCell, so set jCell=iCell to start
- jCell = iCell
- ! write(6,*) 'setting jCell: ', jCell
+ oceanMask = 0
+ addedCells = 0
- 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?
+ iCellStart = maxloc(kmt, dim=1)
- do while (.not.atBoundary)
+ write(6,*) 'Starting index. ', iCellStart
+ write(6,*) 'lat, lon: ', latCell(iCellStart), lonCell(iCellStart)
+ write(6,*) 'Starting kmt: ', kmt(iCellStart)
- ! 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
+ nStack = 1
+ cellStack(nStack) = iCellStart
+ oceanMask(iCellStart) = 1
+ addedCells = 1
- ! 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
+ do while(nStack > 0)
+ oCell = cellStack(nStack)
+ nStack = nStack - 1
+ write(6,*) ' Working on cell ', oCell, addedCells, nStack
- ! if(moveSouth.and..not.atBoundary) write(6,*) ' pushing on to the south : ', jCell
-
- enddo ! .not.atBoundary
+ do i = 1, nEdgesOnCell(oCell)
+ iCell = cellsOnCell(i, oCell)
- ! 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)
+ if(kmt(iCell) > 0 .and. oceanMask(iCell) == 0) then
+ nStack = nStack + 1
+ cellStack(nStack) = iCell
+ oceanMask(iCell) = 1
+ addedCells = addedCells + 1
+ end if
+ end do
+ end do
- ! start the counter at 1
- iEdgeCounter = 1
- edgeList(:) = 0
+ where(oceanMask == 0) kmt(:) = 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
+ write(6,*) addedCells, ' total cells have been in the stack.'
+ write(6,*) 'Done culling inland seas.....'
- connected = .false.
- LeftTurns = 0; RightTurns = 0
- do while (.not.connected)
+ deallocate(cellStack)
+ deallocate(oceanMask)
- call moveAhead(xCell,yCell,zCell,xVertex,yVertex,zVertex, &
- oCell,lCell,iSharedEdge,nCells,nEdges,nVertices,maxEdges,vertexDegree, &
- verticesOnEdge,cellsOnVertex,iCellAhead)
+ end subroutine eliminateLoops
- ! 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), v3(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)
-
- v1 = v1 - ocean
- v2 = v2 - ocean
- v3 = land - ocean
-
- ocean(:) = ocean(:) / sqrt( ocean(1)**2 + ocean(2)**2 + ocean(3)**2)
- land(:) = land(:) / sqrt( land(1)**2 + land(2)**2 + land(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)
- v3(:) = v3(:) / sqrt( v3(1)**2 + v3(2)**2 + v3(3)**2)
-
- call cross_product_in_R3(v3,v1,cross1)
- call cross_product_in_R3(v3,v2,cross2)
-
- d1 = ocean(1)*cross1(1) + ocean(2)*cross1(2) + ocean(3)*cross1(3)
- d2 = ocean(1)*cross2(1) + ocean(2)*cross2(2) + ocean(3)*cross2(3)
-
- kCell = 0
-
- ! write(6,*) lCell, oCell
- ! write(6,*) cellsOnVertex(:,iVertex1), d1
- ! write(6,*) cellsOnVertex(:,iVertex2), d2
- ! write(6,*) v1(:)
- ! write(6,*) v2(:)
-
- if(d1.gt.d2) then
- do i=1,vertexDegree
- kCell=cellsOnVertex(i,iVertex1)
- if(kCell.ne.oCell.and.kCell.ne.lCell) iCellAhead=kCell
- enddo
- endif
-
- if(d2.gt.d1) 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
Modified: branches/ocean_projects/basin/src/module_write_netcdf.F
===================================================================
--- branches/ocean_projects/basin/src/module_write_netcdf.F 2012-10-23 04:36:40 UTC (rev 2258)
+++ branches/ocean_projects/basin/src/module_write_netcdf.F 2012-10-23 13:31:10 UTC (rev 2259)
@@ -50,13 +50,12 @@
integer :: wrVarIDkiteAreasOnVertex
integer :: wrVarIDfEdge
integer :: wrVarIDfVertex
- integer :: wrVarIDh_s
+ integer :: wrVarIDbottomDepth
integer :: wrVarIDu
integer :: wrVarIDboundaryEdge
integer :: wrVarIDboundaryVertex
integer :: wrVarIDu_src
integer :: wrVarIDwindStressMonthly
- integer :: wrVarIDv
integer :: wrVarIDh
integer :: wrVarIDrho
integer :: wrVarIDtemperature
@@ -67,7 +66,7 @@
integer :: wrVarIDtemperatureRestoreMonthly
integer :: wrVarIDsalinityRestoreMonthly
integer :: wrVarIDhZLevel
- integer :: wrVarIDreferenceBottomDepth
+ integer :: wrVarIDrefBottomDepth
integer :: wrLocalnCells
integer :: wrLocalnEdges
@@ -103,7 +102,7 @@
integer, intent(in) :: vertexDegree
integer, intent(in) :: nMonths
character (len=16) :: on_a_sphere
- double precision :: sphere_radius
+ real*8 :: sphere_radius
integer :: nferr
@@ -131,7 +130,7 @@
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, 'nMonths', 12, wrDimIDnMonths)
+ nferr = nf_def_dim(wr_ncid, 'nMonths', nMonths, wrDimIDnMonths)
nferr = nf_def_dim(wr_ncid, 'Time', NF_UNLIMITED, wrDimIDTime)
!
@@ -226,7 +225,7 @@
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)
+ nferr = nf_def_var(wr_ncid, 'bottomDepth', NF_DOUBLE, 1, dimlist, wrVarIDbottomDepth)
dimlist( 1) = wrDimIDnCells
nferr = nf_def_var(wr_ncid, 'temperatureRestore', NF_DOUBLE, 1, dimlist, wrVarIDtemperatureRestore)
dimlist( 1) = wrDimIDnCells
@@ -242,7 +241,7 @@
dimlist( 1) = wrDimIDnVertLevels
nferr = nf_def_var(wr_ncid, 'hZLevel', NF_DOUBLE, 1, dimlist, wrVarIDhZLevel)
dimlist( 1) = wrDimIDnVertLevels
- nferr = nf_def_var(wr_ncid, 'referenceBottomDepth', NF_DOUBLE, 1, dimlist, wrVarIDreferenceBottomDepth)
+ nferr = nf_def_var(wr_ncid, 'refBottomDepth', NF_DOUBLE, 1, dimlist, wrVarIDrefBottomDepth)
dimlist( 1) = wrDimIDnVertLevels
dimlist( 2) = wrDimIDnEdges
dimlist( 3) = wrDimIDTime
@@ -262,10 +261,6 @@
nferr = nf_def_var(wr_ncid, 'windStressMonthly', NF_DOUBLE, 2, dimlist, wrVarIDwindStressMonthly)
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)
@@ -337,13 +332,12 @@
kiteAreasOnVertex, &
fEdge, &
fVertex, &
- h_s, &
+ bottomDepth, &
boundaryEdge, &
boundaryVertex, &
u_src, &
windStressMonthly, &
u, &
- v, &
h, &
rho, &
temperature, &
@@ -354,7 +348,7 @@
temperatureRestoreMonthly, &
salinityRestoreMonthly, &
hZLevel, &
- referenceBottomDepth &
+ refBottomDepth &
)
implicit none
@@ -401,13 +395,12 @@
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
+ real (kind=8), dimension(:), intent(in) :: bottomDepth
integer, dimension(:,:), intent(in) :: boundaryEdge
integer, dimension(:,:), intent(in) :: boundaryVertex
real (kind=8), dimension(:,:), intent(in) :: u_src
real (kind=8), dimension(:,:), intent(in) :: windStressMonthly
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
@@ -418,7 +411,7 @@
real (kind=8), dimension(:,:), intent(in) :: temperatureRestoreMonthly
real (kind=8), dimension(:,:), intent(in) :: salinityRestoreMonthly
real (kind=8), dimension(:), intent(in) :: hZLevel
- real (kind=8), dimension(:), intent(in) :: referenceBottomDepth
+ real (kind=8), dimension(:), intent(in) :: refBottomDepth
integer :: nferr
@@ -609,7 +602,7 @@
start1(1) = 1
count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDh_s, start1, count1, h_s)
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDbottomDepth, start1, count1, bottomDepth)
start2(2) = 1
count2( 1) = wrLocalnVertLevels
@@ -662,7 +655,7 @@
start1(1) = 1
count1( 1) = wrLocalnVertLevels
- nferr = nf_put_vara_double(wr_ncid, wrVarIDreferenceBottomDepth, start1, count1, referenceBottomDepth)
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDrefBottomDepth, start1, count1, refBottomDepth)
start3(3) = time
count3( 1) = wrLocalnVertLevels
@@ -672,12 +665,6 @@
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)
</font>
</pre>