<p><b>duda</b> 2011-02-11 14:26:40 -0700 (Fri, 11 Feb 2011)</p><p>BRANCH COMMIT<br>
<br>
Updates to nonhydrostatic real-data initialization:<br>
<br>
- Add new modules: bitarray, queue, and hinterp<br>
<br>
- In module_test_cases.F: ensure that vegetation category, soil type, <br>
and land mask are all consistent; fix missing soil category over <br>
Antarctica; add fields tmn, seaice, dzs, dz, and surface_pressure;<br>
use new hinterp module to perform masked interpolation of land<br>
surface fields<br>
<br>
- In module_llxy.F, switch to RKIND for real kind, rather than using<br>
a hard-wired kind of 4<br>
<br>
<br>
M src/core_init_nhyd_atmos/module_llxy.F<br>
A src/core_init_nhyd_atmos/module_hinterp.F<br>
M src/core_init_nhyd_atmos/module_test_cases.F<br>
A src/core_init_nhyd_atmos/module_bitarray.F<br>
M src/core_init_nhyd_atmos/Registry<br>
M src/core_init_nhyd_atmos/Makefile<br>
A src/core_init_nhyd_atmos/module_queue.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_init_nhyd_atmos/Makefile
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/Makefile 2011-02-10 22:05:23 UTC (rev 737)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/Makefile 2011-02-11 21:26:40 UTC (rev 738)
@@ -5,6 +5,9 @@
module_advection.o \
module_read_met.o \
module_llxy.o \
+ module_bitarray.o \
+ module_queue.o \
+ module_hinterp.o \
read_geogrid.o
all: core_hyd
@@ -12,8 +15,10 @@
core_hyd: $(OBJS)
ar -ru libdycore.a $(OBJS)
-module_test_cases.o: module_advection.o module_read_met.o read_geogrid.o module_llxy.o
+module_test_cases.o: module_advection.o module_read_met.o read_geogrid.o module_llxy.o module_hinterp.o
+module_hinterp.o: module_queue.o module_bitarray.o
+
module_advection.o:
module_read_met.o:
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-02-10 22:05:23 UTC (rev 737)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-02-11 21:26:40 UTC (rev 738)
@@ -131,14 +131,13 @@
# static terrestrial fields
var persistent real ter ( nCells ) 0 io ter mesh - -
var persistent integer landmask ( nCells ) 0 io landmask mesh - -
-var persistent integer lu_index ( nCells ) 0 io lu_index mesh - -
-var persistent integer soilcat_top ( nCells ) 0 io soilcat_top mesh - -
+var persistent integer ivgtyp ( nCells ) 0 io lu_index mesh - -
+var persistent integer isltyp ( nCells ) 0 io soilcat_top mesh - -
var persistent integer soilcat_bot ( nCells ) 0 io soilcat_bot mesh - -
var persistent real snoalb ( nCells ) 0 io snoalb mesh - -
var persistent real greenfrac ( nMonths nCells ) 0 io greenfrac mesh - -
var persistent real shdmin ( nCells ) 0 io shdmin mesh - -
var persistent real shdmax ( nCells ) 0 io shdmax mesh - -
-var persistent real vegfra ( nCells ) 0 io vegfra mesh - -
var persistent real albedo12m ( nMonths nCells ) 0 io albedo12m mesh - -
# description of the vertical grid structure
@@ -162,37 +161,41 @@
var persistent real dss ( nVertLevels nCells ) 0 io dss mesh - -
# Horizontally interpolated from first-guess data
-var persistent real u_fg ( nFGLevels nEdges ) 0 - u fg - -
-var persistent real v_fg ( nFGLevels nEdges ) 0 - v fg - -
-var persistent real t_fg ( nFGLevels nCells ) 0 - t fg - -
-var persistent real p_fg ( nFGLevels nCells ) 0 - p fg - -
-var persistent real z_fg ( nFGLevels nCells ) 0 - z fg - -
-var persistent real rh_fg ( nFGLevels nCells ) 0 - rh fg - -
-var persistent real soilz_fg ( nCells ) 0 - soilz fg - -
-var persistent real psfc_fg ( nCells ) 0 - psfc fg - -
-var persistent real pmsl_fg ( nCells ) 0 - pmsl fg - -
+var persistent real u_fg ( nFGLevels nEdges Time ) 1 - u fg - -
+var persistent real v_fg ( nFGLevels nEdges Time ) 1 - v fg - -
+var persistent real t_fg ( nFGLevels nCells Time ) 1 - t fg - -
+var persistent real p_fg ( nFGLevels nCells Time ) 1 - p fg - -
+var persistent real z_fg ( nFGLevels nCells Time ) 1 - z fg - -
+var persistent real rh_fg ( nFGLevels nCells Time ) 1 - rh fg - -
+var persistent real soilz_fg ( nCells Time ) 1 - soilz fg - -
+var persistent real psfc_fg ( nCells Time ) 1 - psfc fg - -
+var persistent real pmsl_fg ( nCells Time ) 1 - pmsl fg - -
# Horizontally interpolated from first-guess data
# and should be read in by model
-var persistent real dzs ( nSoilLevels nCells ) 0 o dzs fg - -
-var persistent real sh2o ( nSoilLevels nCells ) 0 o sh2o fg - -
-var persistent real smois ( nSoilLevels nCells ) 0 o smois fg - -
-var persistent real tslb ( nSoilLevels nCells ) 0 o tslb fg - -
-var persistent real skintemp ( nCells ) 0 o skintemp fg - -
-var persistent real sst ( nCells ) 0 o sst fg - -
-var persistent real snow ( nCells ) 0 o snow fg - -
-var persistent real snowc ( nCells ) 0 o snowc fg - -
-var persistent real xice ( nCells ) 0 o xice fg - -
-var persistent real gfs_z ( nVertLevels nCells ) 0 - gfs_z fg - -
+var persistent real dz ( nSoilLevels nCells Time ) 1 o dz fg - -
+var persistent real dzs ( nSoilLevels nCells Time ) 1 o dzs fg - -
+var persistent real sh2o ( nSoilLevels nCells Time ) 1 o sh2o fg - -
+var persistent real smois ( nSoilLevels nCells Time ) 1 o smois fg - -
+var persistent real tslb ( nSoilLevels nCells Time ) 1 o tslb fg - -
+var persistent real tmn ( nCells Time ) 1 o tmn fg - -
+var persistent real skintemp ( nCells Time ) 1 o skintemp fg - -
+var persistent real sst ( nCells Time ) 1 o sst fg - -
+var persistent real snow ( nCells Time ) 1 o snow fg - -
+var persistent real snowc ( nCells Time ) 1 o snowc fg - -
+var persistent real xice ( nCells Time ) 1 o xice fg - -
+var persistent real seaice ( nCells Time ) 1 o seaice fg - -
+var persistent real gfs_z ( nVertLevels nCells Time ) 1 - gfs_z fg - -
+var persistent real vegfra ( nCells Time ) 1 io vegfra mesh - -
# Prognostic variables: read from input, saved in restart, and written to output
-var persistent real u ( nVertLevels nEdges Time ) 2 io u state - -
-var persistent real w ( nVertLevelsP1 nCells Time ) 2 io w state - -
-var persistent real rho ( nVertLevels nCells Time ) 2 io rho state - -
-var persistent real theta ( nVertLevels nCells Time ) 2 io theta state - -
-var persistent real qv ( nVertLevels nCells Time ) 2 io qv state scalars moist
-var persistent real qc ( nVertLevels nCells Time ) 2 io qc state scalars moist
-var persistent real qr ( nVertLevels nCells Time ) 2 io qr state scalars moist
+var persistent real u ( nVertLevels nEdges Time ) 2 o u state - -
+var persistent real w ( nVertLevelsP1 nCells Time ) 2 o w state - -
+var persistent real rho ( nVertLevels nCells Time ) 2 o rho state - -
+var persistent real theta ( nVertLevels nCells Time ) 2 o theta state - -
+var persistent real qv ( nVertLevels nCells Time ) 2 o qv state scalars moist
+var persistent real qc ( nVertLevels nCells Time ) 2 o qc state scalars moist
+var persistent real qr ( nVertLevels nCells Time ) 2 o qr state scalars moist
# state variables diagnosed from prognostic state
var persistent real pressure_p ( nVertLevels nCells Time ) 1 - pressure_p diag - -
@@ -212,13 +215,15 @@
var persistent real exner ( nVertLevels nCells Time ) 1 - exner diag - -
var persistent real exner_base ( nVertLevels nCells Time ) 1 io exner_base diag - -
var persistent real rtheta_base ( nVertLevels nCells Time ) 1 - rtheta_base diag - -
-var persistent real pressure ( nVertLevels nCells Time ) 1 io pressure diag - -
+var persistent real pressure ( nVertLevels nCells Time ) 1 - pressure diag - -
var persistent real pressure_base ( nVertLevels nCells Time ) 1 io pressure_base diag - -
var persistent real rho_base ( nVertLevels nCells Time ) 1 io rho_base diag - -
var persistent real theta_base ( nVertLevels nCells Time ) 1 io theta_base diag - -
var persistent real cqw ( nVertLevels nCells Time ) 1 - cqw diag - -
+var persistent real surface_pressure ( nCells Time ) 1 o surface_pressure diag - -
+
# coupled variables needed by the solver, but not output...
var persistent real ru ( nVertLevels nEdges Time ) 1 - ru diag - -
var persistent real rw ( nVertLevelsP1 nCells Time ) 1 - rw diag - -
Added: branches/atmos_physics/src/core_init_nhyd_atmos/module_bitarray.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_bitarray.F (rev 0)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_bitarray.F 2011-02-11 21:26:40 UTC (rev 738)
@@ -0,0 +1,197 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Module: bitarray_module
+!
+! Purpose: This module provides a two-dimensional bit array and a set of
+! routines to manipulate and examine the bits of the array.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+module bitarray
+
+ type bitarray
+ integer, pointer, dimension(:,:) :: iarray ! Storage array
+ integer :: nx, ny ! Number of bits in the x and y directions
+ integer :: x_int_dim, y_int_dim ! Number of integers in the x and y directions
+ integer :: integer_size ! Number of bits in an integer
+ end type bitarray
+
+ contains
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: bitarray_create
+ !
+ ! Purpose: Allocate and initialize a bit array so that all bits are FALSE
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine bitarray_create(b, i, j)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: i, j
+ type (bitarray), intent(out) :: b
+
+ b%integer_size = bit_size(b%integer_size)
+
+ b%nx = i
+ b%ny = j
+
+ b%x_int_dim = ceiling(real(b%nx)/real(b%integer_size))
+ b%y_int_dim = b%ny
+
+ nullify(b%iarray)
+ allocate(b%iarray(b%x_int_dim, b%y_int_dim))
+ b%iarray = 0
+
+ end subroutine bitarray_create
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: bitarray_copy
+ !
+ ! Purpose: Duplicate a bitarray.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine bitarray_copy(src, dst)
+
+ implicit none
+
+ ! Arguments
+ type (bitarray), intent(in) :: src
+ type (bitarray), intent(out) :: dst
+
+ dst%integer_size = src%integer_size
+
+ dst%nx = src%nx
+ dst%ny = src%ny
+
+ dst%x_int_dim = src%x_int_dim
+ dst%y_int_dim = src%y_int_dim
+
+ allocate(dst%iarray(dst%x_int_dim, dst%y_int_dim))
+ dst%iarray = src%iarray
+
+ end subroutine bitarray_copy
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: bitarray_set
+ !
+ ! Purpose: Set the bit located at (i,j) to TRUE
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine bitarray_set(b, i, j)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: i, j
+ type (bitarray), intent(inout) :: b
+
+ ! Local variables
+ integer :: n_integer, n_bit
+
+ n_integer = ((i-1) / b%integer_size) + 1
+ n_bit = mod((i-1), b%integer_size)
+
+ b%iarray(n_integer, j) = ibset(b%iarray(n_integer, j), n_bit)
+
+ end subroutine bitarray_set
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: bitarray_clear
+ !
+ ! Purpose: Set the bit located at (i,j) to FALSE
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine bitarray_clear(b, i, j)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: i, j
+ type (bitarray), intent(inout) :: b
+
+ ! Local variables
+ integer :: n_integer, n_bit
+
+ n_integer = ((i-1) / b%integer_size) + 1
+ n_bit = mod((i-1), b%integer_size)
+
+ b%iarray(n_integer, j) = ibclr(b%iarray(n_integer, j), n_bit)
+
+ end subroutine bitarray_clear
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: bitarray_test
+ !
+ ! Purpose: To return the value of the bit located at (i,j)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ function bitarray_test(b, i, j)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: i, j
+ type (bitarray), intent(in) :: b
+
+ ! Local variables
+ logical :: bitarray_test
+ integer :: n_integer, n_bit
+
+ n_integer = ((i-1) / b%integer_size) + 1
+ n_bit = mod((i-1), b%integer_size)
+
+ bitarray_test = btest(b%iarray(n_integer,j), n_bit)
+
+ end function bitarray_test
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: bitarray_merge
+ !
+ ! Purpose: The first bitarray argument, b1, is set to the union of the .TRUE.
+ ! bits in b1 and b2. That is, after returning, a bit x in b1 is set if
+ ! either x was set in b1 or x was set in b2. Thus, b1 AND b2 MUST BE BIT
+ ! ARRAYS OF THE SAME DIMENSIONS.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine bitarray_merge(b1, b2)
+
+ implicit none
+
+ ! Arguments
+ type (bitarray), intent(inout) :: b1, b2
+
+ ! Local variables
+ integer :: i, j
+
+ if (b1%x_int_dim /= b2%x_int_dim .or. b1%y_int_dim /= b2%y_int_dim) then
+ write(0,*) 'ERROR: In bitarray_merge(), b1 and b2 have different dimensions.'
+ end if
+
+ do i=1,b1%x_int_dim
+ do j=1,b1%y_int_dim
+ b1%iarray(i,j) = ior(b1%iarray(i,j), b2%iarray(i,j))
+ end do
+ end do
+
+ end subroutine bitarray_merge
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: bitarray_destroy
+ !
+ ! Purpose: To deallocate all allocated memory associated with the bit array
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine bitarray_destroy(b)
+
+ implicit none
+
+ ! Arguments
+ type (bitarray), intent(inout) :: b
+
+ if (associated(b%iarray)) then
+ deallocate(b%iarray)
+ else
+ write(0,*) 'WARNING: In bitarray_destroy(), b is not allocated.'
+ end if
+
+ end subroutine bitarray_destroy
+
+end module bitarray
Added: branches/atmos_physics/src/core_init_nhyd_atmos/module_hinterp.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_hinterp.F (rev 0)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_hinterp.F 2011-02-11 21:26:40 UTC (rev 738)
@@ -0,0 +1,1050 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! MODULE INTERP_MODULE
+!
+! This module provides routines for interpolation.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+module hinterp
+
+ use bitarray
+ use queue
+
+ integer, parameter :: SIXTEEN_POINT=1, FOUR_POINT=2, N_NEIGHBOR=3, &
+ AVERAGE4=4, AVERAGE16=5, W_AVERAGE4=6, W_AVERAGE16=7, &
+ SEARCH=8
+
+ contains
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: interp_array_from_string
+ !
+ ! Purpose:
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ function interp_array_from_string(interp_string)
+
+ implicit none
+
+ ! Arguments
+ character (len=*), intent(in) :: interp_string
+
+ ! Local variables
+ integer :: j, p1, p2, iend, num_methods
+
+ ! Return value
+ integer, pointer, dimension(:) :: interp_array_from_string
+
+ ! Get an idea of how many interpolation methods are in the string
+ ! so we can allocate an appropriately sized array
+ num_methods = 1
+ do j=1,len_trim(interp_string)
+ if (interp_string(j:j) == '+') num_methods = num_methods + 1
+ end do
+
+ allocate(interp_array_from_string(num_methods+1))
+ interp_array_from_string = 0
+
+ iend = len_trim(interp_string)
+
+ p1 = 1
+ p2 = index(interp_string(1:iend),'+')
+ j = 1
+ do while(p2 >= p1)
+ if (index(interp_string(p1:p2-1),'nearest_neighbor') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('nearest_neighbor')) then
+ interp_array_from_string(j) = N_NEIGHBOR
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'average_4pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('average_4pt')) then
+ interp_array_from_string(j) = AVERAGE4
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'average_16pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('average_16pt')) then
+ interp_array_from_string(j) = AVERAGE16
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'wt_average_4pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_4pt')) then
+ interp_array_from_string(j) = W_AVERAGE4
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'wt_average_16pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_16pt')) then
+ interp_array_from_string(j) = W_AVERAGE16
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'four_pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('four_pt')) then
+ interp_array_from_string(j) = FOUR_POINT
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'sixteen_pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('sixteen_pt')) then
+ interp_array_from_string(j) = SIXTEEN_POINT
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'search') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('search')) then
+ interp_array_from_string(j) = SEARCH
+ j = j + 1
+ else
+ if (index(interp_string(p1:p2-1),'average_gcell') == 0) &
+ write(0,*) 'WARNING: Unrecognized interpolation method '//interp_string(p1:p2-1)
+ end if
+ p1 = p2 + 1
+ p2 = index(interp_string(p1:iend),'+') + p1 - 1
+ end do
+
+ p2 = iend+1
+ if (p1 < iend) then
+ if (index(interp_string(p1:p2-1),'nearest_neighbor') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('nearest_neighbor')) then
+ interp_array_from_string(j) = N_NEIGHBOR
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'average_4pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('average_4pt')) then
+ interp_array_from_string(j) = AVERAGE4
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'average_16pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('average_16pt')) then
+ interp_array_from_string(j) = AVERAGE16
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'wt_average_4pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_4pt')) then
+ interp_array_from_string(j) = W_AVERAGE4
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'wt_average_16pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_16pt')) then
+ interp_array_from_string(j) = W_AVERAGE16
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'four_pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('four_pt')) then
+ interp_array_from_string(j) = FOUR_POINT
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'sixteen_pt') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('sixteen_pt')) then
+ interp_array_from_string(j) = SIXTEEN_POINT
+ j = j + 1
+ else if (index(interp_string(p1:p2-1),'search') /= 0 .and. &
+ len_trim(interp_string(p1:p2-1)) == len_trim('search')) then
+ interp_array_from_string(j) = SEARCH
+ j = j + 1
+ else
+ if (index(interp_string(p1:p2-1),'average_gcell') == 0) &
+ write(0,*) 'WARNING: Unrecognized interpolation method '//interp_string(p1:p2-1)
+ end if
+ end if
+
+ return
+
+ end function interp_array_from_string
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: interp_sequence
+ !
+ ! Purpose: Delegates the actual task of interpolation to specific
+ ! interpolation routines defined in the module.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ recursive function interp_sequence(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: start_x, start_y, start_z
+ integer, intent(in) :: end_x, end_y, end_z
+ integer, intent(in) :: izz ! The z-index of the 2d-array to
+ ! interpolate within
+ real (kind=RKIND), intent(in) :: xx , yy ! The location to interpolate to
+ real (kind=RKIND), intent(in) :: msgval
+ real (kind=RKIND), intent(in), optional :: maskval
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array
+ integer, intent(in) :: idx
+ integer, dimension(:), intent(in) :: interp_list
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array
+
+ ! Return value
+ real (kind=RKIND) :: interp_sequence
+
+ ! No more interpolation methods to try
+ if (interp_list(idx) == 0) then
+!write(0,*) 'Returning missing value'
+ interp_sequence = msgval
+ return
+ end if
+
+ if (interp_list(idx) == FOUR_POINT) then
+!write(0,*) 'Calling four_pt'
+ interp_sequence = four_pt(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, &
+ msgval, interp_list, idx+1, maskval, mask_array)
+
+ else if (interp_list(idx) == AVERAGE4) then
+!write(0,*) 'Calling four_pt_average'
+ interp_sequence = four_pt_average(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, &
+ msgval, interp_list, idx+1, maskval, mask_array)
+
+ else if (interp_list(idx) == W_AVERAGE4) then
+!write(0,*) 'Calling wt_four_pt_average'
+ interp_sequence = wt_four_pt_average(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, &
+ msgval, interp_list, idx+1, maskval, mask_array)
+
+ else if (interp_list(idx) == N_NEIGHBOR) then
+!write(0,*) 'Calling nearest_neighbor'
+ interp_sequence = nearest_neighbor(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, &
+ msgval, interp_list, idx+1, maskval, mask_array)
+
+ else if (interp_list(idx) == SIXTEEN_POINT) then
+!write(0,*) 'Calling sixteen_pt'
+ interp_sequence = sixteen_pt(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, &
+ msgval, interp_list, idx+1, maskval, mask_array)
+
+ else if (interp_list(idx) == SEARCH) then
+!write(0,*) 'Calling search_extrap'
+ interp_sequence = search_extrap(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, &
+ msgval, interp_list, idx+1, maskval, mask_array)
+
+ else if (interp_list(idx) == AVERAGE16) then
+!write(0,*) 'Calling sixteen_pt_average'
+ interp_sequence = sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, &
+ msgval, interp_list, idx+1, maskval, mask_array)
+
+ else if (interp_list(idx) == W_AVERAGE16) then
+!write(0,*) 'Calling wt_sixteen_pt_average'
+ interp_sequence = wt_sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, &
+ msgval, interp_list, idx+1, maskval, mask_array)
+ end if
+
+ end function interp_sequence
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: nearest_neighbor
+ !
+ ! Purpose: Returns the point nearest to (xx,yy). If (xx,yy) is outside of the
+ ! array, the point on the edge of the array nearest to (xx,yy) is returned.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ recursive function nearest_neighbor(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: start_x, start_y, start_z
+ integer, intent(in) :: end_x, end_y, end_z
+ integer, intent(in) :: izz ! The z-index of the 2d-array to
+ ! interpolate within
+ real (kind=RKIND), intent(in) :: xx , yy ! The location to interpolate to
+ real (kind=RKIND), intent(in) :: msgval
+ real (kind=RKIND), intent(in), optional :: maskval
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array
+ integer, dimension(:), intent(in) :: interp_list
+ integer, intent(in) :: idx
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array
+
+ ! Return value
+ real (kind=RKIND) :: nearest_neighbor
+
+ ! Local variables
+ integer :: ix, iy
+
+ ix = nint(xx)
+ iy = nint(yy)
+
+ ! The first thing to do is to ensure that the point (xx,yy) is within the array
+ if (ix < start_x .or. ix > end_x) then
+ nearest_neighbor = msgval
+ return
+ end if
+
+ if (iy < start_y .or. iy > end_y) then
+ nearest_neighbor = msgval
+ return
+ end if
+
+ if (present(mask_array) .and. present(maskval)) then
+ if (mask_array(ix,iy) == maskval) then
+ nearest_neighbor = msgval
+ else
+ nearest_neighbor = array(ix,iy,izz)
+ end if
+ else
+ nearest_neighbor = array(ix,iy,izz)
+ end if
+
+ ! If we have a missing value, try the next interpolation method in the sequence
+ if (nearest_neighbor == msgval) then
+ nearest_neighbor = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ end if
+
+ end function nearest_neighbor
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: search_extrap
+ !
+ ! Purpose: Returns the point nearest to (xx,yy) that has a non-missing value.
+ ! If no valid value can be found in the array, msgval is returned.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ recursive function search_extrap(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: start_x, start_y, start_z
+ integer, intent(in) :: end_x, end_y, end_z
+ integer, intent(in) :: izz ! The z-index of the 2d-array to search
+ real (kind=RKIND), intent(in) :: xx , yy ! The location of the search origin
+ real (kind=RKIND), intent(in) :: msgval
+ real (kind=RKIND), intent(in), optional :: maskval
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
+ intent(in) :: array
+ integer, dimension(:), intent(in) :: interp_list
+ integer, intent(in) :: idx
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y), &
+ intent(in), optional :: mask_array
+
+ ! Return value
+ real (kind=RKIND) :: search_extrap
+
+ ! Local variables
+ integer :: i, j
+ real (kind=RKIND) :: distance
+ logical :: found_valid
+ type (bitarray) :: b
+ type (queue) :: q
+ type (q_data) :: qdata
+
+ ! We only search if the starting point is within the array
+ if (nint(xx) < start_x .or. nint(xx) > end_x .or. &
+ nint(yy) < start_y .or. nint(yy) > end_y) then
+ search_extrap = msgval
+ return
+ end if
+
+ call bitarray_create(b, (end_x-start_x+1), (end_y-start_y+1))
+ call q_init(q)
+
+ found_valid = .false.
+ qdata%x = nint(xx)
+ qdata%y = nint(yy)
+ call q_insert(q, qdata)
+ call bitarray_set(b, qdata%x-start_x+1, qdata%y-start_y+1)
+
+ do while (q_isdata(q) .and. (.not. found_valid))
+ qdata = q_remove(q)
+ i = qdata%x
+ j = qdata%y
+
+ if (present(mask_array) .and. present(maskval)) then
+ if (array(i,j,izz) /= msgval .and. mask_array(i,j) /= maskval) found_valid = .true.
+ else
+ if (array(i,j,izz) /= msgval) found_valid = .true.
+ end if
+
+ if (i-1 >= start_x) then
+ if (.not. bitarray_test(b, (i-1)-start_x+1, j-start_y+1)) then
+ qdata%x = i-1
+ qdata%y = j
+ call q_insert(q, qdata)
+ call bitarray_set(b, (i-1)-start_x+1, j-start_y+1)
+ end if
+ end if
+ if (i+1 <= end_x) then
+ if (.not. bitarray_test(b, (i+1)-start_x+1, j-start_y+1)) then
+ qdata%x = i+1
+ qdata%y = j
+ call q_insert(q, qdata)
+ call bitarray_set(b, (i+1)-start_x+1, j-start_y+1)
+ end if
+ end if
+ if (j-1 >= start_y) then
+ if (.not. bitarray_test(b, i-start_x+1, (j-1)-start_y+1)) then
+ qdata%x = i
+ qdata%y = j-1
+ call q_insert(q, qdata)
+ call bitarray_set(b, i-start_x+1, (j-1)-start_y+1)
+ end if
+ end if
+ if (j+1 <= end_y) then
+ if (.not. bitarray_test(b, i-start_x+1, (j+1)-start_y+1)) then
+ qdata%x = i
+ qdata%y = j+1
+ call q_insert(q, qdata)
+ call bitarray_set(b, i-start_x+1, (j+1)-start_y+1)
+ end if
+ end if
+ end do
+
+ if (found_valid) then
+ distance = (real(i)-xx)*(real(i)-xx)+(real(j)-yy)*(real(j)-yy)
+ search_extrap = array(i,j,izz)
+ do while (q_isdata(q))
+ qdata = q_remove(q)
+ if (present(mask_array) .and. present(maskval)) then
+ if (array(qdata%x,qdata%y,izz) /= msgval .and. mask_array(qdata%x,qdata%y) /= maskval) then
+ if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
+ distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
+ search_extrap = array(qdata%x, qdata%y, izz)
+ end if
+ end if
+
+ else
+ if (array(qdata%x,qdata%y,izz) /= msgval) then
+ if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
+ distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
+ search_extrap = array(qdata%x, qdata%y, izz)
+ end if
+ end if
+ end if
+ end do
+ else
+ search_extrap = msgval
+ end if
+
+ call q_destroy(q)
+ call bitarray_destroy(b)
+
+ end function search_extrap
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: four_pt_average
+ !
+ ! Purpose: Average of four surrounding grid point values
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ recursive function four_pt_average(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: start_x, start_y, start_z
+ integer, intent(in) :: end_x, end_y, end_z
+ integer, intent(in) :: izz ! The z-index of the 2d-array to
+ ! interpolate within
+ real (kind=RKIND), intent(in) :: xx, yy ! The location to interpolate to
+ real (kind=RKIND), intent(in) :: msgval
+ real (kind=RKIND), intent(in), optional :: maskval
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
+ intent(in) :: array
+ integer, dimension(:), intent(in) :: interp_list
+ integer, intent(in) :: idx
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y), &
+ intent(in), optional :: mask_array
+
+ ! Return value
+ real (kind=RKIND) :: four_pt_average
+
+ ! Local variables
+ integer :: ifx, ify, icx, icy
+ real (kind=RKIND) :: fxfy, fxcy, cxfy, cxcy
+
+ fxfy = 1.0
+ fxcy = 1.0
+ cxfy = 1.0
+ cxcy = 1.0
+
+ ifx = floor(xx)
+ icx = ceiling(xx)
+ ify = floor(yy)
+ icy = ceiling(yy)
+
+ ! First, make sure that the point is contained in the source array
+ if (ifx < start_x .or. icx > end_x .or. &
+ ify < start_y .or. icy > end_y) then
+
+ ! But if the point is at most half a grid point out, we can
+ ! still proceed with modified ifx, icx, ify, and icy.
+ if (xx > real(start_x)-0.5 .and. ifx < start_x) then
+ ifx = start_x
+ icx = start_x
+ else if (xx < real(end_x)+0.5 .and. icx > end_x) then
+ ifx = end_x
+ icx = end_x
+ end if
+
+ if (yy > real(start_y)-0.5 .and. ify < start_y) then
+ ify = start_y
+ icy = start_y
+ else if (yy < real(end_y)+0.5 .and. icy > end_y) then
+ ify = end_y
+ icy = end_y
+ end if
+
+ if (ifx < start_x .or. icx > end_x .or. &
+ ify < start_y .or. icy > end_y) then
+ four_pt_average = msgval
+ return
+ end if
+ end if
+
+ if (present(mask_array) .and. present(maskval)) then
+ if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0
+ if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0
+ if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0
+ if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0
+ else
+ if (array(ifx, ify, izz) == msgval) fxfy = 0.0
+ if (array(ifx, icy, izz) == msgval) fxcy = 0.0
+ if (array(icx, ify, izz) == msgval) cxfy = 0.0
+ if (array(icx, icy, izz) == msgval) cxcy = 0.0
+ end if
+
+ ! If all four points are missing, try the next interpolation method in the sequence
+ if (fxfy == 0.0 .and. fxcy == 0.0 .and. cxfy == 0.0 .and. cxcy == 0.0) then
+ four_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ else
+ four_pt_average = (fxfy * array(ifx, ify, izz) + &
+ fxcy * array(ifx, icy, izz) + &
+ cxfy * array(icx, ify, izz) + &
+ cxcy * array(icx, icy, izz) ) / (fxfy + fxcy + cxfy + cxcy)
+ end if
+
+ end function four_pt_average
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: wt_four_pt_average
+ !
+ ! Purpose: Weighted average of four surrounding grid point values
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ recursive function wt_four_pt_average(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: start_x, start_y, start_z
+ integer, intent(in) :: end_x, end_y, end_z
+ integer, intent(in) :: izz ! The z-index of the 2d-array to
+ ! interpolate within
+ real (kind=RKIND), intent(in) :: xx, yy ! The location to interpolate to
+ real (kind=RKIND), intent(in) :: msgval
+ real (kind=RKIND), intent(in), optional :: maskval
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
+ intent(in) :: array
+ integer, dimension(:), intent(in) :: interp_list
+ integer, intent(in) :: idx
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y), &
+ intent(in), optional :: mask_array
+
+ ! Return value
+ real (kind=RKIND) :: wt_four_pt_average
+
+ ! Local variables
+ integer :: ifx, ify, icx, icy
+ real (kind=RKIND) :: fxfy, fxcy, cxfy, cxcy
+
+ ifx = floor(xx)
+ icx = ceiling(xx)
+ ify = floor(yy)
+ icy = ceiling(yy)
+
+ fxfy = max(0., 1.0 - sqrt((xx-real(ifx))**2+(yy-real(ify))**2))
+ fxcy = max(0., 1.0 - sqrt((xx-real(ifx))**2+(yy-real(icy))**2))
+ cxfy = max(0., 1.0 - sqrt((xx-real(icx))**2+(yy-real(ify))**2))
+ cxcy = max(0., 1.0 - sqrt((xx-real(icx))**2+(yy-real(icy))**2))
+
+ ! First, make sure that the point is contained in the source array
+ if (ifx < start_x .or. icx > end_x .or. &
+ ify < start_y .or. icy > end_y) then
+
+ ! But if the point is at most half a grid point out, we can
+ ! still proceed with modified ifx, icx, ify, and icy.
+ if (xx > real(start_x)-0.5 .and. ifx < start_x) then
+ ifx = start_x
+ icx = start_x
+ else if (xx < real(end_x)+0.5 .and. icx > end_x) then
+ ifx = end_x
+ icx = end_x
+ end if
+
+ if (yy > real(start_y)-0.5 .and. ifx < start_y) then
+ ify = start_y
+ icy = start_y
+ else if (yy < real(end_y)+0.5 .and. icy > end_y) then
+ ify = end_y
+ icy = end_y
+ end if
+
+ if (ifx < start_x .or. icx > end_x .or. &
+ ify < start_y .or. icy > end_y) then
+ wt_four_pt_average = msgval
+ return
+ end if
+ end if
+
+ if (present(mask_array) .and. present(maskval)) then
+ if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0
+ if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0
+ if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0
+ if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0
+ else
+ if (array(ifx, ify, izz) == msgval) fxfy = 0.0
+ if (array(ifx, icy, izz) == msgval) fxcy = 0.0
+ if (array(icx, ify, izz) == msgval) cxfy = 0.0
+ if (array(icx, icy, izz) == msgval) cxcy = 0.0
+ end if
+
+ ! If all four points are missing, try the next interpolation method in the sequence
+ if (fxfy == 0.0 .and. fxcy == 0.0 .and. cxfy == 0.0 .and. cxcy == 0.0) then
+ wt_four_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ else
+ wt_four_pt_average = (fxfy * array(ifx, ify, izz) + &
+ fxcy * array(ifx, icy, izz) + &
+ cxfy * array(icx, ify, izz) + &
+ cxcy * array(icx, icy, izz) ) / (fxfy + fxcy + cxfy + cxcy)
+ end if
+
+ end function wt_four_pt_average
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: sixteen_pt_average
+ !
+ ! Purpose: Average of sixteen surrounding grid point values
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ recursive function sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: start_x, start_y, start_z
+ integer, intent(in) :: end_x, end_y, end_z
+ integer, intent(in) :: izz ! The z-index of the 2d-array to
+ ! interpolate within
+ real (kind=RKIND), intent(in) :: xx , yy ! The location to interpolate to
+ real (kind=RKIND), intent(in) :: msgval
+ real (kind=RKIND), intent(in), optional :: maskval
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
+ intent(in) :: array
+ integer, dimension(:), intent(in) :: interp_list
+ integer, intent(in) :: idx
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y), &
+ intent(in), optional :: mask_array
+
+ ! Return value
+ real (kind=RKIND) :: sixteen_pt_average
+
+ ! Local variables
+ integer :: i, j, ifx, ify
+ real (kind=RKIND) :: sum, sum_weight
+ real (kind=RKIND), dimension(4,4) :: weights
+
+ ifx = floor(xx)
+ ify = floor(yy)
+
+ ! First see whether the point is far enough within the array to
+ ! allow for a sixteen point average.
+ if (ifx < start_x+1 .or. ifx > end_x-2 .or. &
+ ify < start_y+1 .or. ify > end_y-2) then
+ sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ return
+ end if
+
+ sum_weight = 0.0
+ do i=1,4
+ do j=1,4
+
+ if (present(mask_array) .and. present(maskval)) then
+ if (array(ifx+3-i, ify+3-j, izz) == msgval .or. mask_array(ifx+3-i, ify+3-j) == maskval) then
+ weights(i,j) = 0.0
+ else
+ weights(i,j) = 1.0
+ end if
+ else
+ if (array(ifx+3-i, ify+3-j, izz) == msgval) then
+ weights(i,j) = 0.0
+ else
+ weights(i,j) = 1.0
+ end if
+ end if
+
+ sum_weight = sum_weight + weights(i,j)
+
+ end do
+ end do
+
+ ! If all points are missing, try the next interpolation method in the sequence
+ if (sum_weight == 0.0) then
+ sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ else
+ sum = 0.0
+ do i=1,4
+ do j=1,4
+ sum = sum + weights(i,j) * array(ifx+3-i, ify+3-j, izz)
+ end do
+ end do
+ sixteen_pt_average = sum / sum_weight
+ end if
+
+ end function sixteen_pt_average
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: wt_sixteen_pt_average
+ !
+ ! Purpose: Weighted average of sixteen surrounding grid point values
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ recursive function wt_sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
+ start_y, end_y, start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: start_x, start_y, start_z
+ integer, intent(in) :: end_x, end_y, end_z
+ integer, intent(in) :: izz ! The z-index of the 2d-array to
+ ! interpolate within
+ real (kind=RKIND), intent(in) :: xx , yy ! The location to interpolate to
+ real (kind=RKIND), intent(in) :: msgval
+ real (kind=RKIND), intent(in), optional :: maskval
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
+ intent(in) :: array
+ integer, dimension(:), intent(in) :: interp_list
+ integer, intent(in) :: idx
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y), &
+ intent(in), optional :: mask_array
+
+ ! Return value
+ real (kind=RKIND) :: wt_sixteen_pt_average
+
+ ! Local variables
+ integer :: i, j, ifx, ify
+ real (kind=RKIND) :: sum, sum_weight
+ real (kind=RKIND), dimension(4,4) :: weights
+
+ ifx = floor(xx)
+ ify = floor(yy)
+
+ ! First see whether the point is far enough within the array to
+ ! allow for a sixteen point average.
+ if (ifx < start_x+1 .or. ifx > end_x-2 .or. &
+ ify < start_y+1 .or. ify > end_y-2) then
+ wt_sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ return
+ end if
+
+ sum_weight = 0.0
+ do i=1,4
+ do j=1,4
+
+ if (present(mask_array) .and. present(maskval)) then
+ if (array(ifx+3-i, ify+3-j, izz) == msgval .or. mask_array(ifx+3-i, ify+3-j) == maskval) then
+ weights(i,j) = 0.0
+ else
+ weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
+ end if
+ else
+ if (array(ifx+3-i, ify+3-j, izz) == msgval) then
+ weights(i,j) = 0.0
+ else
+ weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
+ end if
+ end if
+
+ sum_weight = sum_weight + weights(i,j)
+
+ end do
+ end do
+
+ ! If all points are missing, try the next interpolation method in the sequence
+ if (sum_weight == 0.0) then
+ wt_sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ else
+ sum = 0.0
+ do i=1,4
+ do j=1,4
+ sum = sum + weights(i,j) * array(ifx+3-i, ify+3-j, izz)
+ end do
+ end do
+ wt_sixteen_pt_average = sum / sum_weight
+ end if
+
+ end function wt_sixteen_pt_average
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: four_pt
+ !
+ ! Purpose: Bilinear interpolation among four grid values
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ recursive function four_pt(xx, yy, izz, array, start_x, end_x, start_y, &
+ end_y, start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: start_x, start_y, start_z
+ integer, intent(in) :: end_x, end_y, end_z
+ integer, intent(in) :: izz ! The z-index of the 2d-array to
+ ! interpolate within
+ real (kind=RKIND), intent(in) :: xx , yy ! The location to interpolate to
+ real (kind=RKIND), intent(in) :: msgval
+ real (kind=RKIND), intent(in), optional :: maskval
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
+ intent(in) :: array
+ integer, dimension(:), intent(in) :: interp_list
+ integer, intent(in) :: idx
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y), &
+ intent(in), optional :: mask_array
+
+ ! Return value
+ real (kind=RKIND) :: four_pt
+
+ ! Local variables
+ integer :: min_x, min_y, max_x, max_y
+
+ min_x = floor(xx)
+ min_y = floor(yy)
+ max_x = ceiling(xx)
+ max_y = ceiling(yy)
+
+ if (min_x < start_x .or. max_x > end_x) then
+ four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ return
+ end if
+
+ if (min_y < start_y .or. max_y > end_y) then
+ four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ return
+ end if
+
+ ! If we have a missing value, try the next interpolation method in the sequence
+ if (present(mask_array) .and. present(maskval)) then
+ if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) == maskval .or. &
+ array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) == maskval .or. &
+ array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) == maskval .or. &
+ array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) == maskval) then
+ four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ return
+ end if
+ else
+ if (array(min_x,min_y,izz) == msgval .or. &
+ array(max_x,min_y,izz) == msgval .or. &
+ array(min_x,max_y,izz) == msgval .or. &
+ array(max_x,max_y,izz) == msgval ) then
+ four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx)
+ return
+ end if
+ end if
+
+ if (min_x == max_x) then
+ if (min_y == max_y) then
+ four_pt = array(min_x,min_y,izz)
+ else
+ four_pt = array(min_x,min_y,izz)*(real(max_y)-yy) + &
+ array(min_x,max_y,izz)*(yy-real(min_y))
+ end if
+ else if (min_y == max_y) then
+ if (min_x == max_x) then
+ four_pt = array(min_x,min_y,izz)
+ else
+ four_pt = array(min_x,min_y,izz)*(real(max_x)-xx) + &
+ array(max_x,min_y,izz)*(xx-real(min_x))
+ end if
+ else
+ four_pt = (yy - min_y) * (array(min_x,max_y,izz)*(real(max_x)-xx) + &
+ array(max_x,max_y,izz)*(xx-real(min_x))) + &
+ (max_y - yy) * (array(min_x,min_y,izz)*(real(max_x)-xx) + &
+ array(max_x,min_y,izz)*(xx-real(min_x)));
+ end if
+
+ end function four_pt
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: sixteen_pt
+ !
+ ! Purpose: Overlapping parabolic interpolation among sixteen grid values
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ recursive function sixteen_pt(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+
+ implicit none
+
+ ! Arguments
+ integer, intent(in) :: izz ! z-index of 2d-array to interpolate within
+ integer, intent(in) :: start_x, start_y, start_z
+ integer, intent(in) :: end_x, end_y, end_z
+ real (kind=RKIND), intent(in) :: xx , yy ! The location to interpolate to
+ real (kind=RKIND), intent(in) :: msgval
+ real (kind=RKIND), intent(in), optional :: maskval
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
+ intent(in) :: array
+ integer, dimension(:), intent(in) :: interp_list
+ integer, intent(in) :: idx
+ real (kind=RKIND), dimension(start_x:end_x, start_y:end_y), &
+ intent(in), optional :: mask_array
+
+ ! Return value
+ real (kind=RKIND) :: sixteen_pt
+
+ ! Local variables
+ integer :: n , i , j , k , kk , l , ll
+ real (kind=RKIND) :: x , y , a , b , c , d , e , f , g , h
+ real (kind=RKIND), dimension(4,4) :: stl
+ logical :: is_masked
+
+ is_masked = .false.
+
+ if (int(xx) < start_x .or. int(xx) > end_x .or. &
+ int(yy) < start_y .or. int(yy) > end_y) then
+ sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ return
+ end if
+
+ sixteen_pt = 0.0
+ n = 0
+ i = int(xx + 0.00001)
+ j = int(yy + 0.00001)
+ x = xx - i
+ y = yy - j
+
+ if ( ( abs(x) > 0.0001 ) .or. ( abs(y) > 0.0001 ) ) then
+
+ loop_1 : do k = 1,4
+ kk = i + k - 2
+ if ( kk < start_x) then
+ kk = start_x
+ else if ( kk > end_x) then
+ kk = end_x
+ end if
+ loop_2 : do l = 1,4
+ stl(k,l) = 0.
+ ll = j + l - 2
+ if ( ll < start_y ) then
+ ll = start_y
+ else if ( ll > end_y) then
+ ll = end_y
+ end if
+ stl(k,l) = array(kk,ll,izz)
+ n = n + 1
+ if (present(mask_array) .and. present(maskval)) then
+ if (mask_array(kk,ll) == maskval) is_masked = .true.
+ end if
+ if ( stl(k,l) == 0. .and. msgval /= 0.) then
+ stl(k,l) = 1.E-20
+ end if
+ end do loop_2
+ end do loop_1
+
+ ! If we have a missing value, try the next interpolation method in the sequence
+ if (present(mask_array) .and. present(maskval)) then
+ do k=1,4
+ do l=1,4
+ if (stl(k,l) == msgval .or. is_masked) then
+ sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ return
+ end if
+ end do
+ end do
+ else
+ do k=1,4
+ do l=1,4
+ if (stl(k,l) == msgval) then
+ sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx)
+ return
+ end if
+ end do
+ end do
+ end if
+
+ a = oned(x,stl(1,1),stl(2,1),stl(3,1),stl(4,1))
+ b = oned(x,stl(1,2),stl(2,2),stl(3,2),stl(4,2))
+ c = oned(x,stl(1,3),stl(2,3),stl(3,3),stl(4,3))
+ d = oned(x,stl(1,4),stl(2,4),stl(3,4),stl(4,4))
+ sixteen_pt = oned(y,a,b,c,d)
+
+ if (n /= 16) then
+ e = oned(y,stl(1,1),stl(1,2),stl(1,3),stl(1,4))
+ f = oned(y,stl(2,1),stl(2,2),stl(2,3),stl(2,4))
+ g = oned(y,stl(3,1),stl(3,2),stl(3,3),stl(3,4))
+ h = oned(y,stl(4,1),stl(4,2),stl(4,3),stl(4,4))
+ sixteen_pt = (sixteen_pt+oned(x,e,f,g,h)) * 0.5
+ end if
+
+ if (sixteen_pt == 1.E-20) sixteen_pt = 0.
+
+ else
+ if (present(mask_array) .and. present(maskval)) then
+ if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
+ mask_array(i,j) /= maskval .and. array(i,j,izz) /= msgval) then
+ sixteen_pt = array(i,j,izz)
+ else
+ sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ end if
+ else
+ if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. array(i,j,izz) /= msgval) then
+ sixteen_pt = array(i,j,izz)
+ else
+ sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
+ start_z, end_z, msgval, interp_list, idx, maskval, mask_array)
+ end if
+ end if
+ end if
+
+ end function sixteen_pt
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: oned
+ !
+ ! Purpose: 1-dimensional overlapping parabolic interpolation
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ function oned(x,a,b,c,d)
+
+ implicit none
+
+ ! Arguments
+ real (kind=RKIND), intent(in) :: x,a,b,c,d
+
+ ! Return value
+ real (kind=RKIND) :: oned
+
+ oned = 0.
+
+ if ( x == 0. ) then
+ oned = b
+ else if ( x == 1. ) then
+ oned = c
+ end if
+
+ if (b*c /= 0.) then
+ if ( a*d == 0. ) then
+ if ( ( a == 0 ) .and. ( d == 0 ) ) then
+ oned = b*(1.0-x)+c*x
+ else if ( a /= 0. ) then
+ oned = b+x*(0.5*(c-a)+x*(0.5*(c+a)-b))
+ else if ( d /= 0. ) then
+ oned = c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c))
+ end if
+ else
+ oned = (1.0-x)*(b+x*(0.5*(c-a)+x*(0.5*(c+a)-b)))+x*(c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c)))
+ end if
+ end if
+
+ end function oned
+
+end module hinterp
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/module_llxy.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_llxy.F 2011-02-10 22:05:23 UTC (rev 737)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_llxy.F 2011-02-11 21:26:40 UTC (rev 738)
@@ -132,22 +132,22 @@
INTEGER, PARAMETER :: HH=4, VV=5
- REAL (KIND=4), PARAMETER :: PI = 3.141592653589793
+ REAL (KIND=RKIND), PARAMETER :: PI = 3.141592653589793
- REAL (KIND=4), PARAMETER :: DEG_PER_RAD = 180./PI
- REAL (KIND=4), PARAMETER :: RAD_PER_DEG = PI/180.
+ REAL (KIND=RKIND), PARAMETER :: DEG_PER_RAD = 180./PI
+ REAL (KIND=RKIND), PARAMETER :: RAD_PER_DEG = PI/180.
- REAL (KIND=4), PARAMETER :: A_WGS84 = 6378137.
- REAL (KIND=4), PARAMETER :: B_WGS84 = 6356752.314
- REAL (KIND=4), PARAMETER :: RE_WGS84 = A_WGS84
- REAL (KIND=4), PARAMETER :: E_WGS84 = 0.081819192
+ REAL (KIND=RKIND), PARAMETER :: A_WGS84 = 6378137.
+ REAL (KIND=RKIND), PARAMETER :: B_WGS84 = 6356752.314
+ REAL (KIND=RKIND), PARAMETER :: RE_WGS84 = A_WGS84
+ REAL (KIND=RKIND), PARAMETER :: E_WGS84 = 0.081819192
- REAL (KIND=4), PARAMETER :: A_NAD83 = 6378137.
- REAL (KIND=4), PARAMETER :: RE_NAD83 = A_NAD83
- REAL (KIND=4), PARAMETER :: E_NAD83 = 0.0818187034
+ REAL (KIND=RKIND), PARAMETER :: A_NAD83 = 6378137.
+ REAL (KIND=RKIND), PARAMETER :: RE_NAD83 = A_NAD83
+ REAL (KIND=RKIND), PARAMETER :: E_NAD83 = 0.0818187034
- REAL (KIND=4), PARAMETER :: EARTH_RADIUS_M = 6370000.
- REAL (KIND=4), PARAMETER :: EARTH_CIRC_M = 2.*PI*EARTH_RADIUS_M
+ REAL (KIND=RKIND), PARAMETER :: EARTH_RADIUS_M = 6370000.
+ REAL (KIND=RKIND), PARAMETER :: EARTH_CIRC_M = 2.*PI*EARTH_RADIUS_M
INTEGER, PUBLIC, PARAMETER :: PROJ_LATLON = 0
INTEGER, PUBLIC, PARAMETER :: PROJ_LC = 1
@@ -174,43 +174,43 @@
! in an odd row
INTEGER :: jydim ! For Rotated Lat/Lon -- number of rows
INTEGER :: stagger ! For Rotated Lat/Lon -- mass or velocity grid
- REAL (KIND=4) :: phi ! For Rotated Lat/Lon -- domain half-extent in
+ REAL (KIND=RKIND) :: phi ! For Rotated Lat/Lon -- domain half-extent in
! degrees latitude
- REAL (KIND=4) :: lambda ! For Rotated Lat/Lon -- domain half-extend in
+ REAL (KIND=RKIND) :: lambda ! For Rotated Lat/Lon -- domain half-extend in
! degrees longitude
- REAL (KIND=4) :: lat1 ! SW latitude (1,1) in degrees (-90->90N)
- REAL (KIND=4) :: lon1 ! SW longitude (1,1) in degrees (-180->180E)
- REAL (KIND=4) :: lat0 ! For Cassini, latitude of projection pole
- REAL (KIND=4) :: lon0 ! For Cassini, longitude of projection pole
- REAL (KIND=4) :: dx ! Grid spacing in meters at truelats, used
+ REAL (KIND=RKIND) :: lat1 ! SW latitude (1,1) in degrees (-90->90N)
+ REAL (KIND=RKIND) :: lon1 ! SW longitude (1,1) in degrees (-180->180E)
+ REAL (KIND=RKIND) :: lat0 ! For Cassini, latitude of projection pole
+ REAL (KIND=RKIND) :: lon0 ! For Cassini, longitude of projection pole
+ REAL (KIND=RKIND) :: dx ! Grid spacing in meters at truelats, used
! only for ps, lc, and merc projections
- REAL (KIND=4) :: dy ! Grid spacing in meters at truelats, used
+ REAL (KIND=RKIND) :: dy ! Grid spacing in meters at truelats, used
! only for ps, lc, and merc projections
- REAL (KIND=4) :: latinc ! Latitude increment for cylindrical lat/lon
- REAL (KIND=4) :: loninc ! Longitude increment for cylindrical lat/lon
+ REAL (KIND=RKIND) :: latinc ! Latitude increment for cylindrical lat/lon
+ REAL (KIND=RKIND) :: loninc ! Longitude increment for cylindrical lat/lon
! also the lon increment for Gaussian grid
- REAL (KIND=4) :: dlat ! Lat increment for lat/lon grids
- REAL (KIND=4) :: dlon ! Lon increment for lat/lon grids
- REAL (KIND=4) :: stdlon ! Longitude parallel to y-axis (-180->180E)
- REAL (KIND=4) :: truelat1 ! First true latitude (all projections)
- REAL (KIND=4) :: truelat2 ! Second true lat (LC only)
- REAL (KIND=4) :: hemi ! 1 for NH, -1 for SH
- REAL (KIND=4) :: cone ! Cone factor for LC projections
- REAL (KIND=4) :: polei ! Computed i-location of pole point
- REAL (KIND=4) :: polej ! Computed j-location of pole point
- REAL (KIND=4) :: rsw ! Computed radius to SW corner
- REAL (KIND=4) :: rebydx ! Earth radius divided by dx
- REAL (KIND=4) :: knowni ! X-location of known lat/lon
- REAL (KIND=4) :: knownj ! Y-location of known lat/lon
- REAL (KIND=4) :: re_m ! Radius of spherical earth, meters
- REAL (KIND=4) :: rho0 ! For Albers equal area
- REAL (KIND=4) :: nc ! For Albers equal area
- REAL (KIND=4) :: bigc ! For Albers equal area
+ REAL (KIND=RKIND) :: dlat ! Lat increment for lat/lon grids
+ REAL (KIND=RKIND) :: dlon ! Lon increment for lat/lon grids
+ REAL (KIND=RKIND) :: stdlon ! Longitude parallel to y-axis (-180->180E)
+ REAL (KIND=RKIND) :: truelat1 ! First true latitude (all projections)
+ REAL (KIND=RKIND) :: truelat2 ! Second true lat (LC only)
+ REAL (KIND=RKIND) :: hemi ! 1 for NH, -1 for SH
+ REAL (KIND=RKIND) :: cone ! Cone factor for LC projections
+ REAL (KIND=RKIND) :: polei ! Computed i-location of pole point
+ REAL (KIND=RKIND) :: polej ! Computed j-location of pole point
+ REAL (KIND=RKIND) :: rsw ! Computed radius to SW corner
+ REAL (KIND=RKIND) :: rebydx ! Earth radius divided by dx
+ REAL (KIND=RKIND) :: knowni ! X-location of known lat/lon
+ REAL (KIND=RKIND) :: knownj ! Y-location of known lat/lon
+ REAL (KIND=RKIND) :: re_m ! Radius of spherical earth, meters
+ REAL (KIND=RKIND) :: rho0 ! For Albers equal area
+ REAL (KIND=RKIND) :: nc ! For Albers equal area
+ REAL (KIND=RKIND) :: bigc ! For Albers equal area
LOGICAL :: init ! Flag to indicate if this struct is
! ready for use
LOGICAL :: wrap ! For Gaussian -- flag to indicate wrapping
! around globe?
- REAL (KIND=4), POINTER, DIMENSION(:) :: gauss_lat ! Latitude array for Gaussian grid
+ REAL (KIND=RKIND), POINTER, DIMENSION(:) :: gauss_lat ! Latitude array for Gaussian grid
END TYPE proj_info
@@ -282,27 +282,27 @@
INTEGER, INTENT(IN), OPTIONAL :: ixdim
INTEGER, INTENT(IN), OPTIONAL :: jydim
INTEGER, INTENT(IN), OPTIONAL :: stagger
- REAL (KIND=4), INTENT(IN), OPTIONAL :: latinc
- REAL (KIND=4), INTENT(IN), OPTIONAL :: loninc
- REAL (KIND=4), INTENT(IN), OPTIONAL :: lat1
- REAL (KIND=4), INTENT(IN), OPTIONAL :: lon1
- REAL (KIND=4), INTENT(IN), OPTIONAL :: lat0
- REAL (KIND=4), INTENT(IN), OPTIONAL :: lon0
- REAL (KIND=4), INTENT(IN), OPTIONAL :: dx
- REAL (KIND=4), INTENT(IN), OPTIONAL :: stdlon
- REAL (KIND=4), INTENT(IN), OPTIONAL :: truelat1
- REAL (KIND=4), INTENT(IN), OPTIONAL :: truelat2
- REAL (KIND=4), INTENT(IN), OPTIONAL :: knowni
- REAL (KIND=4), INTENT(IN), OPTIONAL :: knownj
- REAL (KIND=4), INTENT(IN), OPTIONAL :: phi
- REAL (KIND=4), INTENT(IN), OPTIONAL :: lambda
- REAL (KIND=4), INTENT(IN), OPTIONAL :: r_earth
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: latinc
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: loninc
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: lat1
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: lon1
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: lat0
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: lon0
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: dx
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: stdlon
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: truelat1
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: truelat2
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: knowni
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: knownj
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: phi
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: lambda
+ REAL (KIND=RKIND), INTENT(IN), OPTIONAL :: r_earth
TYPE(proj_info), INTENT(OUT) :: proj
INTEGER :: iter
- REAL (KIND=4) :: dummy_lon1
- REAL (KIND=4) :: dummy_lon0
- REAL (KIND=4) :: dummy_stdlon
+ REAL (KIND=RKIND) :: dummy_lon1
+ REAL (KIND=RKIND) :: dummy_lon0
+ REAL (KIND=RKIND) :: dummy_stdlon
! First, verify that mandatory parameters are present for the specified proj_code
IF ( proj_code == PROJ_LC ) THEN
@@ -583,10 +583,10 @@
IMPLICIT NONE
TYPE(proj_info), INTENT(IN) :: proj
- REAL (KIND=4), INTENT(IN) :: lat
- REAL (KIND=4), INTENT(IN) :: lon
- REAL (KIND=4), INTENT(OUT) :: i
- REAL (KIND=4), INTENT(OUT) :: j
+ REAL (KIND=RKIND), INTENT(IN) :: lat
+ REAL (KIND=RKIND), INTENT(IN) :: lon
+ REAL (KIND=RKIND), INTENT(OUT) :: i
+ REAL (KIND=RKIND), INTENT(OUT) :: j
IF (.NOT.proj%init) THEN
PRINT '(A)', 'You have not called map_set for this projection'
@@ -642,10 +642,10 @@
IMPLICIT NONE
TYPE(proj_info),INTENT(IN) :: proj
- REAL (KIND=4), INTENT(IN) :: i
- REAL (KIND=4), INTENT(IN) :: j
- REAL (KIND=4), INTENT(OUT) :: lat
- REAL (KIND=4), INTENT(OUT) :: lon
+ REAL (KIND=RKIND), INTENT(IN) :: i
+ REAL (KIND=RKIND), INTENT(IN) :: j
+ REAL (KIND=RKIND), INTENT(OUT) :: lat
+ REAL (KIND=RKIND), INTENT(OUT) :: lon
IF (.NOT.proj%init) THEN
PRINT '(A)', 'You have not called map_set for this projection'
@@ -700,10 +700,10 @@
TYPE(proj_info), INTENT(INOUT) :: proj
! Local vars
- REAL (KIND=4) :: ala1
- REAL (KIND=4) :: alo1
- REAL (KIND=4) :: reflon
- REAL (KIND=4) :: scale_top
+ REAL (KIND=RKIND) :: ala1
+ REAL (KIND=RKIND) :: alo1
+ REAL (KIND=RKIND) :: reflon
+ REAL (KIND=RKIND) :: scale_top
! Executable code
reflon = proj%stdlon + 90.
@@ -734,21 +734,21 @@
IMPLICIT NONE
! Delcare input arguments
- REAL (KIND=4), INTENT(IN) :: lat
- REAL (KIND=4), INTENT(IN) :: lon
+ REAL (KIND=RKIND), INTENT(IN) :: lat
+ REAL (KIND=RKIND), INTENT(IN) :: lon
TYPE(proj_info),INTENT(IN) :: proj
! Declare output arguments
- REAL (KIND=4), INTENT(OUT) :: i !(x-index)
- REAL (KIND=4), INTENT(OUT) :: j !(y-index)
+ REAL (KIND=RKIND), INTENT(OUT) :: i !(x-index)
+ REAL (KIND=RKIND), INTENT(OUT) :: j !(y-index)
! Declare local variables
- REAL (KIND=4) :: reflon
- REAL (KIND=4) :: scale_top
- REAL (KIND=4) :: ala
- REAL (KIND=4) :: alo
- REAL (KIND=4) :: rm
+ REAL (KIND=RKIND) :: reflon
+ REAL (KIND=RKIND) :: scale_top
+ REAL (KIND=RKIND) :: ala
+ REAL (KIND=RKIND) :: alo
+ REAL (KIND=RKIND) :: rm
! BEGIN CODE
@@ -779,20 +779,20 @@
IMPLICIT NONE
! Declare input arguments
- REAL (KIND=4), INTENT(IN) :: i ! Column
- REAL (KIND=4), INTENT(IN) :: j ! Row
+ REAL (KIND=RKIND), INTENT(IN) :: i ! Column
+ REAL (KIND=RKIND), INTENT(IN) :: j ! Row
TYPE (proj_info), INTENT(IN) :: proj
! Declare output arguments
- REAL (KIND=4), INTENT(OUT) :: lat ! -90 -> 90 north
- REAL (KIND=4), INTENT(OUT) :: lon ! -180 -> 180 East
+ REAL (KIND=RKIND), INTENT(OUT) :: lat ! -90 -> 90 north
+ REAL (KIND=RKIND), INTENT(OUT) :: lon ! -180 -> 180 East
! Local variables
- REAL (KIND=4) :: reflon
- REAL (KIND=4) :: scale_top
- REAL (KIND=4) :: xx,yy
- REAL (KIND=4) :: gi2, r2
- REAL (KIND=4) :: arccos
+ REAL (KIND=RKIND) :: reflon
+ REAL (KIND=RKIND) :: scale_top
+ REAL (KIND=RKIND) :: xx,yy
+ REAL (KIND=RKIND) :: gi2, r2
+ REAL (KIND=RKIND) :: arccos
! Begin Code
@@ -844,7 +844,7 @@
TYPE(proj_info), INTENT(INOUT) :: proj
! Local variables
- REAL (KIND=4) :: h, mc, tc, t, rho
+ REAL (KIND=RKIND) :: h, mc, tc, t, rho
h = proj%hemi
@@ -873,14 +873,14 @@
IMPLICIT NONE
! Arguments
- REAL (KIND=4), INTENT(IN) :: lat
- REAL (KIND=4), INTENT(IN) :: lon
- REAL (KIND=4), INTENT(OUT) :: i !(x-index)
- REAL (KIND=4), INTENT(OUT) :: j !(y-index)
+ REAL (KIND=RKIND), INTENT(IN) :: lat
+ REAL (KIND=RKIND), INTENT(IN) :: lon
+ REAL (KIND=RKIND), INTENT(OUT) :: i !(x-index)
+ REAL (KIND=RKIND), INTENT(OUT) :: j !(y-index)
TYPE(proj_info),INTENT(IN) :: proj
! Local variables
- REAL (KIND=4) :: h, mc, tc, t, rho
+ REAL (KIND=RKIND) :: h, mc, tc, t, rho
h = proj%hemi
@@ -914,15 +914,15 @@
implicit none
! Arguments
- REAL (KIND=4), INTENT(IN) :: i ! Column
- REAL (KIND=4), INTENT(IN) :: j ! Row
- REAL (KIND=4), INTENT(OUT) :: lat ! -90 -> 90 north
- REAL (KIND=4), INTENT(OUT) :: lon ! -180 -> 180 East
+ REAL (KIND=RKIND), INTENT(IN) :: i ! Column
+ REAL (KIND=RKIND), INTENT(IN) :: j ! Row
+ REAL (KIND=RKIND), INTENT(OUT) :: lat ! -90 -> 90 north
+ REAL (KIND=RKIND), INTENT(OUT) :: lon ! -180 -> 180 East
TYPE (proj_info), INTENT(IN) :: proj
! Local variables
- REAL (KIND=4) :: h, mc, tc, t, rho, x, y
- REAL (KIND=4) :: chi, a, b, c, d
+ REAL (KIND=RKIND) :: h, mc, tc, t, rho, x, y
+ REAL (KIND=RKIND) :: chi, a, b, c, d
h = proj%hemi
x = (i - proj%knowni + proj%polei)
@@ -966,7 +966,7 @@
TYPE(proj_info), INTENT(INOUT) :: proj
! Local variables
- REAL (KIND=4) :: h, m1, m2, q1, q2, theta, q, sinphi
+ REAL (KIND=RKIND) :: h, m1, m2, q1, q2, theta, q, sinphi
h = proj%hemi
@@ -1014,14 +1014,14 @@
IMPLICIT NONE
! Arguments
- REAL (KIND=4), INTENT(IN) :: lat
- REAL (KIND=4), INTENT(IN) :: lon
- REAL (KIND=4), INTENT(OUT) :: i !(x-index)
- REAL (KIND=4), INTENT(OUT) :: j !(y-index)
+ REAL (KIND=RKIND), INTENT(IN) :: lat
+ REAL (KIND=RKIND), INTENT(IN) :: lon
+ REAL (KIND=RKIND), INTENT(OUT) :: i !(x-index)
+ REAL (KIND=RKIND), INTENT(OUT) :: j !(y-index)
TYPE(proj_info),INTENT(IN) :: proj
! Local variables
- REAL (KIND=4) :: h, q, rho, theta, sinphi
+ REAL (KIND=RKIND) :: h, q, rho, theta, sinphi
h = proj%hemi
@@ -1055,15 +1055,15 @@
implicit none
! Arguments
- REAL (KIND=4), INTENT(IN) :: i ! Column
- REAL (KIND=4), INTENT(IN) :: j ! Row
- REAL (KIND=4), INTENT(OUT) :: lat ! -90 -> 90 north
- REAL (KIND=4), INTENT(OUT) :: lon ! -180 -> 180 East
+ REAL (KIND=RKIND), INTENT(IN) :: i ! Column
+ REAL (KIND=RKIND), INTENT(IN) :: j ! Row
+ REAL (KIND=RKIND), INTENT(OUT) :: lat ! -90 -> 90 north
+ REAL (KIND=RKIND), INTENT(OUT) :: lon ! -180 -> 180 East
TYPE (proj_info), INTENT(IN) :: proj
! Local variables
- REAL (KIND=4) :: h, q, rho, theta, beta, x, y
- REAL (KIND=4) :: a, b, c
+ REAL (KIND=RKIND) :: h, q, rho, theta, beta, x, y
+ REAL (KIND=RKIND) :: a, b, c
h = proj%hemi
@@ -1098,10 +1098,10 @@
TYPE(proj_info), INTENT(INOUT) :: proj
- REAL (KIND=4) :: arg
- REAL (KIND=4) :: deltalon1
- REAL (KIND=4) :: tl1r
- REAL (KIND=4) :: ctl1r
+ REAL (KIND=RKIND) :: arg
+ REAL (KIND=RKIND) :: deltalon1
+ REAL (KIND=RKIND) :: tl1r
+ REAL (KIND=RKIND) :: ctl1r
! Compute cone factor
CALL lc_cone(proj%truelat1, proj%truelat2, proj%cone)
@@ -1138,11 +1138,11 @@
IMPLICIT NONE
! Input Args
- REAL (KIND=4), INTENT(IN) :: truelat1 ! (-90 -> 90 degrees N)
- REAL (KIND=4), INTENT(IN) :: truelat2 ! " " " " "
+ REAL (KIND=RKIND), INTENT(IN) :: truelat1 ! (-90 -> 90 degrees N)
+ REAL (KIND=RKIND), INTENT(IN) :: truelat2 ! " " " " "
! Output Args
- REAL (KIND=4), INTENT(OUT) :: cone
+ REAL (KIND=RKIND), INTENT(OUT) :: cone
! Locals
@@ -1178,22 +1178,22 @@
IMPLICIT NONE
! Input Args
- REAL (KIND=4), INTENT(IN) :: i ! Cartesian X coordinate
- REAL (KIND=4), INTENT(IN) :: j ! Cartesian Y coordinate
+ REAL (KIND=RKIND), INTENT(IN) :: i ! Cartesian X coordinate
+ REAL (KIND=RKIND), INTENT(IN) :: j ! Cartesian Y coordinate
TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure
! Output Args
- REAL (KIND=4), INTENT(OUT) :: lat ! Latitude (-90->90 deg N)
- REAL (KIND=4), INTENT(OUT) :: lon ! Longitude (-180->180 E)
+ REAL (KIND=RKIND), INTENT(OUT) :: lat ! Latitude (-90->90 deg N)
+ REAL (KIND=RKIND), INTENT(OUT) :: lon ! Longitude (-180->180 E)
! Locals
- REAL (KIND=4) :: inew
- REAL (KIND=4) :: jnew
- REAL (KIND=4) :: r
- REAL (KIND=4) :: chi,chi1,chi2
- REAL (KIND=4) :: r2
- REAL (KIND=4) :: xx
- REAL (KIND=4) :: yy
+ REAL (KIND=RKIND) :: inew
+ REAL (KIND=RKIND) :: jnew
+ REAL (KIND=RKIND) :: r
+ REAL (KIND=RKIND) :: chi,chi1,chi2
+ REAL (KIND=RKIND) :: r2
+ REAL (KIND=RKIND) :: xx
+ REAL (KIND=RKIND) :: yy
! BEGIN CODE
@@ -1255,20 +1255,20 @@
IMPLICIT NONE
! Input Args
- REAL (KIND=4), INTENT(IN) :: lat ! Latitude (-90->90 deg N)
- REAL (KIND=4), INTENT(IN) :: lon ! Longitude (-180->180 E)
+ REAL (KIND=RKIND), INTENT(IN) :: lat ! Latitude (-90->90 deg N)
+ REAL (KIND=RKIND), INTENT(IN) :: lon ! Longitude (-180->180 E)
TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure
! Output Args
- REAL (KIND=4), INTENT(OUT) :: i ! Cartesian X coordinate
- REAL (KIND=4), INTENT(OUT) :: j ! Cartesian Y coordinate
+ REAL (KIND=RKIND), INTENT(OUT) :: i ! Cartesian X coordinate
+ REAL (KIND=RKIND), INTENT(OUT) :: j ! Cartesian Y coordinate
! Locals
- REAL (KIND=4) :: arg
- REAL (KIND=4) :: deltalon
- REAL (KIND=4) :: tl1r
- REAL (KIND=4) :: rm
- REAL (KIND=4) :: ctl1r
+ REAL (KIND=RKIND) :: arg
+ REAL (KIND=RKIND) :: deltalon
+ REAL (KIND=RKIND) :: tl1r
+ REAL (KIND=RKIND) :: rm
+ REAL (KIND=RKIND) :: ctl1r
! BEGIN CODE
@@ -1310,7 +1310,7 @@
IMPLICIT NONE
TYPE(proj_info), INTENT(INOUT) :: proj
- REAL (KIND=4) :: clain
+ REAL (KIND=RKIND) :: clain
! Preliminary variables
@@ -1336,12 +1336,12 @@
! Compute i/j coordinate from lat lon for mercator projection
IMPLICIT NONE
- REAL (KIND=4), INTENT(IN) :: lat
- REAL (KIND=4), INTENT(IN) :: lon
+ REAL (KIND=RKIND), INTENT(IN) :: lat
+ REAL (KIND=RKIND), INTENT(IN) :: lon
TYPE(proj_info),INTENT(IN) :: proj
- REAL (KIND=4),INTENT(OUT) :: i
- REAL (KIND=4),INTENT(OUT) :: j
- REAL (KIND=4) :: deltalon
+ REAL (KIND=RKIND),INTENT(OUT) :: i
+ REAL (KIND=RKIND),INTENT(OUT) :: j
+ REAL (KIND=RKIND) :: deltalon
deltalon = lon - proj%lon1
IF (deltalon .LT. -180.) deltalon = deltalon + 360.
@@ -1360,11 +1360,11 @@
! Compute the lat/lon from i/j for mercator projection
IMPLICIT NONE
- REAL (KIND=4),INTENT(IN) :: i
- REAL (KIND=4),INTENT(IN) :: j
+ REAL (KIND=RKIND),INTENT(IN) :: i
+ REAL (KIND=RKIND),INTENT(IN) :: j
TYPE(proj_info),INTENT(IN) :: proj
- REAL (KIND=4), INTENT(OUT) :: lat
- REAL (KIND=4), INTENT(OUT) :: lon
+ REAL (KIND=RKIND), INTENT(OUT) :: lat
+ REAL (KIND=RKIND), INTENT(OUT) :: lon
lat = 2.0*ATAN(EXP(proj%dlon*(proj%rsw + j-proj%knownj)))*deg_per_rad - 90.
@@ -1380,14 +1380,14 @@
! Compute the i/j location of a lat/lon on a LATLON grid.
IMPLICIT NONE
- REAL (KIND=4), INTENT(IN) :: lat
- REAL (KIND=4), INTENT(IN) :: lon
+ REAL (KIND=RKIND), INTENT(IN) :: lat
+ REAL (KIND=RKIND), INTENT(IN) :: lon
TYPE(proj_info), INTENT(IN) :: proj
- REAL (KIND=4), INTENT(OUT) :: i
- REAL (KIND=4), INTENT(OUT) :: j
+ REAL (KIND=RKIND), INTENT(OUT) :: i
+ REAL (KIND=RKIND), INTENT(OUT) :: j
- REAL (KIND=4) :: deltalat
- REAL (KIND=4) :: deltalon
+ REAL (KIND=RKIND) :: deltalat
+ REAL (KIND=RKIND) :: deltalon
! Compute deltalat and deltalon as the difference between the input
! lat/lon and the origin lat/lon
@@ -1410,15 +1410,15 @@
! Compute the lat/lon location of an i/j on a LATLON grid.
IMPLICIT NONE
- REAL (KIND=4), INTENT(IN) :: i
- REAL (KIND=4), INTENT(IN) :: j
+ REAL (KIND=RKIND), INTENT(IN) :: i
+ REAL (KIND=RKIND), INTENT(IN) :: j
TYPE(proj_info), INTENT(IN) :: proj
- REAL (KIND=4), INTENT(OUT) :: lat
- REAL (KIND=4), INTENT(OUT) :: lon
+ REAL (KIND=RKIND), INTENT(OUT) :: lat
+ REAL (KIND=RKIND), INTENT(OUT) :: lon
- REAL (KIND=4) :: i_work, j_work
- REAL (KIND=4) :: deltalat
- REAL (KIND=4) :: deltalon
+ REAL (KIND=RKIND) :: i_work, j_work
+ REAL (KIND=RKIND) :: deltalat
+ REAL (KIND=RKIND) :: deltalon
i_work = i - proj%knowni
j_work = j - proj%knownj
@@ -1452,13 +1452,13 @@
implicit none
! Arguments
- REAL (KIND=4), intent(in) :: lat, lon
- REAL (KIND=4), intent(out) :: i, j
+ REAL (KIND=RKIND), intent(in) :: lat, lon
+ REAL (KIND=RKIND), intent(out) :: i, j
type(proj_info), intent(in) :: proj
! Local variables
- REAL (KIND=4) :: deltalat
- REAL (KIND=4) :: deltalon
+ REAL (KIND=RKIND) :: deltalat
+ REAL (KIND=RKIND) :: deltalon
! Compute deltalat and deltalon as the difference between the input
! lat/lon and the origin lat/lon
@@ -1487,14 +1487,14 @@
implicit none
! Arguments
- REAL (KIND=4), intent(in) :: i, j
- REAL (KIND=4), intent(out) :: lat, lon
+ REAL (KIND=RKIND), intent(in) :: i, j
+ REAL (KIND=RKIND), intent(out) :: lat, lon
type(proj_info), intent(in) :: proj
! Local variables
- REAL (KIND=4) :: deltalat
- REAL (KIND=4) :: deltalon
- REAL (KIND=4) :: i_work, j_work
+ REAL (KIND=RKIND) :: deltalat
+ REAL (KIND=RKIND) :: deltalon
+ REAL (KIND=RKIND) :: i_work, j_work
i_work = i - proj%knowni
j_work = j - proj%knownj
@@ -1524,7 +1524,7 @@
type(proj_info), intent(inout) :: proj
! Local variables
- REAL (KIND=4) :: comp_lat, comp_lon
+ REAL (KIND=RKIND) :: comp_lat, comp_lon
logical :: global_domain
proj%hemi = 1.0
@@ -1551,12 +1551,12 @@
implicit none
! Arguments
- REAL (KIND=4), intent(in) :: lat, lon
- REAL (KIND=4), intent(out) :: i, j
+ REAL (KIND=RKIND), intent(in) :: lat, lon
+ REAL (KIND=RKIND), intent(out) :: i, j
type(proj_info), intent(in) :: proj
! Local variables
- REAL (KIND=4) :: comp_lat, comp_lon
+ REAL (KIND=RKIND) :: comp_lat, comp_lon
! Convert geographic to computational lat/lon
if (abs(proj%lat0) /= 90.) then
@@ -1577,12 +1577,12 @@
implicit none
! Arguments
- REAL (KIND=4), intent(in) :: i, j
- REAL (KIND=4), intent(out) :: lat, lon
+ REAL (KIND=RKIND), intent(in) :: i, j
+ REAL (KIND=RKIND), intent(out) :: lat, lon
type(proj_info), intent(in) :: proj
! Local variables
- REAL (KIND=4) :: comp_lat, comp_lon
+ REAL (KIND=RKIND) :: comp_lat, comp_lon
! Convert i/j to computational lat/lon
call ijll_cyl(i, j, proj, comp_lat, comp_lon)
@@ -1607,16 +1607,16 @@
IMPLICIT NONE
- REAL (KIND=4), INTENT(IN ) :: ilat, ilon
- REAL (KIND=4), INTENT( OUT) :: olat, olon
- REAL (KIND=4), INTENT(IN ) :: lat_np, lon_np, lon_0
+ REAL (KIND=RKIND), INTENT(IN ) :: ilat, ilon
+ REAL (KIND=RKIND), INTENT( OUT) :: olat, olon
+ REAL (KIND=RKIND), INTENT(IN ) :: lat_np, lon_np, lon_0
INTEGER, INTENT(IN ), OPTIONAL :: direction
! >=0, default : computational -> geographical
! < 0 : geographical -> computational
- REAL (KIND=4) :: rlat, rlon
- REAL (KIND=4) :: phi_np, lam_np, lam_0, dlam
- REAL (KIND=4) :: sinphi, cosphi, coslam, sinlam
+ REAL (KIND=RKIND) :: rlat, rlon
+ REAL (KIND=RKIND) :: phi_np, lam_np, lam_0, dlam
+ REAL (KIND=RKIND) :: sinphi, cosphi, coslam, sinlam
! Convert all angles to radians
phi_np = lat_np * rad_per_deg
@@ -1661,9 +1661,9 @@
IMPLICIT NONE
! Arguments
- REAL (KIND=4), INTENT(IN) :: lat, lon
- REAL (KIND=4) :: i, j
- REAL (KIND=4), INTENT(OUT) :: i_real, j_real
+ REAL (KIND=RKIND), INTENT(IN) :: lat, lon
+ REAL (KIND=RKIND) :: i, j
+ REAL (KIND=RKIND), INTENT(OUT) :: i_real, j_real
TYPE (proj_info), INTENT(IN) :: proj
! Local variables
@@ -1832,18 +1832,18 @@
IMPLICIT NONE
! Arguments
- REAL (KIND=4), INTENT(IN) :: i, j
- REAL (KIND=4), INTENT(OUT) :: lat, lon
+ REAL (KIND=RKIND), INTENT(IN) :: i, j
+ REAL (KIND=RKIND), INTENT(OUT) :: lat, lon
TYPE (proj_info), INTENT(IN) :: proj
! Local variables
INTEGER :: ih,jh
INTEGER :: midcol,midrow,ncol,iadd1,iadd2,imt,jh2,knrow,krem,kv,nrow
- REAL (KIND=4) :: i_work, j_work
- REAL (KIND=4) :: dphd,dlmd !Grid increments, degrees
+ REAL (KIND=RKIND) :: i_work, j_work
+ REAL (KIND=RKIND) :: dphd,dlmd !Grid increments, degrees
REAL(KIND=HIGH) :: arg1,arg2,d2r,fctr,glatr,glatd,glond,pi, &
r2d,tlatd,tlond,tlatr,tlonr,tlm0,tph0
- REAL (KIND=4) :: col
+ REAL (KIND=RKIND) :: col
i_work = i
j_work = j
@@ -1950,7 +1950,7 @@
INTEGER :: nlat , i
REAL (KIND=HIGH) , PARAMETER :: pi = 3.141592653589793
REAL (KIND=HIGH) , DIMENSION(nlat) :: cosc , gwt , sinc , colat , wos2 , lat
- REAL (KIND=4) , DIMENSION(nlat) :: lat_sp
+ REAL (KIND=RKIND) , DIMENSION(nlat) :: lat_sp
CALL lggaus(nlat, cosc, gwt, sinc, colat, wos2)
@@ -1987,7 +1987,7 @@
! Convergence criterion for iteration of cos latitude
- REAL (KIND=4) , PARAMETER :: xlim = 1.0E-14
+ REAL (KIND=RKIND) , PARAMETER :: xlim = 1.0E-14
INTEGER :: nzero, i, j
REAL (KIND=HIGH) :: fi, fi1, a, b, g, gm, gp, gt, delta, c, d
@@ -2134,13 +2134,13 @@
IMPLICIT NONE
- REAL (KIND=4) , INTENT(IN) :: lat, lon
- REAL (KIND=4) , INTENT(OUT) :: i, j
+ REAL (KIND=RKIND) , INTENT(IN) :: lat, lon
+ REAL (KIND=RKIND) , INTENT(OUT) :: i, j
TYPE (proj_info), INTENT(IN) :: proj
INTEGER :: n , n_low
LOGICAL :: found = .FALSE.
- REAL (KIND=4) :: diff_1 , diff_nlat
+ REAL (KIND=RKIND) :: diff_1 , diff_nlat
! The easy one first, get the x location. The calling routine has already made
! sure that the necessary assumptions concerning the sign of the deltalon and the
Added: branches/atmos_physics/src/core_init_nhyd_atmos/module_queue.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_queue.F (rev 0)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_queue.F 2011-02-11 21:26:40 UTC (rev 738)
@@ -0,0 +1,226 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Module: queue_module
+!
+! Description: This module implements a queue of user-defined data types and
+! a set of routines related to the maintenance and manipulation of the queue.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+module queue
+
+ type q_data ! The user-defined datatype to store in the queue
+ real (kind=RKIND) :: lat, lon
+ integer :: x, y
+ end type q_data
+
+ type q_item ! Wrapper for item to be stored in the queue
+ type (q_data) :: data
+ type (q_item), pointer :: next
+ end type q_item
+
+ type queue ! The queue object, defined by a head and tail pointer
+ type (q_item), pointer :: head, tail
+ integer :: length
+ end type queue
+
+
+ contains
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: q_init
+ !
+ ! Purpose: To initialize a queue
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine q_init(q)
+
+ implicit none
+
+ ! Arguments
+ type (queue), intent(inout) :: q
+
+ nullify(q%head)
+ nullify(q%tail)
+ q%length = 0
+
+ end subroutine q_init
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: q_insert
+ !
+ ! Purpose: To insert an item in the tail of the queue
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine q_insert(q, qitem)
+
+ implicit none
+
+ ! Arguments
+ type (queue), intent(inout) :: q
+ type (q_data), intent(in) :: qitem
+
+ ! Local variables
+ type (q_item), pointer :: newitem
+
+ allocate(newitem)
+ newitem%data = qitem
+ nullify(newitem%next)
+ if (.not.associated(q%tail)) then
+ q%head=>newitem
+ q%tail=>newitem
+ else
+ q%tail%next=>newitem
+ q%tail=>newitem
+ end if
+
+ q%length = q%length + 1
+
+ end subroutine q_insert
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: q_isdata
+ !
+ ! Purpose: This function returns FALSE if the queue is empty and TRUE otherwise
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ function q_isdata(q)
+
+ implicit none
+
+ ! Arguments
+ type (queue), intent(in) :: q
+
+ ! Local variables
+ logical :: q_isdata
+
+ q_isdata = .false.
+
+ if (associated(q%head) .and. (q%length >= 1)) then
+ q_isdata = .true.
+ end if
+
+ end function q_isdata
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: q_peek
+ !
+ ! Purpose: To return the item in the head of the queue, without
+ ! actually removing the item
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ function q_peek(q)
+
+ implicit none
+
+ ! Arguments
+ type (queue), intent(in) :: q
+
+ ! Local variables
+ type (q_data) :: q_peek
+
+ if (associated(q%head)) then
+ q_peek = q%head%data
+ else
+ write(0,*) 'ERROR: q_peek(): Trying to peek at an empty queue'
+ end if
+
+ end function q_peek
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: q_length
+ !
+ ! Purpose: To return the number of items currently in the queue
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ function q_length(q)
+
+ implicit none
+
+ ! Arguments
+ type (queue), intent(in) :: q
+
+ ! Local variables
+ ! type (q_item), pointer :: cursor
+ integer :: q_length
+
+ q_length = q%length
+
+ ! USE THE FOLLOWING TO COUNT THE LENGTH BY ACTUALLY TRAVERSING THE LINKED LIST
+ ! REPRESENTATION OF THE QUEUE
+ ! if (associated(q%head)) then
+ ! q_length = q_length + 1
+ ! cursor=>q%head
+ ! do while(associated(cursor%next))
+ ! cursor=>cursor%next
+ ! q_length = q_length + 1
+ ! end do
+ ! end if
+
+ end function q_length
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: q_remove
+ !
+ ! Purpose: To return the item stored at the head of the queue
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ function q_remove(q)
+
+ implicit none
+
+ ! Arguments
+ type (queue), intent(inout) :: q
+
+ ! Local variables
+ type (q_data) :: q_remove
+ type (q_item), pointer :: cursor
+
+ if (associated(q%head)) then
+ if (associated(q%head%next)) then
+ cursor=>q%head%next
+ q_remove = q%head%data
+ deallocate(q%head)
+ q%head=>cursor
+ else
+ q_remove = q%head%data
+ deallocate(q%head)
+ nullify(q%head)
+ nullify(q%tail)
+ end if
+ q%length = q%length - 1
+ else
+ write(0,*) 'ERROR: q_remove(): Trying to remove from an empty queue'
+ end if
+
+ end function q_remove
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Name: q_destroy
+ !
+ ! Purpose: To free all memory allocated by the queue, thus destroying any
+ ! items that have not been removed
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine q_destroy(q)
+
+ implicit none
+
+ ! Arguments
+ type (queue), intent(inout) :: q
+
+ ! Local variables
+ type (q_item), pointer :: cursor
+
+ q%length = 0
+
+ if (associated(q%head)) then
+ do while(associated(q%head%next))
+ cursor=>q%head
+ q%head=>q%head%next
+ deallocate(cursor)
+ end do
+ deallocate(q%head)
+ end if
+
+ end subroutine q_destroy
+
+end module queue
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-02-10 22:05:23 UTC (rev 737)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-02-11 21:26:40 UTC (rev 738)
@@ -1968,6 +1968,7 @@
use read_met
use llxy
+ use hinterp
use dmpar
implicit none
@@ -2006,6 +2007,12 @@
integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
integer :: nInterpPoints, ndims
+ integer, dimension(5) :: interp_list
+ real (kind=RKIND) :: maskval
+ real (kind=RKIND) :: msgval
+ real (kind=RKIND) :: fillval
+ integer :: masked
+
!This is temporary variable here. It just need when calculate tangential velocity v.
integer :: eoe, j
integer, dimension(:), pointer :: nEdgesOnEdge, nEdgesOnCell
@@ -2033,12 +2040,15 @@
real (kind=4) :: scalefactor ! NB: this should be a single-precision real
real (kind=RKIND) :: lat_pt, lon_pt, lon_pt_o
real (kind=4), allocatable, dimension(:,:,:) :: rarray ! NB: this should be a single-precision real array
- real (kind=4), allocatable, dimension(:,:) :: maxsnowalb ! NB: this should be a single-precision real array
- real (kind=4), allocatable, dimension(:,:,:) :: vegfra ! NB: this should be a single-precision real array
+ real (kind=RKIND), allocatable, dimension(:,:) :: rslab, maskslab
+ real (kind=RKIND), allocatable, dimension(:,:) :: maxsnowalb
+ real (kind=RKIND), allocatable, dimension(:,:,:) :: vegfra
+ integer, dimension(:), pointer :: mask_array
+ integer, dimension(grid % nEdges), target :: edge_mask
character (len=1024) :: fname
real (kind=RKIND) :: u, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
- real (kind=4) :: lat, lon, x, y
+ real (kind=RKIND) :: lat, lon, x, y
real (kind=RKIND) :: ptop, p0, phi
real (kind=RKIND) :: lon_Edge
@@ -2127,7 +2137,11 @@
r_earth = a
p0 = 1.e+05
+ interp_list(1) = FOUR_POINT
+ interp_list(2) = SEARCH
+ interp_list(3) = 0
+
!
! Scale all distances and areas from a unit sphere to one with radius a
!
@@ -2281,15 +2295,7 @@
deallocate(rarray)
deallocate(ncat)
- !
- ! Derive LANDMASK
- !
- grid % landmask % array(:) = 0
- do iCell=1, grid % nCells
- if (grid % lu_index % array(iCell) /= 16) grid % landmask % array(iCell) = 1
- end do
-
!
! Interpolate SOILCAT_TOP
!
@@ -2414,7 +2420,45 @@
deallocate(ncat)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! KLUDGE TO FIX SOIL TYPE OVER ANTARCTICA
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ where (grid % lu_index % array == 24) grid % soilcat_top % array = 16
+ where (grid % lu_index % array == 24) grid % soilcat_bot % array = 16
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! CORRECT INCONSISTENT SOIL AND LAND USE DATA
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ do iCell = 1,grid % nCells
+ if (grid % lu_index % array(iCell) == 16 .or. &
+ grid % soilcat_top % array(iCell) == 14 .or. &
+ grid % soilcat_bot % array(iCell) == 14) then
+ if (grid % lu_index % array(iCell) /= 16) then
+ write(0,*) 'Turning lu_index into water at ', iCell
+ grid % lu_index % array(iCell) = 16
+ end if
+ if (grid % soilcat_top % array(iCell) /= 14) then
+ write(0,*) 'Turning soilcat_top into water at ', iCell
+ grid % soilcat_top % array(iCell) = 14
+ end if
+ if (grid % soilcat_bot % array(iCell) /= 14) then
+ write(0,*) 'Turning soilcat_bot into water at ', iCell
+ grid % soilcat_bot % array(iCell) = 14
+ end if
+ end if
+ end do
+
+
!
+ ! Derive LANDMASK
+ !
+ grid % landmask % array(:) = 0
+ do iCell=1, grid % nCells
+ if (grid % lu_index % array(iCell) /= 16) grid % landmask % array(iCell) = 1
+ end do
+
+
+ !
! Interpolate SNOALB
!
nx = 186
@@ -2429,12 +2473,12 @@
grid % snoalb % array(:) = 0.0
call map_set(PROJ_LATLON, proj, &
- latinc = 1.0_4, &
- loninc = 1.0_4, &
- knowni = 1.0_4, &
- knownj = 1.0_4, &
- lat1 = -89.5_4, &
- lon1 = -179.5_4)
+ latinc = 1.0, &
+ loninc = 1.0, &
+ knowni = 1.0, &
+ knownj = 1.0, &
+ lat1 = -89.5, &
+ lon1 = -179.5)
write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'maxsnowalb/',1,'-',180,'.',1,'-',180
write(0,*) trim(fname)
@@ -2456,15 +2500,28 @@
write(0,*) 'min/max = ',minval(rarray(:,:,:)),maxval(rarray(:,:,:))
maxsnowalb(181:360,1:180) = rarray(4:183,4:183,1)
+ interp_list(1) = FOUR_POINT
+ interp_list(2) = W_AVERAGE4
+ interp_list(3) = W_AVERAGE16
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
do iCell=1,grid%nCells
- lat = grid % latCell % array(iCell)*DEG_PER_RAD
- lon = grid % lonCell % array(iCell)*DEG_PER_RAD
- call latlon_to_ij(proj, lat, lon, x, y)
- if (x < 0.0) then
- lon = lon + 360.0
+
+ if (grid % landmask % array(iCell) == 1) then
+ lat = grid % latCell % array(iCell)*DEG_PER_RAD
+ lon = grid % lonCell % array(iCell)*DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
+ if (x < 0.5) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+! grid % snoalb % array(iCell) = interp_sequence(x, y, 1, maxsnowalb, 1, 360, 1, 180, 1, 1, -1.e30, interp_list, 1)
+ grid % snoalb % array(iCell) = interp_sequence(x, y, 1, maxsnowalb, 1, 360, 1, 180, 1, 1, 0., interp_list, 1)
+ else
+ grid % snoalb % array(iCell) = 0.0
end if
- grid % snoalb % array(iCell) = four_pt(360, 180, maxsnowalb, x, y)
+
end do
grid % snoalb % array(:) = grid % snoalb % array(:) / 100.0
@@ -2489,12 +2546,12 @@
grid % greenfrac % array(:,:) = 0.0
call map_set(PROJ_LATLON, proj, &
- latinc = 0.144_4, &
- loninc = 0.144_4, &
- knowni = 1.0_4, &
- knownj = 1.0_4, &
- lat1 = -89.928_4, &
- lon1 = -179.928_4)
+ latinc = 0.144, &
+ loninc = 0.144, &
+ knowni = 1.0, &
+ knownj = 1.0, &
+ lat1 = -89.928, &
+ lon1 = -179.928)
write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'greenfrac/',1,'-',1250,'.',1,'-',1250
write(0,*) trim(fname)
@@ -2517,27 +2574,30 @@
vegfra(1251:2500,1:1250,1:12) = rarray(4:1253,4:1253,1:12)
do iCell=1,grid%nCells
- lat = grid % latCell % array(iCell)*DEG_PER_RAD
- lon = grid % lonCell % array(iCell)*DEG_PER_RAD
- call latlon_to_ij(proj, lat, lon, x, y)
- if (x < 0.0) then
- lon = lon + 360.0
+ if (grid % landmask % array(iCell) == 1) then
+ lat = grid % latCell % array(iCell)*DEG_PER_RAD
+ lon = grid % lonCell % array(iCell)*DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
+ if (x < 0.5) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ do k=1,12
+ grid % greenfrac % array(k,iCell) = interp_sequence(x, y, k, vegfra, 1, 2500, 1, 1250, 1, 12, -1.e30, interp_list, 1)
+ end do
+ else
+ grid % greenfrac % array(:,iCell) = 0.0
end if
- do k=1,12
- grid % greenfrac % array(k,iCell) = four_pt(2500, 1250, vegfra(:,:,k), x, y)
- end do
grid % shdmin % array(iCell) = minval(grid % greenfrac % array(:,iCell))
grid % shdmax % array(iCell) = maxval(grid % greenfrac % array(:,iCell))
+
end do
- grid % greenfrac % array(:,:) = grid % greenfrac % array(:,:) / 100.0
-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!TO-DO: The vegfra field should be interpolated in time from the greenfrac field
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- grid % vegfra % array(:) = grid % greenfrac % array(3,:)
+ grid % vegfra % array(:) = grid % greenfrac % array(1,:)
deallocate(rarray)
deallocate(vegfra)
@@ -2558,12 +2618,12 @@
grid % albedo12m % array(:,:) = 0.0
call map_set(PROJ_LATLON, proj, &
- latinc = 0.144_4, &
- loninc = 0.144_4, &
- knowni = 1.0_4, &
- knownj = 1.0_4, &
- lat1 = -89.928_4, &
- lon1 = -179.928_4)
+ latinc = 0.144, &
+ loninc = 0.144, &
+ knowni = 1.0, &
+ knownj = 1.0, &
+ lat1 = -89.928, &
+ lon1 = -179.928)
write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'albedo_ncep/',1,'-',1250,'.',1,'-',1250
write(0,*) trim(fname)
@@ -2586,22 +2646,25 @@
vegfra(1251:2500,1:1250,1:12) = rarray(4:1253,4:1253,1:12)
do iCell=1,grid%nCells
- lat = grid % latCell % array(iCell)*DEG_PER_RAD
- lon = grid % lonCell % array(iCell)*DEG_PER_RAD
- call latlon_to_ij(proj, lat, lon, x, y)
- if (x < 0.0) then
- lon = lon + 360.0
+ if (grid % landmask % array(iCell) == 1) then
+ lat = grid % latCell % array(iCell)*DEG_PER_RAD
+ lon = grid % lonCell % array(iCell)*DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
+ if (x < 0.5) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ do k=1,12
+ grid % albedo12m % array(k,iCell) = interp_sequence(x, y, k, vegfra, 1, 2500, 1, 1250, 1, 12, 0.0, interp_list, 1)
+ end do
+ else
+ grid % albedo12m % array(:,iCell) = 8.0
end if
- do k=1,12
- grid % albedo12m % array(k,iCell) = four_pt(2500, 1250, vegfra(:,:,k), x, y)
- end do
end do
- grid % albedo12m % array(:,:) = grid % albedo12m % array(:,:) / 100.0
-
deallocate(rarray)
deallocate(vegfra)
+
end if ! config_static_interp
@@ -2628,12 +2691,12 @@
if (field % iproj == PROJ_LATLON) then
call map_set(PROJ_LATLON, proj, &
- latinc = field % deltalat, &
- loninc = field % deltalon, &
- knowni = 1.0_4, &
- knownj = 1.0_4, &
- lat1 = field % startlat, &
- lon1 = field % startlon)
+ latinc = real(field % deltalat), &
+ loninc = real(field % deltalon), &
+ knowni = 1.0, &
+ knownj = 1.0, &
+ lat1 = real(field % startlat), &
+ lon1 = real(field % startlon))
end if
@@ -2649,14 +2712,14 @@
lat = latPoints(i)*DEG_PER_RAD
lon = lonPoints(i)*DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- if (x < 0.0) then
+ if (x < 0.5) then
lon = lon + 360.0
call latlon_to_ij(proj, lat, lon, x, y)
end if
if (ndims == 1) then
- destField1d(i) = four_pt(field % nx, field % ny, field % slab, x, y)
+ destField1d(i) = interp_sequence(x, y, 1, field % slab, 1, field % nx, 1, field % ny, 1, 1, -1.e30, interp_list, 1)
else if (ndims == 2) then
- destField2d(k,i) = four_pt(field % nx, field % ny, field % slab, x, y)
+ destField2d(k,i) = interp_sequence(x, y, 1, field % slab, 1, field % nx, 1, field % ny, 1, 1, -1.e30, interp_list, 1)
end if
end do
end if
@@ -2683,8 +2746,8 @@
! Fourth order smoother for terrain
-! nsmterrain = 2
- nsmterrain = 0 ! NO SMOOTHING WHEN USING GFS TERRAIN HEIGHT DIRECTLY
+ nsmterrain = 2
+! nsmterrain = 0 ! NO SMOOTHING WHEN USING GFS TERRAIN HEIGHT DIRECTLY
do i=1,nsmterrain
@@ -2965,6 +3028,42 @@
if (config_met_interp) then
!
+ ! First, try to locate the LANDSEA field for use as an interpolation mask
+ !
+ call read_met_init(trim(config_met_prefix), .false., trim(config_init_date), istatus)
+
+ if (istatus /= 0) then
+ write(0,*) 'Error reading initial met data'
+ return
+ end if
+
+ call read_next_met_field(field, istatus)
+
+ do while (istatus == 0)
+ if (index(field % field, 'LANDSEA') /= 0) then
+
+ allocate(maskslab(-3:field % nx+3, field % ny))
+ maskslab(1:field % nx, 1:field % ny) = field % slab(1:field % nx, 1:field % ny)
+ maskslab(0, 1:field % ny) = field % slab(field % nx, 1:field % ny)
+ maskslab(-1, 1:field % ny) = field % slab(field % nx-1, 1:field % ny)
+ maskslab(-2, 1:field % ny) = field % slab(field % nx-2, 1:field % ny)
+ maskslab(field % nx+1, 1:field % ny) = field % slab(1, 1:field % ny)
+ maskslab(field % nx+2, 1:field % ny) = field % slab(2, 1:field % ny)
+ maskslab(field % nx+3, 1:field % ny) = field % slab(3, 1:field % ny)
+write(0,*) 'minval, maxval of LANDSEA = ', minval(maskslab), maxval(maskslab)
+
+ end if
+
+ deallocate(field % slab)
+ call read_next_met_field(field, istatus)
+ end do
+
+ call read_met_close()
+
+ edge_mask(:) = 1
+
+
+ !
! Horizontally interpolate meteorological data
!
allocate(vert_level(config_nfglevels))
@@ -2978,7 +3077,20 @@
end if
call read_next_met_field(field, istatus)
+
do while (istatus == 0)
+
+ interp_list(1) = FOUR_POINT
+ interp_list(2) = SEARCH
+ interp_list(3) = 0
+
+ maskval = -1.0
+ masked = -1
+ fillval = 0.0
+ msgval = -1.e30
+
+ mask_array => grid % landmask % array
+
if (index(field % field, 'UU') /= 0 .or. &
index(field % field, 'VV') /= 0 .or. &
index(field % field, 'TT') /= 0 .or. &
@@ -3031,12 +3143,12 @@
if (field % iproj == PROJ_LATLON) then
call map_set(PROJ_LATLON, proj, &
- latinc = field % deltalat, &
- loninc = field % deltalon, &
- knowni = 1.0_4, &
- knownj = 1.0_4, &
- lat1 = field % startlat, &
- lon1 = field % startlon)
+ latinc = real(field % deltalat), &
+ loninc = real(field % deltalon), &
+ knowni = 1.0, &
+ knownj = 1.0, &
+ lat1 = real(field % startlat), &
+ lon1 = real(field % startlon))
end if
@@ -3045,6 +3157,9 @@
!
if (index(field % field, 'UU') /= 0) then
write(0,*) 'Interpolating U at ', k, vert_level(k)
+
+ mask_array => edge_mask
+
nInterpPoints = grid % nEdges
latPoints => grid % latEdge % array
lonPoints => grid % lonEdge % array
@@ -3052,6 +3167,9 @@
ndims = 2
else if (index(field % field, 'VV') /= 0) then
write(0,*) 'Interpolating V at ', k, vert_level(k)
+
+ mask_array => edge_mask
+
nInterpPoints = grid % nEdges
latPoints => grid % latEdge % array
lonPoints => grid % lonEdge % array
@@ -3108,6 +3226,17 @@
ndims = 1
else if (index(field % field, 'SM000010') /= 0) then
write(0,*) 'Interpolating SM000010'
+
+ interp_list(1) = SIXTEEN_POINT
+ interp_list(2) = FOUR_POINT
+ interp_list(3) = W_AVERAGE4
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ maskval = 0.0
+ masked = 0
+ fillval = 1.0
+
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
@@ -3116,6 +3245,17 @@
ndims = 2
else if (index(field % field, 'SM010040') /= 0) then
write(0,*) 'Interpolating SM010040'
+
+ interp_list(1) = SIXTEEN_POINT
+ interp_list(2) = FOUR_POINT
+ interp_list(3) = W_AVERAGE4
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ maskval = 0.0
+ masked = 0
+ fillval = 1.0
+
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
@@ -3124,6 +3264,17 @@
ndims = 2
else if (index(field % field, 'SM040100') /= 0) then
write(0,*) 'Interpolating SM040100'
+
+ interp_list(1) = SIXTEEN_POINT
+ interp_list(2) = FOUR_POINT
+ interp_list(3) = W_AVERAGE4
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ maskval = 0.0
+ masked = 0
+ fillval = 1.0
+
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
@@ -3132,6 +3283,17 @@
ndims = 2
else if (index(field % field, 'SM100200') /= 0) then
write(0,*) 'Interpolating SM100200'
+
+ interp_list(1) = SIXTEEN_POINT
+ interp_list(2) = FOUR_POINT
+ interp_list(3) = W_AVERAGE4
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ maskval = 0.0
+ masked = 0
+ fillval = 1.0
+
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
@@ -3140,6 +3302,17 @@
ndims = 2
else if (index(field % field, 'ST000010') /= 0) then
write(0,*) 'Interpolating ST000010'
+
+ interp_list(1) = SIXTEEN_POINT
+ interp_list(2) = FOUR_POINT
+ interp_list(3) = W_AVERAGE4
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ maskval = 0.0
+ masked = 0
+ fillval = 285.0
+
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
@@ -3148,6 +3321,17 @@
ndims = 2
else if (index(field % field, 'ST010040') /= 0) then
write(0,*) 'Interpolating ST010040'
+
+ interp_list(1) = SIXTEEN_POINT
+ interp_list(2) = FOUR_POINT
+ interp_list(3) = W_AVERAGE4
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ maskval = 0.0
+ masked = 0
+ fillval = 285.0
+
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
@@ -3156,6 +3340,17 @@
ndims = 2
else if (index(field % field, 'ST040100') /= 0) then
write(0,*) 'Interpolating ST040100'
+
+ interp_list(1) = SIXTEEN_POINT
+ interp_list(2) = FOUR_POINT
+ interp_list(3) = W_AVERAGE4
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ maskval = 0.0
+ masked = 0
+ fillval = 285.0
+
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
@@ -3164,6 +3359,17 @@
ndims = 2
else if (index(field % field, 'ST100200') /= 0) then
write(0,*) 'Interpolating ST100200'
+
+ interp_list(1) = SIXTEEN_POINT
+ interp_list(2) = FOUR_POINT
+ interp_list(3) = W_AVERAGE4
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ maskval = 0.0
+ masked = 0
+ fillval = 285.0
+
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
@@ -3172,6 +3378,14 @@
ndims = 2
else if (index(field % field, 'SNOW') /= 0) then
write(0,*) 'Interpolating SNOW'
+
+ interp_list(1) = FOUR_POINT
+ interp_list(2) = W_AVERAGE4
+ interp_list(3) = 0
+
+ masked = 0
+ fillval = 0.0
+
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
@@ -3179,6 +3393,17 @@
ndims = 1
else if (index(field % field, 'SEAICE') /= 0) then
write(0,*) 'Interpolating SEAICE'
+
+ interp_list(1) = SIXTEEN_POINT
+ interp_list(2) = FOUR_POINT
+ interp_list(3) = W_AVERAGE4
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ maskval = 1.0
+ masked = 1
+ fillval = 0.0
+
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
@@ -3193,23 +3418,39 @@
ndims = 1
end if
+ allocate(rslab(-3:field % nx+3, field % ny))
+ rslab(1:field % nx, 1:field % ny) = field % slab(1:field % nx, 1:field % ny)
+ rslab(0, 1:field % ny) = field % slab(field % nx, 1:field % ny)
+ rslab(-1, 1:field % ny) = field % slab(field % nx-1, 1:field % ny)
+ rslab(-2, 1:field % ny) = field % slab(field % nx-2, 1:field % ny)
+ rslab(field % nx+1, 1:field % ny) = field % slab(1, 1:field % ny)
+ rslab(field % nx+2, 1:field % ny) = field % slab(2, 1:field % ny)
+ rslab(field % nx+3, 1:field % ny) = field % slab(3, 1:field % ny)
+
do i=1,nInterpPoints
- lat = latPoints(i)*DEG_PER_RAD
- lon = lonPoints(i)*DEG_PER_RAD
- call latlon_to_ij(proj, lat, lon, x, y)
- if (x < 0.0) then
- lon = lon + 360.0
+ if (mask_array(i) /= masked) then
+ lat = latPoints(i)*DEG_PER_RAD
+ lon = lonPoints(i)*DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
+ if (x < 0.5) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ if (ndims == 1) then
+ destField1d(i) = interp_sequence(x, y, 1, rslab, -3, field % nx + 3, 1, field % ny, 1, 1, msgval, interp_list, 1, maskval=maskval, mask_array=maskslab)
+ else if (ndims == 2) then
+ destField2d(k,i) = interp_sequence(x, y, 1, rslab, -3, field % nx + 3, 1, field % ny, 1, 1, msgval, interp_list, 1, maskval=maskval, mask_array=maskslab)
+ end if
+ else
+ if (ndims == 1) then
+ destField1d(i) = fillval
+ else if (ndims == 2) then
+ destField2d(k,i) = fillval
+ end if
end if
- if (ndims == 1) then
- destField1d(i) = four_pt(field % nx, field % ny, field % slab, x, y)
- else if (ndims == 2) then
- destField2d(k,i) = four_pt(field % nx, field % ny, field % slab, x, y)
- end if
end do
-if (index(field % field, 'SEAICE') /= 0) then
- write(0,*) minval(destfield1d(:)), maxval(destfield1d(:))
-end if
+
+ deallocate(rslab)
end if
@@ -3241,7 +3482,140 @@
fg % snowc % array(:) = 0.0
where (fg % snow % array > 0.0) fg % snowc % array = 1.0
+ ! Set TMN based on TSLB
+ do iCell=1,grid%nCells
+ fg % tmn % array(iCell) = fg % tslb % array(4,iCell)
+ end do
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!MGD CHECK
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+do iCell=1,grid%nCells
+ if (grid % landmask % array(iCell) == 1) then
+ if (fg % tslb % array(1,iCell) <= 0.0) write(0,*) 'Bad tslb ', 1, iCell
+ if (fg % tslb % array(2,iCell) <= 0.0) write(0,*) 'Bad tslb ', 2, iCell
+ if (fg % tslb % array(3,iCell) <= 0.0) write(0,*) 'Bad tslb ', 3, iCell
+ if (fg % tslb % array(4,iCell) <= 0.0) write(0,*) 'Bad tslb ', 4, iCell
+
+ if (fg % smois % array(1,iCell) <= 0.0) write(0,*) 'Bad smois ', 1, iCell
+ if (fg % smois % array(2,iCell) <= 0.0) write(0,*) 'Bad smois ', 2, iCell
+ if (fg % smois % array(3,iCell) <= 0.0) write(0,*) 'Bad smois ', 3, iCell
+ if (fg % smois % array(4,iCell) <= 0.0) write(0,*) 'Bad smois ', 4, iCell
+ end if
+end do
+write(0,*) 'Done with soil consistency check'
+
+
+ !
+ ! Get SEAICE from a separate file
+ !
+ call read_met_init('SEAICE_FRACTIONAL', .true., trim(config_init_date), istatus)
+
+ if (istatus /= 0) then
+ write(0,*) 'Error reading initial met data'
+ return
+ end if
+
+ call read_next_met_field(field, istatus)
+ do while (istatus == 0)
+ if (index(field % field, 'SEAICE') /= 0) then
+
+write(0,*) 'PROCESSING SEAICE'
+
+ !
+ ! Set up projection
+ !
+ call map_init(proj)
+
+ if (field % iproj == PROJ_PS) then
+ call map_set(PROJ_PS, proj, &
+ dx = real(field % dx * 1000.0), &
+ truelat1 = real(field % truelat1), &
+ stdlon = real(field % xlonc), &
+ knowni = real(field % nx / 2.0), &
+ knownj = real(field % ny / 2.0), &
+ lat1 = real(field % startlat), &
+ lon1 = real(field % startlon))
+ end if
+
+ if (index(field % field, 'SEAICE') /= 0) then
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % xice % array
+ ndims = 1
+ end if
+
+ interp_list(1) = FOUR_POINT
+ interp_list(2) = W_AVERAGE4
+ interp_list(3) = 0
+
+ masked = 1
+ fillval = 0.0
+ msgval = 1.01
+ mask_array => grid % landmask % array
+
+
+ allocate(rslab(field % nx, field % ny))
+ rslab(:,:) = field % slab(:,:)
+ do i=1,nInterpPoints
+ if (mask_array(i) /= masked) then
+ lat = latPoints(i)*DEG_PER_RAD
+ lon = lonPoints(i)*DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ if (x < 0.5) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ if (ndims == 1) then
+ destField1d(i) = interp_sequence(x, y, 1, rslab, 1, field % nx, 1, field % ny, 1, 1, msgval, interp_list, 1)
+ if (destField1d(i) == msgval) destField1d(i) = fillval
+ else if (ndims == 2) then
+ destField2d(k,i) = interp_sequence(x, y, 1, rslab, 1, field % nx, 1, field % ny, 1, 1, msgval, interp_list, 1)
+ if (destField2d(k,i) == msgval) destField2d(k,i) = fillval
+ end if
+ else
+ if (ndims == 1) then
+ destField1d(i) = fillval
+ else if (ndims == 2) then
+ destField2d(k,i) = fillval
+ end if
+ end if
+ end do
+ deallocate(rslab)
+
+ end if
+
+ deallocate(field % slab)
+ call read_next_met_field(field, istatus)
+ end do
+
+ call read_met_close()
+
+ if (allocated(maskslab)) deallocate(maskslab)
+
+ ! Freeze really cold ocean
+ where (fg % sst % array < 271.0 .and. grid % landmask % array == 0) fg % xice % array = 1.0
+
+ ! Set SEAICE (0/1 flag) based on XICE (fractional ice coverage)
+ fg % seaice % array(:) = 0.0
+ where (fg % xice % array >= 0.5) fg % seaice % array = 1.0
+
+
+ !
+ ! For now, hard-wire soil layer depths and thicknesses
+ !
+ fg % dzs % array(1,:) = 0.10
+ fg % dzs % array(2,:) = 0.30
+ fg % dzs % array(3,:) = 0.60
+ fg % dzs % array(4,:) = 1.00
+
+ fg % dz % array(1,:) = 0.05
+ fg % dz % array(2,:) = 0.25
+ fg % dz % array(3,:) = 0.70
+ fg % dz % array(4,:) = 1.50
+
+
!
! Compute normal wind component and store in fg%u
!
@@ -3519,10 +3893,19 @@
end if ! config_met_interp
+
+ ! Calculate surface pressure
+ do iCell=1,grid%nCells
+ diag % surface_pressure % array(iCell) = 0.5*gravity/rdzw(1) &
+ * (1.25* rho(1,iCell) * (1. + scalars(state % index_qv, 1, iCell)) &
+ - 0.25* rho(2,iCell) * (1. + scalars(state % index_qv, 2, iCell)))
+ diag % surface_pressure % array(iCell) = diag % surface_pressure % array(iCell) + pp(1,iCell) + ppb(1,iCell)
+ end do
end subroutine nhyd_test_case_gfs
+#if 0
real function four_pt(nx, ny, array, xx, yy)
implicit none
@@ -3573,6 +3956,7 @@
return
end function four_pt
+#endif
real function sphere_distance(lat1, lon1, lat2, lon2, radius)
</font>
</pre>