Skip to content

Commit

Permalink
Merge pull request #274 from russfiedler/optimise_rebinnig
Browse files Browse the repository at this point in the history 
Optimise rebinnig
  • Loading branch information
Marshall Ward committed Feb 26, 2019
2 parents 83a2feb + 6fb47f6 commit 9f5866b
Show file tree
Hide file tree
Showing 5 changed files with 263 additions and 178 deletions.
202 changes: 103 additions & 99 deletions src/mom5/ocean_diag/ocean_adv_vel_diag.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -840,71 +840,47 @@ subroutine transport_on_s(Time, Adv_vel)
! mass transports leaving faces of grid cells
wrk1=0.0

if (id_tx_trans > 0) then
if (id_tx_trans > 0 .or. id_tx_trans_int_z > 0) then
wrk1_2d(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Adv_vel%uhrho_et(i,j,k)*Grd%dyte(i,j)*transport_convert
wrk1_2d(i,j) = wrk1_2d(i,j) + wrk1(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_tx_trans, wrk1(:,:,:))
if (id_tx_trans > 0) call diagnose_3d(Time, Grd, id_tx_trans, wrk1(:,:,:))
if (id_tx_trans_int_z > 0) call diagnose_2d(Time, Grd, id_tx_trans_int_z, wrk1_2d(:,:))
endif
if (id_tx_trans_int_z > 0) then
wrk1_2d(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = wrk1_2d(i,j) + Adv_vel%uhrho_et(i,j,k)*Grd%dyte(i,j)*transport_convert
enddo
enddo
enddo
call diagnose_2d(Time, Grd, id_tx_trans_int_z, wrk1_2d(:,:))
endif

if (id_ty_trans > 0) then
if (id_ty_trans > 0 .or. id_ty_trans_int_z > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Adv_vel%vhrho_nt(i,j,k)*Grd%dxtn(i,j)*transport_convert
wrk1_2d(i,j) = wrk1_2d(i,j) + wrk1(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_ty_trans, wrk1(:,:,:))
if (id_ty_trans > 0) call diagnose_3d(Time, Grd, id_ty_trans, wrk1(:,:,:))
if (id_ty_trans_int_z > 0) call diagnose_2d(Time, Grd, id_ty_trans_int_z, wrk1_2d(:,:))
endif
if (id_ty_trans_int_z > 0) then
wrk1_2d(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = wrk1_2d(i,j) + Adv_vel%vhrho_nt(i,j,k)*Grd%dxtn(i,j)*transport_convert
enddo
enddo
enddo
call diagnose_2d(Time, Grd, id_ty_trans_int_z, wrk1_2d(:,:))
endif

if (id_tz_trans > 0) then
if (id_tz_trans > 0 .or. id_tz_trans_int_z > 0) then
wrk1_2d(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Adv_vel%wrho_bt(i,j,k)*Grd%dat(i,j)*transport_convert
wrk1_2d(i,j) = wrk1_2d(i,j) + wrk1(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_tz_trans, wrk1(:,:,:))
if (id_tz_trans > 0) call diagnose_3d(Time, Grd, id_tz_trans, wrk1(:,:,:))
if (id_tz_trans_int_z > 0) call diagnose_2d(TIme, Grd, id_tz_trans_int_z, wrk1_2d(:,:))
endif
if (id_tz_trans_int_z > 0) then
wrk1_2d(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = wrk1_2d(i,j) + Adv_vel%wrho_bt(i,j,k)*Grd%dat(i,j)*transport_convert
enddo
enddo
enddo
call diagnose_2d(TIme, Grd, id_tz_trans_int_z, wrk1_2d(:,:))
endif

if (id_tz_trans_sq > 0) then
do k=1,nk
do j=jsc,jec
Expand Down Expand Up @@ -1012,25 +988,25 @@ subroutine transport_on_nrho (Time, Dens, Adv_vel)
work1(:,:,:) = 0.0
work2(:,:,:) = 0.0

do k_rho=1,neutralrho_nk
do k_rho = 1,neutralrho_nk

tmp(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
do k = 1,nk
do j = jsc,jec
do i = isc,iec
if (k_rho == 1) then
if(Dens%neutralrho(i,j,k) < Dens%neutralrho_bounds(k_rho+1)) then
if (Dens%neutralrho(i,j,k) < Dens%neutralrho_bounds(k_rho+1)) then
tmp(1,i,j) = tmp(1,i,j) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,j) = tmp(2,i,j) + Adv_vel%vhrho_nt(i,j,k)
endif
elseif(k_rho < neutralrho_nk) then
if( (Dens%neutralrho_bounds(k_rho) <= Dens%neutralrho(i,j,k)) .and. &
elseif (k_rho < neutralrho_nk) then
if ((Dens%neutralrho_bounds(k_rho) <= Dens%neutralrho(i,j,k)) .and. &
(Dens%neutralrho(i,j,k) < Dens%neutralrho_bounds(k_rho+1)) ) then
tmp(1,i,j) = tmp(1,i,j) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,j) = tmp(2,i,j) + Adv_vel%vhrho_nt(i,j,k)
endif
else ! if (k_rho == neutralrho_nk) then
if(Dens%neutralrho_bounds(k_rho) <= Dens%neutralrho(i,j,k)) then
if (Dens%neutralrho_bounds(k_rho) <= Dens%neutralrho(i,j,k)) then
tmp(1,i,j) = tmp(1,i,j) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,j) = tmp(2,i,j) + Adv_vel%vhrho_nt(i,j,k)
endif
Expand All @@ -1039,8 +1015,8 @@ subroutine transport_on_nrho (Time, Dens, Adv_vel)
enddo
enddo

do j=jsc,jec
do i=isc,iec
do j = jsc,jec
do i = isc,iec
work1(i,j,k_rho) = (tmp(1,i,j)+work1(i,j,k_rho))*Grd%dyte(i,j)*transport_convert*Grd%tmask(i,j,1)
work2(i,j,k_rho) = (tmp(2,i,j)+work2(i,j,k_rho))*Grd%dxtn(i,j)*transport_convert*Grd%tmask(i,j,1)
enddo
Expand Down Expand Up @@ -1094,6 +1070,7 @@ end subroutine transport_on_nrho
! approach. The weighting approach was
! unnecessary, and added more cost to the scheme.
!
! 2018: Rearranged looping to avoid scanning "impossible" levels by Russ Fiedler.
! </DESCRIPTION>
!
subroutine transport_on_rho (Time, Dens, Adv_vel)
Expand All @@ -1106,7 +1083,9 @@ subroutine transport_on_rho (Time, Dens, Adv_vel)
integer :: i, j, k, k_rho, potrho_nk
real :: work1(isc:iec,jsc:jec,size(Dens%potrho_ref))
real :: work2(isc:iec,jsc:jec,size(Dens%potrho_ref))
real :: tmp(2,isc:iec,jsc:jec)
real :: tmp(2,isc:iec,size(Dens%potrho_ref))

real :: rho_max, rho_min

if (.not.module_is_initialized) then
call mpp_error(FATAL, &
Expand All @@ -1121,36 +1100,47 @@ subroutine transport_on_rho (Time, Dens, Adv_vel)
work1(:,:,:) = 0.0
work2(:,:,:) = 0.0

do k_rho=1,potrho_nk
do j = jsc,jec
tmp(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
if (k_rho == 1) then
if(Dens%potrho(i,j,k) < Dens%potrho_bounds(k_rho+1)) then
tmp(1,i,j) = tmp(1,i,j) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,j) = tmp(2,i,j) + Adv_vel%vhrho_nt(i,j,k)
endif
elseif(k_rho < potrho_nk) then
if( (Dens%potrho_bounds(k_rho) <= Dens%potrho(i,j,k)) .and. &
(Dens%potrho(i,j,k) < Dens%potrho_bounds(k_rho+1)) ) then
tmp(1,i,j) = tmp(1,i,j) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,j) = tmp(2,i,j) + Adv_vel%vhrho_nt(i,j,k)
endif
else ! if (k_rho == potrho_nk) then
if(Dens%potrho_bounds(k_rho) <= Dens%potrho(i,j,k)) then
tmp(1,i,j) = tmp(1,i,j) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,j) = tmp(2,i,j) + Adv_vel%vhrho_nt(i,j,k)
endif
do k = 1,nk
rho_max = maxval(Dens%potrho(isc:iec,j,k),mask=Grd%tmask(isc:iec,j,k)==1)
if (rho_max == -huge(rho_max)) exit ! All land in this row and lower
rho_min = minval(Dens%potrho(isc:iec,j,k),mask=Grd%tmask(isc:iec,j,k)==1)

k_rho = 1
if (rho_min < Dens%potrho_bounds(k_rho+1)) then
do i = isc,iec
if (Dens%potrho(i,j,k) < Dens%potrho_bounds(k_rho+1)) then
tmp(1,i,k_rho) = tmp(1,i,k_rho) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,k_rho) = tmp(2,i,k_rho) + Adv_vel%vhrho_nt(i,j,k)
endif
enddo
endif
do k_rho = 2,potrho_nk-1
if (Dens%potrho_bounds(k_rho) > rho_max .or. Dens%potrho_bounds(k_rho+1) <= rho_min) cycle
do i = isc,iec
if ((Dens%potrho_bounds(k_rho) <= Dens%potrho(i,j,k)) .and. &
(Dens%potrho(i,j,k) < Dens%potrho_bounds(k_rho+1)) ) then
tmp(1,i,k_rho) = tmp(1,i,k_rho) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,k_rho) = tmp(2,i,k_rho) + Adv_vel%vhrho_nt(i,j,k)
endif
enddo
enddo
k_rho = potrho_nk
if (rho_max >= Dens%potrho_bounds(k_rho)) then
do i = isc,iec
if (Dens%potrho_bounds(k_rho) <= Dens%potrho(i,j,k)) then
tmp(1,i,k_rho) = tmp(1,i,k_rho) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,k_rho) = tmp(2,i,k_rho) + Adv_vel%vhrho_nt(i,j,k)
endif
enddo
endif
enddo

do j=jsc,jec
do i=isc,iec
work1(i,j,k_rho) = (tmp(1,i,j)+work1(i,j,k_rho))*Grd%dyte(i,j)*transport_convert*Grd%tmask(i,j,1)
work2(i,j,k_rho) = (tmp(2,i,j)+work2(i,j,k_rho))*Grd%dxtn(i,j)*transport_convert*Grd%tmask(i,j,1)
do k_rho = 1,potrho_nk
do i = isc,iec
work1(i,j,k_rho) = (tmp(1,i,k_rho)+work1(i,j,k_rho))*Grd%dyte(i,j)*transport_convert*Grd%tmask(i,j,1)
work2(i,j,k_rho) = (tmp(2,i,k_rho)+work2(i,j,k_rho))*Grd%dxtn(i,j)*transport_convert*Grd%tmask(i,j,1)
enddo
enddo
enddo
Expand Down Expand Up @@ -1204,6 +1194,7 @@ end subroutine transport_on_rho
! approach. The weighting approach was
! unnecessary, and added more cost to the scheme.
!
! 2018: Rearranged looping to avoid scanning "impossible" levels by Russ Fiedler.
! </DESCRIPTION>
!
subroutine transport_on_theta (Time, Dens, Theta, Adv_vel)
Expand All @@ -1218,7 +1209,9 @@ subroutine transport_on_theta (Time, Dens, Theta, Adv_vel)
integer :: i, j, k, k_theta, theta_nk, tau
real :: work1(isc:iec,jsc:jec,size(Dens%theta_ref))
real :: work2(isc:iec,jsc:jec,size(Dens%theta_ref))
real :: tmp(2,isc:iec,jsc:jec)
real :: tmp(2,isc:iec,size(Dens%theta_ref))

real :: theta_max, theta_min

tau = Time%tau

Expand All @@ -1230,36 +1223,47 @@ subroutine transport_on_theta (Time, Dens, Theta, Adv_vel)
work1(:,:,:) = 0.0
work2(:,:,:) = 0.0

do k_theta=1,theta_nk
do j = jsc,jec
tmp(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
if (k_theta == 1) then
if(Theta%field(i,j,k,tau) < Dens%theta_bounds(k_theta+1)) then
tmp(1,i,j) = tmp(1,i,j) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,j) = tmp(2,i,j) + Adv_vel%vhrho_nt(i,j,k)
endif
elseif(k_theta < theta_nk) then
if( (Dens%theta_bounds(k_theta) <= Theta%field(i,j,k,tau)) .and. &
(Theta%field(i,j,k,tau) < Dens%theta_bounds(k_theta+1)) ) then
tmp(1,i,j) = tmp(1,i,j) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,j) = tmp(2,i,j) + Adv_vel%vhrho_nt(i,j,k)
endif
else ! if (k_theta == theta_nk) then
if(Dens%theta_bounds(k_theta) <= Theta%field(i,j,k,tau)) then
tmp(1,i,j) = tmp(1,i,j) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,j) = tmp(2,i,j) + Adv_vel%vhrho_nt(i,j,k)
endif
do k = 1,nk
theta_max = maxval(Theta%field(isc:iec,j,k,tau),mask=Grd%tmask(isc:iec,j,k)==1)
if (theta_max == -huge(theta_max)) exit ! All land in this row and lower
theta_min = minval(Theta%field(isc:iec,j,k,tau),mask=Grd%tmask(isc:iec,j,k)==1)

