Clean up non-standard syntax and whitespace

src/ALE/MOM_ALE.F90

@@ -869,7 +869,7 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg,
         h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i+1,j,:) )
       endif
       if (associated(OBC)) then
-        if (OBC%segnum_u(I,j) .ne. 0) then
+        if (OBC%segnum_u(I,j) /= 0) then
           if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then
             h1(:) = h_old(i,j,:)
             h2(:) = h_new(i,j,:)
@@ -902,7 +902,7 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg,
         h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i,j+1,:) )
       endif
       if (associated(OBC)) then
-        if (OBC%segnum_v(i,J) .ne. 0) then
+        if (OBC%segnum_v(i,J) /= 0) then
           if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then
             h1(:) = h_old(i,j,:)
             h2(:) = h_new(i,j,:)

src/core/MOM.F90

@@ -1138,11 +1138,9 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
 
   endif ! -------------------------------------------------- end SPLIT
 
-   if (CS%do_dynamics) then!run particles whether or not stepping is split
-     if (CS%use_particles) then
-       call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, CS%tv) ! Run the particles model
-     endif
-   endif
+  if (CS%use_particles .and. CS%do_dynamics) then ! Run particles whether or not stepping is split
+    call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, CS%tv) ! Run the particles model
+  endif
 
 
   if (CS%thickness_diffuse .and. .not.CS%thickness_diffuse_first) then
@@ -3721,8 +3719,8 @@ subroutine MOM_end(CS)
   endif
 
   if (CS%use_particles) then
-     call particles_end(CS%particles)
-     deallocate(CS%particles)
+    call particles_end(CS%particles)
+    deallocate(CS%particles)
   endif
 
   call thickness_diffuse_end(CS%thickness_diffuse_CSp, CS%CDp)

src/core/MOM_continuity_PPM.F90

@@ -1282,7 +1282,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, por_fac
             l_seg = OBC%segnum_v(i,J)
 
             do_I(I) = .false.
-            if(l_seg /= OBC_NONE) &
+            if (l_seg /= OBC_NONE) &
               do_I(i) = (OBC%segment(l_seg)%specified)
 
             if (do_I(i)) FAvi(i) = GV%H_subroundoff*G%dx_Cv(i,J)

src/core/MOM_dynamics_split_RK2.F90

