[Dart-dev] [5916] DART/branches/development/models/POP/model_mod.f90: Most important change: state vector actually contains potential temperature,
nancy at ucar.edu
nancy at ucar.edu
Thu Nov 1 15:01:38 MDT 2012
Revision: 5916
Author: nancy
Date: 2012-11-01 15:01:37 -0600 (Thu, 01 Nov 2012)
Log Message:
-----------
Most important change: state vector actually contains potential temperature,
not in-situ (sensible) temperature. all previous occurrences of KIND_TEMPERATURE
were changed to KIND_POTENTIAL_TEMPERATURE. a new case was added in the
model_interpolate() code to call a new routine if an observation of
KIND_TEMPERATURE was requested. the code that used to be at the end of
the model_interpolate (2 horizontal interpolations and 1 vertical) is
now in a do_interp() routine so it can also be called twice for an in-situ
temp observation: 1) to get the potential temperature and 2) the salinity
at the obs location. then those values plus pressure are used to
convert potential temperature to in-situ temperature. the difference
is in the hundredths of a degree C for shallow obs; for deeper obs
(4000-5000m) the difference is still small (e.g. ~ 0.2 degree).
added a depth -> pressure converter (from POP code) which is
location-independent so can be computed up front for all the
model levels. the potential_temp to in-situ temp converter
is also from the POP code, modified to be more compliant F90
code and removing obsolete language features.
also took the time to make all the subroutine/function declarations
and formatting the same, and made the three routines 'init_conditions',
'init_time', and 'adv_1step' throw fatal errors since they do not actually
perform the intended function. these are only called from perfect_model_obs
when the namelist is set so there is no initial state vector file.
Modified Paths:
--------------
DART/branches/development/models/POP/model_mod.f90
-------------- next part --------------
Modified: DART/branches/development/models/POP/model_mod.f90
===================================================================
--- DART/branches/development/models/POP/model_mod.f90 2012-10-30 20:00:29 UTC (rev 5915)
+++ DART/branches/development/models/POP/model_mod.f90 2012-11-01 21:01:37 UTC (rev 5916)
@@ -29,11 +29,11 @@
nc_check, do_output, to_upper, &
find_namelist_in_file, check_namelist_read, &
open_file, file_exist, find_textfile_dims, &
- file_to_text
+ file_to_text, do_output
use obs_kind_mod, only : KIND_TEMPERATURE, KIND_SALINITY, KIND_DRY_LAND, &
KIND_U_CURRENT_COMPONENT, &
KIND_V_CURRENT_COMPONENT, KIND_SEA_SURFACE_HEIGHT, &
- KIND_SEA_SURFACE_PRESSURE
+ KIND_SEA_SURFACE_PRESSURE, KIND_POTENTIAL_TEMPERATURE
use mpi_utilities_mod, only: my_task_id
use random_seq_mod, only: random_seq_type, init_random_seq, random_gaussian
use dart_pop_mod, only: set_model_time_step, &
@@ -157,6 +157,9 @@
! real, depth of lowest valid cell (0 = land). use only if KMT/KMU not avail.
real(r8), allocatable :: HT(:,:), HU(:,:)
+! compute pressure based on depth - can do once upfront.
+real(r8), allocatable :: pressure(:)
+
real(r8) :: endTime
real(r8) :: ocean_dynamics_timestep = 900.0_r4
integer :: timestepcount = 0
@@ -216,13 +219,11 @@
contains
-!==================================================================
+!------------------------------------------------------------------
+!------------------------------------------------------------------
+subroutine static_init_model()
-
-subroutine static_init_model()
-!------------------------------------------------------------------
-!
! Called to do one time initialization of the model. In this case,
! it reads in the grid information.
@@ -262,7 +263,7 @@
if (do_output()) write(logfileunit, nml=model_nml)
if (do_output()) write( * , nml=model_nml)
-!---------------------------------------------------------------
+
! Set the time step ... causes POP namelists to be read.
! Ensures model_timestep is multiple of 'ocean_dynamics_timestep'
@@ -273,7 +274,7 @@
write(msgstring,*)'assimilation period is ',dd,' days ',ss,' seconds'
call error_handler(E_MSG,'static_init_model',msgstring,source,revision,revdate)
-!---------------------------------------------------------------
+
! get data dimensions, then allocate space, then open the files
! and actually fill in the arrays.
@@ -286,6 +287,7 @@
allocate( HT(Nx,Ny), HU(Nx,Ny))
allocate( ZC(Nz), ZG(Nz))
+
! Fill them in.
! horiz grid initializes ULAT/LON, TLAT/LON as well.
! kmt initializes HT/HU if present in input file.
@@ -293,10 +295,10 @@
call read_topography(Nx, Ny, KMT, KMU)
call read_vert_grid( Nz, ZC, ZG)
-if (debug > 0) call write_grid_netcdf() ! DEBUG only
-if (debug > 0) call write_grid_interptest() ! DEBUG only
+if (debug > 2) call write_grid_netcdf() ! DEBUG only
+if (debug > 2) call write_grid_interptest() ! DEBUG only
-!---------------------------------------------------------------
+
! compute the offsets into the state vector for the start of each
! different variable type.
@@ -322,13 +324,15 @@
model_size = (n3dfields * (Nx * Ny * Nz)) + (n2dfields * (Nx * Ny))
if (do_output()) write(*,*) 'model_size = ', model_size
+! initialize the pressure array - pressure in bars
+allocate(pressure(Nz))
+call dpth2pres(Nz, ZC, pressure)
+
! Initialize the interpolation routines
call init_interp()
end subroutine static_init_model
-
-
!------------------------------------------------------------
subroutine init_interp()
@@ -355,10 +359,8 @@
end subroutine init_interp
-
!------------------------------------------------------------
-
subroutine init_dipole_interp()
! Build the data structure for interpolation for a dipole grid.
@@ -489,13 +491,12 @@
!------------------------------------------------------------
subroutine get_reg_box_indices(lon, lat, x_ind, y_ind)
+ real(r8), intent(in) :: lon, lat
+ integer, intent(out) :: x_ind, y_ind
! Given a longitude and latitude in degrees returns the index of the regular
! lon-lat box that contains the point.
-real(r8), intent(in) :: lon, lat
-integer, intent(out) :: x_ind, y_ind
-
call get_reg_lon_box(lon, x_ind)
call get_reg_lat_box(lat, y_ind)
@@ -504,12 +505,11 @@
!------------------------------------------------------------
subroutine get_reg_lon_box(lon, x_ind)
+ real(r8), intent(in) :: lon
+ integer, intent(out) :: x_ind
! Determine which regular longitude box a longitude is in.
-real(r8), intent(in) :: lon
-integer, intent(out) :: x_ind
-
x_ind = int(num_reg_x * lon / 360.0_r8) + 1
! Watch out for exactly at top; assume all lats and lons in legal range
@@ -520,12 +520,11 @@
!------------------------------------------------------------
subroutine get_reg_lat_box(lat, y_ind)
+ real(r8), intent(in) :: lat
+ integer, intent(out) :: y_ind
! Determine which regular latitude box a latitude is in.
-real(r8), intent(in) :: lat
-integer, intent(out) :: y_ind
-
y_ind = int(num_reg_y * (lat + 90.0_r8) / 180.0_r8) + 1
! Watch out for exactly at top; assume all lats and lons in legal range
@@ -536,11 +535,10 @@
!------------------------------------------------------------
subroutine reg_box_overlap(x_corners, y_corners, is_pole, reg_lon_ind, reg_lat_ind)
+ real(r8), intent(in) :: x_corners(4), y_corners(4)
+ logical, intent(in) :: is_pole
+ integer, intent(out) :: reg_lon_ind(2), reg_lat_ind(2)
-real(r8), intent(in) :: x_corners(4), y_corners(4)
-logical, intent(in) :: is_pole
-integer, intent(out) :: reg_lon_ind(2), reg_lat_ind(2)
-
! Find a set of regular lat lon boxes that covers all of the area covered by
! a dipole grid qaud whose corners are given by the dimension four x_corners
! and y_corners arrays. The two dimensional arrays reg_lon_ind and reg_lat_ind
@@ -609,11 +607,10 @@
!------------------------------------------------------------
subroutine get_quad_corners(x, i, j, corners)
+ real(r8), intent(in) :: x(:, :)
+ integer, intent(in) :: i, j
+ real(r8), intent(out) :: corners(4)
-real(r8), intent(in) :: x(:, :)
-integer, intent(in) :: i, j
-real(r8), intent(out) :: corners(4)
-
! Grabs the corners for a given quadrilateral from the global array of lower
! right corners. Note that corners go counterclockwise around the quad.
@@ -633,11 +630,11 @@
!------------------------------------------------------------
subroutine update_reg_list(reg_list_num, reg_list_lon, reg_list_lat, &
- reg_lon_ind, reg_lat_ind, dipole_lon_index, dipole_lat_index)
+ reg_lon_ind, reg_lat_ind, dipole_lon_index, dipole_lat_index)
-integer, intent(inout) :: reg_list_num(:, :), reg_list_lon(:, :, :), reg_list_lat(:, :, :)
-integer, intent(inout) :: reg_lon_ind(2), reg_lat_ind(2)
-integer, intent(in) :: dipole_lon_index, dipole_lat_index
+ integer, intent(inout) :: reg_list_num(:, :), reg_list_lon(:, :, :), reg_list_lat(:, :, :)
+ integer, intent(inout) :: reg_lon_ind(2), reg_lat_ind(2)
+ integer, intent(in) :: dipole_lon_index, dipole_lat_index
! Updates the data structure listing dipole quads that are in a given regular box
integer :: ind_x, index_x, ind_y
@@ -670,68 +667,53 @@
!------------------------------------------------------------
+subroutine init_conditions(x)
+ real(r8), intent(out) :: x(:)
-
-
-
-
-subroutine init_conditions(x)
-!------------------------------------------------------------------
-!
! Returns a model state vector, x, that is some sort of appropriate
! initial condition for starting up a long integration of the model.
! At present, this is only used if the namelist parameter
! start_from_restart is set to .false. in the program perfect_model_obs.
-! If this option is not to be used in perfect_model_obs, or if no
-! synthetic data experiments using perfect_model_obs are planned,
-! this can be a NULL INTERFACE.
-real(r8), intent(out) :: x(:)
+character(len=128) :: msgstring2, msgstring3
-if ( .not. module_initialized ) call static_init_model
-
+msgstring2 = "cannot run perfect_model_obs with 'start_from_restart = .false.' "
+msgstring3 = 'use pop_to_dart to generate an initial state'
+call error_handler(E_ERR,'init_conditions', &
+ 'WARNING!! POP model has no built-in default state', &
+ source, revision, revdate, &
+ text2=msgstring2, text3=msgstring3)
+
+! this code never reached - just here to avoid compiler warnings
+! about an intent(out) variable not being set to a value.
x = 0.0_r8
end subroutine init_conditions
+!------------------------------------------------------------------
-
subroutine adv_1step(x, time)
-!------------------------------------------------------------------
-!
-! Does a single timestep advance of the model. The input value of
-! the vector x is the starting condition and x is updated to reflect
-! the changed state after a timestep. The time argument is intent
-! in and is used for models that need to know the date/time to
-! compute a timestep, for instance for radiation computations.
-! This interface is only called IF the namelist parameter
-! async is set to 0 in perfect_model_obs or filter -OR- if the
-! program integrate_model is to be used to advance the model
-! state as a separate executable. If none of these options
-! are used (the model will only be advanced as a separate
-! model-specific executable), this can be a NULL INTERFACE.
+ real(r8), intent(inout) :: x(:)
+ type(time_type), intent(in) :: time
-real(r8), intent(inout) :: x(:)
-type(time_type), intent(in) :: time
+! If the model could be called as a subroutine, does a single
+! timestep advance. POP cannot be called this way, so fatal error
+! if this routine is called.
-if ( .not. module_initialized ) call static_init_model
+call error_handler(E_ERR,'adv_1step', &
+ 'POP model cannot be called as a subroutine; async cannot = 0', &
+ source, revision, revdate)
-if (do_output()) then
- call print_time(time,'NULL interface adv_1step (no advance) DART time is')
- call print_time(time,'NULL interface adv_1step (no advance) DART time is',logfileunit)
-endif
-
end subroutine adv_1step
+!------------------------------------------------------------------
+function get_model_size()
+ integer :: get_model_size
-function get_model_size()
-!------------------------------------------------------------------
-!
! Returns the size of the model as an integer. Required for all
! applications.
-integer :: get_model_size
if ( .not. module_initialized ) call static_init_model
@@ -739,40 +721,40 @@
end function get_model_size
+!------------------------------------------------------------------
+subroutine init_time(time)
+ type(time_type), intent(out) :: time
-subroutine init_time(time)
-!------------------------------------------------------------------
-!
-! Companion interface to init_conditions. Returns a time that is somehow
+! Companion interface to init_conditions. Returns a time that is
! appropriate for starting up a long integration of the model.
! At present, this is only used if the namelist parameter
! start_from_restart is set to .false. in the program perfect_model_obs.
-! If this option is not to be used in perfect_model_obs, or if no
-! synthetic data experiments using perfect_model_obs are planned,
-! this can be a NULL INTERFACE.
-type(time_type), intent(out) :: time
+character(len=128) :: msgstring2, msgstring3
-if ( .not. module_initialized ) call static_init_model
+msgstring2 = "cannot run perfect_model_obs with 'start_from_restart = .false.' "
+msgstring3 = 'use pop_to_dart to generate an initial state which contains a timestamp'
+call error_handler(E_ERR,'init_time', &
+ 'WARNING!! POP model has no built-in default time', &
+ source, revision, revdate, &
+ text2=msgstring2, text3=msgstring3)
-! for now, just set to 0
+! this code never reached - just here to avoid compiler warnings
+! about an intent(out) variable not being set to a value.
time = set_time(0,0)
end subroutine init_time
+!------------------------------------------------------------------
-
subroutine model_interpolate(x, location, obs_type, interp_val, istatus)
-!=======================================================================
-!
+ real(r8), intent(in) :: x(:)
+ type(location_type), intent(in) :: location
+ integer, intent(in) :: obs_type
+ real(r8), intent(out) :: interp_val
+ integer, intent(out) :: istatus
-real(r8), intent(in) :: x(:)
-type(location_type), intent(in) :: location
-integer, intent(in) :: obs_type
-real(r8), intent(out) :: interp_val
-integer, intent(out) :: istatus
-
! Model interpolate will interpolate any state variable (S, T, U, V, PSURF) to
! the given location given a state vector. The type of the variable being
! interpolated is obs_type since normally this is used to find the expected
@@ -790,6 +772,9 @@
if ( .not. module_initialized ) call static_init_model
+! print data min/max values
+if (debug > 2) call print_ranges(x)
+
! Let's assume failure. Set return val to missing, then the code can
! just set istatus to something indicating why it failed, and return.
! If the interpolation is good, the interp_val will be set to the
@@ -805,7 +790,7 @@
llat = loc_array(2)
lheight = loc_array(3)
-if (debug > 1) print *, 'requesting interpolation at ', llon, llat, lheight
+if (debug > 1) print *, 'requesting interpolation of ', obs_type, ' at ', llon, llat, lheight
if( vert_is_height(location) ) then
! Nothing to do
@@ -825,23 +810,37 @@
return
endif
-! Do horizontal interpolations for the appropriate levels
-! Find the basic offset of this field
+! kind (in-situ) temperature is a combination of potential temp,
+! salinity, and pressure based on depth. call a routine that
+! interpolates all three, does the conversion, and returns the
+! sensible/in-situ temperature.
+if(obs_type == KIND_TEMPERATURE) then
+ ! we know how to interpolate this from potential temp,
+ ! salinity, and pressure based on depth.
+ call compute_temperature(x, llon, llat, lheight, interp_val, istatus)
+ if (debug > 1) print *, 'interp val, istatus = ', interp_val, istatus
+ return
+endif
+
+! The following kinds are either in the state vector (so you
+! can simply interpolate to find the value) or they are a simple
+! transformation of something in the state vector.
+
convert_to_ssh = .FALSE.
if(obs_type == KIND_SALINITY) then
- base_offset = start_index(1)
-else if(obs_type == KIND_TEMPERATURE) then
- base_offset = start_index(2)
+ base_offset = start_index(S_index)
+else if(obs_type == KIND_POTENTIAL_TEMPERATURE) then
+ base_offset = start_index(T_index)
else if(obs_type == KIND_U_CURRENT_COMPONENT) then
- base_offset = start_index(3)
+ base_offset = start_index(U_index)
else if(obs_type == KIND_V_CURRENT_COMPONENT) then
- base_offset = start_index(4)
+ base_offset = start_index(V_index)
else if(obs_type == KIND_SEA_SURFACE_PRESSURE) then
- base_offset = start_index(5)
+ base_offset = start_index(PSURF_index)
else if(obs_type == KIND_SEA_SURFACE_HEIGHT) then
- base_offset = start_index(5) ! simple linear transform of PSURF
+ base_offset = start_index(PSURF_index) ! simple linear transform of PSURF
convert_to_ssh = .TRUE.
else
! Not a legal type for interpolation, return istatus error
@@ -849,54 +848,34 @@
return
endif
-if (debug > 1) print *, 'base offset now ', base_offset
-
-! For Sea Surface Height,Pressure don't need the vertical coordinate
+! For Sea Surface Height or Pressure don't need the vertical coordinate
! SSP needs to be converted to a SSH if height is required.
if( vert_is_surface(location) ) then
call lon_lat_interpolate(x(base_offset:), llon, llat, obs_type, 1, interp_val, istatus)
if (convert_to_ssh .and. (istatus == 0)) then
interp_val = interp_val / 980.6_r8 ! POP uses CGS units
endif
+ if (debug > 1) print *, 'interp val, istatus = ', interp_val, istatus
return
endif
! Get the bounding vertical levels and the fraction between bottom and top
-call height_bounds(lheight, nz, zc, hgt_bot, hgt_top, hgt_fract, hstatus)
+call height_bounds(lheight, Nz, ZC, hgt_bot, hgt_top, hgt_fract, hstatus)
if(hstatus /= 0) then
istatus = 12
return
endif
-! Find the base location for the bottom height and interpolate horizontally
-! on this level. Do bottom first in case it is below the ocean floor; can
-! avoid the second horizontal interpolation.
-offset = base_offset + (hgt_bot - 1) * nx * ny
-if (debug > 1) print *, 'relative bot height offset = ', offset - base_offset
-if (debug > 1) print *, 'absolute bot height offset = ', offset
-call lon_lat_interpolate(x(offset:), llon, llat, obs_type, hgt_bot, bot_val, istatus)
-! Failed istatus from interpolate means give up
-if(istatus /= 0) return
+! do a 2d interpolation for the value at the bottom level, then again for
+! the top level, then do a linear interpolation in the vertical to get the
+! final value. this sets both interp_val and istatus.
+call do_interp(x, base_offset, hgt_bot, hgt_top, hgt_fract, &
+ llon, llat, obs_type, interp_val, istatus)
+if (debug > 1) print *, 'interp val, istatus = ', interp_val, istatus
-! Find the base location for the top height and interpolate horizontally
-! on this level.
-offset = base_offset + (hgt_top - 1) * nx * ny
-if (debug > 1) print *, 'relative top height offset = ', offset - base_offset
-if (debug > 1) print *, 'absolute top height offset = ', offset
-call lon_lat_interpolate(x(offset:), llon, llat, obs_type, hgt_top, top_val, istatus)
-! Failed istatus from interpolate means give up
-if(istatus /= 0) return
-
-
-! Then weight them by the vertical fraction and return
-interp_val = bot_val + hgt_fract * (top_val - bot_val)
-if (debug > 1) print *, 'model_interp: interp val = ',interp_val
-
-! All good.
-istatus = 0
-
end subroutine model_interpolate
+!------------------------------------------------------------------
! Three different types of grids are used here. The POP dipole
! grid is referred to as a dipole grid and each region is
@@ -916,10 +895,14 @@
! The irregular grid is also assumed to be global east
! west for all applications.
+!------------------------------------------------------------------
subroutine lon_lat_interpolate(x, lon, lat, var_type, height, interp_val, istatus)
-!=======================================================================
-!
+ real(r8), intent(in) :: x(:)
+ real(r8), intent(in) :: lon, lat
+ integer, intent(in) :: var_type, height
+ real(r8), intent(out) :: interp_val
+ integer, intent(out) :: istatus
! Subroutine to interpolate to a lon lat location given the state vector
! for that level, x. This works just on one horizontal slice.
@@ -929,12 +912,6 @@
! This routine works for either the dipole or a regular lat-lon grid.
! Successful interpolation returns istatus=0.
-real(r8), intent(in) :: x(:)
-real(r8), intent(in) :: lon, lat
-integer, intent(in) :: var_type, height
-real(r8), intent(out) :: interp_val
-integer, intent(out) :: istatus
-
! Local storage
integer :: lat_bot, lat_top, lon_bot, lon_top, num_inds, start_ind
integer :: x_ind, y_ind
@@ -1073,19 +1050,16 @@
!------------------------------------------------------------
+function get_val(lon_index, lat_index, nlon, x, var_type, height, masked)
+ integer, intent(in) :: lon_index, lat_index, nlon, var_type, height
+ real(r8), intent(in) :: x(:)
+ logical, intent(out) :: masked
+ real(r8) :: get_val
-function get_val(lon_index, lat_index, nlon, x, var_type, height, masked)
-!=======================================================================
-!
! Returns the value from a single level array given the lat and lon indices
! 'masked' returns true if this is NOT a valid grid location (e.g. land, or
! below the ocean floor in shallower areas).
-integer, intent(in) :: lon_index, lat_index, nlon, var_type, height
-real(r8), intent(in) :: x(:)
-logical, intent(out) :: masked
-real(r8) :: get_val
-
if ( .not. module_initialized ) call static_init_model
! check the land/ocean bottom map and return if not valid water cell.
@@ -1103,22 +1077,21 @@
end function get_val
-
!------------------------------------------------------------
subroutine get_irreg_box(lon, lat, lon_array, lat_array, &
- found_x, found_y, lon_fract, lat_fract, istatus)
-!
+ found_x, found_y, lon_fract, lat_fract, istatus)
+
+ real(r8), intent(in) :: lon, lat
+ real(r8), intent(in) :: lon_array(nx, ny), lat_array(nx, ny)
+ real(r8), intent(out) :: lon_fract, lat_fract
+ integer, intent(out) :: found_x, found_y, istatus
+
! Given a longitude and latitude of a point (lon and lat) and the
! longitudes and latitudes of the lower left corner of the regular grid
! boxes, gets the indices of the grid box that contains the point and
! the fractions along each directrion for interpolation.
-real(r8), intent(in) :: lon, lat
-real(r8), intent(in) :: lon_array(nx, ny), lat_array(nx, ny)
-real(r8), intent(out) :: lon_fract, lat_fract
-integer, intent(out) :: found_x, found_y, istatus
-
! Local storage
integer :: lat_status, lon_top, lat_top
@@ -1141,12 +1114,13 @@
!------------------------------------------------------------
-
subroutine lon_bounds(lon, nlons, lon_array, bot, top, fract)
+ real(r8), intent(in) :: lon
+ integer, intent(in) :: nlons
+ real(r8), intent(in) :: lon_array(:, :)
+ integer, intent(out) :: bot, top
+ real(r8), intent(out) :: fract
-!=======================================================================
-!
-
! Given a longitude lon, the array of longitudes for grid boundaries, and the
! number of longitudes in the grid, returns the indices of the longitude
! below and above the location longitude and the fraction of the distance
@@ -1154,12 +1128,6 @@
! Since longitude grids are going to be regularly spaced, this could be made more efficient.
! Algorithm fails for a silly grid that has only two longitudes separated by 180 degrees.
-real(r8), intent(in) :: lon
-integer, intent(in) :: nlons
-real(r8), intent(in) :: lon_array(:, :)
-integer, intent(out) :: bot, top
-real(r8), intent(out) :: fract
-
! Local storage
integer :: i
real(r8) :: dist_bot, dist_top
@@ -1188,14 +1156,16 @@
end subroutine lon_bounds
-
!-------------------------------------------------------------
subroutine lat_bounds(lat, nlats, lat_array, bot, top, fract, istatus)
+ real(r8), intent(in) :: lat
+ integer, intent(in) :: nlats
+ real(r8), intent(in) :: lat_array(:, :)
+ integer, intent(out) :: bot, top
+ real(r8), intent(out) :: fract
+ integer, intent(out) :: istatus
-!=======================================================================
-!
-
! Given a latitude lat, the array of latitudes for grid boundaries, and the
! number of latitudes in the grid, returns the indices of the latitude
! below and above the location latitude and the fraction of the distance
@@ -1204,13 +1174,6 @@
! northernmost (2 returned). If one really had lots of polar obs would
! want to worry about interpolating around poles.
-real(r8), intent(in) :: lat
-integer, intent(in) :: nlats
-real(r8), intent(in) :: lat_array(:, :)
-integer, intent(out) :: bot, top
-real(r8), intent(out) :: fract
-integer, intent(out) :: istatus
-
! Local storage
integer :: i
@@ -1243,19 +1206,16 @@
end subroutine lat_bounds
+!------------------------------------------------------------------
-
function lon_dist(lon1, lon2)
-!=======================================================================
-!
+ real(r8), intent(in) :: lon1, lon2
+ real(r8) :: lon_dist
! Returns the smallest signed distance between lon1 and lon2 on the sphere
! If lon1 is less than 180 degrees east of lon2 the distance is negative
! If lon1 is less than 180 degrees west of lon2 the distance is positive
-real(r8), intent(in) :: lon1, lon2
-real(r8) :: lon_dist
-
if ( .not. module_initialized ) call static_init_model
lon_dist = lon2 - lon1
@@ -1269,16 +1229,14 @@
end function lon_dist
-
!------------------------------------------------------------
-
subroutine get_dipole_quad(lon, lat, qlons, qlats, num_inds, start_ind, &
- x_inds, y_inds, found_x, found_y, istatus)
+ x_inds, y_inds, found_x, found_y, istatus)
-real(r8), intent(in) :: lon, lat, qlons(:, :), qlats(:, :)
-integer, intent(in) :: num_inds, start_ind, x_inds(:), y_inds(:)
-integer, intent(out) :: found_x, found_y, istatus
+ real(r8), intent(in) :: lon, lat, qlons(:, :), qlats(:, :)
+ integer, intent(in) :: num_inds, start_ind, x_inds(:), y_inds(:)
+ integer, intent(out) :: found_x, found_y, istatus
! Given the lon and lat of a point, and a list of the
! indices of the quads that might contain a point at (lon, lat), determines
@@ -1309,13 +1267,12 @@
end subroutine get_dipole_quad
+!------------------------------------------------------------
-!------------------------------------------------------------
function in_quad(lon, lat, x_corners, y_corners)
+ real(r8), intent(in) :: lon, lat, x_corners(4), y_corners(4)
+ logical :: in_quad
-real(r8), intent(in) :: lon, lat, x_corners(4), y_corners(4)
-logical :: in_quad
-
! Return in_quad true if the point (lon, lat) is in the quad with
! the given corners.
@@ -1371,15 +1328,14 @@
end function in_quad
-
!------------------------------------------------------------
subroutine line_intercept(side_x_in, side_y, x_point_in, y_point, &
- cant_be_in_box, in_box, intercept_above, intercept_below)
+ cant_be_in_box, in_box, intercept_above, intercept_below)
-real(r8), intent(in) :: side_x_in(2), side_y(2), x_point_in, y_point
-logical, intent(out) :: cant_be_in_box, in_box
-integer, intent(out) :: intercept_above, intercept_below
+ real(r8), intent(in) :: side_x_in(2), side_y(2), x_point_in, y_point
+ logical, intent(out) :: cant_be_in_box, in_box
+ integer, intent(out) :: intercept_above, intercept_below
! Find the intercept of a vertical line from point (x_point, y_point) and
! a line segment with endpoints side_x and side_y.
@@ -1466,10 +1422,10 @@
!------------------------------------------------------------
subroutine quad_bilinear_interp(lon_in, lat, x_corners_in, y_corners, &
- p, interp_val)
+ p, interp_val)
-real(r8), intent(in) :: lon_in, lat, x_corners_in(4), y_corners(4), p(4)
-real(r8), intent(out) :: interp_val
+ real(r8), intent(in) :: lon_in, lat, x_corners_in(4), y_corners(4), p(4)
+ real(r8), intent(out) :: interp_val
! Given a longitude and latitude (lon_in, lat), the longitude and
! latitude of the 4 corners of a quadrilateral and the values at the
@@ -1566,10 +1522,9 @@
!------------------------------------------------------------
subroutine mat3x3(m, v, r)
+ real(r8), intent(in) :: m(3, 3), v(3)
+ real(r8), intent(out) :: r(3)
-real(r8), intent(in) :: m(3, 3), v(3)
-real(r8), intent(out) :: r(3)
-
! Solves rank 3 linear system mr = v for r
! using Cramer's rule. This isn't the best choice
! for speed or numerical stability so might want to replace
@@ -1592,11 +1547,11 @@
end subroutine mat3x3
!------------------------------------------------------------
+
function deter3(m)
+ real(r8), intent(in) :: m(3, 3)
+ real(r8) :: deter3
-real(r8) :: deter3
-real(r8), intent(in) :: m(3, 3)
-
! Computes determinant of 3x3 matrix m
deter3 = m(1,1)*m(2,2)*m(3,3) + m(1,2)*m(2,3)*m(3,1) + &
@@ -1604,19 +1559,16 @@
m(1,1)*m(2,3)*m(3,2) - m(3,3)*m(2,1)*m(1,2)
end function deter3
+
!------------------------------------------------------------
-
-
subroutine height_bounds(lheight, nheights, hgt_array, bot, top, fract, istatus)
-!=======================================================================
-!
-real(r8), intent(in) :: lheight
-integer, intent(in) :: nheights
-real(r8), intent(in) :: hgt_array(nheights)
-integer, intent(out) :: bot, top
-real(r8), intent(out) :: fract
-integer, intent(out) :: istatus
+ real(r8), intent(in) :: lheight
+ integer, intent(in) :: nheights
+ real(r8), intent(in) :: hgt_array(nheights)
+ integer, intent(out) :: bot, top
+ real(r8), intent(out) :: fract
+ integer, intent(out) :: istatus
! Local variables
integer :: i
@@ -1641,8 +1593,8 @@
bot = 2
! NOTE: the fract definition is the relative distance from bottom to top
fract = 1.0_r8
-if (debug > 1) print *, 'above first level in height'
-if (debug > 1) print *, 'hgt_array, top, bot, fract=', hgt_array(1), top, bot, fract
+if (debug > 7) print *, 'above first level in height'
+if (debug > 7) print *, 'hgt_array, top, bot, fract=', hgt_array(1), top, bot, fract
return
endif
@@ -1653,7 +1605,7 @@
top = i -1
bot = i
fract = (hgt_array(bot) - lheight) / (hgt_array(bot) - hgt_array(top))
-if (debug > 1) print *, 'i, hgt_array, top, bot, fract=', i, hgt_array(i), top, bot, fract
+if (debug > 7) print *, 'i, hgt_array, top, bot, fract=', i, hgt_array(i), top, bot, fract
return
endif
enddo
@@ -1663,28 +1615,28 @@
end subroutine height_bounds
+!------------------------------------------------------------------
+function get_model_time_step()
+ type(time_type) :: get_model_time_step
-function get_model_time_step()
-!------------------------------------------------------------------
-!
! Returns the the time step of the model; the smallest increment
! in time that the model is capable of advancing the state in a given
! implementation. This interface is required for all applications.
-type(time_type) :: get_model_time_step
-
if ( .not. module_initialized ) call static_init_model
get_model_time_step = model_timestep
end function get_model_time_step
+!------------------------------------------------------------------
+subroutine get_state_meta_data(index_in, location, var_type)
+ integer, intent(in) :: index_in
+ type(location_type), intent(out) :: location
+ integer, intent(out), optional :: var_type
-subroutine get_state_meta_data(index_in, location, var_type)
-!------------------------------------------------------------------
-!
! Given an integer index into the state vector structure, returns the
! associated location. A second intent(out) optional argument kind
! can be returned if the model has more than one type of field (for
@@ -1692,10 +1644,6 @@
! required for all filter applications as it is required for computing
! the distance between observations and state variables.
-integer, intent(in) :: index_in
-type(location_type), intent(out) :: location
-integer, intent(out), optional :: var_type
-
real(r8) :: lat, lon, depth
integer :: lon_index, lat_index, depth_index, local_var
@@ -1728,17 +1676,16 @@
end subroutine get_state_meta_data
+!------------------------------------------------------------------
subroutine get_state_indices(index_in, lat_index, lon_index, depth_index, var_type)
-!------------------------------------------------------------------
-!
+ integer, intent(in) :: index_in
+ integer, intent(out) :: lat_index, lon_index, depth_index
+ integer, intent(out) :: var_type
+
! Given an integer index into the state vector structure, returns the
! associated array indices for lat, lon, and depth, as well as the type.
-integer, intent(in) :: index_in
-integer, intent(out) :: lat_index, lon_index, depth_index
-integer, intent(out) :: var_type
-
integer :: startind, offset
if ( .not. module_initialized ) call static_init_model
@@ -1760,18 +1707,16 @@
end subroutine get_state_indices
+!------------------------------------------------------------------
subroutine get_state_kind(index_in, var_type, startind, offset)
-!------------------------------------------------------------------
-!
+ integer, intent(in) :: index_in
+ integer, intent(out) :: var_type, startind, offset
+
! Given an integer index into the state vector structure, returns the kind,
! and both the starting offset for this kind, as well as the offset into
! the block of this kind.
-integer, intent(in) :: index_in
-integer, intent(out) :: var_type, startind, offset
-
-
if ( .not. module_initialized ) call static_init_model
if (debug > 5) print *, 'asking for meta data about index ', index_in
@@ -1780,7 +1725,7 @@
var_type = KIND_SALINITY
startind = start_index(S_index)
else if (index_in < start_index(T_index+1)) then
- var_type = KIND_TEMPERATURE
+ var_type = KIND_POTENTIAL_TEMPERATURE
startind = start_index(T_index)
else if (index_in < start_index(U_index+1)) then
var_type = KIND_U_CURRENT_COMPONENT
@@ -1802,17 +1747,15 @@
end subroutine get_state_kind
+!------------------------------------------------------------------
+subroutine get_state_kind_inc_dry(index_in, var_type)
+ integer, intent(in) :: index_in
+ integer, intent(out) :: var_type
-subroutine get_state_kind_inc_dry(index_in, var_type)
-!------------------------------------------------------------------
-!
! Given an integer index into the state vector structure, returns the
! type, taking into account the ocean bottom and dry land.
-integer, intent(in) :: index_in
-integer, intent(out) :: var_type
-
integer :: lon_index, lat_index, depth_index, startind, offset
if ( .not. module_initialized ) call static_init_model
@@ -1835,24 +1778,26 @@
end subroutine get_state_kind_inc_dry
+!------------------------------------------------------------------
subroutine end_model()
-!------------------------------------------------------------------
-!
-! Does any shutdown and clean-up needed for model. Can be a NULL
-! INTERFACE if the model has no need to clean up storage, etc.
-! if ( .not. module_initialized ) call static_init_model
+! Shutdown and clean-up.
-deallocate(ULAT, ULON, TLAT, TLON, KMT, KMU, HT, HU)
-deallocate(ZC, ZG)
+! assume if one is allocated, they all were. if no one ever
+! called the init routine, don't try to dealloc something that
+! was never alloc'd.
+if (allocated(ULAT)) deallocate(ULAT, ULON, TLAT, TLON, KMT, KMU, HT, HU)
+if (allocated(ZC)) deallocate(ZC, ZG, pressure)
end subroutine end_model
+!------------------------------------------------------------------
+function nc_write_model_atts( ncFileID ) result (ierr)
+ integer, intent(in) :: ncFileID ! netCDF file identifier
+ integer :: ierr ! return value of function
-function nc_write_model_atts( ncFileID ) result (ierr)
-!------------------------------------------------------------------
! TJH -- Writes the model-specific attributes to a netCDF file.
! This includes coordinate variables and some metadata, but NOT
! the model state vector.
@@ -1871,9 +1816,6 @@
! NF90_put_var ! provide values for variable
! NF90_CLOSE ! close: save updated netCDF dataset
-integer, intent(in) :: ncFileID ! netCDF file identifier
-integer :: ierr ! return value of function
-
integer :: nDimensions, nVariables, nAttributes, unlimitedDimID
!----------------------------------------------------------------------
@@ -2311,10 +2253,15 @@
end function nc_write_model_atts
+!------------------------------------------------------------------
+function nc_write_model_vars( ncFileID, statevec, copyindex, timeindex ) result (ierr)
+ integer, intent(in) :: ncFileID ! netCDF file identifier
+ real(r8), dimension(:), intent(in) :: statevec
+ integer, intent(in) :: copyindex
+ integer, intent(in) :: timeindex
+ integer :: ierr ! return value of function
-function nc_write_model_vars( ncFileID, statevec, copyindex, timeindex ) result (ierr)
-!------------------------------------------------------------------
! TJH 24 Oct 2006 -- Writes the model variables to a netCDF file.
!
! TJH 29 Jul 2003 -- for the moment, all errors are fatal, so the
@@ -2337,12 +2284,6 @@
! NF90_put_var ! provide values for variable
! NF90_CLOSE ! close: save updated netCDF dataset
-integer, intent(in) :: ncFileID ! netCDF file identifier
-real(r8), dimension(:), intent(in) :: statevec
-integer, intent(in) :: copyindex
-integer, intent(in) :: timeindex
-integer :: ierr ! return value of function
-
integer :: nDimensions, nVariables, nAttributes, unlimitedDimID
integer :: VarID
@@ -2432,11 +2373,13 @@
end function nc_write_model_vars
+!------------------------------------------------------------------
+subroutine pert_model_state(state, pert_state, interf_provided)
+ real(r8), intent(in) :: state(:)
+ real(r8), intent(out) :: pert_state(:)
+ logical, intent(out) :: interf_provided
-subroutine pert_model_state(state, pert_state, interf_provided)
-!------------------------------------------------------------------
-!
! Perturbs a model state for generating initial ensembles.
! The perturbed state is returned in pert_state.
! A model may choose to provide a NULL INTERFACE by returning
@@ -2447,10 +2390,6 @@
! should be returned as .true. if the model wants to do its own
! perturbing of states.
-real(r8), intent(in) :: state(:)
-real(r8), intent(out) :: pert_state(:)
-logical, intent(out) :: interf_provided
-
integer :: i, var_type
logical, save :: random_seq_init = .false.
@@ -2479,34 +2418,75 @@
end subroutine pert_model_state
+!------------------------------------------------------------------
+subroutine ens_mean_for_model(ens_mean)
+ real(r8), intent(in) :: ens_mean(:)
-
-subroutine ens_mean_for_model(ens_mean)
-!------------------------------------------------------------------
! If needed by the model interface, this is the current mean
! for all state vector items across all ensembles. It is up to this
! code to allocate space and save a copy if it is going to be used
! later on. For now, we are ignoring it.
-real(r8), intent(in) :: ens_mean(:)
-
if ( .not. module_initialized ) call static_init_model
+
end subroutine ens_mean_for_model
+!------------------------------------------------------------------
+subroutine print_ranges(x)
+ real(r8), intent(in) :: x(:)
+! intended for debugging use = print out the data min/max for each
+! field in the state vector, along with the starting and ending
+! indices for each field.
+integer :: s, e
+
+! S_index = 1
+! T_index = 2
+! U_index = 3
+! V_index = 4
+! PSURF_index = 5
+
+s = 1
+e = start_index(T_index)-1
+print *, 'min/max salinity: ', minval(x(s:e)), maxval(x(s:e))
+print *, 's/e salinity: ', s, e
+
+s = start_index(T_index)
+e = start_index(U_index)-1
+print *, 'min/max pot temp: ', minval(x(s:e)), maxval(x(s:e))
+print *, 's/e pot temp: ', s, e
+
+s = start_index(U_index)
+e = start_index(V_index)-1
+print *, 'min/max U current: ', minval(x(s:e)), maxval(x(s:e))
+print *, 's/e U current: ', s, e
+
+s = start_index(V_index)
+e = start_index(PSURF_index)-1
+print *, 'min/max V current: ', minval(x(s:e)), maxval(x(s:e))
+print *, 's/e V current: ', s, e
+
+s = start_index(PSURF_index)
+e = size(x)
+print *, 'min/max Surf Pres: ', minval(x(s:e)), maxval(x(s:e))
+print *, 's/e Surf Pres: ', s, e
+
+end subroutine print_ranges
+
+!------------------------------------------------------------------
+
subroutine restart_file_to_sv(filename, state_vector, model_time)
-!------------------------------------------------------------------
+ character(len=*), intent(in) :: filename
+ real(r8), intent(inout) :: state_vector(:)
+ type(time_type), intent(out) :: model_time
+
! Reads the current time and state variables from a POP restart
! file and packs them into a dart state vector.
-character(len=*), intent(in) :: filename
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list