k_theta = 1
if (theta_min < Dens%theta_bounds(k_theta+1)) then
do i = isc,iec
if (Theta%field(i,j,k,tau) < Dens%theta_bounds(k_theta+1)) then
tmp(1,i,k_theta) = tmp(1,i,k_theta) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,k_theta) = tmp(2,i,k_theta) + Adv_vel%vhrho_nt(i,j,k)
endif
enddo
endif
do k_theta = 2,theta_nk-1
if (Dens%theta_bounds(k_theta) > theta_max .or. Dens%theta_bounds(k_theta+1) <= theta_min ) cycle
do i = isc,iec
if ((Dens%theta_bounds(k_theta) <= Theta%field(i,j,k,tau)) .and. &
(Theta%field(i,j,k,tau) < Dens%theta_bounds(k_theta+1)) ) then
tmp(1,i,k_theta) = tmp(1,i,k_theta) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,k_theta) = tmp(2,i,k_theta) + Adv_vel%vhrho_nt(i,j,k)
endif
enddo
enddo
k_theta = theta_nk
if (theta_max >= Dens%theta_bounds(k_theta)) then
do i = isc,iec
if (Dens%theta_bounds(k_theta) <= Theta%field(i,j,k,tau)) then
tmp(1,i,k_theta) = tmp(1,i,k_theta) + Adv_vel%uhrho_et(i,j,k)
tmp(2,i,k_theta) = tmp(2,i,k_theta) + Adv_vel%vhrho_nt(i,j,k)
endif
enddo
endif
enddo