@@ -1368,12 +1368,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_h_PFu = register_diag_field('ocean_model', 'h_PFu', diag%axesCuL, Time, &
       'Thickness Multiplied Zonal Pressure Force Acceleration', &
       'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if(CS%id_h_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+  if (CS%id_h_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_h_PFv = register_diag_field('ocean_model', 'h_PFv', diag%axesCvL, Time, &
       'Thickness Multiplied Meridional Pressure Force Acceleration', &
       'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if(CS%id_h_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
+  if (CS%id_h_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, Time, &
       'Depth-integral of Zonal Pressure Force Acceleration', &
@@ -1398,12 +1398,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_h_CAu = register_diag_field('ocean_model', 'h_CAu', diag%axesCuL, Time, &
       'Thickness Multiplied Zonal Coriolis and Advective Acceleration', &
       'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if(CS%id_h_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+  if (CS%id_h_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_h_CAv = register_diag_field('ocean_model', 'h_CAv', diag%axesCvL, Time, &
       'Thickness Multiplied Meridional Coriolis and Advective Acceleration', &
       'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if(CS%id_h_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
+  if (CS%id_h_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, Time, &
       'Depth-integral of Zonal Coriolis and Advective Acceleration', &
@@ -1448,12 +1448,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_h_u_BT_accel = register_diag_field('ocean_model', 'h_u_BT_accel', diag%axesCuL, Time, &
       'Thickness Multiplied Barotropic Anomaly Zonal Acceleration', &
       'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if(CS%id_h_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+  if (CS%id_h_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_h_v_BT_accel = register_diag_field('ocean_model', 'h_v_BT_accel', diag%axesCvL, Time, &
       'Thickness Multiplied Barotropic Anomaly Meridional Acceleration', &
       'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if(CS%id_h_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
+  if (CS%id_h_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, Time, &
       'Depth-integral of Barotropic Anomaly Zonal Acceleration', &
@@ -1472,7 +1472,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, &
       'Meridional Pressure Force Acceleration multiplied by the viscous remnant', &
       'm s-2', conversion=US%L_T2_to_m_s2)
-  if(CS%id_PFv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
+  if (CS%id_PFv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, &
       'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', &
@@ -1481,7 +1481,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, &
       'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', &
       'm s-2', conversion=US%L_T2_to_m_s2)
-  if(CS%id_CAv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
+  if (CS%id_CAv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, &
       'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', &
@@ -1490,7 +1490,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, Time, &
       'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', &
       'm s-2', conversion=US%L_T2_to_m_s2)
-  if(CS%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
+  if (CS%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
 
   id_clock_Cor        = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
   id_clock_continuity = cpu_clock_id('(Ocean continuity equation)',      grain=CLOCK_MODULE)

src/core/MOM_dynamics_unsplit.F90

@@ -239,8 +239,6 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &
   h_av(:,:,:) = 0; hp(:,:,:) = 0
   up(:,:,:) = 0; upp(:,:,:) = 0
   vp(:,:,:) = 0; vpp(:,:,:) = 0
-  if (CS%id_ueffA > 0) ueffA(:,:,:) = 0
-  if (CS%id_veffA > 0) veffA(:,:,:) = 0
 
   dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
   if (dyn_p_surf) then
@@ -431,22 +429,23 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &
   call disable_averaging(CS%diag)
   call enable_averages(dt, Time_local, CS%diag)
 
-! Calculate effective areas and post data
+  ! Calculate effective areas and post data
   if (CS%id_ueffA > 0) then
-     do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-        if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k)/up(I,j,k)
-     enddo ; enddo ; enddo
-     call post_data(CS%id_ueffA, ueffA, CS%diag)
+    ueffA(:,:,:) = 0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k)/up(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_ueffA, ueffA, CS%diag)
   endif
 
   if (CS%id_veffA > 0) then
-     do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-        if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k)/vp(i,J,k)
-     enddo ; enddo ; enddo
-     call post_data(CS%id_veffA, veffA, CS%diag)
+    veffA(:,:,:) = 0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k)/vp(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_veffA, veffA, CS%diag)
   endif
 
-
 ! h_av = (h + hp)/2
   do k=1,nz
     do j=js-2,je+2 ; do i=is-2,ie+2

src/core/MOM_dynamics_unsplit_RK2.F90

@@ -250,8 +250,6 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
   h_av(:,:,:) = 0; hp(:,:,:) = 0
   up(:,:,:) = 0
   vp(:,:,:) = 0
-  if (CS%id_ueffA > 0) ueffA(:,:,:) = 0
-  if (CS%id_veffA > 0) veffA(:,:,:) = 0
 
   dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
   if (dyn_p_surf) then
@@ -452,17 +450,19 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
 
 ! Calculate effective areas and post data
   if (CS%id_ueffA > 0) then
-     do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-        if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k)/up(I,j,k)
-     enddo ; enddo ; enddo
-     call post_data(CS%id_ueffA, ueffA, CS%diag)
+    ueffA(:,:,:) = 0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k)/up(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_ueffA, ueffA, CS%diag)
   endif
 
   if (CS%id_veffA > 0) then
-     do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-        if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k)/vp(i,J,k)
-     enddo ; enddo ; enddo
-     call post_data(CS%id_veffA, veffA, CS%diag)
+    veffA(:,:,:) = 0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k)/vp(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_veffA, veffA, CS%diag)
   endif
 
 

src/core/MOM_open_boundary.F90

@@ -859,7 +859,7 @@ subroutine initialize_segment_data(G, OBC, PF)
               ! siz(3) is constituent for tidal variables
               call field_size(filename, 'constituent', siz, no_domain=.true.)
               ! expect third dimension to be number of constituents in MOM_input
-              if (siz(3) .ne. OBC%n_tide_constituents .and. OBC%add_tide_constituents) then
+              if (siz(3) /= OBC%n_tide_constituents .and. OBC%add_tide_constituents) then
                 call MOM_error(FATAL, 'Number of constituents in input data is not '//&
                     'the same as the number specified')
               endif
@@ -897,7 +897,7 @@ subroutine initialize_segment_data(G, OBC, PF)
         ! Check if this is a tidal field. If so, the number
         ! of expected constituents must be 1.
         if ((index(segment%field(m)%name, 'phase') > 0) .or. (index(segment%field(m)%name, 'amp') > 0)) then
-          if (OBC%n_tide_constituents .gt. 1 .and. OBC%add_tide_constituents) then
+          if (OBC%n_tide_constituents > 1 .and. OBC%add_tide_constituents) then
             call MOM_error(FATAL, 'Only one constituent is supported when specifying '//&
                 'tidal boundary conditions by value rather than file.')
           endif
@@ -997,7 +997,7 @@ subroutine initialize_obc_tides(OBC, US, param_file)
   ! If the nodal correction is based on a different time, initialize that.
   ! Otherwise, it can use N from the time reference.
   if (OBC%add_nodal_terms) then
-    if (sum(nodal_ref_date) .ne. 0) then
+    if (sum(nodal_ref_date) /= 0) then
       ! A reference date was provided for the nodal correction
       nodal_time = set_date(nodal_ref_date(1), nodal_ref_date(2), nodal_ref_date(3))
       call astro_longitudes_init(nodal_time, nodal_longitudes)
@@ -3939,7 +3939,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
         endif
         ! no dz for tidal variables
         if (segment%field(m)%nk_src > 1 .and.&
-            (index(segment%field(m)%name, 'phase') .le. 0 .and. index(segment%field(m)%name, 'amp') .le. 0)) then
+            (index(segment%field(m)%name, 'phase') <= 0 .and. index(segment%field(m)%name, 'amp') <= 0)) then
           call time_interp_external(segment%field(m)%fid_dz,Time, tmp_buffer_in)
           if (turns /= 0) then
             ! TODO: This is hardcoded for 90 degrees, and needs to be generalized.

