[Dart-dev] [6043] DART/branches/development/models: add enough files to allow the quickbuild to finish
nancy at ucar.edu
nancy at ucar.edu
Thu Apr 4 09:10:08 MDT 2013
Revision: 6043
Author: nancy
Date: 2013-04-04 09:10:06 -0600 (Thu, 04 Apr 2013)
Log Message:
-----------
add enough files to allow the quickbuild to finish
for the CESM dir. totally untested.
Modified Paths:
--------------
DART/branches/development/models/CESM/dart_to_cesm.f90
DART/branches/development/models/CESM/model_mod.f90
DART/branches/development/models/CESM/work/path_names_filter
DART/branches/development/models/CESM/work/path_names_wakeup_filter
Added Paths:
-----------
DART/branches/development/models/POP/pop_model_mod.f90
DART/branches/development/models/cam/cam_model_mod.f90
DART/branches/development/models/clm/clm_model_mod.f90
-------------- next part --------------
Modified: DART/branches/development/models/CESM/dart_to_cesm.f90
===================================================================
--- DART/branches/development/models/CESM/dart_to_cesm.f90 2013-04-03 21:50:48 UTC (rev 6042)
+++ DART/branches/development/models/CESM/dart_to_cesm.f90 2013-04-04 15:10:06 UTC (rev 6043)
@@ -34,7 +34,7 @@
use time_manager_mod, only : time_type, print_time, print_date, operator(-)
use model_mod, only : static_init_model, sv_to_restart_file, &
get_model_size, get_cesm_restart_filename
-use dart_cesm_mod, only : write_cesm_namelist
+!use dart_cesm_mod, only : write_cesm_namelist
implicit none
@@ -110,7 +110,7 @@
call sv_to_restart_file(statevector, cesm_restart_filename, model_time)
if ( advance_time_present ) then
- call write_cesm_namelist(model_time, adv_to_time)
+ ! call write_cesm_namelist(model_time, adv_to_time)
endif
!----------------------------------------------------------------------
Modified: DART/branches/development/models/CESM/model_mod.f90
===================================================================
--- DART/branches/development/models/CESM/model_mod.f90 2013-04-03 21:50:48 UTC (rev 6042)
+++ DART/branches/development/models/CESM/model_mod.f90 2013-04-04 15:10:06 UTC (rev 6043)
@@ -60,7 +60,10 @@
get_close_maxdist_init, &
get_close_obs_init, &
get_close_obs, &
- ens_mean_for_model
+ ens_mean_for_model, &
+ restart_file_to_sv, &
+ sv_to_restart_file, &
+ get_cesm_restart_filename
! version controlled file description for error handling, do not edit
character(len=128), parameter :: &
@@ -73,18 +76,18 @@
character(len=256) :: msgstring
logical, save :: module_initialized = .false.
-logical :: output_state_vector
-integer :: assimilation_period_days, assimilation_period_seconds
-real(r8) :: model_perturbation_amplitude
+! FIXME: for now make cam the only default
+logical :: include_CAM = .true.
+logical :: include_POP = .false.
+logical :: include_CLM = .false.
integer :: debug = 0 ! turn up for more and more debug messages
namelist /model_nml/ &
- output_state_vector, &
- assimilation_period_days, & ! for now, this is the timestep
- assimilation_period_seconds, &
- model_perturbation_amplitude,&
+ include_CAM, &
+ include_POP, &
+ include_CLM, &
debug
type(time_type) :: model_time, model_timestep
@@ -120,9 +123,9 @@
if (do_output()) write( * , nml=model_nml)
-call cam_static_init_model()
-call pop_static_init_model()
-call clm_static_init_model()
+if (include_CAM) call cam_static_init_model()
+if (include_POP) call pop_static_init_model()
+if (include_CLM) call clm_static_init_model()
model_timestep = get_model_time_step()
@@ -187,10 +190,14 @@
if ( .not. module_initialized ) call static_init_model
-pop_model_size = pop_get_model_size()
-clm_model_size = clm_get_model_size()
-cam_model_size = cam_get_model_size()
+cam_model_size = 0
+pop_model_size = 0
+clm_model_size = 0
+if (include_CAM) cam_model_size = cam_get_model_size()
+if (include_POP) pop_model_size = pop_get_model_size()
+if (include_CLM) clm_model_size = clm_get_model_size()
+
get_model_size = cam_model_size + clm_model_size + pop_model_size
end function get_model_size
@@ -289,9 +296,9 @@
if ( .not. module_initialized ) call static_init_model
-cam_time = cam_get_model_time_step()
-pop_time = pop_get_model_time_step()
-clm_time = clm_get_model_time_step()
+if (include_CAM) cam_time = cam_get_model_time_step()
+if (include_POP) pop_time = pop_get_model_time_step()
+if (include_CLM) clm_time = clm_get_model_time_step()
! FIXME:
! make sure they are compatible here
@@ -347,9 +354,9 @@
! Shutdown and clean-up.
-call cam_end_model()
-call pop_end_model()
-call clm_end_model()
+if (include_CAM) call cam_end_model()
+if (include_POP) call pop_end_model()
+if (include_CLM) call clm_end_model()
end subroutine end_model
@@ -363,22 +370,28 @@
if ( .not. module_initialized ) call static_init_model
-rc = cam_nc_write_model_atts(ncFileID)
-if (rc /= 0) then
- nc_write_model_atts = rc
- return
+if (include_CAM) then
+ rc = cam_nc_write_model_atts(ncFileID)
+ if (rc /= 0) then
+ nc_write_model_atts = rc
+ return
+ endif
endif
-rc = pop_nc_write_model_atts(ncFileID)
-if (rc /= 0) then
- nc_write_model_atts = rc
- return
+if (include_POP) then
+ rc = pop_nc_write_model_atts(ncFileID)
+ if (rc /= 0) then
+ nc_write_model_atts = rc
+ return
+ endif
endif
-rc = clm_nc_write_model_atts(ncFileID)
-if (rc /= 0) then
- nc_write_model_atts = rc
- return
+if (include_CLM) then
+ rc = clm_nc_write_model_atts(ncFileID)
+ if (rc /= 0) then
+ nc_write_model_atts = rc
+ return
+ endif
endif
nc_write_model_atts = 0 ! If we got here, things went well.
@@ -394,26 +407,35 @@
integer, intent(in) :: timeindex
integer :: nc_write_model_vars ! function return value
-integer :: rc
+integer :: rc, x_start, x_end
if ( .not. module_initialized ) call static_init_model
-rc = cam_nc_write_model_vars(ncFileID, statevec, copyindex, timeindex)
-if (rc /= 0) then
- nc_write_model_vars = rc
- return
+if (include_CAM) then
+ call set_start_end('CAM', x_start, x_end)
+ rc = cam_nc_write_model_vars(ncFileID, statevec(x_start:x_end), copyindex, timeindex)
+ if (rc /= 0) then
+ nc_write_model_vars = rc
+ return
+ endif
endif
-rc = pop_nc_write_model_vars(ncFileID, statevec+cam_model_size, copyindex, timeindex)
-if (rc /= 0) then
- nc_write_model_vars = rc
- return
+if (include_POP) then
+ call set_start_end('POP', x_start, x_end)
+ rc = pop_nc_write_model_vars(ncFileID, statevec(x_start:x_end), copyindex, timeindex)
+ if (rc /= 0) then
+ nc_write_model_vars = rc
+ return
+ endif
endif
-rc = clm_nc_write_model_vars(ncFileID, statevec+cam_model_size+pop_model_size, copyindex, timeindex)
-if (rc /= 0) then
- nc_write_model_vars = rc
- return
+if (include_CLM) then
+ call set_start_end('CLM', x_start, x_end)
+ rc = clm_nc_write_model_vars(ncFileID, statevec(x_start:x_end), copyindex, timeindex)
+ if (rc /= 0) then
+ nc_write_model_vars = rc
+ return
+ endif
endif
nc_write_model_vars = 0 ! If we got here, things went well.
@@ -441,14 +463,20 @@
if ( .not. module_initialized ) call static_init_model
-call set_start_end('CAM', x_start, x_end)
-call cam_pert_model_state(state(x_start:x_end), pert_state(x_start:x_end), interf_provided)
+if (include_CAM) then
+ call set_start_end('CAM', x_start, x_end)
+ call cam_pert_model_state(state(x_start:x_end), pert_state(x_start:x_end), interf_provided)
+endif
-call set_start_end('POP', x_start, x_end)
-call pop_pert_model_state(state(x_start:x_end), pert_state(x_start:x_end), interf_provided)
+if (include_POP) then
+ call set_start_end('POP', x_start, x_end)
+ call pop_pert_model_state(state(x_start:x_end), pert_state(x_start:x_end), interf_provided)
+endif
-call set_start_end('CLM', x_start, x_end)
-call clm_pert_model_state(state(x_start:x_end), pert_state(x_start:x_end), interf_provided)
+if (include_CLM) then
+ call set_start_end('CLM', x_start, x_end)
+ call clm_pert_model_state(state(x_start:x_end), pert_state(x_start:x_end), interf_provided)
+endif
end subroutine pert_model_state
@@ -466,14 +494,20 @@
if ( .not. module_initialized ) call static_init_model
-call set_start_end('CAM', x_start, x_end)
-call cam_ens_mean_for_model(ens_mean(x_start:x_end))
+if (include_CAM) then
+ call set_start_end('CAM', x_start, x_end)
+ call cam_ens_mean_for_model(ens_mean(x_start:x_end))
+endif
-call set_start_end('POP', x_start, x_end)
-call pop_ens_mean_for_model(ens_mean(x_start:x_end))
+if (include_POP) then
+ call set_start_end('POP', x_start, x_end)
+ call pop_ens_mean_for_model(ens_mean(x_start:x_end))
+endif
-call set_start_end('CLM', x_start, x_end)
-call clm_ens_mean_for_model(ens_mean(x_start:x_end))
+if (include_CLM) then
+ call set_start_end('CLM', x_start, x_end)
+ call clm_ens_mean_for_model(ens_mean(x_start:x_end))
+endif
end subroutine ens_mean_for_model
@@ -493,19 +527,35 @@
state_vector = MISSING_R8
-call set_start_end('CAM', x_start, x_end)
-call cam_restart_file_to_sv(filename, state_vector(x_start:x_end), model_time)
+if (include_CAM) then
+ call set_start_end('CAM', x_start, x_end)
+ call cam_restart_file_to_sv(filename, state_vector(x_start:x_end), model_time)
+endif
-call set_start_end('POP', x_start, x_end)
-call pop_restart_file_to_sv(filename, state_vector(x_start:x_end), model_time)
+if (include_POP) then
+ call set_start_end('POP', x_start, x_end)
+ call pop_restart_file_to_sv(filename, state_vector(x_start:x_end), model_time)
+endif
-call set_start_end('CLM', x_start, x_end)
-call clm_restart_file_to_sv(filename, state_vector(x_start:x_end), model_time)
+if (include_CLM) then
+ call set_start_end('CLM', x_start, x_end)
+ call clm_restart_file_to_sv(filename, state_vector(x_start:x_end), model_time)
+endif
end subroutine restart_file_to_sv
!------------------------------------------------------------------
+subroutine get_cesm_restart_filename(filename)
+ character(len=*), intent(out) :: filename
+
+! FIXME:
+filename = 'dummy'
+
+end subroutine get_cesm_restart_filename
+
+!------------------------------------------------------------------
+
subroutine sv_to_restart_file(state_vector, filename, statedate)
real(r8), intent(in) :: state_vector(:)
character(len=*), intent(in) :: filename
Modified: DART/branches/development/models/CESM/work/path_names_filter
===================================================================
--- DART/branches/development/models/CESM/work/path_names_filter 2013-04-03 21:50:48 UTC (rev 6042)
+++ DART/branches/development/models/CESM/work/path_names_filter 2013-04-04 15:10:06 UTC (rev 6043)
@@ -11,7 +11,7 @@
models/POP/pop_model_mod.f90
models/POP/dart_pop_mod.f90
models/clm/clm_model_mod.f90
-mpi_utilities/null_mpi_utilities_mod.f90
+mpi_utilities/mpi_utilities_mod.f90
obs_def/obs_def_mod.f90
obs_kind/obs_kind_mod.f90
obs_model/obs_model_mod.f90
Modified: DART/branches/development/models/CESM/work/path_names_wakeup_filter
===================================================================
--- DART/branches/development/models/CESM/work/path_names_wakeup_filter 2013-04-03 21:50:48 UTC (rev 6042)
+++ DART/branches/development/models/CESM/work/path_names_wakeup_filter 2013-04-04 15:10:06 UTC (rev 6043)
@@ -1,5 +1,5 @@
common/types_mod.f90
filter/wakeup_filter.f90
-mpi_utilities/null_mpi_utilities_mod.f90
+mpi_utilities/mpi_utilities_mod.f90
time_manager/time_manager_mod.f90
utilities/utilities_mod.f90
Added: DART/branches/development/models/POP/pop_model_mod.f90
===================================================================
--- DART/branches/development/models/POP/pop_model_mod.f90 (rev 0)
+++ DART/branches/development/models/POP/pop_model_mod.f90 2013-04-04 15:10:06 UTC (rev 6043)
@@ -0,0 +1,3630 @@
+! DART software - Copyright 2004 - 2011 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
+
+module pop_model_mod
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+! 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=128), parameter :: &
+ source = '$URL$', &
+ revision = '$Revision$', &
+ 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(:)
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list