[Dart-dev] [6897] DART/trunk/models: These were a test for a fully-coupled assimilation with three model_mod() modules.
nancy at ucar.edu
nancy at ucar.edu
Thu Apr 17 15:08:02 MDT 2014
Revision: 6897
Author: thoar
Date: 2014-04-17 15:08:01 -0600 (Thu, 17 Apr 2014)
Log Message:
-----------
These were a test for a fully-coupled assimilation with three model_mod() modules.
The test is over, and these should not be in common circulation.
Removed Paths:
-------------
DART/trunk/models/POP/pop_model_mod.f90
DART/trunk/models/cam/cam_model_mod.f90
DART/trunk/models/clm/clm_model_mod.f90
-------------- next part --------------
Deleted: DART/trunk/models/POP/pop_model_mod.f90
===================================================================
--- DART/trunk/models/POP/pop_model_mod.f90 2014-04-17 21:05:44 UTC (rev 6896)
+++ DART/trunk/models/POP/pop_model_mod.f90 2014-04-17 21:08:01 UTC (rev 6897)
@@ -1,3630 +0,0 @@
-! DART software - Copyright 2004 - 2013 UCAR. This open source software is
-! provided by UCAR, "as is", without charge, subject to all terms of use at
-! http://www.image.ucar.edu/DAReS/DART/DART_download
-!
-! $Id$
-
-module pop_model_mod
-
-! This is the interface between the POP ocean model and DART.
-
-! Modules that are absolutely required for use are listed
-use types_mod, only : r4, r8, SECPERDAY, MISSING_R8, rad2deg, PI
-use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time,&
- print_time, print_date, &
- operator(*), operator(+), operator(-), &
- operator(>), operator(<), operator(/), &
- operator(/=), operator(<=)
-use location_mod, only : location_type, get_dist, &
- get_close_maxdist_init, get_close_obs_init, &
- set_location, &
- VERTISHEIGHT, get_location, vert_is_height, &
- vert_is_level, vert_is_surface, &
- loc_get_close_obs => get_close_obs, get_close_type
-use utilities_mod, only : register_module, error_handler, &
- E_ERR, E_WARN, E_MSG, logfileunit, get_unit, &
- nc_check, do_output, to_upper, &
- find_namelist_in_file, check_namelist_read, &
- open_file, file_exist, find_textfile_dims, &
- 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_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, &
- get_horiz_grid_dims, get_vert_grid_dim, &
- read_horiz_grid, read_topography, read_vert_grid, &
- get_pop_restart_filename
-
-use typesizes
-use netcdf
-
-implicit none
-private
-
-! these routines must be public and you cannot change
-! the arguments - they will be called *from* the DART code.
-public :: pop_get_model_size, &
- pop_adv_1step, &
- pop_get_state_meta_data, &
- pop_model_interpolate, &
- pop_get_model_time_step, &
- pop_static_init_model, &
- pop_end_model, &
- pop_init_time, &
- pop_init_conditions, &
- pop_nc_write_model_atts, &
- pop_nc_write_model_vars, &
- pop_pert_model_state, &
- get_close_maxdist_init, &
- get_close_obs_init, &
- pop_get_close_obs, &
- pop_ens_mean_for_model
-
-! generally useful routines for various support purposes.
-! the interfaces here can be changed as appropriate.
-public :: pop_get_gridsize, pop_restart_file_to_sv, pop_sv_to_restart_file, &
- get_pop_restart_filename, pop_test_interpolation
-
-! version controlled file description for error handling, do not edit
-character(len=256), parameter :: source = &
- "$URL$"
-character(len=32 ), parameter :: revision = "$Revision$"
-character(len=128), parameter :: revdate = "$Date$"
-
-character(len=256) :: msgstring
-logical, save :: module_initialized = .false.
-
-! Storage for a random sequence for perturbing a single initial state
-type(random_seq_type) :: random_seq
-
-! things which can/should be in the model_nml
-logical :: output_state_vector = .true.
-integer :: assimilation_period_days = 1
-integer :: assimilation_period_seconds = 0
-real(r8) :: model_perturbation_amplitude = 0.2
-logical :: update_dry_cell_walls = .false.
-integer :: debug = 0 ! turn up for more and more debug messages
-
-! FIXME: currently the update_dry_cell_walls namelist value DOES
-! NOTHING. it needs additional code to detect the cells which are
-! wet, but within 1 cell of the bottom/sides/etc.
-
-namelist /pop_model_nml/ &
- output_state_vector, &
- assimilation_period_days, & ! for now, this is the timestep
- assimilation_period_seconds, &
- model_perturbation_amplitude,&
- update_dry_cell_walls, &
- debug
-
-!------------------------------------------------------------------
-!
-! The DART state vector (control vector) will consist of: S, T, U, V, PSURF
-! (Salinity, Temperature, U velocity, V velocity, Sea Surface Height).
-! S, T are 3D arrays, located at cell centers. U,V are at grid cell corners.
-! PSURF is a 2D field (X,Y only). The Z direction is downward.
-!
-! FIXME: proposed change 1: we put SSH first, then T,U,V, then S, then
-! any optional tracers, since SSH is the only 2D
-! field; all tracers are 3D. this simplifies the
-! mapping to and from the vars to state vector.
-!
-! FIXME: proposed change 2: we make this completely namelist driven,
-! both contents and order of vars. this should
-! wait until restart files are in netcdf format,
-! to avoid problems with incompatible namelist
-! and IC files. it also complicates the mapping
-! to and from the vars to state vector.
-!------------------------------------------------------------------
-
-integer, parameter :: n3dfields = 4
-integer, parameter :: n2dfields = 1
-integer, parameter :: nfields = n3dfields + n2dfields
-
-! (the absoft compiler likes them to all be the same length during declaration)
-! we trim the blanks off before use anyway, so ...
-character(len=128) :: progvarnames(nfields) = (/'SALT ','TEMP ','UVEL ','VVEL ','PSURF'/)
-
-integer, parameter :: S_index = 1
-integer, parameter :: T_index = 2
-integer, parameter :: U_index = 3
-integer, parameter :: V_index = 4
-integer, parameter :: PSURF_index = 5
-
-integer :: start_index(nfields)
-
-! Grid parameters - the values will be read from a
-! standard POP namelist and filled in here.
-
-! nx, ny and nz are the size of the dipole (or irregular) grids.
-integer :: Nx=-1, Ny=-1, Nz=-1 ! grid counts for each field
-
-! locations of cell centers (C) and edges (G) for each axis.
-real(r8), allocatable :: ZC(:), ZG(:)
-
-! These arrays store the longitude and latitude of the lower left corner of
-! each of the dipole u quadrilaterals and t quadrilaterals.
-real(r8), allocatable :: ULAT(:,:), ULON(:,:), TLAT(:,:), TLON(:,:)
-
-! integer, lowest valid cell number in the vertical
-integer, allocatable :: KMT(:, :), KMU(:, :)
-! 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
-type(time_type) :: model_time, model_timestep
-
-integer :: model_size ! the state vector length
-
-
-INTERFACE vector_to_prog_var
- MODULE PROCEDURE vector_to_2d_prog_var
- MODULE PROCEDURE vector_to_3d_prog_var
-END INTERFACE
-
-
-!------------------------------------------------
-
-! The regular grid used for dipole interpolation divides the sphere into
-! a set of regularly spaced lon-lat boxes. The number of boxes in
-! longitude and latitude are set by num_reg_x and num_reg_y. Making the
-! number of regular boxes smaller decreases the computation required for
-! doing each interpolation but increases the static storage requirements
-! and the initialization computation (which seems to be pretty small).
-integer, parameter :: num_reg_x = 90, num_reg_y = 90
-
-! The max_reg_list_num controls the size of temporary storage used for
-! initializing the regular grid. Four arrays
-! of size num_reg_x*num_reg_y*max_reg_list_num are needed. The initialization
-! fails and returns an error if max_reg_list_num is too small. A value of
-! 30 is sufficient for the x3 POP grid with 180 regular lon and lat boxes
-! and a value of 80 is sufficient for for the x1 grid.
-integer, parameter :: max_reg_list_num = 80
-
-! The dipole interpolation keeps a list of how many and which dipole quads
-! overlap each regular lon-lat box. The number for the u and t grids are
-! stored in u_dipole_num and t_dipole_num. The allocatable arrays
-! u_dipole_lon(lat)_list and t_dipole_lon(lat)_list list the longitude
-! and latitude indices for the overlapping dipole quads. The entry in
-! u_dipole_start and t_dipole_start for a given regular lon-lat box indicates
-! where the list of dipole quads begins in the u_dipole_lon(lat)_list and
-! t_dipole_lon(lat)_list arrays.
-
-integer :: u_dipole_start(num_reg_x, num_reg_y)
-integer :: u_dipole_num (num_reg_x, num_reg_y) = 0
-integer :: t_dipole_start(num_reg_x, num_reg_y)
-integer :: t_dipole_num (num_reg_x, num_reg_y) = 0
-integer, allocatable :: u_dipole_lon_list(:), t_dipole_lon_list(:)
-integer, allocatable :: u_dipole_lat_list(:), t_dipole_lat_list(:)
-
-! Need to check for pole quads: for now we are not interpolating in them
-integer :: pole_x, t_pole_y, u_pole_y
-
-
-
-! Have a global variable saying whether this is dipole or regular lon-lat grid
-! This should be initialized static_init_model. Code to do this is below.
-logical :: dipole_grid
-
-contains
-
-!------------------------------------------------------------------
-!------------------------------------------------------------------
-
-subroutine pop_static_init_model()
-
-! Called to do one time initialization of the model. In this case,
-! it reads in the grid information.
-
-integer :: iunit, io
-integer :: ss, dd
-
-! The Plan:
-!
-! read in the grid sizes from the horiz grid file and the vert grid file
-! horiz is netcdf, vert is ascii
-!
-! allocate space, and read in actual grid values
-!
-! figure out model timestep. FIXME: from where?
-!
-! Compute the model size.
-!
-! set the index numbers where the field types change
-!
-
-if ( module_initialized ) return ! only need to do this once.
-
-! Print module information to log file and stdout.
-call register_module(source, revision, revdate)
-
-! Since this routine calls other routines that could call this routine
-! we'll say we've been initialized pretty dang early.
-module_initialized = .true.
-
-! Read the DART namelist for this model
-call find_namelist_in_file('input.nml', 'pop_model_nml', iunit)
-read(iunit, nml = pop_model_nml, iostat = io)
-call check_namelist_read(iunit, io, 'pop_model_nml')
-
-! Record the namelist values used for the run
-call error_handler(E_MSG,'static_init_model','pop_model_nml values are',' ',' ',' ')
-if (do_output()) write(logfileunit, nml=pop_model_nml)
-if (do_output()) write( * , nml=pop_model_nml)
-
-
-! Set the time step ... causes POP namelists to be read.
-! Ensures model_timestep is multiple of 'ocean_dynamics_timestep'
-
-model_timestep = set_model_time_step()
-
-call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400)
-
-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.
-
-call get_horiz_grid_dims(Nx, Ny)
-call get_vert_grid_dim(Nz)
-
-! Allocate space for grid variables.
-allocate(ULAT(Nx,Ny), ULON(Nx,Ny), TLAT(Nx,Ny), TLON(Nx,Ny))
-allocate( KMT(Nx,Ny), KMU(Nx,Ny))
-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.
-call read_horiz_grid(Nx, Ny, ULAT, ULON, TLAT, TLON)
-call read_topography(Nx, Ny, KMT, KMU)
-call read_vert_grid( Nz, ZC, ZG)
-
-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.
-
-! record where in the state vector the data type changes
-! from one type to another, by computing the starting
-! index for each block of data.
-start_index(S_index) = 1
-start_index(T_index) = start_index(S_index) + (Nx * Ny * Nz)
-start_index(U_index) = start_index(T_index) + (Nx * Ny * Nz)
-start_index(V_index) = start_index(U_index) + (Nx * Ny * Nz)
-start_index(PSURF_index) = start_index(V_index) + (Nx * Ny * Nz)
-
-! in spite of the staggering, all grids are the same size
-! and offset by half a grid cell. 4 are 3D and 1 is 2D.
-! e.g. S,T,U,V = 256 x 225 x 70
-! e.g. PSURF = 256 x 225
-
-if (do_output()) write(logfileunit, *) 'Using grid : Nx, Ny, Nz = ', &
- Nx, Ny, Nz
-if (do_output()) write( * , *) 'Using grid : Nx, Ny, Nz = ', &
- Nx, Ny, Nz
-
-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 pop_static_init_model
-
-!------------------------------------------------------------
-
-subroutine init_interp()
-
-! Initializes data structures needed for POP interpolation for
-! either dipole or irregular grid.
-! This should be called at static_init_model time to avoid
-! having all this temporary storage in the middle of a run.
-
-integer :: i
-
-! Determine whether this is a irregular lon-lat grid or a dipole.
-! Do this by seeing if the lons have the same values at both
-! the first and last latitude row; this is not the case for dipole.
-
-dipole_grid = .false.
-do i = 1, nx
- if(ulon(i, 1) /= ulon(i, ny)) then
- dipole_grid = .true.
- call init_dipole_interp()
- return
- endif
-enddo
-
-end subroutine init_interp
-
-!------------------------------------------------------------
-
-subroutine init_dipole_interp()
-
-! Build the data structure for interpolation for a dipole grid.
-
-! Need a temporary data structure to build this.
-! These arrays keep a list of the x and y indices of dipole quads
-! that potentially overlap the regular boxes. Need one for the u
-! and one for the t grid.
-integer :: ureg_list_lon(num_reg_x, num_reg_y, max_reg_list_num)
-integer :: ureg_list_lat(num_reg_x, num_reg_y, max_reg_list_num)
-integer :: treg_list_lon(num_reg_x, num_reg_y, max_reg_list_num)
-integer :: treg_list_lat(num_reg_x, num_reg_y, max_reg_list_num)
-
-
-real(r8) :: u_c_lons(4), u_c_lats(4), t_c_lons(4), t_c_lats(4), pole_row_lon
-integer :: i, j, k, pindex
-integer :: reg_lon_ind(2), reg_lat_ind(2), u_total, t_total, u_index, t_index
-logical :: is_pole
-
-! Begin by finding the quad that contains the pole for the dipole t_grid.
-! To do this locate the u quad with the pole on its right boundary. This is on
-! the row that is opposite the shifted pole and exactly follows a lon circle.
-pole_x = nx / 2;
-! Search for the row at which the longitude flips over the pole
-pole_row_lon = ulon(pole_x, 1);
-do i = 1, ny
- pindex = i
- if(ulon(pole_x, i) /= pole_row_lon) exit
-enddo
-
-! Pole boxes for u have indices pole_x or pole_x-1 and index - 1;
-! (it's right before the flip).
-u_pole_y = pindex - 1;
-
-! Locate the T dipole quad that contains the pole.
-! We know it is in either the same lat quad as the u pole edge or one higher.
-! Figure out if the pole is more or less than halfway along
-! the u quad edge to find the right one.
-if(ulat(pole_x, u_pole_y) > ulat(pole_x, u_pole_y + 1)) then
- t_pole_y = u_pole_y;
-else
- t_pole_y = u_pole_y + 1;
-endif
-
-! Loop through each of the dipole grid quads
-do i = 1, nx
- ! There's no wraparound in y, one box less than grid boundaries
- do j = 1, ny - 1
- ! Is this the pole quad for the T grid?
- is_pole = (i == pole_x .and. j == t_pole_y)
-
- ! Set up array of lons and lats for the corners of these u and t quads
- call get_quad_corners(ulon, i, j, u_c_lons)
- call get_quad_corners(ulat, i, j, u_c_lats)
- call get_quad_corners(tlon, i, j, t_c_lons)
- call get_quad_corners(tlat, i, j, t_c_lats)
-
- ! Get list of regular boxes that cover this u dipole quad
- ! false indicates that for the u grid there's nothing special about pole
- call reg_box_overlap(u_c_lons, u_c_lats, .false., reg_lon_ind, reg_lat_ind)
-
- ! Update the temporary data structures for the u quad
- call update_reg_list(u_dipole_num, ureg_list_lon, &
- ureg_list_lat, reg_lon_ind, reg_lat_ind, i, j)
-
- ! Repeat for t dipole quads
- call reg_box_overlap(t_c_lons, t_c_lats, is_pole, reg_lon_ind, reg_lat_ind)
- call update_reg_list(t_dipole_num, treg_list_lon, &
- treg_list_lat, reg_lon_ind, reg_lat_ind, i, j)
- enddo
-enddo
-
-if (do_output()) write(*,*)'to determine (minimum) max_reg_list_num values for new grids ...'
-if (do_output()) write(*,*)'u_dipole_num is ',maxval(u_dipole_num)
-if (do_output()) write(*,*)'t_dipole_num is ',maxval(t_dipole_num)
-
-! Invert the temporary data structure. The total number of entries will be
-! the sum of the number of dipole cells for each regular cell.
-u_total = sum(u_dipole_num)
-t_total = sum(t_dipole_num)
-
-! Allocate storage for the final structures in module storage
-allocate(u_dipole_lon_list(u_total), u_dipole_lat_list(u_total))
-allocate(t_dipole_lon_list(t_total), t_dipole_lat_list(t_total))
-
-! Fill up the long list by traversing the temporary structure. Need indices
-! to keep track of where to put the next entry.
-u_index = 1
-t_index = 1
-
-! Loop through each regular grid box
-do i = 1, num_reg_x
- do j = 1, num_reg_y
-
- ! The list for this regular box starts at the current indices.
- u_dipole_start(i, j) = u_index
- t_dipole_start(i, j) = t_index
-
- ! Copy all the close dipole quads for regular u box(i, j)
- do k = 1, u_dipole_num(i, j)
- u_dipole_lon_list(u_index) = ureg_list_lon(i, j, k)
- u_dipole_lat_list(u_index) = ureg_list_lat(i, j, k)
- u_index = u_index + 1
- enddo
-
- ! Copy all the close dipoles for regular t box (i, j)
- do k = 1, t_dipole_num(i, j)
- t_dipole_lon_list(t_index) = treg_list_lon(i, j, k)
- t_dipole_lat_list(t_index) = treg_list_lat(i, j, k)
- t_index = t_index + 1
- enddo
-
- enddo
-enddo
-
-! Confirm that the indices come out okay as debug
-if(u_index /= u_total + 1) then
- msgstring = 'Storage indices did not balance for U grid: : contact DART developers'
- call error_handler(E_ERR, 'init_dipole_interp', msgstring, source, revision, revdate)
-endif
-if(t_index /= t_total + 1) then
- msgstring = 'Storage indices did not balance for T grid: : contact DART developers'
- call error_handler(E_ERR, 'init_dipole_interp', msgstring, source, revision, revdate)
-endif
-
-end subroutine init_dipole_interp
-
-!------------------------------------------------------------
-
-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.
-
-call get_reg_lon_box(lon, x_ind)
-call get_reg_lat_box(lat, y_ind)
-
-end subroutine get_reg_box_indices
-
-!------------------------------------------------------------
-
-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.
-
-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
-if(lon == 360.0_r8) x_ind = num_reg_x
-
-end subroutine get_reg_lon_box
-
-!------------------------------------------------------------
-
-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.
-
-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
-if(lat == 90.0_r8) y_ind = num_reg_y
-
-end subroutine get_reg_lat_box
-
-!------------------------------------------------------------
-
-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)
-
-! 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
-! return the first and last indices of the regular boxes in latitude and
-! longitude respectively. These indices may wraparound for reg_lon_ind.
-! A special computation is needed for a dipole quad that has the true north
-! pole in its interior. The logical is_pole is set to true if this is the case.
-! This can only happen for the t grid. If the longitude boxes overlap 0
-! degrees, the indices returned are adjusted by adding the total number of
-! boxes to the second index (e.g. the indices might be 88 and 93 for a case
-! with 90 longitude boxes).
-
-real(r8) :: lat_min, lat_max, lon_min, lon_max
-integer :: i
-
-! A quad containing the pole is fundamentally different
-if(is_pole) then
- ! Need all longitude boxes
- reg_lon_ind(1) = 1
- reg_lon_ind(2) = num_reg_x
- ! Need to cover from lowest latitude to top box
- lat_min = minval(y_corners)
- reg_lat_ind(1) = int(num_reg_y * (lat_min + 90.0_r8) / 180.0_r8) + 1
- call get_reg_lat_box(lat_min, reg_lat_ind(1))
- reg_lat_ind(2) = num_reg_y
-else
- ! All other quads do not contain pole (pole could be on edge but no problem)
- ! This is specific to the dipole POP grids that do not go to the south pole
- ! Finding the range of latitudes is cake
- lat_min = minval(y_corners)
- lat_max = maxval(y_corners)
-
- ! Figure out the indices of the regular boxes for min and max lats
- call get_reg_lat_box(lat_min, reg_lat_ind(1))
- call get_reg_lat_box(lat_max, reg_lat_ind(2))
-
- ! Lons are much trickier. Need to make sure to wraparound the
- ! right way. There is no guarantee on direction of lons in the
- ! high latitude dipole rows.
- ! All longitudes for non-pole rows have to be within 180 degrees
- ! of one another.
- lon_min = minval(x_corners)
- lon_max = maxval(x_corners)
- if((lon_max - lon_min) > 180.0_r8) then
- ! If the max longitude value is more than 180
- ! degrees larger than the min, then there must be wraparound.
- ! Then, find the smallest value > 180 and the largest < 180 to get range.
- lon_min = 360.0_r8
- lon_max = 0.0_r8
- do i=1, 4
- if(x_corners(i) > 180.0_r8 .and. x_corners(i) < lon_min) lon_min = x_corners(i)
- if(x_corners(i) < 180.0_r8 .and. x_corners(i) > lon_max) lon_max = x_corners(i)
- enddo
- endif
-
- ! Get the indices for the extreme longitudes
- call get_reg_lon_box(lon_min, reg_lon_ind(1))
- call get_reg_lon_box(lon_max, reg_lon_ind(2))
-
- ! Watch for wraparound again; make sure that second index is greater than first
- if(reg_lon_ind(2) < reg_lon_ind(1)) reg_lon_ind(2) = reg_lon_ind(2) + num_reg_x
-endif
-
-end subroutine reg_box_overlap
-
-!------------------------------------------------------------
-
-subroutine get_quad_corners(x, i, j, corners)
- 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.
-
-integer :: ip1
-
-! Have to worry about wrapping in longitude but not in latitude
-ip1 = i + 1
-if(ip1 > nx) ip1 = 1
-
-corners(1) = x(i, j )
-corners(2) = x(ip1, j )
-corners(3) = x(ip1, j+1)
-corners(4) = x(i, j+1)
-
-end subroutine get_quad_corners
-
-!------------------------------------------------------------
-
-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)
-
- 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
-
-! Loop through indices for each possible regular cell
-! Have to watch for wraparound in longitude
-if(reg_lon_ind(2) < reg_lon_ind(1)) reg_lon_ind(2) = reg_lon_ind(2) + num_reg_x
-
-do ind_x = reg_lon_ind(1), reg_lon_ind(2)
- ! Inside loop, need to go back to wraparound indices to find right box
- index_x = ind_x
- if(index_x > num_reg_x) index_x = index_x - num_reg_x
-
- do ind_y = reg_lat_ind(1), reg_lat_ind(2)
- ! Make sure the list storage isn't full
- if(reg_list_num(index_x, ind_y) >= max_reg_list_num) then
- write(msgstring,*) 'max_reg_list_num (',max_reg_list_num,') is too small ... increase'
- call error_handler(E_ERR, 'update_reg_list', msgstring, source, revision, revdate)
- endif
-
- ! Increment the count
- reg_list_num(index_x, ind_y) = reg_list_num(index_x, ind_y) + 1
- ! Store this quad in the list for this regular box
- reg_list_lon(index_x, ind_y, reg_list_num(index_x, ind_y)) = dipole_lon_index
- reg_list_lat(index_x, ind_y, reg_list_num(index_x, ind_y)) = dipole_lat_index
- enddo
-enddo
-
-end subroutine update_reg_list
-
-!------------------------------------------------------------
-
-subroutine pop_init_conditions(x)
- real(r8), intent(out) :: 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.
-
-character(len=128) :: msgstring2, msgstring3
-
-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 pop_init_conditions
-
-!------------------------------------------------------------------
-
-subroutine pop_adv_1step(x, 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.
-
-call error_handler(E_ERR,'adv_1step', &
- 'POP model cannot be called as a subroutine; async cannot = 0', &
- source, revision, revdate)
-
-end subroutine pop_adv_1step
-
-!------------------------------------------------------------------
-
-function pop_get_model_size()
- integer :: pop_get_model_size
-
-! Returns the size of the model as an integer. Required for all
-! applications.
-
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-pop_get_model_size = model_size
-
-end function pop_get_model_size
-
-!------------------------------------------------------------------
-
-subroutine pop_init_time(time)
- type(time_type), intent(out) :: time
-
-! 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.
-
-character(len=128) :: msgstring2, msgstring3
-
-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)
-
-! 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 pop_init_time
-
-!------------------------------------------------------------------
-
-subroutine pop_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
-
-! 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
-! value of an observation at some location. The interpolated value is
-! returned in interp_val and istatus is 0 for success.
-
-! Local storage
-real(r8) :: loc_array(3), llon, llat, lheight
-integer :: base_offset, offset, ind
-integer :: hgt_bot, hgt_top
-real(r8) :: hgt_fract
-real(r8) :: top_val, bot_val
-integer :: hstatus
-logical :: convert_to_ssh
-
-if ( .not. module_initialized ) call pop_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
-! good value, and the last line here sets istatus to 0.
-! make any error codes set here be in the 10s
-
-interp_val = MISSING_R8 ! the DART bad value flag
-istatus = 99 ! unknown error
-
-! Get the individual locations values
-loc_array = get_location(location)
-llon = loc_array(1)
-llat = loc_array(2)
-lheight = loc_array(3)
-
-if (debug > 1) print *, 'requesting interpolation of ', obs_type, ' at ', llon, llat, lheight
-
-if( vert_is_height(location) ) then
- ! Nothing to do
-elseif ( vert_is_surface(location) ) then
- ! Nothing to do
-elseif (vert_is_level(location)) then
- ! convert the level index to an actual depth
- ind = nint(loc_array(3))
- if ( (ind < 1) .or. (ind > size(zc)) ) then
- istatus = 11
- return
- else
- lheight = zc(ind)
- endif
-else ! if pressure or undefined, we don't know what to do
- istatus = 17
- return
-endif
-
-! 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(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(U_index)
-else if(obs_type == KIND_V_CURRENT_COMPONENT) then
- base_offset = start_index(V_index)
-else if(obs_type == KIND_SEA_SURFACE_PRESSURE) then
- base_offset = start_index(PSURF_index)
-else if(obs_type == KIND_SEA_SURFACE_HEIGHT) then
- 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
- istatus = 15
- return
-endif
-
-! 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)
-if(hstatus /= 0) then
- istatus = 12
- return
-endif
-
-! 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
-
-end subroutine pop_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
-! referred to as a quad, short for quadrilateral.
-! The longitude latitude rectangular grid with possibly irregular
-! spacing in latitude used for some POP applications and testing
-! is referred to as the irregular grid and each region is
-! called a box.
-! Finally, a regularly spaced longitude latitude grid is used
-! as a computational tool for interpolating from the dipole
-! grid. This is referred to as the regular grid and each region
-! is called a box.
-! All grids are referenced by the index of the lower left corner
-! of the quad or box.
-
-! The dipole grid is assumed to be global for all applications.
-! 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.
-! NOTE: Using array sections to pass in the x array may be inefficient on some
-! compiler/platform setups. Might want to pass in the entire array with a base
-! offset value instead of the section if this is an issue.
-! This routine works for either the dipole or a regular lat-lon grid.
-! Successful interpolation returns istatus=0.
-
-! Local storage
-integer :: lat_bot, lat_top, lon_bot, lon_top, num_inds, start_ind
-integer :: x_ind, y_ind
-real(r8) :: p(4), x_corners(4), y_corners(4), xbot, xtop
-real(r8) :: lon_fract, lat_fract
-logical :: masked
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-! Succesful return has istatus of 0
-istatus = 0
-
-! Get the lower left corner for either grid type
-if(dipole_grid) then
- ! Figure out which of the regular grid boxes this is in
- call get_reg_box_indices(lon, lat, x_ind, y_ind)
-
- ! Is this on the U or T grid?
- if(is_on_ugrid(var_type)) then
- ! On U grid
- num_inds = u_dipole_num (x_ind, y_ind)
- start_ind = u_dipole_start(x_ind, y_ind)
-
- ! If there are no quads overlapping, can't do interpolation
- if(num_inds == 0) then
- istatus = 1
- return
- endif
-
- ! Search the list of quads to see if (lon, lat) is in one
- call get_dipole_quad(lon, lat, ulon, ulat, num_inds, start_ind, &
- u_dipole_lon_list, u_dipole_lat_list, lon_bot, lat_bot, istatus)
- ! Fail on bad istatus return
- if(istatus /= 0) return
-
- ! Getting corners for accurate interpolation
- call get_quad_corners(ulon, lon_bot, lat_bot, x_corners)
- call get_quad_corners(ulat, lon_bot, lat_bot, y_corners)
-
- ! Fail if point is in one of the U boxes that go through the
- ! pole (this could be fixed up if necessary)
- if(lat_bot == u_pole_y .and. (lon_bot == pole_x -1 .or. &
- lon_bot == pole_x)) then
- istatus = 4
- return
- endif
-
- else
- ! On T grid
- num_inds = t_dipole_num (x_ind, y_ind)
- start_ind = t_dipole_start(x_ind, y_ind)
- call get_dipole_quad(lon, lat, tlon, tlat, num_inds, start_ind, &
- t_dipole_lon_list, t_dipole_lat_list, lon_bot, lat_bot, istatus)
- ! Fail on bad istatus return
- if(istatus /= 0) return
-
- ! Fail if point is in T box that covers pole
- if(lon_bot == pole_x .and. lat_bot == t_pole_y) then
- istatus = 5
- return
- endif
-
- ! Getting corners for accurate interpolation
- call get_quad_corners(tlon, lon_bot, lat_bot, x_corners)
- call get_quad_corners(tlat, lon_bot, lat_bot, y_corners)
- endif
-
-else
- ! This is an irregular grid
- ! U and V are on velocity grid
- if (is_on_ugrid(var_type)) then
- ! Get the corner indices and the fraction of the distance between
- call get_irreg_box(lon, lat, ulon, ulat, &
- lon_bot, lat_bot, lon_fract, lat_fract, istatus)
- else
- ! Eta, T and S are on the T grid
- ! Get the corner indices
- call get_irreg_box(lon, lat, tlon, tlat, &
- lon_bot, lat_bot, lon_fract, lat_fract, istatus)
- endif
-
- ! Return passing through error status
- if(istatus /= 0) return
-
-endif
-
-! Find the indices to get the values for interpolating
-lat_top = lat_bot + 1
-if(lat_top > ny) then
- istatus = 2
- return
-endif
-
-! Watch for wraparound in longitude
-lon_top = lon_bot + 1
-if(lon_top > nx) lon_top = 1
-
-! Get the values at the four corners of the box or quad
-! Corners go around counterclockwise from lower left
-p(1) = get_val(lon_bot, lat_bot, nx, x, var_type, height, masked)
-if(masked) then
- istatus = 3
- return
-endif
-
-p(2) = get_val(lon_top, lat_bot, nx, x, var_type, height, masked)
-if(masked) then
- istatus = 3
- return
-endif
-
-p(3) = get_val(lon_top, lat_top, nx, x, var_type, height, masked)
-if(masked) then
- istatus = 3
- return
-endif
-
-p(4) = get_val(lon_bot, lat_top, nx, x, var_type, height, masked)
-if(masked) then
- istatus = 3
- return
-endif
-
-! Full bilinear interpolation for quads
-if(dipole_grid) then
- call quad_bilinear_interp(lon, lat, x_corners, y_corners, p, interp_val)
-else
- ! Rectangular biliear interpolation
- xbot = p(1) + lon_fract * (p(2) - p(1))
- xtop = p(4) + lon_fract * (p(3) - p(4))
- ! Now interpolate in latitude
- interp_val = xbot + lat_fract * (xtop - xbot)
-endif
-
-end subroutine lon_lat_interpolate
-
-!------------------------------------------------------------
-
-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
-
-! Returns the value from a single level array given the lat and lon indices
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list