src/core/MOM_variables.F90

@@ -316,10 +316,10 @@ module MOM_variables
 
 !> pointers to grids modifying cell metric at porous barriers
 type, public :: porous_barrier_ptrs
-   real, pointer, dimension(:,:,:) :: por_face_areaU => NULL() !< fractional open area of U-faces [nondim]
-   real, pointer, dimension(:,:,:) :: por_face_areaV => NULL() !< fractional open area of V-faces [nondim]
-   real, pointer, dimension(:,:,:) :: por_layer_widthU => NULL() !< fractional open width of U-faces [nondim]
-   real, pointer, dimension(:,:,:) :: por_layer_widthV => NULL() !< fractional open width of V-faces [nondim]
+  real, pointer, dimension(:,:,:) :: por_face_areaU => NULL() !< fractional open area of U-faces [nondim]
+  real, pointer, dimension(:,:,:) :: por_face_areaV => NULL() !< fractional open area of V-faces [nondim]
+  real, pointer, dimension(:,:,:) :: por_layer_widthU => NULL() !< fractional open width of U-faces [nondim]
+  real, pointer, dimension(:,:,:) :: por_layer_widthV => NULL() !< fractional open width of V-faces [nondim]
 end type porous_barrier_ptrs
 
 

src/framework/MOM_checksums.F90

