Skip to content

Commit

Permalink
memory for Txx instead hTxx
Browse files Browse the repository at this point in the history
  • Loading branch information
NoraLoose committed Dec 2, 2023
1 parent 3ec78a8 commit 2069dcf
Showing 1 changed file with 38 additions and 41 deletions.
79 changes: 38 additions & 41 deletions src/parameterizations/lateral/MOM_Zanna_Bolton.F90
Original file line number Diff line number Diff line change
Expand Up @@ -76,9 +76,9 @@ module MOM_Zanna_Bolton
Txy !< Subgrid stress xy component in q [L2 T-2 ~> m2 s-2]

real, dimension(:,:,:), allocatable :: &
hTxx_with_memory, & !< Memory-enhanced thickness-weighted subgrid stress xx component [L4 T-2 ~> m4 s-2]
hTyy_with_memory, & !< Memory-enhanced thickness-weighted subgrid stress yy component [L4 T-2 ~> m4 s-2]
hTxy_with_memory, & !< Memory-enhanced thickness-weighted subgrid stress xy component [L4 T-2 ~> m4 s-2]
Txx_with_memory, & !< Memory-enhanced subgrid stress xx component [L3 T-2 ~> m3 s-2]
Tyy_with_memory, & !< Memory-enhanced subgrid stress yy component [L3 T-2 ~> m3 s-2]
Txy_with_memory, & !< Memory-enhanced subgrid stress xy component [L3 T-2 ~> m3 s-2]
ZBtau !< Memory timescale for ZB forcing [T ~> s]

real, dimension(:,:), allocatable :: &
Expand Down Expand Up @@ -315,9 +315,9 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020)
allocate(CS%kappa_q(SZIB_(G),SZJB_(G))); CS%kappa_q(:,:) = 0.

if (CS%ZB_memory_type > 0) then
allocate(CS%hTxx_with_memory(SZI_(G),SZJ_(G),SZK_(GV))); CS%hTxx_with_memory(:,:,:) = 0.
allocate(CS%hTyy_with_memory(SZI_(G),SZJ_(G),SZK_(GV))); CS%hTyy_with_memory(:,:,:) = 0.
allocate(CS%hTxy_with_memory(SZIB_(G),SZJB_(G),SZK_(GV))); CS%hTxy_with_memory(:,:,:) = 0.
allocate(CS%Txx_with_memory(SZI_(G),SZJ_(G),SZK_(GV))); CS%Txx_with_memory(:,:,:) = 0.
allocate(CS%Tyy_with_memory(SZI_(G),SZJ_(G),SZK_(GV))); CS%Tyy_with_memory(:,:,:) = 0.
allocate(CS%Txy_with_memory(SZIB_(G),SZJB_(G),SZK_(GV))); CS%Txy_with_memory(:,:,:) = 0.
allocate(CS%ZBtau(SZI_(G),SZJ_(G),SZK_(GV))); CS%ZBtau(:,:,:) = 0.
endif

Expand Down Expand Up @@ -407,9 +407,9 @@ subroutine ZB_2020_end(CS)
endif

if (CS%ZB_memory_type > 0) then
deallocate(CS%hTxx_with_memory)
deallocate(CS%hTyy_with_memory)
deallocate(CS%hTxy_with_memory)
deallocate(CS%Txx_with_memory)
deallocate(CS%Tyy_with_memory)
deallocate(CS%Txy_with_memory)
deallocate(CS%ZBtau)
endif
end subroutine ZB_2020_end
Expand Down Expand Up @@ -602,26 +602,23 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS)
real :: sh_xx_q !< Horizontal tension interpolated to q point [T-1 ~> s-1]
real :: D !< deformation rate [T-1 ~> s-1]
real :: tau !< eddy memory time scale tau [T ~> s]
real :: h_neglect !< Thickness so small it can be lost in
! roundoff and so neglected [H ~> m or kg m-2]
integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq
integer :: i, j, k
character(len=160) :: mesg ! The text of an error message

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB

h_neglect = GV%H_subroundoff

call pass_var(CS%Txx, G%Domain)
call pass_var(CS%Tyy, G%Domain)
call pass_var(CS%Txy, G%Domain, position=CORNER)

