Merge branch 'cbegeman/ocn/enable-e3sm-shared-pi' (PR #6485)

Enable E3SM shared pi in MPAS-Ocean This PR enables the use of E3SM shared constant pi in MPAS-Ocean when run in coupled configurations. [non-BFB]
E3SM-Project · Oct 23, 2024 · accee4c · accee4c
2 parents 4532660 + b3c8342
commit accee4c
Show file tree

Hide file tree

Showing 24 changed files with 56 additions and 45 deletions.
diff --git a/components/mpas-ocean/src/analysis_members/Makefile b/components/mpas-ocean/src/analysis_members/Makefile
@@ -2,6 +2,9 @@
 
 OBJS = mpas_ocn_analysis_driver.o
 
+UTILS = shr_kind_mod.o \
+        shr_const_mod.o
+
 MEMBERS = mpas_ocn_global_stats.o \
           mpas_ocn_okubo_weiss.o \
           mpas_ocn_layer_volume_weighted_averages.o \
@@ -35,7 +38,7 @@ MEMBERS = mpas_ocn_global_stats.o \
 
 all: $(OBJS)
 
-mpas_ocn_analysis_driver.o: $(MEMBERS)
+mpas_ocn_analysis_driver.o: $(UTILS) $(MEMBERS)
 
 mpas_ocn_okubo_weiss.o: mpas_ocn_okubo_weiss_eigenvalues.o
 

diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F
@@ -19,6 +19,9 @@
 
 module ocn_conservation_check
 
+   use shr_kind_mod, only: SHR_KIND_R8
+   use shr_const_mod
+
    use mpas_derived_types
    use mpas_pool_routines
    use mpas_dmpar
@@ -158,10 +161,9 @@ subroutine ocn_init_conservation_check(domain, err)!{{{
       !-----------------------------------------------------------------
 
       ! taking PI from SHR_CONST_PI in share/util/shr_const_mod.F90 to match coupler
-      !real (kind=RKIND), parameter :: piE3SM = 3.14159265358979323846_RKIND  ! pi
-      real (kind=RKIND), parameter :: piE3SM = 3.141592653589793_RKIND
+      real (kind=RKIND), parameter :: piE3SM = SHR_CONST_PI
       ! taking earth radius from SHR_CONST_REARTH in share/util/shr_const_mod.F90 to match coupler
-      real (kind=RKIND), parameter :: earthRadiusE3SM = 6.371229e6_RKIND ! radius of earth, m
+      real (kind=RKIND), parameter :: earthRadiusE3SM = SHR_CONST_REARTH ! radius of earth, m
 
       err = 0
 

diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F
@@ -681,7 +681,7 @@ SUBROUTINE harmonic_analysis_solve(MNP,nfreq,hmat,GLOELV,haff,haface,emagt,phase
       REAL(kind=RKIND),ALLOCATABLE :: hap(:),hax(:) 
       REAL(kind=RKIND),ALLOCATABLE :: ha(:,:)
 
-      convrd=180.0_RKIND/pii
+      convrd=180.0_RKIND/pi
 
       mm = 2*nfreq
       ALLOCATE ( PHASEE(nfreq),EMAG(nfreq) )

diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_lagrangian_particle_tracking.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_lagrangian_particle_tracking.F
@@ -3833,7 +3833,7 @@ subroutine convert_latlon_from_xyz(lat, lon, x, y, z) !{{{
 
     ! ensure range in 0, 2*pi
     if (lon < 0.0_RKIND) then
-      lon = 2*pii + lon
+      lon = 2*pi + lon
     end if
 
   end subroutine convert_latlon_from_xyz !}}}

diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_okubo_weiss.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_okubo_weiss.F
@@ -946,8 +946,8 @@ subroutine ocn_compute_eddy_stats(dminfo, block, nVertLevels, nCells, nCellsSolv
 
       ! for lat/lon coordinates, convert from radians to degrees for output
       if (config_AM_okuboWeiss_use_lat_lon_coords) then
-         wsPosX = wsPosX *180.0_RKIND/pii
-         wsPosY = wsPosY *180.0_RKIND/pii
+         wsPosX = wsPosX *180.0_RKIND/pi
+         wsPosY = wsPosY *180.0_RKIND/pi
       end if
 
       call mpas_timer_stop("stats per proc")

diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_time_filters.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_time_filters.F
@@ -62,7 +62,7 @@ module ocn_time_filters
    !--------------------------------------------------------------------
 #ifdef TIME_FILTERS_DEBUG
    integer :: iEdgeOutput = 0, iBlockOutput = 0, iklevel = 1
-   real (kind=RKIND) :: lonEdgePoint = 10.0_RKIND*pii/180.0_RKIND, latEdgePoint = 30.0_RKIND*pii/180.0_RKIND
+   real (kind=RKIND) :: lonEdgePoint = 10.0_RKIND*pi/180.0_RKIND, latEdgePoint = 30.0_RKIND*pi/180.0_RKIND
 #endif
 
 !***********************************************************************
@@ -182,7 +182,7 @@ subroutine ocn_init_time_filters(domain, err)!{{{
       call mpas_pool_get_subpool(block % structs, 'mesh', statePool)
       call mpas_pool_get_array(statePool, 'latEdge', latEdge)
       call mpas_pool_get_array(statePool, 'lonEdge', lonEdge)
-      print *,  'lon = ', 180.0_RKIND/pii*lonEdge(iEdgeOutput), ' lat = ', 180.0_RKIND/pii*latEdge(iEdgeOutput), &
+      print *,  'lon = ', 180.0_RKIND/pi*lonEdge(iEdgeOutput), ' lat = ', 180.0_RKIND/pi*latEdge(iEdgeOutput), &
                           ' iklevel=',iklevel, ' iEdgeOutput=',iEdgeOutput, ' iBlockOutput = ', iBlockOutput
 #endif
 

diff --git a/components/mpas-ocean/src/analysis_members/shr_const_mod.F b/components/mpas-ocean/src/analysis_members/shr_const_mod.F
@@ -0,0 +1 @@
+../../../../share/util/shr_const_mod.F90
diff --git a/components/mpas-ocean/src/analysis_members/shr_kind_mod.F b/components/mpas-ocean/src/analysis_members/shr_kind_mod.F
@@ -0,0 +1 @@
+../../../../share/util/shr_kind_mod.F90
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F
@@ -2587,7 +2587,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
 
                            ! Reset tracer2 to 2 in top n layers
                            ! in zonal bands, and 1 outside
-                           lat = latCell(iCell)*180./3.1415
+                           lat = latCell(iCell)*180.0_RKIND/pi
                            if (     lat>-60.0.and.lat<-55.0 &
                                 .or.lat>-40.0.and.lat<-35.0 &
                                 .or.lat>- 2.5.and.lat<  2.5 &
@@ -2604,7 +2604,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
 
                            ! Reset tracer3 to 2 in top n layers
                            ! in zonal bands, and 1 outside
-                           lat = latCell(iCell)*180./3.1415
+                           lat = latCell(iCell)*180.0_RKIND/pi
                            if (     lat>-55.0.and.lat<-50.0 &
                                 .or.lat>-35.0.and.lat<-30.0 &
                                 .or.lat>-15.0.and.lat<-10.0 &

diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
@@ -2225,7 +2225,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 
                            ! Reset tracer2 to 2 in top n layers
                            ! in zonal bands, and 1 outside
-                           lat = latCell(iCell)*180./3.1415
+                           lat = latCell(iCell)*180.0_RKIND/pi
                            if (     lat>-60.0.and.lat<-55.0 &
                                 .or.lat>-40.0.and.lat<-35.0 &
                                 .or.lat>- 2.5.and.lat<  2.5 &
@@ -2242,7 +2242,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 
                            ! Reset tracer3 to 2 in top n layers
                            ! in zonal bands, and 1 outside
-                           lat = latCell(iCell)*180./3.1415
+                           lat = latCell(iCell)*180.0_RKIND/pi
                            if (     lat>-55.0.and.lat<-50.0 &
                                 .or.lat>-35.0.and.lat<-30.0 &
                                 .or.lat>-15.0.and.lat<-10.0 &

diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F
@@ -2405,7 +2405,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{
 
                            ! Reset tracer2 to 2 in top n layers
                            ! in zonal bands, and 1 outside
-                           lat = latCell(iCell)*180./3.1415
+                           lat = latCell(iCell)*180.0_RKIND/pi
                            if (     lat>-60.0.and.lat<-55.0 &
                                 .or.lat>-40.0.and.lat<-35.0 &
                                 .or.lat>- 2.5.and.lat<  2.5 &
@@ -2422,7 +2422,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{
 
                            ! Reset tracer3 to 2 in top n layers
                            ! in zonal bands, and 1 outside
-                           lat = latCell(iCell)*180./3.1415
+                           lat = latCell(iCell)*180.0_RKIND/pi
                            if (     lat>-55.0.and.lat<-50.0 &
                                 .or.lat>-35.0.and.lat<-30.0 &
                                 .or.lat>-15.0.and.lat<-10.0 &

diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile
@@ -143,7 +143,7 @@ mpas_ocn_tracer_advection_std.o: mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_trac
 
 mpas_ocn_tracer_advection_vert.o: mpas_ocn_mesh.o mpas_ocn_config.o
 
-mpas_ocn_tracer_advection_shared.o: mpas_ocn_mesh.o mpas_ocn_config.o
+mpas_ocn_tracer_advection_shared.o: mpas_ocn_constants.o mpas_ocn_mesh.o mpas_ocn_config.o
 
 mpas_ocn_tracer_hmix_redi.o: mpas_ocn_constants.o mpas_ocn_config.o
 

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_constants.F b/components/mpas-ocean/src/shared/mpas_ocn_constants.F
@@ -65,6 +65,7 @@ module ocn_constants
 
    real (kind=RKIND), public :: &
       T0_Kelvin        ,&! zero point for Celsius
+      pi               ,&! pi
       mpercm           ,&! meters per m
       cmperm           ,&! m per meter
       days_per_second  ,&! days per second
@@ -117,6 +118,7 @@ subroutine ocn_constants_init(configPool, packagePool)!{{{
        !-----------------------------------------------------------------------
 
        T0_Kelvin = 273.16_RKIND               ! zero point for Celsius
+       pi        = 3.14159265358979323846_RKIND ! pi
        rho_air   = 1.2_RKIND                  ! ambient air density (kg/m^3)
        rho_sw    = config_density0            ! density of salt water (kg/m^3)
        rho_fw    = 1.0e3_RKIND                ! avg. water density (kg/m^3)
@@ -159,6 +161,7 @@ subroutine ocn_constants_init(configPool, packagePool)!{{{
 
 #ifdef MPAS_ESM_SHR_CONST
        T0_Kelvin              = SHR_CONST_TKFRZ         ! zero point for Celsius
+       pi                     = SHR_CONST_PI            ! zero point for Celsius
        cp_sw                  = SHR_CONST_CPSW          ! erg/g/K
        cp_air                 = SHR_CONST_CPDAIR        ! J/kg/K
        rho_air                = SHR_CONST_RHODAIR       ! kg/m^3

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_gm.F b/components/mpas-ocean/src/shared/mpas_ocn_gm.F
@@ -484,7 +484,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
 
                ! Compute the speed of the first baroclinic mode from the Brunt-Vaisala frequency.
                cGMphaseSpeed(iEdge) = max(c_min,                          &
-                                          sum_hN/(config_GM_spatially_variable_baroclinic_mode*3.141592_RKIND))
+                                          sum_hN/(config_GM_spatially_variable_baroclinic_mode*pi))
 
             end do
             !$omp end do
@@ -615,7 +615,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
                ! for the first-mode (m=1) baroclinic gravity wave speed in m/s
                ! c_m = 1/pi integral(N dz )
 
-               c_mode1 = max(1.0E-5_RKIND, sum_hN/3.141592_RKIND)
+               c_mode1 = max(1.0E-5_RKIND, sum_hN/pi)
 
                ! See Hallberg (2013) https://doi.org/10.1016/j.ocemod.2013.08.007
                ! just after eqn 4 the defines the deformation

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F b/components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F
@@ -217,8 +217,8 @@ subroutine ocn_manufactured_solution_init(domain, err)!{{{
 
       if (.not. config_use_manufactured_solution) return
 
-      kx = 2.0_RKIND*pii / config_manufactured_solution_wavelength_x
-      ky = 2.0_RKIND*pii / config_manufactured_solution_wavelength_y
+      kx = 2.0_RKIND*pi / config_manufactured_solution_wavelength_x
+      ky = 2.0_RKIND*pi / config_manufactured_solution_wavelength_y
 
       eta0 = config_manufactured_solution_amplitude
 

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_mesh.F b/components/mpas-ocean/src/shared/mpas_ocn_mesh.F
@@ -33,6 +33,7 @@ module ocn_mesh
    use mpas_log
 
    use ocn_config
+   use ocn_constants
 
    implicit none
    private
@@ -1792,7 +1793,7 @@ subroutine ocn_meshScaling() !{{{
                ! neighboring cells are circles for this calculation.
                cellWidth = 2.0_RKIND* &
                            sqrt((areaCell(cell1) + areaCell(cell2))/ &
-                                2.0_RKIND/pii)
+                                2.0_RKIND/pi)
                meshScalingDel2(iEdge) =  cellWidth/ &
                                          config_hmix_ref_cell_width
                meshScalingDel4(iEdge) = (cellWidth/ &

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F
@@ -1527,15 +1527,15 @@ subroutine ocn_tend_freq_filtered_thickness(tendPool, statePool, &
          do k = 1, maxLevelCell(iCell)
 
             tend_lowFreqDivergence(k,iCell) = &
-               -2.0*pii/thickness_filter_timescale_sec &
+               -2.0*pi/thickness_filter_timescale_sec &
               *(lowFreqDivergence(k,iCell)  - div_hu(k) &
                 + div_hu_btr * layerThickness(k,iCell)/totalThickness)
 
             tend_highFreqThickness(k,iCell) = - div_hu(k) + &
                   div_hu_btr*layerThickness(k,iCell)/totalThickness + &
                   lowFreqDivergence(k,iCell) + &
                   use_highFreqThick_restore* &
-                  (-2.0*pii/highFreqThick_restore_time_sec* &
+                  (-2.0*pi/highFreqThick_restore_time_sec* &
                    highFreqThickness(k,iCell) )
 
          end do ! vert (k) loop

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F b/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F
@@ -243,8 +243,8 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo
         ! compute the tidalHeight
         if (trim(config_tidal_forcing_model) == 'monochromatic') then
           tidalHeight = config_tidal_forcing_monochromatic_amp * &
-            SIN(2.0_RKIND*pii/config_tidal_forcing_monochromatic_period * daysSinceStartOfSim - &
-            pii*config_tidal_forcing_monochromatic_phaseLag/180.0_RKIND) - &
+            SIN(2.0_RKIND*pi/config_tidal_forcing_monochromatic_period * daysSinceStartOfSim - &
+            pi*config_tidal_forcing_monochromatic_phaseLag/180.0_RKIND) - &
             config_tidal_forcing_monochromatic_baseline
         elseif (trim(config_tidal_forcing_model) == 'linear') then
            tidalHeight = max(config_tidal_forcing_linear_min, &

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_shared.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_shared.F
@@ -24,10 +24,12 @@ module ocn_tracer_advection_shared
    use mpas_derived_types
    use mpas_hash
    use mpas_sort
+   use mpas_constants
    use mpas_geometry_utils
    use mpas_log
 
    use ocn_config
+   use ocn_constants
    use ocn_mesh
 
    implicit none
@@ -370,7 +372,6 @@ subroutine computeDerivTwo(derivTwo, err)!{{{
          xec, yec, zec,    &! arc bisection coords
          thetae_tmp,       &! angle
          xv1, xv2, yv1, yv2, zv1, zv2, &! vertex cart coords
-         pii,              &! pi math constant
          length_scale,     &! length scale
          cos2t, costsint, sin2t  ! trig function temps
 
@@ -413,7 +414,6 @@ subroutine computeDerivTwo(derivTwo, err)!{{{
 
       ! Initialize derivTwo and pi
 
-      pii = 2.*asin(1.0)
       derivTwo(:,:,:) = 0.0_RKIND
 
       do iCell = 1, nCellsAll
@@ -461,9 +461,9 @@ subroutine computeDerivTwo(derivTwo, err)!{{{
             end do
 
             if ( zc(1) == 1.0_RKIND) then
-               theta_abs(iCell) = pii/2.0_RKIND
+               theta_abs(iCell) = pi/2.0_RKIND
             else
-               theta_abs(iCell) = pii/2.0_RKIND &
+               theta_abs(iCell) = pi/2.0_RKIND &
                         - mpas_sphere_angle( xc(1), yc(1), zc(1), &
                                              xc(2), yc(2), zc(2), &
                                    0.0_RKIND, 0.0_RKIND, 1.0_RKIND) 
@@ -509,7 +509,7 @@ subroutine computeDerivTwo(derivTwo, err)!{{{
                angle_2d(i) = angleEdge(edgesOnCell(i,iCell))
                iEdge = edgesOnCell(i,iCell)
                if ( iCell .ne. cellsOnEdge(1,iEdge)) &
-                  angle_2d(i) = angle_2d(i) - pii
+                  angle_2d(i) = angle_2d(i) - pi
 
                xp(i) = dcEdge(edgesOnCell(i,iCell)) * cos(angle_2d(i))
                yp(i) = dcEdge(edgesOnCell(i,iCell)) * sin(angle_2d(i))

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F
@@ -199,7 +199,7 @@ subroutine ocn_vel_forcing_topographic_wave_drag_init(err)!{{{
               call mpas_log_write(" Topographic wave drag scheme is: Jayne and St. Laurent")
               call mpas_log_write("")
               tensorScheme = .false.
-              kappa = pii/10000.0_RKIND
+              kappa = pi/10000.0_RKIND
 
               do iEdge = 1, nEdgesAll
                   kmax =  maxLevelEdgeTop(iEdge)
@@ -241,7 +241,7 @@ subroutine ocn_vel_forcing_topographic_wave_drag_init(err)!{{{
 
                       topographicWaveDrag(iEdge) = topographicWaveDragCoeff * gam  &
                                                  * bed_slope_edges(iEdge)**2 * Nbar * Nb &
-                                                 / (8.0_RKIND * omegaM2 * pii**2)
+                                                 / (8.0_RKIND * omegaM2 * pi**2)
                   endif
               enddo
           else if (config_topographic_wave_drag_scheme.EQ."LGF") then
@@ -256,7 +256,7 @@ subroutine ocn_vel_forcing_topographic_wave_drag_init(err)!{{{
                       topographicWaveDrag(iEdge) = topographicWaveDragCoeff  &
                                                  * (sqrt((topo_buoyancy_N1B(iEdge)**2 - omegaM2**2)   &
                                                        * (topo_buoyancy_N1V(iEdge)**2 - omegaM2**2))) &
-                                                 / (4.0_RKIND*pii*omegaM2)
+                                                 / (4.0_RKIND*pi*omegaM2)
                       normalCoeffTWD(iEdge) = (lonGradEdge(iEdge)**2) &
                                          * (cos(angleEdge(iEdge))**2) + (latGradEdge(iEdge)**2)*(sin(angleEdge(iEdge))**2) &
                                          - 2.0_RKIND*latGradEdge(iEdge)*lonGradEdge(iEdge)*sin(angleEdge(iEdge))*(cos(angleEdge(iEdge)))

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_leith.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_leith.F
@@ -159,7 +159,7 @@ subroutine ocn_vel_hmix_leith_tend(div, relVort, tend, err)!{{{
          dcEdgeInv = 1.0_RKIND / dcEdge(iEdge)
          dvEdgeInv = 1.0_RKIND / dvEdge(iEdge)
 
-         visc2tmp = (leithParam*dxLeith*meshScalingDel2(iEdge)/pii)**3
+         visc2tmp = (leithParam*dxLeith*meshScalingDel2(iEdge)/pi)**3
 
          do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
 

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F
@@ -863,8 +863,8 @@ subroutine ocn_vel_self_attraction_loading_init(domain,err)!{{{
             ! Pre-compute sin and cos of latCell (co-latitude) values
             allocate(sinLatCell(nCellsAll), cosLatCell(nCellsAll))
             do iCell = 1,nCellsAll
-                sinLatCell(iCell) = sin(0.5_RKIND*pii-latCell(iCell))
-                cosLatCell(iCell) = cos(0.5_RKIND*pii-latCell(iCell))
+                sinLatCell(iCell) = sin(0.5_RKIND*pi-latCell(iCell))
+                cosLatCell(iCell) = cos(0.5_RKIND*pi-latCell(iCell))
             enddo
 
             ! Calculate blocking indices
@@ -2526,7 +2526,7 @@ subroutine associatedLegendrePolynomials(n, m, startIdx, endIdx, l, pmnm2, pmnm1
         if (n == m) then
 
             do iCell = startIdx,endIdx
-                pmnm2(iCell) = sqrt(1.0_RKIND/(4.0_RKIND*pii))*sinLatCell(iCell)**m
+                pmnm2(iCell) = sqrt(1.0_RKIND/(4.0_RKIND*pi))*sinLatCell(iCell)**m
                 do i = 1,m
                     pmnm2(iCell) = pmnm2(iCell)*sqrt(real(2*i+1,RKIND)/real(2*i,RKIND))
                 enddo