do j=jsc,jec
do i=isc,iec
work1(i,j,k_theta) = (tmp(1,i,j)+work1(i,j,k_theta))*Grd%dyte(i,j)*transport_convert*Grd%tmask(i,j,1)
work2(i,j,k_theta) = (tmp(2,i,j)+work2(i,j,k_theta))*Grd%dxtn(i,j)*transport_convert*Grd%tmask(i,j,1)
do k_theta = 1,theta_nk
do i = isc,iec
work1(i,j,k_theta) = (tmp(1,i,k_theta)+work1(i,j,k_theta))*Grd%dyte(i,j)*transport_convert*Grd%tmask(i,j,1)
work2(i,j,k_theta) = (tmp(2,i,k_theta)+work2(i,j,k_theta))*Grd%dxtn(i,j)*transport_convert*Grd%tmask(i,j,1)
enddo
enddo
enddo
Expand Down
28 changes: 21 additions & 7 deletions src/mom5/ocean_diag/ocean_tracer_diag.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -550,7 +550,7 @@ subroutine ocean_tracer_diag_init(Grid, Domain, Time, Time_steps, Thickness, T_p
id_tracer_on_rho(n) = register_diag_field ('ocean_model', 'temp_on_rho', &
Dens%potrho_axes(1:3), Time%model_time, &
'temperature on potential density surface', 'C', &
missing_value=missing_value, range=(/0.0,1.e3/))
missing_value=missing_value, range=(/-5.0,100./))
id_tracer_zrho_on_rho(n) = register_diag_field ('ocean_model', 'temp_zrho_on_rho', &
Dens%potrho_axes(1:3), Time%model_time, &
'temperature*dz/drho on potential density surface', 'degC * m/(kg/m^3)', &
Expand Down Expand Up @@ -4053,6 +4053,8 @@ subroutine diagnose_tracer_on_rho(Time, Dens, Tracer, ntracer)
integer :: potrho_nk
integer :: i, j, k, k_rho, tau
real, dimension(isd:ied,jsd:jed,size(Dens%potrho_ref(:))) :: tracer_on_rho
real, dimension(jsc:jec) :: rho_minj,rho_maxj
real :: rho_min,rho_max

potrho_nk = size(Dens%potrho_ref(:))
tau = Time%tau
Expand All @@ -4070,9 +4072,21 @@ subroutine diagnose_tracer_on_rho(Time, Dens, Tracer, ntracer)
! since the initial value for tracer_on_rho is 0.

! interpolate tracer field onto rho-surface
do k_rho=1,potrho_nk
do k=1,nk-1
do j=jsc,jec

do k = 1,nk-1
do j = jsc,jec
rho_maxj(j) = maxval(Dens%potrho(isc:iec,j,k+1),mask=Grd%tmask(isc:iec,j,k+1)==1.)
rho_minj(j) = minval(Dens%potrho(isc:iec,j,k),mask=Grd%tmask(isc:iec,j,k)==1.)
enddo
rho_max = maxval(rho_maxj)
rho_min = minval(rho_minj)
if (rho_max == -huge(rho_max)) exit ! only rock below this level
do k_rho = 1,potrho_nk
if (rho_max < Dens%potrho_ref(k_rho)) cycle
if (rho_min > Dens%potrho_ref(k_rho)) cycle
do j = jsc,jec
if (rho_maxj(j) < Dens%potrho_ref(k_rho)) cycle
if (rho_minj(j) > Dens%potrho_ref(k_rho)) cycle
do i=isc,iec
if( Dens%potrho_ref(k_rho) > Dens%potrho(i,j,k) ) then
if( Dens%potrho_ref(k_rho) <= Dens%potrho(i,j,k+1)) then
Expand All @@ -4089,9 +4103,9 @@ subroutine diagnose_tracer_on_rho(Time, Dens, Tracer, ntracer)
enddo

! ensure masking is applied to interpolated field
do k_rho=1,potrho_nk
do j=jsc,jec
do i=isc,iec
do k_rho = 1,potrho_nk
do j = jsc,jec
do i = isc,iec
tracer_on_rho(i,j,k_rho) = tracer_on_rho(i,j,k_rho)*Grd%tmask(i,j,1)
enddo
enddo
Expand Down

0 comments on commit 9f5866b

Please sign in to comment.