<p><b>dwj07@fsu.edu</b> 2011-12-22 10:18:00 -0700 (Thu, 22 Dec 2011)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Merging trunk to time_averaging branch<br>
</p><hr noshade><pre><font color="gray">
Property changes on: branches/ocean_projects/time_averaging
___________________________________________________________________
Modified: svn:mergeinfo
- /branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
+ /branches/cam_mpas_nh:1260-1270
/branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
/trunk/mpas:1271-1274
Modified: branches/ocean_projects/time_averaging/Makefile
===================================================================
--- branches/ocean_projects/time_averaging/Makefile 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/Makefile 2011-12-22 17:18:00 UTC (rev 1275)
@@ -31,7 +31,7 @@
"CFLAGS = -g" \
"LDFLAGS = -g -C" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
xlf:
( make all \
@@ -43,7 +43,7 @@
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
ftn:
( make all \
@@ -55,7 +55,7 @@
"CFLAGS = -fast" \
"LDFLAGS = " \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
pgi:
( make all \
@@ -67,7 +67,7 @@
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
pgi-nersc:
( make all \
@@ -79,7 +79,7 @@
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
pgi-llnl:
( make all \
@@ -91,7 +91,7 @@
"CFLAGS = -fast" \
"LDFLAGS = " \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
pgi-serial:
( make all \
@@ -103,7 +103,7 @@
"CFLAGS = -O0 -g" \
"LDFLAGS = -O0 -g -Mbounds -Mchkptr" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
ifort-serial:
( make all \
@@ -115,7 +115,7 @@
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
ifort-papi:
( make all \
@@ -127,7 +127,7 @@
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_PAPI -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" \
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_PAPI -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" \
"PAPILIBS = -L$(PAPI)/lib -lpapi" )
ifort-papi-serial:
@@ -140,7 +140,7 @@
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_PAPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" \
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_PAPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" \
"PAPILIBS = -L$(PAPI)/lib -lpapi" )
ifort:
@@ -153,7 +153,7 @@
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
gfortran:
( make all \
@@ -165,7 +165,7 @@
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3 -m64" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
gfortran-serial:
( make all \
@@ -177,7 +177,7 @@
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3 -m64" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
g95:
( make all \
@@ -189,7 +189,7 @@
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
g95-serial:
( make all \
@@ -201,7 +201,7 @@
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
pathscale-nersc:
( make all \
@@ -213,7 +213,7 @@
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
cray-nersc:
( make all \
@@ -225,7 +225,7 @@
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
intel-nersc:
( make all \
@@ -237,7 +237,7 @@
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
CPPINCLUDES = -I../inc -I$(NETCDF)/include -I$(PAPI)/include
FCINCLUDES = -I../inc -I$(NETCDF)/include -I$(PAPI)/include
Modified: branches/ocean_projects/time_averaging/src/core_atmos_physics/mpas_atmphys_date_time.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_atmos_physics/mpas_atmphys_date_time.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_atmos_physics/mpas_atmphys_date_time.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,6 +1,8 @@
!=============================================================================================
module mpas_atmphys_date_time
+ use mpas_kind_types
+
implicit none
private
public:: get_julgmt, &
Modified: branches/ocean_projects/time_averaging/src/core_atmos_physics/mpas_atmphys_initialize_real.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_atmos_physics/mpas_atmphys_initialize_real.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_atmos_physics/mpas_atmphys_initialize_real.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,6 @@
!=============================================================================================
module mpas_atmphys_initialize_real
+ use mpas_kind_types
use mpas_configure, only: config_met_prefix, &
config_frac_seaice, &
config_input_sst, &
@@ -82,12 +83,12 @@
if(field % iproj == PROJ_LATLON) then
call map_set(PROJ_LATLON, proj, &
- latinc = real(field % deltalat), &
- loninc = real(field % deltalon), &
- knowni = 1.0, &
- knownj = 1.0, &
- lat1 = real(field % startlat), &
- lon1 = real(field % startlon))
+ latinc = real(field % deltalat,RKIND), &
+ loninc = real(field % deltalon,RKIND), &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
end if
!Interpolate field to each MPAS grid cell:
@@ -110,10 +111,10 @@
if(index(field % field,'SST') /= 0) then
fg % sst % array(iCell) = interp_sequence(x,y,1,slab_r8,1,field%nx, &
- 1,field%ny,1,1,-1.e30,interp_list,1)
+ 1,field%ny,1,1,-1.e30_RKIND,interp_list,1)
elseif(index(field % field,'SEAICE') /= 0) then
fg % xice % array(iCell) = interp_sequence(x,y,1,slab_r8,1,field%nx, &
- 1,field%ny,1,1,-1.e30,interp_list,1)
+ 1,field%ny,1,1,-1.e30_RKIND,interp_list,1)
endif
end do
Modified: branches/ocean_projects/time_averaging/src/core_atmos_physics/mpas_atmphys_vars.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_atmos_physics/mpas_atmphys_vars.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_atmos_physics/mpas_atmphys_vars.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,7 @@
!=============================================================================================
module mpas_atmphys_vars
+
+ use mpas_kind_types
implicit none
public
Modified: branches/ocean_projects/time_averaging/src/core_hyd_atmos/mpas_atmh_advection.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_hyd_atmos/mpas_atmh_advection.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_hyd_atmos/mpas_atmh_advection.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,6 @@
module atmh_advection
+ use mpas_kind_types
use mpas_grid_types
use mpas_configure
use mpas_constants
@@ -116,7 +117,7 @@
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
xc(2), yc(2), zc(2), &
- 0., 0., 1. )
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND )
! angles from cell center to neighbor centers (thetav)
@@ -335,7 +336,7 @@
! Computes the angle between arcs AB and AC, given points A, B, and C
! Equation numbers w.r.t. http://mathworld.wolfram.com/SphericalTrigonometry.html
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
+ real (kind=RKIND) function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
implicit none
@@ -355,9 +356,9 @@
real (kind=RKIND) :: s ! Semiperimeter of the triangle
real (kind=RKIND) :: sin_angle
- a = acos(max(min(bx*cx + by*cy + bz*cz,1.0),-1.0)) ! Eqn. (3)
- b = acos(max(min(ax*cx + ay*cy + az*cz,1.0),-1.0)) ! Eqn. (2)
- c = acos(max(min(ax*bx + ay*by + az*bz,1.0),-1.0)) ! Eqn. (1)
+ a = acos(max(min(bx*cx + by*cy + bz*cz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (3)
+ b = acos(max(min(ax*cx + ay*cy + az*cz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (2)
+ c = acos(max(min(ax*bx + ay*by + az*bz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (1)
ABx = bx - ax
ABy = by - ay
@@ -373,12 +374,12 @@
s = 0.5*(a + b + c)
! sin_angle = sqrt((sin(s-b)*sin(s-c))/(sin(b)*sin(c))) ! Eqn. (28)
- sin_angle = sqrt(min(1.,max(0.,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
+ sin_angle = sqrt(min(1.0_RKIND,max(0.0_RKIND,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
if ((Dx*ax + Dy*ay + Dz*az) >= 0.0) then
- sphere_angle = 2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ sphere_angle = 2.0 * asin(max(min(sin_angle,1.0_RKIND),-1.0_RKIND))
else
- sphere_angle = -2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ sphere_angle = -2.0 * asin(max(min(sin_angle,1.0_RKIND),-1.0_RKIND))
end if
end function sphere_angle
@@ -390,7 +391,7 @@
! Computes the angle between vectors AB and AC, given points A, B, and C, and
! a vector (u,v,w) normal to the plane.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
+ real (kind=RKIND) function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
implicit none
@@ -425,9 +426,9 @@
cos_angle = (ABx*ACx + ABy*ACy + ABz*ACz) / (mAB * mAC)
if ((Dx*u + Dy*v + Dz*w) >= 0.0) then
- plane_angle = acos(max(min(cos_angle,1.0),-1.0))
+ plane_angle = acos(max(min(cos_angle,1.0_RKIND),-1.0_RKIND))
else
- plane_angle = -acos(max(min(cos_angle,1.0),-1.0))
+ plane_angle = -acos(max(min(cos_angle,1.0_RKIND),-1.0_RKIND))
end if
end function plane_angle
@@ -440,7 +441,7 @@
! B=(bx, by, bz). It is assumed that both A and B lie on the surface of the
! same sphere centered at the origin.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function arc_length(ax, ay, az, bx, by, bz)
+ real (kind=RKIND) function arc_length(ax, ay, az, bx, by, bz)
implicit none
@@ -575,7 +576,7 @@
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
-SUBROUTine MIGS (A,N,X,INDX)
+SUBROUTINE MIGS (A,N,X,INDX)
!
! Subroutine to invert matrix A(N,N) with the inverse stored
! in X(N,N) in the output. Copyright (c) Tao Pang 2001.
@@ -617,7 +618,7 @@
X(J,I) = X(J,I)/A(INDX(J),J)
END DO
END DO
-END SUBROUTine MIGS
+END SUBROUTINE MIGS
SUBROUTINE ELGS (A,N,INDX)
@@ -646,7 +647,7 @@
DO I = 1, N
C1= 0.0
DO J = 1, N
- C1 = AMAX1(C1,ABS(A(I,J)))
+ C1 = MAX(C1,ABS(A(I,J)))
END DO
C(I) = C1
END DO
Modified: branches/ocean_projects/time_averaging/src/core_hyd_atmos/mpas_atmh_test_cases.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_hyd_atmos/mpas_atmh_test_cases.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_hyd_atmos/mpas_atmh_test_cases.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -258,7 +258,7 @@
if (config_test_case == 2) then
r_pert = sphere_distance( grid % latEdge % array (iEdge), grid % lonEdge % array (iEdge), &
- lat_pert, lon_pert, 1.)/(pert_radius)
+ lat_pert, lon_pert, 1.0_RKIND)/(pert_radius)
u_pert = u_perturbation*exp(-r_pert**2)*(lat2-lat1)*a/grid % dvEdge % array(iEdge)
else if (config_test_case == 3) then
@@ -468,7 +468,7 @@
end subroutine atmh_test_case_1
- real function sphere_distance(lat1, lon1, lat2, lon2, radius)
+ real (kind=RKIND) function sphere_distance(lat1, lon1, lat2, lon2, radius)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
! sphere with given radius.
@@ -487,7 +487,7 @@
end function sphere_distance
- real function AA(theta)
+ real (kind=RKIND) function AA(theta)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! A, used in height field computation for Rossby-Haurwitz wave
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -506,7 +506,7 @@
end function AA
- real function BB(theta)
+ real (kind=RKIND) function BB(theta)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! B, used in height field computation for Rossby-Haurwitz wave
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -524,7 +524,7 @@
end function BB
- real function CC(theta)
+ real (kind=RKIND) function CC(theta)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! C, used in height field computation for Rossby-Haurwitz wave
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Modified: branches/ocean_projects/time_averaging/src/core_hyd_atmos/mpas_atmh_time_integration.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_hyd_atmos/mpas_atmh_time_integration.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_hyd_atmos/mpas_atmh_time_integration.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1652,9 +1652,9 @@
! add in vertical flux to get max and min estimate
s_max_update(iScalar) = s_max_update(iScalar) &
- - rdnw(k) * (max(0.,v_flux(iScalar,iCell,km0)) - min(0.,v_flux(iScalar,iCell,km1)))
+ - rdnw(k) * (max(0.0_RKIND,v_flux(iScalar,iCell,km0)) - min(0.0_RKIND,v_flux(iScalar,iCell,km1)))
s_min_update(iScalar) = s_min_update(iScalar) &
- - rdnw(k) * (min(0.,v_flux(iScalar,iCell,km0)) - max(0.,v_flux(iScalar,iCell,km1)))
+ - rdnw(k) * (min(0.0_RKIND,v_flux(iScalar,iCell,km0)) - max(0.0_RKIND,v_flux(iScalar,iCell,km1)))
end do
@@ -1671,8 +1671,8 @@
fdir = -1.0
end if
flux = -fdir * h_flux(iScalar,iEdge)/grid % areaCell % array(iCell)
- s_max_update(iScalar) = s_max_update(iScalar) + max(0.,flux)
- s_min_update(iScalar) = s_min_update(iScalar) + min(0.,flux)
+ s_max_update(iScalar) = s_max_update(iScalar) + max(0.0_RKIND,flux)
+ s_min_update(iScalar) = s_min_update(iScalar) + min(0.0_RKIND,flux)
end do
@@ -1687,9 +1687,9 @@
s_min_update (iScalar) = s_min_update (iScalar) / h_new (k,iCell)
s_upwind = s_update(iScalar,iCell,km0) / h_new(k,iCell)
if ( s_max_update(iScalar) > s_max(iScalar) .and. config_monotonic) &
- scale_in (iScalar,iCell,km0) = max(0.,(s_max(iScalar)-s_upwind)/(s_max_update(iScalar)-s_upwind+eps))
+ scale_in (iScalar,iCell,km0) = max(0.0_RKIND,(s_max(iScalar)-s_upwind)/(s_max_update(iScalar)-s_upwind+eps))
if ( s_min_update(iScalar) < s_min(iScalar) ) &
- scale_out (iScalar,iCell,km0) = max(0.,(s_upwind-s_min(iScalar))/(s_upwind-s_min_update(iScalar)+eps))
+ scale_out (iScalar,iCell,km0) = max(0.0_RKIND,(s_upwind-s_min(iScalar))/(s_upwind-s_min_update(iScalar)+eps))
end do
end do ! end loop over cells to compute scale factor
Modified: branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_hinterp.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_hinterp.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_hinterp.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -545,10 +545,10 @@
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))
+ fxfy = max(0.0_RKIND, 1.0 - sqrt((xx-real(ifx))**2+(yy-real(ify))**2))
+ fxcy = max(0.0_RKIND, 1.0 - sqrt((xx-real(ifx))**2+(yy-real(icy))**2))
+ cxfy = max(0.0_RKIND, 1.0 - sqrt((xx-real(icx))**2+(yy-real(ify))**2))
+ cxcy = max(0.0_RKIND, 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. &
@@ -744,13 +744,13 @@
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))
+ weights(i,j) = max(0.0_RKIND, 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))
+ weights(i,j) = max(0.0_RKIND, 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
end if
end if
Modified: branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_llxy.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_llxy.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_llxy.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -130,6 +130,8 @@
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ USE MPAS_KIND_TYPES
+
INTEGER, PARAMETER :: HH=4, VV=5
REAL (KIND=RKIND), PARAMETER :: PI = 3.141592653589793
@@ -1154,10 +1156,10 @@
! intersects the Earth's surface at each of the distinctly different
! latitudes
IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
- cone = ALOG10(COS(truelat1*rad_per_deg)) - &
- ALOG10(COS(truelat2*rad_per_deg))
- cone = cone /(ALOG10(TAN((45.0 - ABS(truelat1)/2.0) * rad_per_deg)) - &
- ALOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg)))
+ cone = LOG10(COS(truelat1*rad_per_deg)) - &
+ LOG10(COS(truelat2*rad_per_deg))
+ cone = cone /(LOG10(TAN((45.0 - ABS(truelat1)/2.0) * rad_per_deg)) - &
+ LOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg)))
ELSE
cone = SIN(ABS(truelat1)*rad_per_deg )
ENDIF
@@ -1220,9 +1222,9 @@
! Longitude
lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone
# if ( defined (G95) && ( DA_CORE == 1 ) )
- lon = DMOD(lon+360., 360.)
+ lon = DMOD(lon+360., 360.0_RKIND)
# else
- lon = AMOD(lon+360., 360.)
+ lon = MOD(lon+360., 360.0_RKIND)
# endif
! Latitude. Latitude determined by solving an equation adapted
@@ -1323,7 +1325,7 @@
proj%rsw = 0.
IF (proj%lat1 .NE. 0.) THEN
- proj%rsw = (ALOG(TAN(0.5*((proj%lat1+90.)*rad_per_deg))))/proj%dlon
+ proj%rsw = (LOG(TAN(0.5*((proj%lat1+90.)*rad_per_deg))))/proj%dlon
ENDIF
RETURN
@@ -1347,7 +1349,7 @@
IF (deltalon .LT. -180.) deltalon = deltalon + 360.
IF (deltalon .GT. 180.) deltalon = deltalon - 360.
i = proj%knowni + (deltalon/(proj%dlon*deg_per_rad))
- j = proj%knownj + (ALOG(TAN(0.5*((lat + 90.) * rad_per_deg)))) / &
+ j = proj%knownj + (LOG(TAN(0.5*((lat + 90.) * rad_per_deg)))) / &
proj%dlon - proj%rsw
RETURN
@@ -1531,7 +1533,7 @@
! Try to determine whether this domain has global coverage
if (abs(proj%lat1 - proj%latinc/2. + 90.) < 0.001 .and. &
- abs(mod(proj%lon1 - proj%loninc/2. - proj%stdlon,360.)) < 0.001) then
+ abs(mod(proj%lon1 - proj%loninc/2. - proj%stdlon,360.0_RKIND)) < 0.001) then
global_domain = .true.
else
global_domain = .false.
Modified: branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_queue.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_queue.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_queue.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -7,6 +7,9 @@
module init_atm_queue
+ use mpas_kind_types
+
+
type q_data ! The user-defined datatype to store in the queue
real (kind=RKIND) :: lat, lon
integer :: x, y
Modified: branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -166,8 +166,8 @@
integer :: eoe, j
integer, dimension(:), pointer :: nEdgesOnEdge
integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge
- real, dimension(:), pointer :: dvEdge, AreaCell
- real, dimension(:,:), pointer :: weightsOnEdge
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, AreaCell
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
real (kind=RKIND) :: u, v, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
@@ -626,7 +626,7 @@
if (config_test_case == 2) then
r_pert = sphere_distance( grid % latEdge % array (iEdge), grid % lonEdge % array (iEdge), &
- lat_pert, lon_pert, 1.)/(pert_radius)
+ lat_pert, lon_pert, 1.0_RKIND)/(pert_radius)
u_pert = u_perturbation*exp(-r_pert**2)*(lat2-lat1)*a/grid % dvEdge % array(iEdge)
else if (config_test_case == 3) then
@@ -757,9 +757,9 @@
if (config_theta_adv_order ==3) then
diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
- - sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ - sign(1.0_RKIND,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
- + sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ + sign(1.0_RKIND,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
end if
end do
@@ -855,7 +855,7 @@
! renormalize for setting cell-face fluxes
do k=1,nz1
- flux_zonal(k) = sign(1.,lat2_in-lat1_in)*flux_zonal(k)*dlat*a/dvEdge/u0
+ flux_zonal(k) = sign(1.0_RKIND,lat2_in-lat1_in)*flux_zonal(k)*dlat*a/dvEdge/u0
end do
end subroutine init_atm_calc_flux_zonal
@@ -885,7 +885,7 @@
integer :: eoe, j
integer, dimension(:), pointer :: nEdgesOnEdge
integer, dimension(:,:), pointer :: edgesOnEdge
- real, dimension(:,:), pointer :: weightsOnEdge
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, cell1, cell2, nCellsSolve
integer :: index_qv
@@ -1200,7 +1200,7 @@
temp = p(k,i)*thi(k,i)
pres = p0*p(k,i)**(1./rcp)
qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
- scalars(index_qv,k,i) = amin1(0.014,rh(k,i)*qvs)
+ scalars(index_qv,k,i) = min(0.014_RKIND,rh(k,i)*qvs)
end do
end do
@@ -1422,8 +1422,8 @@
integer :: eoe, j
integer, dimension(:), pointer :: nEdgesOnEdge
integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge
- real, dimension(:), pointer :: dvEdge, AreaCell, xCell, yCell
- real, dimension(:,:), pointer :: weightsOnEdge
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, AreaCell, xCell, yCell
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
integer :: index_qv
@@ -1649,7 +1649,7 @@
! smoothing grid for the upper level >> but not propoer for parallel programing
dzmin=.7
do k=2,nz1
- sm = .25*min((zc(k)-zc(k-1))/dz,1.)
+ sm = .25*min((zc(k)-zc(k-1))/dz,1.0_RKIND)
do i=1,grid % nCells
hx(k,i) = hx(k-1,i)
end do
@@ -1826,7 +1826,7 @@
temp = p(k,i)*t(k,i)
pres = p0*p(k,i)**(1./rcp)
qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
- scalars(index_qv,k,i) = amin1(0.014,rh(k,i)*qvs)
+ scalars(index_qv,k,i) = min(0.014_RKIND,rh(k,i)*qvs)
end do
do k=1,nz1
@@ -1949,9 +1949,9 @@
if (config_theta_adv_order ==3) then
diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
- - sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ - sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
- + sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ + sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
end if
end do
@@ -2520,12 +2520,12 @@
grid % soiltemp % array(:) = 0.0
call map_set(PROJ_LATLON, proj, &
- latinc = 1.0, &
- loninc = 1.0, &
- knowni = 1.0, &
- knownj = 1.0, &
- lat1 = -89.5, &
- lon1 = -179.5)
+ latinc = 1.0_RKIND, &
+ loninc = 1.0_RKIND, &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = -89.5_RKIND, &
+ lon1 = -179.5_RKIND)
write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'soiltemp_1deg/',1,'-',180,'.',1,'-',180
write(0,*) trim(fname)
@@ -2568,8 +2568,8 @@
end if
if (y < 1.0) y = 1.0
if (y > 179.0) y = 179.0
-! grid % soiltemp % array(iCell) = interp_sequence(x, y, 1, soiltemp_1deg, 1, 360, 1, 180, 1, 1, -1.e30, interp_list, 1)
- grid % soiltemp % array(iCell) = interp_sequence(x, y, 1, soiltemp_1deg, 1, 360, 1, 180, 1, 1, 0., interp_list, 1)
+! grid % soiltemp % array(iCell) = interp_sequence(x, y, 1, soiltemp_1deg, 1, 360, 1, 180, 1, 1, -1.e30_RKIND, interp_list, 1)
+ grid % soiltemp % array(iCell) = interp_sequence(x, y, 1, soiltemp_1deg, 1, 360, 1, 180, 1, 1, 0.0_RKIND, interp_list, 1)
else
grid % soiltemp % array(iCell) = 0.0
end if
@@ -2595,12 +2595,12 @@
grid % snoalb % array(:) = 0.0
call map_set(PROJ_LATLON, proj, &
- latinc = 1.0, &
- loninc = 1.0, &
- knowni = 1.0, &
- knownj = 1.0, &
- lat1 = -89.5, &
- lon1 = -179.5)
+ latinc = 1.0_RKIND, &
+ loninc = 1.0_RKIND, &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = -89.5_RKIND, &
+ lon1 = -179.5_RKIND)
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)
@@ -2643,8 +2643,8 @@
end if
if (y < 1.0) y = 1.0
if (y > 179.0) y = 179.0
-! 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)
+! grid % snoalb % array(iCell) = interp_sequence(x, y, 1, maxsnowalb, 1, 360, 1, 180, 1, 1, -1.e30_RKIND, interp_list, 1)
+ grid % snoalb % array(iCell) = interp_sequence(x, y, 1, maxsnowalb, 1, 360, 1, 180, 1, 1, 0.0_RKIND, interp_list, 1)
else
grid % snoalb % array(iCell) = 0.0
end if
@@ -2673,12 +2673,12 @@
grid % greenfrac % array(:,:) = 0.0
call map_set(PROJ_LATLON, proj, &
- latinc = 0.144, &
- loninc = 0.144, &
- knowni = 1.0, &
- knownj = 1.0, &
- lat1 = -89.928, &
- lon1 = -179.928)
+ latinc = 0.144_RKIND, &
+ loninc = 0.144_RKIND, &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = -89.928_RKIND, &
+ lon1 = -179.928_RKIND)
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)
@@ -2715,7 +2715,7 @@
if (y < 1.0) y = 1.0
if (y > 1249.0) y = 1249.0
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)
+ grid % greenfrac % array(k,iCell) = interp_sequence(x, y, k, vegfra, 1, 2500, 1, 1250, 1, 12, -1.e30_RKIND, interp_list, 1)
end do
else
grid % greenfrac % array(:,iCell) = 0.0
@@ -2744,12 +2744,12 @@
grid % albedo12m % array(:,:) = 0.0
call map_set(PROJ_LATLON, proj, &
- latinc = 0.144, &
- loninc = 0.144, &
- knowni = 1.0, &
- knownj = 1.0, &
- lat1 = -89.928, &
- lon1 = -179.928)
+ latinc = 0.144_RKIND, &
+ loninc = 0.144_RKIND, &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = -89.928_RKIND, &
+ lon1 = -179.928_RKIND)
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)
@@ -2786,7 +2786,7 @@
if (y < 1.0) y = 1.0
if (y > 1249.0) y = 1249.0
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)
+ grid % albedo12m % array(k,iCell) = interp_sequence(x, y, k, vegfra, 1, 2500, 1, 1250, 1, 12, 0.0_RKIND, interp_list, 1)
end do
else
grid % albedo12m % array(:,iCell) = 8.0
@@ -2822,12 +2822,12 @@
if (field % iproj == PROJ_LATLON) then
call map_set(PROJ_LATLON, proj, &
- latinc = real(field % deltalat), &
- loninc = real(field % deltalon), &
- knowni = 1.0, &
- knownj = 1.0, &
- lat1 = real(field % startlat), &
- lon1 = real(field % startlon))
+ latinc = real(field % deltalat,RKIND), &
+ loninc = real(field % deltalon,RKIND), &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
end if
@@ -2848,9 +2848,9 @@
call latlon_to_ij(proj, lat, lon, x, y)
end if
if (ndims == 1) then
- destField1d(i) = interp_sequence(x, y, 1, field % slab, 1, field % nx, 1, field % ny, 1, 1, -1.e30, interp_list, 1)
+ destField1d(i) = interp_sequence(x, y, 1, field % slab, 1, field % nx, 1, field % ny, 1, 1, -1.e30_RKIND, interp_list, 1)
else if (ndims == 2) then
- destField2d(k,i) = interp_sequence(x, y, 1, field % slab, 1, field % nx, 1, field % ny, 1, 1, -1.e30, interp_list, 1)
+ destField2d(k,i) = interp_sequence(x, y, 1, field % slab, 1, field % nx, 1, field % ny, 1, 1, -1.e30_RKIND, interp_list, 1)
end if
end do
end if
@@ -3014,9 +3014,9 @@
hx(k,:) = hx(k-1,:)
dzminf = zw(k)-zw(k-1)
-! dzmin = max(.5,1.-.5*zw(k)/hm)
+! dzmin = max(0.5_RKIND,1.-.5*zw(k)/hm)
- sm = .05*min(.5*zw(k)/hm,1.)
+ sm = .05*min(0.5_RKIND*zw(k)/hm,1.0_RKIND)
do i=1,50
do iCell=1,grid %nCells
@@ -3290,18 +3290,18 @@
if (field % iproj == PROJ_LATLON) then
call map_set(PROJ_LATLON, proj, &
- latinc = real(field % deltalat), &
- loninc = real(field % deltalon), &
- knowni = 1.0, &
- knownj = 1.0, &
- lat1 = real(field % startlat), &
- lon1 = real(field % startlon))
+ latinc = real(field % deltalat,RKIND), &
+ loninc = real(field % deltalon,RKIND), &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
else if (field % iproj == PROJ_GAUSS) then
call map_set(PROJ_GAUSS, proj, &
nlat = nint(field % deltalat), &
- loninc = real(field % deltalon), &
- lat1 = real(field % startlat), &
- lon1 = real(field % startlon))
+ loninc = real(field % deltalon,RKIND), &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
! nxmax = nint(360.0 / field % deltalon), &
end if
@@ -3737,13 +3737,13 @@
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))
+ dx = real(field % dx * 1000.0,RKIND), &
+ truelat1 = real(field % truelat1,RKIND), &
+ stdlon = real(field % xlonc,RKIND), &
+ knowni = real(field % nx / 2.0,RKIND), &
+ knownj = real(field % ny / 2.0,RKIND), &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
end if
if (index(field % field, 'SEAICE') /= 0) then
@@ -4084,9 +4084,9 @@
if (config_theta_adv_order ==3) then
diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
- - sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ - sign(1.0_RKIND,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
- + sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ + sign(1.0_RKIND,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
end if
end do
@@ -4215,18 +4215,18 @@
if (field % iproj == PROJ_LATLON) then
call map_set(PROJ_LATLON, proj, &
- latinc = real(field % deltalat), &
- loninc = real(field % deltalon), &
- knowni = 1.0, &
- knownj = 1.0, &
- lat1 = real(field % startlat), &
- lon1 = real(field % startlon))
+ latinc = real(field % deltalat,RKIND), &
+ loninc = real(field % deltalon,RKIND), &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
else if (field % iproj == PROJ_GAUSS) then
call map_set(PROJ_GAUSS, proj, &
nlat = nint(field % deltalat), &
- loninc = real(field % deltalon), &
- lat1 = real(field % startlat), &
- lon1 = real(field % startlon))
+ loninc = real(field % deltalon,RKIND), &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
! nxmax = nint(360.0 / field % deltalon), &
end if
@@ -4247,7 +4247,7 @@
lon = lon - 360.0
call latlon_to_ij(proj, lat, lon, x, y)
end if
- fg % sst % array(iCell) = interp_sequence(x, y, 1, slab_r8, 1, field % nx, 1, field % ny, 1, 1, -1.e30, interp_list, 1)
+ fg % sst % array(iCell) = interp_sequence(x, y, 1, slab_r8, 1, field % nx, 1, field % ny, 1, 1, -1.e30_RKIND, interp_list, 1)
end do
deallocate(slab_r8)
@@ -4272,18 +4272,18 @@
if (field % iproj == PROJ_LATLON) then
call map_set(PROJ_LATLON, proj, &
- latinc = real(field % deltalat), &
- loninc = real(field % deltalon), &
- knowni = 1.0, &
- knownj = 1.0, &
- lat1 = real(field % startlat), &
- lon1 = real(field % startlon))
+ latinc = real(field % deltalat,RKIND), &
+ loninc = real(field % deltalon,RKIND), &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
else if (field % iproj == PROJ_GAUSS) then
call map_set(PROJ_GAUSS, proj, &
nlat = nint(field % deltalat), &
- loninc = real(field % deltalon), &
- lat1 = real(field % startlat), &
- lon1 = real(field % startlon))
+ loninc = real(field % deltalon,RKIND), &
+ lat1 = real(field % startlat,RKIND), &
+ lon1 = real(field % startlon,RKIND))
! nxmax = nint(360.0 / field % deltalon), &
end if
@@ -4304,7 +4304,7 @@
lon = lon - 360.0
call latlon_to_ij(proj, lat, lon, x, y)
end if
- fg % xice % array(iCell) = interp_sequence(x, y, 1, slab_r8, 1, field % nx, 1, field % ny, 1, 1, -1.e30, interp_list, 1)
+ fg % xice % array(iCell) = interp_sequence(x, y, 1, slab_r8, 1, field % nx, 1, field % ny, 1, 1, -1.e30_RKIND, interp_list, 1)
end do
@@ -4390,7 +4390,7 @@
#endif
- real function sphere_distance(lat1, lon1, lat2, lon2, radius)
+ real (kind=RKIND) function sphere_distance(lat1, lon1, lat2, lon2, radius)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
! sphere with given radius.
@@ -4433,12 +4433,12 @@
do while (nearest_cell /= current_cell)
current_cell = nearest_cell
- current_distance = sphere_distance(latCell(current_cell), lonCell(current_cell), target_lat, target_lon, 1.0)
+ current_distance = sphere_distance(latCell(current_cell), lonCell(current_cell), target_lat, target_lon, 1.0_RKIND)
nearest_cell = current_cell
nearest_distance = current_distance
do i = 1, nEdgesOnCell(current_cell)
iCell = cellsOnCell(i,current_cell)
- d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0)
+ d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0_RKIND)
if (d < nearest_distance) then
nearest_cell = iCell
nearest_distance = d
@@ -4476,13 +4476,13 @@
do while (nearest_edge /= current_edge)
current_edge = nearest_edge
- current_distance = sphere_distance(latEdge(current_edge), lonEdge(current_edge), target_lat, target_lon, 1.0)
+ current_distance = sphere_distance(latEdge(current_edge), lonEdge(current_edge), target_lat, target_lon, 1.0_RKIND)
nearest_edge = current_edge
nearest_distance = current_distance
cell1 = cellsOnEdge(1,current_edge)
cell2 = cellsOnEdge(2,current_edge)
- cell1_dist = sphere_distance(latCell(cell1), lonCell(cell1), target_lat, target_lon, 1.0)
- cell2_dist = sphere_distance(latCell(cell2), lonCell(cell2), target_lat, target_lon, 1.0)
+ cell1_dist = sphere_distance(latCell(cell1), lonCell(cell1), target_lat, target_lon, 1.0_RKIND)
+ cell2_dist = sphere_distance(latCell(cell2), lonCell(cell2), target_lat, target_lon, 1.0_RKIND)
if (cell1_dist < cell2_dist) then
iCell = cell1
else
@@ -4490,7 +4490,7 @@
end if
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i,iCell)
- d = sphere_distance(latEdge(iEdge), lonEdge(iEdge), target_lat, target_lon, 1.0)
+ d = sphere_distance(latEdge(iEdge), lonEdge(iEdge), target_lat, target_lon, 1.0_RKIND)
if (d < nearest_distance) then
nearest_edge = iEdge
nearest_distance = d
@@ -4501,7 +4501,7 @@
end function nearest_edge
- real function vertical_interp(target_z, nz, zf, order, extrap, surface_val, sealev_val)
+ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surface_val, sealev_val)
implicit none
Modified: branches/ocean_projects/time_averaging/src/core_nhyd_atmos/mpas_atm_advection.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_nhyd_atmos/mpas_atm_advection.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_nhyd_atmos/mpas_atm_advection.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,6 @@
module atm_advection
+ use mpas_kind_types
use mpas_grid_types
use mpas_configure
use mpas_constants
@@ -117,7 +118,7 @@
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
xc(2), yc(2), zc(2), &
- 0., 0., 1. )
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND )
! angles from cell center to neighbor centers (thetav)
@@ -390,7 +391,7 @@
! Computes the angle between arcs AB and AC, given points A, B, and C
! Equation numbers w.r.t. http://mathworld.wolfram.com/SphericalTrigonometry.html
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
+ real (kind=RKIND) function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
implicit none
@@ -410,9 +411,9 @@
real (kind=RKIND) :: s ! Semiperimeter of the triangle
real (kind=RKIND) :: sin_angle
- a = acos(max(min(bx*cx + by*cy + bz*cz,1.0),-1.0)) ! Eqn. (3)
- b = acos(max(min(ax*cx + ay*cy + az*cz,1.0),-1.0)) ! Eqn. (2)
- c = acos(max(min(ax*bx + ay*by + az*bz,1.0),-1.0)) ! Eqn. (1)
+ a = acos(max(min(bx*cx + by*cy + bz*cz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (3)
+ b = acos(max(min(ax*cx + ay*cy + az*cz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (2)
+ c = acos(max(min(ax*bx + ay*by + az*bz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (1)
ABx = bx - ax
ABy = by - ay
@@ -428,12 +429,12 @@
s = 0.5*(a + b + c)
! sin_angle = sqrt((sin(s-b)*sin(s-c))/(sin(b)*sin(c))) ! Eqn. (28)
- sin_angle = sqrt(min(1.,max(0.,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
+ sin_angle = sqrt(min(1.0_RKIND,max(0.0_RKIND,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
if ((Dx*ax + Dy*ay + Dz*az) >= 0.0) then
- sphere_angle = 2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ sphere_angle = 2.0 * asin(max(min(sin_angle,1.0_RKIND),-1.0_RKIND))
else
- sphere_angle = -2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ sphere_angle = -2.0 * asin(max(min(sin_angle,1.0_RKIND),-1.0_RKIND))
end if
end function sphere_angle
@@ -445,7 +446,7 @@
! Computes the angle between vectors AB and AC, given points A, B, and C, and
! a vector (u,v,w) normal to the plane.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
+ real (kind=RKIND) function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
implicit none
@@ -480,9 +481,9 @@
cos_angle = (ABx*ACx + ABy*ACy + ABz*ACz) / (mAB * mAC)
if ((Dx*u + Dy*v + Dz*w) >= 0.0) then
- plane_angle = acos(max(min(cos_angle,1.0),-1.0))
+ plane_angle = acos(max(min(cos_angle,1.0_RKIND),-1.0_RKIND))
else
- plane_angle = -acos(max(min(cos_angle,1.0),-1.0))
+ plane_angle = -acos(max(min(cos_angle,1.0_RKIND),-1.0_RKIND))
end if
end function plane_angle
@@ -495,7 +496,7 @@
! B=(bx, by, bz). It is assumed that both A and B lie on the surface of the
! same sphere centered at the origin.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function arc_length(ax, ay, az, bx, by, bz)
+ real (kind=RKIND) function arc_length(ax, ay, az, bx, by, bz)
implicit none
@@ -701,7 +702,7 @@
DO I = 1, N
C1= 0.0
DO J = 1, N
- C1 = AMAX1(C1,ABS(A(I,J)))
+ C1 = MAX(C1,ABS(A(I,J)))
END DO
C(I) = C1
END DO
@@ -838,7 +839,7 @@
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
xc(2), yc(2), zc(2), &
- 0., 0., 1. )
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND )
! angles from cell center to neighbor centers (thetav)
@@ -890,10 +891,10 @@
do i=2,n-1
ip1 = i+1
if (ip1 == n) ip1 = 1
- thetat(i) = plane_angle( 0.,0.,0., &
- xp(i)-xp(i-1), yp(i)-yp(i-1), 0., &
- xp(ip1)-xp(i), yp(ip1)-yp(i), 0., &
- 0., 0., 1.)
+ thetat(i) = plane_angle( 0.0_RKIND, 0.0_RKIND, 0.0_RKIND, &
+ xp(i)-xp(i-1), yp(i)-yp(i-1), 0.0_RKIND, &
+ xp(ip1)-xp(i), yp(ip1)-yp(i), 0.0_RKIND, &
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND)
thetat(i) = thetat(i) + thetat(i-1)
end do
Modified: branches/ocean_projects/time_averaging/src/core_nhyd_atmos/mpas_atm_test_cases.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_nhyd_atmos/mpas_atm_test_cases.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_nhyd_atmos/mpas_atm_test_cases.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -5,7 +5,9 @@
use mpas_constants
use mpas_dmpar
use atm_advection
+#ifdef DO_PHYSICS
use mpas_atmphys_control
+#endif
contains
@@ -94,6 +96,8 @@
stop
end if
+
+#ifdef DO_PHYSICS
!initialization of surface input variables technically not needed to run our current set of
!idealized test cases:
if (config_test_case > 0) then
@@ -105,6 +109,7 @@
end do
endif
+#endif
end subroutine atm_setup_test_case
@@ -153,8 +158,8 @@
integer :: eoe, j
integer, dimension(:), pointer :: nEdgesOnEdge
integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge
- real, dimension(:), pointer :: dvEdge, AreaCell
- real, dimension(:,:), pointer :: weightsOnEdge
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, AreaCell
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
real (kind=RKIND) :: u, v, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
@@ -666,7 +671,7 @@
if (config_test_case == 2) then
r_pert = sphere_distance( grid % latEdge % array (iEdge), grid % lonEdge % array (iEdge), &
- lat_pert, lon_pert, 1.)/(pert_radius)
+ lat_pert, lon_pert, 1.0_RKIND)/(pert_radius)
u_pert = u_perturbation*exp(-r_pert**2)*(lat2-lat1)*a/grid % dvEdge % array(iEdge)
else if (config_test_case == 3) then
@@ -786,9 +791,9 @@
if (config_theta_adv_order ==3) then
diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
- - sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ - sign(1.0_RKIND,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
- + sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ + sign(1.0_RKIND,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
end if
end do
@@ -881,7 +886,7 @@
! renormalize for setting cell-face fluxes
do k=1,nz1
- flux_zonal(k) = sign(1.,lat2_in-lat1_in)*flux_zonal(k)*dlat*a/dvEdge/u0
+ flux_zonal(k) = sign(1.0_RKIND,lat2_in-lat1_in)*flux_zonal(k)*dlat*a/dvEdge/u0
end do
end subroutine atm_calc_flux_zonal
@@ -1004,7 +1009,7 @@
integer :: eoe, j
integer, dimension(:), pointer :: nEdgesOnEdge
integer, dimension(:,:), pointer :: edgesOnEdge
- real, dimension(:,:), pointer :: weightsOnEdge
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, cell1, cell2, nCellsSolve
integer :: index_qv
@@ -1320,7 +1325,7 @@
temp = p(k,i)*thi(k,i)
pres = p0*p(k,i)**(1./rcp)
qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
- scalars(index_qv,k,i) = amin1(0.014,rh(k,i)*qvs)
+ scalars(index_qv,k,i) = min(0.014_RKIND,rh(k,i)*qvs)
end do
end do
@@ -1542,8 +1547,8 @@
integer :: eoe, j
integer, dimension(:), pointer :: nEdgesOnEdge
integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge
- real, dimension(:), pointer :: dvEdge, AreaCell, xCell, yCell
- real, dimension(:,:), pointer :: weightsOnEdge
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, AreaCell, xCell, yCell
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
integer :: index_qv
@@ -1769,7 +1774,7 @@
! smoothing grid for the upper level >> but not propoer for parallel programing
dzmin=.7
do k=2,nz1
- sm = .25*min((zc(k)-zc(k-1))/dz,1.)
+ sm = .25*min((zc(k)-zc(k-1))/dz,1.0_RKIND)
do i=1,grid % nCells
hx(k,i) = hx(k-1,i)
end do
@@ -1946,7 +1951,7 @@
temp = p(k,i)*t(k,i)
pres = p0*p(k,i)**(1./rcp)
qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
- scalars(index_qv,k,i) = amin1(0.014,rh(k,i)*qvs)
+ scalars(index_qv,k,i) = min(0.014_RKIND,rh(k,i)*qvs)
end do
do k=1,nz1
@@ -2069,9 +2074,9 @@
if (config_theta_adv_order ==3) then
diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
- - sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ - sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
- + sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ + sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
end if
end do
@@ -2128,7 +2133,7 @@
!----------------------------------------------------------------------------------------------------------
- real function sphere_distance(lat1, lon1, lat2, lon2, radius)
+ real (kind=RKIND) function sphere_distance(lat1, lon1, lat2, lon2, radius)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
! sphere with given radius.
@@ -2147,10 +2152,10 @@
end function sphere_distance
!--------------------------------------------------------------------
- real function env_qv( z, temperature, pressure, rh_max )
+ real (kind=RKIND) function env_qv( z, temperature, pressure, rh_max )
implicit none
- real z, temperature, pressure, ztr, es, qvs, p0, rh_max
+ real (kind=RKIND) :: z, temperature, pressure, ztr, es, qvs, p0, rh_max
p0 = 100000.
Modified: branches/ocean_projects/time_averaging/src/core_nhyd_atmos/mpas_atm_time_integration.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_nhyd_atmos/mpas_atm_time_integration.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_nhyd_atmos/mpas_atm_time_integration.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -680,7 +680,7 @@
type (mesh_type) :: grid
integer :: iCell, iEdge, k, cell1, cell2
integer, dimension(:,:), pointer :: cellsOnEdge
- real, dimension(:,:,:), pointer :: zf, zf3
+ real (kind=RKIND), dimension(:,:,:), pointer :: zf, zf3
real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, areaCell
real (kind=RKIND) :: flux
@@ -718,9 +718,9 @@
tend % w % array(k,cell1) = tend % w % array(k,cell1) - zf(k,1,iEdge)*flux
!3rd order stencil
if (config_theta_adv_order == 3) then
- tend % w % array(k,cell2) = tend % w % array(k,cell2) + sign(1.,tend % u % array(k,iEdge)) &
+ tend % w % array(k,cell2) = tend % w % array(k,cell2) + sign(1.0_RKIND,tend % u % array(k,iEdge)) &
*config_coef_3rd_order*zf3(k,2,iEdge)*flux
- tend % w % array(k,cell1) = tend % w % array(k,cell1) - sign(1.,tend % u % array(k,iEdge)) &
+ tend % w % array(k,cell1) = tend % w % array(k,cell1) - sign(1.0_RKIND,tend % u % array(k,iEdge)) &
*config_coef_3rd_order*zf3(k,1,iEdge)*flux
end if
@@ -869,7 +869,7 @@
do k=2,nVertLevels
- kr = min(nVertLevels,k+ nint(.5-sign(.5,zx(k,iEdge)+zx(k+1,iEdge))))
+ kr = min(nVertLevels,k+ nint(.5-sign(0.5_RKIND,zx(k,iEdge)+zx(k+1,iEdge))))
kl = min(nVertLevels,2*k+1-kr)
pr = zz(k,cell2)*rtheta_pp_old(k ,cell2)+.5*(zgrid(k ,cell1) +zgrid(k +1,cell1) &
-zgrid(k ,cell2) -zgrid(k +1,cell2)) &
@@ -1167,16 +1167,16 @@
!SHP-mtn
flux = cf1*ru(1,iEdge) + cf2*ru(2,iEdge) + cf3*ru(3,iEdge)
- w(1,cell2) = w(1,cell2) - (zb(1,2,iEdge) + sign(1.,flux)*coef_3rd_order*zb3(1,2,iEdge)) &
+ w(1,cell2) = w(1,cell2) - (zb(1,2,iEdge) + sign(1.0_RKIND,flux)*coef_3rd_order*zb3(1,2,iEdge)) &
*flux/(cf1*rho_zz(1,cell2)+cf2*rho_zz(2,cell2)+cf3*rho_zz(3,cell2))
- w(1,cell1) = w(1,cell1) + (zb(1,1,iEdge) + sign(1.,flux)*coef_3rd_order*zb3(1,1,iEdge)) &
+ w(1,cell1) = w(1,cell1) + (zb(1,1,iEdge) + sign(1.0_RKIND,flux)*coef_3rd_order*zb3(1,1,iEdge)) &
*flux/(cf1*rho_zz(1,cell1)+cf2*rho_zz(2,cell1)+cf3*rho_zz(3,cell1))
do k = 2, nVertLevels
flux = (fzm(k)*ru(k,iEdge)+fzp(k)*ru(k-1,iEdge))
- w(k,cell2) = w(k,cell2) - (zb(k,2,iEdge)+sign(1.,flux)*coef_3rd_order*zb3(k,2,iEdge)) &
+ w(k,cell2) = w(k,cell2) - (zb(k,2,iEdge)+sign(1.0_RKIND,flux)*coef_3rd_order*zb3(k,2,iEdge)) &
*flux/(fzm(k)*rho_zz(k,cell2)+fzp(k)*rho_zz(k-1,cell2))
- w(k,cell1) = w(k,cell1) + (zb(k,1,iEdge)+sign(1.,flux)*coef_3rd_order*zb3(k,1,iEdge)) &
+ w(k,cell1) = w(k,cell1) + (zb(k,1,iEdge)+sign(1.0_RKIND,flux)*coef_3rd_order*zb3(k,1,iEdge)) &
*flux/(fzm(k)*rho_zz(k,cell1)+fzp(k)*rho_zz(k-1,cell1))
enddo
@@ -1298,7 +1298,7 @@
do i=1,nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
do k=1,grid % nVertLevels
- scalar_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uhAvg(k,iEdge))*adv_coefs_3rd(i,iEdge)
+ scalar_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,uhAvg(k,iEdge))*adv_coefs_3rd(i,iEdge)
do iScalar=1,s_old % num_scalars
flux_arr(iScalar,k) = flux_arr(iScalar,k) + scalar_weight* scalar_new(iScalar,k,iCell)
end do
@@ -1633,7 +1633,7 @@
do i=1,nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
do k=1,grid % nVertLevels
- scalar_weight = uhAvg(k,iEdge)*(adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uhAvg(k,iEdge))*adv_coefs_3rd(i,iEdge))
+ scalar_weight = uhAvg(k,iEdge)*(adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,uhAvg(k,iEdge))*adv_coefs_3rd(i,iEdge))
flux_arr(k,iEdge) = flux_arr(k,iEdge) + scalar_weight* scalar_new(k,iCell)
end do
end do
@@ -1651,7 +1651,7 @@
do k = 2, nVertLevels
scalar_new(k,iCell) = scalar_old(k,iCell)*h_old(k,iCell)
- flux_upwind = dt*(max(0.,wwAvg(k,iCell))*scalar_old(k-1,iCell) + min(0.,wwAvg(k,iCell))*scalar_old(k,iCell))
+ flux_upwind = dt*(max(0.0_RKIND,wwAvg(k,iCell))*scalar_old(k-1,iCell) + min(0.0_RKIND,wwAvg(k,iCell))*scalar_old(k,iCell))
scalar_new(k-1,iCell) = scalar_new(k-1,iCell) - flux_upwind*rdnw(k-1)
scalar_new(k ,iCell) = scalar_new(k ,iCell) + flux_upwind*rdnw(k)
wdtn(k,iCell) = dt*wdtn(k,iCell) - flux_upwind
@@ -1661,8 +1661,8 @@
! contributions to the update: first the vertical flux component, then the horizontal
do k=1,nVertLevels
- scale_in (k,iCell) = - rdnw(k)*(min(0.,wdtn(k+1,iCell))-max(0.,wdtn(k,iCell)))
- scale_out(k,iCell) = - rdnw(k)*(max(0.,wdtn(k+1,iCell))-min(0.,wdtn(k,iCell)))
+ scale_in (k,iCell) = - rdnw(k)*(min(0.0_RKIND,wdtn(k+1,iCell))-max(0.0_RKIND,wdtn(k,iCell)))
+ scale_out(k,iCell) = - rdnw(k)*(max(0.0_RKIND,wdtn(k+1,iCell))-min(0.0_RKIND,wdtn(k,iCell)))
end do
end do
@@ -1677,15 +1677,15 @@
if (cell1 <= grid%nCellsSolve .or. cell2 <= grid%nCellsSolve) then ! only for owned cells
do k=1,grid % nVertLevels
flux_upwind = grid % dvEdge % array(iEdge) * dt * &
- (max(0.,uhAvg(k,iEdge))*scalar_old(k,cell1) + min(0.,uhAvg(k,iEdge))*scalar_old(k,cell2))
+ (max(0.0_RKIND,uhAvg(k,iEdge))*scalar_old(k,cell1) + min(0.0_RKIND,uhAvg(k,iEdge))*scalar_old(k,cell2))
flux_arr(k,iEdge) = dt*flux_arr(k,iEdge) - flux_upwind
scalar_new(k,cell1) = scalar_new(k,cell1) - flux_upwind / areaCell(cell1)
scalar_new(k,cell2) = scalar_new(k,cell2) + flux_upwind / areaCell(cell2)
- scale_out(k,cell1) = scale_out(k,cell1) - max(0.,flux_arr(k,iEdge)) / areaCell(cell1)
- scale_in (k,cell1) = scale_in (k,cell1) - min(0.,flux_arr(k,iEdge)) / areaCell(cell1)
- scale_out(k,cell2) = scale_out(k,cell2) + min(0.,flux_arr(k,iEdge)) / areaCell(cell2)
- scale_in (k,cell2) = scale_in (k,cell2) + max(0.,flux_arr(k,iEdge)) / areaCell(cell2)
+ scale_out(k,cell1) = scale_out(k,cell1) - max(0.0_RKIND,flux_arr(k,iEdge)) / areaCell(cell1)
+ scale_in (k,cell1) = scale_in (k,cell1) - min(0.0_RKIND,flux_arr(k,iEdge)) / areaCell(cell1)
+ scale_out(k,cell2) = scale_out(k,cell2) + min(0.0_RKIND,flux_arr(k,iEdge)) / areaCell(cell2)
+ scale_in (k,cell2) = scale_in (k,cell2) + max(0.0_RKIND,flux_arr(k,iEdge)) / areaCell(cell2)
end do
end if
@@ -1700,10 +1700,10 @@
s_upwind = scalar_new(k,iCell)/h_new(k,iCell)
scale_factor = (s_max(k,iCell)-s_upwind)/(s_max_update-s_upwind+eps)
- scale_in(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
+ scale_in(k,iCell) = min( 1.0_RKIND, max( 0.0_RKIND, scale_factor) )
scale_factor = (s_upwind-s_min(k,iCell))/(s_upwind-s_min_update+eps)
- scale_out(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
+ scale_out(k,iCell) = min( 1.0_RKIND, max( 0.0_RKIND, scale_factor) )
end do
end do
@@ -1729,8 +1729,8 @@
if (cell1 <= grid%nCellsSolve .or. cell2 <= grid%nCellsSolve) then
do k = 1, nVertLevels
flux = flux_arr(k,iEdge)
- flux = max(0.,flux) * min(scale_out(k,cell1), scale_in(k,cell2)) &
- + min(0.,flux) * min(scale_in(k,cell1), scale_out(k,cell2))
+ flux = max(0.0_RKIND,flux) * min(scale_out(k,cell1), scale_in(k,cell2)) &
+ + min(0.0_RKIND,flux) * min(scale_in(k,cell1), scale_out(k,cell2))
flux_arr(k,iEdge) = flux
end do
end if
@@ -1741,8 +1741,8 @@
do iCell=1,grid % nCells
do k = 2, nVertLevels
flux = wdtn(k,iCell)
- flux = max(0.,flux) * min(scale_out(k-1,iCell), scale_in(k ,iCell)) &
- + min(0.,flux) * min(scale_out(k ,iCell), scale_in(k-1,iCell))
+ flux = max(0.0_RKIND,flux) * min(scale_out(k-1,iCell), scale_in(k ,iCell)) &
+ + min(0.0_RKIND,flux) * min(scale_out(k ,iCell), scale_in(k-1,iCell))
wdtn(k,iCell) = flux
end do
end do
@@ -1790,7 +1790,7 @@
do iCell = 1, grid%nCells
do k=1, grid%nVertLevels
- scalar_new_in(iScalar,k,iCell) = max(0.,scalar_new(k,iCell))
+ scalar_new_in(iScalar,k,iCell) = max(0.0_RKIND,scalar_new(k,iCell))
end do
end do
@@ -2072,7 +2072,7 @@
k = 2
wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
do k=3,nVertLevels-1
- wduz(k) = flux3( u(k-2,iEdge),u(k-1,iEdge),u(k,iEdge),u(k+1,iEdge),0.5*(rw(k,cell1)+rw(k,cell2)), 1. )
+ wduz(k) = flux3( u(k-2,iEdge),u(k-1,iEdge),u(k,iEdge),u(k+1,iEdge),0.5*(rw(k,cell1)+rw(k,cell2)), 1.0_RKIND )
end do
k = nVertLevels
wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
@@ -2368,7 +2368,7 @@
do i=1,nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
do k=2,grid % nVertLevels
- scalar_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,ru_edge_w(k))*adv_coefs_3rd(i,iEdge)
+ scalar_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,ru_edge_w(k))*adv_coefs_3rd(i,iEdge)
flux_arr(k) = flux_arr(k) + scalar_weight* w(k,iCell)
end do
end do
@@ -2582,7 +2582,7 @@
k = 2
wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
do k=3,nVertLevels-1
- wdwz(k) = flux3( w(k-2,iCell),w(k-1,iCell),w(k,iCell),w(k+1,iCell),0.5*(rw(k,iCell)+rw(k-1,iCell)), 1. )
+ wdwz(k) = flux3( w(k-2,iCell),w(k-1,iCell),w(k,iCell),w(k+1,iCell),0.5*(rw(k,iCell)+rw(k-1,iCell)), 1.0_RKIND )
end do
k = nVertLevels
wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
@@ -2678,7 +2678,7 @@
do i=1,nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
do k=1,grid % nVertLevels
- scalar_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,ru(k,iEdge))*adv_coefs_3rd(i,iEdge)
+ scalar_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,ru(k,iEdge))*adv_coefs_3rd(i,iEdge)
flux_arr(k) = flux_arr(k) + scalar_weight* theta_m(k,iCell)
end do
end do
@@ -3324,10 +3324,10 @@
* (grid % fzp % array(k) * grid % zz % array(k-1,cell1) + grid % fzm % array(k) * grid % zz % array(k,cell1))
!3rd order! stencil
if (config_theta_adv_order ==3) then
- diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + sign(1.,flux)*config_coef_3rd_order &
+ diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + sign(1.0_RKIND,flux)*config_coef_3rd_order &
* grid % zb3 % array(k,2,iEdge)*flux &
* (grid % fzp % array(k) * grid % zz % array(k-1,cell2) + grid % fzm % array(k) * grid % zz % array(k,cell2))
- diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - sign(1.,flux)*config_coef_3rd_order &
+ diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - sign(1.0_RKIND,flux)*config_coef_3rd_order &
* grid % zb3 % array(k,1,iEdge)*flux &
* (grid % fzp % array(k) * grid % zz % array(k-1,cell1) + grid % fzm % array(k) * grid % zz % array(k,cell1))
end if
@@ -3405,11 +3405,11 @@
t(k) = state_new % theta_m % array(k,iCell)/(1. + 1.61*state_new % scalars % array(state_new % index_qv,k,iCell))
rho(k) = grid % zz % array(k,iCell)*state_new % rho_zz % array(k,iCell)
p(k) = diag % exner % array(k,iCell)
- qv(k) = max(0.,state_new % scalars % array(state_new % index_qv,k,iCell))
- qc(k) = max(0.,state_new % scalars % array(state_new % index_qc,k,iCell))
- qr(k) = max(0.,state_new % scalars % array(state_new % index_qr,k,iCell))
- qc1(k) = max(0.,state_old % scalars % array(state_old % index_qc,k,iCell))
- qr1(k) = max(0.,state_old % scalars % array(state_old % index_qr,k,iCell))
+ qv(k) = max(0.0_RKIND,state_new % scalars % array(state_new % index_qv,k,iCell))
+ qc(k) = max(0.0_RKIND,state_new % scalars % array(state_new % index_qc,k,iCell))
+ qr(k) = max(0.0_RKIND,state_new % scalars % array(state_new % index_qr,k,iCell))
+ qc1(k) = max(0.0_RKIND,state_old % scalars % array(state_old % index_qc,k,iCell))
+ qr1(k) = max(0.0_RKIND,state_old % scalars % array(state_old % index_qr,k,iCell))
dzu(k) = grid % dzu % array(k)
end do
@@ -3486,8 +3486,8 @@
do i=1,nx
do k=1,nz1
qrprod(k) = qc1t(k,i) &
- -(qc1t(k,i)-dt*amax1(ackess*(qc1(k,i)-.001), &
- 0.))/(1.+dt*ckess*qr1(k,i)**.875)
+ -(qc1t(k,i)-dt*max(ackess*(qc1(k,i)-.001), &
+ 0.0_RKIND))/(1.+dt*ckess*qr1(k,i)**.875)
velqr(k) = (qr1(k,i)*r(k))**1.1364*rhalf(k)
qvs(k) = pc(k)*exp(f2x*(pk(k)*t1t(k,i)-273.) &
/(pk(k)*t1t(k,i)- 36.))
@@ -3505,23 +3505,23 @@
artemp = 36340.*(.5*(velqr(2)+velqr(1))+veld-velu)
artot = artot+dt*artemp
do k=1,nz1
- qc1t(k,i) = amax1(qc1t(k,i)-qrprod(k),0.)
- qr1t(k,i) = amax1(qr1t(k,i)+qrprod(k),0.)
+ qc1t(k,i) = max(qc1t(k,i)-qrprod(k),0.0_RKIND)
+ qr1t(k,i) = max(qr1t(k,i)+qrprod(k),0.0_RKIND)
prod(k) = (qv1t(k,i)-qvs(k))/(1.+qvs(k)*f5 &
/(pk(k)*t1t(k,i)-36.)**2)
end do
do k=1,nz1
- ern(k) = amin1(dt*(((1.6+124.9*(r(k)*qr1t(k,i))**.2046) &
+ ern(k) = min(dt*(((1.6+124.9*(r(k)*qr1t(k,i))**.2046) &
*(r(k)*qr1t(k,i))**.525)/(2.55e6*pc(k) &
/(3.8 *qvs(k))+5.4e5))*(dim(qvs(k),qv1t(k,i)) &
/(r(k)*qvs(k))), &
- amax1(-prod(k)-qc1t(k,i),0.),qr1t(k,i))
+ max(-prod(k)-qc1t(k,i),0.0_RKIND),qr1t(k,i))
end do
do k=1,nz1
- buoycy(k) = f0(k)*(amax1(prod(k),-qc1t(k,i))-ern(k))
- qv1t(k,i) = amax1(qv1t(k,i) &
- -amax1(prod(k),-qc1t(k,i))+ern(k),0.)
- qc1t(k,i) = qc1t(k,i)+amax1(prod(k),-qc1t(k,i))
+ buoycy(k) = f0(k)*(max(prod(k),-qc1t(k,i))-ern(k))
+ qv1t(k,i) = max(qv1t(k,i) &
+ -max(prod(k),-qc1t(k,i))+ern(k),0.0_RKIND)
+ qc1t(k,i) = qc1t(k,i)+max(prod(k),-qc1t(k,i))
qr1t(k,i) = qr1t(k,i)-ern(k)
t1t (k,i) = t1t (k,i)+buoycy(k)
end do
Modified: branches/ocean_projects/time_averaging/src/core_ocean/mpas_ocn_advection.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_ocean/mpas_ocn_advection.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_ocean/mpas_ocn_advection.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,6 @@
module ocn_advection
+ use mpas_kind_types
use mpas_grid_types
use mpas_configure
use mpas_constants
@@ -117,7 +118,7 @@
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
xc(2), yc(2), zc(2), &
- 0., 0., 1. )
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND )
! angles from cell center to neighbor centers (thetav)
@@ -390,7 +391,7 @@
! Computes the angle between arcs AB and AC, given points A, B, and C
! Equation numbers w.r.t. http://mathworld.wolfram.com/SphericalTrigonometry.html
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
+ real (kind=RKIND) function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
implicit none
@@ -410,9 +411,9 @@
real (kind=RKIND) :: s ! Semiperimeter of the triangle
real (kind=RKIND) :: sin_angle
- a = acos(max(min(bx*cx + by*cy + bz*cz,1.0),-1.0)) ! Eqn. (3)
- b = acos(max(min(ax*cx + ay*cy + az*cz,1.0),-1.0)) ! Eqn. (2)
- c = acos(max(min(ax*bx + ay*by + az*bz,1.0),-1.0)) ! Eqn. (1)
+ a = acos(max(min(bx*cx + by*cy + bz*cz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (3)
+ b = acos(max(min(ax*cx + ay*cy + az*cz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (2)
+ c = acos(max(min(ax*bx + ay*by + az*bz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (1)
ABx = bx - ax
ABy = by - ay
@@ -428,12 +429,12 @@
s = 0.5*(a + b + c)
! sin_angle = sqrt((sin(s-b)*sin(s-c))/(sin(b)*sin(c))) ! Eqn. (28)
- sin_angle = sqrt(min(1.,max(0.,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
+ sin_angle = sqrt(min(1.0_RKIND,max(0.0_RKIND,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
if ((Dx*ax + Dy*ay + Dz*az) >= 0.0) then
- sphere_angle = 2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ sphere_angle = 2.0 * asin(max(min(sin_angle,1.0_RKIND),-1.0_RKIND))
else
- sphere_angle = -2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ sphere_angle = -2.0 * asin(max(min(sin_angle,1.0_RKIND),-1.0_RKIND))
end if
end function sphere_angle
@@ -445,7 +446,7 @@
! Computes the angle between vectors AB and AC, given points A, B, and C, and
! a vector (u,v,w) normal to the plane.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
+ real (kind=RKIND) function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
implicit none
@@ -480,9 +481,9 @@
cos_angle = (ABx*ACx + ABy*ACy + ABz*ACz) / (mAB * mAC)
if ((Dx*u + Dy*v + Dz*w) >= 0.0) then
- plane_angle = acos(max(min(cos_angle,1.0),-1.0))
+ plane_angle = acos(max(min(cos_angle,1.0_RKIND),-1.0_RKIND))
else
- plane_angle = -acos(max(min(cos_angle,1.0),-1.0))
+ plane_angle = -acos(max(min(cos_angle,1.0_RKIND),-1.0_RKIND))
end if
end function plane_angle
@@ -495,7 +496,7 @@
! B=(bx, by, bz). It is assumed that both A and B lie on the surface of the
! same sphere centered at the origin.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function arc_length(ax, ay, az, bx, by, bz)
+ real (kind=RKIND) function arc_length(ax, ay, az, bx, by, bz)
implicit none
@@ -839,7 +840,7 @@
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
xc(2), yc(2), zc(2), &
- 0., 0., 1. )
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND )
! angles from cell center to neighbor centers (thetav)
@@ -891,10 +892,10 @@
do i=2,n-1
ip1 = i+1
if (ip1 == n) ip1 = 1
- thetat(i) = plane_angle( 0.,0.,0., &
- xp(i)-xp(i-1), yp(i)-yp(i-1), 0., &
- xp(ip1)-xp(i), yp(ip1)-yp(i), 0., &
- 0., 0., 1.)
+ thetat(i) = plane_angle( 0.0_RKIND, 0.0_RKIND, 0.0_RKIND, &
+ xp(i)-xp(i-1), yp(i)-yp(i-1), 0.0_RKIND, &
+ xp(ip1)-xp(i), yp(ip1)-yp(i), 0.0_RKIND, &
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND)
thetat(i) = thetat(i) + thetat(i-1)
end do
Modified: branches/ocean_projects/time_averaging/src/core_ocean/mpas_ocn_vel_vadv.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_ocean/mpas_ocn_vel_vadv.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_ocean/mpas_ocn_vel_vadv.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -108,9 +108,9 @@
integer, dimension(:), pointer :: maxLevelEdgeTop
integer, dimension(:,:), pointer :: cellsOnEdge
- real :: wTopEdge
- real, dimension(:), allocatable :: w_dudzTopEdge
- real, dimension(:), pointer :: zMidZLevel
+ real (kind=RKIND) :: wTopEdge
+ real (kind=RKIND), dimension(:), allocatable :: w_dudzTopEdge
+ real (kind=RKIND), dimension(:), pointer :: zMidZLevel
if(.not.velVadvOn) return
Modified: branches/ocean_projects/time_averaging/src/core_sw/mpas_sw_advection.F
===================================================================
--- branches/ocean_projects/time_averaging/src/core_sw/mpas_sw_advection.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/core_sw/mpas_sw_advection.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,6 @@
module sw_advection
+ use mpas_kind_types
use mpas_grid_types
use mpas_configure
use mpas_constants
@@ -117,7 +118,7 @@
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
xc(2), yc(2), zc(2), &
- 0., 0., 1. )
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND )
! angles from cell center to neighbor centers (thetav)
@@ -390,7 +391,7 @@
! Computes the angle between arcs AB and AC, given points A, B, and C
! Equation numbers w.r.t. http://mathworld.wolfram.com/SphericalTrigonometry.html
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
+ real (kind=RKIND) function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
implicit none
@@ -410,9 +411,9 @@
real (kind=RKIND) :: s ! Semiperimeter of the triangle
real (kind=RKIND) :: sin_angle
- a = acos(max(min(bx*cx + by*cy + bz*cz,1.0),-1.0)) ! Eqn. (3)
- b = acos(max(min(ax*cx + ay*cy + az*cz,1.0),-1.0)) ! Eqn. (2)
- c = acos(max(min(ax*bx + ay*by + az*bz,1.0),-1.0)) ! Eqn. (1)
+ a = acos(max(min(bx*cx + by*cy + bz*cz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (3)
+ b = acos(max(min(ax*cx + ay*cy + az*cz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (2)
+ c = acos(max(min(ax*bx + ay*by + az*bz,1.0_RKIND),-1.0_RKIND)) ! Eqn. (1)
ABx = bx - ax
ABy = by - ay
@@ -428,12 +429,12 @@
s = 0.5*(a + b + c)
! sin_angle = sqrt((sin(s-b)*sin(s-c))/(sin(b)*sin(c))) ! Eqn. (28)
- sin_angle = sqrt(min(1.,max(0.,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
+ sin_angle = sqrt(min(1.0_RKIND,max(0.0_RKIND,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
if ((Dx*ax + Dy*ay + Dz*az) >= 0.0) then
- sphere_angle = 2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ sphere_angle = 2.0 * asin(max(min(sin_angle,1.0_RKIND),-1.0_RKIND))
else
- sphere_angle = -2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ sphere_angle = -2.0 * asin(max(min(sin_angle,1.0_RKIND),-1.0_RKIND))
end if
end function sphere_angle
@@ -445,7 +446,7 @@
! Computes the angle between vectors AB and AC, given points A, B, and C, and
! a vector (u,v,w) normal to the plane.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
+ real (kind=RKIND) function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
implicit none
@@ -480,9 +481,9 @@
cos_angle = (ABx*ACx + ABy*ACy + ABz*ACz) / (mAB * mAC)
if ((Dx*u + Dy*v + Dz*w) >= 0.0) then
- plane_angle = acos(max(min(cos_angle,1.0),-1.0))
+ plane_angle = acos(max(min(cos_angle,1.0_RKIND),-1.0_RKIND))
else
- plane_angle = -acos(max(min(cos_angle,1.0),-1.0))
+ plane_angle = -acos(max(min(cos_angle,1.0_RKIND),-1.0_RKIND))
end if
end function plane_angle
@@ -495,7 +496,7 @@
! B=(bx, by, bz). It is assumed that both A and B lie on the surface of the
! same sphere centered at the origin.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function arc_length(ax, ay, az, bx, by, bz)
+ real (kind=RKIND) function arc_length(ax, ay, az, bx, by, bz)
implicit none
@@ -838,7 +839,7 @@
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
xc(2), yc(2), zc(2), &
- 0., 0., 1. )
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND )
! angles from cell center to neighbor centers (thetav)
@@ -890,10 +891,10 @@
do i=2,n-1
ip1 = i+1
if (ip1 == n) ip1 = 1
- thetat(i) = plane_angle( 0.,0.,0., &
- xp(i)-xp(i-1), yp(i)-yp(i-1), 0., &
- xp(ip1)-xp(i), yp(ip1)-yp(i), 0., &
- 0., 0., 1.)
+ thetat(i) = plane_angle( 0.0_RKIND, 0.0_RKIND, 0.0_RKIND, &
+ xp(i)-xp(i-1), yp(i)-yp(i-1), 0.0_RKIND, &
+ xp(ip1)-xp(i), yp(ip1)-yp(i), 0.0_RKIND, &
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND)
thetat(i) = thetat(i) + thetat(i-1)
end do
Modified: branches/ocean_projects/time_averaging/src/framework/Makefile
===================================================================
--- branches/ocean_projects/time_averaging/src/framework/Makefile 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/framework/Makefile 2011-12-22 17:18:00 UTC (rev 1275)
@@ -4,7 +4,8 @@
ZOLTANOBJ = mpas_zoltan_interface.o
endif
-OBJS = mpas_framework.o \
+OBJS = mpas_kind_types.o \
+ mpas_framework.o \
mpas_timer.o \
mpas_timekeeping.o \
mpas_configure.o \
@@ -28,10 +29,18 @@
mpas_configure.o: mpas_dmpar.o
+mpas_constants.o: mpas_kind_types.o
+
mpas_grid_types.o: mpas_dmpar.o
-mpas_dmpar.o: mpas_sort.o streams.o
+mpas_dmpar.o: mpas_sort.o streams.o mpas_kind_types.o
+mpas_sort.o: mpas_kind_types.o
+
+mpas_timekeeping.o: mpas_kind_types.o
+
+mpas_timer.o: mpas_kind_types.o
+
mpas_block_decomp.o: mpas_grid_types.o mpas_hash.o mpas_configure.o
mpas_io_input.o: mpas_grid_types.o mpas_dmpar.o mpas_block_decomp.o mpas_sort.o mpas_configure.o mpas_timekeeping.o $(ZOLTANOBJ)
Modified: branches/ocean_projects/time_averaging/src/framework/mpas_constants.F
===================================================================
--- branches/ocean_projects/time_averaging/src/framework/mpas_constants.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/framework/mpas_constants.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,7 @@
module mpas_constants
+ use mpas_kind_types
+
real (kind=RKIND), parameter :: pii = 3.141592653589793
real (kind=RKIND), parameter :: a = 6371229.0
real (kind=RKIND), parameter :: omega = 7.29212e-5
Modified: branches/ocean_projects/time_averaging/src/framework/mpas_dmpar.F
===================================================================
--- branches/ocean_projects/time_averaging/src/framework/mpas_dmpar.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/framework/mpas_dmpar.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,15 +1,16 @@
module mpas_dmpar
+ use mpas_kind_types
use mpas_sort
#ifdef _MPI
include 'mpif.h'
integer, parameter :: MPI_INTEGERKIND = MPI_INTEGER
-#if (RKIND == 8)
+#ifdef SINGLE_PRECISION
+ integer, parameter :: MPI_REALKIND = MPI_REAL
+#else
integer, parameter :: MPI_REALKIND = MPI_DOUBLE_PRECISION
-#else
- integer, parameter :: MPI_REALKIND = MPI_REAL
#endif
#endif
@@ -19,6 +20,7 @@
type dm_info
integer :: nprocs, my_proc_id, comm, info
+ logical :: using_external_comm
end type dm_info
@@ -45,23 +47,30 @@
contains
- subroutine mpas_dmpar_init(dminfo)
+ subroutine mpas_dmpar_init(dminfo, mpi_comm)
implicit none
type (dm_info), intent(inout) :: dminfo
+ integer, intent(in), optional :: mpi_comm ! Optional: externally-supplied MPI communicator
#ifdef _MPI
integer :: mpi_rank, mpi_size
integer :: mpi_ierr
+ if (present(mpi_comm)) then
+ dminfo % comm = mpi_comm
+ dminfo % using_external_comm = .true.
+ else
+ call MPI_Init(mpi_ierr)
+ dminfo % comm = MPI_COMM_WORLD
+ dminfo % using_external_comm = .false.
+ end if
+
! Find out our rank and the total number of processors
- call MPI_Init(mpi_ierr)
- call MPI_Comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_ierr)
- call MPI_Comm_size(MPI_COMM_WORLD, mpi_size, mpi_ierr)
+ call MPI_Comm_rank(dminfo % comm, mpi_rank, mpi_ierr)
+ call MPI_Comm_size(dminfo % comm, mpi_size, mpi_ierr)
- dminfo % comm = MPI_COMM_WORLD
-
dminfo % nprocs = mpi_size
dminfo % my_proc_id = mpi_rank
@@ -75,6 +84,7 @@
dminfo % comm = 0
dminfo % my_proc_id = IO_NODE
dminfo % nprocs = 1
+ dminfo % using_external_comm = .false.
#endif
end subroutine mpas_dmpar_init
@@ -89,7 +99,9 @@
#ifdef _MPI
integer :: mpi_ierr
- call MPI_Finalize(mpi_ierr)
+ if (.not. dminfo % using_external_comm) then
+ call MPI_Finalize(mpi_ierr)
+ end if
#endif
end subroutine mpas_dmpar_finalize
Modified: branches/ocean_projects/time_averaging/src/framework/mpas_io_input.F
===================================================================
--- branches/ocean_projects/time_averaging/src/framework/mpas_io_input.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/framework/mpas_io_input.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1226,10 +1226,10 @@
#include "input_field0dreal.inc"
-#if (RKIND == 8)
+#if SINGLE_PRECISION
+ nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % scalar)
+#else
nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start1, count1, field % scalar)
-#else
- nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % scalar)
#endif
end subroutine mpas_io_input_field0d_real
@@ -1261,10 +1261,10 @@
#include "input_field1dreal.inc"
-#if (RKIND == 8)
+#if SINGLE_PRECISION
+ nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % array)
+#else
nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start1, count1, field % array)
-#else
- nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % array)
#endif
end subroutine mpas_io_input_field1d_real
@@ -1290,10 +1290,10 @@
#include "input_field2dreal.inc"
-#if (RKIND == 8)
+#if SINGLE_PRECISION
+ nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start2, count2, field % array)
+#else
nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start2, count2, field % array)
-#else
- nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start2, count2, field % array)
#endif
end subroutine mpas_io_input_field2d_real
@@ -1321,10 +1321,10 @@
#include "input_field3dreal.inc"
-#if (RKIND == 8)
+#if SINGLE_PRECISION
+ nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start3, count3, field % array)
+#else
nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start3, count3, field % array)
-#else
- nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start3, count3, field % array)
#endif
end subroutine mpas_io_input_field3d_real
@@ -1348,10 +1348,10 @@
#include "input_field0dreal_time.inc"
-#if (RKIND == 8)
+#if SINGLE_PRECISION
+ nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % scalar)
+#else
nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start1, count1, field % scalar)
-#else
- nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % scalar)
#endif
end subroutine mpas_io_input_field0d_real_time
@@ -1377,10 +1377,10 @@
#include "input_field1dreal_time.inc"
-#if (RKIND == 8)
+#if SINGLE_PRECISION
+ nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start2, count2, field % array)
+#else
nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start2, count2, field % array)
-#else
- nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start2, count2, field % array)
#endif
end subroutine mpas_io_input_field1d_real_time
@@ -1408,10 +1408,10 @@
#include "input_field2dreal_time.inc"
-#if (RKIND == 8)
+#if SINGLE_PRECISION
+ nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start3, count3, field % array)
+#else
nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start3, count3, field % array)
-#else
- nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start3, count3, field % array)
#endif
end subroutine mpas_io_input_field2d_real_time
@@ -1441,10 +1441,10 @@
#include "input_field3dreal_time.inc"
-#if (RKIND == 8)
+#if SINGLE_PRECISION
+ nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start4, count4, field % array)
+#else
nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start4, count4, field % array)
-#else
- nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start4, count4, field % array)
#endif
end subroutine mpas_io_input_field3d_real_time
Modified: branches/ocean_projects/time_averaging/src/framework/mpas_io_output.F
===================================================================
--- branches/ocean_projects/time_averaging/src/framework/mpas_io_output.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/framework/mpas_io_output.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -429,10 +429,10 @@
#include "output_field0dreal.inc"
-#if (RKIND == 8)
+#ifdef SINGLE_PRECISION
+ nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start1, count1, field % scalar)
+#else
nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start1, count1, field % scalar)
-#else
- nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start1, count1, field % scalar)
#endif
nferr = nf_sync(output_obj % wr_ncid)
@@ -458,10 +458,10 @@
#include "output_field1dreal.inc"
-#if (RKIND == 8)
+#ifdef SINGLE_PRECISION
+ nferr = nf_put_vara_real(output_obj % wr_ncid, VarID, start1, count1, field % array)
+#else
nferr = nf_put_vara_double(output_obj % wr_ncid, VarID, start1, count1, field % array)
-#else
- nferr = nf_put_vara_real(output_obj % wr_ncid, VarID, start1, count1, field % array)
#endif
nferr = nf_sync(output_obj % wr_ncid)
@@ -489,10 +489,10 @@
#include "output_field2dreal.inc"
-#if (RKIND == 8)
+#ifdef SINGLE_PRECISION
+ nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start2, count2, field % array)
+#else
nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start2, count2, field % array)
-#else
- nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start2, count2, field % array)
#endif
nferr = nf_sync(output_obj % wr_ncid)
@@ -522,10 +522,10 @@
#include "output_field3dreal.inc"
-#if (RKIND == 8)
+#ifdef SINGLE_PRECISION
+ nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start3, count3, field % array)
+#else
nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start3, count3, field % array)
-#else
- nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start3, count3, field % array)
#endif
nferr = nf_sync(output_obj % wr_ncid)
@@ -551,10 +551,10 @@
#include "output_field0dreal_time.inc"
-#if (RKIND == 8)
+#ifdef SINGLE_PRECISION
+ nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start1, count1, field % scalar)
+#else
nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start1, count1, field % scalar)
-#else
- nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start1, count1, field % scalar)
#endif
nferr = nf_sync(output_obj % wr_ncid)
@@ -582,10 +582,10 @@
#include "output_field1dreal_time.inc"
-#if (RKIND == 8)
+#ifdef SINGLE_PRECISION
+ nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start2, count2, field % array)
+#else
nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start2, count2, field % array)
-#else
- nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start2, count2, field % array)
#endif
nferr = nf_sync(output_obj % wr_ncid)
@@ -615,10 +615,10 @@
#include "output_field2dreal_time.inc"
-#if (RKIND == 8)
+#ifdef SINGLE_PRECISION
+ nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start3, count3, field % array)
+#else
nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start3, count3, field % array)
-#else
- nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start3, count3, field % array)
#endif
nferr = nf_sync(output_obj % wr_ncid)
@@ -650,10 +650,10 @@
#include "output_field3dreal_time.inc"
-#if (RKIND == 8)
+#ifdef SINGLE_PRECISION
+ nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start4, count4, field % array)
+#else
nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start4, count4, field % array)
-#else
- nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start4, count4, field % array)
#endif
nferr = nf_sync(output_obj % wr_ncid)
Copied: branches/ocean_projects/time_averaging/src/framework/mpas_kind_types.F (from rev 1274, trunk/mpas/src/framework/mpas_kind_types.F)
===================================================================
--- branches/ocean_projects/time_averaging/src/framework/mpas_kind_types.F (rev 0)
+++ branches/ocean_projects/time_averaging/src/framework/mpas_kind_types.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -0,0 +1,15 @@
+module mpas_kind_types
+
+#ifdef SINGLE_PRECISION
+ integer, parameter :: RKIND = selected_real_kind(6)
+#else
+ integer, parameter :: RKIND = selected_real_kind(12)
+#endif
+
+ contains
+
+ subroutine dummy_kinds()
+
+ end subroutine dummy_kinds
+
+end module mpas_kind_types
Modified: branches/ocean_projects/time_averaging/src/framework/mpas_sort.F
===================================================================
--- branches/ocean_projects/time_averaging/src/framework/mpas_sort.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/framework/mpas_sort.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,7 @@
module mpas_sort
+ use mpas_kind_types
+
interface quicksort
module procedure mpas_quicksort_int
module procedure mpas_quicksort_real
Modified: branches/ocean_projects/time_averaging/src/framework/mpas_timekeeping.F
===================================================================
--- branches/ocean_projects/time_averaging/src/framework/mpas_timekeeping.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/framework/mpas_timekeeping.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,6 @@
module mpas_timekeeping
+ use mpas_kind_types
use ESMF_BaseMod
use ESMF_Stubs
use ESMF_CalendarMod
Modified: branches/ocean_projects/time_averaging/src/framework/mpas_timer.F
===================================================================
--- branches/ocean_projects/time_averaging/src/framework/mpas_timer.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/framework/mpas_timer.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,7 @@
module mpas_timer
+ use mpas_kind_types
+
implicit none
save
! private
Modified: branches/ocean_projects/time_averaging/src/operators/mpas_spline_interpolation.F
===================================================================
--- branches/ocean_projects/time_averaging/src/operators/mpas_spline_interpolation.F 2011-12-22 15:30:21 UTC (rev 1274)
+++ branches/ocean_projects/time_averaging/src/operators/mpas_spline_interpolation.F 2011-12-22 17:18:00 UTC (rev 1275)
@@ -1,5 +1,7 @@
module mpas_spline_interpolation
+ use mpas_kind_types
+
implicit none
private
</font>
</pre>