!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module diag_bsf 1,16
!BOP
! !MODULE: diag_bsf
!
! !DESCRIPTION:
! This module contains code to compute the barotropic stream function
! !REVISION HISTORY:
! SVN:$Id$
! !USES:
use POP_KindsMod
use POP_ErrorMod
use POP_IOUnitsMod
use POP_GridHorzMod
use POP_FieldMod
use POP_HaloMod
use kinds_mod
use broadcast
use domain_size
use domain
use constants
use exit_mod
use global_reductions
use gather_scatter
use grid
use io
use registry
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: init_diag_bsf, &
pcg_diag_bsf_solver
!EOP
!BOC
integer (int_kind), parameter :: &
mxisle = 200 ! max number of islands
real (r8) :: &
r_norm ! normalization constant
integer (int_kind) :: &
nisle ! number of islands
integer (int_kind), dimension(mxisle) :: &
npts_is ! number of island boundary points
integer (int_kind), dimension(nx_block,ny_block,max_blocks_tropic) :: &
ISMASK_B ! island mask on barotropic decomposition
! coefficients for the 9-point operator
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
BSFC, BSFN, BSFS, BSFE, BSFW, BSFNE, BSFSE, BSFNW, BSFSW
real (r8),dimension(nx_block,ny_block,max_blocks_tropic) :: &
BSF_AN, &
BSF_AE, &
BSF_ANE, &
BSF_A0, &
BSF_RCALCT_B
real (r8), dimension(mxisle) :: &
aislandr ! island diagonal coefficient
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_diag_bsf
! !INTERFACE:
!
subroutine init_diag_bsf (BSF_in_contents_file, errorCode) 1,36
!
! !DESCRIPTION:
! This subroutine initializes variables associated with the computation of
! the barotropic stream function
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
logical (POP_logical), intent(in) :: &
BSF_in_contents_file ! true if BFS is in tavg_contents
! !INPUT PARAMETERS:
integer (POP_i4), intent(out) :: &
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
!
! initialize the necessary fields for the diagnostic barotropic
! streamfunction calculation
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
isle, &
iblock, &
i,j, &
bsf_error_flag, &! error flag for the missing TAVG_2D fields, etc.
nml_error ! namelist i/o error flag
integer (int_kind), dimension(:,:,:), allocatable :: &
ISMASK, &! island mask on baroclinic decomposition
WORK1
logical (log_kind), dimension(:,:,:), allocatable :: &
MASKI
logical (log_kind) :: &
maskt
real (r8), dimension(:,:,:), allocatable :: &
BSFC, BSFN, BSFS, BSFE, BSFW, BSFNE, BSFSE, BSFNW, BSFSW, &
WORK0,WORK2,RCALC_TMP ! coefficients for the 9-point operator
real (r8) :: &
xne,xse,xnw,xsw, &! contribution to coefficients from x,y
yne,yse,ynw,ysw, &! components of divergence
ase,anw,asw
logical (log_kind) :: ldiag_bsf
namelist /bsf_diagnostic_nml/ ldiag_bsf
!-----------------------------------------------------------------------
!
! read bsf diagnostic namelist
!
!-----------------------------------------------------------------------
errorCode = POP_Success
ldiag_bsf = .false.
if (my_task == master_task) then
open (nml_in, file=nml_filename, status='old', iostat=nml_error)
if (nml_error /= 0) then
nml_error = -1
else
nml_error = 1
endif
!*** keep reading until find right namelist
do while (nml_error > 0)
read(nml_in, nml=bsf_diagnostic_nml,iostat=nml_error)
end do
if (nml_error == 0) close(nml_in)
end if
call broadcast_scalar(nml_error, master_task)
if (nml_error /= 0) then
call exit_POP(sigAbort,'ERROR reading bsf_diagnostic_nml namelist')
endif
call broadcast_scalar(ldiag_bsf, master_task)
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,ndelim_fmt)
write(stdout,blank_fmt)
write(stdout,*) ' Barotropic Streamfunction Diagnostic:'
write(stdout,blank_fmt)
write(stdout,*) ' bsf_diagnostic_nml namelist settings:'
write(stdout,blank_fmt)
write(stdout,bsf_diagnostic_nml)
write(stdout,blank_fmt)
endif
!-----------------------------------------------------------------------
!
! consistency check (BSF must be requested in tavg_contents file)
!
!-----------------------------------------------------------------------
if (ldiag_bsf .and. .not. BSF_in_contents_file) &
call exit_POP(sigAbort,'ERROR: BSF must be requested in tavg_contents file' )
if (BSF_in_contents_file .and. .not. ldiag_bsf) &
call exit_POP(sigAbort,'ERROR: BSF is in tavg_contents file, but ldiag_bsf is false.')
if (.not. ldiag_bsf) then
return
else
call register_string ('ldiag_bsf')
endif
if (my_task == master_task) then
write (stdout,*) 'Initializing diagnostic BSF variables ....'
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
endif
allocate(WORK0 (nx_block,ny_block,nblocks_clinic), &
WORK2 (nx_block,ny_block,nblocks_clinic), &
BSFC (nx_block,ny_block,nblocks_clinic), &
BSFN (nx_block,ny_block,nblocks_clinic), &
BSFNE (nx_block,ny_block,nblocks_clinic), &
BSFE (nx_block,ny_block,nblocks_clinic), &
BSFSE (nx_block,ny_block,nblocks_clinic), &
BSFS (nx_block,ny_block,nblocks_clinic), &
BSFSW (nx_block,ny_block,nblocks_clinic), &
BSFW (nx_block,ny_block,nblocks_clinic), &
BSFNW (nx_block,ny_block,nblocks_clinic), &
RCALC_TMP(nx_block,ny_block,nblocks_clinic))
allocate(ISMASK (nx_block,ny_block,nblocks_clinic), &
WORK1 (nx_block,ny_block,nblocks_clinic))
allocate(MASKI (nx_block,ny_block,nblocks_clinic))
r_norm = c0
ISMASK_B = 0
WORK1 = 0
do iblock = 1, nblocks_clinic
WORK2(:,:,iblock) = p25*merge( 1, 0, KMU(:,:,iblock) > 0 )
enddo
!-----------------------------------------------------------------------
!
! 9-point coefficients
!
!-----------------------------------------------------------------------
do iblock = 1,nblocks_clinic
WORK0(:,:,iblock) = c0
BSFC (:,:,iblock) = c0
BSFN (:,:,iblock) = c0
BSFNE(:,:,iblock) = c0
BSFE (:,:,iblock) = c0
BSFSE(:,:,iblock) = c0
BSFS (:,:,iblock) = c0
BSFSW(:,:,iblock) = c0
BSFW (:,:,iblock) = c0
BSFNW(:,:,iblock) = c0
RCALC_TMP(:,:,iblock) = c0
do j=2,ny_block
do i=2,nx_block
xne = WORK2(i ,j ,iblock)*DXUR(i ,j ,iblock)*DYU(i ,j ,iblock)
xse = WORK2(i ,j-1,iblock)*DXUR(i ,j-1,iblock)*DYU(i ,j-1,iblock)
xnw = WORK2(i-1,j ,iblock)*DXUR(i-1,j ,iblock)*DYU(i-1,j ,iblock)
xsw = WORK2(i-1,j-1,iblock)*DXUR(i-1,j-1,iblock)*DYU(i-1,j-1,iblock)
yne = WORK2(i ,j ,iblock)*DYUR(i ,j ,iblock)*DXU(i ,j ,iblock)
yse = WORK2(i ,j-1,iblock)*DYUR(i ,j-1,iblock)*DXU(i ,j-1,iblock)
ynw = WORK2(i-1,j ,iblock)*DYUR(i-1,j ,iblock)*DXU(i-1,j ,iblock)
ysw = WORK2(i-1,j-1,iblock)*DYUR(i-1,j-1,iblock)*DXU(i-1,j-1,iblock)
BSFNE(i,j,iblock) = xne + yne
BSFSE(i,j,iblock) = xse + yse
BSFNW(i,j,iblock) = xnw + ynw
BSFSW(i,j,iblock) = xsw + ysw
BSFE(i,j,iblock) = xne + xse - yne - yse
BSFW(i,j,iblock) = xnw + xsw - ynw - ysw
BSFN(i,j,iblock) = yne + ynw - xne - xnw
BSFS(i,j,iblock) = yse + ysw - xse - xsw
!*** diagnoal coefficient
BSFC(i,j,iblock) = -(BSFNE(i,j,iblock) + BSFSE(i,j,iblock) &
+ BSFNW(i,j,iblock) + BSFSW(i,j,iblock))
WORK0 (i,j,iblock) = TAREA(i,j,iblock)**2
enddo ! i
enddo ! j
enddo ! iblock
call POP_HaloUpdate (BSFN, POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'init_diag_bsf: error updating halo for BSFN')
return
endif
call POP_HaloUpdate (BSFNE,POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'init_diag_bsf: error updating halo for BSFNE')
return
endif
call POP_HaloUpdate (BSFE, POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'init_diag_bsf: error updating halo for BSFE')
return
endif
call POP_HaloUpdate (BSFSE,POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'init_diag_bsf: error updating halo for BSFSE')
return
endif
call POP_HaloUpdate (BSFS, POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'init_diag_bsf: error updating halo for BSFS')
return
endif
call POP_HaloUpdate (BSFSW,POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'init_diag_bsf: error updating halo for BSFSW')
return
endif
call POP_HaloUpdate (BSFW, POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'init_diag_bsf: error updating halo for BSFW')
return
endif
call POP_HaloUpdate (BSFNW,POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'init_diag_bsf: error updating halo for BSFNW')
return
endif
call redistribute_blocks(BSF_AN, distrb_tropic, BSFN, distrb_clinic)
call redistribute_blocks(BSF_AE, distrb_tropic, BSFE, distrb_clinic)
call redistribute_blocks(BSF_ANE, distrb_tropic, BSFNE, distrb_clinic)
call redistribute_blocks(BSF_A0, distrb_tropic, BSFC , distrb_clinic)
r_norm = c1/global_sum(WORK0,distrb_clinic,field_loc_center,RCALCT)
call redistribute_blocks(BSF_RCALCT_B, distrb_tropic, RCALCT, distrb_clinic)
!-----------------------------------------------------------------------
!
! calculate island mask and redistribute ISMASK after it has been computed
!
!-----------------------------------------------------------------------
call island_mask_diag_bsf (ISMASK, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'init_diag_bsf: error in island_mask_diag_bsf')
return
endif
call redistribute_blocks(ISMASK_B, distrb_tropic, ISMASK, distrb_clinic)
!-----------------------------------------------------------------------
! calculate island coefficients
!-----------------------------------------------------------------------
do isle=1,mxisle
aislandr(isle) = c0
enddo
if ( my_task == master_task ) then
write (stdout,*) ' '
write (stdout,*) ' island coefficients: '
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
endif
do isle=1,nisle
WORK0 = c0
MASKI(:,:,:) = (ISMASK(:,:,:) == isle)
do iblock = 1,nblocks_clinic
do j=2,ny_block-1
do i=2,nx_block-1
!*** north faces
maskt = MASKI(i,j,iblock) .and. (ISMASK(i,j+1,iblock) == 0)
if (maskt) WORK0(i,j,iblock) = BSFN(i,j,iblock)
!*** northeast faces
maskt = MASKI(i,j,iblock) .and. (ISMASK(i+1,j+1,iblock) == 0)
if (maskt) WORK0(i,j,iblock) = WORK0(i,j,iblock) + BSFNE(i,j,iblock)
!*** northwest faces
maskt = MASKI(i,j,iblock) .and. (ISMASK(i-1,j+1,iblock) == 0)
if (maskt) WORK0(i,j,iblock) = WORK0(i,j,iblock) + BSFNW(i,j,iblock)
!*** south faces
maskt = MASKI(i,j,iblock) .and. (ISMASK(i,j-1,iblock) == 0)
if (maskt) WORK0(i,j,iblock) = WORK0(i,j,iblock) + BSFS(i,j,iblock)
!*** southeast faces
maskt = MASKI(i,j,iblock) .and. (ISMASK(i+1,j-1,iblock) == 0)
if (maskt) WORK0(i,j,iblock) = WORK0(i,j,iblock) + BSFSE(i,j,iblock)
!*** southwest faces
maskt = MASKI(i,j,iblock) .and. (ISMASK(i-1,j-1,iblock) == 0)
if (maskt) WORK0(i,j,iblock) = WORK0(i,j,iblock) + BSFSW(i,j,iblock)
!*** east faces
maskt = MASKI(i,j,iblock) .and. (ISMASK(i+1,j,iblock) == 0)
if (maskt) WORK0(i,j,iblock) = WORK0(i,j,iblock) + BSFE(i,j,iblock)
!*** west faces
maskt = MASKI(i,j,iblock) .and. (ISMASK(i-1,j,iblock) == 0)
if (maskt) WORK0(i,j,iblock) = WORK0(i,j,iblock) + BSFW(i,j,iblock)
end do ! i
end do ! j
end do ! iblock
!*** diagonal coefficient on island
aislandr(isle) = -c1/global_sum(WORK0, distrb_clinic, field_loc_center)
if ( my_task == master_task ) then
write (stdout,*) ' island ', isle, ' coefficient = ', c1/aislandr(isle)
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
endif
enddo ! isle
deallocate(WORK0,WORK2,BSFC,BSFN,BSFNE,BSFE,BSFSE,BSFS,BSFSW, &
BSFW,BSFNW,RCALC_TMP)
deallocate(ISMASK,WORK1)
deallocate(MASKI)
!-----------------------------------------------------------------------
!EOC
end subroutine init_diag_bsf
!***********************************************************************
!BOP
! !IROUTINE: pcg_diag_bsf
! !INTERFACE:
subroutine pcg_diag_bsf ( X, B, errorCode ) 1,23
! !DESCRIPTION:
!
! Conjugate-gradient solver for diagnosed barotropic streamfunction
! in free-surface formulation.
! Based upon pcg in module solvers.F90
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8),dimension(nx_block,ny_block,max_blocks_tropic),intent(in) :: &
B
! !INPUT/OUTPUT PARAMETERS:
real (r8),dimension(nx_block,ny_block,max_blocks_tropic),intent(inout) :: &
X
! !OUTPUT PARAMETERS:
integer(POP_i4), intent(out) :: &
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind), parameter :: mxscan =1000, ncheck = 100
real (r8) :: &
err = 1.0e-6_r8
integer (int_kind) :: &
m, isle, mscan, iblock
real (r8) :: &
eta0, eta1, eta2, rr, risle, rms_residual
logical (log_kind), dimension(nx_block,ny_block,max_blocks_tropic) :: &
CALCZ
real (r8), dimension(nx_block,ny_block,max_blocks_tropic) :: &
R, S, Q, WORK1, WORK2, WORK3, BSFCR
type (block) :: &
this_block ! block information for current block
!-----------------------------------------------------------------------
!
! initialization
!
!-----------------------------------------------------------------------
errorCode = POP_Success
call POP_HaloUpdate (X, POP_haloTropic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'pcg_diag_bsf: error updating halo for X')
return
endif
R = c0
Q = c0
WORK1 = c0
WORK2 = c0
CALCZ = (ISMASK_B == 0) ! mask for interior ocean points
!-----------------------------------------------------------------------
!
! diagonal preconditioner
!
!-----------------------------------------------------------------------
where (BSF_A0 /= c0)
BSFCR = c1/BSF_A0
elsewhere
BSFCR = c0
end where
where ( .not. CALCZ ) BSFCR = c0
do isle=1,nisle
where ( abs(ISMASK_B) == isle ) BSFCR = aislandr(isle)
enddo
!-----------------------------------------------------------------------
!
! compute initial residual
!
!-----------------------------------------------------------------------
call POP_HaloUpdate (X, POP_haloTropic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'pcg_diag_bsf: error updating halo for X')
return
endif
!$OMP PARALLEL DO PRIVATE(iblock,this_block)
do iblock=1,nblocks_tropic
this_block = get_block(blocks_tropic(iblock),iblock)
call btrop_operator_diag(WORK1,X,this_block,iblock)
WORK1(:,:,iblock) = B(:,:,iblock) - WORK1(:,:,iblock)
where (CALCZ(:,:,iblock))
R (:,:,iblock) = WORK1(:,:,iblock)
WORK2(:,:,iblock) = WORK1(:,:,iblock)
endwhere
end do ! block loop
!$OMP END PARALLEL DO
do isle=1,nisle ! island loop
WORK3 = merge(c1, c0, ISMASK_B == isle)
risle = global_sum( WORK1*WORK3,distrb_tropic,field_loc_center)
where ( abs(ISMASK_B) == isle ) R = risle
where ( ISMASK_B == isle ) WORK2 = risle/npts_is(isle)
enddo
rms_residual = sqrt(global_sum(WORK2*R,distrb_tropic, &
field_loc_center, BSF_RCALCT_B)*r_norm)
if ( my_task == master_task ) then
write (stdout,*) ' '
write (stdout,*) ' convergence info from pcg_diag_bsf: '
write (stdout,*) ' iter = ', 0, ' rms_resid = ', rms_residual
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
endif
!-----------------------------------------------------------------------
!
! initialize fields and scalars
!
!-----------------------------------------------------------------------
call POP_HaloUpdate (R, POP_haloTropic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'pcg_diag_bsf: error updating halo for R')
return
endif
S = c0
eta0 = c1
mscan = mxscan
!-----------------------------------------------------------------------
!
! iterate
!
!-----------------------------------------------------------------------
iter_loop: do m=1,mxscan
!$OMP PARALLEL DO PRIVATE(iblock,this_block,isle)
do iblock=1,nblocks_tropic
this_block = get_block(blocks_tropic(iblock),iblock)
!*** use diagnoal preconditioner
WORK1(:,:,iblock) = R(:,:,iblock) * BSFCR (:,:,iblock)
!-----------------------------------------------------------------------
!
! update conjugate direction vector s
!
!-----------------------------------------------------------------------
WORK2(:,:,iblock) = R(:,:,iblock)
do isle=1,nisle ! island loop
where ( ISMASK_B(:,:,iblock) == isle ) &
WORK2(:,:,iblock) = R(:,:,iblock)/npts_is(isle)
enddo
end do ! block loop
!$OMP END PARALLEL DO
!*** (r,(PC)r)
eta1 = global_sum( WORK2*WORK1, distrb_tropic, &
field_loc_center,BSF_RCALCT_B )
eta2 = eta1/eta0
call POP_HaloUpdate (WORK1, POP_haloTropic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'pcg_diag_bsf: error updating halo for WORK1')
return
endif
call POP_HaloUpdate (S,POP_haloTropic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'pcg_diag_bsf: error updating halo for S')
return
endif
!$OMP PARALLEL DO PRIVATE(iblock,this_block)
do iblock=1,nblocks_tropic
this_block = get_block(blocks_tropic(iblock),iblock)
S(:,:,iblock) = WORK1(:,:,iblock) + S(:,:,iblock)*eta2
!-----------------------------------------------------------------------
!
! compute As
!
!-----------------------------------------------------------------------
call btrop_operator_diag ( WORK1,S,this_block, iblock)
end do ! block loop
!$OMP END PARALLEL DO
where (CALCZ)
Q = WORK1
WORK2 = WORK1
endwhere
do isle=1,nisle ! island loop
WORK3 = merge(c1, c0, ISMASK_B == isle)
risle = global_sum( WORK1, distrb_tropic, field_loc_center,WORK3 )
where ( abs(ISMASK_B) == isle ) Q = risle
where ( ISMASK_B == isle ) WORK2 = risle/npts_is(isle)
enddo
!-----------------------------------------------------------------------
!
! compute next solution and residual
!
!-----------------------------------------------------------------------
call POP_HaloUpdate (Q,POP_haloTropic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'pcg_diag_bsf: error updating halo for Q')
return
endif
call POP_HaloUpdate (S,POP_haloTropic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'pcg_diag_bsf: error updating halo for S')
return
endif
eta0 = eta1
eta1 = eta0/global_sum( WORK2*S, distrb_tropic, &
field_loc_center, BSF_RCALCT_B )
!$OMP PARALLEL DO PRIVATE(iblock,this_block)
do iblock=1,nblocks_tropic
this_block = get_block(blocks_tropic(iblock),iblock)
X(:,:,iblock) = X(:,:,iblock) + eta1*S(:,:,iblock)
R(:,:,iblock) = R(:,:,iblock) - eta1*Q(:,:,iblock)
if (mod(m,ncheck) == 0) then
call btrop_operator_diag ( WORK1, X, this_block, iblock)
WORK1(:,:,iblock) = B(:,:,iblock) - WORK1(:,:,iblock)
where (CALCZ(:,:,iblock))
R(:,:,iblock) = WORK1(:,:,iblock)
WORK2(:,:,iblock) = WORK1(:,:,iblock)
endwhere
endif ! mod(m,ncheck)
end do ! block loop
!$OMP END PARALLEL DO
if (mod(m,ncheck) == 0) then
do isle=1,nisle ! island loop
WORK3 = merge(c1, c0, ISMASK_B == isle)
risle = global_sum( WORK1, distrb_tropic, &
field_loc_center, WORK3 )
where ( abs(ISMASK_B) == isle ) R = risle
where ( ISMASK_B == isle ) WORK2 = risle/npts_is(isle)
enddo
!-----------------------------------------------------------------------
!
! test for convergence
!
!-----------------------------------------------------------------------
rr = global_sum( WORK2*R, distrb_tropic, &
field_loc_center,BSF_RCALCT_B ) ! (r,r)
rms_residual = sqrt(rr*r_norm)
if ( my_task == master_task ) then
write (stdout,*) ' iter = ', m, ' rms_resid = ', rms_residual
endif
if (rms_residual < err) then
mscan = m
exit iter_loop
endif
endif ! mod(m,ncheck)
enddo iter_loop
if ( mscan == mxscan ) then
if ( my_task == master_task ) then
write (stdout,*) &
' WARNING (pcg_diag_bsf): No convergence after ', &
mxscan, ' scans'
endif
endif
!EOC
!-----------------------------------------------------------------------
end subroutine pcg_diag_bsf
!***********************************************************************
!BOP
! !IROUTINE: island_mask_diag_bsf
! !INTERFACE:
subroutine island_mask_diag_bsf (ISMASK, errorCode) 1,8
!
! !DESCRIPTION:
! This subroutine finds island mask for the barotropic streamfunction
! diagnostic computation
! !REVISION HISTORY:
! same as module
! !INPUT/OUTPUT PARAMETERS:
integer (int_kind), dimension(:,:,:), intent(inout) :: &
ISMASK ! island mask on baroclinic decomposition
! !OUTPUT PARAMETERS:
integer (POP_i4), intent(out) :: &
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
!
! output variables: nisle, npts_is, ISMASK
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: i, j, n, imax, error_code
integer (int_kind), dimension(:,:), allocatable :: &
KALL, KINT, I1, I2, ISMASK_G, WORK1, WORK2
logical (log_kind), dimension(:,:), allocatable :: &
LAND, LANDLEFT
!-----------------------------------------------------------------------
!
! allocate space
!
!-----------------------------------------------------------------------
errorCode = POP_Success
allocate(KALL (nx_global,ny_global), &
KINT (nx_global,ny_global), &
I1 (nx_global,ny_global), &
I2 (nx_global,ny_global), &
ISMASK_G (nx_global,ny_global), &
WORK1 (nx_global,ny_global), &
WORK2 (nx_global,ny_global) )
allocate(LAND (nx_global,ny_global), &
LANDLEFT (nx_global,ny_global) )
KALL = 0; KINT = 0; I1 = 0; I2 = 0; ISMASK_G = 0; WORK1 = 0; WORK2 = 0
call gather_global (WORK1, KMT, master_task,distrb_clinic)
error_code = 0
if ( my_task == master_task ) then
LAND = .false.
where ( WORK1 == 0 ) LAND = .true.
do j=1,ny_global
do i=1,nx_global
I1(i,j) = i
I2(i,j) = j
enddo
enddo
where ( .not. LAND )
KINT = nx_global * ny_global + 1
elsewhere
KINT = I1 + (I2-1) * nx_global
endwhere
KALL = KINT
imax = nx_global + 2 * ny_global
do i=1,imax
WORK1 = cshift(KALL, +1, 2)
WORK1(:,ny_global) = KALL(:,ny_global)
WORK2 = cshift(KALL, -1, 2)
WORK2(:,1) = KALL(:,1)
KALL = min(WORK1, WORK2, KALL)
WORK1 = cshift(KALL, +1, 2)
WORK1(:,ny_global) = KALL(:,ny_global)
WORK2 = cshift(KALL, -1, 2)
WORK2(:,1) = KALL(:,1)
KALL = min(WORK1, WORK2, KALL)
WORK1 = cshift(KALL, +1, 1)
WORK2 = cshift(KALL, -1, 1)
KALL = min(WORK1, WORK2, KALL)
WORK1 = cshift(KALL, +1, 1)
WORK2 = cshift(KALL, -1, 1)
KALL = min(WORK1, WORK2, KALL)
where (LAND) KINT = KALL
KALL = KINT
enddo
ISMASK_G = 0
LANDLEFT = LAND
where ( KINT == KINT(1,ny_global) ) ! north continent
ISMASK_G = 999 ! NOTE: nisle must be less than 999
LANDLEFT = .false.
endwhere
n = 0
do i=1,nx_global
do j=1,ny_global
if (LANDLEFT(i,j)) then
n = n + 1
where ( KINT == KINT(i,j) )
ISMASK_G = n
LANDLEFT = .false.
endwhere
endif
enddo
enddo
nisle = n
write (stdout,*) ' '
write (stdout,*) ' number of islands = ', nisle
if ( nisle > mxisle ) error_code = 1
endif ! if master_task
call broadcast_scalar (error_code, master_task)
if ( error_code /= 0) then
call exit_POP(sigAbort, &
'ERROR (island_mask_diag_bsf): nisle > mxisle')
endif
if ( my_task == master_task ) then
WORK1 = cshift(ISMASK_G, +1, 1)
WORK2 = cshift(ISMASK_G, -1, 1)
KALL = max(ISMASK_G, WORK1, WORK2)
WORK1 = cshift(ISMASK_G, +1, 2)
WORK2 = cshift(ISMASK_G, -1, 2)
KALL = max(KALL, WORK1, WORK2)
WORK1 = cshift(ISMASK_G, +1, 2)
WORK2 = cshift(WORK1, +1, 1)
KALL = max(KALL, WORK2)
WORK1 = cshift(ISMASK_G, +1, 2)
WORK2 = cshift(WORK1, -1, 1)
KALL = max(KALL, WORK2)
WORK1 = cshift(ISMASK_G, -1, 2)
WORK2 = cshift(WORK1, +1, 1)
KALL = max(KALL, WORK2)
WORK1 = cshift(ISMASK_G, -1, 2)
WORK2 = cshift(WORK1, -1, 1)
KALL = max(KALL, WORK2)
where (.not.LAND)
ISMASK_G = -KALL
endwhere
ISMASK_G = -ISMASK_G
write (stdout,*) ' '
do i=1,nisle
npts_is(i) = count(ISMASK_G == i)
write (stdout,*) ' island ', i,' # points on boundary = ', npts_is(i)
enddo
endif ! if master_task
call scatter_global(ISMASK, ISMASK_G, master_task, distrb_clinic, &
field_loc_NEcorner, field_type_scalar)
call POP_HaloUpdate (ISMASK, POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'island_mask_diag_bsf: error updating halo for ISMASK')
return
endif
call broadcast_scalar ( nisle, master_task )
call broadcast_array ( npts_is, master_task )
deallocate(KALL,KINT,I1,I2,ISMASK_G,WORK1,WORK2)
deallocate(LAND,LANDLEFT)
!-----------------------------------------------------------------------
!EOC
end subroutine island_mask_diag_bsf
!***********************************************************************
!BOP
! !IROUTINE: btrop_operator_diag
! !INTERFACE:
subroutine btrop_operator_diag(AX,X,this_block,bid) 3
! !DESCRIPTION:
! This routine applies the nine-point stencil operator for the
! barotropic solver, using weights local to this diagnostics
! module. It takes advantage of some 9pt weights being
! shifted versions of others.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,max_blocks_tropic), &
intent(in) :: &
X ! array to be operated on
type (block), intent(in) :: &
this_block ! block info for this block
integer (int_kind), intent(in) :: &
bid ! local block address for this block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,max_blocks_tropic), &
intent(out) :: &
AX ! nine point operator result (Ax)
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j ! dummy counters
!-----------------------------------------------------------------------
AX(:,:,bid) = c0
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
AX(i,j,bid) = BSF_A0 (i ,j ,bid)*X(i ,j ,bid) + &
BSF_AN (i ,j ,bid)*X(i ,j+1,bid) + &
BSF_AN (i ,j-1,bid)*X(i ,j-1,bid) + &
BSF_AE (i ,j ,bid)*X(i+1,j ,bid) + &
BSF_AE (i-1,j ,bid)*X(i-1,j ,bid) + &
BSF_ANE(i ,j ,bid)*X(i+1,j+1,bid) + &
BSF_ANE(i ,j-1,bid)*X(i+1,j-1,bid) + &
BSF_ANE(i-1,j ,bid)*X(i-1,j+1,bid) + &
BSF_ANE(i-1,j-1,bid)*X(i-1,j-1,bid)
end do
end do
!-----------------------------------------------------------------------
!EOC
end subroutine btrop_operator_diag
!BOP
! !IROUTINE: pcg_diag_bsf_solver
! !INTERFACE:
subroutine pcg_diag_bsf_solver(SOLN, RHS, errorCode) 1,5
! !DESCRIPTION:
! Solves the elliptic equation for bsf computation.
! First redistributes arrays to the barotropic block distribution for
! better performance of the solver.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,max_blocks_clinic), &
intent(in) :: &
RHS ! right-hand-side of linear system
! for blocks in baroclinic distribution
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,max_blocks_clinic), &
intent(inout) :: &
SOLN ! on input, initial guess
! on output, final solution for sfc pressure
! !OUTPUT PARAMETERS:
integer (POP_i4), intent(out) :: &
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
real (r8), dimension(nx_block,ny_block,max_blocks_tropic) :: &
S_TROPIC, &! surface pressure in barotropic distribution
RHS_TROPIC ! right hand side in barotropic distribution
!-----------------------------------------------------------------------
!
! switch to the barotropic distribution for iterative solvers
!
!-----------------------------------------------------------------------
errorCode = POP_Success
call redistribute_blocks(S_TROPIC, distrb_tropic, &
SOLN, distrb_clinic)
call redistribute_blocks(RHS_TROPIC, distrb_tropic, &
RHS, distrb_clinic)
!-----------------------------------------------------------------------
!
! call bsf diagnostic solver
!
!-----------------------------------------------------------------------
call pcg_diag_bsf(S_TROPIC,RHS_TROPIC, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'pcg_diag_bsf_solver: error in pcg_diag_bsf')
return
endif
!-----------------------------------------------------------------------
!
! switch solution back to the baroclinic distribution
!
!-----------------------------------------------------------------------
call redistribute_blocks(SOLN, distrb_clinic, &
S_TROPIC, distrb_tropic)
!-----------------------------------------------------------------------
!EOC
end subroutine pcg_diag_bsf_solver
!***********************************************************************
end module diag_bsf
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||