diff --git a/config_src/coupled_driver/MOM_surface_forcing.F90 b/config_src/coupled_driver/MOM_surface_forcing.F90 index a76d37cd6e..7d2af296e0 100644 --- a/config_src/coupled_driver/MOM_surface_forcing.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing.F90 @@ -305,7 +305,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & do j=js-2,je+2 ; do i=is-2,ie+2 fluxes%TKE_tidal(i,j) = CS%TKE_tidal(i,j) fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j) - enddo; enddo + enddo ; enddo if (restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) @@ -353,7 +353,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice do j=js,je ; do i=is,ie if (sfc_state%SST(i,j) <= -0.0539*sfc_state%SSS(i,j)) open_ocn_mask(i,j)=0.0 - enddo; enddo + enddo ; enddo endif if (CS%salt_restore_as_sflux) then do j=js,je ; do i=is,ie @@ -361,7 +361,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore) fluxes%salt_flux(i,j) = 1.e-3*G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const)* & (CS%basin_mask(i,j)*open_ocn_mask(i,j)*CS%srestore_mask(i,j)) *delta_sss ! kg Salt m-2 s-1 - enddo; enddo + enddo ; enddo if (CS%adjust_net_srestore_to_zero) then if (CS%adjust_net_srestore_by_scaling) then call adjust_area_mean_to_zero(fluxes%salt_flux, G, fluxes%saltFluxGlobalScl) @@ -382,7 +382,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & (CS%Rho0*CS%Flux_const) * & delta_sss / (0.5*(sfc_state%SSS(i,j) + data_restore(i,j))) endif - enddo; enddo + enddo ; enddo if (CS%adjust_net_srestore_to_zero) then if (CS%adjust_net_srestore_by_scaling) then call adjust_area_mean_to_zero(fluxes%vprec, G, fluxes%vPrecGlobalScl) @@ -392,7 +392,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & fluxes%vPrecGlobalAdj = reproducing_sum(work_sum(:,:), isr, ier, jsr, jer) / CS%area_surf do j=js,je ; do i=is,ie fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - fluxes%vPrecGlobalAdj ) * G%mask2dT(i,j) - enddo; enddo + enddo ; enddo endif endif endif @@ -536,12 +536,12 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & call adjust_area_mean_to_zero(net_FW2, G, fluxes%netFWGlobalScl) do j=js,je ; do i=is,ie fluxes%vprec(i,j) = fluxes%vprec(i,j) + (net_FW2(i,j) - net_FW(i,j)/G%areaT(i,j)) * G%mask2dT(i,j) - enddo; enddo + enddo ; enddo else fluxes%netFWGlobalAdj = reproducing_sum(net_FW(:,:), isr, ier, jsr, jer) / CS%area_surf do j=js,je ; do i=is,ie fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - fluxes%netFWGlobalAdj ) * G%mask2dT(i,j) - enddo; enddo + enddo ; enddo endif endif diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index dca6b8a837..395a4d3abb 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -1079,9 +1079,9 @@ subroutine ocean_model_data2D_get(OS,Ocean, name, array2D,isc,jsc) case('mask') array2D(isc:,jsc:) = OS%grid%mask2dT(g_isc:g_iec,g_jsc:g_jec) !OR same result -! do j=g_jsc,g_jec; do i=g_isc,g_iec +! do j=g_jsc,g_jec ; do i=g_isc,g_iec ! array2D(isc+i-g_isc,jsc+j-g_jsc) = OS%grid%mask2dT(i,j) -! enddo; enddo +! enddo ; enddo case('t_surf') array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET case('t_pme') diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 2727f42e1f..80a622b5ec 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -273,11 +273,11 @@ program MOM_main else calendar = uppercase(calendar) if (calendar(1:6) == 'JULIAN') then ; calendar_type = JULIAN - else if (calendar(1:9) == 'GREGORIAN') then ; calendar_type = GREGORIAN - else if (calendar(1:6) == 'NOLEAP') then ; calendar_type = NOLEAP - else if (calendar(1:10)=='THIRTY_DAY') then ; calendar_type = THIRTY_DAY_MONTHS - else if (calendar(1:11)=='NO_CALENDAR') then; calendar_type = NO_CALENDAR - else if (calendar(1:1) /= ' ') then + elseif (calendar(1:9) == 'GREGORIAN') then ; calendar_type = GREGORIAN + elseif (calendar(1:6) == 'NOLEAP') then ; calendar_type = NOLEAP + elseif (calendar(1:10)=='THIRTY_DAY') then ; calendar_type = THIRTY_DAY_MONTHS + elseif (calendar(1:11)=='NO_CALENDAR') then; calendar_type = NO_CALENDAR + elseif (calendar(1:1) /= ' ') then call MOM_error(FATAL,'MOM_driver: Invalid namelist value '//trim(calendar)//' for calendar') else call MOM_error(FATAL,'MOM_driver: No namelist value for calendar') @@ -641,7 +641,7 @@ program MOM_main call get_date(Time, yr, mon, day, hr, mins, sec) write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & 'Current model time: year, month, day, hour, minute, second' - end if + endif call close_file(unit) endif diff --git a/config_src/solo_driver/MOM_surface_forcing.F90 b/config_src/solo_driver/MOM_surface_forcing.F90 index 37bcaea17e..9a25b4b11a 100644 --- a/config_src/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/solo_driver/MOM_surface_forcing.F90 @@ -599,16 +599,16 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, CS) time_lev_daily = days - 365*floor(real(days) / 365.0) if (time_lev_daily < 31) then ; time_lev_monthly = 0 - else if (time_lev_daily < 59) then ; time_lev_monthly = 1 - else if (time_lev_daily < 90) then ; time_lev_monthly = 2 - else if (time_lev_daily < 120) then ; time_lev_monthly = 3 - else if (time_lev_daily < 151) then ; time_lev_monthly = 4 - else if (time_lev_daily < 181) then ; time_lev_monthly = 5 - else if (time_lev_daily < 212) then ; time_lev_monthly = 6 - else if (time_lev_daily < 243) then ; time_lev_monthly = 7 - else if (time_lev_daily < 273) then ; time_lev_monthly = 8 - else if (time_lev_daily < 304) then ; time_lev_monthly = 9 - else if (time_lev_daily < 334) then ; time_lev_monthly = 10 + elseif (time_lev_daily < 59) then ; time_lev_monthly = 1 + elseif (time_lev_daily < 90) then ; time_lev_monthly = 2 + elseif (time_lev_daily < 120) then ; time_lev_monthly = 3 + elseif (time_lev_daily < 151) then ; time_lev_monthly = 4 + elseif (time_lev_daily < 181) then ; time_lev_monthly = 5 + elseif (time_lev_daily < 212) then ; time_lev_monthly = 6 + elseif (time_lev_daily < 243) then ; time_lev_monthly = 7 + elseif (time_lev_daily < 273) then ; time_lev_monthly = 8 + elseif (time_lev_daily < 304) then ; time_lev_monthly = 9 + elseif (time_lev_daily < 334) then ; time_lev_monthly = 10 else ; time_lev_monthly = 11 endif @@ -847,16 +847,16 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, CS) time_lev_daily = days - 365*floor(real(days) / 365.0) if (time_lev_daily < 31) then ; time_lev_monthly = 0 - else if (time_lev_daily < 59) then ; time_lev_monthly = 1 - else if (time_lev_daily < 90) then ; time_lev_monthly = 2 - else if (time_lev_daily < 120) then ; time_lev_monthly = 3 - else if (time_lev_daily < 151) then ; time_lev_monthly = 4 - else if (time_lev_daily < 181) then ; time_lev_monthly = 5 - else if (time_lev_daily < 212) then ; time_lev_monthly = 6 - else if (time_lev_daily < 243) then ; time_lev_monthly = 7 - else if (time_lev_daily < 273) then ; time_lev_monthly = 8 - else if (time_lev_daily < 304) then ; time_lev_monthly = 9 - else if (time_lev_daily < 334) then ; time_lev_monthly = 10 + elseif (time_lev_daily < 59) then ; time_lev_monthly = 1 + elseif (time_lev_daily < 90) then ; time_lev_monthly = 2 + elseif (time_lev_daily < 120) then ; time_lev_monthly = 3 + elseif (time_lev_daily < 151) then ; time_lev_monthly = 4 + elseif (time_lev_daily < 181) then ; time_lev_monthly = 5 + elseif (time_lev_daily < 212) then ; time_lev_monthly = 6 + elseif (time_lev_daily < 243) then ; time_lev_monthly = 7 + elseif (time_lev_daily < 273) then ; time_lev_monthly = 8 + elseif (time_lev_daily < 304) then ; time_lev_monthly = 9 + elseif (time_lev_daily < 334) then ; time_lev_monthly = 10 else ; time_lev_monthly = 11 endif @@ -1153,7 +1153,7 @@ subroutine buoyancy_forcing_from_data_override(sfc_state, fluxes, day, dt, G, CS ! but evap is normally a positive quantity in the files fluxes%latent(i,j) = CS%latent_heat_vapor*fluxes%evap(i,j) fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j) - enddo; enddo + enddo ; enddo call data_override('OCN', 'sens', fluxes%sens(:,:), day, & is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in) @@ -1162,7 +1162,7 @@ subroutine buoyancy_forcing_from_data_override(sfc_state, fluxes, day, dt, G, CS do j=js,je ; do i=is,ie fluxes%sens(i,j) = -fluxes%sens(i,j) ! Normal convention is positive into the ocean ! but sensible is normally a positive quantity in the files - enddo; enddo + enddo ; enddo call data_override('OCN', 'sw', fluxes%sw(:,:), day, & is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index e4d297ddc8..c56d8a3fc3 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -1225,7 +1225,7 @@ subroutine ALE_initThicknessToCoord( CS, G, GV, h ) do j = G%jsd,G%jed ; do i = G%isd,G%ied h(i,j,:) = GV%m_to_H * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j) ) - enddo; enddo + enddo ; enddo end subroutine ALE_initThicknessToCoord diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 1f3488a7bc..ebe8b93bf6 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -1147,7 +1147,7 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) totalThickness = 0.0 do k = 1,nz totalThickness = totalThickness + h(i,j,k) - end do + enddo zOld(nz+1) = - nominalDepth do k = nz,1,-1 @@ -1190,8 +1190,8 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) call adjust_interface_motion( CS, nz, h(i,j,:), dzInterface(i,j,:) ) - end do - end do + enddo + enddo end subroutine build_zstar_grid @@ -1236,7 +1236,7 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) totalThickness = 0.0 do k = 1,nz totalThickness = totalThickness + h(i,j,k) - end do + enddo call build_sigma_column(CS%sigma_CS, nominalDepth, totalThickness, zNew) @@ -1244,7 +1244,7 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) zOld(nz+1) = -nominalDepth do k = nz,1,-1 zOld(k) = zOld(k+1) + h(i, j, k) - end do + enddo call filtered_grid_motion( CS, nz, zOld, zNew, dzInterface(i,j,:) ) @@ -1267,8 +1267,8 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) dzInterface(i,j,CS%nk+1) = 0. #endif - end do - end do + enddo + enddo end subroutine build_sigma_grid @@ -1393,8 +1393,8 @@ subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS ) endif #endif - end do ! end loop on i - end do ! end loop on j + enddo ! end loop on i + enddo ! end loop on j end subroutine build_rho_grid @@ -1466,7 +1466,7 @@ subroutine build_grid_HyCOM1( G, GV, h, tv, h_new, dzInterface, CS ) else ! on land dzInterface(i,j,:) = 0. endif ! mask2dT - enddo; enddo ! i,j + enddo ; enddo ! i,j call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) @@ -1597,7 +1597,7 @@ subroutine build_grid_SLight(G, GV, h, tv, dzInterface, CS) else ! on land dzInterface(i,j,:) = 0. endif ! mask2dT - enddo; enddo ! i,j + enddo ; enddo ! i,j end subroutine build_grid_SLight @@ -1704,7 +1704,7 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) total_height = 0.0 do k = 1,nz total_height = total_height + h(i,j,k) - end do + enddo eta = total_height - local_depth @@ -1715,7 +1715,7 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) z_inter(1) = eta do k = 1,nz z_inter(k+1) = z_inter(k) - delta_h - end do + enddo ! Refine grid in the middle do k = 1,nz+1 @@ -1725,15 +1725,15 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) if ( x <= x1 ) then t = y1*x/x1 - else if ( (x > x1 ) .and. ( x < x2 )) then + elseif ( (x > x1 ) .and. ( x < x2 )) then t = y1 + (y2-y1) * (x-x1) / (x2-x1) else t = y2 + (1.0-y2) * (x-x2) / (1.0-x2) - end if + endif z_inter(k) = -t * max_depth + eta - end do + enddo ! Modify interface heights to account for topography z_inter(nz+1) = - local_depth @@ -1742,8 +1742,8 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) do k = nz,1,-1 if ( z_inter(k) < (z_inter(k+1) + min_thickness) ) then z_inter(k) = z_inter(k+1) + min_thickness - end if - end do + endif + enddo ! Chnage in interface position x = 0. ! Left boundary at x=0 @@ -1751,11 +1751,11 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) do k = 2,nz x = x + h(i,j,k) dzInterface(i,j,k) = z_inter(k) - x - end do + enddo dzInterface(i,j,nz+1) = 0. - end do - end do + enddo + enddo stop 'OOOOOOPS' ! For some reason the gnu compiler will not let me delete this ! routine???? @@ -1792,17 +1792,17 @@ subroutine inflate_vanished_layers_old( CS, G, GV, h ) ! Build grid for current column do k = 1,GV%ke hTmp(k) = h(i,j,k) - end do + enddo call old_inflate_layers_1d( CS%min_thickness, GV%ke, hTmp ) ! Save modified grid do k = 1,GV%ke h(i,j,k) = hTmp(k) - end do + enddo - end do - end do + enddo + enddo end subroutine inflate_vanished_layers_old @@ -1859,7 +1859,7 @@ subroutine convective_adjustment(G, GV, h, tv) call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), & densities(k+1), tv%eqn_of_state ) stratified = .false. - end if + endif enddo ! k if ( stratified ) exit @@ -1962,7 +1962,7 @@ subroutine set_target_densities_from_GV( GV, CS ) CS%target_density(nz+1) = GV%Rlay(nz)+0.5*(GV%Rlay(nz)-GV%Rlay(nz-1)) do k = 2,nz CS%target_density(k) = CS%target_density(k-1) + CS%coordinateResolution(k) - end do + enddo CS%target_density_set = .true. end subroutine set_target_densities_from_GV @@ -2080,7 +2080,7 @@ function getCoordinateInterfaces( CS ) ! The following line has an "abs()" to allow ferret users to reference ! data by index. It is a temporary work around... :( -AJA getCoordinateInterfaces(:) = abs( getCoordinateInterfaces(:) ) - end if + endif end function getCoordinateInterfaces diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 10ba747d14..c0620122c1 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -140,7 +140,7 @@ subroutine buildGridFromH(nz, h, x) x(1) = 0.0 do k = 1,nz x(k+1) = x(k) + h(k) - end do + enddo end subroutine buildGridFromH @@ -389,21 +389,21 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & call PLM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) if ( CS%boundary_extrapolation ) then call PLM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect) - end if + endif iMethod = INTEGRATION_PLM case ( REMAPPING_PPM_H4 ) call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) if ( CS%boundary_extrapolation ) then call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) - end if + endif iMethod = INTEGRATION_PPM case ( REMAPPING_PPM_IH4 ) call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) if ( CS%boundary_extrapolation ) then call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) - end if + endif iMethod = INTEGRATION_PPM case ( REMAPPING_PQM_IH4IH3 ) call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge ) @@ -412,7 +412,7 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & if ( CS%boundary_extrapolation ) then call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, & ppoly_r_coefs, h_neglect ) - end if + endif iMethod = INTEGRATION_PQM case ( REMAPPING_PQM_IH6IH5 ) call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E, h_neglect_edge ) @@ -421,7 +421,7 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & if ( CS%boundary_extrapolation ) then call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, & ppoly_r_coefs, h_neglect ) - end if + endif iMethod = INTEGRATION_PQM case default call MOM_error( FATAL, 'MOM_remapping, build_reconstructions_1d: '//& @@ -1119,7 +1119,7 @@ subroutine remapByProjection( n0, h0, u0, ppoly0_E, ppoly0_coefs, & call integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, & xL, xR, h1(iTarget), u1(iTarget), jStart, xStart, h_neglect ) - end do ! end iTarget loop on target grid cells + enddo ! end iTarget loop on target grid cells end subroutine remapByProjection @@ -1206,7 +1206,7 @@ subroutine remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, dx1, & if (present(h1)) h1(iTarget) = hNew endif - end do ! end iTarget loop on target grid cells + enddo ! end iTarget loop on target grid cells end subroutine remapByDeltaZ @@ -1321,7 +1321,7 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, call MOM_error( FATAL,'The selected integration method is invalid' ) end select - end if ! end checking whether source cell is vanished + endif ! end checking whether source cell is vanished ! 2. Cell is not vanished else @@ -1454,8 +1454,8 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, do k = jL+1,jR-1 q = q + h0(k) * u0(k) hAct = hAct + h0(k) - end do - end if + enddo + endif ! Integrate from left boundary of cell jR up to xR xi0 = 0.0 @@ -1494,7 +1494,7 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, call MOM_error( FATAL,'The selected integration method is invalid' ) end select - end if ! end integration for non-vanished cells + endif ! end integration for non-vanished cells ! The cell average is the integrated value divided by the cell width #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__ @@ -1507,7 +1507,7 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, uAve = q / hC #endif - end if ! end if clause to check if cell is vanished + endif ! endif clause to check if cell is vanished end subroutine integrateReconOnInterval diff --git a/src/ALE/P1M_functions.F90 b/src/ALE/P1M_functions.F90 index 8590a7297f..75490bee9f 100644 --- a/src/ALE/P1M_functions.F90 +++ b/src/ALE/P1M_functions.F90 @@ -87,7 +87,7 @@ subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coef, h_neglect ) ppoly_coef(k,1) = u0_l ppoly_coef(k,2) = u0_r - u0_l - end do ! end loop on interior cells + enddo ! end loop on interior cells end subroutine P1M_interpolation @@ -147,7 +147,7 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef ) if ( (u1 - u0) * (ppoly_E(2,1) - u0_r) < 0.0 ) then slope = 2.0 * ( ppoly_E(2,1) - u0 ) - end if + endif ! Using the limited slope, the left edge value is reevaluated and ! the interpolant coefficients recomputed @@ -155,7 +155,7 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef ) ppoly_E(1,1) = u0 - 0.5 * slope else ppoly_E(1,1) = u0 - end if + endif ppoly_coef(1,1) = ppoly_E(1,1) ppoly_coef(1,2) = ppoly_E(1,2) - ppoly_E(1,1) @@ -175,13 +175,13 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef ) if ( (u1 - u0) * (u0_l - ppoly_E(N-1,2)) < 0.0 ) then slope = 2.0 * ( u1 - ppoly_E(N-1,2) ) - end if + endif if ( h1 /= 0.0 ) then ppoly_E(N,2) = u1 + 0.5 * slope else ppoly_E(N,2) = u1 - end if + endif ppoly_coef(N,1) = ppoly_E(N,1) ppoly_coef(N,2) = ppoly_E(N,2) - ppoly_E(N,1) diff --git a/src/ALE/P3M_functions.F90 b/src/ALE/P3M_functions.F90 index acc3e064ce..3034d2a8b4 100644 --- a/src/ALE/P3M_functions.F90 +++ b/src/ALE/P3M_functions.F90 @@ -144,7 +144,7 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) else h_l = h(k-1) u_l = u(k-1) - end if + endif if ( k == N ) then h_r = h(k) @@ -152,7 +152,7 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) else h_r = h(k+1) u_r = u(k+1) - end if + endif ! Compute limited slope sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hNeglect ) @@ -163,28 +163,28 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c ) else slope = 0.0 - end if + endif ! If the slopes are close to zero in machine precision and in absolute ! value, we set the slope to zero. This prevents asymmetric representation ! near extrema. These expressions are both nondimensional. if ( abs(u1_l*h_c) < eps ) then u1_l = 0.0 - end if + endif if ( abs(u1_r*h_c) < eps ) then u1_r = 0.0 - end if + endif ! The edge slopes are limited from above by the respective ! one-sided slopes if ( abs(u1_l) > abs(sigma_l) ) then u1_l = sigma_l - end if + endif if ( abs(u1_r) > abs(sigma_r) ) then u1_r = sigma_r - end if + endif ! Build cubic interpolant (compute the coefficients) call build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coef ) @@ -197,7 +197,7 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) ! cubic coefficients if ( monotonic == 0 ) then call monotonize_cubic( h_c, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ) - end if + endif ! Store edge slopes ppoly_S(k,1) = u1_l @@ -206,7 +206,7 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) ! Recompute coefficients of cubic call build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coef ) - end do ! loop on cells + enddo ! loop on cells end subroutine P3M_limiter @@ -278,7 +278,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & slope = 2.0 * ( u1 - u0 ) / ( h0 + hNeglect ) if ( abs(u1_r) > abs(slope) ) then u1_r = slope - end if + endif ! The right edge value in the boundary cell is taken to be the left ! edge value in the neighboring cell @@ -297,7 +297,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & u0_l = u0_r u1_l = 0.0 u1_r = 0.0 - end if + endif ! Store edge values and slope, build cubic and check monotonicity ppoly_E(i0,1) = u0_l @@ -317,7 +317,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & ppoly_S(i0,2) = u1_r call build_cubic_interpolant( h, i0, ppoly_E, ppoly_S, ppoly_coef ) - end if + endif ! ----- Right boundary ----- i0 = N-1 @@ -338,7 +338,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & slope = 2.0 * ( u1 - u0 ) / ( h1 + hNeglect ) if ( abs(u1_l) > abs(slope) ) then u1_l = slope - end if + endif ! The left edge value in the boundary cell is taken to be the right ! edge value in the neighboring cell @@ -357,7 +357,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & u0_r = u0_l u1_l = 0.0 u1_r = 0.0 - end if + endif ! Store edge values and slope, build cubic and check monotonicity ppoly_E(i1,1) = u0_l @@ -376,7 +376,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & ppoly_S(i1,2) = u1_r call build_cubic_interpolant( h, i1, ppoly_E, ppoly_S, ppoly_coef ) - end if + endif end subroutine P3M_boundary_extrapolation @@ -474,10 +474,10 @@ integer function is_cubic_monotonic( ppoly_coef, k ) if ( abs(c) > 1e-15 ) then xi_0 = 0.5 * ( -b - sqrt( rho ) ) / c xi_1 = 0.5 * ( -b + sqrt( rho ) ) / c - else if ( abs(b) > 1e-15 ) then + elseif ( abs(b) > 1e-15 ) then xi_0 = - a / b xi_1 = - a / b - end if + endif ! If one of the roots of the first derivative lies in (0,1), ! the cubic is not monotonic. @@ -486,11 +486,11 @@ integer function is_cubic_monotonic( ppoly_coef, k ) monotonic = 0 else monotonic = 1 - end if + endif else ! there are no real roots --> cubic is monotonic monotonic = 1 - end if + endif ! Set the return value is_cubic_monotonic = monotonic @@ -560,11 +560,11 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ! set them to zero if ( u1_l*slope <= 0.0 ) then u1_l = 0.0 - end if + endif if ( u1_r*slope <= 0.0 ) then u1_r = 0.0 - end if + endif ! Compute the location of the inflexion point, which is the root ! of the second derivative @@ -582,8 +582,8 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ! If the inflexion point lies in [0,1], change boolean value if ( (xi_ip >= 0.0) .AND. (xi_ip <= 1.0) ) then found_ip = 1 - end if - end if + endif + endif ! When there is an inflexion point within [0,1], check the slope ! to see if it is consistent with the limited PLM slope. If not, @@ -599,9 +599,9 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r inflexion_l = 1 else inflexion_r = 1 - end if - end if - end if ! found_ip + endif + endif + endif ! found_ip ! At this point, if the cubic is not monotonic, we know where the ! inflexion point should lie. When the cubic is monotonic, both @@ -618,12 +618,12 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r u1_l = 0.0 u1_r = 3.0 * (u0_r - u0_l) / h - else if (u1_l_tmp*slope < 0.0) then + elseif (u1_l_tmp*slope < 0.0) then u1_r = u1_r_tmp u1_l = 1.5*(u0_r - u0_l)/h - 0.5*u1_r - else if (u1_r_tmp*slope < 0.0) then + elseif (u1_r_tmp*slope < 0.0) then u1_l = u1_l_tmp u1_r = 3.0*(u0_r - u0_l)/h - 2.0*u1_l @@ -633,9 +633,9 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r u1_l = u1_l_tmp u1_r = u1_r_tmp - end if + endif - end if ! end treating case with inflexion point on the left + endif ! end treating case with inflexion point on the left ! Move inflexion point on the right if ( inflexion_r == 1 ) then @@ -648,12 +648,12 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r u1_l = 3.0 * (u0_r - u0_l) / h u1_r = 0.0 - else if (u1_l_tmp*slope < 0.0) then + elseif (u1_l_tmp*slope < 0.0) then u1_r = u1_r_tmp u1_l = 3.0*(u0_r - u0_l)/h - 2.0*u1_r - else if (u1_r_tmp*slope < 0.0) then + elseif (u1_r_tmp*slope < 0.0) then u1_l = u1_l_tmp u1_r = 1.5*(u0_r - u0_l)/h - 0.5*u1_l @@ -663,17 +663,17 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r u1_l = u1_l_tmp u1_r = u1_r_tmp - end if + endif - end if ! end treating case with inflexion point on the right + endif ! end treating case with inflexion point on the right if ( abs(u1_l*h) < eps ) then u1_l = 0.0 - end if + endif if ( abs(u1_r*h) < eps ) then u1_r = 0.0 - end if + endif end subroutine monotonize_cubic diff --git a/src/ALE/PCM_functions.F90 b/src/ALE/PCM_functions.F90 index bcb963faa6..6d407b0cc5 100644 --- a/src/ALE/PCM_functions.F90 +++ b/src/ALE/PCM_functions.F90 @@ -56,7 +56,7 @@ subroutine PCM_reconstruction( N, u, ppoly_E, ppoly_coef ) ! The edge values are equal to the cell average do k = 1,N ppoly_E(k,:) = u(k) - end do + enddo end subroutine PCM_reconstruction diff --git a/src/ALE/PLM_functions.F90 b/src/ALE/PLM_functions.F90 index 73f9206c21..12cd558e60 100644 --- a/src/ALE/PLM_functions.F90 +++ b/src/ALE/PLM_functions.F90 @@ -102,7 +102,7 @@ subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coef, h_neglect ) ! Extrema in the mean values require a PCM reconstruction avoid generating ! larger extreme values. slope = 0.0 - end if + endif ! This block tests to see if roundoff causes edge values to be out of bounds u_min = min( u_l, u_c, u_r ) @@ -130,7 +130,7 @@ subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coef, h_neglect ) ppoly_E(k,1) = u_c - 0.5 * slope ppoly_E(k,2) = u_c + 0.5 * slope - end do ! end loop on interior cells + enddo ! end loop on interior cells ! Boundary cells use PCM. Extrapolation is handled in a later routine. slp(1) = 0. diff --git a/src/ALE/PPM_functions.F90 b/src/ALE/PPM_functions.F90 index d0eb8325ad..11dabad684 100644 --- a/src/ALE/PPM_functions.F90 +++ b/src/ALE/PPM_functions.F90 @@ -198,7 +198,7 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef, h_neglect) slope = 2.0 * ( u1 - u0 ) if ( abs(u1_r) > abs(slope) ) then u1_r = slope - end if + endif ! The right edge value in the boundary cell is taken to be the left ! edge value in the neighboring cell @@ -215,11 +215,11 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef, h_neglect) if ( exp1 > exp2 ) then u0_l = 3.0 * u0 - 2.0 * u0_r - end if + endif if ( exp1 < -exp2 ) then u0_r = 3.0 * u0 - 2.0 * u0_l - end if + endif ppoly_E(i0,1) = u0_l ppoly_E(i0,2) = u0_r @@ -251,7 +251,7 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef, h_neglect) slope = 2.0 * ( u1 - u0 ) if ( abs(u1_l) > abs(slope) ) then u1_l = slope - end if + endif ! The left edge value in the boundary cell is taken to be the right ! edge value in the neighboring cell @@ -268,11 +268,11 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef, h_neglect) if ( exp1 > exp2 ) then u0_l = 3.0 * u1 - 2.0 * u0_r - end if + endif if ( exp1 < -exp2 ) then u0_r = 3.0 * u1 - 2.0 * u0_l - end if + endif ppoly_E(i1,1) = u0_l ppoly_E(i1,2) = u0_r diff --git a/src/ALE/PQM_functions.F90 b/src/ALE/PQM_functions.F90 index 6c89c7ac10..3a4e517e57 100644 --- a/src/ALE/PQM_functions.F90 +++ b/src/ALE/PQM_functions.F90 @@ -83,7 +83,7 @@ subroutine PQM_reconstruction( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ppoly_coef(k,4) = d ppoly_coef(k,5) = e - end do ! end loop on cells + enddo ! end loop on cells end subroutine PQM_reconstruction @@ -171,7 +171,7 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c ) else slope = 0.0 - end if + endif ! If one of the slopes has the wrong sign compared with the ! limited PLM slope, it is set equal to the limited PLM slope @@ -186,7 +186,7 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) u1_r = 0.0 inflexion_l = -1 inflexion_r = -1 - end if + endif ! Edge values are bounded and averaged when discontinuous and not ! monotonic, edge slopes are consistent and the cell is not an extremum. @@ -232,12 +232,12 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) inflexion_l = 1 else inflexion_r = 1 - end if - end if + endif + endif ! If both x1 and x2 do not lie in [0,1], check whether ! only x1 lies in [0,1] - else if ( (x1 >= 0.0) .AND. (x1 <= 1.0) ) then + elseif ( (x1 >= 0.0) .AND. (x1 <= 1.0) ) then gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b @@ -249,11 +249,11 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) inflexion_l = 1 else inflexion_r = 1 - end if - end if + endif + endif ! If x1 does not lie in [0,1], check whether x2 lies in [0,1] - else if ( (x2 >= 0.0) .AND. (x2 <= 1.0) ) then + elseif ( (x2 >= 0.0) .AND. (x2 <= 1.0) ) then gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b @@ -265,12 +265,12 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) inflexion_l = 1 else inflexion_r = 1 - end if - end if + endif + endif - end if ! end checking where the inflexion points lie + endif ! end checking where the inflexion points lie - end if ! end checking if alpha1 != 0 AND rho >= 0 + endif ! end checking if alpha1 != 0 AND rho >= 0 ! If alpha1 is zero, the second derivative of the quartic reduces ! to a straight line @@ -289,14 +289,14 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) inflexion_l = 1 else inflexion_r = 1 - end if - end if ! check slope consistency + endif + endif ! check slope consistency - end if + endif - end if ! end check whether we can find the root of the straight line + endif ! end check whether we can find the root of the straight line - end if ! end checking whether to shift inflexion points + endif ! end checking whether to shift inflexion points ! At this point, we know onto which edge to shift inflexion points if ( inflexion_l == 1 ) then @@ -316,15 +316,15 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) u0_r = 5.0 * u_c - 4.0 * u0_l u1_r = 20.0 * (u_c - u0_l) / ( h_c + hNeglect ) - else if ( u1_r * slope < 0.0 ) then + elseif ( u1_r * slope < 0.0 ) then u1_r = 0.0 u0_l = (5.0*u_c - 3.0*u0_r) / 2.0 u1_l = 10.0 * (-u_c + u0_r) / (3.0 * h_c + hNeglect) - end if + endif - else if ( inflexion_r == 1 ) then + elseif ( inflexion_r == 1 ) then ! We modify the edge slopes so that both inflexion points ! collapse onto the right edge @@ -341,15 +341,15 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) u0_r = ( 5.0 * u_c - 3.0 * u0_l ) / 2.0 u1_r = 10.0 * (u_c - u0_l) / (3.0 * h_c + hNeglect) - else if ( u1_r * slope < 0.0 ) then + elseif ( u1_r * slope < 0.0 ) then u1_r = 0.0 u0_l = 5.0 * u_c - 4.0 * u0_r u1_l = 20.0 * ( -u_c + u0_r ) / (h_c + hNeglect) - end if + endif - end if ! clause to check where to collapse inflexion points + endif ! clause to check where to collapse inflexion points ! Save edge values and edge slopes for reconstruction ppoly_E(k,1) = u0_l @@ -357,7 +357,7 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) ppoly_S(k,1) = u1_l ppoly_S(k,2) = u1_r - end do ! end loop on interior cells + enddo ! end loop on interior cells ! Constant reconstruction within boundary cells ppoly_E(1,:) = u(1) @@ -431,7 +431,7 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef ) slope = 2.0 * ( u1 - u0 ) if ( abs(u1_r) > abs(slope) ) then u1_r = slope - end if + endif ! The right edge value in the boundary cell is taken to be the left ! edge value in the neighboring cell @@ -448,11 +448,11 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef ) if ( exp1 > exp2 ) then u0_l = 3.0 * u0 - 2.0 * u0_r - end if + endif if ( exp1 < -exp2 ) then u0_r = 3.0 * u0 - 2.0 * u0_l - end if + endif ppoly_E(i0,1) = u0_l ppoly_E(i0,2) = u0_r @@ -489,7 +489,7 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef ) slope = 2.0 * ( u1 - u0 ) if ( abs(u1_l) > abs(slope) ) then u1_l = slope - end if + endif ! The left edge value in the boundary cell is taken to be the right ! edge value in the neighboring cell @@ -506,11 +506,11 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef ) if ( exp1 > exp2 ) then u0_l = 3.0 * u1 - 2.0 * u0_r - end if + endif if ( exp1 < -exp2 ) then u0_r = 3.0 * u1 - 2.0 * u0_l - end if + endif ppoly_E(i1,1) = u0_l ppoly_E(i1,2) = u0_r @@ -636,7 +636,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coef, else u0_l = u_plm u1_l = slope / (h0 + hNeglect) - end if + endif ! Monotonize quartic inflexion_l = 0 @@ -664,18 +664,18 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coef, gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b if ( gradient1 * slope < 0.0 ) then inflexion_l = 1 - end if - end if + endif + endif x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1 if ( (x2 > 0.0) .and. (x2 < 1.0) ) then gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b if ( gradient2 * slope < 0.0 ) then inflexion_l = 1 - end if - end if + endif + endif - end if + endif if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then @@ -684,10 +684,10 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coef, gradient1 = 3.0 * d * (x1**2) + 2.0 * c * x1 + b if ( gradient1 * slope < 0.0 ) then inflexion_l = 1 - end if - end if + endif + endif - end if + endif if ( inflexion_l == 1 ) then @@ -706,15 +706,15 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coef, u0_r = 5.0 * um - 4.0 * u0_l u1_r = 20.0 * (um - u0_l) / ( h0 + hNeglect ) - else if ( u1_r * slope < 0.0 ) then + elseif ( u1_r * slope < 0.0 ) then u1_r = 0.0 u0_l = (5.0*um - 3.0*u0_r) / 2.0 u1_l = 10.0 * (-um + u0_r) / (3.0 * h0 + hNeglect ) - end if + endif - end if + endif ! Store edge values, edge slopes and coefficients ppoly_E(i0,1) = u0_l @@ -789,7 +789,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coef, else u0_r = u_plm u1_r = slope / h1 - end if + endif ! Monotonize quartic inflexion_r = 0 @@ -817,18 +817,18 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coef, gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b if ( gradient1 * slope < 0.0 ) then inflexion_r = 1 - end if - end if + endif + endif x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1 if ( (x2 > 0.0) .and. (x2 < 1.0) ) then gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b if ( gradient2 * slope < 0.0 ) then inflexion_r = 1 - end if - end if + endif + endif - end if + endif if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then @@ -837,10 +837,10 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coef, gradient1 = 3.0 * d * (x1**2) + 2.0 * c * x1 + b if ( gradient1 * slope < 0.0 ) then inflexion_r = 1 - end if - end if + endif + endif - end if + endif if ( inflexion_r == 1 ) then @@ -859,15 +859,15 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coef, u0_r = ( 5.0 * um - 3.0 * u0_l ) / 2.0 u1_r = 10.0 * (um - u0_l) / (3.0 * h1) - else if ( u1_r * slope < 0.0 ) then + elseif ( u1_r * slope < 0.0 ) then u1_r = 0.0 u0_l = 5.0 * um - 4.0 * u0_r u1_l = 20.0 * ( -um + u0_r ) / h1 - end if + endif - end if + endif ! Store edge values, edge slopes and coefficients ppoly_E(i1,1) = u0_l diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index d3141cfd2d..84bb9e5518 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -123,14 +123,14 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & xTmp(1) = 0.0 do k = 1,count_nonzero_layers xTmp(k+1) = xTmp(k) + h_nv(k) - end do + enddo ! Compute densities on source column p(:) = CS%ref_pressure call calculate_density(T, S, p, densities, 1, nz, eqn_of_state) do k = 1,count_nonzero_layers densities(k) = densities(mapping(k)) - end do + enddo ! Based on source column density profile, interpolate to generate a new grid call build_and_interpolate_grid(CS%interp_CS, densities, count_nonzero_layers, & @@ -141,10 +141,10 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & call old_inflate_layers_1d(CS%min_thickness, CS%nk, h_new) ! Comment: The following adjustment of h_new, and re-calculation of h_new via x1 needs to be removed - x1(1) = 0.0 ; do k = 1,CS%nk ; x1(k+1) = x1(k) + h_new(k) ; end do + x1(1) = 0.0 ; do k = 1,CS%nk ; x1(k+1) = x1(k) + h_new(k) ; enddo do k = 1,CS%nk h_new(k) = x1(k+1) - x1(k) - end do + enddo else ! count_nonzero_layers <= 1 if (nz == CS%nk) then @@ -231,12 +231,12 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ if ( count_nonzero_layers <= 1 ) then h1(:) = h0(:) exit ! stop iterations here - end if + endif xTmp(1) = 0.0 do k = 1,count_nonzero_layers xTmp(k+1) = xTmp(k) + hTmp(k) - end do + enddo ! Compute densities within current water column call calculate_density( T_tmp, S_tmp, p, densities,& @@ -244,7 +244,7 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ do k = 1,count_nonzero_layers densities(k) = densities(mapping(k)) - end do + enddo ! One regridding iteration ! Based on global density profile, interpolate to generate a new grid @@ -252,12 +252,12 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ hTmp, xTmp, CS%target_density, nz, h1, x1, h_neglect, h_neglect_edge) call old_inflate_layers_1d( CS%min_thickness, nz, h1 ) - x1(1) = 0.0 ; do k = 1,nz ; x1(k+1) = x1(k) + h1(k) ; end do + x1(1) = 0.0 ; do k = 1,nz ; x1(k+1) = x1(k) + h1(k) ; enddo ! Remap T and S from previous grid to new grid do k = 1,nz h1(k) = x1(k+1) - x1(k) - end do + enddo call remapping_core_h(remapCS, nz, h0, S, nz, h1, Tmp, h_neglect, h_neglect_edge) S_tmp(:) = Tmp(:) @@ -273,7 +273,7 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ x0(k) = x0(k-1) + h0(k-1) x1(k) = x1(k-1) + h1(k-1) deviation = deviation + (x0(k)-x1(k))**2 - end do + enddo deviation = sqrt( deviation / (nz-1) ) m = m + 1 @@ -281,7 +281,7 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ ! Copy final grid onto start grid for next iteration h0(:) = h1(:) - end do ! end regridding iterations + enddo ! end regridding iterations if (CS%integrate_downward_for_e) then zInterface(1) = 0. @@ -330,12 +330,12 @@ subroutine copy_finite_thicknesses(nk, h_in, threshold, nout, h_out, mapping) if (h_out(nout) > thickest_h_out) then thickest_h_out = h_out(nout) k_thickest = nout - end if + endif else ! Add up mass in vanished layers thickness_in_vanished = thickness_in_vanished + h_in(k) - end if - end do + endif + enddo ! No finite layers if (nout <= 1) return @@ -367,8 +367,8 @@ subroutine old_inflate_layers_1d( min_thickness, nk, h ) do k = 1,nk if ( h(k) > min_thickness ) then count_nonzero_layers = count_nonzero_layers + 1 - end if - end do + endif + enddo ! If all layer thicknesses are greater than the threshold, exit routine if ( count_nonzero_layers == nk ) return @@ -377,9 +377,9 @@ subroutine old_inflate_layers_1d( min_thickness, nk, h ) if ( count_nonzero_layers == 0 ) then do k = 1,nk h(k) = min_thickness - end do + enddo return - end if + endif ! Inflate zero layers correction = 0.0 @@ -388,8 +388,8 @@ subroutine old_inflate_layers_1d( min_thickness, nk, h ) delta = min_thickness - h(k) correction = correction + delta h(k) = h(k) + delta - end if - end do + endif + enddo ! Modify thicknesses of nonzero layers to ensure volume conservation maxThickness = h(1) @@ -398,8 +398,8 @@ subroutine old_inflate_layers_1d( min_thickness, nk, h ) if ( h(k) > maxThickness ) then maxThickness = h(k) k_found = k - end if - end do + endif + enddo h(k_found) = h(k_found) - correction diff --git a/src/ALE/polynomial_functions.F90 b/src/ALE/polynomial_functions.F90 index 0cc4eb0b71..78c75f53a0 100644 --- a/src/ALE/polynomial_functions.F90 +++ b/src/ALE/polynomial_functions.F90 @@ -45,7 +45,7 @@ real function evaluation_polynomial( coeff, ncoef, x ) f = 0.0 do k = 1,ncoef f = f + coeff(k) * ( x**(k-1) ) - end do + enddo evaluation_polynomial = f @@ -73,7 +73,7 @@ real function first_derivative_polynomial( coeff, ncoef, x ) f = 0.0 do k = 2,ncoef f = f + REAL(k-1)*coeff(k) * ( x**(k-2) ) - end do + enddo first_derivative_polynomial = f @@ -99,7 +99,7 @@ real function integration_polynomial( xi0, xi1, Coeff, npoly ) do k = 1,npoly+1 integral = integral + Coeff(k) * (xi1**k - xi0**k) / real(k) - end do + enddo ! !One non-answer-changing way of unrolling the above is: ! k=1 diff --git a/src/ALE/regrid_edge_slopes.F90 b/src/ALE/regrid_edge_slopes.F90 index e07f3c3bd5..59d36e3e0e 100644 --- a/src/ALE/regrid_edge_slopes.F90 +++ b/src/ALE/regrid_edge_slopes.F90 @@ -116,23 +116,23 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect ) tri_b(i+1) = a * u(i) + b * u(i+1) - end do ! end loop on cells + enddo ! end loop on cells ! Boundary conditions: left boundary x(1) = 0.0 do i = 2,5 x(i) = x(i-1) + h(i-1) - end do + enddo do i = 1,4 do j = 1,4 Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - end do + enddo Bsys(i) = u(i) * ( h(i) ) - end do + enddo call solve_linear_system( Asys, Bsys, Csys, 4 ) @@ -148,17 +148,17 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect ) x(1) = 0.0 do i = 2,5 x(i) = x(i-1) + h(N-5+i) - end do + enddo do i = 1,4 do j = 1,4 Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - end do + enddo Bsys(i) = u(N-4+i) * ( h(N-4+i) ) - end do + enddo call solve_linear_system( Asys, Bsys, Csys, 4 ) @@ -176,7 +176,7 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect ) do i = 2,N edge_slopes(i,1) = tri_x(i) edge_slopes(i-1,2) = tri_x(i) - end do + enddo edge_slopes(1,1) = tri_x(1) edge_slopes(N,2) = tri_x(N+1) @@ -364,7 +364,7 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) tri_u(k+1) = beta tri_b(k+1) = a * u(k-1) + b * u(k) + c * u(k+1) + d * u(k+2) - end do ! end loop on cells + enddo ! end loop on cells ! Use a right-biased stencil for the second row @@ -481,17 +481,17 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) x(1) = 0.0 do i = 2,7 x(i) = x(i-1) + h(i-1) - end do + enddo do i = 1,6 do j = 1,6 Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - end do + enddo Bsys(i) = u(i) * h(i) - end do + enddo call solve_linear_system( Asys, Bsys, Csys, 6 ) @@ -621,17 +621,17 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) x(1) = 0.0 do i = 2,7 x(i) = x(i-1) + h(N-7+i) - end do + enddo do i = 1,6 do j = 1,6 Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - end do + enddo Bsys(i) = u(N-6+i) * h(N-6+i) - end do + enddo call solve_linear_system( Asys, Bsys, Csys, 6 ) @@ -652,7 +652,7 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) do i = 2,N edge_slopes(i,1) = tri_x(i) edge_slopes(i-1,2) = tri_x(i) - end do + enddo edge_slopes(1,1) = tri_x(1) edge_slopes(N,2) = tri_x(N+1) diff --git a/src/ALE/regrid_edge_values.F90 b/src/ALE/regrid_edge_values.F90 index d43cf5cc36..5fe4700c38 100644 --- a/src/ALE/regrid_edge_values.F90 +++ b/src/ALE/regrid_edge_values.F90 @@ -89,7 +89,7 @@ subroutine bound_edge_values( N, h, u, edge_val, h_neglect ) k0 = 1 k1 = 1 k2 = 2 - else if ( k == N ) then + elseif ( k == N ) then k0 = N-1 k1 = N k2 = N @@ -97,7 +97,7 @@ subroutine bound_edge_values( N, h, u, edge_val, h_neglect ) k0 = k-1 k1 = k k2 = k+1 - end if + endif ! All cells can now be treated equally h_l = h(k0) @@ -119,7 +119,7 @@ subroutine bound_edge_values( N, h, u, edge_val, h_neglect ) slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c ) else slope = 0.0 - end if + endif ! The limiter must be used in the local coordinate system to each cell. ! Hence, we must multiply the slope by h1. The multiplication by 0.5 is @@ -129,11 +129,11 @@ subroutine bound_edge_values( N, h, u, edge_val, h_neglect ) if ( (u_l-u0_l)*(u0_l-u_c) < 0.0 ) then u0_l = u_c - sign( min( abs(slope), abs(u0_l-u_c) ), slope ) - end if + endif if ( (u_r-u0_r)*(u0_r-u_c) < 0.0 ) then u0_r = u_c + sign( min( abs(slope), abs(u0_r-u_c) ), slope ) - end if + endif ! Finally bound by neighboring cell means in case of round off u0_l = max( min( u0_l, max(u_l, u_c) ), min(u_l, u_c) ) @@ -143,7 +143,7 @@ subroutine bound_edge_values( N, h, u, edge_val, h_neglect ) edge_val(k,1) = u0_l edge_val(k,2) = u0_r - end do ! loop on interior edges + enddo ! loop on interior edges end subroutine bound_edge_values @@ -178,9 +178,9 @@ subroutine average_discontinuous_edge_values( N, edge_val ) u0_avg = 0.5 * ( u0_minus + u0_plus ) edge_val(k,2) = u0_avg edge_val(k+1,1) = u0_avg - end if + endif - end do ! end loop on interior edges + enddo ! end loop on interior edges end subroutine average_discontinuous_edge_values @@ -224,9 +224,9 @@ subroutine check_discontinuous_edge_values( N, u, edge_val ) u0_avg = max( min( u0_avg, max(um_minus, um_plus) ), min(um_minus, um_plus) ) edge_val(k,2) = u0_avg edge_val(k+1,1) = u0_avg - end if + endif - end do ! end loop on interior edges + enddo ! end loop on interior edges end subroutine check_discontinuous_edge_values @@ -284,7 +284,7 @@ subroutine edge_values_explicit_h2( N, h, u, edge_val, h_neglect ) ! value of left cell edge_val(k-1,2) = edge_val(k,1) - end do ! end loop on interior cells + enddo ! end loop on interior cells ! Boundary edge values are simply equal to the boundary cell averages edge_val(1,1) = u(1) @@ -388,24 +388,24 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect ) endif #endif - end do ! end loop on interior cells + enddo ! end loop on interior cells ! Determine first two edge values f1 = max( hNeglect, hMinFrac*sum(h(1:4)) ) x(1) = 0.0 do i = 2,5 x(i) = x(i-1) + max(f1, h(i-1)) - end do + enddo do i = 1,4 do j = 1,4 A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) - end do + enddo B(i) = u(i) * max(f1, h(i) ) - end do + enddo call solve_linear_system( A, B, C, 4 ) @@ -433,17 +433,17 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect ) x(1) = 0.0 do i = 2,5 x(i) = x(i-1) + max(f1, h(N-5+i)) - end do + enddo do i = 1,4 do j = 1,4 A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) - end do + enddo B(i) = u(N-4+i) * max(f1, h(N-4+i) ) - end do + enddo call solve_linear_system( A, B, C, 4 ) @@ -461,10 +461,10 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect ) do i = 1,4 do j = 1,4 A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) - end do + enddo write(0,*) A(i,:) B(i) = u(N-4+i) * ( h(N-4+i) ) - end do + enddo write(0,*) 'B=',B write(0,*) 'C=',C write(0,*) 'h(:N)=',h(N-3:N) @@ -561,24 +561,24 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect ) tri_b(i+1) = a * u(i) + b * u(i+1) - end do ! end loop on cells + enddo ! end loop on cells ! Boundary conditions: left boundary h0 = max( hNeglect, hMinFrac*sum(h(1:4)) ) x(1) = 0.0 do i = 2,5 x(i) = x(i-1) + max( h0, h(i-1) ) - end do + enddo do i = 1,4 do j = 1,4 Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - end do + enddo Bsys(i) = u(i) * max( h0, h(i) ) - end do + enddo call solve_linear_system( Asys, Bsys, Csys, 4 ) @@ -591,17 +591,17 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect ) x(1) = 0.0 do i = 2,5 x(i) = x(i-1) + max( h0, h(N-5+i) ) - end do + enddo do i = 1,4 do j = 1,4 Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - end do + enddo Bsys(i) = u(N-4+i) * max( h0, h(N-4+i) ) - end do + enddo call solve_linear_system( Asys, Bsys, Csys, 4 ) @@ -615,7 +615,7 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect ) do i = 2,N edge_val(i,1) = tri_x(i) edge_val(i-1,2) = tri_x(i) - end do + enddo edge_val(1,1) = tri_x(1) edge_val(N,2) = tri_x(N+1) @@ -812,7 +812,7 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect ) tri_u(k+1) = beta tri_b(k+1) = a * u(k-1) + b * u(k) + c * u(k+1) + d * u(k+2) - end do ! end loop on cells + enddo ! end loop on cells ! Use a right-biased stencil for the second row @@ -940,17 +940,17 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect ) x(1) = 0.0 do i = 2,7 x(i) = x(i-1) + max( g, h(i-1) ) - end do + enddo do i = 1,6 do j = 1,6 Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - end do + enddo Bsys(i) = u(i) * max( g, h(i) ) - end do + enddo call solve_linear_system( Asys, Bsys, Csys, 6 ) @@ -1085,17 +1085,17 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect ) x(1) = 0.0 do i = 2,7 x(i) = x(i-1) + max( g, h(N-7+i) ) - end do + enddo do i = 1,6 do j = 1,6 Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - end do + enddo Bsys(i) = u(N-6+i) * max( g, h(N-6+i) ) - end do + enddo call solve_linear_system( Asys, Bsys, Csys, 6 ) @@ -1110,7 +1110,7 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect ) do i = 2,N edge_val(i,1) = tri_x(i) edge_val(i-1,2) = tri_x(i) - end do + enddo edge_val(1,1) = tri_x(1) edge_val(N,2) = tri_x(N+1) diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index d9d2a19228..fd445e7318 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -104,7 +104,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs ) - end if + endif case ( INTERPOLATION_P1M_H4 ) degree = DEGREE_1 @@ -112,11 +112,11 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge ) else call edge_values_explicit_h2( n0, h0, densities, ppoly0_E, h_neglect_edge ) - end if + endif call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs ) - end if + endif case ( INTERPOLATION_P1M_IH4 ) degree = DEGREE_1 @@ -124,18 +124,18 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge ) else call edge_values_explicit_h2( n0, h0, densities, ppoly0_E, h_neglect_edge ) - end if + endif call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs ) - end if + endif case ( INTERPOLATION_PLM ) degree = DEGREE_1 call PLM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call PLM_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) - end if + endif case ( INTERPOLATION_PPM_H4 ) if ( n0 >= 4 ) then @@ -145,15 +145,15 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & if (extrapolate) then call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, & ppoly0_coefs, h_neglect ) - end if + endif else degree = DEGREE_1 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E, h_neglect_edge ) call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs ) - end if - end if + endif + endif case ( INTERPOLATION_PPM_IH4 ) if ( n0 >= 4 ) then @@ -163,15 +163,15 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & if (extrapolate) then call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, & ppoly0_coefs, h_neglect ) - end if + endif else degree = DEGREE_1 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E, h_neglect_edge ) call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs ) - end if - end if + endif + endif case ( INTERPOLATION_P3M_IH4IH3 ) if ( n0 >= 4 ) then @@ -183,15 +183,15 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & if (extrapolate) then call P3M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect, h_neglect_edge ) - end if + endif else degree = DEGREE_1 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E, h_neglect_edge ) call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs ) - end if - end if + endif + endif case ( INTERPOLATION_P3M_IH6IH5 ) if ( n0 >= 6 ) then @@ -203,15 +203,15 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & if (extrapolate) then call P3M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect, h_neglect_edge ) - end if + endif else degree = DEGREE_1 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E, h_neglect_edge ) call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs ) - end if - end if + endif + endif case ( INTERPOLATION_PQM_IH4IH3 ) if ( n0 >= 4 ) then @@ -223,15 +223,15 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & if (extrapolate) then call PQM_boundary_extrapolation_v1( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect ) - end if + endif else degree = DEGREE_1 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E, h_neglect_edge ) call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs ) - end if - end if + endif + endif case ( INTERPOLATION_PQM_IH6IH5 ) if ( n0 >= 6 ) then @@ -243,15 +243,15 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & if (extrapolate) then call PQM_boundary_extrapolation_v1( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect ) - end if + endif else degree = DEGREE_1 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E, h_neglect_edge ) call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs ) - end if - end if + endif + endif end select end subroutine regridding_set_ppolys @@ -288,7 +288,7 @@ subroutine interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefs, & t = target_values(k) x1(k) = get_polynomial_coordinate ( n0, h0, x0, ppoly0_E, ppoly0_coefs, t, degree ) h1(k-1) = x1(k) - x1(k-1) - end do + enddo h1(n1) = x1(n1+1) - x1(n1) end subroutine interpolate_grid @@ -373,7 +373,7 @@ function get_polynomial_coordinate ( N, h, x_g, ppoly_E, ppoly_coefs, & if ( target_value <= ppoly_E(1,1) ) then x_tgt = x_g(1) return ! return because there is no need to look further - end if + endif ! Since discontinuous edge values are allowed, we check whether the target ! value lies between two discontinuous edge values at interior interfaces @@ -383,8 +383,8 @@ function get_polynomial_coordinate ( N, h, x_g, ppoly_E, ppoly_coefs, & x_tgt = x_g(k) return ! return because there is no need to look further exit - end if - end do + endif + enddo ! If the target value is outside the range of all values, we ! force the target coordinate to be equal to the lowest or @@ -392,7 +392,7 @@ function get_polynomial_coordinate ( N, h, x_g, ppoly_E, ppoly_coefs, & if ( target_value >= ppoly_E(N,2) ) then x_tgt = x_g(N+1) return ! return because there is no need to look further - end if + endif ! At this point, we know that the target value is bounded and does not ! lie between discontinuous, monotonic edge values. Therefore, @@ -404,8 +404,8 @@ function get_polynomial_coordinate ( N, h, x_g, ppoly_E, ppoly_coefs, & ( target_value < ppoly_E(k,2) ) ) then k_found = k exit - end if - end do + endif + enddo ! At this point, 'k_found' should be strictly positive. If not, this is ! a major failure because it means we could not find any target cell @@ -419,14 +419,14 @@ function get_polynomial_coordinate ( N, h, x_g, ppoly_E, ppoly_coefs, & 'inconsistent interpolant (perhaps not monotonically '//& 'increasing)' call MOM_error( FATAL, 'Aborting execution' ) - end if + endif ! Reset all polynomial coefficients to 0 and copy those pertaining to ! the found cell a(:) = 0.0 do i = 1,degree+1 a(i) = ppoly_coefs(k_found,i) - end do + enddo ! Guess value to start Newton-Raphson iterations (middle of cell) xi0 = 0.5 @@ -439,7 +439,7 @@ function get_polynomial_coordinate ( N, h, x_g, ppoly_E, ppoly_coefs, & if ( ( iter > NR_ITERATIONS ) .OR. & ( abs(delta) < NR_TOLERANCE ) ) then exit - end if + endif numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + & a(5)*xi0*xi0*xi0*xi0 - target_value @@ -459,16 +459,16 @@ function get_polynomial_coordinate ( N, h, x_g, ppoly_E, ppoly_coefs, & xi0 = 0.0 grad = a(2) if ( grad == 0.0 ) xi0 = xi0 + eps - end if + endif if ( xi0 > 1.0 ) then xi0 = 1.0 grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5) if ( grad == 0.0 ) xi0 = xi0 - eps - end if + endif iter = iter + 1 - end do ! end Newton-Raphson iterations + enddo ! end Newton-Raphson iterations x_tgt = x_g(k_found) + xi0 * h(k_found) end function get_polynomial_coordinate diff --git a/src/ALE/regrid_solvers.F90 b/src/ALE/regrid_solvers.F90 index 7e44039831..18ef1e5e0b 100644 --- a/src/ALE/regrid_solvers.F90 +++ b/src/ALE/regrid_solvers.F90 @@ -63,16 +63,16 @@ subroutine solve_linear_system( A, B, X, system_size ) else ! Go to the next row to see ! if there is a valid pivot there k = k + 1 - end if + endif - end do ! end loop to find pivot + enddo ! end loop to find pivot ! If no pivot could be found, the system is singular and we need ! to end the execution if ( .NOT. found_pivot ) then write(0,*) ' A=',A call MOM_error( FATAL, 'The linear system is singular !' ) - end if + endif ! If the pivot is in a row that is different than row i, that is if ! k is different than i, we need to swap those two rows @@ -81,18 +81,18 @@ subroutine solve_linear_system( A, B, X, system_size ) swap_a = A(i,j) A(i,j) = A(k,j) A(k,j) = swap_a - end do + enddo swap_b = B(i) B(i) = B(k) B(k) = swap_b - end if + endif ! Transform pivot to 1 by dividing the entire row ! (right-hand side included) by the pivot pivot = A(i,i) do j = i,system_size A(i,j) = A(i,j) / pivot - end do + enddo B(i) = B(i) / pivot ! #INV: At this point, A(i,i) is a suitable pivot and it is equal to 1 @@ -103,11 +103,11 @@ subroutine solve_linear_system( A, B, X, system_size ) factor = A(k,i) do j = (i+1),system_size ! j is the column index A(k,j) = A(k,j) - factor * A(i,j) - end do + enddo B(k) = B(k) - factor * B(i) - end do + enddo - end do ! end loop on i + enddo ! end loop on i ! Solve system by back substituting @@ -116,9 +116,9 @@ subroutine solve_linear_system( A, B, X, system_size ) X(i) = B(i) do j = (i+1),system_size X(i) = X(i) - A(i,j) * X(j) - end do + enddo X(i) = X(i) / A(i,i) - end do + enddo end subroutine solve_linear_system @@ -147,18 +147,18 @@ subroutine solve_tridiagonal_system( Al, Ad, Au, B, X, system_size ) do k = 1,N-1 Al(k+1) = Al(k+1) / Ad(k) Ad(k+1) = Ad(k+1) - Al(k+1) * Au(k) - end do + enddo ! Forward sweep do k = 2,N B(k) = B(k) - Al(k) * B(k-1) - end do + enddo ! Backward sweep X(N) = B(N) / Ad(N) do k = N-1,1,-1 X(k) = ( B(k) - Au(k)*X(k+1) ) / Ad(k) - end do + enddo end subroutine solve_tridiagonal_system diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9fca715e42..c1dcf4cf33 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2882,7 +2882,7 @@ subroutine extract_surface_state(CS, sfc_state) endif ! numberOfErrors endif ! localError endif ! mask2dT - enddo; enddo + enddo ; enddo call sum_across_PEs(numberOfErrors) if (numberOfErrors>0) then write(msg(1:240),'(3(a,i9,x))') 'There were a total of ',numberOfErrors, & diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index 9d01f108d1..690fcb42e9 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -549,7 +549,7 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) vhm = 10.0*vhc elseif (abs(vhc) > c1*abs(vhm)) then if (abs(vhc) < c2*abs(vhm)) then ; vhc = (3.0*vhc+(1.0-c2*3.0)*vhm) - else if (abs(vhc) <= c3*abs(vhm)) then ; vhc = vhm + elseif (abs(vhc) <= c3*abs(vhm)) then ; vhc = vhm else ; vhc = slope*vhc+(1.0-c3*slope)*vhm endif endif diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index c430179917..a54d7bb01f 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -665,7 +665,7 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, LB, vol_CFL, & if (marginal) then ; h_u(I,j,k) = h_marg else ; h_u(I,j,k) = h_avg ; endif - enddo; enddo ; enddo + enddo ; enddo ; enddo if (present(visc_rem_u)) then !$OMP parallel do default(shared) do k=1,nz ; do j=jsh,jeh ; do I=ish-1,ieh @@ -1948,7 +1948,7 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dMx, dMn)) ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j)) endif - enddo; enddo + enddo ; enddo if (local_open_BC) then do n=1, OBC%number_of_segments @@ -1975,7 +1975,7 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ ! Left/right values following Eq. B2 in Lin 1994, MWR (132) h_L(i,j) = 0.5*( h_im1 + h_in(i,j) ) + oneSixth*( slp(i-1,j) - slp(i,j) ) h_R(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + oneSixth*( slp(i,j) - slp(i+1,j) ) - enddo; enddo + enddo ; enddo endif if (local_open_BC) then diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 91f9f6546b..3668861db7 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -763,7 +763,7 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg, PF) if (Je_obc>Js_obc) then OBC%segment(l_seg)%direction = OBC_DIRECTION_E - else if (Je_obcIs_obc) then OBC%segment(l_seg)%direction = OBC_DIRECTION_S - else if (Ie_obc0.) then + elseif (G%mask2dCu(I,j)>0.) then h_stack(:) = h(i+ishift,j,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & G%ke, h_stack, segment%field(m)%buffer_dst(I,J,:)) - else if (G%mask2dCu(I,j+1)>0.) then + elseif (G%mask2dCu(I,j+1)>0.) then h_stack(:) = h(i+ishift,j+1,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src,segment%field(m)%dz_src(I,j,:), & @@ -2462,10 +2462,10 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) if (segment%field(m)%name == 'V') then allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,G%ke)) allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc:je_obc)) - else if (segment%field(m)%name == 'U') then + elseif (segment%field(m)%name == 'U') then allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,G%ke)) allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc)) - else if (segment%field(m)%name == 'DVDX') then + elseif (segment%field(m)%name == 'DVDX') then allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,G%ke)) else allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,G%ke)) @@ -2474,10 +2474,10 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) if (segment%field(m)%name == 'U') then allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,G%ke)) allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc:je_obc)) - else if (segment%field(m)%name == 'V') then + elseif (segment%field(m)%name == 'V') then allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,G%ke)) allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc)) - else if (segment%field(m)%name == 'DUDY') then + elseif (segment%field(m)%name == 'DUDY') then allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,G%ke)) else allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,G%ke)) @@ -2583,12 +2583,12 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) if (associated(segment%field(m)%buffer_dst)) then do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc segment%tr_Reg%Tr(1)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k) - enddo; enddo; enddo + enddo ; enddo ; enddo if (.not. segment%tr_Reg%Tr(1)%is_initialized) then ! if the tracer reservoir has not yet been initialized, then set to external value. do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc segment%tr_Reg%Tr(1)%tres(i,j,k) = segment%tr_Reg%Tr(1)%t(i,j,k) - enddo; enddo; enddo + enddo ; enddo ; enddo segment%tr_Reg%Tr(1)%is_initialized=.true. endif else @@ -2598,12 +2598,12 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) if (associated(segment%field(m)%buffer_dst)) then do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc segment%tr_Reg%Tr(2)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k) - enddo; enddo; enddo + enddo ; enddo ; enddo if (.not. segment%tr_Reg%Tr(1)%is_initialized) then !if the tracer reservoir has not yet been initialized, then set to external value. do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc segment%tr_Reg%Tr(2)%tres(i,j,k) = segment%tr_Reg%Tr(2)%t(i,j,k) - enddo; enddo; enddo + enddo ; enddo ; enddo segment%tr_Reg%Tr(1)%is_initialized=.true. endif else diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 7add057e0e..943ed8eb83 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -386,7 +386,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if ((CS%id_Tpot > 0) .or. (CS%id_tob > 0)) then do k=1,nz ; do j=js,je ; do i=is,ie work_3d(i,j,k) = gsw_pt_from_ct(tv%S(i,j,k),tv%T(i,j,k)) - enddo; enddo ; enddo + enddo ; enddo ; enddo if (CS%id_Tpot > 0) call post_data(CS%id_Tpot, work_3d, CS%diag) if (CS%id_tob > 0) call post_data(CS%id_tob, work_3d(:,:,nz), CS%diag, mask=G%mask2dT) endif @@ -403,7 +403,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if ((CS%id_Sprac > 0) .or. (CS%id_sob > 0)) then do k=1,nz ; do j=js,je ; do i=is,ie work_3d(i,j,k) = gsw_sp_from_sr(tv%S(i,j,k)) - enddo; enddo ; enddo + enddo ; enddo ; enddo if (CS%id_Sprac > 0) call post_data(CS%id_Sprac, work_3d, CS%diag) if (CS%id_sob > 0) call post_data(CS%id_sob, work_3d(:,:,nz), CS%diag, mask=G%mask2dT) endif @@ -718,19 +718,19 @@ subroutine find_weights(Rlist, R_in, k, nz, wt, wt_p) if ((k_lower == 1) .or. (R_in >= Rlist(k_lower))) exit k_upper = k_lower inc = inc*2 - end do + enddo else do k_upper = min(k_upper+inc, nz) if ((k_upper == nz) .or. (R_in < Rlist(k_upper))) exit k_lower = k_upper inc = inc*2 - end do + enddo endif if ((k_lower == 1) .and. (R_in <= Rlist(k_lower))) then k = 1 ; wt = 1.0 ; wt_p = 0.0 - else if ((k_upper == nz) .and. (R_in >= Rlist(k_upper))) then + elseif ((k_upper == nz) .and. (R_in >= Rlist(k_upper))) then k = nz-1 ; wt = 0.0 ; wt_p = 1.0 else do @@ -741,7 +741,7 @@ subroutine find_weights(Rlist, R_in, k, nz, wt, wt_p) else k_lower = k_new endif - end do + enddo ! Uncomment this as a code check ! if ((R_in < Rlist(k_lower)) .or. (R_in >= Rlist(k_upper)) .or. (k_upper-k_lower /= 1)) & diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 4db4d30c18..17427fb80f 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -234,7 +234,7 @@ subroutine MOM_sum_output_init(G, param_file, directory, ntrnc, & call get_filename_appendix(filename_appendix) if (len_trim(filename_appendix) > 0) then energyfile = trim(energyfile) //'.'//trim(filename_appendix) - end if + endif CS%energyfile = trim(slasher(directory))//trim(energyfile) call log_param(param_file, mdl, "output_path/ENERGYFILE", CS%energyfile) @@ -606,11 +606,11 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, CS, tracer_CSp, OBC, dt_forc else if ((CS%timeunit >= 0.99) .and. (CS%timeunit < 1.01)) then time_units = " [seconds] " - else if ((CS%timeunit >= 3599.0) .and. (CS%timeunit < 3601.0)) then + elseif ((CS%timeunit >= 3599.0) .and. (CS%timeunit < 3601.0)) then time_units = " [hours] " - else if ((CS%timeunit >= 86399.0) .and. (CS%timeunit < 86401.0)) then + elseif ((CS%timeunit >= 86399.0) .and. (CS%timeunit < 86401.0)) then time_units = " [days] " - else if ((CS%timeunit >= 3.0e7) .and. (CS%timeunit < 3.2e7)) then + elseif ((CS%timeunit >= 3.0e7) .and. (CS%timeunit < 3.2e7)) then time_units = " [years] " else write(time_units,'(9x,"[",es8.2," s] ")') CS%timeunit diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index ec4a78fc7b..9244b33738 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -120,7 +120,7 @@ subroutine wave_speed(h, tv, G, GV, cg1, CS, full_halos, use_ebt_mode, & if (calc_modal_structure) then do k=1,nz; do j=js,je; do i=is,ie modal_structure(i,j,k) = 0.0 - enddo; enddo; enddo + enddo ; enddo ; enddo endif S => tv%S ; T => tv%T diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index dceed058f2..2df645c338 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -2341,7 +2341,7 @@ subroutine convert_temp_salt_for_TEOS10(T, S, press, G, kd, mask_z, EOS) ! S(i,j,k) = gsw_sa_from_sp(S(i,j,k),p,G%geoLonT(i,j),G%geoLatT(i,j)) T(i,j,k) = gsw_ct_from_pt(S(i,j,k),T(i,j,k)) endif - enddo; enddo; enddo + enddo ; enddo ; enddo end subroutine convert_temp_salt_for_TEOS10 ! Extractor routine for the EOS type if the members need to be accessed outside this module diff --git a/src/framework/MOM_checksums.F90 b/src/framework/MOM_checksums.F90 index fd880f6656..da90ef1ad7 100644 --- a/src/framework/MOM_checksums.F90 +++ b/src/framework/MOM_checksums.F90 @@ -198,7 +198,7 @@ integer function subchk(array, HI, di, dj, scale) do j=HI%jsc+dj,HI%jec+dj; do i=HI%isc+di,HI%iec+di bc = bitcount(abs(scale*array(i,j))) subchk = subchk + bc - enddo; enddo + enddo ; enddo call sum_across_PEs(subchk) subchk=mod(subchk,1000000000) end function subchk @@ -385,7 +385,7 @@ integer function subchk(array, HI, di, dj, scale) do J=HI%jsc+dj,HI%jec+dj; do I=HI%isc+di,HI%iec+di bc = bitcount(abs(scale*array(I,J))) subchk = subchk + bc - enddo; enddo + enddo ; enddo call sum_across_PEs(subchk) subchk=mod(subchk,1000000000) end function subchk @@ -573,7 +573,7 @@ integer function subchk(array, HI, di, dj, scale) do j=HI%jsc+dj,HI%jec+dj; do I=HI%isc+di,HI%iec+di bc = bitcount(abs(scale*array(I,j))) subchk = subchk + bc - enddo; enddo + enddo ; enddo call sum_across_PEs(subchk) subchk=mod(subchk,1000000000) end function subchk @@ -718,7 +718,7 @@ integer function subchk(array, HI, di, dj, scale) do J=HI%jsc+dj,HI%jec+dj; do i=HI%isc+di,HI%iec+di bc = bitcount(abs(scale*array(i,J))) subchk = subchk + bc - enddo; enddo + enddo ; enddo call sum_across_PEs(subchk) subchk=mod(subchk,1000000000) end function subchk diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 176e6e6d13..d4f8dbff57 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -225,7 +225,7 @@ subroutine fill_miss_2d(aout,good,fill,prev,G,smooth,num_pass,relc,crit,keep_bug endif enddo enddo - else if (nfill == nfill_prev) then + elseif (nfill == nfill_prev) then print *,& 'Unable to fill missing points using either data at the same vertical level from a connected basin'//& 'or using a point from a previous vertical level. Make sure that the original data has some valid'//& @@ -236,7 +236,7 @@ subroutine fill_miss_2d(aout,good,fill,prev,G,smooth,num_pass,relc,crit,keep_bug nfill = sum(fill_pts(is:ie,js:je)) call sum_across_PEs(nfill) - end do + enddo if (do_smooth) then do k=1,npass @@ -1010,7 +1010,7 @@ subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) zi(:,:)=mp(1:ni,1:nj) mp = fill_boundaries(zi,cyclic_x,tripolar_n) -end do +enddo diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 079ac6ba3a..17ba449715 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -216,7 +216,7 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit call MOM_error(WARNING, "MOM_io create_file: "//trim(vars(k)%name)//& " has unrecognized t_grid "//trim(vars(k)%t_grid)) end select - end do + enddo if ((use_lath .or. use_lonh .or. use_latq .or. use_lonq)) then if (.not.domain_set) call MOM_error(FATAL, "create_file: "//& @@ -260,13 +260,13 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit ! Set appropriate units, depending on the value. if (timeunit < 0.0) then time_units = "days" ! The default value. - else if ((timeunit >= 0.99) .and. (timeunit < 1.01)) then + elseif ((timeunit >= 0.99) .and. (timeunit < 1.01)) then time_units = "seconds" - else if ((timeunit >= 3599.0) .and. (timeunit < 3601.0)) then + elseif ((timeunit >= 3599.0) .and. (timeunit < 3601.0)) then time_units = "hours" - else if ((timeunit >= 86399.0) .and. (timeunit < 86401.0)) then + elseif ((timeunit >= 86399.0) .and. (timeunit < 86401.0)) then time_units = "days" - else if ((timeunit >= 3.0e7) .and. (timeunit < 3.2e7)) then + elseif ((timeunit >= 3.0e7) .and. (timeunit < 3.2e7)) then time_units = "years" else write(time_units,'(es8.2," s")') timeunit diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 6944647008..bb34fe4985 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -869,7 +869,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV) seconds = seconds + 60*minute + 3600*hour if (year <= 9999) then write(restartname,'("_Y",I4.4,"_D",I3.3,"_S",I5.5)') year, days, seconds - else if (year <= 99999) then + elseif (year <= 99999) then write(restartname,'("_Y",I5.5,"_D",I3.3,"_S",I5.5)') year, days, seconds else write(restartname,'("_Y",I10.10,"_D",I3.3,"_S",I5.5)') year, days, seconds @@ -920,8 +920,8 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV) restartname = restartname(1:length-3)//'.'//trim(filename_appendix)//'.nc' else restartname = restartname(1:length) //'.'//trim(filename_appendix) - end if - end if + endif + endif restartpath = trim(directory)// trim(restartname) @@ -1453,8 +1453,8 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & restartname = restartname(1:length-3)//'.'//trim(filename_appendix)//'.nc' else restartname = restartname(1:length) //'.'//trim(filename_appendix) - end if - end if + endif + endif filepath = trim(directory) // trim(restartname) if (num_restart < 10) then diff --git a/src/framework/MOM_string_functions.F90 b/src/framework/MOM_string_functions.F90 index f56834a8f6..142134ddc5 100644 --- a/src/framework/MOM_string_functions.F90 +++ b/src/framework/MOM_string_functions.F90 @@ -42,7 +42,7 @@ function lowercase(input_string) do k=1, len_trim(input_string) if (lowercase(k:k) >= 'A' .and. lowercase(k:k) <= 'Z') & lowercase(k:k) = achar(ichar(lowercase(k:k))+co) - end do + enddo end function lowercase function uppercase(input_string) @@ -59,7 +59,7 @@ function uppercase(input_string) do k=1, len_trim(input_string) if (uppercase(k:k) >= 'a' .and. uppercase(k:k) <= 'z') & uppercase(k:k) = achar(ichar(uppercase(k:k))+co) - end do + enddo end function uppercase function left_int(i) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 94b628221c..6be4f7d0d3 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1089,7 +1089,7 @@ subroutine add_shelf_flux(G, CS, state, forces, fluxes) if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then sponge_area = sponge_area + G%areaT(i,j) endif - enddo; enddo + enddo ; enddo ! take into account changes in mass (or thickness) when imposing ice shelf mass if (CS%shelf_mass_is_dynamic .and. CS%override_shelf_movement .and. & @@ -1125,7 +1125,7 @@ subroutine add_shelf_flux(G, CS, state, forces, fluxes) shelf_mass1 = shelf_mass1 + (CS%mass_shelf(i,j) * CS%area_shelf_h(i,j)) endif - enddo; enddo + enddo ; enddo call mpp_sum(shelf_mass0); call mpp_sum(shelf_mass1) delta_mass_shelf = (shelf_mass1 - shelf_mass0)/CS%time_step ! delta_mass_shelf = (shelf_mass1 - shelf_mass0)* & @@ -1152,7 +1152,7 @@ subroutine add_shelf_flux(G, CS, state, forces, fluxes) fluxes%sens(i,j) = fluxes%vprec(i,j) * CS%Cp * CS%T0 ! W /m^2 fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * CS%S0*1.0e-3 ! kg (salt)/(m^2 s) endif - enddo; enddo + enddo ; enddo if (CS%DEBUG) then if (is_root_pe()) write(*,*)'Mean melt flux (kg/(m^2 s)),dt',mean_melt_flux,CS%time_step @@ -1679,12 +1679,12 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl endif - ! else if (CS%shelf_mass_is_dynamic) then + ! elseif (CS%shelf_mass_is_dynamic) then ! call initialize_ice_shelf_boundary ( CS%u_face_mask_boundary, CS%v_face_mask_boundary, & ! CS%u_flux_boundary_values, CS%v_flux_boundary_values, & ! CS%u_boundary_values, CS%v_boundary_values, CS%h_boundary_values, & ! CS%hmask, G, param_file) - end if + endif if (CS%shelf_mass_is_dynamic .and. .not. CS%override_shelf_movement) then ! the only reason to initialize boundary conds is if the shelf is dynamic @@ -1694,7 +1694,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl !MJH CS%u_boundary_values, CS%v_boundary_values, CS%h_boundary_values, & !MJH CS%hmask, G, param_file) - end if + endif if (new_sim .and. (.not. (CS%override_shelf_movement .and. CS%mass_from_file))) then @@ -5856,14 +5856,14 @@ subroutine savearray2(fname,A,flag) WRITE(fin) 'SECOND DIMENSION TOO LARGE' CLOSE(fin) RETURN -END IF +ENDIF DO i=1,M WRITE(ln,'(E17.9)') A(i,1) DO j=2,N WRITE(sing,'(E17.9)') A(i,j) ln = TRIM(ln) // ' ' // TRIM(sing) - END DO + ENDDO if (i == 1) THEN @@ -5889,14 +5889,14 @@ subroutine savearray2(fname,A,flag) FMT1 = TRIM(FMT1) // ')' - END IF + ENDIF WRITE(UNIT=fin,IOSTAT=iock,FMT=TRIM(FMT1)) TRIM(ln) if (iock /= 0) THEN - PRINT*,iock - END IF -END DO + PRINT *,iock + ENDIF +ENDDO CLOSE(FIN) diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index ba84d55763..f59728af8f 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -538,7 +538,7 @@ subroutine set_grid_metrics_cartesian(G, param_file) if (units_temp(1:1) == 'k') then ! Axes are measured in km. dx_everywhere = 1000.0 * G%len_lon / (REAL(niglobal)) dy_everywhere = 1000.0 * G%len_lat / (REAL(njglobal)) - else if (units_temp(1:1) == 'm') then ! Axes are measured in m. + elseif (units_temp(1:1) == 'm') then ! Axes are measured in m. dx_everywhere = G%len_lon / (REAL(niglobal)) dy_everywhere = G%len_lat / (REAL(njglobal)) else ! Axes are measured in degrees of latitude and longitude. @@ -679,7 +679,7 @@ subroutine set_grid_metrics_spherical(G, param_file) ! G%dxBu(I,J) = G%Rad_Earth * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 ) G%dyBu(I,J) = G%Rad_Earth * dLat*PI_180 G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) - enddo; enddo + enddo ; enddo do J=JsdB,JedB ; do i=isd,ied G%geoLonCv(i,J) = grid_LonT(i) @@ -690,7 +690,7 @@ subroutine set_grid_metrics_spherical(G, param_file) G%dxCv(i,J) = G%Rad_Earth * COS( G%geoLatCv(i,J)*PI_180 ) * dL_di ! G%dxCv(i,J) = G%Rad_Earth * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 ) G%dyCv(i,J) = G%Rad_Earth * dLat*PI_180 - enddo; enddo + enddo ; enddo do j=jsd,jed ; do I=IsdB,IedB G%geoLonCu(I,j) = grid_lonB(I) @@ -701,7 +701,7 @@ subroutine set_grid_metrics_spherical(G, param_file) G%dxCu(I,j) = G%Rad_Earth * COS( G%geoLatCu(I,j)*PI_180 ) * dL_di ! G%dxCu(I,j) = G%Rad_Earth * dLon*PI_180 * COS( latitude ) G%dyCu(I,j) = G%Rad_Earth * dLat*PI_180 - enddo; enddo + enddo ; enddo do j=jsd,jed ; do i=isd,ied G%geoLonT(i,j) = grid_LonT(i) @@ -717,7 +717,7 @@ subroutine set_grid_metrics_spherical(G, param_file) ! dL_di = G%geoLatCv(i,max(jsd,J-1))*PI_180 ! In radians ! G%areaT(i,j) = Rad_Earth**2*dLon*dLat*ABS(SIN(latitude)-SIN(dL_di)) G%areaT(i,j) = G%dxT(i,j) * G%dyT(i,j) - enddo; enddo + enddo ; enddo call callTree_leave("set_grid_metrics_spherical()") end subroutine set_grid_metrics_spherical @@ -1207,7 +1207,7 @@ function Int_dj_dy(y, GP) if (y >= y_eq_enhance) then r = r + I_C0*0.5*(GP%lat_enhance_factor - 1.0)*y_eq_enhance - else if (y <= -y_eq_enhance) then + elseif (y <= -y_eq_enhance) then r = r - I_C0*0.5*(GP%lat_enhance_factor - 1.0)*y_eq_enhance else r = r + I_C0*0.5*(GP%lat_enhance_factor - 1.0) * & diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index b150b8c4ad..e818c33acd 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -1187,9 +1187,9 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file) call create_file(unit, trim(filepath), vars, nFlds_used, fields, & file_threading, dG=G) - do J=Jsq,Jeq; do I=Isq,Ieq; out_q(I,J) = G%geoLatBu(I,J); enddo; enddo + do J=Jsq,Jeq; do I=Isq,Ieq; out_q(I,J) = G%geoLatBu(I,J); enddo ; enddo call write_field(unit, fields(1), G%Domain%mpp_domain, out_q) - do J=Jsq,Jeq; do I=Isq,Ieq; out_q(I,J) = G%geoLonBu(I,J); enddo; enddo + do J=Jsq,Jeq; do I=Isq,Ieq; out_q(I,J) = G%geoLonBu(I,J); enddo ; enddo call write_field(unit, fields(2), G%Domain%mpp_domain, out_q) call write_field(unit, fields(3), G%Domain%mpp_domain, G%geoLatT) call write_field(unit, fields(4), G%Domain%mpp_domain, G%geoLonT) @@ -1210,7 +1210,7 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file) do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = G%dyCv(i,J) ; enddo ; enddo call write_field(unit, fields(10), G%Domain%mpp_domain, out_v) - do j=js,je ; do i=is,ie ; out_h(i,j) = G%dxT(i,j); enddo; enddo + do j=js,je ; do i=is,ie ; out_h(i,j) = G%dxT(i,j); enddo ; enddo call write_field(unit, fields(11), G%Domain%mpp_domain, out_h) do j=js,je ; do i=is,ie ; out_h(i,j) = G%dyT(i,j) ; enddo ; enddo call write_field(unit, fields(12), G%Domain%mpp_domain, out_h) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 49153586b7..491c806a6b 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1709,18 +1709,18 @@ subroutine initialize_temp_salt_linear(T, S, G, param_file, just_read_params) ! S(:,:,1) = S_top ! do k = 2,G%ke ! S(:,:,k) = S(:,:,k-1) + delta_S -! end do +! enddo do k = 1,G%ke S(:,:,k) = S_top - S_range*((real(k)-0.5)/real(G%ke)) T(:,:,k) = T_top - T_range*((real(k)-0.5)/real(G%ke)) - end do + enddo ! ! Prescribe temperature ! delta_T = T_range / ( G%ke - 1.0 ) ! T(:,:,1) = T_top ! do k = 2,G%ke ! T(:,:,k) = T(:,:,k-1) + delta_T -! end do +! enddo ! delta = 1 ! T(:,:,G%ke/2 - (delta-1):G%ke/2 + delta) = 1.0 @@ -1848,7 +1848,7 @@ subroutine initialize_sponges_file(G, GV, use_temperature, tv, param_file, CSp, ! apply the sponges, along with the interface heights. ! call initialize_sponge(Idamp, eta, G, param_file, CSp) deallocate(eta) - else if (.not. new_sponges) then ! ALE mode + elseif (.not. new_sponges) then ! ALE mode call field_size(filename,eta_var,siz,no_domain=.true.) if (siz(1) /= G%ieg-G%isg+1 .or. siz(2) /= G%jeg-G%jsg+1) & @@ -1872,7 +1872,7 @@ subroutine initialize_sponges_file(G, GV, use_temperature, tv, param_file, CSp, enddo ; enddo ; enddo do k=1,nz; do j=js,je ; do i=is,ie h(i,j,k) = eta(i,j,k)-eta(i,j,k+1) - enddo ; enddo; enddo + enddo ; enddo ; enddo call initialize_ALE_sponge(Idamp, G, param_file, ALE_CSp, h, nz_data) deallocate(eta) deallocate(h) @@ -1910,7 +1910,7 @@ subroutine initialize_sponges_file(G, GV, use_temperature, tv, param_file, CSp, call set_up_sponge_field(tmp, tv%T, G, nz, CSp) call MOM_read_data(filename, salin_var, tmp(:,:,:), G%Domain) call set_up_sponge_field(tmp, tv%S, G, nz, CSp) - else if (use_temperature) then + elseif (use_temperature) then call set_up_ALE_sponge_field(filename, potemp_var, Time, G, tv%T, ALE_CSp) call set_up_ALE_sponge_field(filename, salin_var, Time, G, tv%S, ALE_CSp) endif diff --git a/src/initialization/midas_vertmap.F90 b/src/initialization/midas_vertmap.F90 index 373062ffc3..8d022d97cc 100644 --- a/src/initialization/midas_vertmap.F90 +++ b/src/initialization/midas_vertmap.F90 @@ -408,7 +408,7 @@ function bisect_fast(a, x, lo, hi) result(bi_r) if (PRESENT(lo)) then where (lo>0) lo_=lo -end if +endif if (PRESENT(hi)) then where (hi>0) hi_=hi endif @@ -950,7 +950,7 @@ subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) zi(:,:)=mp(1:ni,1:nj) mp = fill_boundaries(zi,cyclic_x,tripolar_n) -end do +enddo return diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 2672308fd7..de5a97363b 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -308,15 +308,15 @@ subroutine init_oda(Time, G, GV, CS) do i=1, CS%ni; do j=1, CS%nj if ( global2D(i,j) > 1 ) then T_grid%mask(i,j,k) = 1.0 - end if - end do; end do + endif + enddo ; enddo if (k == 1) then T_grid%z(:,:,k) = global2D/2 else T_grid%z(:,:,k) = T_grid%z(:,:,k-1) + (global2D + global2D_old)/2 - end if + endif global2D_old = global2D - end do + enddo call ocean_da_core_init(CS%mpp_domain, T_grid, CS%Profiles, Time) @@ -363,7 +363,7 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) CS%nk, CS%h(i,j,:), T(i,j,:)) call remapping_core_h(CS%remapCS, GV%ke, h(i,j,:), tv%S(i,j,:), & CS%nk, CS%h(i,j,:), S(i,j,:)) - enddo; enddo + enddo ; enddo do m=1,CS%ensemble_size call mpp_redistribute(CS%domains(m)%mpp_domain, T,& @@ -449,7 +449,7 @@ subroutine get_posterior_tracer(Time, CS, G, GV, h, tv, increment) used=send_data(CS%Ocean_posterior%id_s(m), CS%Ocean_posterior%S(isc:iec,jsc:jec,:,m), CS%Time) endif endif - end do + enddo tv => CS%tv h => CS%h @@ -479,7 +479,7 @@ subroutine oda(Time, CS) !! switch back to ensemble member pelist call set_current_pelist(CS%ensemble_pelist(CS%ensemble_id,:)) - end if + endif return end subroutine oda diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index c4e771375c..11798d3bdb 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -433,7 +433,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, else do j=js-2,je+2 ; do I=Isq-1,Ieq+1 h_u(I,j) = 0.5 * (h(i,j,k) + h(i+1,j,k)) - enddo; enddo + enddo ; enddo do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2 h_v(i,J) = 0.5 * (h(i,j,k) + h(i,j+1,k)) enddo ; enddo diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 076dab7b56..3be1ae6192 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -284,7 +284,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & do j=jsd,jed ; do i=isd,ied ; test(i,j,1) = 1.0 ; test(i,j,2) = 0.0 ; enddo ; enddo do m=1,CS%nMode ; do fr=1,CS%nFreq call create_group_pass(pass_En, CS%En(:,:,:,fr,m), G%domain) - enddo; enddo + enddo ; enddo call create_group_pass(pass_test, test(:,:,1), test(:,:,2), G%domain, stagger=AGRID) call start_group_pass(pass_test, G%domain) @@ -920,7 +920,7 @@ subroutine refract(En, cn, freq, dt, G, NAngle, use_PPMang) call MOM_error(WARNING, "refract: CFL exceeds 1.", .true.) if (CFL_ang(i,j,A) > 0.0) then ; CFL_ang(i,j,A) = 1.0 ; else ; CFL_ang(i,j,A) = -1.0 ; endif endif - enddo; enddo + enddo ; enddo ! Advect in angular space if (.not.use_PPMang) then @@ -931,7 +931,7 @@ subroutine refract(En, cn, freq, dt, G, NAngle, use_PPMang) else Flux_E(i,A) = CFL_ang(i,j,A) * En2d(i,A+1) endif - enddo; enddo + enddo ; enddo else ! Use PPM do i=is,ie @@ -1109,7 +1109,7 @@ subroutine propagate(En, cn, freq, dt, G, CS, NAngle) ! Apply the propagation WITH CORNER ADVECTION/FINITE VOLUME APPROACH. LB%jsh = js ; LB%jeh = je ; LB%ish = is ; LB%ieh = ie call propagate_corner_spread(En(:,:,a), a, NAngle, speed, dt, G, CS, LB) - end do ! a-loop + enddo ! a-loop else ! IMPLEMENT PPM ADVECTION IN HORIZONTAL----------------------- ! These could be in the control structure, as they do not vary. @@ -1436,7 +1436,7 @@ subroutine propagate_corner_spread(En, energized_wedge, NAngle, speed, dt, G, CS enddo ! m-loop ! update energy in cell En(i,j) = sum(E_new)/Nsubrays - enddo; enddo + enddo ; enddo end subroutine propagate_corner_spread ! #@# This subroutine needs a doxygen description @@ -2069,7 +2069,7 @@ subroutine PPM_reconstruction_x(h_in, h_l, h_r, G, LB, simple_2nd) slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dMx, dMn)) ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j)) endif - enddo; enddo + enddo ; enddo do j=jsl,jel ; do i=isl,iel ! Neighboring values should take into account any boundaries. The 3 @@ -2081,7 +2081,7 @@ subroutine PPM_reconstruction_x(h_in, h_l, h_r, G, LB, simple_2nd) ! Left/right values following Eq. B2 in Lin 1994, MWR (132) h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + oneSixth*( slp(i-1,j) - slp(i,j) ) h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + oneSixth*( slp(i,j) - slp(i+1,j) ) - enddo; enddo + enddo ; enddo endif call PPM_limit_pos(h_in, h_l, h_r, 0.0, G, isl, iel, jsl, jel) @@ -2502,7 +2502,7 @@ subroutine internal_tides_init(Time, G, GV, param_file, diag, CS) ! will be multiplied by N and En to get into [W m-2] CS%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*GV%Rho0*& kappa_itides * h2(i,j) - enddo; enddo + enddo ; enddo ! Read in prescribed coast/ridge/shelf angles from file call get_param(param_file, mdl, "REFL_ANGLE_FILE", refl_angle_file, & @@ -2556,7 +2556,7 @@ subroutine internal_tides_init(Time, G, GV, param_file, diag, CS) do i=isd,ied; do j=jsd,jed if (ridge_temp(i,j) == 1) then; CS%refl_dbl(i,j) = .true. else ; CS%refl_dbl(i,j) = .false. ; endif - enddo; enddo + enddo ; enddo ! Read in prescribed land mask from file (if overwriting -BDM). ! This should be done in MOM_initialize_topography subroutine diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 308b7ca9b6..93aeb6f750 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -198,7 +198,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ allocate(CS%Ref_h%p(CS%nz_data,CS%num_col)) do col=1,CS%num_col ; do K=1,CS%nz_data CS%Ref_h%p(K,col) = data_h(CS%col_i(col),CS%col_j(col),K) - enddo; enddo + enddo ; enddo endif total_sponge_cols = CS%num_col @@ -224,7 +224,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & CS%num_col_u = CS%num_col_u + 1 - enddo; enddo + enddo ; enddo if (CS%num_col_u > 0) then @@ -247,7 +247,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ allocate(CS%Ref_hu%p(CS%nz_data,CS%num_col_u)) do col=1,CS%num_col_u ; do K=1,CS%nz_data CS%Ref_hu%p(K,col) = data_hu(CS%col_i_u(col),CS%col_j_u(col),K) - enddo; enddo + enddo ; enddo endif total_sponge_cols_u = CS%num_col_u call sum_across_PEs(total_sponge_cols_u) @@ -261,7 +261,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ Iresttime_v(i,J) = 0.5 * (Iresttime(i,j) + Iresttime(i,j+1)) if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) & CS%num_col_v = CS%num_col_v + 1 - enddo; enddo + enddo ; enddo if (CS%num_col_v > 0) then @@ -400,7 +400,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & CS%num_col_u = CS%num_col_u + 1 - enddo; enddo + enddo ; enddo if (CS%num_col_u > 0) then @@ -432,7 +432,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) Iresttime_v(i,J) = 0.5 * (Iresttime(i,j) + Iresttime(i,j+1)) if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) & CS%num_col_v = CS%num_col_v + 1 - enddo; enddo + enddo ; enddo if (CS%num_col_v > 0) then @@ -863,7 +863,7 @@ subroutine apply_ALE_sponge(h, dt, G, CS, Time) ! u points do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB; do k=1,nz hu(I,j,k) = 0.5 * (h(i,j,k) + h(i+1,j,k)) - enddo; enddo; enddo + enddo ; enddo ; enddo if (CS%new_sponges) then if (.not. present(Time)) & @@ -935,7 +935,7 @@ subroutine apply_ALE_sponge(h, dt, G, CS, Time) ! v points do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec; do k=1,nz hv(i,J,k) = 0.5 * (h(i,j,k) + h(i,j+1,k)) - enddo; enddo; enddo + enddo ; enddo ; enddo do c=1,CS%num_col_v i = CS%col_i_v(c) ; j = CS%col_j_v(c) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index c0289bbd79..2861218128 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -401,7 +401,7 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) 'Constant value to enhance VT2 in KPP.', & default=1.0) endif - end if + endif call closeParameterBlock(paramFile) call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 201588a2c2..7eafb011bd 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -1343,7 +1343,7 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & h_prev = h_ent ; h_ent = h_prev+dh_Newt if (h_ent > h_max) then h_ent = 0.5*(h_prev+h_max) - else if (h_ent < h_min) then + elseif (h_ent < h_min) then h_ent = 0.5*(h_prev+h_min) endif diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 index afdebe4ae5..2678b18e1a 100644 --- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 +++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 @@ -1091,7 +1091,7 @@ subroutine find_PE_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, & ColHt_chg = ColHt_core * y1 if (ColHt_chg < 0.0) PE_chg = PE_chg - pres * ColHt_chg if (present(ColHt_cor)) ColHt_cor = -pres * min(ColHt_chg, 0.0) - else if (present(ColHt_cor)) then + elseif (present(ColHt_cor)) then y1 = dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) ColHt_cor = -pres * min(ColHt_core * y1, 0.0) endif diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 89a11217fa..a58773d066 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -1710,7 +1710,7 @@ subroutine find_PE_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, & ColHt_chg = ColHt_core * y1 if (ColHt_chg < 0.0) PE_chg = PE_chg - pres * ColHt_chg if (present(ColHt_cor)) ColHt_cor = -pres * min(ColHt_chg, 0.0) - else if (present(ColHt_cor)) then + elseif (present(ColHt_cor)) then y1 = dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) ColHt_cor = -pres * min(ColHt_core * y1, 0.0) endif diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 index c9f10826db..5f3f982dd1 100644 --- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 +++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 @@ -736,7 +736,7 @@ subroutine entrainment_diffusive(u, v, h, tv, fluxes, dt, G, GV, CS, ea, eb, & F(i,k) = MIN(F(i,k), ds_dsp1(i,k)*( ((F(i,k-1) + & dsp1_ds(i,k-1)*F(i,k-1)) - F(i,k-2)) + (h(i,j,k-1) - Angstrom))) F(i,k) = MAX(F(i,k),MIN(minF(i,k),0.0)) - else if (k == kb(i)+1) then + elseif (k == kb(i)+1) then F(i,k) = MIN(F(i,k), ds_dsp1(i,k)*( ((F(i,k-1) + eakb(i)) - & eb_kmb(i)) + (h(i,j,k-1) - Angstrom))) F(i,k) = MAX(F(i,k),MIN(minF(i,k),0.0)) @@ -791,7 +791,7 @@ subroutine entrainment_diffusive(u, v, h, tv, fluxes, dt, G, GV, CS, ea, eb, & ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*F_cor eb(i,j,k) = eb(i,j,k) + F_cor - else if ((k==kb(i)) .and. (F(i,k) > 0.0)) then + elseif ((k==kb(i)) .and. (F(i,k) > 0.0)) then ! Rho_cor is the density anomaly that needs to be corrected, ! taking into account that the true potential density of the ! deepest buffer layer is not exactly what is returned as dS_kb. @@ -817,7 +817,7 @@ subroutine entrainment_diffusive(u, v, h, tv, fluxes, dt, G, GV, CS, ea, eb, & ea(i,j,k) = ea(i,j,k) + ea_cor eb(i,j,k) = eb(i,j,k) - (dS_kb(i) * I_dSkbp1(i)) * ea_cor - else if (k < kb(i)) then + elseif (k < kb(i)) then ! Repetative, unless ea(kb) has been corrected. ea(i,j,k) = ea(i,j,k+1) endif @@ -1007,7 +1007,7 @@ subroutine F_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb, ! elsewhere, so F should always be nonnegative. ea(i,j,k) = dsp1_ds(i,k)*F(i,k) eb(i,j,k) = F(i,k) - else if (k == kb(i)) then + elseif (k == kb(i)) then ea(i,j,k) = eakb(i) eb(i,j,k) = F(i,k) elseif (k == kb(i)-1) then diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index e65af9183c..3c9188b6bb 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -372,7 +372,7 @@ subroutine int_tide_input_init(Time, G, GV, param_file, diag, CS, itide) ! Compute the fixed part of internal tidal forcing; units are [J m-2] here. CS%TKE_itidal_coef(i,j) = 0.5*kappa_h2_factor*GV%Rho0*& kappa_itides * itide%h2(i,j) * itide%tideamp(i,j)**2 - enddo; enddo + enddo ; enddo CS%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal_itide',diag%axesT1,Time, & diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index 2952d9ac9b..ef6c160f9f 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -295,7 +295,7 @@ subroutine opacity_from_chl(optics, fluxes, G, CS, chl_in) chl_in(i,j,k), i, j, k, G%geoLonT(i,j), G%geoLatT(i,j) call MOM_error(FATAL,"MOM_opacity opacity_from_chl: "//trim(mesg)) endif - enddo; enddo; enddo + enddo ; enddo ; enddo else ! Only the 2-d surface chlorophyll can be read in from a file. The ! same value is assumed for all layers. diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index cc772bdb53..9906083597 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -391,7 +391,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & visc%Kd_extra_T(i,j,k) = 0.0 visc%Kd_extra_S(i,j,k) = 0.0 endif - enddo; enddo + enddo ; enddo if (associated(dd%KT_extra)) then ; do K=1,nz+1 ; do i=is,ie dd%KT_extra(i,j,K) = KT_extra(i,K) enddo ; enddo ; endif diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index e50d5db614..ec0b5a80b3 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -643,7 +643,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, CS, symmetrize) if (oldfn >= ustarsq) then cycle - else if ((oldfn + Dfn) <= ustarsq) then + elseif ((oldfn + Dfn) <= ustarsq) then Dh = h_at_vel(i,k) else Dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/Dfn) @@ -659,7 +659,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, CS, symmetrize) if (oldfn >= ustarsq) then cycle - else if ((oldfn + Dfn) <= ustarsq) then + elseif ((oldfn + Dfn) <= ustarsq) then Dh = h_at_vel(i,k) else Dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/Dfn) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 2fc99c48fc..1bb8eb48dd 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -273,14 +273,14 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, " "//trim(int_tide_profile_str)//" unavailable in CVMix. Available "//& "profiles in CVMix are "//trim(SIMMONS_PROFILE_STRING)//" and "//& trim(SCHMITTNER_PROFILE_STRING)//".") - else if (.not.CS%use_CVMix_tidal .and. (CS%int_tide_profile == SIMMONS_04.or. & + elseif (.not.CS%use_CVMix_tidal .and. (CS%int_tide_profile == SIMMONS_04.or. & CS%int_tide_profile == SCHMITTNER)) then call MOM_error(FATAL, "tidal_mixing_init: Tidal mixing profiles "// & trim(SIMMONS_PROFILE_STRING)//" and "//trim(SCHMITTNER_PROFILE_STRING)//& " are available only when USE_CVMix_TIDAL is True.") endif - else if (CS%use_CVMix_tidal) then + elseif (CS%use_CVMix_tidal) then call MOM_error(FATAL, "tidal_mixing_init: Cannot set INT_TIDE_DISSIPATION to False "// & "when USE_CVMix_TIDAL is set to True.") endif @@ -294,7 +294,7 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, if (CS%use_CVMix_tidal) then call MOM_error(FATAL, "tidal_mixing_init: Lee wave driven dissipation scheme cannot "// & "be used when CVMix tidal mixing scheme is active.") - end if + endif call get_param(param_file, mdl, "LEE_WAVE_PROFILE", tmpstr, & "LEE_WAVE_PROFILE selects the vertical profile of energy \n"//& "dissipation with LEE_WAVE_DISSIPATION. Valid values are:\n"//& @@ -325,7 +325,7 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, if (CS%use_CVMix_tidal) then call MOM_error(FATAL, "tidal_mixing_init: Polzin scheme cannot "// & "be used when CVMix tidal mixing scheme is active.") - end if + endif call get_param(param_file, mdl, "NU_POLZIN", CS%Nu_Polzin, & "When the Polzin decay profile is used, this is a \n"//& "non-dimensional constant in the expression for the \n"//& @@ -407,7 +407,7 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, if (CS%use_CVMix_tidal) then call MOM_error(FATAL, "tidal_mixing_init: Tidal amplitude files are "// & "not compatible with CVMix tidal mixing. ") - end if + endif call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, & "The path to the file containing the spatially varying \n"//& "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc") @@ -438,7 +438,7 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, ! Compute the fixed part of internal tidal forcing; units are [kg s-2] here. CS%TKE_itidal(i,j) = 0.5*CS%kappa_h2_factor*GV%Rho0*& CS%kappa_itides*CS%h2(i,j)*utide*utide - enddo; enddo + enddo ; enddo endif diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 4226e4fa8c..5f9a8f8281 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1373,7 +1373,7 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, CS) do k=1,nz ; do I=Isq,Ieq ; if (abs(u(I,j,k)) > maxvel) then u(I,j,k) = SIGN(truncvel,u(I,j,k)) if (h(i,j,k) + h(i+1,j,k) > H_report) CS%ntrunc = CS%ntrunc + 1 - endif ; enddo ; enddo + endif ; enddo ; enddo endif ; endif enddo ! j-loop else ! Do not report accelerations leading to large velocities. @@ -1408,7 +1408,7 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, CS) call write_u_accel(I, j, u_old, h, ADp, CDp, dt, G, GV, CS%PointAccel_CSp, & vel_report(I,j), -vel_report(I,j), forces%taux(I,j)*dt_Rho0, & a=CS%a_u(:,j,:), hv=CS%h_u(:,j,:)) - endif ; enddo; enddo + endif ; enddo ; enddo endif if (len_trim(CS%v_trunc_file) > 0) then @@ -1459,7 +1459,7 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, CS) do k=1,nz ; do i=is,ie ; if (abs(v(i,J,k)) > maxvel) then v(i,J,k) = SIGN(truncvel,v(i,J,k)) if (h(i,j,k) + h(i,j+1,k) > H_report) CS%ntrunc = CS%ntrunc + 1 - endif ; enddo ; enddo + endif ; enddo ; enddo endif ; endif enddo ! J-loop else ! Do not report accelerations leading to large velocities. @@ -1494,7 +1494,7 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, CS) call write_v_accel(i, J, v_old, h, ADp, CDp, dt, G, GV, CS%PointAccel_CSp, & vel_report(i,J), -vel_report(i,J), forces%tauy(i,J)*dt_Rho0, & a=CS%a_v(:,J,:),hv=CS%h_v(:,J,:)) - endif ; enddo; enddo + endif ; enddo ; enddo endif end subroutine vertvisc_limit_vel diff --git a/src/tracer/DOME_tracer.F90 b/src/tracer/DOME_tracer.F90 index 38f3f4ee57..d8f6b9a972 100644 --- a/src/tracer/DOME_tracer.F90 +++ b/src/tracer/DOME_tracer.F90 @@ -214,7 +214,7 @@ subroutine initialize_DOME_tracer(restart, day, G, GV, h, diag, OBC, CS, & ! This adds the stripes of tracer to every layer. CS%tr(i,j,k,m) = CS%tr(i,j,k,m) + tr_y enddo - enddo; enddo; enddo + enddo ; enddo ; enddo if (NTR > 7) then do j=js,je ; do i=is,ie diff --git a/src/tracer/MOM_OCMIP2_CO2calc.F90 b/src/tracer/MOM_OCMIP2_CO2calc.F90 index 8c2809418d..aa87c19e73 100644 --- a/src/tracer/MOM_OCMIP2_CO2calc.F90 +++ b/src/tracer/MOM_OCMIP2_CO2calc.F90 @@ -421,7 +421,7 @@ function drtsafe(k0, k1, k2, kb, k1p, k2p, k3p, ksi, kw, ks, kf, & swap=fl fl=fh fh=swap -end if +endif drtsafe=0.5*(x1+x2) dxold=abs(x2-x1) dx=dxold @@ -446,7 +446,7 @@ function drtsafe(k0, k1, k2, kb, k1p, k2p, k3p, ksi, kw, ks, kf, & ! write (6,*) 'Exiting drtsafe at B on iteration ', j, ', ph = ', -log10(drtsafe) return endif - end if + endif if (abs(dx) < xacc) then ! write (6,*) 'Exiting drtsafe at C on iteration ', j, ', ph = ', -log10(drtsafe) return @@ -459,7 +459,7 @@ function drtsafe(k0, k1, k2, kb, k1p, k2p, k3p, ksi, kw, ks, kf, & else xh=drtsafe fh=f - end if + endif enddo !} j return diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 65000627ef..91d96bcc65 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -298,7 +298,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia !Jasmin does not want to apply the maximum for now !if (tr_ptr(i,j,k) > g_tracer%src_var_valid_max) tr_ptr(i,j,k) = g_tracer%src_var_valid_max endif - enddo; enddo ; enddo + enddo ; enddo ; enddo !jgj: Reset CASED to 0 below K=1 if (trim(g_tracer_name) == 'cased') then @@ -306,7 +306,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia if (tr_ptr(i,j,k) /= CS%tracer_land_val) then tr_ptr(i,j,k) = 0.0 endif - enddo; enddo ; enddo + enddo ; enddo ; enddo endif elseif(.not. g_tracer%requires_restart) then !Do nothing for this tracer, it is initialized by the tracer package @@ -521,12 +521,12 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, rho_dzt(:,:,:) = GV%H_to_kg_m2 * GV%Angstrom do k = 1, nk ; do j = jsc, jec ; do i = isc, iec !{ rho_dzt(i,j,k) = GV%H_to_kg_m2 * h_old(i,j,k) - enddo; enddo ; enddo !} + enddo ; enddo ; enddo !} dzt(:,:,:) = 1.0 do k = 1, nk ; do j = jsc, jec ; do i = isc, iec !{ dzt(i,j,k) = GV%H_to_m * h_old(i,j,k) - enddo; enddo ; enddo !} + enddo ; enddo ; enddo !} do j=jsc,jec ; do i=isc,iec surface_field(i,j) = tv%S(i,j,1) diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index b3232c1bca..00e97deb8a 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -394,14 +394,14 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) hEff_sum(:,:) = 0. do k = 1,CS%nsurf-1 ; do j=G%jsc,G%jec ; do i=G%isc-1,G%iec hEff_sum(i,j) = hEff_sum(i,j) + CS%uhEff(i,j,k) - enddo ; enddo; enddo + enddo ; enddo ; enddo call post_data(CS%id_uhEff_2d, hEff_sum, CS%diag) endif if (CS%id_vhEff_2d>0) then hEff_sum(:,:) = 0. do k = 1,CS%nsurf-1 ; do j=G%jsc-1,G%jec ; do i=G%isc,G%iec hEff_sum(i,j) = hEff_sum(i,j) + CS%vhEff(i,j,k) - enddo ; enddo; enddo + enddo ; enddo ; enddo call post_data(CS%id_vhEff_2d, hEff_sum, CS%diag) endif diff --git a/src/tracer/MOM_neutral_diffusion_aux.F90 b/src/tracer/MOM_neutral_diffusion_aux.F90 index 1ecfd7a25a..94ffe5234c 100644 --- a/src/tracer/MOM_neutral_diffusion_aux.F90 +++ b/src/tracer/MOM_neutral_diffusion_aux.F90 @@ -507,7 +507,7 @@ real function refine_nondim_position(CS, T_ref, S_ref, alpha_ref, beta_ref, P_to fa = fb fb = fc fc = fa - end if + endif tol = 2. * machep * abs ( sb ) + CS%xtol m = 0.5 * ( c - sb ) if ( abs ( m ) <= tol .or. fb == 0. ) then @@ -526,12 +526,12 @@ real function refine_nondim_position(CS, T_ref, S_ref, alpha_ref, beta_ref, P_to r = fb / fc p = s0 * ( 2. * m * q * ( q - r ) - ( sb - sa ) * ( r - 1. ) ) q = ( q - 1. ) * ( r - 1. ) * ( s0 - 1. ) - end if + endif if ( 0. < p ) then q = - q else p = - p - end if + endif s0 = e e = d if ( 2. * p < 3. * m * q - abs ( tol * q ) .and. & @@ -540,17 +540,17 @@ real function refine_nondim_position(CS, T_ref, S_ref, alpha_ref, beta_ref, P_to else e = m d = e - end if - end if + endif + endif sa = sb fa = fb if ( tol < abs ( d ) ) then sb = sb + d - else if ( 0. < m ) then + elseif ( 0. < m ) then sb = sb + tol else sb = sb - tol - end if + endif call drho_at_pos(CS, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, & sb, fb) if ( ( 0. < fb .and. 0. < fc ) .or. & @@ -559,7 +559,7 @@ real function refine_nondim_position(CS, T_ref, S_ref, alpha_ref, beta_ref, P_to fc = fa e = sb - sa d = e - end if + endif enddo ! Modified from original to ensure that the minimum is found fa = ABS(fa) ; fb = ABS(fb) ; fc = ABS(fc) diff --git a/src/tracer/MOM_offline_aux.F90 b/src/tracer/MOM_offline_aux.F90 index 4c63ea2b33..deb395bc4a 100644 --- a/src/tracer/MOM_offline_aux.F90 +++ b/src/tracer/MOM_offline_aux.F90 @@ -252,7 +252,7 @@ subroutine distribute_residual_uh_barotropic(G, GV, hvol, uh) else h2d(i,k) = GV%H_subroundoff endif - enddo; enddo + enddo ; enddo ! Distribute flux. Note min/max is intended to make sure that the mass transport ! does not deplete a cell @@ -320,7 +320,7 @@ subroutine distribute_residual_vh_barotropic(G, GV, hvol, vh) else h2d(j,k) = GV%H_subroundoff endif - enddo; enddo + enddo ; enddo ! Distribute flux evenly throughout a column do j=js-1,je diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index 8da247186e..a821219cd5 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -607,7 +607,7 @@ real function remaining_transport_sum(CS, uhtr, vhtr) if (ABS(vhtr(i,J,k))>vh_neglect) then remaining_transport_sum = remaining_transport_sum + ABS(vhtr(i,J,k)) endif - enddo; enddo; enddo + enddo ; enddo ; enddo call sum_across_PEs(remaining_transport_sum) end function remaining_transport_sum @@ -852,15 +852,15 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, CS, h_pre, do k = 1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 eatr_sub(i,j,k) = eatr(i,j,k) ebtr_sub(i,j,k) = ebtr(i,j,k) - enddo; enddo ; enddo + enddo ; enddo ; enddo do k = 1, nz ; do j=js-1,je+1 ; do i=is-2,ie+1 uhtr_sub(I,j,k) = uhtr(I,j,k) - enddo; enddo ; enddo + enddo ; enddo ; enddo do k = 1, nz ; do j=js-2,je+1 ; do i=is-1,ie+1 vhtr_sub(i,J,k) = vhtr(i,J,k) - enddo; enddo ; enddo + enddo ; enddo ; enddo ! Calculate 3d mass transports to be used in this iteration @@ -881,7 +881,7 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, CS, h_pre, call update_h_horizontal_flux(G, GV, uhtr_sub, vhtr_sub, h_pre, h_new) do k = 1, nz ; do i = is-1, ie+1 ; do j=js-1, je+1 h_vol(i,j,k) = h_pre(i,j,k)*G%areaT(i,j) - enddo; enddo; enddo + enddo ; enddo ; enddo call advect_tracer(h_pre, uhtr_sub, vhtr_sub, CS%OBC, dt_iter, G, GV, & CS%tracer_adv_CSp, CS%tracer_Reg, h_vol, max_iter_in=30, x_first_in=x_before_y) @@ -898,7 +898,7 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, CS, h_pre, call update_h_horizontal_flux(G, GV, uhtr_sub, vhtr_sub, h_pre, h_new) do k = 1, nz ; do i = is-1, ie+1 ; do j=js-1, je+1 h_vol(i,j,k) = h_pre(i,j,k)*G%areaT(i,j) - enddo; enddo; enddo + enddo ; enddo ; enddo call advect_tracer(h_pre, uhtr_sub, vhtr_sub, CS%OBC, dt_iter, G, GV, & CS%tracer_adv_CSp, CS%tracer_Reg, h_vol, max_iter_in=30, x_first_in=x_before_y) @@ -922,15 +922,15 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, CS, h_pre, do k = 1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 eatr(i,j,k) = eatr(i,j,k) - eatr_sub(i,j,k) ebtr(i,j,k) = ebtr(i,j,k) - ebtr_sub(i,j,k) - enddo; enddo ; enddo + enddo ; enddo ; enddo do k = 1, nz ; do j=js-1,je+1 ; do i=is-2,ie+1 uhtr(I,j,k) = uhtr(I,j,k) - uhtr_sub(I,j,k) - enddo; enddo ; enddo + enddo ; enddo ; enddo do k = 1, nz ; do j=js-2,je+1 ; do i=is-1,ie+1 vhtr(i,J,k) = vhtr(i,J,k) - vhtr_sub(i,J,k) - enddo; enddo ; enddo + enddo ; enddo ; enddo call pass_var(eatr,G%Domain) call pass_var(ebtr,G%Domain) @@ -946,7 +946,7 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, CS, h_pre, sum_v = sum_v + abs(vhtr(i,J-1,k))+abs(vhtr(I,J,k)) sum_abs_fluxes = sum_abs_fluxes + abs(eatr(i,j,k)) + abs(ebtr(i,j,k)) + abs(uhtr(I-1,j,k)) + & abs(uhtr(I,j,k)) + abs(vhtr(i,J-1,k)) + abs(vhtr(i,J,k)) - enddo; enddo; enddo + enddo ; enddo ; enddo call sum_across_PEs(sum_abs_fluxes) print *, "Remaining u-flux, v-flux:", sum_u, sum_v @@ -958,7 +958,7 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, CS, h_pre, ! Switch order of Strang split every iteration z_first = .not. z_first x_before_y = .not. x_before_y - end do + enddo end subroutine offline_advection_layer @@ -1025,7 +1025,7 @@ subroutine update_offline_fields(CS, h, fluxes, do_ale) if (CS%G%mask2dT(i,j)<1.0) then CS%h_end(i,j,k) = CS%GV%Angstrom endif - enddo; enddo ; enddo + enddo ; enddo ; enddo do k=1,nz+1 ; do j=js,je ; do i=is,ie CS%Kd(i,j,k) = max(0.0, CS%Kd(i,j,k)) @@ -1038,13 +1038,13 @@ subroutine update_offline_fields(CS, h, fluxes, do_ale) if (CS%G%mask2dCv(i,J)<1.0) then CS%vhtr(i,J,k) = 0.0 endif - enddo; enddo ; enddo + enddo ; enddo ; enddo do k=1,nz ; do j=js,je ; do I=is-1,ie if (CS%G%mask2dCu(I,j)<1.0) then CS%uhtr(I,j,k) = 0.0 endif - enddo; enddo ; enddo + enddo ; enddo ; enddo if (CS%debug) then call uvchksum("[uv]htr_sub after update_offline_fields", CS%uhtr, CS%vhtr, CS%G%HI) diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index f61b5a6a5e..7f9b975863 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -85,18 +85,18 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & !$OMP h_old,convert_flux,h_neglect,eb,tr) & !$OMP private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1) !$OMP do - do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo; enddo + do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo ; enddo if (present(sfc_flux)) then if (convert_flux) then !$OMP do do j = js, je; do i = is,ie sfc_src(i,j) = (sfc_flux(i,j)*dt) * GV%kg_m2_to_H - enddo; enddo + enddo ; enddo else !$OMP do do j = js, je; do i = is,ie sfc_src(i,j) = sfc_flux(i,j) - enddo; enddo + enddo ; enddo endif endif if (present(btm_flux)) then @@ -104,12 +104,12 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & !$OMP do do j = js, je; do i = is,ie btm_src(i,j) = (btm_flux(i,j)*dt) * GV%kg_m2_to_H - enddo; enddo + enddo ; enddo else !$OMP do do j = js, je; do i = is,ie btm_src(i,j) = btm_flux(i,j) - enddo; enddo + enddo ; enddo endif endif @@ -271,7 +271,7 @@ subroutine applyTracerBoundaryFluxesInOut(G, GV, Tr, dt, fluxes, h, evap_CFL_lim if (present(in_flux_optional)) then do j=js,je ; do i=is,ie in_flux(i,j) = in_flux_optional(i,j) - enddo; enddo + enddo ; enddo endif if (present(out_flux_optional)) then do j=js,je ; do i=is,ie diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 491803c4e5..0edc123d26 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -723,7 +723,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & tmp = h_srt(i,k2-1,j) ; h_srt(i,k2-1,j) = h_srt(i,k2,j) ; h_srt(i,k2,j) = tmp enddo endif ; enddo - enddo; enddo + enddo ; enddo !$OMP do do j=js-1,je+1 max_srt(j) = 0 diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90 index d2cc4dafbb..68fbd3fdd0 100644 --- a/src/tracer/dye_example.F90 +++ b/src/tracer/dye_example.F90 @@ -233,7 +233,7 @@ subroutine initialize_dye_tracer(restart, day, G, GV, h, diag, OBC, CS, sponge_C z_bot = z_bot + h(i,j,k)*GV%H_to_m enddo endif - enddo; enddo + enddo ; enddo enddo end subroutine initialize_dye_tracer @@ -312,7 +312,7 @@ subroutine dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS z_bot = z_bot + h_new(i,j,k)*GV%H_to_m enddo endif - enddo; enddo + enddo ; enddo enddo end subroutine dye_tracer_column_physics diff --git a/src/tracer/oil_tracer.F90 b/src/tracer/oil_tracer.F90 index b3f595f175..36e503e10c 100644 --- a/src/tracer/oil_tracer.F90 +++ b/src/tracer/oil_tracer.F90 @@ -301,7 +301,7 @@ subroutine initialize_oil_tracer(restart, day, G, GV, h, diag, OBC, CS, & CS%oil_source_i=i CS%oil_source_j=j endif - enddo; enddo + enddo ; enddo CS%Time => day CS%diag => diag diff --git a/src/user/BFB_initialization.F90 b/src/user/BFB_initialization.F90 index 8e6443ae4a..b4d317d289 100644 --- a/src/user/BFB_initialization.F90 +++ b/src/user/BFB_initialization.F90 @@ -75,10 +75,10 @@ subroutine BFB_set_coord(Rlay, g_prime, GV, param_file, eqn_of_state) g_prime(k) = (Rlay(k) - Rlay(k-1))*GV%g_earth/GV%rho0 else g_prime(k) = GV%g_earth - end if + endif !Rlay(:) = 0.0 !g_prime(:) = 0.0 - end do + enddo if (first_call) call write_BFB_log(param_file) diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90 index 7aa2943ff0..e3aa923179 100644 --- a/src/user/BFB_surface_forcing.F90 +++ b/src/user/BFB_surface_forcing.F90 @@ -202,12 +202,12 @@ subroutine BFB_buoyancy_forcing(state, fluxes, day, dt, G, CS) ! density in kg m-3 that is being restored toward. if (G%geoLatT(i,j) < CS%lfrslat) then Temp_restore = CS%SST_s - else if (G%geoLatT(i,j) > CS%lfrnlat) then + elseif (G%geoLatT(i,j) > CS%lfrnlat) then Temp_restore = CS%SST_n else Temp_restore = (CS%SST_s - CS%SST_n)/(CS%lfrslat - CS%lfrnlat) * & (G%geoLatT(i,j) - CS%lfrslat) + CS%SST_s - end if + endif density_restore = Temp_restore*CS%drho_dt + CS%Rho0 diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index 0e9a18ffad..3b30e2ee31 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -71,7 +71,7 @@ subroutine DOME2d_initialize_topography ( D, G, param_file, max_depth ) if ( x <= l1 ) then D(i,j) = bay_depth * max_depth - else if (( x > l1 ) .and. ( x < l2 )) then + elseif (( x > l1 ) .and. ( x < l2 )) then D(i,j) = bay_depth * max_depth + (1.0-bay_depth) * max_depth * & ( x - l1 ) / (l2 - l1) else @@ -453,7 +453,7 @@ subroutine DOME2d_initialize_sponges(G, GV, tv, param_file, use_ALE, CSp, ACSp) h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo - enddo; enddo + enddo ; enddo ! Store the grid on which the T/S sponge data will reside call initialize_ALE_sponge(Idamp, G, param_file, ACSp, h, nz) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index b8d46798e4..34ef50b8cb 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -238,7 +238,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv, just_read_par do j=js,je ; do i=is,ie delta_h = G%bathyT(i,j) / dfloat(nz) h(i,j,:) = GV%m_to_H * delta_h - end do ; end do + enddo ; enddo case default call MOM_error(FATAL,"isomip_initialize: "// & @@ -570,7 +570,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, tv, PF, use_ALE, CSp, ACSp) do j=js,je ; do i=is,ie delta_h = G%bathyT(i,j) / dfloat(nz) h(i,j,:) = delta_h - end do ; end do + enddo ; enddo case default call MOM_error(FATAL,"ISOMIP_initialize_sponges: "// & diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index a41e3b55a2..c464a2b1f6 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -798,14 +798,14 @@ subroutine Surface_Bands_by_data_override(day_center,G,GV,CS) temp_x(i,j)=0.0 temp_y(i,j)=0.0 endif - enddo; enddo + enddo ; enddo ! Interpolate to u/v grids do j = G%jsc,G%jec ; do I = G%IscB,G%IecB CS%STKx0(I,j,b) = 0.5 * (temp_x(i,j) + temp_x(i+1,j)) - enddo; enddo + enddo ; enddo do j = G%JscB,G%JecB ; do i = G%isc,G%iec CS%STKy0(i,J,b) = 0.5 * (temp_y(i,j) + temp_y(i,j+1)) - enddo; enddo + enddo ; enddo ! Disperse into halo on u/v grids call pass_vector(CS%STKx0(:,:,b),CS%STKy0(:,:,b), G%Domain, To_ALL) enddo @@ -895,7 +895,7 @@ subroutine get_Langmuir_Number( LA, G, GV, HBL, USTAR, I, J, & if (.not.(WaveMethod==LF17)) then LA = max(0.1,sqrt(USTAR/(LA_STK+1.e-8))) - end if + endif if (LA_Misalignment) then WaveDirection = atan2(LA_STKy,LA_STKx) diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90 index 736e6d662b..83740a1d61 100644 --- a/src/user/Phillips_initialization.F90 +++ b/src/user/Phillips_initialization.F90 @@ -310,14 +310,14 @@ subroutine Phillips_initialize_topography(D, G, param_file, max_depth) D(i,j) = Htop*sin(PI*(G%geoLonT(i,j)-x1)/(x2-x1))**2 if (G%geoLatT(i,j)>y1 .and. G%geoLatT(i,j)x3 .and. G%geoLonT(i,j)x3 .and. G%geoLonT(i,j)y1 .and. G%geoLatT(i,j) 0. .and. x > (1. - east_sponge_width)) then + elseif (east_sponge_time_scale > 0. .and. x > (1. - east_sponge_width)) then dist = 1. - (1. - x) / east_sponge_width Idamp(i,j) = 1. / east_sponge_time_scale * max(0., min(1., dist)) endif diff --git a/src/user/dumbbell_initialization.F90 b/src/user/dumbbell_initialization.F90 index e2bc9b5869..b19afec76c 100644 --- a/src/user/dumbbell_initialization.F90 +++ b/src/user/dumbbell_initialization.F90 @@ -177,7 +177,7 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, param_file, just_read_param do j=js,je ; do i=is,ie delta_h = G%bathyT(i,j) / dfloat(nz) h(i,j,:) = delta_h - end do ; end do + enddo ; enddo end select diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90 index 2eeda73243..839918270d 100644 --- a/src/user/dumbbell_surface_forcing.F90 +++ b/src/user/dumbbell_surface_forcing.F90 @@ -162,7 +162,7 @@ subroutine dumbbell_buoyancy_forcing(state, fluxes, day, dt, G, CS) ((CS%S_restore(i,j) - state%SSS(i,j)) / & (0.5 * (CS%S_restore(i,j) + state%SSS(i,j)))) - end if + endif enddo ; enddo endif ! end RESTOREBUOY @@ -234,7 +234,7 @@ subroutine dumbbell_dynamic_forcing(state, fluxes, day, dt, G, CS) G%mask2dT(i,j) * sin(deg_rad*(rdays/CS%slp_period)) fluxes%p_surf_full(i,j) = CS%forcing_mask(i,j) * CS%slp_amplitude * & G%mask2dT(i,j) * sin(deg_rad*(rdays/CS%slp_period)) - enddo; enddo + enddo ; enddo @@ -339,7 +339,7 @@ subroutine dumbbell_surface_forcing_init(Time, G, param_file, diag, CS) if ((x>0.25)) then CS%forcing_mask(i,j) = 1 CS%S_restore(i,j) = CS%S_surf + CS%S_range - else if ((x<-0.25)) then + elseif ((x<-0.25)) then CS%forcing_mask(i,j) = 1 CS%S_restore(i,j) = CS%S_surf - CS%S_range endif diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90 index 790185d0ee..017a36bc9a 100644 --- a/src/user/seamount_initialization.F90 +++ b/src/user/seamount_initialization.F90 @@ -178,7 +178,7 @@ subroutine seamount_initialize_thickness ( h, G, GV, param_file, just_read_param do j=js,je ; do i=is,ie delta_h = G%bathyT(i,j) / dfloat(nz) h(i,j,:) = GV%m_to_H * delta_h - end do ; end do + enddo ; enddo end select diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90 index a33718b243..14f31e6916 100644 --- a/src/user/sloshing_initialization.F90 +++ b/src/user/sloshing_initialization.F90 @@ -101,7 +101,7 @@ subroutine sloshing_initialize_thickness ( h, G, GV, param_file, just_read_param ! Define uniform interfaces do k = 0,nz z_unif(k+1) = -real(k)/real(nz) - end do + enddo ! 1. Define stratification n = 3 @@ -117,17 +117,17 @@ subroutine sloshing_initialize_thickness ( h, G, GV, param_file, just_read_param if ( x <= x1 ) then t = y1*x/x1 - else if ( (x > x1 ) .and. ( x < x2 )) then + elseif ( (x > x1 ) .and. ( x < x2 )) then t = y1 + (y2-y1) * (x-x1) / (x2-x1) else t = y2 + (1.0-y2) * (x-x2) / (1.0-x2) - end if + endif t = - z_unif(k) z_inter(k) = -t * G%max_depth - end do + enddo ! 2. Define displacement a0 = 75.0; ! Displacement amplitude (meters) @@ -140,15 +140,15 @@ subroutine sloshing_initialize_thickness ( h, G, GV, param_file, just_read_param if ( k == 1 ) then displ(k) = 0.0 - end if + endif if ( k == nz+1 ) then displ(k) = 0.0 - end if + endif z_inter(k) = z_inter(k) + displ(k) - end do + enddo ! 3. The last interface must coincide with the seabed z_inter(nz+1) = -G%bathyT(i,j) @@ -159,9 +159,9 @@ subroutine sloshing_initialize_thickness ( h, G, GV, param_file, just_read_param if ( z_inter(k) < (z_inter(k+1) + GV%Angstrom_Z) ) then z_inter(k) = z_inter(k+1) + GV%Angstrom_Z - end if + endif - end do + enddo ! 4. Define layers total_height = 0.0 @@ -169,7 +169,7 @@ subroutine sloshing_initialize_thickness ( h, G, GV, param_file, just_read_param h(i,j,k) = GV%m_to_H * (z_inter(k) - z_inter(k+1)) total_height = total_height + h(i,j,k) - end do + enddo enddo ; enddo @@ -234,7 +234,7 @@ subroutine sloshing_initialize_temperature_salinity ( T, S, h, G, GV, param_file !S(:,:,1) = S_ref !do k = 2,G%ke ! S(:,:,k) = S(:,:,k-1) + delta_S - !end do + !enddo deltah = G%max_depth / nz do j=js,je ; do i=is,ie @@ -252,7 +252,7 @@ subroutine sloshing_initialize_temperature_salinity ( T, S, h, G, GV, param_file T(:,:,1) = T_ref do k = 2,G%ke T(:,:,k) = T(:,:,k-1) + delta_T - end do + enddo kdelta = 2 T(:,:,G%ke/2 - (kdelta-1):G%ke/2 + kdelta) = 1.0 diff --git a/src/user/soliton_initialization.F90 b/src/user/soliton_initialization.F90 index 1d4981a003..c9e7eec40e 100644 --- a/src/user/soliton_initialization.F90 +++ b/src/user/soliton_initialization.F90 @@ -55,7 +55,7 @@ subroutine soliton_initialize_thickness(h, G, GV) val4 = val2*((2.0*val3/(1.0+(val3*val3)))**2) h(i,j,k) = GV%m_to_H * (0.25*val4 * (6.0*y*y+3.0) * exp(-0.5*y*y)) enddo - end do ; end do + enddo ; enddo end subroutine soliton_initialize_thickness