!if (CS%hTxx_with_memory(20,20,1) == 0.) then
!if (CS%Txx_with_memory(20,20,1) == 0.) then
! ! This is lazy initial condition
! do k = 1,nz
! CS%hTxx_with_memory(:,:,k) = CS%Txx(:,:,k) * G%mask2dT(:,:)
! CS%hTyy_with_memory(:,:,k) = CS%Tyy(:,:,k) * G%mask2dT(:,:)
! CS%Txx_with_memory(:,:,k) = CS%Txx(:,:,k) * G%mask2dT(:,:)
! CS%Tyy_with_memory(:,:,k) = CS%Tyy(:,:,k) * G%mask2dT(:,:)
! enddo
!endif

Expand All @@ -638,17 +635,17 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS)
+ (CS%sh_xy(I-1,J,k) + CS%sh_xy(I,J-1,k)) )
D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h)
tau = CS%ZB_memory_const / (D + 1.0/CS%ZB_max_memory)
tau = 7776000.0 ! 90 days
!tau = 7776000.0 ! 90 days
CS%ZBtau(i,j,k) = tau

! Eulerian memory with implicit time stepping
CS%hTxx_with_memory(i,j,k) = tau / (tau + CS%dt) * CS%hTxx_with_memory(i,J,k) &
+ CS%dt / (tau + CS%dt) * h(i,j,k) * CS%Txx(i,j,k)
CS%hTyy_with_memory(i,j,k) = tau / (tau + CS%dt) * CS%hTyy_with_memory(i,J,k) &
+ CS%dt / (tau + CS%dt) * h(i,j,k) * CS%Tyy(i,j,k)
CS%Txx_with_memory(i,j,k) = tau / (tau + CS%dt) * CS%Txx_with_memory(i,J,k) &
+ CS%dt / (tau + CS%dt) * CS%Txx(i,j,k)
CS%Tyy_with_memory(i,j,k) = tau / (tau + CS%dt) * CS%Tyy_with_memory(i,J,k) &
+ CS%dt / (tau + CS%dt) * CS%Tyy(i,j,k)

CS%Txx(i,j,k) = CS%hTxx_with_memory(i,j,k) / (h(i,j,k) + h_neglect)
CS%Tyy(i,j,k) = CS%hTyy_with_memory(i,j,k) / (h(i,j,k) + h_neglect)
CS%Txx(i,j,k) = CS%Txx_with_memory(i,j,k)
CS%Tyy(i,j,k) = CS%Tyy_with_memory(i,j,k)
enddo ; enddo


Expand All @@ -658,12 +655,12 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS)
+ (CS%sh_xx(i+1,j,k) + CS%sh_xx(i,j+1,k)))
D = sqrt(CS%sh_xy(I,J,k) * CS%sh_xy(I,J,k) + sh_xx_q * sh_xx_q)
tau = CS%ZB_memory_const / (D + 1.0/CS%ZB_max_memory)
tau = 7776000.0 ! 90 days
!tau = 7776000.0 ! 90 days

CS%hTxy_with_memory(I,J,k) = tau / (tau + CS%dt) * CS%hTxy_with_memory(i,J,k) &
+ CS%dt / (tau + CS%dt) * CS%hq(I,J,k) * CS%Txy(I,J,k)
CS%Txy_with_memory(I,J,k) = tau / (tau + CS%dt) * CS%Txy_with_memory(i,J,k) &
+ CS%dt / (tau + CS%dt) * CS%Txy(I,J,k)

CS%Txy(I,J,k) = CS%hTxy_with_memory(I,J,k) / (CS%hq(I,J,k) + h_neglect)
CS%Txy(I,J,k) = CS%Txy_with_memory(I,J,k)

enddo ; enddo

Expand All @@ -688,10 +685,10 @@ subroutine tensor_advection(u, v, G, GV, CS)
! Local variables

real, dimension(SZI_(G),SZJ_(G)) :: &
hTxx_with_memory, & !< thickness-weighted and advected tensor-component Txx [L4 T-2 ~> m4 s-2]
hTyy_with_memory !< thickness-weighted and advected tensor-component Tyy [L4 T-2 ~> m4 s-2]
Txx_with_memory, & !< advected tensor-component Txx [L3 T-2 ~> m3 s-2]
Tyy_with_memory !< advected tensor-component Tyy [L3 T-2 ~> m3 s-2]
real, dimension(SZIB_(G),SZJB_(G)) :: &
hTxy_with_memory !< thickness-weighted and advected tensor-component Txy [L4 T-2 ~> m4 s-2]
Txy_with_memory !< advected tensor-component Txy [L3 T-2 ~> m3 s-2]