@@ -110,7 +110,7 @@ subroutine chksum0(scalar, mesg, scale, logunit)
     call chksum_error(FATAL, 'NaN detected: '//trim(mesg))
 
   scaling = 1.0 ; if (present(scale)) scaling = scale
-  iounit = error_unit; if(present(logunit)) iounit = logunit
+  iounit = error_unit ; if (present(logunit)) iounit = logunit
 
   if (calculateStatistics) then
     rs = scaling * scalar
@@ -147,7 +147,7 @@ subroutine zchksum(array, mesg, scale, logunit)
   endif
 
   scaling = 1.0 ; if (present(scale)) scaling = scale
-  iounit = error_unit; if(present(logunit)) iounit = logunit
+  iounit = error_unit ; if (present(logunit)) iounit = logunit
 
   if (calculateStatistics) then
     if (present(scale)) then
@@ -352,7 +352,7 @@ subroutine chksum_h_2d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logu
   endif
 
   scaling = 1.0 ; if (present(scale)) scaling = scale
-  iounit = error_unit; if(present(logunit)) iounit = logunit
+  iounit = error_unit ; if (present(logunit)) iounit = logunit
 
   if (calculateStatistics) then
     if (present(scale)) then
@@ -618,7 +618,7 @@ subroutine chksum_B_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners,
   endif
 
   scaling = 1.0 ; if (present(scale)) scaling = scale
-  iounit = error_unit; if(present(logunit)) iounit = logunit
+  iounit = error_unit ; if (present(logunit)) iounit = logunit
   sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
   if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
@@ -901,7 +901,7 @@ subroutine chksum_u_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners,
   endif
 
   scaling = 1.0 ; if (present(scale)) scaling = scale
-  iounit = error_unit; if(present(logunit)) iounit = logunit
+  iounit = error_unit ; if (present(logunit)) iounit = logunit
   sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
   if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
@@ -1079,7 +1079,7 @@ subroutine chksum_v_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners,
   endif
 
   scaling = 1.0 ; if (present(scale)) scaling = scale
-  iounit = error_unit; if(present(logunit)) iounit = logunit
+  iounit = error_unit ; if (present(logunit)) iounit = logunit
   sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
   if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
@@ -1246,7 +1246,7 @@ subroutine chksum_h_3d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logu
   endif
 
   scaling = 1.0 ; if (present(scale)) scaling = scale
-  iounit = error_unit; if(present(logunit)) iounit = logunit
+  iounit = error_unit ; if (present(logunit)) iounit = logunit
 
   if (calculateStatistics) then
     if (present(scale)) then
@@ -1397,7 +1397,7 @@ subroutine chksum_B_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners,
   endif
 
   scaling = 1.0 ; if (present(scale)) scaling = scale
-  iounit = error_unit; if(present(logunit)) iounit = logunit
+  iounit = error_unit ; if (present(logunit)) iounit = logunit
   sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
   if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
@@ -1576,7 +1576,7 @@ subroutine chksum_u_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners,
   endif
 
   scaling = 1.0 ; if (present(scale)) scaling = scale
-  iounit = error_unit; if(present(logunit)) iounit = logunit
+  iounit = error_unit ; if (present(logunit)) iounit = logunit
   sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
   if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
@@ -1754,7 +1754,7 @@ subroutine chksum_v_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners,
   endif
 
   scaling = 1.0 ; if (present(scale)) scaling = scale
-  iounit = error_unit; if(present(logunit)) iounit = logunit
+  iounit = error_unit ; if (present(logunit)) iounit = logunit
   sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
   if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 

src/framework/MOM_document.F90

@@ -672,22 +672,22 @@ function real_array_string(vals, sep)
   integer :: j, n, ns
   logical :: doWrite
   character(len=10) :: separator
-  n=1 ; doWrite=.true. ; real_array_string=''
+  n = 1 ; doWrite = .true. ; real_array_string = ''
   if (present(sep)) then
-    separator=sep ; ns=len(sep)
+    separator = sep ; ns = len(sep)
   else
-    separator=', ' ; ns=2
+    separator = ', ' ; ns = 2
   endif
   do j=1,size(vals)
-    doWrite=.true.
-    if (j<size(vals)) then
-      if (vals(j)==vals(j+1)) then
-        n=n+1
-        doWrite=.false.
+    doWrite = .true.
+    if (j < size(vals)) then
+      if (vals(j) == vals(j+1)) then
+        n = n+1
+        doWrite = .false.
       endif
     endif
     if (doWrite) then
-      if(len(real_array_string)>0) then ! Write separator if a number has already been written
+      if (len(real_array_string) > 0) then ! Write separator if a number has already been written
         real_array_string = real_array_string // separator(1:ns)
       endif
       if (n>1) then