From f02deaefb9a9e71eacbe9187e784dad382421f3b Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 5 May 2018 11:16:38 -0400 Subject: [PATCH] Replaced '.eq.' with '==' and '.ne.' with '/=' Replace the older Fortran syntax '.eq.', and '.ne.' with the clearer and more succinct syntax '==' and '/='. All answers are bitwise identical. --- config_src/solo_driver/coupler_types.F90 | 12 +- src/ALE/P1M_functions.F90 | 4 +- src/ALE/P3M_functions.F90 | 18 +- src/ALE/PQM_functions.F90 | 28 +- src/ALE/regrid_edge_values.F90 | 6 +- src/ALE/regrid_solvers.F90 | 2 +- src/core/MOM.F90 | 2 +- src/core/MOM_CoriolisAdv.F90 | 6 +- src/core/MOM_open_boundary.F90 | 2 +- src/diagnostics/MOM_wave_speed.F90 | 22 +- src/diagnostics/MOM_wave_structure.F90 | 6 +- src/equation_of_state/MOM_EOS.F90 | 2 +- src/framework/MOM_diag_mediator.F90 | 16 +- src/framework/MOM_horizontal_regridding.F90 | 50 +- src/ice_shelf/MOM_ice_shelf.F90 | 580 +++++++++--------- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 16 +- src/ice_shelf/shelf_triangular_FEstuff.F90 | 44 +- src/ice_shelf/user_shelf_init.F90 | 2 +- .../MOM_fixed_initialization.F90 | 2 +- src/initialization/midas_vertmap.F90 | 18 +- src/ocean_data_assim/MOM_oda_driver.F90 | 6 +- .../lateral/MOM_internal_tides.F90 | 12 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 4 +- src/parameterizations/vertical/MOM_KPP.F90 | 10 +- .../vertical/MOM_opacity.F90 | 4 +- .../vertical/MOM_tidal_mixing.F90 | 8 +- src/tracer/MOM_OCMIP2_CO2calc.F90 | 6 +- src/tracer/MOM_generic_tracer.F90 | 10 +- src/tracer/MOM_neutral_diffusion.F90 | 4 +- src/tracer/MOM_tracer_advect.F90 | 2 +- src/user/DOME2d_initialization.F90 | 8 +- src/user/DOME_initialization.F90 | 2 +- src/user/MOM_wave_interface.F90 | 32 +- src/user/adjustment_initialization.F90 | 8 +- src/user/dumbbell_initialization.F90 | 2 +- src/user/sloshing_initialization.F90 | 4 +- 36 files changed, 480 insertions(+), 480 deletions(-) diff --git a/config_src/solo_driver/coupler_types.F90 b/config_src/solo_driver/coupler_types.F90 index ba4ce0d3fa..99a74e085c 100644 --- a/config_src/solo_driver/coupler_types.F90 +++ b/config_src/solo_driver/coupler_types.F90 @@ -314,7 +314,7 @@ subroutine coupler_type_copy_1d_2d(var_in, var_out, is, ie, js, je, & if (var_in%num_bcs >= 0) & call CT_spawn_1d_2d(var_in, var_out, (/ is, is, ie, ie /), (/ js, js, je, je /), suffix) - if ((var_out%num_bcs > 0) .and. (diag_name .ne. ' ')) & + if ((var_out%num_bcs > 0) .and. (diag_name /= ' ')) & call CT_set_diags_2d(var_out, diag_name, axes, time) end subroutine coupler_type_copy_1d_2d @@ -365,7 +365,7 @@ subroutine coupler_type_copy_1d_3d(var_in, var_out, is, ie, js, je, kd, & if (var_in%num_bcs >= 0) & call CT_spawn_1d_3d(var_in, var_out, (/ is, is, ie, ie /), (/ js, js, je, je /), (/1, kd/), suffix) - if ((var_out%num_bcs > 0) .and. (diag_name .ne. ' ')) & + if ((var_out%num_bcs > 0) .and. (diag_name /= ' ')) & call CT_set_diags_3d(var_out, diag_name, axes, time) end subroutine coupler_type_copy_1d_3d @@ -408,7 +408,7 @@ subroutine coupler_type_copy_2d_2d(var_in, var_out, is, ie, js, je, & if (var_in%num_bcs >= 0) & call CT_spawn_2d_2d(var_in, var_out, (/ is, is, ie, ie /), (/ js, js, je, je /), suffix) - if ((var_out%num_bcs > 0) .and. (diag_name .ne. ' ')) & + if ((var_out%num_bcs > 0) .and. (diag_name /= ' ')) & call CT_set_diags_2d(var_out, diag_name, axes, time) end subroutine coupler_type_copy_2d_2d @@ -459,7 +459,7 @@ subroutine coupler_type_copy_2d_3d(var_in, var_out, is, ie, js, je, kd, & if (var_in%num_bcs >= 0) & call CT_spawn_2d_3d(var_in, var_out, (/ is, is, ie, ie /), (/ js, js, je, je /), (/1, kd/), suffix) - if ((var_out%num_bcs > 0) .and. (diag_name .ne. ' ')) & + if ((var_out%num_bcs > 0) .and. (diag_name /= ' ')) & call CT_set_diags_3d(var_out, diag_name, axes, time) end subroutine coupler_type_copy_2d_3d @@ -502,7 +502,7 @@ subroutine coupler_type_copy_3d_2d(var_in, var_out, is, ie, js, je, & if (var_in%num_bcs >= 0) & call CT_spawn_3d_2d(var_in, var_out, (/ is, is, ie, ie /), (/ js, js, je, je /), suffix) - if ((var_out%num_bcs > 0) .and. (diag_name .ne. ' ')) & + if ((var_out%num_bcs > 0) .and. (diag_name /= ' ')) & call CT_set_diags_2d(var_out, diag_name, axes, time) end subroutine coupler_type_copy_3d_2d @@ -553,7 +553,7 @@ subroutine coupler_type_copy_3d_3d(var_in, var_out, is, ie, js, je, kd, & if (var_in%num_bcs >= 0) & call CT_spawn_3d_3d(var_in, var_out, (/ is, is, ie, ie /), (/ js, js, je, je /), (/1, kd/), suffix) - if ((var_out%num_bcs > 0) .and. (diag_name .ne. ' ')) & + if ((var_out%num_bcs > 0) .and. (diag_name /= ' ')) & call CT_set_diags_3d(var_out, diag_name, axes, time) end subroutine coupler_type_copy_3d_3d diff --git a/src/ALE/P1M_functions.F90 b/src/ALE/P1M_functions.F90 index 14190bd3ea..aafaec2580 100644 --- a/src/ALE/P1M_functions.F90 +++ b/src/ALE/P1M_functions.F90 @@ -151,7 +151,7 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) ! Using the limited slope, the left edge value is reevaluated and ! the interpolant coefficients recomputed - if ( h0 .NE. 0.0 ) then + if ( h0 /= 0.0 ) then ppoly_E(1,1) = u0 - 0.5 * slope else ppoly_E(1,1) = u0 @@ -177,7 +177,7 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) slope = 2.0 * ( u1 - ppoly_E(N-1,2) ) end if - if ( h1 .NE. 0.0 ) then + if ( h1 /= 0.0 ) then ppoly_E(N,2) = u1 + 0.5 * slope else ppoly_E(N,2) = u1 diff --git a/src/ALE/P3M_functions.F90 b/src/ALE/P3M_functions.F90 index 7ea9f9283b..a532ca7003 100644 --- a/src/ALE/P3M_functions.F90 +++ b/src/ALE/P3M_functions.F90 @@ -133,7 +133,7 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect u_c = u(k) h_c = h(k) - if ( k .EQ. 1 ) then + if ( k == 1 ) then h_l = h(k) u_l = u(k) else @@ -141,7 +141,7 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect u_l = u(k-1) end if - if ( k .EQ. N ) then + if ( k == N ) then h_r = h(k) u_r = u(k) else @@ -190,7 +190,7 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients, h_neglect ! If cubic is not monotonic, monotonize it by modifiying the ! edge slopes, store the new edge slopes and recompute the ! cubic coefficients - if ( monotonic .EQ. 0 ) then + if ( monotonic == 0 ) then call monotonize_cubic( h_c, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ) end if @@ -301,7 +301,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici call build_cubic_interpolant( h, i0, ppoly_E, ppoly_S, ppoly_coefficients ) monotonic = is_cubic_monotonic( ppoly_coefficients, i0 ) - if ( monotonic .EQ. 0 ) then + if ( monotonic == 0 ) then call monotonize_cubic( h0, u0_l, u0_r, 0.0, slope, slope, u1_l, u1_r ) ! Rebuild cubic after monotonization @@ -360,7 +360,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coeffici call build_cubic_interpolant( h, i1, ppoly_E, ppoly_S, ppoly_coefficients ) monotonic = is_cubic_monotonic( ppoly_coefficients, i1 ) - if ( monotonic .EQ. 0 ) then + if ( monotonic == 0 ) then call monotonize_cubic( h1, u0_l, u0_r, slope, 0.0, slope, u1_l, u1_r ) ! Rebuild cubic after monotonization @@ -564,7 +564,7 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ! There is a possible root (and inflexion point) only if a3 is nonzero. ! When a3 is zero, the second derivative of the cubic is constant (the ! cubic degenerates into a parabola) and no inflexion point exists. - if ( a3 .NE. 0.0 ) then + if ( a3 /= 0.0 ) then ! Location of inflexion point xi_ip = - a2 / (3.0 * a3) @@ -579,7 +579,7 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ! decide on which side we want to collapse the inflexion point. ! If the inflexion point lies on one of the edges, the cubic is ! guaranteed to be monotonic - if ( found_ip .EQ. 1 ) then + if ( found_ip == 1 ) then slope_ip = a1 + 2.0*a2*xi_ip + 3.0*a3*xi_ip*xi_ip ! Check whether slope is consistent @@ -597,7 +597,7 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ! 'inflexion_l' and 'inflexion_r' are set to 0 and nothing is to be done. ! Move inflexion point on the left - if ( inflexion_l .EQ. 1 ) then + if ( inflexion_l == 1 ) then u1_l_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_r u1_r_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_l @@ -627,7 +627,7 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r end if ! end treating case with inflexion point on the left ! Move inflexion point on the right - if ( inflexion_r .EQ. 1 ) then + if ( inflexion_r == 1 ) then u1_l_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_r u1_r_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_l diff --git a/src/ALE/PQM_functions.F90 b/src/ALE/PQM_functions.F90 index b91b52b437..8d15e6dd98 100644 --- a/src/ALE/PQM_functions.F90 +++ b/src/ALE/PQM_functions.F90 @@ -191,7 +191,7 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) ! monotonic, edge slopes are consistent and the cell is not an extremum. ! We now need to check and encorce the monotonicity of the quartic within ! the cell - if ( (inflexion_l .EQ. 0) .AND. (inflexion_r .EQ. 0) ) then + if ( (inflexion_l == 0) .AND. (inflexion_r == 0) ) then a = u0_l b = h_c * u1_l @@ -208,7 +208,7 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3 ! Check whether inflexion points exist - if (( alpha1 .ne. 0.0 ) .and. ( rho >= 0.0 )) then + if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) then sqrt_rho = sqrt( rho ) @@ -273,7 +273,7 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) ! If alpha1 is zero, the second derivative of the quartic reduces ! to a straight line - if (( alpha1 .eq. 0.0 ) .and. ( alpha2 .ne. 0.0 )) then + if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then x1 = - alpha3 / alpha2 if ( (x1 >= 0.0) .AND. (x1 <= 1.0) ) then @@ -298,7 +298,7 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) end if ! end checking whether to shift inflexion points ! At this point, we know onto which edge to shift inflexion points - if ( inflexion_l .EQ. 1 ) then + if ( inflexion_l == 1 ) then ! We modify the edge slopes so that both inflexion points ! collapse onto the left edge @@ -323,7 +323,7 @@ subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect ) end if - else if ( inflexion_r .EQ. 1 ) then + else if ( inflexion_r == 1 ) then ! We modify the edge slopes so that both inflexion points ! collapse onto the right edge @@ -608,7 +608,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff ! Compute coefficient for rational function based on mean and right ! edge value and slope - if (u1_r.ne.0.) then ! HACK by AJA + if (u1_r /= 0.) then ! HACK by AJA beta = 2.0 * ( u0_r - um ) / ( (h0 + hNeglect)*u1_r) - 1.0 else beta = 0. @@ -651,7 +651,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff ! Check whether inflexion points exist. If so, transform the quartic ! so that both inflexion points coalesce on the left edge. - if (( alpha1 .ne. 0.0 ) .and. ( rho >= 0.0 )) then + if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) then sqrt_rho = sqrt( rho ) @@ -673,7 +673,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff end if - if (( alpha1 .eq. 0.0 ) .and. ( alpha2 .ne. 0.0 )) then + if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then x1 = - alpha3 / alpha2 if ( (x1 >= 0.0) .and. (x1 <= 1.0) ) then @@ -685,7 +685,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff end if - if ( inflexion_l .eq. 1 ) then + if ( inflexion_l == 1 ) then ! We modify the edge slopes so that both inflexion points ! collapse onto the left edge @@ -757,7 +757,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff ! Compute coefficient for rational function based on mean and left ! edge value and slope - if (um-u0_l.ne.0.) then ! HACK by AJA + if (um-u0_l /= 0.) then ! HACK by AJA beta = 0.5*h1*u1_l / (um-u0_l) - 1.0 else beta = 0. @@ -766,7 +766,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff ar = u0_l ! Right edge value estimate based on rational function - if (1+beta.ne.0.) then ! HACK by AJA + if (1+beta /= 0.) then ! HACK by AJA u0_r = (ar + 2*br + beta*br ) / ((1+beta)*(1+beta)) else u0_r = um + 0.5 * slope ! PLM @@ -804,7 +804,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff ! Check whether inflexion points exist. If so, transform the quartic ! so that both inflexion points coalesce on the right edge. - if (( alpha1 .ne. 0.0 ) .and. ( rho >= 0.0 )) then + if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) then sqrt_rho = sqrt( rho ) @@ -826,7 +826,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff end if - if (( alpha1 .eq. 0.0 ) .and. ( alpha2 .ne. 0.0 )) then + if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then x1 = - alpha3 / alpha2 if ( (x1 >= 0.0) .and. (x1 <= 1.0) ) then @@ -838,7 +838,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coeff end if - if ( inflexion_r .eq. 1 ) then + if ( inflexion_r == 1 ) then ! We modify the edge slopes so that both inflexion points ! collapse onto the right edge diff --git a/src/ALE/regrid_edge_values.F90 b/src/ALE/regrid_edge_values.F90 index d7e8ee54b5..1df7c6ec69 100644 --- a/src/ALE/regrid_edge_values.F90 +++ b/src/ALE/regrid_edge_values.F90 @@ -88,11 +88,11 @@ subroutine bound_edge_values( N, h, u, edge_values, h_neglect ) ! boundary cell and the right neighbor of the right boundary cell ! is assumed to be the same as the right boundary cell. This ! effectively makes boundary cells look like extrema. - if ( k .EQ. 1 ) then + if ( k == 1 ) then k0 = 1 k1 = 1 k2 = 2 - else if ( k .EQ. N ) then + else if ( k == N ) then k0 = N-1 k1 = N k2 = N @@ -179,7 +179,7 @@ subroutine average_discontinuous_edge_values( N, edge_values ) ! Edge value on the right of the edge u0_plus = edge_values(k+1,1) - if ( u0_minus .NE. u0_plus ) then + if ( u0_minus /= u0_plus ) then u0_avg = 0.5 * ( u0_minus + u0_plus ) edge_values(k,2) = u0_avg edge_values(k+1,1) = u0_avg diff --git a/src/ALE/regrid_solvers.F90 b/src/ALE/regrid_solvers.F90 index 1fb85651cf..1f15f97d1b 100644 --- a/src/ALE/regrid_solvers.F90 +++ b/src/ALE/regrid_solvers.F90 @@ -79,7 +79,7 @@ subroutine solve_linear_system( A, B, X, system_size ) ! 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 - if ( k .NE. i ) then + if ( k /= i ) then do j = 1,system_size swap_a = A(i,j) A(i,j) = A(k,j) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e13f52c854..ea22a1832a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1408,7 +1408,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS ! Note that for the layer mode case, the calls to tracer sources and sinks is embedded in ! main_offline_advection_layer. Warning: this may not be appropriate for tracers that ! exchange with the atmosphere - if (time_interval .NE. dt_offline) then + if (time_interval /= dt_offline) then call MOM_error(FATAL, & "For offline tracer mode in a non-ALE configuration, dt_offline must equal time_interval") endif diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index 4019752728..9d01f108d1 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -860,7 +860,7 @@ subroutine gradKE(u, v, h, KE, KEx, KEy, k, OBC, G, CS) ! Calculate KE (Kinetic energy for use in the -grad(KE) acceleration term). - if (CS%KE_Scheme.eq.KE_ARAKAWA) then + if (CS%KE_Scheme == KE_ARAKAWA) then ! The following calculation of Kinetic energy includes the metric terms ! identified in Arakawa & Lamb 1982 as important for KE conservation. It ! also includes the possibility of partially-blocked tracer cell faces. @@ -871,7 +871,7 @@ subroutine gradKE(u, v, h, KE, KEx, KEy, k, OBC, G, CS) +G%areaCv(i,J-1)*(v(i,J-1,k)*v(i,J-1,k)) ) & )*0.25*G%IareaT(i,j) enddo ; enddo - elseif (CS%KE_Scheme.eq.KE_SIMPLE_GUDONOV) then + elseif (CS%KE_Scheme == KE_SIMPLE_GUDONOV) then ! The following discretization of KE is based on the one-dimensinal Gudonov ! scheme which does not take into account any geometric factors do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -881,7 +881,7 @@ subroutine gradKE(u, v, h, KE, KEx, KEy, k, OBC, G, CS) vm = 0.5*( v(i, J ,k) - ABS( v(i, J ,k) ) ) ; vm2 = vm*vm KE(i,j) = ( max(up2,um2) + max(vp2,vm2) ) *0.5 enddo ; enddo - elseif (CS%KE_Scheme.eq.KE_GUDONOV) then + elseif (CS%KE_Scheme == KE_GUDONOV) then ! The following discretization of KE is based on the one-dimensinal Gudonov ! scheme but has been adapted to take horizontal grid factors into account do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 9f5d79ef4e..91f9f6546b 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -303,7 +303,7 @@ subroutine open_boundary_config(G, param_file, OBC) call get_param(param_file, mdl, "NK", OBC%ke, & "The number of model layers", default=0, do_not_log=.true.) - if (config1 .ne. "none") OBC%user_BCs_set_globally = .true. + if (config1 /= "none") OBC%user_BCs_set_globally = .true. if (OBC%number_of_segments > 0) then call get_param(param_file, mdl, "OBC_ZERO_VORTICITY", OBC%zero_vorticity, & diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index bb0c4395f7..ec4a78fc7b 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -849,7 +849,7 @@ subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) enddo ! print resutls (for debugging only) - !if (ig .eq. 83 .and. jg .eq. 2) then + !if (ig == 83 .and. jg == 2) then ! if (nmodes>1)then ! print *, "Results after finding first mode:" ! print *, "first guess at lam_1=", 1./speed2_tot @@ -878,7 +878,7 @@ subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) ! set number of intervals within search range numint = nint((lamMax - lamMin)/lamInc) - !if (ig .eq. 144 .and. jg .eq. 5) then + !if (ig == 144 .and. jg == 5) then ! print *, 'Looking for other eigenvalues at', ig, jg ! print *, 'Wave_speed: lamMin=', lamMin ! print *, 'Wave_speed: cnMax=', 1/sqrt(lamMin) @@ -899,7 +899,7 @@ subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) xl = xr - lamInc call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), & nrows,xr,det_r,ddet_r) - !if (ig .eq. 83 .and. jg .eq. 2) then + !if (ig == 83 .and. jg == 2) then ! print *, "Move interval" ! print *, "iint=",iint ! print *, "@ xr=",xr @@ -911,7 +911,7 @@ subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) nrootsfound = nrootsfound + 1 xbl(nrootsfound) = xl xbr(nrootsfound) = xr - !if (ig .eq. 144 .and. jg .eq. 5) then + !if (ig == 144 .and. jg == 5) then ! print *, "Root located without subdivision!" ! print *, "between xbl=",xl,"and xbr=",xr !endif @@ -939,7 +939,7 @@ subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) nrootsfound = nrootsfound + 1 xbl(nrootsfound) = xl_sub xbr(nrootsfound) = xr - !if (ig .eq. 144 .and. jg .eq. 5) then + !if (ig == 144 .and. jg == 5) then ! print *, "Root located after subdiving",sub_it," times!" ! print *, "between xbl=",xl_sub,"and xbr=",xr !endif @@ -954,7 +954,7 @@ subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) call MOM_error(WARNING, "wave_speed: root not found "// & " after sub_it_max subdivisions of original"// & " interval.") - !if (ig .eq. 144 .and. jg .eq. 5) then + !if (ig == 144 .and. jg == 5) then !print *, "xbl=",xbl !print *, "xbr=",xbr !print *, "Wave_speed: kc=",kc @@ -979,7 +979,7 @@ subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) ! oops, lamMax not large enough - could add code to increase (BDM) ! set unfound modes to zero for now (BDM) cn(i,j,nrootsfound+2:nmodes) = 0.0 - !if (ig .eq. 83 .and. jg .eq. 2) then + !if (ig == 83 .and. jg == 2) then ! call MOM_error(WARNING, "wave_speed: not all modes found "// & ! " within search range: increase numint.") ! print *, "Increase lamMax at ig=",ig," jg=",jg @@ -1030,7 +1030,7 @@ subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) ! ----- Spot check - comment out later (BDM) ---------- !ig = G%idg_offset + i !jg = G%jdg_offset + j - !if (ig .eq. 83 .and. jg .eq. 2) then + !if (ig == 83 .and. jg == 2) then !! print *, "nmodes=",nmodes ! print *, "lam_1=",lam_1 ! print *, "lamMin=",lamMin @@ -1065,9 +1065,9 @@ subroutine tridiag_det(a,b,c,nrows,lam,det_out,ddet_out) real :: I_rescale ! inverse of rescale integer :: n ! row (layer interface) index - if (size(b) .ne. nrows) call MOM_error(WARNING, "Diagonal b must be same length as nrows.") - if (size(a) .ne. nrows) call MOM_error(WARNING, "Diagonal a must be same length as nrows.") - if (size(c) .ne. nrows) call MOM_error(WARNING, "Diagonal c must be same length as nrows.") + if (size(b) /= nrows) call MOM_error(WARNING, "Diagonal b must be same length as nrows.") + if (size(a) /= nrows) call MOM_error(WARNING, "Diagonal a must be same length as nrows.") + if (size(c) /= nrows) call MOM_error(WARNING, "Diagonal c must be same length as nrows.") I_rescale = 1.0/rescale diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index 3286cae5d8..88f5bc06d5 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -275,7 +275,7 @@ subroutine wave_structure(h, tv, G, GV, cn, ModeNum, freq, CS, En, full_halos) do i=is,ie ; if (cn(i,j)>0.0)then !----for debugging, remove later---- ig = i + G%idg_offset ; jg = j + G%jdg_offset - !if (ig .eq. CS%int_tide_source_x .and. jg .eq. CS%int_tide_source_y) then + !if (ig == CS%int_tide_source_x .and. jg == CS%int_tide_source_y) then !----------------------------------- if (G%mask2dT(i,j) > 0.5) then @@ -534,7 +534,7 @@ subroutine wave_structure(h, tv, G, GV, cn, ModeNum, freq, CS, En, full_halos) CS%num_intfaces(i,j) = nzm !----for debugging; delete later---- - !if (ig .eq. ig_stop .and. jg .eq. jg_stop) then + !if (ig == ig_stop .and. jg == jg_stop) then !print *, 'cn(ig,jg)=', cn(i,j) !print *, "e_guess=", e_guess(1:kc-1) !print *, "|e_guess|=", sqrt(sum(e_guess(1:kc-1)**2)) @@ -680,7 +680,7 @@ subroutine tridiag_solver(a,b,c,h,y,method,x) !enddo ; enddo !print *, 'A(2,1),A(2,2),A(1,2)=', A_check(2,1), A_check(2,2), A_check(1,2) !y_check = matmul(A_check,x) - !if (all(y_check .ne. y))then + !if (all(y_check /= y))then ! print *, "tridiag_solver: Uh oh, something's not right!" ! print *, "y=", y ! print *, "y_check=", y_check diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index f3747ff33b..dc6a9869da 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -2332,7 +2332,7 @@ subroutine convert_temp_salt_for_TEOS10(T, S, press, G, kd, mask_z, EOS) if (.not.associated(EOS)) call MOM_error(FATAL, & "convert_temp_salt_to_TEOS10 called with an unassociated EOS_type EOS.") - if ((EOS%form_of_EOS .ne. EOS_TEOS10) .and. (EOS%form_of_EOS .ne. EOS_NEMO)) return + if ((EOS%form_of_EOS /= EOS_TEOS10) .and. (EOS%form_of_EOS /= EOS_NEMO)) return do k=1,kd ; do j=G%jsc,G%jec ; do i=G%isc,G%iec if (mask_z(i,j,k) >= 1.0) then diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index aae2b9d9b1..34bde56f02 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -1357,21 +1357,21 @@ integer function register_diag_field(module_name, field_name, axes, init_time, & ! Register diagnostics remapped to z vertical coordinate if (axes%rank == 3) then remap_axes => null() - if ((axes%id .eq. diag_cs%axesTL%id)) then + if ((axes%id == diag_cs%axesTL%id)) then remap_axes => diag_cs%remap_axesTL(i) - elseif (axes%id .eq. diag_cs%axesBL%id) then + elseif (axes%id == diag_cs%axesBL%id) then remap_axes => diag_cs%remap_axesBL(i) - elseif (axes%id .eq. diag_cs%axesCuL%id ) then + elseif (axes%id == diag_cs%axesCuL%id ) then remap_axes => diag_cs%remap_axesCuL(i) - elseif (axes%id .eq. diag_cs%axesCvL%id) then + elseif (axes%id == diag_cs%axesCvL%id) then remap_axes => diag_cs%remap_axesCvL(i) - elseif (axes%id .eq. diag_cs%axesTi%id) then + elseif (axes%id == diag_cs%axesTi%id) then remap_axes => diag_cs%remap_axesTi(i) - elseif (axes%id .eq. diag_cs%axesBi%id) then + elseif (axes%id == diag_cs%axesBi%id) then remap_axes => diag_cs%remap_axesBi(i) - elseif (axes%id .eq. diag_cs%axesCui%id ) then + elseif (axes%id == diag_cs%axesCui%id ) then remap_axes => diag_cs%remap_axesCui(i) - elseif (axes%id .eq. diag_cs%axesCvi%id) then + elseif (axes%id == diag_cs%axesCvi%id) then remap_axes => diag_cs%remap_axesCvi(i) endif ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 8726ce6770..79e3ebb60d 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -191,15 +191,15 @@ subroutine fill_miss_2d(aout,good,fill,prev,G,smooth,num_pass,relc,crit,keep_bug do j=js,je i_loop: do i=is,ie - if (good_(i,j) .eq. 1.0 .or. fill(i,j) .eq. 0.) cycle i_loop + if (good_(i,j) == 1.0 .or. fill(i,j) == 0.) cycle i_loop ge=good_(i+1,j);gw=good_(i-1,j) gn=good_(i,j+1);gs=good_(i,j-1) east=0.0;west=0.0;north=0.0;south=0.0 - if (ge.eq.1.0) east=aout(i+1,j)*ge - if (gw.eq.1.0) west=aout(i-1,j)*gw - if (gn.eq.1.0) north=aout(i,j+1)*gn - if (gs.eq.1.0) south=aout(i,j-1)*gs + if (ge == 1.0) east=aout(i+1,j)*ge + if (gw == 1.0) west=aout(i-1,j)*gw + if (gn == 1.0) north=aout(i,j+1)*gn + if (gs == 1.0) south=aout(i,j-1)*gs ngood = ge+gw+gn+gs if (ngood > 0.) then @@ -219,13 +219,13 @@ subroutine fill_miss_2d(aout,good,fill,prev,G,smooth,num_pass,relc,crit,keep_bug if (nfill == nfill_prev .and. PRESENT(prev)) then do j=js,je do i=is,ie - if (fill_pts(i,j).eq.1.0) then + if (fill_pts(i,j) == 1.0) then aout(i,j)=prev(i,j) fill_pts(i,j)=0.0 endif enddo enddo - else if (nfill .eq. nfill_prev) then + else if (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'//& @@ -243,7 +243,7 @@ subroutine fill_miss_2d(aout,good,fill,prev,G,smooth,num_pass,relc,crit,keep_bug call pass_var(aout,G%Domain) do j=js,je do i=is,ie - if (fill(i,j) .eq. 1) then + if (fill(i,j) == 1) then east=max(good(i+1,j),fill(i+1,j)) ; west=max(good(i-1,j),fill(i-1,j)) north=max(good(i,j+1),fill(i,j+1)) ; south=max(good(i,j-1),fill(i,j-1)) !### Appropriate parentheses should be added here, but they will change answers. @@ -264,7 +264,7 @@ subroutine fill_miss_2d(aout,good,fill,prev,G,smooth,num_pass,relc,crit,keep_bug do j=js,je do i=is,ie - if (good_(i,j).eq.0.0 .and. fill_pts(i,j) .eq. 1.0) then + if (good_(i,j) == 0.0 .and. fill_pts(i,j) == 1.0) then print *,'in fill_miss, fill, good,i,j= ',fill_pts(i,j),good_(i,j),i,j call MOM_error(FATAL,"MOM_initialize: "// & "fill is true and good is false after fill_miss, how did this happen? ") @@ -348,40 +348,40 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, rcode = NF90_OPEN(filename, NF90_NOWRITE, ncid) - if (rcode .ne. 0) call MOM_error(FATAL,"error opening file "//trim(filename)//& + if (rcode /= 0) call MOM_error(FATAL,"error opening file "//trim(filename)//& " in hinterp_extrap") rcode = NF90_INQ_VARID(ncid, varnam, varid) - if (rcode .ne. 0) call MOM_error(FATAL,"error finding variable "//trim(varnam)//& + if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(varnam)//& " in file "//trim(filename)//" in hinterp_extrap") rcode = NF90_INQUIRE_VARIABLE(ncid, varid, ndims=ndims, dimids=dims) - if (rcode .ne. 0) call MOM_error(FATAL,'error inquiring dimensions hinterp_extrap') + if (rcode /= 0) call MOM_error(FATAL,'error inquiring dimensions hinterp_extrap') if (ndims < 3) call MOM_error(FATAL,"Variable "//trim(varnam)//" in file "// & trim(filename)//" has too few dimensions.") rcode = NF90_INQUIRE_DIMENSION(ncid, dims(1), dim_name(1), len=id) - if (rcode .ne. 0) call MOM_error(FATAL,"error reading dimension 1 data for "// & + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 1 data for "// & trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap") rcode = NF90_INQ_VARID(ncid, dim_name(1), dim_id(1)) - if (rcode .ne. 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(1))//& + if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(1))//& " in file "//trim(filename)//" in hinterp_extrap") rcode = NF90_INQUIRE_DIMENSION(ncid, dims(2), dim_name(2), len=jd) - if (rcode .ne. 0) call MOM_error(FATAL,"error reading dimension 2 data for "// & + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 2 data for "// & trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap") rcode = NF90_INQ_VARID(ncid, dim_name(2), dim_id(2)) - if (rcode .ne. 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(2))//& + if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(2))//& " in file "//trim(filename)//" in hinterp_extrap") rcode = NF90_INQUIRE_DIMENSION(ncid, dims(3), dim_name(3), len=kd) - if (rcode .ne. 0) call MOM_error(FATAL,"error reading dimension 3 data for "// & + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 3 data for "// & trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap") rcode = NF90_INQ_VARID(ncid, dim_name(3), dim_id(3)) - if (rcode .ne. 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(3))//& + if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(3))//& " in file "//trim(filename)//" in hinterp_extrap") missing_value=0.0 rcode = NF90_GET_ATT(ncid, varid, "_FillValue", missing_value) - if (rcode .ne. 0) call MOM_error(FATAL,"error finding missing value for "//& + if (rcode /= 0) call MOM_error(FATAL,"error finding missing value for "//& trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap") if (allocated(lon_in)) deallocate(lon_in) @@ -397,15 +397,15 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, start = 1; count = 1; count(1) = id rcode = NF90_GET_VAR(ncid, dim_id(1), lon_in, start, count) - if (rcode .ne. 0) call MOM_error(FATAL,"error reading dimension 1 values for var_name "// & + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 1 values for var_name "// & trim(varnam)//",dim_name "//trim(dim_name(1))//" in file "// trim(filename)//" in hinterp_extrap") start = 1; count = 1; count(1) = jd rcode = NF90_GET_VAR(ncid, dim_id(2), lat_in, start, count) - if (rcode .ne. 0) call MOM_error(FATAL,"error reading dimension 2 values for var_name "// & + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 2 values for var_name "// & trim(varnam)//",dim_name "//trim(dim_name(2))//" in file "// trim(filename)//" in hinterp_extrap") start = 1; count = 1; count(1) = kd rcode = NF90_GET_VAR(ncid, dim_id(3), z_in, start, count) - if (rcode .ne. 0) call MOM_error(FATAL,"error reading dimension 3 values for var_name "// & + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 3 values for var_name "// & trim(varnam//",dim_name "//trim(dim_name(3)))//" in file "// trim(filename)//" in hinterp_extrap") call cpu_clock_end(id_clock_read) @@ -470,7 +470,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, if (is_root_pe()) then start = 1; start(3) = k; count = 1; count(1) = id; count(2) = jd rcode = NF90_GET_VAR(ncid,varid, tr_in, start, count) - if (rcode .ne. 0) call MOM_error(FATAL,"hinterp_and_extract_from_Fie: "//& + if (rcode /= 0) call MOM_error(FATAL,"hinterp_and_extract_from_Fie: "//& "error reading level "//trim(laynum)//" of variable "//& trim(varnam)//" in file "// trim(filename)) @@ -982,7 +982,7 @@ subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) do j=1,nj do i=1,ni - if (fill(i,j) .eq. 1) then + if (fill(i,j) == 1) then B(i,j,1)=1-nm(i+1,j);B(i,j,2)=1-nm(i-1,j) B(i,j,3)=1-nm(i,j+1);B(i,j,4)=1-nm(i,j-1) endif @@ -992,7 +992,7 @@ subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) do n=1,niter do j=1,nj do i=1,ni - if (fill(i,j) .eq. 1) then + if (fill(i,j) == 1) then bsum = real(B(i,j,1)+B(i,j,2)+B(i,j,3)+B(i,j,4)) Isum = 1.0/bsum res(i,j)=Isum*(B(i,j,1)*mp(i+1,j)+B(i,j,2)*mp(i-1,j)+& diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index d26ca9961a..60f0c688d7 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -313,7 +313,7 @@ function slope_limiter (num, denom) real :: slope_limiter real :: r - if (denom .eq. 0) then + if (denom == 0) then slope_limiter = 0 elseif (num*denom <= 0) then slope_limiter = 0 @@ -820,7 +820,7 @@ subroutine shelf_calc_flux(state, forces, fluxes, Time, time_step, CS) call update_OD_ffrac_uncoupled (CS) endif - if (CS%velocity_update_sub_counter .eq. CS%nstep_velocity) then + if (CS%velocity_update_sub_counter == CS%nstep_velocity) then if (is_root_pe()) write(*,*) "ABOUT TO CALL VELOCITY SOLVER" @@ -877,7 +877,7 @@ subroutine change_thickness_using_melt(CS,G,time_step, fluxes) do j=G%jsc,G%jec do i=G%isc,G%iec - if ((CS%hmask(i,j) .eq. 1) .or. (CS%hmask(i,j) .eq. 2)) then + if ((CS%hmask(i,j) == 1) .or. (CS%hmask(i,j) == 2)) then ! first, zero out fluxes applied during previous time step if (associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0 if (associated(fluxes%sens)) fluxes%sens(i,j) = 0.0 @@ -908,7 +908,7 @@ subroutine change_thickness_using_melt(CS,G,time_step, fluxes) do j=G%jsd,G%jed do i=G%isd,G%ied - if ((CS%hmask(i,j) .eq. 1) .or. (CS%hmask(i,j) .eq. 2)) then + if ((CS%hmask(i,j) == 1) .or. (CS%hmask(i,j) == 2)) then CS%mass_shelf(i,j) = CS%h_shelf(i,j)*CS%density_ice endif enddo @@ -1268,7 +1268,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl call get_param(param_file, mdl, "GROUNDING_LINE_COUPLE", CS%GL_couple, & "THIS PARAMETER NEEDS A DESCRIPTION.", default=.false.) if (CS%GL_regularize) CS%GL_couple = .false. - if (CS%GL_regularize .and. (CS%n_sub_regularize.eq.0)) call MOM_error (FATAL, & + if (CS%GL_regularize .and. (CS%n_sub_regularize == 0)) call MOM_error (FATAL, & "GROUNDING_LINE_INTERP_SUBGRID_N must be a positive integer if GL regularization is used") endif call get_param(param_file, mdl, "SHELF_THERMO", CS%isthermo, & @@ -1667,7 +1667,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl ! next make sure mass is consistent with thickness do j=G%jsd,G%jed do i=G%isd,G%ied - if ((CS%hmask(i,j) .eq. 1) .or. (CS%hmask(i,j) .eq. 2)) then + if ((CS%hmask(i,j) == 1) .or. (CS%hmask(i,j) == 2)) then CS%mass_shelf(i,j) = CS%h_shelf(i,j)*CS%density_ice endif enddo @@ -1704,7 +1704,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl ! next make sure mass is consistent with thickness do j=G%jsd,G%jed do i=G%isd,G%ied - if ((CS%hmask(i,j) .eq. 1) .or. (CS%hmask(i,j) .eq. 2)) then + if ((CS%hmask(i,j) == 1) .or. (CS%hmask(i,j) == 2)) then CS%mass_shelf(i,j) = CS%h_shelf(i,j)*CS%density_ice endif enddo @@ -1729,11 +1729,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl if (.not. G%symmetric) then do j=G%jsd,G%jed do i=G%isd,G%ied - if (((i+G%idg_offset) .eq. (G%domain%nihalo+1)).and.(CS%u_face_mask(i-1,j).eq.3)) then + if (((i+G%idg_offset) == (G%domain%nihalo+1)).and.(CS%u_face_mask(i-1,j) == 3)) then CS%u_shelf (i-1,j-1) = CS%u_boundary_values (i-1,j-1) CS%u_shelf (i-1,j) = CS%u_boundary_values (i-1,j) endif - if (((j+G%jdg_offset) .eq. (G%domain%njhalo+1)).and.(CS%v_face_mask(i,j-1).eq.3)) then + if (((j+G%jdg_offset) == (G%domain%njhalo+1)).and.(CS%v_face_mask(i,j-1) == 3)) then CS%u_shelf (i-1,j-1) = CS%u_boundary_values (i-1,j-1) CS%u_shelf (i,j-1) = CS%u_boundary_values (i,j-1) endif @@ -2240,7 +2240,7 @@ subroutine ice_shelf_advect(CS, time_step, melt_rate, Time) do j=jsd,jed do i=isd,ied thick_bd = CS%thickness_boundary_values(i,j) - if (thick_bd .ne. 0.0) then + if (thick_bd /= 0.0) then CS%h_shelf(i,j) = CS%thickness_boundary_values(i,j) endif enddo @@ -2262,7 +2262,7 @@ subroutine ice_shelf_advect(CS, time_step, melt_rate, Time) do j=jsd,jed do i=isd,ied - if (CS%hmask(i,j) .eq. 1) then + if (CS%hmask(i,j) == 1) then CS%h_shelf (i,j) = h_after_vflux(i,j) endif enddo @@ -2373,7 +2373,7 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) nodefloat = 0 do k=0,1 do l=0,1 - if ((CS%hmask(i,j) .eq. 1) .and. & + if ((CS%hmask(i,j) == 1) .and. & (rhoi/rhow * H_node(i-1+k,j-1+l) - G%bathyT(i,j) <= 0)) then nodefloat = nodefloat + 1 endif @@ -2404,13 +2404,13 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) isym=0 ! must prepare phi - if (FE .eq. 1) then + if (FE == 1) then allocate (Phi (isd:ied,jsd:jed,1:8,1:4)) ; Phi(:,:,:,:)=0 do j=jsd,jed do i=isd,ied - if (((i > isd) .and. (j > jsd)) .or. (isym .eq. 1)) then + if (((i > isd) .and. (j > jsd)) .or. (isym == 1)) then X(:,:) = geolonq (i-1:i,j-1:j)*1000 Y(:,:) = geolatq (i-1:i,j-1:j)*1000 else @@ -2427,7 +2427,7 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) enddo endif - if (FE .eq. 1) then + if (FE == 1) then call calc_shelf_visc_bilinear (CS, u, v) call pass_var (CS%ice_visc_bilinear, G%domain) @@ -2445,7 +2445,7 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) do j=G%jsd,G%jed do i=G%isd,G%ied - if (FE .eq. 1) then + if (FE == 1) then CS%taub_beta_eff_bilinear (i,j) = CS%taub_beta_eff_bilinear (i,j) * CS%float_frac (i,j) else CS%taub_beta_eff_upper_tri (i,j) = CS%taub_beta_eff_upper_tri (i,j) * CS%float_frac (i,j) @@ -2454,20 +2454,20 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) enddo enddo - if (FE .eq. 1) then + if (FE == 1) then call apply_boundary_values_bilinear (CS, time, Phisub, H_node, float_cond, & rhoi/rhow, u_bdry_cont, v_bdry_cont) - elseif (FE .eq. 2) then + elseif (FE == 2) then call apply_boundary_values_triangle (CS, time, u_bdry_cont, v_bdry_cont) endif Au(:,:) = 0.0 ; Av(:,:) = 0.0 - if (FE .eq. 1) then + if (FE == 1) then call CG_action_bilinear (Au, Av, u, v, Phi, Phisub, CS%umask, CS%vmask, CS%hmask, H_node, & CS%ice_visc_bilinear, float_cond, G%bathyT, CS%taub_beta_eff_bilinear, G%areaT, & G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi/rhow) - elseif (FE .eq. 2) then + elseif (FE == 2) then call CG_action_triangular (Au, Av, u, v, CS%umask, CS%vmask, CS%hmask, CS%ice_visc_upper_tri, & CS%ice_visc_lower_tri, CS%taub_beta_eff_upper_tri, CS%taub_beta_eff_lower_tri, & G%dxT, G%dyT, G%areaT, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, isym) @@ -2479,10 +2479,10 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) err_init = 0 ; err_tempu = 0; err_tempv = 0 do j=jsumstart,G%jecB do i=isumstart,G%iecB - if (CS%umask(i,j) .eq. 1) then + if (CS%umask(i,j) == 1) then err_tempu = ABS (Au(i,j) + u_bdry_cont(i,j) - TAUDX(i,j)) endif - if (CS%vmask(i,j) .eq. 1) then + if (CS%vmask(i,j) == 1) then err_tempv = MAX(ABS (Av(i,j) + v_bdry_cont(i,j) - TAUDY(i,j)), err_tempu) endif if (err_tempv >= err_init) then @@ -2513,7 +2513,7 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) if (is_root_pe()) print *,"linear solve done",iters," iterations" - if (FE .eq. 1) then + if (FE == 1) then call calc_shelf_visc_bilinear (CS,u,v) call pass_var (CS%ice_visc_bilinear, G%domain) call pass_var (CS%taub_beta_eff_bilinear, G%domain) @@ -2525,7 +2525,7 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) call pass_var (CS%taub_beta_eff_lower_tri, G%domain) endif - if (iter .eq. 1) then + if (iter == 1) then ! call savearray2 ("visc1",CS%ice_visc_bilinear,CS%write_output_to_file) endif @@ -2533,7 +2533,7 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) do j=G%jsd,G%jed do i=G%isd,G%ied - if (FE .eq. 1) then + if (FE == 1) then CS%taub_beta_eff_bilinear (i,j) = CS%taub_beta_eff_bilinear (i,j) * CS%float_frac (i,j) else CS%taub_beta_eff_upper_tri (i,j) = CS%taub_beta_eff_upper_tri (i,j) * CS%float_frac (i,j) @@ -2544,20 +2544,20 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) u_bdry_cont (:,:) = 0 ; v_bdry_cont (:,:) = 0 - if (FE .eq. 1) then + if (FE == 1) then call apply_boundary_values_bilinear (CS, time, Phisub, H_node, float_cond, & rhoi/rhow, u_bdry_cont, v_bdry_cont) - elseif (FE .eq. 2) then + elseif (FE == 2) then call apply_boundary_values_triangle (CS, time, u_bdry_cont, v_bdry_cont) endif Au(:,:) = 0 ; Av(:,:) = 0 - if (FE .eq. 1) then + if (FE == 1) then call CG_action_bilinear (Au, Av, u, v, Phi, Phisub, CS%umask, CS%vmask, CS%hmask, H_node, & CS%ice_visc_bilinear, float_cond, G%bathyT, CS%taub_beta_eff_bilinear, G%areaT, G%isc-1, & G%iec+1, G%jsc-1, G%jec+1, rhoi/rhow) - elseif (FE .eq. 2) then + elseif (FE == 2) then call CG_action_triangular (Au, Av, u, v, CS%umask, CS%vmask, CS%hmask, CS%ice_visc_upper_tri, & CS%ice_visc_lower_tri, CS%taub_beta_eff_upper_tri, CS%taub_beta_eff_lower_tri, & G%dxT, G%dyT, G%areaT, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, isym) @@ -2565,14 +2565,14 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) err_max = 0 - if (CS%nonlin_solve_err_mode .eq. 1) then + if (CS%nonlin_solve_err_mode == 1) then do j=jsumstart,G%jecB do i=isumstart,G%iecB - if (CS%umask(i,j) .eq. 1) then + if (CS%umask(i,j) == 1) then err_tempu = ABS (Au(i,j) + u_bdry_cont(i,j) - TAUDX(i,j)) endif - if (CS%vmask(i,j) .eq. 1) then + if (CS%vmask(i,j) == 1) then err_tempv = MAX(ABS (Av(i,j) + v_bdry_cont(i,j) - TAUDY(i,j)), err_tempu) endif if (err_tempv >= err_max) then @@ -2583,17 +2583,17 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) call mpp_max (err_max) - elseif (CS%nonlin_solve_err_mode .eq. 2) then + elseif (CS%nonlin_solve_err_mode == 2) then max_vel = 0 ; tempu = 0 ; tempv = 0 do j=jsumstart,G%jecB do i=isumstart,G%iecB - if (CS%umask(i,j) .eq. 1) then + if (CS%umask(i,j) == 1) then err_tempu = ABS (u_last(i,j)-u(i,j)) tempu = u(i,j) endif - if (CS%vmask(i,j) .eq. 1) then + if (CS%vmask(i,j) == 1) then err_tempv = MAX(ABS (v_last(i,j)- v(i,j)), err_tempu) tempv = SQRT(v(i,j)**2+tempu**2) endif @@ -2730,20 +2730,20 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE jsumstart = JSUMSTART_INT_ endif - if (FE .eq. 1) then + if (FE == 1) then visc => CS%ice_visc_bilinear beta => CS%taub_beta_eff_bilinear - elseif (FE .eq. 2) then + elseif (FE == 2) then visc => CS%ice_visc_upper_tri visc_lo => CS%ice_visc_lower_tri beta => CS%taub_beta_eff_upper_tri beta_lo => CS%taub_beta_eff_lower_tri endif - if (FE .eq. 1) then + if (FE == 1) then call apply_boundary_values_bilinear (CS, time, Phisub, H_node, float_cond, & CS%density_ice/CS%density_ocean_avg, ubd, vbd) - elseif (FE .eq. 2) then + elseif (FE == 2) then call apply_boundary_values_triangle (CS, time, ubd, vbd) endif @@ -2754,11 +2754,11 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE call pass_vector(RHSu, RHSv, G%domain, TO_ALL, BGRID_NE) - if (FE .eq. 1) then + if (FE == 1) then call matrix_diagonal_bilinear(CS, float_cond, H_node, & CS%density_ice/CS%density_ocean_avg, Phisub, DIAGu, DIAGv) ! DIAGu(:,:) = 1 ; DIAGv(:,:) = 1 - elseif (FE .eq. 2) then + elseif (FE == 2) then call matrix_diagonal_triangle (CS, DIAGu, DIAGv) DIAGu(:,:) = 1 ; DIAGv(:,:) = 1 endif @@ -2767,11 +2767,11 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE - if (FE .eq. 1) then + if (FE == 1) then call CG_action_bilinear (Au, Av, u, v, Phi, Phisub, umask, vmask, hmask, & H_node, visc, float_cond, G%bathyT, beta, G%areaT, isc-1, iec+1, jsc-1, & jec+1, CS%density_ice/CS%density_ocean_avg) - elseif (FE .eq. 2) then + elseif (FE == 2) then call CG_action_triangular (Au, Av, u, v, umask, vmask, hmask, visc, visc_lo, & beta, beta_lo, G%dxT, G%dyT, G%areaT, isc-1, iec+1, jsc-1, jec+1, isym) endif @@ -2784,8 +2784,8 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=jsumstart,jecq do i=isumstart,iecq - if (umask(i,j) .eq. 1) dot_p1 = dot_p1 + Ru(i,j)**2 - if (vmask(i,j) .eq. 1) dot_p1 = dot_p1 + Rv(i,j)**2 + if (umask(i,j) == 1) dot_p1 = dot_p1 + Ru(i,j)**2 + if (vmask(i,j) == 1) dot_p1 = dot_p1 + Rv(i,j)**2 enddo enddo @@ -2797,8 +2797,8 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=JSUMSTART_INT_,jecq do i=ISUMSTART_INT_,iecq - if (umask(i,j) .eq. 1) sum_vec(i,j) = Ru(i,j)**2 - if (vmask(i,j) .eq. 1) sum_vec(i,j) = sum_vec(i,j) + Rv(i,j)**2 + if (umask(i,j) == 1) sum_vec(i,j) = Ru(i,j)**2 + if (vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + Rv(i,j)**2 enddo enddo @@ -2811,8 +2811,8 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=jsdq,jedq do i=isdq,iedq - if (umask(i,j) .eq. 1) Zu(i,j) = Ru (i,j) / DIAGu (i,j) - if (vmask(i,j) .eq. 1) Zv(i,j) = Rv (i,j) / DIAGv (i,j) + if (umask(i,j) == 1) Zu(i,j) = Ru (i,j) / DIAGu (i,j) + if (vmask(i,j) == 1) Zv(i,j) = Rv (i,j) / DIAGv (i,j) enddo enddo @@ -2843,13 +2843,13 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE Au(:,:) = 0 ; Av(:,:) = 0 - if (FE .eq. 1) then + if (FE == 1) then call CG_action_bilinear (Au, Av, Du, Dv, Phi, Phisub, umask, vmask, hmask, & H_node, visc, float_cond, G%bathyT, beta, G%areaT, is, ie, js, & je, CS%density_ice/CS%density_ocean_avg) - elseif (FE .eq. 2) then + elseif (FE == 2) then call CG_action_triangular (Au, Av, Du, Dv, umask, vmask, hmask, visc, visc_lo, & beta, beta_lo, G%dxT, G%dyT, G%areaT, is, ie, js, je, isym) @@ -2865,11 +2865,11 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE dot_p1 = 0 ; dot_p2 = 0 do j=jsumstart,jecq do i=isumstart,iecq - if (umask(i,j) .eq. 1) then + if (umask(i,j) == 1) then dot_p1 = dot_p1 + Zu(i,j)*Ru(i,j) dot_p2 = dot_p2 + Du(i,j)*Au(i,j) endif - if (vmask(i,j) .eq. 1) then + if (vmask(i,j) == 1) then dot_p1 = dot_p1 + Zv(i,j)*Rv(i,j) dot_p2 = dot_p2 + Dv(i,j)*Av(i,j) endif @@ -2882,12 +2882,12 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=jscq,jecq do i=iscq,iecq - if (umask(i,j) .eq. 1) sum_vec(i,j) = Zu(i,j) * Ru(i,j) - if (vmask(i,j) .eq. 1) sum_vec(i,j) = sum_vec(i,j) + & + if (umask(i,j) == 1) sum_vec(i,j) = Zu(i,j) * Ru(i,j) + if (vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + & Zv(i,j) * Rv(i,j) - if (umask(i,j) .eq. 1) sum_vec_2(i,j) = Du(i,j) * Au(i,j) - if (vmask(i,j) .eq. 1) sum_vec_2(i,j) = sum_vec_2(i,j) + & + if (umask(i,j) == 1) sum_vec_2(i,j) = Du(i,j) * Au(i,j) + if (vmask(i,j) == 1) sum_vec_2(i,j) = sum_vec_2(i,j) + & Dv(i,j) * Av(i,j) enddo enddo @@ -2910,17 +2910,17 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=jsd,jed do i=isd,ied - if (umask(i,j) .eq. 1) u(i,j) = u(i,j) + alpha_k * Du(i,j) - if (vmask(i,j) .eq. 1) v(i,j) = v(i,j) + alpha_k * Dv(i,j) + if (umask(i,j) == 1) u(i,j) = u(i,j) + alpha_k * Du(i,j) + if (vmask(i,j) == 1) v(i,j) = v(i,j) + alpha_k * Dv(i,j) enddo enddo do j=jsd,jed do i=isd,ied - if (umask(i,j) .eq. 1) then + if (umask(i,j) == 1) then Ru_old(i,j) = Ru(i,j) ; Zu_old(i,j) = Zu(i,j) endif - if (vmask(i,j) .eq. 1) then + if (vmask(i,j) == 1) then Rv_old(i,j) = Rv(i,j) ; Zv_old(i,j) = Zv(i,j) endif enddo @@ -2931,18 +2931,18 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=jsd,jed do i=isd,ied - if (umask(i,j) .eq. 1) Ru(i,j) = Ru(i,j) - alpha_k * Au(i,j) - if (vmask(i,j) .eq. 1) Rv(i,j) = Rv(i,j) - alpha_k * Av(i,j) + if (umask(i,j) == 1) Ru(i,j) = Ru(i,j) - alpha_k * Au(i,j) + if (vmask(i,j) == 1) Rv(i,j) = Rv(i,j) - alpha_k * Av(i,j) enddo enddo do j=jsdq,jedq do i=isdq,iedq - if (umask(i,j) .eq. 1) then + if (umask(i,j) == 1) then Zu(i,j) = Ru (i,j) / DIAGu (i,j) endif - if (vmask(i,j) .eq. 1) then + if (vmask(i,j) == 1) then Zv(i,j) = Rv (i,j) / DIAGv (i,j) endif enddo @@ -2956,11 +2956,11 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE dot_p1 = 0 ; dot_p2 = 0 do j=jsumstart,jecq do i=isumstart,iecq - if (umask(i,j) .eq. 1) then + if (umask(i,j) == 1) then dot_p1 = dot_p1 + Zu(i,j)*Ru(i,j) dot_p2 = dot_p2 + Zu_old(i,j)*Ru_old(i,j) endif - if (vmask(i,j) .eq. 1) then + if (vmask(i,j) == 1) then dot_p1 = dot_p1 + Zv(i,j)*Rv(i,j) dot_p2 = dot_p2 + Zv_old(i,j)*Rv_old(i,j) endif @@ -2975,12 +2975,12 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=JSUMSTART_INT_,jecq do i=ISUMSTART_INT_,iecq - if (umask(i,j) .eq. 1) sum_vec(i,j) = Zu(i,j) * Ru(i,j) - if (vmask(i,j) .eq. 1) sum_vec(i,j) = sum_vec(i,j) + & + if (umask(i,j) == 1) sum_vec(i,j) = Zu(i,j) * Ru(i,j) + if (vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + & Zv(i,j) * Rv(i,j) - if (umask(i,j) .eq. 1) sum_vec_2(i,j) = Zu_old(i,j) * Ru_old(i,j) - if (vmask(i,j) .eq. 1) sum_vec_2(i,j) = sum_vec_2(i,j) + & + if (umask(i,j) == 1) sum_vec_2(i,j) = Zu_old(i,j) * Ru_old(i,j) + if (vmask(i,j) == 1) sum_vec_2(i,j) = sum_vec_2(i,j) + & Zv_old(i,j) * Rv_old(i,j) enddo enddo @@ -3002,8 +3002,8 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=jsd,jed do i=isd,ied - if (umask(i,j) .eq. 1) Du(i,j) = Zu(i,j) + beta_k * Du(i,j) - if (vmask(i,j) .eq. 1) Dv(i,j) = Zv(i,j) + beta_k * Dv(i,j) + if (umask(i,j) == 1) Du(i,j) = Zu(i,j) + beta_k * Du(i,j) + if (vmask(i,j) == 1) Dv(i,j) = Zv(i,j) + beta_k * Dv(i,j) enddo enddo @@ -3015,10 +3015,10 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=jsumstart,jecq do i=isumstart,iecq - if (umask(i,j) .eq. 1) then + if (umask(i,j) == 1) then dot_p1 = dot_p1 + Ru(i,j)**2 endif - if (vmask(i,j) .eq. 1) then + if (vmask(i,j) == 1) then dot_p1 = dot_p1 + Rv(i,j)**2 endif enddo @@ -3031,8 +3031,8 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=JSUMSTART_INT_,jecq do i=ISUMSTART_INT_,iecq - if (umask(i,j) .eq. 1) sum_vec(i,j) = Ru(i,j)**2 - if (vmask(i,j) .eq. 1) sum_vec(i,j) = sum_vec(i,j) + Rv(i,j)**2 + if (umask(i,j) == 1) sum_vec(i,j) = Ru(i,j)**2 + if (vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + Rv(i,j)**2 enddo enddo @@ -3058,7 +3058,7 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE cg_halo = cg_halo - 1 - if (cg_halo .eq. 0) then + if (cg_halo == 0) then ! pass vectors call pass_vector(Du, Dv, G%domain, TO_ALL, BGRID_NE) call pass_vector(u, v, G%domain, TO_ALL, BGRID_NE) @@ -3070,15 +3070,15 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE do j=jsdq,jedq do i=isdq,iedq - if (umask(i,j) .eq. 3) then + if (umask(i,j) == 3) then u(i,j) = u_bdry(i,j) - elseif (umask(i,j) .eq. 0) then + elseif (umask(i,j) == 0) then u(i,j) = 0 endif - if (vmask(i,j) .eq. 3) then + if (vmask(i,j) == 3) then v(i,j) = v_bdry(i,j) - elseif (vmask(i,j) .eq. 0) then + elseif (vmask(i,j) == 0) then v(i,j) = 0 endif enddo @@ -3086,7 +3086,7 @@ subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE call pass_vector (u,v, G%domain, TO_ALL, BGRID_NE) - if (conv_flag .eq. 0) then + if (conv_flag == 0) then iters = CS%cg_max_iterations endif @@ -3148,25 +3148,25 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ ((j+j_off) >= G%domain%njhalo+1)) then ! based on mehmet's code - only if btw north & south boundaries stencil(:) = -1 -! if (i+i_off .eq. G%domain%nihalo+G%domain%nihalo) +! if (i+i_off == G%domain%nihalo+G%domain%nihalo) do i=is,ie if (((i+i_off) <= G%domain%niglobal+G%domain%nihalo) .AND. & ((i+i_off) >= G%domain%nihalo+1)) then - if (i+i_off .eq. G%domain%nihalo+1) then + if (i+i_off == G%domain%nihalo+1) then at_west_bdry=.true. else at_west_bdry=.false. endif - if (i+i_off .eq. G%domain%niglobal+G%domain%nihalo) then + if (i+i_off == G%domain%niglobal+G%domain%nihalo) then at_east_bdry=.true. else at_east_bdry=.false. endif - if (hmask(i,j) .eq. 1) then + if (hmask(i,j) == 1) then dxh = G%dxT(i,j) ; dyh = G%dyT(i,j) ; dxdyh = G%areaT(i,j) @@ -3178,7 +3178,7 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ ! 1ST DO LEFT FACE - if (u_face_mask (i-1,j) .eq. 4.) then + if (u_face_mask (i-1,j) == 4.) then flux_diff_cell = flux_diff_cell + dyh * time_step * u_flux_boundary_values (i-1,j) / dxdyh @@ -3187,7 +3187,7 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ ! get u-velocity at center of left face u_face = 0.5 * (CS%u_shelf(i-1,j-1) + CS%u_shelf(i-1,j)) - ! if (at_west_bdry .and. (i .eq. G%isc)) then + ! if (at_west_bdry .and. (i == G%isc)) then ! print *, j, u_face, stencil(-1) ! endif @@ -3195,12 +3195,12 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ ! i may not cover all the cases.. but i cover the realistic ones - if (at_west_bdry .AND. (hmask(i-1,j).eq.3)) then ! at western bdry but there is a + if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then ! at western bdry but there is a ! thickness bdry condition, and the stencil contains it stencil (-1) = CS%thickness_boundary_values(i-1,j) flux_diff_cell = flux_diff_cell + ABS(u_face) * dyh * time_step * stencil(-1) / dxdyh - elseif (hmask(i-1,j) * hmask(i-2,j) .eq. 1) then ! h(i-2) and h(i-1) are valid + elseif (hmask(i-1,j) * hmask(i-2,j) == 1) then ! h(i-2) and h(i-1) are valid phi = slope_limiter (stencil(-1)-stencil(-2), stencil(0)-stencil(-1)) flux_diff_cell = flux_diff_cell + ABS(u_face) * dyh* time_step / dxdyh * & (stencil(-1) - phi * (stencil(-1)-stencil(0))/2) @@ -3214,7 +3214,7 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ endif elseif (u_face < 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available - if (hmask(i-1,j) * hmask(i+1,j) .eq. 1) then ! h(i-1) and h(i+1) are both valid + if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid phi = slope_limiter (stencil(0)-stencil(1), stencil(-1)-stencil(0)) flux_diff_cell = flux_diff_cell - ABS(u_face) * dyh * time_step / dxdyh * & (stencil(0) - phi * (stencil(0)-stencil(-1))/2) @@ -3222,7 +3222,7 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ else flux_diff_cell = flux_diff_cell - ABS(u_face) * dyh * time_step / dxdyh * stencil(0) - if ((hmask(i-1,j) .eq. 0) .OR. (hmask(i-1,j) .eq. 2)) then + if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2)) then flux_enter(i-1,j,2) = ABS(u_face) * dyh * time_step * stencil(0) endif endif @@ -3233,7 +3233,7 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ ! get u-velocity at center of right face - if (u_face_mask (i+1,j) .eq. 4.) then + if (u_face_mask (i+1,j) == 4.) then flux_diff_cell = flux_diff_cell + dyh * time_step * u_flux_boundary_values (i+1,j) / dxdyh @@ -3243,12 +3243,12 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ if (u_face < 0) then !flux is into cell - we need info from h(i+2), h(i+1) if available - if (at_east_bdry .AND. (hmask(i+1,j).eq.3)) then ! at eastern bdry but there is a + if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then ! at eastern bdry but there is a ! thickness bdry condition, and the stencil contains it flux_diff_cell = flux_diff_cell + ABS(u_face) * dyh * time_step * stencil(1) / dxdyh - elseif (hmask(i+1,j) * hmask(i+2,j) .eq. 1) then ! h(i+2) and h(i+1) are valid + elseif (hmask(i+1,j) * hmask(i+2,j) == 1) then ! h(i+2) and h(i+1) are valid phi = slope_limiter (stencil(1)-stencil(2), stencil(0)-stencil(1)) flux_diff_cell = flux_diff_cell + ABS(u_face) * dyh * time_step / dxdyh * & @@ -3264,7 +3264,7 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ elseif (u_face > 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available - if (hmask(i-1,j) * hmask(i+1,j) .eq. 1) then ! h(i-1) and h(i+1) are both valid + if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid phi = slope_limiter (stencil(0)-stencil(-1), stencil(1)-stencil(0)) flux_diff_cell = flux_diff_cell - ABS(u_face) * dyh * time_step / dxdyh * & @@ -3276,7 +3276,7 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ flux_diff_cell = flux_diff_cell - ABS(u_face) * dyh * time_step / dxdyh * stencil(0) - if ((hmask(i+1,j) .eq. 0) .OR. (hmask(i+1,j) .eq. 2)) then + if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2)) then flux_enter(i+1,j,1) = ABS(u_face) * dyh * time_step * stencil(0) endif @@ -3288,29 +3288,29 @@ subroutine ice_shelf_advect_thickness_x (CS, time_step, h0, h_after_uflux, flux_ endif - elseif ((hmask(i,j) .eq. 0) .OR. (hmask(i,j) .eq. 2)) then + elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then - if (at_west_bdry .AND. (hmask(i-1,j) .EQ. 3)) then + if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then u_face = 0.5 * (CS%u_shelf(i-1,j-1) + CS%u_shelf(i-1,j)) flux_enter (i,j,1) = ABS(u_face) * G%dyT(i,j) * time_step * CS%thickness_boundary_values(i-1,j) - elseif (u_face_mask (i-1,j) .eq. 4.) then + elseif (u_face_mask (i-1,j) == 4.) then flux_enter (i,j,1) = G%dyT(i,j) * time_step * u_flux_boundary_values (i-1,j) endif - if (at_east_bdry .AND. (hmask(i+1,j) .EQ. 3)) then + if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then u_face = 0.5 * (CS%u_shelf(i,j-1) + CS%u_shelf(i,j)) flux_enter(i,j,2) = ABS(u_face) * G%dyT(i,j) * time_step * CS%thickness_boundary_values(i+1,j) - elseif (u_face_mask (i+1,j) .eq. 4.) then + elseif (u_face_mask (i+1,j) == 4.) then flux_enter (i,j,2) = G%dyT(i,j) * time_step * u_flux_boundary_values (i+1,j) endif - if ((i .eq. is) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i-1,j) .eq. 1)) then + if ((i == is) .AND. (hmask(i,j) == 0) .AND. (hmask(i-1,j) == 1)) then ! this is solely for the purposes of keeping the mask consistent while advancing ! the front without having to call pass_var - if cell is empty and cell to left ! is ice-covered then this cell will become partly covered hmask(i,j) = 2 - elseif ((i .eq. ie) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i+1,j) .eq. 1)) then + elseif ((i == ie) .AND. (hmask(i,j) == 0) .AND. (hmask(i+1,j) == 1)) then ! this is solely for the purposes of keeping the mask consistent while advancing ! the front without having to call pass_var - if cell is empty and cell to left ! is ice-covered then this cell will become partly covered @@ -3394,19 +3394,19 @@ subroutine ice_shelf_advect_thickness_y (CS, time_step, h_after_uflux, h_after_v if (((j+j_off) <= G%domain%njglobal+G%domain%njhalo) .AND. & ((j+j_off) >= G%domain%njhalo+1)) then - if (j+j_off .eq. G%domain%njhalo+1) then + if (j+j_off == G%domain%njhalo+1) then at_south_bdry=.true. else at_south_bdry=.false. endif - if (j+j_off .eq. G%domain%njglobal+G%domain%njhalo) then + if (j+j_off == G%domain%njglobal+G%domain%njhalo) then at_north_bdry=.true. else at_north_bdry=.false. endif - if (hmask(i,j) .eq. 1) then + if (hmask(i,j) == 1) then dxh = G%dxT(i,j) ; dyh = G%dyT(i,j) ; dxdyh = G%areaT(i,j) h_after_vflux (i,j) = h_after_uflux (i,j) @@ -3415,7 +3415,7 @@ subroutine ice_shelf_advect_thickness_y (CS, time_step, h_after_uflux, h_after_v ! 1ST DO south FACE - if (v_face_mask (i,j-1) .eq. 4.) then + if (v_face_mask (i,j-1) == 4.) then flux_diff_cell = flux_diff_cell + dxh * time_step * v_flux_boundary_values (i,j-1) / dxdyh @@ -3428,11 +3428,11 @@ subroutine ice_shelf_advect_thickness_y (CS, time_step, h_after_uflux, h_after_v ! i may not cover all the cases.. but i cover the realistic ones - if (at_south_bdry .AND. (hmask(i,j-1).eq.3)) then ! at western bdry but there is a + if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then ! at western bdry but there is a ! thickness bdry condition, and the stencil contains it flux_diff_cell = flux_diff_cell + ABS(v_face) * dxh * time_step * stencil(-1) / dxdyh - elseif (hmask(i,j-1) * hmask(i,j-2) .eq. 1) then ! h(j-2) and h(j-1) are valid + elseif (hmask(i,j-1) * hmask(i,j-2) == 1) then ! h(j-2) and h(j-1) are valid phi = slope_limiter (stencil(-1)-stencil(-2), stencil(0)-stencil(-1)) flux_diff_cell = flux_diff_cell + ABS(v_face) * dxh * time_step / dxdyh * & @@ -3446,14 +3446,14 @@ subroutine ice_shelf_advect_thickness_y (CS, time_step, h_after_uflux, h_after_v elseif (v_face < 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available - if (hmask(i,j-1) * hmask(i,j+1) .eq. 1) then ! h(j-1) and h(j+1) are both valid + if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid phi = slope_limiter (stencil(0)-stencil(1), stencil(-1)-stencil(0)) flux_diff_cell = flux_diff_cell - ABS(v_face) * dxh * time_step / dxdyh * & (stencil(0) - phi * (stencil(0)-stencil(-1))/2) else flux_diff_cell = flux_diff_cell - ABS(v_face) * dxh * time_step / dxdyh * stencil(0) - if ((hmask(i,j-1) .eq. 0) .OR. (hmask(i,j-1) .eq. 2)) then + if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2)) then flux_enter(i,j-1,4) = ABS(v_face) * dyh * time_step * stencil(0) endif @@ -3465,7 +3465,7 @@ subroutine ice_shelf_advect_thickness_y (CS, time_step, h_after_uflux, h_after_v ! NEXT DO north FACE - if (v_face_mask(i,j+1) .eq. 4.) then + if (v_face_mask(i,j+1) == 4.) then flux_diff_cell = flux_diff_cell + dxh * time_step * v_flux_boundary_values (i,j+1) / dxdyh @@ -3476,10 +3476,10 @@ subroutine ice_shelf_advect_thickness_y (CS, time_step, h_after_uflux, h_after_v if (v_face < 0) then !flux is into cell - we need info from h(j+2), h(j+1) if available - if (at_north_bdry .AND. (hmask(i,j+1).eq.3)) then ! at eastern bdry but there is a + if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then ! at eastern bdry but there is a ! thickness bdry condition, and the stencil contains it flux_diff_cell = flux_diff_cell + ABS(v_face) * dxh * time_step * stencil(1) / dxdyh - elseif (hmask(i,j+1) * hmask(i,j+2) .eq. 1) then ! h(j+2) and h(j+1) are valid + elseif (hmask(i,j+1) * hmask(i,j+2) == 1) then ! h(j+2) and h(j+1) are valid phi = slope_limiter (stencil(1)-stencil(2), stencil(0)-stencil(1)) flux_diff_cell = flux_diff_cell + ABS(v_face) * dxh * time_step / dxdyh * & (stencil(1) - phi * (stencil(1)-stencil(0))/2) @@ -3491,7 +3491,7 @@ subroutine ice_shelf_advect_thickness_y (CS, time_step, h_after_uflux, h_after_v elseif (v_face > 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available - if (hmask(i,j-1) * hmask(i,j+1) .eq. 1) then ! h(j-1) and h(j+1) are both valid + if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid phi = slope_limiter (stencil(0)-stencil(-1), stencil(1)-stencil(0)) flux_diff_cell = flux_diff_cell - ABS(v_face) * dxh * time_step / dxdyh * & (stencil(0) - phi * (stencil(0)-stencil(1))/2) @@ -3499,7 +3499,7 @@ subroutine ice_shelf_advect_thickness_y (CS, time_step, h_after_uflux, h_after_v ! (o.w. flux would most likely be out of cell) ! but h(j+2) is not flux_diff_cell = flux_diff_cell - ABS(v_face) * dxh * time_step / dxdyh * stencil(0) - if ((hmask(i,j+1) .eq. 0) .OR. (hmask(i,j+1) .eq. 2)) then + if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2)) then flux_enter(i,j+1,3) = ABS(v_face) * dxh * time_step * stencil(0) endif endif @@ -3510,28 +3510,28 @@ subroutine ice_shelf_advect_thickness_y (CS, time_step, h_after_uflux, h_after_v h_after_vflux (i,j) = h_after_vflux (i,j) + flux_diff_cell - elseif ((hmask(i,j) .eq. 0) .OR. (hmask(i,j) .eq. 2)) then + elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then - if (at_south_bdry .AND. (hmask(i,j-1) .EQ. 3)) then + if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then v_face = 0.5 * (CS%u_shelf(i-1,j-1) + CS%u_shelf(i,j-1)) flux_enter (i,j,3) = ABS(v_face) * G%dxT(i,j) * time_step * CS%thickness_boundary_values(i,j-1) - elseif (v_face_mask(i,j-1) .eq. 4.) then + elseif (v_face_mask(i,j-1) == 4.) then flux_enter (i,j,3) = G%dxT(i,j) * time_step * v_flux_boundary_values (i,j-1) endif - if (at_north_bdry .AND. (hmask(i,j+1) .EQ. 3)) then + if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then v_face = 0.5 * (CS%u_shelf(i-1,j) + CS%u_shelf(i,j)) flux_enter (i,j,4) = ABS(v_face) * G%dxT(i,j) * time_step * CS%thickness_boundary_values(i,j+1) - elseif (v_face_mask(i,j+1) .eq. 4.) then + elseif (v_face_mask(i,j+1) == 4.) then flux_enter (i,j,4) = G%dxT(i,j) * time_step * v_flux_boundary_values (i,j+1) endif - if ((j .eq. js) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i,j-1) .eq. 1)) then + if ((j == js) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j-1) == 1)) then ! this is solely for the purposes of keeping the mask consistent while advancing ! the front without having to call pass_var - if cell is empty and cell to left ! is ice-covered then this cell will become partly covered hmask (i,j) = 2 - elseif ((j .eq. je) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i,j+1) .eq. 1)) then + elseif ((j == je) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j+1) == 1)) then ! this is solely for the purposes of keeping the mask consistent while advancing ! the front without having to call pass_var - if cell is empty and cell to left ! is ice-covered then this cell will become partly covered @@ -3612,7 +3612,7 @@ subroutine shelf_advance_front (CS, flux_enter) mapi(1) = -1 ; mapi(2) = 1 ; mapi(3:4) = 0 mapj(3) = -1 ; mapj(4) = 1 ; mapj(1:2) = 0 - do while (iter_flag .eq. 1) + do while (iter_flag == 1) iter_flag = 0 @@ -3664,7 +3664,7 @@ subroutine shelf_advance_front (CS, flux_enter) h_reference = h_reference / real(n_flux) partial_vol = h_shelf (i,j) * area_shelf_h (i,j) + tot_flux - if ((partial_vol / dxdyh) .eq. h_reference) then ! cell is exactly covered, no overflow + if ((partial_vol / dxdyh) == h_reference) then ! cell is exactly covered, no overflow hmask (i,j) = 1 h_shelf (i,j) = h_reference area_shelf_h(i,j) = dxdyh @@ -3689,33 +3689,33 @@ subroutine shelf_advance_front (CS, flux_enter) n_flux = 0 ; new_partial (:) = 0 do k=1,2 - if (u_face_mask (i-2+k,j) .eq. 2) then + if (u_face_mask (i-2+k,j) == 2) then n_flux = n_flux + 1 - elseif (hmask (i+2*k-3,j) .eq. 0) then + elseif (hmask (i+2*k-3,j) == 0) then n_flux = n_flux + 1 new_partial (k) = 1 endif enddo do k=1,2 - if (v_face_mask (i,j-2+k) .eq. 2) then + if (v_face_mask (i,j-2+k) == 2) then n_flux = n_flux + 1 - elseif (hmask (i,j+2*k-3) .eq. 0) then + elseif (hmask (i,j+2*k-3) == 0) then n_flux = n_flux + 1 new_partial (k+2) = 1 endif enddo - if (n_flux .eq. 0) then ! there is nowhere to put the extra ice! + if (n_flux == 0) then ! there is nowhere to put the extra ice! h_shelf(i,j) = h_reference + partial_vol / dxdyh else h_shelf(i,j) = h_reference do k=1,2 - if (new_partial(k) .eq. 1) & + if (new_partial(k) == 1) & flux_enter_replace (i+2*k-3,j,3-k) = partial_vol / real(n_flux) enddo do k=1,2 ! ### Combine these two loops? - if (new_partial(k+2) .eq. 1) & + if (new_partial(k+2) == 1) & flux_enter_replace(i,j+2*k-3,5-k) = partial_vol / real(n_flux) enddo endif @@ -3751,8 +3751,8 @@ subroutine ice_shelf_min_thickness_calve (CS, h_shelf, area_shelf_h,hmask) do j=G%jsd,G%jed do i=G%isd,G%ied -! if ((h_shelf(i,j) < CS%min_thickness_simple_calve) .and. (hmask(i,j).eq.1) .and. & -! (CS%float_frac(i,j) .eq. 0.0)) then +! if ((h_shelf(i,j) < CS%min_thickness_simple_calve) .and. (hmask(i,j) == 1) .and. & +! (CS%float_frac(i,j) == 0.0)) then if ((h_shelf(i,j) < CS%min_thickness_simple_calve) .and. (area_shelf_h(i,j) > 0.)) then h_shelf(i,j) = 0.0 area_shelf_h(i,j) = 0.0 @@ -3775,7 +3775,7 @@ subroutine calve_to_mask (CS, h_shelf, area_shelf_h, hmask, calve_mask) if (CS%calve_to_mask) then do j=G%jsc,G%jec do i=G%isc,G%iec - if ((calve_mask(i,j) .eq. 0.0) .and. (hmask(i,j) .ne. 0.0)) then + if ((calve_mask(i,j) == 0.0) .and. (hmask(i,j) /= 0.0)) then h_shelf(i,j) = 0.0 area_shelf_h(i,j) = 0.0 hmask(i,j) = 0.0 @@ -3874,35 +3874,35 @@ subroutine calc_shelf_driving_stress (CS, TAUD_X, TAUD_Y, OD, FE) dxdyh = G%areaT(i,j) ! print *,dxh," ",dyh," ",dxdyh - if (hmask(i,j) .eq. 1) then ! we are inside the global computational bdry, at an ice-filled cell + if (hmask(i,j) == 1) then ! we are inside the global computational bdry, at an ice-filled cell ! calculate sx - if ((i+i_off) .eq. gisc) then ! at left computational bdry - if (hmask(i+1,j) .eq. 1) then + if ((i+i_off) == gisc) then ! at left computational bdry + if (hmask(i+1,j) == 1) then sx = (S(i+1,j)-S(i,j))/dxh else sx = 0 endif - elseif ((i+i_off) .eq. giec) then ! at right computational bdry - if (hmask(i-1,j) .eq. 1) then + elseif ((i+i_off) == giec) then ! at right computational bdry + if (hmask(i-1,j) == 1) then sx = (S(i,j)-S(i-1,j))/dxh else sx=0 endif else ! interior - if (hmask(i+1,j) .eq. 1) then + if (hmask(i+1,j) == 1) then cnt = cnt+1 sx = S(i+1,j) else sx = S(i,j) endif - if (hmask(i-1,j) .eq. 1) then + if (hmask(i-1,j) == 1) then cnt = cnt+1 sx = sx - S(i-1,j) else sx = sx - S(i,j) endif - if (cnt .eq. 0) then + if (cnt == 0) then sx=0 else sx = sx / (cnt * dxh) @@ -3912,32 +3912,32 @@ subroutine calc_shelf_driving_stress (CS, TAUD_X, TAUD_Y, OD, FE) cnt = 0 ! calculate sy, similarly - if ((j+j_off) .eq. gjsc) then ! at south computational bdry - if (hmask(i,j+1) .eq. 1) then + if ((j+j_off) == gjsc) then ! at south computational bdry + if (hmask(i,j+1) == 1) then sy = (S(i,j+1)-S(i,j))/dyh else sy = 0 endif - elseif ((j+j_off) .eq. gjec) then ! at nprth computational bdry - if (hmask(i,j-1) .eq. 1) then + elseif ((j+j_off) == gjec) then ! at nprth computational bdry + if (hmask(i,j-1) == 1) then sy = (S(i,j)-S(i,j-1))/dyh else sy = 0 endif else ! interior - if (hmask(i,j+1) .eq. 1) then + if (hmask(i,j+1) == 1) then cnt = cnt+1 sy = S(i,j+1) else sy = S(i,j) endif - if (hmask(i,j-1) .eq. 1) then + if (hmask(i,j-1) == 1) then cnt = cnt+1 sy = sy - S(i,j-1) else sy = sy - S(i,j) endif - if (cnt .eq. 0) then + if (cnt == 0) then sy=0 else sy = sy / (cnt * dyh) @@ -3945,7 +3945,7 @@ subroutine calc_shelf_driving_stress (CS, TAUD_X, TAUD_Y, OD, FE) endif - if (FE .eq. 1) then + if (FE == 1) then ! SW vertex taud_x (i-1,j-1) = taud_x (i-1,j-1) - .25 * rho * grav * H(i,j) * sx * dxdyh @@ -3984,14 +3984,14 @@ subroutine calc_shelf_driving_stress (CS, TAUD_X, TAUD_Y, OD, FE) endif - if (float_frac(i,j) .eq. 1) then + if (float_frac(i,j) == 1) then neumann_val = .5 * grav * (rho * H (i,j) ** 2 - rhow * D(i,j) ** 2) else neumann_val = .5 * grav * (1-rho/rhow) * rho * H(i,j) ** 2 endif - if ((u_face_mask(i-1,j) .eq. 2) .OR. (hmask(i-1,j) .eq. 0) .OR. (hmask(i-1,j) .eq. 2) ) then + if ((u_face_mask(i-1,j) == 2) .OR. (hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2) ) then ! left face of the cell is at a stress boundary ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated ! pressure on either side of the face @@ -4005,19 +4005,19 @@ subroutine calc_shelf_driving_stress (CS, TAUD_X, TAUD_Y, OD, FE) taud_x(i-1,j) = taud_x(i-1,j) - .5 * dyh * neumann_val endif - if ((u_face_mask(i,j) .eq. 2) .OR. (hmask(i+1,j) .eq. 0) .OR. (hmask(i+1,j) .eq. 2) ) then + if ((u_face_mask(i,j) == 2) .OR. (hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2) ) then ! right face of the cell is at a stress boundary taud_x(i,j-1) = taud_x(i,j-1) + .5 * dyh * neumann_val taud_x(i,j) = taud_x(i,j) + .5 * dyh * neumann_val endif - if ((v_face_mask(i,j-1) .eq. 2) .OR. (hmask(i,j-1) .eq. 0) .OR. (hmask(i,j-1) .eq. 2) ) then + if ((v_face_mask(i,j-1) == 2) .OR. (hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2) ) then ! south face of the cell is at a stress boundary taud_y(i-1,j-1) = taud_y(i-1,j-1) - .5 * dxh * neumann_val taud_y(i,j-1) = taud_y(i,j-1) - .5 * dxh * neumann_val endif - if ((v_face_mask(i,j) .eq. 2) .OR. (hmask(i,j+1) .eq. 0) .OR. (hmask(i,j+1) .eq. 2) ) then + if ((v_face_mask(i,j) == 2) .OR. (hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2) ) then ! north face of the cell is at a stress boundary taud_y(i-1,j) = taud_y(i-1,j) + .5 * dxh * neumann_val ! note negative sign due to direction of normal vector taud_y(i,j) = taud_y(i,j) + .5 * dxh * neumann_val @@ -4084,17 +4084,17 @@ subroutine init_boundary_values (CS, time, input_flux, input_thick, new_sim) do j=jsd,jed do i=isd,ied -! if ((i .eq. 4) .AND. ((mpp_pe() .eq. 0) .or. (mpp_pe() .eq. 6))) then +! if ((i == 4) .AND. ((mpp_pe() == 0) .or. (mpp_pe() == 6))) then ! print *,hmask(i,j),i,j,mpp_pe() ! endif - if (hmask(i,j) .eq. 3) then + if (hmask(i,j) == 3) then thickness_boundary_values (i,j) = input_thick endif - if ((hmask(i,j) .eq. 0) .or. (hmask(i,j) .eq. 1) .or. (hmask(i,j) .eq. 2)) then + if ((hmask(i,j) == 0) .or. (hmask(i,j) == 1) .or. (hmask(i,j) == 2)) then if ((i <= iec).and.(i >= isc)) then - if (u_face_mask (i-1,j) .eq. 3) then + if (u_face_mask (i-1,j) == 3) then u_boundary_values (i-1,j-1) = (1 - ((G%geoLatBu(i-1,j-1) - 0.5*CS%len_lat)*2./CS%len_lat)**2) * & 1.5 * input_flux / input_thick u_boundary_values (i-1,j) = (1 - ((G%geoLatBu(i-1,j) - 0.5*CS%len_lat)*2./CS%len_lat)**2) * & @@ -4105,12 +4105,12 @@ subroutine init_boundary_values (CS, time, input_flux, input_thick, new_sim) if (.not.(new_sim)) then if (.not. G%symmetric) then - if (((i+i_off) .eq. (G%domain%nihalo+1)).and.(u_face_mask(i-1,j).eq.3)) then + if (((i+i_off) == (G%domain%nihalo+1)).and.(u_face_mask(i-1,j) == 3)) then CS%u_shelf (i-1,j-1) = u_boundary_values (i-1,j-1) CS%u_shelf (i-1,j) = u_boundary_values (i-1,j) ! print *, u_boundary_values (i-1,j) endif - if (((j+j_off) .eq. (G%domain%njhalo+1)).and.(v_face_mask(i,j-1).eq.3)) then + if (((j+j_off) == (G%domain%njhalo+1)).and.(v_face_mask(i,j-1) == 3)) then CS%u_shelf (i-1,j-1) = u_boundary_values (i-1,j-1) CS%u_shelf (i,j-1) = u_boundary_values (i,j-1) endif @@ -4145,14 +4145,14 @@ subroutine CG_action_triangular (uret, vret, u, v, umask, vmask, hmask, nu_upper do i=is,ie do j=js,je - if (hmask(i,j) .eq. 1) then ! this cell's vertices contain degrees of freedom + if (hmask(i,j) == 1) then ! this cell's vertices contain degrees of freedom ux = (u(i,j-1)-u(i-1,j-1))/dxh(i,j) vx = (v(i,j-1)-v(i-1,j-1))/dxh(i,j) uy = (u(i-1,j)-u(i-1,j-1))/dyh(i,j) vy = (v(i-1,j)-v(i-1,j-1))/dyh(i,j) - if (umask(i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node + if (umask(i,j-1) == 1) then ! this (bot right) is a degree of freedom node uret(i,j-1) = uret(i,j-1) + & .5 * dxdyh(i,j) * nu_lower (i,j) * ((4*ux+2*vy) * (1./dxh(i,j)) + (uy+vy) * (0./dyh(i,j))) @@ -4169,7 +4169,7 @@ subroutine CG_action_triangular (uret, vret, u, v, umask, vmask, hmask, nu_upper v(i-1,j) + v(i,j-1)) endif - if (umask(i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node + if (umask(i-1,j) == 1) then ! this (top left) is a degree of freedom node uret(i-1,j) = uret(i-1,j) + & .5 * dxdyh(i,j) * nu_lower (i,j) * ((4*ux+2*vy) * (0./dxh(i,j)) + (uy+vy) * (1./dyh(i,j))) @@ -4186,7 +4186,7 @@ subroutine CG_action_triangular (uret, vret, u, v, umask, vmask, hmask, nu_upper v(i-1,j) + v(i,j-1)) endif - if (umask(i-1,j-1) .eq. 1) then ! this (bot left) is a degree of freedom node + if (umask(i-1,j-1) == 1) then ! this (bot left) is a degree of freedom node uret(i-1,j-1) = uret(i-1,j-1) + & .5 * dxdyh(i,j) * nu_upper (i,j) * ((4*ux+2*vy) * (-1./dxh(i,j)) + (uy+vy) * (-1./dyh(i,j))) @@ -4209,7 +4209,7 @@ subroutine CG_action_triangular (uret, vret, u, v, umask, vmask, hmask, nu_upper uy = (u(i,j)-u(i,j-1))/dyh(i,j) vy = (v(i,j)-v(i,j-1))/dyh(i,j) - if (umask(i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node + if (umask(i,j-1) == 1) then ! this (bot right) is a degree of freedom node uret(i,j-1) = uret(i,j-1) + & .5 * dxdyh(i,j) * nu_upper (i,j) * ((4*ux+2*vy) * (0./dxh(i,j)) + (uy+vy) * (-1./dyh(i,j))) @@ -4226,7 +4226,7 @@ subroutine CG_action_triangular (uret, vret, u, v, umask, vmask, hmask, nu_upper u(i-1,j) + u(i,j-1)) endif - if (umask(i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node + if (umask(i-1,j) == 1) then ! this (top left) is a degree of freedom node uret(i-1,j) = uret(i-1,j) + & .5 * dxdyh(i,j) * nu_upper (i,j) * ((4*ux+2*vy) * (-1./dxh(i,j)) + (uy+vy) * (0./dyh(i,j))) @@ -4243,7 +4243,7 @@ subroutine CG_action_triangular (uret, vret, u, v, umask, vmask, hmask, nu_upper u(i-1,j) + u(i,j-1)) endif - if (umask(i,j) .eq. 1) then ! this (top right) is a degree of freedom node + if (umask(i,j) == 1) then ! this (top right) is a degree of freedom node uret(i,j) = uret(i,j) + & .5 * dxdyh(i,j) * nu_upper (i,j) * ((4*ux+2*vy) * (1./dxh(i,j)) + (uy+vy) * (1./dyh(i,j))) @@ -4297,7 +4297,7 @@ subroutine CG_action_bilinear (uret, vret, u, v, Phi, Phisub, umask, vmask, hmas ! Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q ! Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q -! Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear +! Phi_k is equal to 1 at vertex k, and 0 at vertex l /= k, and bilinear real :: ux, vx, uy, vy, uq, vq, area, basel integer :: iq, jq, iphi, jphi, i, j, ilq, jlq @@ -4307,7 +4307,7 @@ subroutine CG_action_bilinear (uret, vret, u, v, Phi, Phisub, umask, vmask, hmas xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3)) do j=js,je - do i=is,ie ; if (hmask(i,j) .eq. 1) then + do i=is,ie ; if (hmask(i,j) == 1) then ! dxh = G%dxh(i,j) ! dyh = G%dyh(i,j) ! @@ -4329,13 +4329,13 @@ subroutine CG_action_bilinear (uret, vret, u, v, Phi, Phisub, umask, vmask, hmas do iq=1,2 ; do jq=1,2 - if (iq .eq. 2) then + if (iq == 2) then ilq = 2 else ilq = 1 endif - if (jq .eq. 2) then + if (jq == 2) then jlq = 2 else jlq = 1 @@ -4372,41 +4372,41 @@ subroutine CG_action_bilinear (uret, vret, u, v, Phi, Phisub, umask, vmask, hmas v(i,j) * Phi(i,j,8,2*(jq-1)+iq) do iphi=1,2 ; do jphi=1,2 - if (umask (i-2+iphi,j-2+jphi) .eq. 1) then + if (umask (i-2+iphi,j-2+jphi) == 1) then uret (i-2+iphi,j-2+jphi) = uret (i-2+iphi,j-2+jphi) + & .25 * area * nu (i,j) * ((4*ux+2*vy) * Phi(i,j,2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + & (uy+vx) * Phi(i,j,2*(2*(jphi-1)+iphi),2*(jq-1)+iq)) endif - if (vmask (i-2+iphi,j-2+jphi) .eq. 1) then + if (vmask (i-2+iphi,j-2+jphi) == 1) then vret (i-2+iphi,j-2+jphi) = vret (i-2+iphi,j-2+jphi) + & .25 * area * nu (i,j) * ((uy+vx) * Phi(i,j,2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + & (4*vy+2*ux) * Phi(i,j,2*(2*(jphi-1)+iphi),2*(jq-1)+iq)) endif - if (iq .eq. iphi) then + if (iq == iphi) then ilq = 2 else ilq = 1 endif - if (jq .eq. jphi) then + if (jq == jphi) then jlq = 2 else jlq = 1 endif - if (float_cond(i,j) .eq. 0) then + if (float_cond(i,j) == 0) then - if (umask (i-2+iphi,j-2+jphi) .eq. 1) then + if (umask (i-2+iphi,j-2+jphi) == 1) then uret (i-2+iphi,j-2+jphi) = uret (i-2+iphi,j-2+jphi) + & .25 * beta(i,j) * area * uq * xquad(ilq) * xquad(jlq) endif - if (vmask (i-2+iphi,j-2+jphi) .eq. 1) then + if (vmask (i-2+iphi,j-2+jphi) == 1) then vret (i-2+iphi,j-2+jphi) = vret (i-2+iphi,j-2+jphi) + & .25 * beta(i,j) * area * vq * xquad(ilq) * xquad(jlq) @@ -4415,25 +4415,25 @@ subroutine CG_action_bilinear (uret, vret, u, v, Phi, Phisub, umask, vmask, hmas endif Ucontr(iphi,jphi) = Ucontr(iphi,jphi) + .25 * area * uq * xquad(ilq) * xquad(jlq) * beta(i,j) -! if ((i.eq.27) .and. (j.eq.8) .and. (iphi.eq.1) .and. (jphi.eq.1)) & +! if ((i == 27) .and. (j == 8) .and. (iphi == 1) .and. (jphi == 1)) & ! print *, "grid", uq, .25 * area * uq * xquad(ilq) * xquad(jlq) !endif enddo ; enddo enddo ; enddo - if (float_cond(i,j) .eq. 1) then + if (float_cond(i,j) == 1) then Usubcontr = 0.0 ; Vsubcontr = 0.0 ; basel = D(i,j) Ucell(:,:) = u(i-1:i,j-1:j) ; Vcell(:,:) = v(i-1:i,j-1:j) ; Hcell(:,:) = H_node(i-1:i,j-1:j) call CG_action_subgrid_basal_bilinear & (Phisub, Hcell, Ucell, Vcell, area, basel, dens_ratio, Usubcontr, Vsubcontr, i, j) do iphi=1,2 ; do jphi=1,2 - if (umask (i-2+iphi,j-2+jphi) .eq. 1) then + if (umask (i-2+iphi,j-2+jphi) == 1) then uret (i-2+iphi,j-2+jphi) = uret (i-2+iphi,j-2+jphi) + Usubcontr (iphi,jphi) * beta(i,j) endif - if (vmask (i-2+iphi,j-2+jphi) .eq. 1) then + if (vmask (i-2+iphi,j-2+jphi) == 1) then vret (i-2+iphi,j-2+jphi) = vret (i-2+iphi,j-2+jphi) + Vsubcontr (iphi,jphi) * beta(i,j) - !if ( (iphi.eq.1) .and. (jphi.eq.1)) 8 + !if ( (iphi == 1) .and. (jphi == 1)) 8 ! print *, i,j, Usubcontr (iphi,jphi) * beta(i,j), " ", Ucontr(iphi,jphi) endif enddo ; enddo @@ -4497,7 +4497,7 @@ subroutine CG_action_subgrid_basal_bilinear (Phisub, H, U, V, DXDYH, D, dens_rat Ucontr (m,n) = Ucontr (m,n) + subarea * 0.25 * Phisub(i,j,m,n,qx,qy) * uq Vcontr (m,n) = Vcontr (m,n) + subarea * 0.25 * Phisub(i,j,m,n,qx,qy) * vq - ! if ((i_m .eq. 27) .and. (j_m .eq. 8) .and. (m.eq.1) .and. (n.eq.1)) & + ! if ((i_m == 27) .and. (j_m == 8) .and. (m == 1) .and. (n == 1)) & print *, "in subgrid", uq, Phisub(i,j,m,n,qx,qy) endif @@ -4540,12 +4540,12 @@ subroutine matrix_diagonal_triangle (CS, u_diagonal, v_diagonal) nu_lower => CS%ice_visc_lower_tri ; nu_upper => CS%ice_visc_upper_tri beta_lower => CS%taub_beta_eff_lower_tri ; beta_upper => CS%taub_beta_eff_upper_tri - do i=isc-1,iec+1 ; do j=jsc-1,jec+1 ; if (hmask(i,j) .eq. 1) then + do i=isc-1,iec+1 ; do j=jsc-1,jec+1 ; if (hmask(i,j) == 1) then dxh = G%dxT(i,j) dyh = G%dyT(i,j) dxdyh = G%areaT(i,j) - if (umask (i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node + if (umask (i,j-1) == 1) then ! this (bot right) is a degree of freedom node ux = 1./dxh ; uy = 0./dyh vx = 0. ; vy = 0. @@ -4585,7 +4585,7 @@ subroutine matrix_diagonal_triangle (CS, u_diagonal, v_diagonal) endif - if (umask (i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node + if (umask (i-1,j) == 1) then ! this (top left) is a degree of freedom node ux = 0./dxh ; uy = 1./dyh vx = 0. ; vy = 0. @@ -4625,7 +4625,7 @@ subroutine matrix_diagonal_triangle (CS, u_diagonal, v_diagonal) endif - if (umask (i-1,j-1) .eq. 1) then ! this (bot left) is a degree of freedom node + if (umask (i-1,j-1) == 1) then ! this (bot left) is a degree of freedom node ux = -1./dxh ; uy = -1./dyh vx = 0. ; vy = 0. @@ -4646,7 +4646,7 @@ subroutine matrix_diagonal_triangle (CS, u_diagonal, v_diagonal) beta_lower(i,j) * dxdyh * 1./24 endif - if (umask (i,j) .eq. 1) then ! this (top right) is a degree of freedom node + if (umask (i,j) == 1) then ! this (top right) is a degree of freedom node ux = 1./ dxh ; uy = 1./dyh vx = 0. ; vy = 0. @@ -4718,7 +4718,7 @@ subroutine matrix_diagonal_bilinear(CS, float_cond, H_node, dens_ratio, Phisub, ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j - do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (hmask(i,j) .eq. 1) then + do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (hmask(i,j) == 1) then dxh = G%dxT(i,j) dyh = G%dyT(i,j) @@ -4742,19 +4742,19 @@ subroutine matrix_diagonal_bilinear(CS, float_cond, H_node, dens_ratio, Phisub, do iphi=1,2 ; do jphi=1,2 - if (iq .eq. iphi) then + if (iq == iphi) then ilq = 2 else ilq = 1 endif - if (jq .eq. jphi) then + if (jq == jphi) then jlq = 2 else jlq = 1 endif - if (umask (i-2+iphi,j-2+jphi) .eq. 1) then + if (umask (i-2+iphi,j-2+jphi) == 1) then ux = Phi (2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq) uy = Phi (2*(2*(jphi-1)+iphi), 2*(jq-1)+iq) @@ -4767,14 +4767,14 @@ subroutine matrix_diagonal_bilinear(CS, float_cond, H_node, dens_ratio, Phisub, uq = xquad(ilq) * xquad(jlq) - if (float_cond(i,j) .eq. 0) then + if (float_cond(i,j) == 0) then u_diagonal (i-2+iphi,j-2+jphi) = u_diagonal (i-2+iphi,j-2+jphi) + & .25 * beta(i,j) * dxdyh * uq * xquad(ilq) * xquad(jlq) endif endif - if (vmask (i-2+iphi,j-2+jphi) .eq. 1) then + if (vmask (i-2+iphi,j-2+jphi) == 1) then vx = Phi (2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq) vy = Phi (2*(2*(jphi-1)+iphi), 2*(jq-1)+iq) @@ -4787,7 +4787,7 @@ subroutine matrix_diagonal_bilinear(CS, float_cond, H_node, dens_ratio, Phisub, vq = xquad(ilq) * xquad(jlq) - if (float_cond(i,j) .eq. 0) then + if (float_cond(i,j) == 0) then v_diagonal (i-2+iphi,j-2+jphi) = v_diagonal (i-2+iphi,j-2+jphi) + & .25 * beta(i,j) * dxdyh * vq * xquad(ilq) * xquad(jlq) endif @@ -4795,13 +4795,13 @@ subroutine matrix_diagonal_bilinear(CS, float_cond, H_node, dens_ratio, Phisub, endif enddo ; enddo enddo ; enddo - if (float_cond(i,j) .eq. 1) then + if (float_cond(i,j) == 1) then Usubcontr = 0.0 ; Vsubcontr = 0.0 ; basel = G%bathyT(i,j) Hcell(:,:) = H_node(i-1:i,j-1:j) call CG_diagonal_subgrid_basal_bilinear & (Phisub, Hcell, dxdyh, basel, dens_ratio, Usubcontr, Vsubcontr) do iphi=1,2 ; do jphi=1,2 - if (umask (i-2+iphi,j-2+jphi) .eq. 1) then + if (umask (i-2+iphi,j-2+jphi) == 1) then u_diagonal (i-2+iphi,j-2+jphi) = u_diagonal (i-2+iphi,j-2+jphi) + Usubcontr (iphi,jphi) * beta(i,j) v_diagonal (i-2+iphi,j-2+jphi) = v_diagonal (i-2+iphi,j-2+jphi) + Vsubcontr (iphi,jphi) * beta(i,j) endif @@ -4888,9 +4888,9 @@ subroutine apply_boundary_values_triangle (CS, time, u_boundary_contr, v_boundar domain_width = CS%len_lat - do i=isc-1,iec+1 ; do j=jsc-1,jec+1 ; if (hmask(i,j) .eq. 1) then + do i=isc-1,iec+1 ; do j=jsc-1,jec+1 ; if (hmask(i,j) == 1) then - if ((umask(i-1,j-1) .eq. 3) .OR. (umask(i,j-1) .eq. 3) .OR. (umask(i-1,j) .eq. 3)) then + if ((umask(i-1,j-1) == 3) .OR. (umask(i,j-1) == 3) .OR. (umask(i-1,j) == 3)) then dxh = G%dxT(i,j) dyh = G%dyT(i,j) @@ -4901,7 +4901,7 @@ subroutine apply_boundary_values_triangle (CS, time, u_boundary_contr, v_boundar uy = (u_boundary_values(i-1,j)-u_boundary_values(i-1,j-1))/dyh vy = (v_boundary_values(i-1,j)-v_boundary_values(i-1,j-1))/dyh - if (umask (i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node + if (umask (i,j-1) == 1) then ! this (bot right) is a degree of freedom node u_boundary_contr (i,j-1) = u_boundary_contr (i,j-1) + & .5 * dxdyh * nu_lower (i,j) * ((4*ux+2*vy) * (1./dxh) + (uy+vy) * (0./dyh)) @@ -4918,7 +4918,7 @@ subroutine apply_boundary_values_triangle (CS, time, u_boundary_contr, v_boundar v_boundary_values(i-1,j) + v_boundary_values(i,j-1)) endif - if (umask (i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node + if (umask (i-1,j) == 1) then ! this (top left) is a degree of freedom node u_boundary_contr (i-1,j) = u_boundary_contr (i-1,j) + & .5 * dxdyh * nu_lower (i,j) * ((4*ux+2*vy) * (0./dxh) + (uy+vy) * (1./dyh)) @@ -4935,7 +4935,7 @@ subroutine apply_boundary_values_triangle (CS, time, u_boundary_contr, v_boundar v_boundary_values(i-1,j) + v_boundary_values(i,j-1)) endif - if (umask (i-1,j-1) .eq. 1) then ! this (bot left) is a degree of freedom node + if (umask (i-1,j-1) == 1) then ! this (bot left) is a degree of freedom node u_boundary_contr (i-1,j-1) = u_boundary_contr (i-1,j-1) + & .5 * dxdyh * nu_upper (i,j) * ((4*ux+2*vy) * (-1./dxh) + (uy+vy) * (-1./dyh)) @@ -4954,7 +4954,7 @@ subroutine apply_boundary_values_triangle (CS, time, u_boundary_contr, v_boundar endif - if ((umask(i,j) .eq. 3) .OR. (umask(i,j-1) .eq. 3) .OR. (umask(i-1,j) .eq. 3)) then + if ((umask(i,j) == 3) .OR. (umask(i,j-1) == 3) .OR. (umask(i-1,j) == 3)) then dxh = G%dxT(i,j) dyh = G%dyT(i,j) @@ -4965,7 +4965,7 @@ subroutine apply_boundary_values_triangle (CS, time, u_boundary_contr, v_boundar uy = (u_boundary_values(i,j)-u_boundary_values(i,j-1))/dyh vy = (v_boundary_values(i,j)-v_boundary_values(i,j-1))/dyh - if (umask (i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node + if (umask (i,j-1) == 1) then ! this (bot right) is a degree of freedom node u_boundary_contr (i,j-1) = u_boundary_contr (i,j-1) + & .5 * dxdyh * nu_upper (i,j) * ((4*ux+2*vy) * (0./dxh) + (uy+vy) * (-1./dyh)) @@ -4984,7 +4984,7 @@ subroutine apply_boundary_values_triangle (CS, time, u_boundary_contr, v_boundar u_boundary_values(i,j-1)) endif - if (umask (i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node + if (umask (i-1,j) == 1) then ! this (top left) is a degree of freedom node u_boundary_contr (i-1,j) = u_boundary_contr (i-1,j) + & .5 * dxdyh * nu_upper (i,j) * ((4*ux+2*vy) * (-1./dxh) + (uy+vy) * (0./dyh)) @@ -5003,7 +5003,7 @@ subroutine apply_boundary_values_triangle (CS, time, u_boundary_contr, v_boundar u_boundary_values(i,j-1)) endif - if (umask (i,j) .eq. 1) then ! this (top right) is a degree of freedom node + if (umask (i,j) == 1) then ! this (top right) is a degree of freedom node u_boundary_contr (i,j) = u_boundary_contr (i,j) + & .5 * dxdyh * nu_upper (i,j) * ((4*ux+2*vy) * (1./dxh) + (uy+vy) * (1./dyh)) @@ -5081,13 +5081,13 @@ subroutine apply_boundary_values_bilinear(CS, time, Phisub, H_node, float_cond, ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j - do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (hmask(i,j) .eq. 1) then + do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (hmask(i,j) == 1) then ! process this cell if any corners have umask set to non-dirichlet bdry. ! NOTE: vmask not considered, probably should be - if ((umask(i-1,j-1) .eq. 3) .OR. (umask(i,j-1) .eq. 3) .OR. & - (umask(i-1,j) .eq. 3) .OR. (umask(i,j) .eq. 3)) then + if ((umask(i-1,j-1) == 3) .OR. (umask(i,j-1) == 3) .OR. & + (umask(i-1,j) == 3) .OR. (umask(i,j) == 3)) then dxh = G%dxT(i,j) @@ -5144,40 +5144,40 @@ subroutine apply_boundary_values_bilinear(CS, time, Phisub, H_node, float_cond, do iphi=1,2 ; do jphi=1,2 - if (iq .eq. iphi) then + if (iq == iphi) then ilq = 2 else ilq = 1 endif - if (jq .eq. jphi) then + if (jq == jphi) then jlq = 2 else jlq = 1 endif - if (umask (i-2+iphi,j-2+jphi) .eq. 1) then + if (umask (i-2+iphi,j-2+jphi) == 1) then u_boundary_contr (i-2+iphi,j-2+jphi) = u_boundary_contr (i-2+iphi,j-2+jphi) + & .25 * dxdyh * nu (i,j) * ( (4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + & (uy+vx) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq) ) - if (float_cond(i,j) .eq. 0) then + if (float_cond(i,j) == 0) then u_boundary_contr (i-2+iphi,j-2+jphi) = u_boundary_contr (i-2+iphi,j-2+jphi) + & .25 * beta(i,j) * dxdyh * uq * xquad(ilq) * xquad(jlq) endif endif - if (vmask (i-2+iphi,j-2+jphi) .eq. 1) then + if (vmask (i-2+iphi,j-2+jphi) == 1) then v_boundary_contr (i-2+iphi,j-2+jphi) = v_boundary_contr (i-2+iphi,j-2+jphi) + & .25 * dxdyh * nu (i,j) * ( (uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + & (4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq)) - if (float_cond(i,j) .eq. 0) then + if (float_cond(i,j) == 0) then v_boundary_contr (i-2+iphi,j-2+jphi) = v_boundary_contr (i-2+iphi,j-2+jphi) + & .25 * beta(i,j) * dxdyh * vq * xquad(ilq) * xquad(jlq) endif @@ -5186,18 +5186,18 @@ subroutine apply_boundary_values_bilinear(CS, time, Phisub, H_node, float_cond, enddo ; enddo enddo ; enddo - if (float_cond(i,j) .eq. 1) then + if (float_cond(i,j) == 1) then Usubcontr = 0.0 ; Vsubcontr = 0.0 ; basel = G%bathyT(i,j) Ucell(:,:) = u_boundary_values(i-1:i,j-1:j) ; Vcell(:,:) = v_boundary_values(i-1:i,j-1:j) Hcell(:,:) = H_node(i-1:i,j-1:j) call CG_action_subgrid_basal_bilinear & (Phisub, Hcell, Ucell, Vcell, dxdyh, basel, dens_ratio, Usubcontr, Vsubcontr) do iphi=1,2 ; do jphi = 1,2 - if (umask (i-2+iphi,j-2+jphi) .eq. 1) then + if (umask (i-2+iphi,j-2+jphi) == 1) then u_boundary_contr (i-2+iphi,j-2+jphi) = u_boundary_contr (i-2+iphi,j-2+jphi) + & Usubcontr(iphi,jphi) * beta (i,j) endif - if (vmask (i-2+iphi,j-2+jphi) .eq. 1) then + if (vmask (i-2+iphi,j-2+jphi) == 1) then v_boundary_contr (i-2+iphi,j-2+jphi) = v_boundary_contr (i-2+iphi,j-2+jphi) + & Vsubcontr(iphi,jphi) * beta (i,j) endif @@ -5265,7 +5265,7 @@ subroutine calc_shelf_visc_triangular (CS,u,v) dyh = G%dyT(i,j) dxdyh = G%areaT(i,j) - if (hmask (i,j) .eq. 1) then + if (hmask (i,j) == 1) then ux = (u(i,j-1)-u(i-1,j-1)) / dxh vx = (v(i,j-1)-v(i-1,j-1)) / dxh uy = (u(i-1,j)-u(i-1,j-1)) / dyh @@ -5341,7 +5341,7 @@ subroutine calc_shelf_visc_bilinear (CS, u, v) dyh = G%dyT(i,j) dxdyh = G%areaT(i,j) - if (hmask (i,j) .eq. 1) then + if (hmask (i,j) == 1) then ux = (u(i,j) + u(i,j-1) - u(i-1,j) - u(i-1,j-1)) / (2*dxh) vx = (v(i,j) + v(i,j-1) - v(i-1,j) - v(i-1,j-1)) / (2*dxh) uy = (u(i,j) - u(i,j-1) + u(i-1,j) - u(i-1,j-1)) / (2*dyh) @@ -5388,7 +5388,7 @@ subroutine update_OD_ffrac (CS, ocean_mass, counter, nstep_velocity, time_step, enddo enddo - if (counter .eq. nstep_velocity) then + if (counter == nstep_velocity) then do j=jsc,jec do i=isc,iec @@ -5464,7 +5464,7 @@ subroutine bilinear_shape_functions (X, Y, Phi, area) ! ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j -! Phi_i is equal to 1 at vertex i, and 0 at vertex k .ne. i, and bilinear +! Phi_i is equal to 1 at vertex i, and 0 at vertex k /= i, and bilinear ! ! This should be a one-off; once per nonlinear solve? once per lifetime? ! ... will all cells have the same shape and dimension? @@ -5487,13 +5487,13 @@ subroutine bilinear_shape_functions (X, Y, Phi, area) xnode = 2-mod(node,2) ; ynode = ceiling(REAL(node)/2) - if (ynode .eq. 1) then + if (ynode == 1) then yexp = 1-yquad(qpoint) else yexp = yquad(qpoint) endif - if (1 .eq. xnode) then + if (1 == xnode) then xexp = 1-xquad(qpoint) else xexp = xquad(qpoint) @@ -5556,12 +5556,12 @@ subroutine bilinear_shape_functions_subgrid (Phisub, nsub) do k=1,2 do l=1,2 val = 1.0 - if (k .eq. 1) then + if (k == 1) then val = val * (1.0-x) else val = val * x endif - if (l .eq. 1) then + if (l == 1) then val = val * (1.0-y) else val = val * y @@ -5633,7 +5633,7 @@ subroutine update_velocity_masks (CS) do j=js,G%jed do i=is,G%ied - if (hmask(i,j) .eq. 1) then + if (hmask(i,j) == 1) then umask(i-1:i,j-1:j) = 1. vmask(i-1:i,j-1:j) = 1. @@ -5690,40 +5690,40 @@ subroutine update_velocity_masks (CS) ! vmask (i-1,j-1:j) = 0. !endif - !if (j_off+j .eq. gjsc+1) then !bot boundary + !if (j_off+j == gjsc+1) then !bot boundary ! v_face_mask (i,j-1) = 0. ! umask (i-1:i,j-1) = 0. ! vmask (i-1:i,j-1) = 0. - !elseif (j_off+j .eq. gjec) then !top boundary + !elseif (j_off+j == gjec) then !top boundary ! v_face_mask (i,j) = 0. ! umask (i-1:i,j) = 0. ! vmask (i-1:i,j) = 0. !endif if (i < G%ied) then - if ((hmask(i+1,j) .eq. 0) & - .OR. (hmask(i+1,j) .eq. 2)) then + if ((hmask(i+1,j) == 0) & + .OR. (hmask(i+1,j) == 2)) then !right boundary or adjacent to unfilled cell u_face_mask (i,j) = 2. endif endif if (i > G%isd) then - if ((hmask(i-1,j) .eq. 0) .OR. (hmask(i-1,j) .eq. 2)) then + if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2)) then !adjacent to unfilled cell u_face_mask (i-1,j) = 2. endif endif if (j > G%jsd) then - if ((hmask(i,j-1) .eq. 0) .OR. (hmask(i,j-1) .eq. 2)) then + if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2)) then !adjacent to unfilled cell v_face_mask (i,j-1) = 2. endif endif if (j < G%jed) then - if ((hmask(i,j+1) .eq. 0) .OR. (hmask(i,j+1) .eq. 2)) then + if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2)) then !adjacent to unfilled cell v_face_mask (i,j) = 2. endif @@ -5768,7 +5768,7 @@ subroutine interpolate_H_to_B (CS, h_shelf, hmask, H_node) num_h = 0 do k=0,1 do l=0,1 - if (hmask (i+k,j+l) .eq. 1.0) then + if (hmask (i+k,j+l) == 1.0) then summ = summ + h_shelf (i+k,j+l) num_h = num_h + 1 endif @@ -5866,7 +5866,7 @@ subroutine savearray2(fname,A,flag) END DO - if (i.eq.1) THEN + if (i == 1) THEN lh = LEN(TRIM(ln)) @@ -5893,7 +5893,7 @@ subroutine savearray2(fname,A,flag) WRITE(UNIT=fin,IOSTAT=iock,FMT=TRIM(FMT1)) TRIM(ln) - if (iock .ne. 0) THEN + if (iock /= 0) THEN PRINT*,iock END IF END DO @@ -5949,7 +5949,7 @@ subroutine solo_time_step (CS, time_step, n, Time, min_time_step_in) local_u_max = 0 ; local_v_max = 0 - if (hmask (i,j) .eq. 1.0) then + if (hmask (i,j) == 1.0) then ! all 4 corners of the cell should have valid velocity values; otherwise something is wrong ! this is done by checking that umask and vmask are nonzero at all 4 corners do ki=1,2 ; do kj = 1,2 @@ -5991,7 +5991,7 @@ subroutine solo_time_step (CS, time_step, n, Time, min_time_step_in) call ice_shelf_advect (CS, time_step_int, CS%lprec, Time) - if (mpp_pe() .eq. 7) then + if (mpp_pe() == 7) then call savearray2 ("hmask",CS%hmask,CS%write_output_to_file) !!! OVS!!! ! call savearray2 ("tshelf",CS%t_shelf,CS%write_output_to_file) @@ -6109,7 +6109,7 @@ subroutine ice_shelf_temp(CS, time_step, melt_rate, Time) do i=isd,ied t_bd = CS%t_boundary_values(i,j) ! if (CS%hmask(i,j) > 1) then - if ((CS%hmask(i,j) .eq. 3) .or. (CS%hmask(i,j) .eq. -2)) then + if ((CS%hmask(i,j) == 3) .or. (CS%hmask(i,j) == -2)) then CS%t_shelf(i,j) = CS%t_boundary_values(i,j) endif enddo @@ -6140,7 +6140,7 @@ subroutine ice_shelf_temp(CS, time_step, melt_rate, Time) do j=jsd,jed do i=isd,ied -! if (CS%hmask(i,j) .eq. 1) then +! if (CS%hmask(i,j) == 1) then if (CS%h_shelf(i,j) > 0.0) then CS%t_shelf (i,j) = th_after_vflux(i,j)/CS%h_shelf (i,j) else @@ -6153,7 +6153,7 @@ subroutine ice_shelf_temp(CS, time_step, melt_rate, Time) do i=isd,ied t_bd = CS%t_boundary_values(i,j) ! if (CS%hmask(i,j) > 1) then - if ((CS%hmask(i,j) .eq. 3) .or. (CS%hmask(i,j) .eq. -2)) then + if ((CS%hmask(i,j) == 3) .or. (CS%hmask(i,j) == -2)) then CS%t_shelf(i,j) = t_bd ! CS%t_shelf(i,j) = -15.0 endif @@ -6162,7 +6162,7 @@ subroutine ice_shelf_temp(CS, time_step, melt_rate, Time) do j=jsc,jec do i=isc,iec - if ((CS%hmask(i,j) .eq. 1) .or. (CS%hmask(i,j) .eq. 2)) then + if ((CS%hmask(i,j) == 1) .or. (CS%hmask(i,j) == 2)) then if (CS%h_shelf(i,j) > 0.0) then ! CS%t_shelf (i,j) = CS%t_shelf (i,j) + time_step*(adot*Tsurf -melt_rate (i,j)*Tbot(i,j))/CS%h_shelf (i,j) CS%t_shelf (i,j) = CS%t_shelf (i,j) + time_step*(adot*Tsurf -3/spy*Tbot(i,j))/CS%h_shelf (i,j) @@ -6248,25 +6248,25 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter ((j+j_off) >= G%domain%njhalo+1)) then ! based on mehmet's code - only if btw north & south boundaries stencil(:) = -1 -! if (i+i_off .eq. G%domain%nihalo+G%domain%nihalo) +! if (i+i_off == G%domain%nihalo+G%domain%nihalo) do i=is,ie if (((i+i_off) <= G%domain%niglobal+G%domain%nihalo) .AND. & ((i+i_off) >= G%domain%nihalo+1)) then - if (i+i_off .eq. G%domain%nihalo+1) then + if (i+i_off == G%domain%nihalo+1) then at_west_bdry=.true. else at_west_bdry=.false. endif - if (i+i_off .eq. G%domain%niglobal+G%domain%nihalo) then + if (i+i_off == G%domain%niglobal+G%domain%nihalo) then at_east_bdry=.true. else at_east_bdry=.false. endif - if (hmask(i,j) .eq. 1) then + if (hmask(i,j) == 1) then dxh = G%dxT(i,j) ; dyh = G%dyT(i,j) ; dxdyh = G%areaT(i,j) @@ -6278,7 +6278,7 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter ! 1ST DO LEFT FACE - if (u_face_mask (i-1,j) .eq. 4.) then + if (u_face_mask (i-1,j) == 4.) then flux_diff_cell = flux_diff_cell + dyh * time_step * u_flux_boundary_values (i-1,j) * & t_boundary(i-1,j) / dxdyh @@ -6290,7 +6290,7 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter ! get u-velocity at center of left face u_face = 0.5 * (CS%u_shelf(i-1,j-1) + CS%u_shelf(i-1,j)) - ! if (at_west_bdry .and. (i .eq. G%isc)) then + ! if (at_west_bdry .and. (i == G%isc)) then ! print *, j, u_face, stencil(-1) ! endif @@ -6298,12 +6298,12 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter ! i may not cover all the cases.. but i cover the realistic ones - if (at_west_bdry .AND. (hmask(i-1,j).eq.3)) then ! at western bdry but there is a + if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then ! at western bdry but there is a ! thickness bdry condition, and the stencil contains it stencil (-1) = CS%t_boundary_values(i-1,j)*CS%h_shelf(i-1,j) flux_diff_cell = flux_diff_cell + ABS(u_face) * dyh * time_step * stencil(-1) / dxdyh - elseif (hmask(i-1,j) * hmask(i-2,j) .eq. 1) then ! h(i-2) and h(i-1) are valid + elseif (hmask(i-1,j) * hmask(i-2,j) == 1) then ! h(i-2) and h(i-1) are valid phi = slope_limiter (stencil(-1)-stencil(-2), stencil(0)-stencil(-1)) flux_diff_cell = flux_diff_cell + ABS(u_face) * dyh* time_step / dxdyh * & (stencil(-1) - phi * (stencil(-1)-stencil(0))/2) @@ -6317,7 +6317,7 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter endif elseif (u_face < 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available - if (hmask(i-1,j) * hmask(i+1,j) .eq. 1) then ! h(i-1) and h(i+1) are both valid + if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid phi = slope_limiter (stencil(0)-stencil(1), stencil(-1)-stencil(0)) flux_diff_cell = flux_diff_cell - ABS(u_face) * dyh * time_step / dxdyh * & (stencil(0) - phi * (stencil(0)-stencil(-1))/2) @@ -6325,7 +6325,7 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter else flux_diff_cell = flux_diff_cell - ABS(u_face) * dyh * time_step / dxdyh * stencil(0) - if ((hmask(i-1,j) .eq. 0) .OR. (hmask(i-1,j) .eq. 2)) then + if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2)) then flux_enter(i-1,j,2) = ABS(u_face) * dyh * time_step * stencil(0) endif endif @@ -6336,7 +6336,7 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter ! get u-velocity at center of right face - if (u_face_mask (i+1,j) .eq. 4.) then + if (u_face_mask (i+1,j) == 4.) then flux_diff_cell = flux_diff_cell + dyh * time_step * u_flux_boundary_values (i+1,j) *& t_boundary(i+1,j)/ dxdyh @@ -6349,12 +6349,12 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter if (u_face < 0) then !flux is into cell - we need info from h(i+2), h(i+1) if available - if (at_east_bdry .AND. (hmask(i+1,j).eq.3)) then ! at eastern bdry but there is a + if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then ! at eastern bdry but there is a ! thickness bdry condition, and the stencil contains it flux_diff_cell = flux_diff_cell + ABS(u_face) * dyh * time_step * stencil(1) / dxdyh - elseif (hmask(i+1,j) * hmask(i+2,j) .eq. 1) then ! h(i+2) and h(i+1) are valid + elseif (hmask(i+1,j) * hmask(i+2,j) == 1) then ! h(i+2) and h(i+1) are valid phi = slope_limiter (stencil(1)-stencil(2), stencil(0)-stencil(1)) flux_diff_cell = flux_diff_cell + ABS(u_face) * dyh * time_step / dxdyh * & @@ -6370,7 +6370,7 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter elseif (u_face > 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available - if (hmask(i-1,j) * hmask(i+1,j) .eq. 1) then ! h(i-1) and h(i+1) are both valid + if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid phi = slope_limiter (stencil(0)-stencil(-1), stencil(1)-stencil(0)) flux_diff_cell = flux_diff_cell - ABS(u_face) * dyh * time_step / dxdyh * & @@ -6382,7 +6382,7 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter flux_diff_cell = flux_diff_cell - ABS(u_face) * dyh * time_step / dxdyh * stencil(0) - if ((hmask(i+1,j) .eq. 0) .OR. (hmask(i+1,j) .eq. 2)) then + if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2)) then flux_enter(i+1,j,1) = ABS(u_face) * dyh * time_step * stencil(0) endif @@ -6395,34 +6395,34 @@ subroutine ice_shelf_advect_temp_x (CS, time_step, h0, h_after_uflux, flux_enter endif - elseif ((hmask(i,j) .eq. 0) .OR. (hmask(i,j) .eq. 2)) then + elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then - if (at_west_bdry .AND. (hmask(i-1,j) .EQ. 3)) then + if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then u_face = 0.5 * (CS%u_shelf(i-1,j-1) + CS%u_shelf(i-1,j)) flux_enter (i,j,1) = ABS(u_face) * G%dyT(i,j) * time_step * t_boundary(i-1,j)* & CS%thickness_boundary_values(i+1,j) - elseif (u_face_mask (i-1,j) .eq. 4.) then + elseif (u_face_mask (i-1,j) == 4.) then flux_enter (i,j,1) = G%dyT(i,j) * time_step * u_flux_boundary_values (i-1,j)*t_boundary(i-1,j) ! flux_enter (i,j,1) = G%dyh(i,j) * time_step * CS%u_shelf(i,j)*t_boundary (i-1,j) ! assume no flux bc for temp endif - if (at_east_bdry .AND. (hmask(i+1,j) .EQ. 3)) then + if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then u_face = 0.5 * (CS%u_shelf(i,j-1) + CS%u_shelf(i,j)) flux_enter(i,j,2) = ABS(u_face) * G%dyT(i,j) * time_step * t_boundary(i+1,j)* & CS%thickness_boundary_values(i+1,j) - elseif (u_face_mask (i+1,j) .eq. 4.) then + elseif (u_face_mask (i+1,j) == 4.) then flux_enter (i,j,2) = G%dyT(i,j) * time_step * u_flux_boundary_values (i+1,j) * t_boundary(i+1,j) ! assume no flux bc for temp ! flux_enter (i,j,2) = G%dyh(i,j) * time_step * CS%u_shelf(i,j)*t_boundary (i+1,j) endif -! if ((i .eq. is) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i-1,j) .eq. 1)) then +! if ((i == is) .AND. (hmask(i,j) == 0) .AND. (hmask(i-1,j) == 1)) then ! this is solely for the purposes of keeping the mask consistent while advancing ! the front without having to call pass_var - if cell is empty and cell to left ! is ice-covered then this cell will become partly covered ! hmask(i,j) = 2 -! elseif ((i .eq. ie) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i+1,j) .eq. 1)) then +! elseif ((i == ie) .AND. (hmask(i,j) == 0) .AND. (hmask(i+1,j) == 1)) then ! this is solely for the purposes of keeping the mask consistent while advancing ! the front without having to call pass_var - if cell is empty and cell to left ! is ice-covered then this cell will become partly covered @@ -6507,18 +6507,18 @@ subroutine ice_shelf_advect_temp_y (CS, time_step, h_after_uflux, h_after_vflux, if (((j+j_off) <= G%domain%njglobal+G%domain%njhalo) .AND. & ((j+j_off) >= G%domain%njhalo+1)) then - if (j+j_off .eq. G%domain%njhalo+1) then + if (j+j_off == G%domain%njhalo+1) then at_south_bdry=.true. else at_south_bdry=.false. endif - if (j+j_off .eq. G%domain%njglobal+G%domain%njhalo) then + if (j+j_off == G%domain%njglobal+G%domain%njhalo) then at_north_bdry=.true. else at_north_bdry=.false. endif - if (hmask(i,j) .eq. 1) then + if (hmask(i,j) == 1) then dxh = G%dxT(i,j) ; dyh = G%dyT(i,j) ; dxdyh = G%areaT(i,j) h_after_vflux (i,j) = h_after_uflux (i,j) @@ -6527,7 +6527,7 @@ subroutine ice_shelf_advect_temp_y (CS, time_step, h_after_uflux, h_after_vflux, ! 1ST DO south FACE - if (v_face_mask (i,j-1) .eq. 4.) then + if (v_face_mask (i,j-1) == 4.) then flux_diff_cell = flux_diff_cell + dxh * time_step * v_flux_boundary_values (i,j-1) * & t_boundary(i,j-1)/ dxdyh @@ -6543,11 +6543,11 @@ subroutine ice_shelf_advect_temp_y (CS, time_step, h_after_uflux, h_after_vflux, ! i may not cover all the cases.. but i cover the realistic ones - if (at_south_bdry .AND. (hmask(i,j-1).eq.3)) then ! at western bdry but there is a + if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then ! at western bdry but there is a ! thickness bdry condition, and the stencil contains it flux_diff_cell = flux_diff_cell + ABS(v_face) * dxh * time_step * stencil(-1) / dxdyh - elseif (hmask(i,j-1) * hmask(i,j-2) .eq. 1) then ! h(j-2) and h(j-1) are valid + elseif (hmask(i,j-1) * hmask(i,j-2) == 1) then ! h(j-2) and h(j-1) are valid phi = slope_limiter (stencil(-1)-stencil(-2), stencil(0)-stencil(-1)) flux_diff_cell = flux_diff_cell + ABS(v_face) * dxh * time_step / dxdyh * & @@ -6561,14 +6561,14 @@ subroutine ice_shelf_advect_temp_y (CS, time_step, h_after_uflux, h_after_vflux, elseif (v_face < 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available - if (hmask(i,j-1) * hmask(i,j+1) .eq. 1) then ! h(j-1) and h(j+1) are both valid + if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid phi = slope_limiter (stencil(0)-stencil(1), stencil(-1)-stencil(0)) flux_diff_cell = flux_diff_cell - ABS(v_face) * dxh * time_step / dxdyh * & (stencil(0) - phi * (stencil(0)-stencil(-1))/2) else flux_diff_cell = flux_diff_cell - ABS(v_face) * dxh * time_step / dxdyh * stencil(0) - if ((hmask(i,j-1) .eq. 0) .OR. (hmask(i,j-1) .eq. 2)) then + if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2)) then flux_enter(i,j-1,4) = ABS(v_face) * dyh * time_step * stencil(0) endif @@ -6580,7 +6580,7 @@ subroutine ice_shelf_advect_temp_y (CS, time_step, h_after_uflux, h_after_vflux, ! NEXT DO north FACE - if (v_face_mask(i,j+1) .eq. 4.) then + if (v_face_mask(i,j+1) == 4.) then flux_diff_cell = flux_diff_cell + dxh * time_step * v_flux_boundary_values (i,j+1) *& t_boundary(i,j+1)/ dxdyh @@ -6594,10 +6594,10 @@ subroutine ice_shelf_advect_temp_y (CS, time_step, h_after_uflux, h_after_vflux, if (v_face < 0) then !flux is into cell - we need info from h(j+2), h(j+1) if available - if (at_north_bdry .AND. (hmask(i,j+1).eq.3)) then ! at eastern bdry but there is a + if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then ! at eastern bdry but there is a ! thickness bdry condition, and the stencil contains it flux_diff_cell = flux_diff_cell + ABS(v_face) * dxh * time_step * stencil(1) / dxdyh - elseif (hmask(i,j+1) * hmask(i,j+2) .eq. 1) then ! h(j+2) and h(j+1) are valid + elseif (hmask(i,j+1) * hmask(i,j+2) == 1) then ! h(j+2) and h(j+1) are valid phi = slope_limiter (stencil(1)-stencil(2), stencil(0)-stencil(1)) flux_diff_cell = flux_diff_cell + ABS(v_face) * dxh * time_step / dxdyh * & (stencil(1) - phi * (stencil(1)-stencil(0))/2) @@ -6609,7 +6609,7 @@ subroutine ice_shelf_advect_temp_y (CS, time_step, h_after_uflux, h_after_vflux, elseif (v_face > 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available - if (hmask(i,j-1) * hmask(i,j+1) .eq. 1) then ! h(j-1) and h(j+1) are both valid + if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid phi = slope_limiter (stencil(0)-stencil(-1), stencil(1)-stencil(0)) flux_diff_cell = flux_diff_cell - ABS(v_face) * dxh * time_step / dxdyh * & (stencil(0) - phi * (stencil(0)-stencil(1))/2) @@ -6617,7 +6617,7 @@ subroutine ice_shelf_advect_temp_y (CS, time_step, h_after_uflux, h_after_vflux, ! (o.w. flux would most likely be out of cell) ! but h(j+2) is not flux_diff_cell = flux_diff_cell - ABS(v_face) * dxh * time_step / dxdyh * stencil(0) - if ((hmask(i,j+1) .eq. 0) .OR. (hmask(i,j+1) .eq. 2)) then + if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2)) then flux_enter(i,j+1,3) = ABS(v_face) * dxh * time_step * stencil(0) endif endif @@ -6628,35 +6628,35 @@ subroutine ice_shelf_advect_temp_y (CS, time_step, h_after_uflux, h_after_vflux, h_after_vflux (i,j) = h_after_vflux (i,j) + flux_diff_cell - elseif ((hmask(i,j) .eq. 0) .OR. (hmask(i,j) .eq. 2)) then + elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then - if (at_south_bdry .AND. (hmask(i,j-1) .EQ. 3)) then + if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then v_face = 0.5 * (CS%v_shelf(i-1,j-1) + CS%v_shelf(i,j-1)) flux_enter (i,j,3) = ABS(v_face) * G%dxT(i,j) * time_step * t_boundary(i,j-1)* & CS%thickness_boundary_values(i,j-1) - elseif (v_face_mask(i,j-1) .eq. 4.) then + elseif (v_face_mask(i,j-1) == 4.) then flux_enter (i,j,3) = G%dxT(i,j) * time_step * v_flux_boundary_values (i,j-1)*t_boundary(i,j-1) ! assume no flux bc for temp ! flux_enter (i,j,3) = G%dxh(i,j) * time_step * CS%v_shelf(i,j)*t_boundary (i,j-1) endif - if (at_north_bdry .AND. (hmask(i,j+1) .EQ. 3)) then + if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then v_face = 0.5 * (CS%v_shelf(i-1,j) + CS%v_shelf(i,j)) flux_enter (i,j,4) = ABS(v_face) * G%dxT(i,j) * time_step * t_boundary(i,j+1)* & CS%thickness_boundary_values(i,j+1) - elseif (v_face_mask(i,j+1) .eq. 4.) then + elseif (v_face_mask(i,j+1) == 4.) then flux_enter (i,j,4) = G%dxT(i,j) * time_step * v_flux_boundary_values (i,j+1)*t_boundary(i,j+1) ! assume no flux bc for temp ! flux_enter (i,j,4) = G%dxh(i,j) * time_step * CS%v_shelf(i,j)*t_boundary (i,j+1) endif -! if ((j .eq. js) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i,j-1) .eq. 1)) then +! if ((j == js) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j-1) == 1)) then ! this is solely for the purposes of keeping the mask consistent while advancing ! the front without having to call pass_var - if cell is empty and cell to left ! is ice-covered then this cell will become partly covered ! hmask (i,j) = 2 - ! elseif ((j .eq. je) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i,j+1) .eq. 1)) then + ! elseif ((j == je) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j+1) == 1)) then ! this is solely for the purposes of keeping the mask consistent while advancing the ! front without having to call pass_var - if cell is empty and cell to left is ! ice-covered then this cell will become partly covered diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 9ce2f37032..01202a013b 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -123,7 +123,7 @@ subroutine initialize_ice_thickness_from_file (h_shelf, area_shelf_h, hmask, G, if (area_shelf_h (i,j) >= G%areaT(i,j)) then hmask(i,j) = 1. - elseif (area_shelf_h (i,j) .eq. 0.0) then + elseif (area_shelf_h (i,j) == 0.0) then hmask(i,j) = 0. elseif ((area_shelf_h(i,j) > 0) .and. (area_shelf_h(i,j) <= G%areaT(i,j))) then hmask(i,j) = 2. @@ -209,7 +209,7 @@ subroutine initialize_ice_thickness_channel (h_shelf, area_shelf_h, hmask, G, PF endif endif - if ((i+G%idg_offset) .eq. G%domain%nihalo+1) then + if ((i+G%idg_offset) == G%domain%nihalo+1) then hmask(i-1,j) = 3.0 endif @@ -311,7 +311,7 @@ end subroutine initialize_ice_thickness_channel ! ! upstream boundary - set either dirichlet or flux condition -! if ((i+G%idg_offset) .eq. G%domain%nihalo+1) then +! if ((i+G%idg_offset) == G%domain%nihalo+1) then ! if (flux_bdry) then ! u_face_mask_boundary (i-1,j) = 4.0 ! u_flux_boundary_values (i-1,j) = input_flux @@ -328,14 +328,14 @@ end subroutine initialize_ice_thickness_channel ! ! side boundaries: no flow -! if (G%jdg_offset+j .eq. gjsc+1) then !bot boundary -! if (len_stress .eq. 0. .OR. G%geoLonCv(i,j-1) <= len_stress) then +! if (G%jdg_offset+j == gjsc+1) then !bot boundary +! if (len_stress == 0. .OR. G%geoLonCv(i,j-1) <= len_stress) then ! v_face_mask_boundary (i,j-1) = 0. ! else ! v_face_mask_boundary (i,j-1) = 1. ! endif -! elseif (G%jdg_offset+j .eq. gjec) then !top boundary -! if (len_stress .eq. 0. .OR. G%geoLonCv(i,j-1) <= len_stress) then +! elseif (G%jdg_offset+j == gjec) then !top boundary +! if (len_stress == 0. .OR. G%geoLonCv(i,j-1) <= len_stress) then ! v_face_mask_boundary (i,j) = 0. ! else ! v_face_mask_boundary (i,j) = 1. @@ -344,7 +344,7 @@ end subroutine initialize_ice_thickness_channel ! ! downstream boundary - CFBC -! if (i+G%idg_offset .eq. giec) then +! if (i+G%idg_offset == giec) then ! u_face_mask_boundary(i,j) = 2.0 ! endif diff --git a/src/ice_shelf/shelf_triangular_FEstuff.F90 b/src/ice_shelf/shelf_triangular_FEstuff.F90 index 6829774386..088ada5507 100644 --- a/src/ice_shelf/shelf_triangular_FEstuff.F90 +++ b/src/ice_shelf/shelf_triangular_FEstuff.F90 @@ -192,12 +192,12 @@ subroutine matrix_diagonal_triangle (CS, u_diagonal, v_diagonal) nu_lower => CS%ice_visc_lower_tri ; nu_upper => CS%ice_visc_upper_tri beta_lower => CS%taub_beta_eff_lower_tri ; beta_upper => CS%taub_beta_eff_upper_tri - do i=isc-1,iec+1 ; do j=jsc-1,jec+1 ; if (hmask(i,j) .eq. 1) then + do i=isc-1,iec+1 ; do j=jsc-1,jec+1 ; if (hmask(i,j) == 1) then dxh = G%dxT(i,j) dyh = G%dyT(i,j) dxdyh = G%areaT(i,j) - if (umask (i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node + if (umask (i,j-1) == 1) then ! this (bot right) is a degree of freedom node ux = 1./dxh ; uy = 0./dyh vx = 0. ; vy = 0. @@ -237,7 +237,7 @@ subroutine matrix_diagonal_triangle (CS, u_diagonal, v_diagonal) endif - if (umask (i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node + if (umask (i-1,j) == 1) then ! this (top left) is a degree of freedom node ux = 0./dxh ; uy = 1./dyh vx = 0. ; vy = 0. @@ -277,7 +277,7 @@ subroutine matrix_diagonal_triangle (CS, u_diagonal, v_diagonal) endif - if (umask (i-1,j-1) .eq. 1) then ! this (bot left) is a degree of freedom node + if (umask (i-1,j-1) == 1) then ! this (bot left) is a degree of freedom node ux = -1./dxh ; uy = -1./dyh vx = 0. ; vy = 0. @@ -298,7 +298,7 @@ subroutine matrix_diagonal_triangle (CS, u_diagonal, v_diagonal) beta_lower(i,j) * dxdyh * 1./24 endif - if (umask (i,j) .eq. 1) then ! this (top right) is a degree of freedom node + if (umask (i,j) == 1) then ! this (top right) is a degree of freedom node ux = 1./ dxh ; uy = 1./dyh vx = 0. ; vy = 0. @@ -360,9 +360,9 @@ end subroutine matrix_diagonal_triangle !~ domain_width = CS%len_lat - !~ do i=isc-1,iec+1 ; do j=jsc-1,jec+1 ; if (hmask(i,j) .eq. 1) then + !~ do i=isc-1,iec+1 ; do j=jsc-1,jec+1 ; if (hmask(i,j) == 1) then - !~ if ((umask(i-1,j-1) .eq. 3) .OR. (umask(i,j-1) .eq. 3) .OR. (umask(i-1,j) .eq. 3)) then + !~ if ((umask(i-1,j-1) == 3) .OR. (umask(i,j-1) == 3) .OR. (umask(i-1,j) == 3)) then !~ dxh = G%dxh(i,j) !~ dyh = G%dyh(i,j) @@ -373,7 +373,7 @@ end subroutine matrix_diagonal_triangle !~ uy = (u_boundary_values(i-1,j)-u_boundary_values(i-1,j-1))/dyh !~ vy = (v_boundary_values(i-1,j)-v_boundary_values(i-1,j-1))/dyh - !~ if (umask (i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node + !~ if (umask (i,j-1) == 1) then ! this (bot right) is a degree of freedom node !~ u_boundary_contr (i,j-1) = u_boundary_contr (i,j-1) + & !~ .5 * dxdyh * nu_lower (i,j) * ((4*ux+2*vy) * (1./dxh) + (uy+vy) * (0./dyh)) @@ -390,7 +390,7 @@ end subroutine matrix_diagonal_triangle !~ v_boundary_values(i-1,j) + v_boundary_values(i,j-1)) !~ endif - !~ if (umask (i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node + !~ if (umask (i-1,j) == 1) then ! this (top left) is a degree of freedom node !~ u_boundary_contr (i-1,j) = u_boundary_contr (i-1,j) + & !~ .5 * dxdyh * nu_lower (i,j) * ((4*ux+2*vy) * (0./dxh) + (uy+vy) * (1./dyh)) @@ -407,7 +407,7 @@ end subroutine matrix_diagonal_triangle !~ v_boundary_values(i-1,j) + v_boundary_values(i,j-1)) !~ endif - !~ if (umask (i-1,j-1) .eq. 1) then ! this (bot left) is a degree of freedom node + !~ if (umask (i-1,j-1) == 1) then ! this (bot left) is a degree of freedom node !~ u_boundary_contr (i-1,j-1) = u_boundary_contr (i-1,j-1) + & !~ .5 * dxdyh * nu_upper (i,j) * ((4*ux+2*vy) * (-1./dxh) + (uy+vy) * (-1./dyh)) @@ -426,7 +426,7 @@ end subroutine matrix_diagonal_triangle !~ endif - !~ if ((umask(i,j) .eq. 3) .OR. (umask(i,j-1) .eq. 3) .OR. (umask(i-1,j) .eq. 3)) then + !~ if ((umask(i,j) == 3) .OR. (umask(i,j-1) == 3) .OR. (umask(i-1,j) == 3)) then !~ dxh = G%dxh(i,j) !~ dyh = G%dyh(i,j) @@ -437,7 +437,7 @@ end subroutine matrix_diagonal_triangle !~ uy = (u_boundary_values(i,j)-u_boundary_values(i,j-1))/dyh !~ vy = (v_boundary_values(i,j)-v_boundary_values(i,j-1))/dyh - !~ if (umask (i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node + !~ if (umask (i,j-1) == 1) then ! this (bot right) is a degree of freedom node !~ u_boundary_contr (i,j-1) = u_boundary_contr (i,j-1) + & !~ .5 * dxdyh * nu_upper (i,j) * ((4*ux+2*vy) * (0./dxh) + (uy+vy) * (-1./dyh)) @@ -456,7 +456,7 @@ end subroutine matrix_diagonal_triangle !~ u_boundary_values(i,j-1)) !~ endif - !~ if (umask (i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node + !~ if (umask (i-1,j) == 1) then ! this (top left) is a degree of freedom node !~ u_boundary_contr (i-1,j) = u_boundary_contr (i-1,j) + & !~ .5 * dxdyh * nu_upper (i,j) * ((4*ux+2*vy) * (-1./dxh) + (uy+vy) * (0./dyh)) @@ -475,7 +475,7 @@ end subroutine matrix_diagonal_triangle !~ u_boundary_values(i,j-1)) !~ endif - !~ if (umask (i,j) .eq. 1) then ! this (top right) is a degree of freedom node + !~ if (umask (i,j) == 1) then ! this (top right) is a degree of freedom node !~ u_boundary_contr (i,j) = u_boundary_contr (i,j) + & !~ .5 * dxdyh * nu_upper (i,j) * ((4*ux+2*vy) * (1./dxh) + (uy+vy) * (1./dyh)) @@ -551,7 +551,7 @@ end subroutine matrix_diagonal_triangle !~ dyh = G%dyh(i,j) !~ dxdyh = G%dxdyh(i,j) - !~ if (hmask (i,j) .eq. 1) then + !~ if (hmask (i,j) == 1) then !~ ux = (u(i,j-1)-u(i-1,j-1)) / dxh !~ vx = (v(i,j-1)-v(i-1,j-1)) / dxh !~ uy = (u(i-1,j)-u(i-1,j-1)) / dyh @@ -605,14 +605,14 @@ end subroutine matrix_diagonal_triangle !~ do i=is,ie !~ do j=js,je - !~ if (hmask(i,j) .eq. 1) then ! this cell's vertices contain degrees of freedom + !~ if (hmask(i,j) == 1) then ! this cell's vertices contain degrees of freedom !~ ux = (u(i,j-1)-u(i-1,j-1))/dxh(i,j) !~ vx = (v(i,j-1)-v(i-1,j-1))/dxh(i,j) !~ uy = (u(i-1,j)-u(i-1,j-1))/dyh(i,j) !~ vy = (v(i-1,j)-v(i-1,j-1))/dyh(i,j) - !~ if (umask(i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node + !~ if (umask(i,j-1) == 1) then ! this (bot right) is a degree of freedom node !~ uret(i,j-1) = uret(i,j-1) + & !~ .5 * dxdyh(i,j) * nu_lower (i,j) * ((4*ux+2*vy) * (1./dxh(i,j)) + (uy+vy) * (0./dyh(i,j))) @@ -629,7 +629,7 @@ end subroutine matrix_diagonal_triangle !~ v(i-1,j) + v(i,j-1)) !~ endif - !~ if (umask(i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node + !~ if (umask(i-1,j) == 1) then ! this (top left) is a degree of freedom node !~ uret(i-1,j) = uret(i-1,j) + & !~ .5 * dxdyh(i,j) * nu_lower (i,j) * ((4*ux+2*vy) * (0./dxh(i,j)) + (uy+vy) * (1./dyh(i,j))) @@ -646,7 +646,7 @@ end subroutine matrix_diagonal_triangle !~ v(i-1,j) + v(i,j-1)) !~ endif - !~ if (umask(i-1,j-1) .eq. 1) then ! this (bot left) is a degree of freedom node + !~ if (umask(i-1,j-1) == 1) then ! this (bot left) is a degree of freedom node !~ uret(i-1,j-1) = uret(i-1,j-1) + & !~ .5 * dxdyh(i,j) * nu_upper (i,j) * ((4*ux+2*vy) * (-1./dxh(i,j)) + (uy+vy) * (-1./dyh(i,j))) @@ -669,7 +669,7 @@ end subroutine matrix_diagonal_triangle !~ uy = (u(i,j)-u(i,j-1))/dyh(i,j) !~ vy = (v(i,j)-v(i,j-1))/dyh(i,j) - !~ if (umask(i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node + !~ if (umask(i,j-1) == 1) then ! this (bot right) is a degree of freedom node !~ uret(i,j-1) = uret(i,j-1) + & !~ .5 * dxdyh(i,j) * nu_upper (i,j) * ((4*ux+2*vy) * (0./dxh(i,j)) + (uy+vy) * (-1./dyh(i,j))) @@ -686,7 +686,7 @@ end subroutine matrix_diagonal_triangle !~ u(i-1,j) + u(i,j-1)) !~ endif - !~ if (umask(i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node + !~ if (umask(i-1,j) == 1) then ! this (top left) is a degree of freedom node !~ uret(i-1,j) = uret(i-1,j) + & !~ .5 * dxdyh(i,j) * nu_upper (i,j) * ((4*ux+2*vy) * (-1./dxh(i,j)) + (uy+vy) * (0./dyh(i,j))) @@ -703,7 +703,7 @@ end subroutine matrix_diagonal_triangle !~ u(i-1,j) + u(i,j-1)) !~ endif - !~ if (umask(i,j) .eq. 1) then ! this (top right) is a degree of freedom node + !~ if (umask(i,j) == 1) then ! this (top right) is a degree of freedom node !~ uret(i,j) = uret(i,j) + & !~ .5 * dxdyh(i,j) * nu_upper (i,j) * ((4*ux+2*vy) * (1./dxh(i,j)) + (uy+vy) * (1./dyh(i,j))) diff --git a/src/ice_shelf/user_shelf_init.F90 b/src/ice_shelf/user_shelf_init.F90 index abe44044ea..650d4a9e5f 100644 --- a/src/ice_shelf/user_shelf_init.F90 +++ b/src/ice_shelf/user_shelf_init.F90 @@ -232,7 +232,7 @@ subroutine USER_update_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, G, C endif ; endif ; endif - if ((i+G%idg_offset) .eq. G%domain%nihalo+1) then + if ((i+G%idg_offset) == G%domain%nihalo+1) then hmask(i-1,j) = 3.0 endif diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 7aff08540a..0275bfc205 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -241,7 +241,7 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF) call log_param(PF, mdl, "!MAXIMUM_DEPTH", max_depth, & "The (diagnosed) maximum depth of the ocean.", units="m") endif - if (trim(config) .ne. "DOME") then + if (trim(config) /= "DOME") then call limit_topography(D, G, PF, max_depth) endif diff --git a/src/initialization/midas_vertmap.F90 b/src/initialization/midas_vertmap.F90 index 7225217b99..d8c30b345c 100644 --- a/src/initialization/midas_vertmap.F90 +++ b/src/initialization/midas_vertmap.F90 @@ -262,7 +262,7 @@ function tracer_z_init(tr_in,z_edges,e,nkml,nkbl,land_fill,wet,nlay,nlevs,debug, do j=1,ny i_loop: do i=1,nx - if (nlevs_data(i,j) .eq. 0 .or. wet(i,j) .eq. 0.) then + if (nlevs_data(i,j) == 0 .or. wet(i,j) == 0.) then tr(i,j,:) = land_fill cycle i_loop endif @@ -297,7 +297,7 @@ function tracer_z_init(tr_in,z_edges,e,nkml,nkbl,land_fill,wet,nlay,nlevs,debug, if (debug_) then if (PRESENT(i_debug)) then - if (i.eq.i_debug.and.j.eq.j_debug) then + if (i == i_debug.and.j == j_debug) then print *,'0001 k,k_top,k_bot,sum(wt),sum(z2-z1) = ',k,k_top,k_bot,sum(wt),sum(z2-z1) endif endif @@ -321,7 +321,7 @@ function tracer_z_init(tr_in,z_edges,e,nkml,nkbl,land_fill,wet,nlay,nlevs,debug, ! endif if (debug_) then if (PRESENT(i_debug)) then - if (i.eq.i_debug.and.j.eq.j_debug) then + if (i == i_debug.and.j == j_debug) then print *,'0002 k,k_top,k_bot,k_bot_prev,sl_tr = ',k,k_top,k_bot,k_bot_prev,sl_tr endif endif @@ -333,7 +333,7 @@ function tracer_z_init(tr_in,z_edges,e,nkml,nkbl,land_fill,wet,nlay,nlevs,debug, if (debug_) then if (PRESENT(i_debug)) then - if (i.eq.i_debug.and.j.eq.j_debug) then + if (i == i_debug.and.j == j_debug) then print *,'0003 k,tr = ',k,tr(i,j,k) endif endif @@ -357,7 +357,7 @@ function tracer_z_init(tr_in,z_edges,e,nkml,nkbl,land_fill,wet,nlay,nlevs,debug, if (debug_) then if (PRESENT(i_debug)) then - if (i.eq.i_debug.and.j.eq.j_debug) then + if (i == i_debug.and.j == j_debug) then print *,'0004 k,kz,nlevs,sl_tr,tr = ',k,kz,nlevs_data(i,j),sl_tr,tr(i,j,k) print *,'0005 k,kz,tr(kz),tr(kz-1),tr(kz+1) = ',k,kz,tr_1d(kz),tr_1d(kz-1),tr_1d(kz+1),z_edges(kz+2) endif @@ -792,7 +792,7 @@ function find_interfaces(rho,zin,Rb,depth,nlevs,nkml,nkbl,hml,debug) result(zi) if (dir == 1) then do k=2,nlevs_data(i,j)-1 if (rho_(i,k) - rho_(i,k-1) < 0.0 ) then - if (k.eq.2) then + if (k == 2) then rho_(i,k-1)=rho_(i,k)-epsln else drhodz = (rho_(i,k+1)-rho_(i,k-1))/(zin(k+1)-zin(k-1)) @@ -807,7 +807,7 @@ function find_interfaces(rho,zin,Rb,depth,nlevs,nkml,nkbl,hml,debug) result(zi) else do k=nlevs_data(i,j)-1,2,-1 if (rho_(i,k+1) - rho_(i,k) < 0.0) then - if (k .eq. nlevs_data(i,j)-1) then + if (k == nlevs_data(i,j)-1) then rho_(i,k+1)=rho_(i,k-1)+epsln else drhodz = (rho_(i,k+1)-rho_(i,k-1))/(zin(k+1)-zin(k-1)) @@ -922,7 +922,7 @@ subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) do j=1,nj do i=1,ni - if (fill(i,j) .eq. 1) then + if (fill(i,j) == 1) then B(i,j,1)=1-nm(i+1,j);B(i,j,2)=1-nm(i-1,j) B(i,j,3)=1-nm(i,j+1);B(i,j,4)=1-nm(i,j-1) endif @@ -932,7 +932,7 @@ subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) do n=1,niter do j=1,nj do i=1,ni - if (fill(i,j) .eq. 1) then + if (fill(i,j) == 1) then bsum = real(B(i,j,1)+B(i,j,2)+B(i,j,3)+B(i,j,4)) Isum = 1.0/bsum res(i,j)=Isum*(B(i,j,1)*mp(i+1,j)+B(i,j,2)*mp(i-1,j)+& diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 37c46407af..2672308fd7 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -310,7 +310,7 @@ subroutine init_oda(Time, G, GV, CS) T_grid%mask(i,j,k) = 1.0 end if end do; end do - if (k .eq. 1) then + if (k == 1) then T_grid%z(:,:,k) = global2D/2 else T_grid%z(:,:,k) = T_grid%z(:,:,k-1) + (global2D + global2D_old)/2 @@ -525,9 +525,9 @@ subroutine set_analysis_time(Time,CS) CS%Time=increment_time(CS%Time,CS%assim_frequency*3600) call get_date(Time, yr, mon, day, hr, min, sec) - if (pe() .eq. mpp_root_pe()) print *, 'Model Time: ', yr, mon, day, hr, min, sec + if (pe() == mpp_root_pe()) print *, 'Model Time: ', yr, mon, day, hr, min, sec call get_date(CS%time, yr, mon, day, hr, min, sec) - if (pe() .eq. mpp_root_pe()) print *, 'Assimilation Time: ', yr, mon, day, hr, min, sec + if (pe() == mpp_root_pe()) print *, 'Assimilation Time: ', yr, mon, day, hr, min, sec endif if (CS%Time < Time) then call MOM_error(FATAL, " set_analysis_time: " // & diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index c133df5d4e..7fcb842bee 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -440,8 +440,8 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & Ub(i,j,fr,m) = CS%wave_structure_CSp%Uavg_profile(i,j,nzm) Umax(i,j,fr,m) = maxval(CS%wave_structure_CSp%Uavg_profile(i,j,1:nzm)) !! for debugging print profile, etc. Delete later - !if (id_g .eq. 260 .and. & - ! jd_g .eq. 50 .and. & + !if (id_g == 260 .and. & + ! jd_g == 50 .and. & ! tot_En_mode(i,j,1,1)>500.0) then ! print *, 'Profiles for mode ',m,' and frequency ',fr ! print *, 'id_g=', id_g, 'jd_g=', jd_g @@ -661,7 +661,7 @@ subroutine sum_En(G, CS, En, label) En_sum = En_sum + tmpForSumming enddo En_sum_diff = En_sum - CS%En_sum - if (CS%En_sum .ne. 0.0) then + if (CS%En_sum /= 0.0) then En_sum_pdiff= (En_sum_diff/CS%En_sum)*100.0 else En_sum_pdiff= 0.0; @@ -1790,7 +1790,7 @@ subroutine reflect(En, NAngle, CS, G, LB) id_g = i + G%idg_offset ! redistribute energy in angular space if ray will hit boundary ! i.e., if energy is in a reflecting cell - if (angle_c(i,j) .ne. CS%nullangle) then + if (angle_c(i,j) /= CS%nullangle) then do a=1,NAngle if (En(i,j,a) > 0.0) then ! if ray is incident, keep specified boundary angle @@ -1818,7 +1818,7 @@ subroutine reflect(En, NAngle, CS, G, LB) endif a_r = nint(angle_r/Angle_size) + 1 do while (a_r > Nangle) ; a_r = a_r - Nangle ; enddo - if (a .ne. a_r) then + if (a /= a_r) then En_reflected(a_r) = part_refl(i,j)*En(i,j,a) En(i,j,a) = (1.0-part_refl(i,j))*En(i,j,a) endif @@ -2536,7 +2536,7 @@ subroutine internal_tides_init(Time, G, GV, param_file, diag, CS) do j=jsd,jed do i=isd,ied ! flag cells with partial reflection - if (CS%refl_angle(i,j) .ne. CS%nullangle .and. & + if (CS%refl_angle(i,j) /= CS%nullangle .and. & CS%refl_pref(i,j) < 1.0 .and. CS%refl_pref(i,j) > 0.0) then CS%refl_pref_logical(i,j) = .true. endif diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 61555090ab..532362f082 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -934,10 +934,10 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) "used which introduced potential restart issues. This flag will be \n"//& "deprecated in a future release.", default=.false.) if (CS%interpolate_Res_fn) then - if (CS%Res_coef_visc .ne. CS%Res_coef_khth) call MOM_error(FATAL, & + if (CS%Res_coef_visc /= CS%Res_coef_khth) call MOM_error(FATAL, & "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//& "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_SCALE_COEF.") - if (CS%Res_fn_power_visc .ne. CS%Res_fn_power_khth) call MOM_error(FATAL, & + if (CS%Res_fn_power_visc /= CS%Res_fn_power_khth) call MOM_error(FATAL, & "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//& "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_FN_POWER.") endif diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index f309fff485..08ed2fb130 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -303,7 +303,7 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves) '\t MatchBoth = match gradient for both diffusivity and NLT\n'// & '\t ParabolicNonLocal = sigma*(1-sigma)^2 for diffusivity; (1-sigma)^2 for NLT', & default='SimpleShapes') - if (CS%MatchTechnique.eq.'ParabolicNonLocal') then + if (CS%MatchTechnique == 'ParabolicNonLocal') then ! This forces Cs2 (Cs in non-local computation) to equal 1 for parabolic non-local option. ! May be used during CVMix initialization. Cs_is_one=.true. @@ -960,16 +960,16 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & ! LMD94 shape function, not matching is equivalent to matching to a zero diffusivity. !BGR/ Add option for use of surface buoyancy flux with total sw flux. - if (CS%SW_METHOD .eq. SW_METHOD_ALL_SW) then + if (CS%SW_METHOD == SW_METHOD_ALL_SW) then surfBuoyFlux = buoyFlux(i,j,1) - elseif (CS%SW_METHOD .eq. SW_METHOD_MXL_SW) then + elseif (CS%SW_METHOD == SW_METHOD_MXL_SW) then surfBuoyFlux = buoyFlux(i,j,1) - buoyFlux(i,j,int(kOBL)+1) ! We know the actual buoyancy flux into the OBL - elseif (CS%SW_METHOD .eq. SW_METHOD_LV1_SW) then + elseif (CS%SW_METHOD == SW_METHOD_LV1_SW) then surfBuoyFlux = buoyFlux(i,j,1) - buoyFlux(i,j,2) endif ! If option "MatchBoth" is selected in CVMix, MOM should be capable of matching. - if (.not. (CS%MatchTechnique.eq.'MatchBoth')) then + if (.not. (CS%MatchTechnique == 'MatchBoth')) then Kdiffusivity(:,:) = 0. ! Diffusivities for heat and salt (m2/s) Kviscosity(:) = 0. ! Viscosity (m2/s) else diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index 3fd52e456d..4a14bc41ba 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -609,12 +609,12 @@ subroutine opacity_init(Time, G, param_file, diag, tracer_flow, CS, optics) "The number of bands of penetrating shortwave radiation.", & default=1) if (CS%Opacity_scheme == DOUBLE_EXP ) then - if (optics%nbands.ne.2) then + if (optics%nbands /= 2) then call MOM_error(FATAL, "set_opacity: "// & "Cannot use a double_exp opacity scheme with nbands!=2.") endif elseif (CS%Opacity_scheme == SINGLE_EXP ) then - if (optics%nbands.ne.1) then + if (optics%nbands /= 1) then call MOM_error(FATAL, "set_opacity: "// & "Cannot use a single_exp opacity scheme with nbands!=1.") endif diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 9ecf1374ef..b786f9c919 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -267,14 +267,14 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, end select ! Check profile consistency - if (CS%use_CVMix_tidal .and. (CS%int_tide_profile.eq.STLAURENT_02 .or. & - CS%int_tide_profile.eq.POLZIN_09)) then + if (CS%use_CVMix_tidal .and. (CS%int_tide_profile == STLAURENT_02 .or. & + CS%int_tide_profile == POLZIN_09)) then call MOM_error(FATAL, "tidal_mixing_init: Tidal mixing profile"// & " "//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.eq.SIMMONS_04.or. & - CS%int_tide_profile.eq.SCHMITTNER)) then + else if (.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.") diff --git a/src/tracer/MOM_OCMIP2_CO2calc.F90 b/src/tracer/MOM_OCMIP2_CO2calc.F90 index b5187d5d1d..8c2809418d 100644 --- a/src/tracer/MOM_OCMIP2_CO2calc.F90 +++ b/src/tracer/MOM_OCMIP2_CO2calc.F90 @@ -336,7 +336,7 @@ subroutine MOM_ocmip2_co2calc(dope_vec, mask, & ! recommended (xacc of 10**-9 drops precision to 2 significant ! figures). ! - if (mask(i,j) .ne. 0.0) then !{ + if (mask(i,j) /= 0.0) then !{ htotal(i,j) = drtsafe(k0, k1, k2, kb, k1p, k2p, k3p, ksi, kw, & ks, kf, bt, dic_in(i,j), ft, pt_in(i,j),& sit_in(i,j), st, ta_in(i,j), & @@ -433,7 +433,7 @@ function drtsafe(k0, k1, k2, kb, k1p, k2p, k3p, ksi, kw, ks, kf, & dxold=dx dx=0.5*(xh-xl) drtsafe=xl+dx - if (xl .eq. drtsafe) then + if (xl == drtsafe) then ! write (6,*) 'Exiting drtsafe at A on iteration ', j, ', ph = ', -log10(drtsafe) return endif @@ -442,7 +442,7 @@ function drtsafe(k0, k1, k2, kb, k1p, k2p, k3p, ksi, kw, ks, kf, & dx=f/df temp=drtsafe drtsafe=drtsafe-dx - if (temp .eq. drtsafe) then + if (temp == drtsafe) then ! write (6,*) 'Exiting drtsafe at B on iteration ', j, ', ph = ', -log10(drtsafe) return endif diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 1e9f5577c3..6ec168b499 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -266,7 +266,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia g_tracer=>CS%g_tracer_list do - if (INDEX(CS%IC_file, '_NULL_') .ne. 0) then + if (INDEX(CS%IC_file, '_NULL_') /= 0) then call MOM_error(WARNING,"The name of the IC_file "//trim(CS%IC_file)//& " indicates no MOM initialization was asked for the generic tracers."//& "Bypassing the MOM initialization of ALL generic tracers!") @@ -293,7 +293,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia !Check/apply the bounds for each g_tracer do k=1,nk ; do j=jsc,jec ; do i=isc,iec - if (tr_ptr(i,j,k) .ne. CS%tracer_land_val) then + if (tr_ptr(i,j,k) /= CS%tracer_land_val) then if (tr_ptr(i,j,k) < g_tracer%src_var_valid_min) tr_ptr(i,j,k) = g_tracer%src_var_valid_min !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 @@ -301,9 +301,9 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia enddo; enddo ; enddo !jgj: Reset CASED to 0 below K=1 - if (trim(g_tracer_name) .eq. 'cased') then + if (trim(g_tracer_name) == 'cased') then do k=2,nk ; do j=jsc,jec ; do i=isc,iec - if (tr_ptr(i,j,k) .ne. CS%tracer_land_val) then + if (tr_ptr(i,j,k) /= CS%tracer_land_val) then tr_ptr(i,j,k) = 0.0 endif enddo; enddo ; enddo @@ -404,7 +404,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia if (associated(CS%g_diag_list)) then g_diag=>CS%g_diag_list do - if (g_diag%Z_diag .ne. 0) & + if (g_diag%Z_diag /= 0) & call register_Z_tracer(g_diag%field_ptr, trim(g_diag%name),g_diag%longname , g_diag%units, & day, G, diag_to_Z_CSp) diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 9e68b543de..f5ce1f10e6 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -1448,7 +1448,7 @@ subroutine neutral_surface_T_eval(nk, ns, k_sub, Ks, Ps, T_mean, T_int, deg, iMe ks_top = k_sub ks_bot = k_sub + 1 - if ( Ks(ks_top) .ne. Ks(ks_bot) ) then + if ( Ks(ks_top) /= Ks(ks_bot) ) then call MOM_error(FATAL, "Neutral surfaces span more than one layer") endif kl = Ks(k_sub) @@ -2231,7 +2231,7 @@ logical function test_rnp(expected_pos, test_pos, title) if (test_rnp) then write(stdunit,'(A, f20.16, " .neq. ", f20.16, " <-- WRONG")') title, expected_pos, test_pos else - write(stdunit,'(A, f20.16, " .eq. ", f20.16)') title, expected_pos, test_pos + write(stdunit,'(A, f20.16, " == ", f20.16)') title, expected_pos, test_pos endif end function test_rnp !> Deallocates neutral_diffusion control structure diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index efb0b6ceed..1389365139 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -566,7 +566,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & segment%tr_Reg%Tr(m)%tres(I,j,k)= (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(I,j,k) + & dt*(u_L_in*Tr(m)%t(I+ishift,j,k) - & u_L_out*segment%tr_Reg%Tr(m)%t(I,j,k))) -! if (j.eq.10 .and. segment%direction==OBC_DIRECTION_E .and. m==2 .and. k.eq.1) & +! if (j == 10 .and. segment%direction==OBC_DIRECTION_E .and. m==2 .and. k == 1) & ! print *,'tres=',segment%tr_Reg%Tr(m)%tres(I,j,k),& ! segment%tr_Reg%Tr(m)%t(I,j,k), fac1 endif diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index 4186c2d34d..0e9a18ffad 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -309,7 +309,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, param_file, end select ! Modify salinity and temperature when z coordinates are used - if ( coordinateMode(verticalCoordinate) .eq. REGRIDDING_ZSTAR ) then + if ( coordinateMode(verticalCoordinate) == REGRIDDING_ZSTAR ) then index_bay_z = Nint ( dome2d_depth_bay * G%ke ) do j = G%jsc,G%jec ; do i = G%isc,G%iec x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon @@ -321,7 +321,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, param_file, endif ! Z initial conditions ! Modify salinity and temperature when sigma coordinates are used - if ( coordinateMode(verticalCoordinate) .eq. REGRIDDING_SIGMA ) then + if ( coordinateMode(verticalCoordinate) == REGRIDDING_SIGMA ) then do i = G%isc,G%iec ; do j = G%jsc,G%jec x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon if ( x <= dome2d_width_bay ) then @@ -333,8 +333,8 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, param_file, ! Modify temperature when rho coordinates are used T(G%isc:G%iec,G%jsc:G%jec,1:G%ke) = 0.0 - if (( coordinateMode(verticalCoordinate) .eq. REGRIDDING_RHO ) .or. & - ( coordinateMode(verticalCoordinate) .eq. REGRIDDING_LAYER )) then + if (( coordinateMode(verticalCoordinate) == REGRIDDING_RHO ) .or. & + ( coordinateMode(verticalCoordinate) == REGRIDDING_LAYER )) then do i = G%isc,G%iec ; do j = G%jsc,G%jec x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon if ( x <= dome2d_width_bay ) then diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 3e6baf1f23..7d6d5644a9 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -279,7 +279,7 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg) Def_Rad = sqrt(D_edge*g_prime_tot) / (1.0e-4*1000.0) tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5e3*Def_Rad) * GV%m_to_H - if (OBC%number_of_segments .ne. 1) then + if (OBC%number_of_segments /= 1) then print *, 'Error in DOME OBC segment setup' return !!! Need a better error message here endif diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 5cc38f8f24..c183b5b7a3 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -688,7 +688,7 @@ subroutine Surface_Bands_by_data_override(day_center,G,GV,CS) varread1 = 'wavenumber' !Old method gives wavenumber varread2 = 'frequency' !New method gives frequency rcode_wn = NF90_OPEN(trim(SurfBandFileName), NF90_NOWRITE, ncid) - if (rcode_wn .ne. 0) then + if (rcode_wn /= 0) then call MOM_error(FATAL,"error opening file "//trim(SurfBandFileName)//& " in MOM_wave_interface.") endif @@ -696,49 +696,49 @@ subroutine Surface_Bands_by_data_override(day_center,G,GV,CS) rcode_wn = NF90_INQ_VARID(ncid, varread1, varid_wn) rcode_fr = NF90_INQ_VARID(ncid, varread2, varid_fr) - if (rcode_wn .ne. 0 .and. rcode_fr .ne. 0) then + if (rcode_wn /= 0 .and. rcode_fr /= 0) then call MOM_error(FATAL,"error finding variable "//trim(varread1)//& " or "//trim(varread2)//" in file "//trim(SurfBandFileName)//" in MOM_wave_interface.") - elseif (rcode_wn.eq.0) then + elseif (rcode_wn == 0) then ! wavenumbers found: PartitionMode=0 rcode_wn = NF90_INQUIRE_VARIABLE(ncid, varid_wn, ndims=ndims, & dimids=dims) - if (rcode_wn .ne. 0) then + if (rcode_wn /= 0) then call MOM_error(FATAL, & 'error inquiring dimensions MOM_wave_interface.') endif rcode_wn = NF90_INQUIRE_DIMENSION(ncid, dims(1), dim_name(1), len=id) - if (rcode_wn .ne. 0) then + if (rcode_wn /= 0) then call MOM_error(FATAL,"error reading dimension 1 data for "// & trim(varread1)//" in file "// trim(SurfBandFileName)// & " in MOM_wave_interface.") endif rcode_wn = NF90_INQ_VARID(ncid, dim_name(1), dim_id(1)) - if (rcode_wn .ne. 0) then + if (rcode_wn /= 0) then call MOM_error(FATAL,"error finding variable "//trim(dim_name(1))//& " in file "//trim(SurfBandFileName)//" in MOM_wave_interace.") endif ! Allocating size of wavenumber bins allocate( CS%WaveNum_Cen(1:id) ) ; CS%WaveNum_Cen(:)=0.0 - elseif (rcode_fr.eq.0) then + elseif (rcode_fr == 0) then ! frequencies found: PartitionMode=1 rcode_fr = NF90_INQUIRE_VARIABLE(ncid, varid_fr, ndims=ndims, & dimids=dims) - if (rcode_fr .ne. 0) then + if (rcode_fr /= 0) then call MOM_error(FATAL,& 'error inquiring dimensions MOM_wave_interface.') endif rcode_fr = NF90_INQUIRE_DIMENSION(ncid, dims(1), dim_name(1), len=id) - if (rcode_fr .ne. 0) then + if (rcode_fr /= 0) then call MOM_error(FATAL,"error reading dimension 1 data for "// & trim(varread2)//" in file "// trim(SurfBandFileName)// & " in MOM_wave_interface.") endif rcode_fr = NF90_INQ_VARID(ncid, dim_name(1), dim_id(1)) - if (rcode_fr .ne. 0) then + if (rcode_fr /= 0) then call MOM_error(FATAL,"error finding variable "//trim(dim_name(1))//& " in file "//trim(SurfBandFileName)//" in MOM_wave_interace.") endif @@ -758,7 +758,7 @@ subroutine Surface_Bands_by_data_override(day_center,G,GV,CS) start = 1; count = 1; count(1) = id if (PartitionMode==0) then rcode_wn = NF90_GET_VAR(ncid, dim_id(1), CS%WaveNum_Cen, start, count) - if (rcode_wn .ne. 0) then + if (rcode_wn /= 0) then call MOM_error(FATAL,& "error reading dimension 1 values for var_name "// & trim(varread1)//",dim_name "//trim(dim_name(1))// & @@ -767,7 +767,7 @@ subroutine Surface_Bands_by_data_override(day_center,G,GV,CS) NUMBANDS = ID elseif (PartitionMode==1) then rcode_fr = NF90_GET_VAR(ncid, dim_id(1), CS%Freq_Cen, start, count) - if (rcode_fr .ne. 0) then + if (rcode_fr /= 0) then call MOM_error(FATAL,& "error reading dimension 1 values for var_name "// & trim(varread2)//",dim_name "//trim(dim_name(1))// & @@ -1142,7 +1142,7 @@ subroutine StokesMixing(G, GV, DT, h, u, v, WAVES ) do k = 1, G%ke do j = G%jscB, G%jecB do i = G%iscB, G%iecB - if (k.eq.1) then + if (k == 1) then dTauUp = 0. dTauDn = 0.5*(WAVES%Kvs(i,j,k+1)+WAVES%Kvs(i+1,j,k+1))*& (waves%us_x(i,j,k)-waves%us_x(i,j,k+1))& @@ -1154,7 +1154,7 @@ subroutine StokesMixing(G, GV, DT, h, u, v, WAVES ) dTauDn = 0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i+1,j,k+1))*& (waves%us_x(i,j,k)-waves%us_x(i,j,k+1))& /(GV%H_to_m *0.5*(h(i,j,k)+h(i,j,k+1)) ) - elseif (k.eq.G%ke) then + elseif (k == G%ke) then dTauUp = 0.5*(waves%Kvs(i,j,k)+waves%Kvs(i+1,j,k))*& (waves%us_x(i,j,k-1)-waves%us_x(i,j,k))& /(GV%H_to_m *0.5*(h(i,j,k-1)+h(i,j,k)) ) @@ -1169,7 +1169,7 @@ subroutine StokesMixing(G, GV, DT, h, u, v, WAVES ) do k = 1, G%ke do j = G%jscB, G%jecB do i = G%iscB, G%iecB - if (k.eq.1) then + if (k == 1) then dTauUp = 0. dTauDn = 0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i,j+1,k+1))& *(waves%us_y(i,j,k)-waves%us_y(i,j,k+1))& @@ -1181,7 +1181,7 @@ subroutine StokesMixing(G, GV, DT, h, u, v, WAVES ) dTauDn = 0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i,j+1,k+1))*& (waves%us_y(i,j,k)-waves%us_y(i,j,k+1))& /(GV%H_to_m *0.5*(h(i,j,k)+h(i,j,k+1)) ) - elseif (k.eq.G%ke) then + elseif (k == G%ke) then dTauUp = 0.5*(waves%Kvs(i,j,k)+waves%Kvs(i,j+1,k))*& (waves%us_y(i,j,k-1)-waves%us_y(i,j,k))& /(GV%H_to_m *0.5*(h(i,j,k-1)+h(i,j,k)) ) diff --git a/src/user/adjustment_initialization.F90 b/src/user/adjustment_initialization.F90 index 68975fc41f..c73c8a12e4 100644 --- a/src/user/adjustment_initialization.F90 +++ b/src/user/adjustment_initialization.F90 @@ -110,7 +110,7 @@ subroutine adjustment_initialize_thickness ( h, G, GV, param_file, just_read_par select case ( coordinateMode(verticalCoordinate) ) case ( REGRIDDING_LAYER, REGRIDDING_RHO ) - if (delta_S_strat.ne.0.) then + if (delta_S_strat /= 0.) then adjustment_delta = adjustment_deltaS / delta_S_strat * G%max_depth do k=1,nz+1 e0(k) = adjustment_delta-(G%max_depth+2*adjustment_delta) * (real(k-1) / real(nz)) @@ -128,7 +128,7 @@ subroutine adjustment_initialize_thickness ( h, G, GV, param_file, just_read_par end do target_values = target_values - 1000. do j=js,je ; do i=is,ie - if (front_wave_length.ne.0.) then + if (front_wave_length /= 0.) then y = ( 0.125 + G%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) ) yy = 2. * ( G%geoLatT(i,j) - 0.5 * G%len_lat ) / adjustment_width yy = min(1.0, yy); yy = max(-1.0, yy) @@ -142,7 +142,7 @@ subroutine adjustment_initialize_thickness ( h, G, GV, param_file, just_read_par x = x * acos( 0. ) delta_S = adjustment_deltaS * 0.5 * (1. - sin( x ) ) do k=2,nz - if (dSdz.ne.0.) then + if (dSdz /= 0.) then eta1D(k) = ( target_values(k) - ( S_ref + delta_S ) ) / dSdz else eta1D(k) = e0(k) - (0.5*adjustment_delta) * sin( x ) @@ -258,7 +258,7 @@ subroutine adjustment_initialize_temperature_salinity ( T, S, h, G, GV, param_fi do k=nz,1,-1 eta1d(k) = eta1d(k+1) + h(i,j,k)*GV%H_to_m enddo - if (front_wave_length.ne.0.) then + if (front_wave_length /= 0.) then y = ( 0.125 + G%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) ) yy = 2. * ( G%geoLatT(i,j) - 0.5 * G%len_lat ) / front_wave_length yy = min(1.0, yy); yy = max(-1.0, yy) diff --git a/src/user/dumbbell_initialization.F90 b/src/user/dumbbell_initialization.F90 index 88b80e84c6..e2bc9b5869 100644 --- a/src/user/dumbbell_initialization.F90 +++ b/src/user/dumbbell_initialization.F90 @@ -344,7 +344,7 @@ subroutine dumbbell_initialize_sponges(G, GV, tv, param_file, use_ALE, CSp, ACSp S(i,j,k)=S_ref - 0.5*S_range enddo endif -! if (j.eq.G%jsc) print *,'i,Sponge S= ',i,S(i,1,1) +! if (j == G%jsc) print *,'i,Sponge S= ',i,S(i,1,1) enddo enddo diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90 index be64ef163e..06b1df3218 100644 --- a/src/user/sloshing_initialization.F90 +++ b/src/user/sloshing_initialization.F90 @@ -138,11 +138,11 @@ subroutine sloshing_initialize_thickness ( h, G, GV, param_file, just_read_param x = G%geoLonT(i,j) / G%len_lon displ(k) = a0 * cos(acos(-1.0)*x) + weight_z; - if ( k .EQ. 1 ) then + if ( k == 1 ) then displ(k) = 0.0 end if - if ( k .EQ. nz+1 ) then + if ( k == nz+1 ) then displ(k) = 0.0 end if