real :: uT, vT ! Interpolation of advective velocity to T points
real :: uQ, vQ ! Interpolation of advective velocity to Q points
Expand All @@ -708,9 +705,9 @@ subroutine tensor_advection(u, v, G, GV, CS)



call pass_var(CS%hTxx_with_memory, G%Domain)
call pass_var(CS%hTyy_with_memory, G%Domain)
call pass_var(CS%hTxy_with_memory, G%Domain, position=CORNER)
call pass_var(CS%Txx_with_memory, G%Domain)
call pass_var(CS%Tyy_with_memory, G%Domain)
call pass_var(CS%Txy_with_memory, G%Domain, position=CORNER)
call pass_vector(u,v,G%Domain)
call pass_var(G%dxCv,G%Domain)
call pass_var(G%dyCu,G%Domain)
Expand Down Expand Up @@ -760,16 +757,16 @@ subroutine tensor_advection(u, v, G, GV, CS)
y = 1. + dy / G%dyBu(i0,j0)
endif

hTxx_with_memory(i,j) = bilin_interp(CS%hTxx_with_memory(i0,j0,k), CS%hTxx_with_memory(i0,j0+1,k), &
CS%hTxx_with_memory(i0+1,j0,k), CS%hTxx_with_memory(i0+1,j0+1,k), x, y)
hTyy_with_memory(i,j) = bilin_interp(CS%hTyy_with_memory(i0,j0,k), CS%hTyy_with_memory(i0,j0+1,k), &
CS%hTyy_with_memory(i0+1,j0,k), CS%hTyy_with_memory(i0+1,j0+1,k), x, y)
Txx_with_memory(i,j) = bilin_interp(CS%Txx_with_memory(i0,j0,k), CS%Txx_with_memory(i0,j0+1,k), &
CS%Txx_with_memory(i0+1,j0,k), CS%Txx_with_memory(i0+1,j0+1,k), x, y)
Tyy_with_memory(i,j) = bilin_interp(CS%Tyy_with_memory(i0,j0,k), CS%Tyy_with_memory(i0,j0+1,k), &
CS%Tyy_with_memory(i0+1,j0,k), CS%Tyy_with_memory(i0+1,j0+1,k), x, y)
enddo; enddo

do j=js-1,je+1 ; do i=is-1,ie+1
! Save computations to storage array and enforcing zero B.C.
CS%hTxx_with_memory(i,j,k) = hTxx_with_memory(i,j) * G%mask2dT(i,j)
CS%hTyy_with_memory(i,j,k) = hTyy_with_memory(i,j) * G%mask2dT(i,j)
CS%Txx_with_memory(i,j,k) = Txx_with_memory(i,j) * G%mask2dT(i,j)
CS%Tyy_with_memory(i,j,k) = Tyy_with_memory(i,j) * G%mask2dT(i,j)
enddo ; enddo

do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1
Expand Down Expand Up @@ -811,14 +808,14 @@ subroutine tensor_advection(u, v, G, GV, CS)
y = 1. + dy / G%dyT(I0+1,j0+1)
endif

hTxy_with_memory(I,J) = bilin_interp(CS%hTxy_with_memory(I0,J0,k), CS%hTxy_with_memory(I0,J0+1,k), &
CS%hTxy_with_memory(I0+1,J0,k), CS%hTxy_with_memory(I0+1,J0+1,k), x, y)
Txy_with_memory(I,J) = bilin_interp(CS%Txy_with_memory(I0,J0,k), CS%Txy_with_memory(I0,J0+1,k), &
CS%Txy_with_memory(I0+1,J0,k), CS%Txy_with_memory(I0+1,J0+1,k), x, y)

enddo; enddo

do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+1
! Save computations to storage array and enforcing zero B.C.
CS%hTxy_with_memory(I,J,k) = hTxy_with_memory(I,J) * G%mask2dBu(I,J)
CS%Txy_with_memory(I,J,k) = Txy_with_memory(I,J) * G%mask2dBu(I,J)
enddo ; enddo

enddo ! k loop
Expand Down

0 comments on commit 2069dcf

Please sign in to comment.