!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module hmix_aniso 1,14
!BOP
! !MODULE: hmix_aniso
!
! !DESCRIPTION:
! This module contains routines for computing horizontal friction
! terms from the divergence of a stress derived from the
! rate-of-strain tensor; the stress is anisotropic in general,
! although the viscous coefficients may be chosen to make the
! stress isotropic.
!
! !REVISION HISTORY:
! SVN:$Id: hmix_aniso.F90 21356 2010-03-01 22:12:38Z njn01 $
! !USES:
use kinds_mod
use blocks
use communicate
use distribution
use domain
use constants
use broadcast
use POP_ReductionsMod
use gather_scatter
use diagnostics
use time_management
use grid
use io
use exit_mod
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: init_aniso, &
hdiffu_aniso
!EOP
!BOC
!-----------------------------------------------------------------------
!
! anisotropic viscosity variables
!
!-----------------------------------------------------------------------
integer (int_kind), parameter :: &! choices for alignment type
hmix_alignment_type_flow = 1, &
hmix_alignment_type_east = 2, &
hmix_alignment_type_grid = 3
integer (int_kind) :: &
hmix_alignment_itype, &! integer choice for alignment type
vconst_5 ! int parameter for variable visc form
logical (log_kind) :: &
lvariable_hmix_aniso, &! spatially varying mixing coeffs
lsmag_aniso ! smagorinski nonlinear viscosity
real (r8) :: &
visc_para, &! viscosity parallel to alignment direction
visc_perp, &! viscosity perpendicular to alignment direction
vconst_1, &! coefficients for variable viscosity form
vconst_2, &! coefficients for variable viscosity form
vconst_3, &! coefficients for variable viscosity form
vconst_4, &! coefficients for variable viscosity form
vconst_6, &! coefficients for variable viscosity form
vconst_7, &! coefficients for variable viscosity form
smag_lat, &! latitude at which to vary perp Smag visc
smag_lat_fact, &! coeff of latitude-depend Smag visc
smag_lat_gauss ! Gaussian width of latitude-dep Smag visc
!*** the following variables are used only with smag viscosity
real (r8) :: &
c_para, &! dimensionless smag coefficient for
c_perp, &! parallel and perpendicular directions
u_para, &! velocities for grid reynolds number viscous
u_perp ! limit in parallel and perpendicular directions
real (r8), dimension (:,:,:), allocatable :: &
K1E,K1W,K2N,K2S, &! metric factors
H1E,H1W,H2N,H2S, &! grid spacings
AMAX_CFL, &! 1/2 maximum cfl-allowed viscosity
DSMIN, &! min(DXU, DYU)
F_PARA_SMAG, &! horizontal variations of the
F_PERP_SMAG ! Smagorinsky coefficients
!-----------------------------------------------------------------------
!
! F_PARA and F_PERP represent the actual anisotropic viscosities
! used when lvariable_hmix_aniso is true. For lsmag_aniso option,
! they contain the time-invariant "numerical" parts of the
! anisotropic viscosity coefficients.
!
!-----------------------------------------------------------------------
real (r8), dimension(:,:,:,:), allocatable :: &
F_PARA, &! spatial dependence of viscosity
! parallel to alignment direction
F_PERP ! spatial dependence of viscosity
! perpendicular to alignment direction
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_aniso
! !INTERFACE:
subroutine init_aniso 1,29
! !DESCRIPTION:
! Initializes various choices, allocates necessary space and computes
! some time-independent arrays for anisotropic viscosity.
!
! !REVISION HISTORY:
! same as module
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
iblock, &! block loop index
i,j,k, &! loop index
cfl_warn1, &! flag for warning about CFL limit
cfl_warn2, &! flag for warning about CFL limit
nml_error ! error flag for namelist
real (r8), dimension(:,:), allocatable :: &
WORKA,WORKB ! local work arrays
real (r8) :: &
f_para_min, f_para_max, &
f_perp_min, f_perp_max
character (char_len) :: &
hmix_alignment_choice, &! choice of direction that breaks isotropy
var_viscosity_infile, &! filename for variable viscosity factor
var_viscosity_infile_fmt, & ! format (bin or nc) of above file
var_viscosity_outfile, &! file for output of internally-computed
var_viscosity_outfile_fmt ! viscosity and format of that file
character (16), parameter :: &
param_fmt = '(a15,2x,1pe12.5)'
integer(int_kind) :: errorCode
!-----------------------------------------------------------------------
!
! input namelist for setting various anisotropic viscosity options
!
!-----------------------------------------------------------------------
namelist /hmix_aniso_nml/ &
hmix_alignment_choice, lvariable_hmix_aniso, lsmag_aniso, &
visc_para, visc_perp, c_para, c_perp, u_para, u_perp, &
vconst_1, vconst_2, vconst_3, vconst_4, vconst_5, vconst_6, &
vconst_7, smag_lat, smag_lat_fact, smag_lat_gauss, &
var_viscosity_infile, var_viscosity_infile_fmt, &
var_viscosity_outfile, var_viscosity_outfile_fmt
!-----------------------------------------------------------------------
!
! read input namelist for anisotropic viscosity
!
! DEFAULT VALUES
! hmix_alignment_choice : flow-aligned
! lvariable_hmix_aniso : false
! lsmag_aniso : false
! visc_para : zero
! visc_perp : zero
! c_para : zero
! c_perp : zero
! u_para : zero
! u_perp : zero
! vconst_1 : 1.e7_r8
! vconst_2 : 24.5_r8
! vconst_3 : 0.2_r8
! vconst_4 : 1.e-8_r8
! vconst_5 : 3
! vconst_6 : 1.e7_r8
! vconst_7 : 45.0_r8
! smag_lat : 20.0_r8
! smag_lat_fact : 0.98_r8
! smag_lat_gauss : 98.0_r8
!
!-----------------------------------------------------------------------
hmix_alignment_choice = 'flow'
lvariable_hmix_aniso = .false.
lsmag_aniso = .false.
visc_para = c0
visc_perp = c0
c_para = c0
c_perp = c0
u_para = c0
u_perp = c0
vconst_1 = 1.e7_r8
vconst_2 = 24.5_r8
vconst_3 = 0.2_r8
vconst_4 = 1.e-8_r8
vconst_5 = 3
vconst_6 = 1.e7_r8
vconst_7 = 45.0_r8
smag_lat = 20.0_r8
smag_lat_fact = 0.98_r8
smag_lat_gauss = 98.0_r8
var_viscosity_infile = 'unknown_var_viscosity_infile'
var_viscosity_infile_fmt = 'bin'
var_viscosity_outfile = 'unknown_var_viscosity_outfile'
var_viscosity_outfile_fmt = 'bin'
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
do while (nml_error > 0)
read(nml_in, nml=hmix_aniso_nml,iostat=nml_error)
end do
if (nml_error == 0) close(nml_in)
endif
call broadcast_scalar(nml_error, master_task)
if (nml_error /= 0) then
call exit_POP(sigAbort, &
'ERROR reading hmix_aniso_nml')
endif
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,ndelim_fmt)
write(stdout,'(a30)') 'Anisotropic viscosity options:'
write(stdout,blank_fmt)
write(stdout,hmix_aniso_nml)
write(stdout,blank_fmt)
select case (hmix_alignment_choice(1:4))
case ('flow')
hmix_alignment_itype = hmix_alignment_type_flow
write(stdout,'(a47)') &
' Anisotropy aligned along local flow direction'
case ('east')
hmix_alignment_itype = hmix_alignment_type_east
write(stdout,'(a46)') &
' Anisotropy aligned along east-west direction'
case ('grid')
hmix_alignment_itype = hmix_alignment_type_grid
write(stdout,'(a43)') &
' Anisotropy aligned along grid coordinates'
case default
hmix_alignment_itype = -1000
end select
if (.not. lvariable_hmix_aniso) then
write(stdout,'(a52)') &
' Spatially variable anisotropic viscosity disabled'
endif
if (.not. lsmag_aniso) then
write(stdout,'(a36)') 'Using input anisotropic viscosities:'
write(stdout,param_fmt) ' visc_para = ',visc_para
write(stdout,param_fmt) ' visc_perp = ',visc_perp
else
write(stdout,'(a38)') 'Using nonlinear Smagorinski viscosites'
write(stdout,'(a37)') ' input anistropic smag coefficients:'
write(stdout,param_fmt) ' c_para = ',c_para
write(stdout,param_fmt) ' c_perp = ',c_perp
write(stdout,param_fmt) ' u_para = ',u_para
write(stdout,param_fmt) ' u_perp = ',u_perp
if (smag_lat_fact /= c0) then
write(stdout,'(a31)') 'Using latitudinal variation in '
write(stdout,'(a35)') 'Smagorinski perpendicular viscosity'
write(stdout,param_fmt) 'smag_lat =',smag_lat
write(stdout,param_fmt) 'smag_lat_fact =',smag_lat_fact
write(stdout,param_fmt) 'smag_lat_gauss=',smag_lat_gauss
endif
endif
if (lvariable_hmix_aniso) then
if (trim(var_viscosity_infile) == 'ccsm-internal') then
write(stdout,'(a38)') &
'Using variable anisotropic parameters:'
write(stdout,param_fmt) ' vconst_1 = ',vconst_1
write(stdout,param_fmt) ' vconst_2 = ',vconst_2
write(stdout,param_fmt) ' vconst_3 = ',vconst_3
write(stdout,param_fmt) ' vconst_4 = ',vconst_4
write(stdout,'(a15,2x,i6)') ' vconst_5 = ',vconst_5
write(stdout,param_fmt) ' vconst_6 = ',vconst_6
write(stdout,param_fmt) ' vconst_7 = ',vconst_7
if (trim(var_viscosity_outfile) /= &
'unknown_var_viscosity_outfile') then
write(stdout,'(a44,a)') &
'Computed variable viscosity output to file: ', &
trim(var_viscosity_outfile)
endif
else
write(stdout,'(a47,a)') &
'Variable anisotropic viscosity read from file: ', &
trim(var_viscosity_infile)
endif
endif
endif
call broadcast_scalar(hmix_alignment_itype, master_task)
call broadcast_scalar(lvariable_hmix_aniso, master_task)
call broadcast_scalar(lsmag_aniso, master_task)
call broadcast_scalar(visc_para, master_task)
call broadcast_scalar(visc_perp, master_task)
call broadcast_scalar(c_para, master_task)
call broadcast_scalar(c_perp, master_task)
call broadcast_scalar(u_para, master_task)
call broadcast_scalar(u_perp, master_task)
call broadcast_scalar(vconst_1, master_task)
call broadcast_scalar(vconst_2, master_task)
call broadcast_scalar(vconst_3, master_task)
call broadcast_scalar(vconst_4, master_task)
call broadcast_scalar(vconst_5, master_task)
call broadcast_scalar(vconst_6, master_task)
call broadcast_scalar(vconst_7, master_task)
call broadcast_scalar(smag_lat, master_task)
call broadcast_scalar(smag_lat_fact, master_task)
call broadcast_scalar(smag_lat_gauss, master_task)
call broadcast_scalar(var_viscosity_infile, master_task)
call broadcast_scalar(var_viscosity_infile_fmt, master_task)
call broadcast_scalar(var_viscosity_outfile, master_task)
call broadcast_scalar(var_viscosity_outfile_fmt, master_task)
if (hmix_alignment_itype == -1000) then
call exit_POP(sigAbort, &
'Unknown type for anisotropic alignment direction')
endif
!-----------------------------------------------------------------------
!
! allocate and initialize various 2D arrays:
! AMAX_CFL is 1/2 the maximum allowed viscosity due to CFL limit
!
!-----------------------------------------------------------------------
allocate(K1E(nx_block,ny_block,nblocks_clinic), &
K1W(nx_block,ny_block,nblocks_clinic), &
K2N(nx_block,ny_block,nblocks_clinic), &
K2S(nx_block,ny_block,nblocks_clinic), &
H1E(nx_block,ny_block,nblocks_clinic), &
H1W(nx_block,ny_block,nblocks_clinic), &
H2N(nx_block,ny_block,nblocks_clinic), &
H2S(nx_block,ny_block,nblocks_clinic), &
AMAX_CFL(nx_block,ny_block,nblocks_clinic), &
WORKA(nx_block,ny_block), &
WORKB(nx_block,ny_block))
if (lsmag_aniso) then
allocate(DSMIN(nx_block,ny_block,nblocks_clinic))
endif
do iblock=1,nblocks_clinic
H2S(:,:,iblock) = HTE(:,:,iblock)
H1W(:,:,iblock) = HTN(:,:,iblock)
H2N(:,:,iblock) = eoshift(H2S(:,:,iblock),dim=2,shift=1)
H1E(:,:,iblock) = eoshift(H1W(:,:,iblock),dim=1,shift=1)
WORKA = H2S(:,:,iblock) + H2N(:,:,iblock)
WORKB = eoshift(WORKA,dim=1,shift=-1)
K1W(:,:,iblock) = c2*(WORKA - WORKB)/ &
(WORKA + WORKB)/H1W(:,:,iblock)
K1E(:,:,iblock) = eoshift(K1W(:,:,iblock),dim=1,shift=1)
WORKA = H1W(:,:,iblock) + H1E(:,:,iblock)
WORKB = eoshift(WORKA,dim=2,shift=-1)
K2S(:,:,iblock) = c2*(WORKA - WORKB)/ &
(WORKA + WORKB)/H2S(:,:,iblock)
K2N(:,:,iblock) = eoshift(K2S(:,:,iblock),dim=2,shift=1)
AMAX_CFL(:,:,iblock) = p125/(dtu*(DXUR(:,:,iblock)**2 + &
DYUR(:,:,iblock)**2))
if (lsmag_aniso) then
DSMIN(:,:,iblock) = min(DXU(:,:,iblock),DYU(:,:,iblock))
endif
end do ! block loop
!-----------------------------------------------------------------------
!
! allocate space and compute spatial dependence of
! anisotropic viscosity coefficients here.
!
!-----------------------------------------------------------------------
if (lvariable_hmix_aniso) then
allocate(F_PARA(nx_block,ny_block,km,nblocks_clinic), &
F_PERP(nx_block,ny_block,km,nblocks_clinic))
if (trim(var_viscosity_infile) == 'ccsm-internal') then
call compute_ccsm_var_viscosity
else
call read_var_viscosity(trim(var_viscosity_infile), &
trim(var_viscosity_infile_fmt))
endif
do k=1,km
f_para_min = POP_GlobalMinval(F_PARA(:,:,k,:), POP_distrbClinic, errorCode)
f_para_max = POP_GlobalMaxval(F_PARA(:,:,k,:), POP_distrbClinic, errorCode)
f_perp_min = POP_GlobalMinval(F_PERP(:,:,k,:), POP_distrbClinic, errorCode)
f_perp_max = POP_GlobalMaxval(F_PERP(:,:,k,:), POP_distrbClinic, errorCode)
if (my_task == master_task) then
write(stdout,"(' vertical level = ',i3)") k
write(stdout,'(a14,1x,1pe12.5,3x,a12,1x,1pe12.5)') &
' Min F_PARA =',f_para_min, &
'Max F_PARA =',f_para_max
write(stdout,'(a14,1x,1pe12.5,3x,a12,1x,1pe12.5)') &
' Min F_PERP =',f_perp_min, &
'Max F_PERP =',f_perp_max
endif
enddo
endif
!-----------------------------------------------------------------------
!
! check that viscosities are less than 1/2 CFL limit
! taper where necessary
!
!-----------------------------------------------------------------------
if (.not. lsmag_aniso) then
if (lvariable_hmix_aniso) then
! Enforce limits through tapering
!$OMP PARALLEL DO PRIVATE(iblock,k)
do iblock=1,nblocks_clinic
do k=1,km
where (F_PARA(:,:,k,iblock) > &
AMAX_CFL(:,:,iblock))
F_PARA(:,:,k,iblock) = AMAX_CFL(:,:,iblock)
endwhere
where (F_PERP(:,:,k,iblock) > &
AMAX_CFL(:,:,iblock))
F_PERP(:,:,k,iblock) = AMAX_CFL(:,:,iblock)
endwhere
enddo
enddo
!$OMP END PARALLEL DO
if (trim(var_viscosity_outfile) /= &
'unknown_var_viscosity_outfile') then
call write_var_viscosity(trim(var_viscosity_outfile), &
trim(var_viscosity_outfile_fmt))
endif
else
cfl_warn1 = 0
cfl_warn2 = 0
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
if (ANY(visc_para > AMAX_CFL(:,:,iblock))) cfl_warn1 = 1
if (ANY(visc_perp > AMAX_CFL(:,:,iblock))) cfl_warn2 = 1
end do
!$OMP END PARALLEL DO
if (POP_GlobalSum(cfl_warn1,POP_distrbClinic,errorCode) /= 0) then
if (my_task == master_task) then
write(stdout,'(a21)') 'WARNING (hmix_aniso):'
write(stdout,'(a51)') &
' parallel viscosity > 1/2 CFL limit at some points'
endif
endif
if (POP_GlobalSum(cfl_warn2,POP_distrbClinic,errorCode) /= 0) then
if (my_task == master_task) then
write(stdout,'(a21)') 'WARNING (hmix_aniso):'
write(stdout,'(a56)') &
' perpendicular viscosity > 1/2 CFL limit at some points'
endif
endif
endif
endif
!-----------------------------------------------------------------------
!
! set up latitudinal variation of Smagorinsky viscosity
!
!-----------------------------------------------------------------------
if ( lsmag_aniso ) then
allocate(F_PARA_SMAG(nx_block,ny_block,max_blocks_clinic), &
F_PERP_SMAG(nx_block,ny_block,max_blocks_clinic))
F_PARA_SMAG = c1
!*** latitude dependence of F_PERP_SMAG for viscous CCSM runs
if (smag_lat_fact /= c0) then
!$OMP PARALLEL DO PRIVATE(iblock,WORKA)
do iblock=1,nblocks_clinic
WORKA = abs(ULAT(:,:,iblock))*radian
where (WORKA >= smag_lat)
F_PERP_SMAG(:,:,iblock) = c1 - smag_lat_fact* &
exp(-(WORKA - smag_lat)**2/smag_lat_gauss)
elsewhere
F_PERP_SMAG(:,:,iblock) = c1 - smag_lat_fact
endwhere
end do
!$OMP END PARALLEL DO
endif
endif
!-----------------------------------------------------------------------
!
! deallocate work space and write viscosity if required
!
!-----------------------------------------------------------------------
deallocate (WORKA,WORKB)
if (my_task == master_task) then
write(stdout,blank_fmt)
endif
!-----------------------------------------------------------------------
!EOC
end subroutine init_aniso
!***********************************************************************
!BOP
! !IROUTINE: hdiffu_aniso
! !INTERFACE:
subroutine hdiffu_aniso(k,HDUK,HDVK,UMIXK,VMIXK,this_block) 1,1
! !DESCRIPTION:
! This routine computes the viscous terms in the momentum equation
! as the divergence of a stress tensor which is linearly related
! to the rate-of-strain tensor with viscous coefficents $\nu_{\|}$
! and $\nu_{\bot}$. These coefficients represent energy dissipation
! in directions parallel and perpendicular to a specified
! alignment direction which breaks the isotropy of the dissipation.
!
! A functional approach is used to derived the discrete operator,
! which ensures positive-definite energy dissipation, provided
! $\nu_{\|} > \nu_{\bot}$.
!
! In the functional approach, each U-cell is subdivided horizontally
! into 4 subcells (or quarter-cells), and the strain and stress
! tensors are evaluated in each subcell.
!
! The viscosities may optionally be evaluated with Smagorinsky-like
! non-linear dependence on the deformation rate, which is
! proportional to the norm of the strain tensor. With the
! smagorinski option, the viscosities are evalutated as
! \begin{eqnarray}
! \nu_{\|} &=& \max[c_{\|} |D|ds^2),u_{\|} ds] \\
! \nu_{\bot} &=& \max[c_{\bot}|D|ds^2),u_{\bot}ds]
! \end{eqnarray}
! where $ds = \min(dx,dy)$, and the $|D|$ is the deformation rate:
! $|D| = \surd{2}|E|$, $c_{\|}$ and $c_{\bot}$ are dimensionless
! coefficients of order 1, and $u_{\|}$ and $u_{\bot}$ are
! velocities associated with the grid Reynolds number which
! determine a minimum background viscosity in regions where the
! nonlinear viscosity is too small to control grid-point noise.
! Typically $u_{\|}$ and $u_{\bot}$ are order 1 cm/s.
! The non-linear viscosities are also automatically limited
! to be no greater than 1/2 the maximum valued allowed by the
! viscous CFL limit.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: k ! depth level index
real (r8), dimension(nx_block,ny_block), intent(in) :: &
UMIXK, &! U velocity at level k and mix time level
VMIXK ! V velocity at level k and mix time level
type (block), intent(in) :: &
this_block ! block info for this sub block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(out) :: &
HDUK, &! Hdiff(Ub) at level k
HDVK ! Hdiff(Vb) at level k
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,iq, &! dummy loop index
bid ! local block address
real (r8), dimension(nx_block,ny_block) :: &
GE,GW,GN,GS, &! depth ratios in each direction
NORM1, NORM2, &! normals in each direction for alignment
HDIFFCFL ! for cfl diagnostics
real (r8) :: &
work1,work2,work3,work4,work5,work6,work7,work8, &
ue,uw,un,us, &! u in (e,w,n,s) mult. by gamma
ve,vw,vn,vs, &! v in (e,w,n,s) mult. by gamma
FX,FY ! area*friction terms
real (r8), dimension(nx_block,ny_block,4) :: &
A,B,C,D, &! viscous coefficients
VTMP1, VTMP2, &! temp viscosity factors
S11,S22,S12, &! stress tensor
E11,E22,E12 ! rate-of-strain tensor
!-----------------------------------------------------------------------
!
! Indices for quarter-cells
! X = U point
! N,S,E,W label north,south,east,west faces of U cell
! 1,2,3,4 label SW,NW,NE,SE quarter cells surrounding a U-point
!
!
! X ------- X ------- X
! | | |
! | 1 | 4 |
! | | |
! ---------------- N ----------------
! | | | | |
! | 3 | 2 | 3 | 2 |
! | | | | |
! X ------ W ------ X ------ E ------ X
! | | | | |
! | 4 | 1 | 4 | 1 |
! | | | | |
! ---------------- S ----------------
! | | |
! | 2 | 3 |
! | | |
! X ------- X ------- X
!
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! find local block address
!
!-----------------------------------------------------------------------
bid = this_block%local_id
HDUK = c0
HDVK = c0
!-----------------------------------------------------------------------
!
! compute gamma (depth ratios) for partial bottom cells
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
do j=this_block%jb-1,this_block%je+1
do i=this_block%ib-1,this_block%ie+1
GW(i,j) = min(DZU(i ,j,k,bid), &
DZU(i-1,j,k,bid))/DZU(i,j,k,bid)
GE(i,j) = min(DZU(i ,j,k,bid), &
DZU(i+1,j,k,bid))/DZU(i,j,k,bid)
GS(i,j) = min(DZU(i,j ,k,bid), &
DZU(i,j-1,k,bid))/DZU(i,j,k,bid)
GN(i,j) = min(DZU(i,j ,k,bid), &
DZU(i,j+1,k,bid))/DZU(i,j,k,bid)
end do
end do
else
GN = c1
GS = c1
GE = c1
GW = c1
endif
!-----------------------------------------------------------------------
!
! compute rate-of-strain tensor in each quarter-cell
!
!-----------------------------------------------------------------------
do j=this_block%jb-1,this_block%je+1
do i=this_block%ib-1,this_block%ie+1
uw = GW(i,j)*UMIXK(i-1,j)
ue = GE(i,j)*UMIXK(i+1,j)
us = GS(i,j)*UMIXK(i,j-1)
un = GN(i,j)*UMIXK(i,j+1)
vw = GW(i,j)*VMIXK(i-1,j)
ve = GE(i,j)*VMIXK(i+1,j)
vs = GS(i,j)*VMIXK(i,j-1)
vn = GN(i,j)*VMIXK(i,j+1)
work1 = (UMIXK(i,j) - uw)/H1W(i,j,bid)
work2 = (ue - UMIXK(i,j))/H1E(i,j,bid)
work3 = p5*K2S(i,j,bid)*(VMIXK(i,j) + vs)
work4 = p5*K2N(i,j,bid)*(VMIXK(i,j) + vn)
E11(i,j,1) = work1 + work3
E11(i,j,2) = work1 + work4
E11(i,j,3) = work2 + work4
E11(i,j,4) = work2 + work3
work1 = (VMIXK(i,j) - vs)/H2S(i,j,bid)
work2 = (vn - VMIXK(i,j))/H2N(i,j,bid)
work3 = p5*K1W(i,j,bid)*(UMIXK(i,j) + uw)
work4 = p5*K1E(i,j,bid)*(UMIXK(i,j) + ue)
E22(i,j,1) = work1 + work3
E22(i,j,2) = work2 + work3
E22(i,j,3) = work2 + work4
E22(i,j,4) = work1 + work4
work1 = (UMIXK(i,j) - us)/H2S(i,j,bid)
work2 = (un - UMIXK(i,j))/H2N(i,j,bid)
work3 = (VMIXK(i,j) - vw)/H1W(i,j,bid)
work4 = (ve - VMIXK(i,j))/H1E(i,j,bid)
work5 = K2S(i,j,bid)*(UMIXK(i,j) + us)
work6 = K2N(i,j,bid)*(UMIXK(i,j) + un)
work7 = K1W(i,j,bid)*(VMIXK(i,j) + vw)
work8 = K1E(i,j,bid)*(VMIXK(i,j) + ve)
E12(i,j,1) = work1 + work3 - p5*(work5 + work7)
E12(i,j,2) = work2 + work3 - p5*(work6 + work7)
E12(i,j,3) = work2 + work4 - p5*(work6 + work8)
E12(i,j,4) = work1 + work4 - p5*(work5 + work8)
end do
end do
!-----------------------------------------------------------------------
!
! calculate normal directions for alignment
!
!-----------------------------------------------------------------------
select case (hmix_alignment_itype)
case (hmix_alignment_type_grid)
case (hmix_alignment_type_east)
NORM1 = cos(ANGLE(:,:,bid))
NORM2 = -sin(ANGLE(:,:,bid))
case (hmix_alignment_type_flow)
do j=this_block%jb-1,this_block%je+1
do i=this_block%ib-1,this_block%ie+1
work5 = sqrt(UMIXK(i,j)**2 + VMIXK(i,j)**2)
if (work5 >= eps) then
NORM1(i,j) = UMIXK(i,j)/work5 ! n1
NORM2(i,j) = VMIXK(i,j)/work5 ! n2
else
NORM1 = c0
NORM2 = c0
endif
end do
end do
end select
!-----------------------------------------------------------------------
!
! calculate some viscous factors based on chosen type
!
!-----------------------------------------------------------------------
if (lsmag_aniso) then ! smagorinsky nonlinear viscosity
do iq = 1,4 ! loop over quarter cells
do j=this_block%jb-1,this_block%je+1
do i=this_block%ib-1,this_block%ie+1
! deformation rate |D| = sqrt(2)|E| (use temp work6)
! note E12 is twice the (1,2) component of the strain tensor
work6 = sqrt(c2*(E11(i,j,iq)**2 + E22(i,j,iq)**2) + &
E12(i,j,iq)**2)
! compute nonlinear viscosities (ds = min(dx,dy))
! max[c_para|D|ds**2),u_para*ds]
! max[c_perp|D|ds**2),u_perp*ds]
VTMP1(i,j,iq) = c_para*F_PARA_SMAG(i,j,bid)*work6* &
DSMIN(i,j,bid)*DSMIN(i,j,bid)
VTMP2(i,j,iq) = c_perp*F_PERP_SMAG(i,j,bid)*work6* &
DSMIN(i,j,bid)*DSMIN(i,j,bid)
end do
end do
end do
if (lvariable_hmix_aniso) then
do iq = 1,4 ! loop over quarter cells
do j=this_block%jb-1,this_block%je+1
do i=this_block%ib-1,this_block%ie+1
VTMP1(i,j,iq) = max( VTMP1(i,j,iq), &
F_PARA(i,j,k,bid))
VTMP2(i,j,iq) = max( VTMP2(i,j,iq), &
F_PERP(i,j,k,bid))
end do
end do
end do
endif
! taper viscosities when greater than 1/2 CFL limit
do iq = 1,4
do j=this_block%jb-1,this_block%je+1
do i=this_block%ib-1,this_block%ie+1
VTMP1(i,j,iq) = min( VTMP1(i,j,iq ), &
AMAX_CFL(i,j,bid)) ! parallel viscosity
VTMP2(i,j,iq) = min( VTMP2(i,j,iq ), &
AMAX_CFL(i,j,bid)) ! perpendicular visc
enddo
enddo
enddo
else ! no smagorinski nonlinearity
if (lvariable_hmix_aniso) then
do iq=1,4
VTMP1(:,:,iq) = F_PARA(:,:,k,bid) ! parallel viscosity
VTMP2(:,:,iq) = F_PERP(:,:,k,bid) ! perpendicular visc
end do
else
VTMP1 = visc_para
VTMP2 = visc_perp
endif
endif
!-----------------------------------------------------------------------
!
! compute viscous coefficients A,B,C,D multiplying strain tensor
! use only one of the following anisotropic forms:
! grid-aligned, east-aligned or flow-aligned
!
! compute stress tensor from strain tensor
!
!-----------------------------------------------------------------------
select case (hmix_alignment_itype)
case (hmix_alignment_type_grid)
A = p5*(VTMP1+VTMP2)
B = p5*(VTMP1+VTMP2)
C = c0
D = VTMP2
case (hmix_alignment_type_east, &
hmix_alignment_type_flow)
do iq = 1,4
do j=this_block%jb-1,this_block%je+1
do i=this_block%ib-1,this_block%ie+1
A(i,j,iq) = p5*(VTMP1(i,j,iq)+VTMP2(i,j,iq)) - &
c2*(VTMP1(i,j,iq)-VTMP2(i,j,iq))* &
(NORM1(i,j)*NORM2(i,j))**2
B(i,j,iq) = p5*(VTMP1(i,j,iq)+VTMP2(i,j,iq)) - &
c2*(VTMP1(i,j,iq)-VTMP2(i,j,iq))* &
(NORM1(i,j)*NORM2(i,j))**2
C(i,j,iq) = (VTMP1(i,j,iq)-VTMP2(i,j,iq))* &
NORM1(i,j)*NORM2(i,j)*(NORM1(i,j)**2-NORM2(i,j)**2)
D(i,j,iq) = VTMP2(i,j,iq) + &
c2*(VTMP1(i,j,iq)-VTMP2(i,j,iq))* &
(NORM1(i,j)*NORM2(i,j))**2
end do
end do
end do
end select
do iq = 1,4
do j=this_block%jb-1,this_block%je+1
do i=this_block%ib-1,this_block%ie+1
S11(i,j,iq) = A(i,j,iq)*E11(i,j,iq) - B(i,j,iq)*E22(i,j,iq) + &
C(i,j,iq)*E12(i,j,iq)
S22(i,j,iq) = - B(i,j,iq)*E11(i,j,iq) + A(i,j,iq)*E22(i,j,iq) - &
C(i,j,iq)*E12(i,j,iq)
S12(i,j,iq) = C(i,j,iq)*(E11(i,j,iq) - E22(i,j,iq)) + &
D(i,j,iq)*E12(i,j,iq)
end do
end do
end do
!-----------------------------------------------------------------------
!
! compute friction terms from stresses
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! x-component
!
!-----------------------------------------------------------------------
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
work1 = H2S(i,j,bid)*S11(i,j,1) + H2N(i,j,bid)*S11(i,j,2)
work2 = H2S(i,j,bid)*S11(i,j,4) + H2N(i,j,bid)*S11(i,j,3)
work3 = (H2S(i+1,j,bid)*S11(i+1,j,1) + &
H2N(i+1,j,bid)*S11(i+1,j,2))*GE(i,j)
work4 = (H2S(i-1,j,bid)*S11(i-1,j,4) + &
H2N(i-1,j,bid)*S11(i-1,j,3))*GW(i,j)
!*** <h2*S11>_e - <h2*S11>_w
FX = p25*(work2 + work3 - work1 - work4)
work1 = H1W(i,j,bid)*S12(i,j,1) + H1E(i,j,bid)*S12(i,j,4)
work2 = H1W(i,j,bid)*S12(i,j,2) + H1E(i,j,bid)*S12(i,j,3)
work3 = (H1W(i,j+1,bid)*S12(i,j+1,1) + &
H1E(i,j+1,bid)*S12(i,j+1,4))*GN(i,j)
work4 = (H1W(i,j-1,bid)*S12(i,j-1,2) + &
H1E(i,j-1,bid)*S12(i,j-1,3))*GS(i,j)
!*** <h1*S12>_n*(1+h2n*k2n/2)
!*** -<h1*S12>_s*(1-h2s*k2s/2)
FX = FX + p25*((work2 + work3)* &
(c1 + p5*H2N(i,j,bid)*K2N(i,j,bid)) &
- (work1 + work4)* &
(c1 - p5*H2S(i,j,bid)*K2S(i,j,bid)))
work1 = H2S(i,j,bid)*S22(i,j,1) + H2N(i,j,bid)*S22(i,j,2)
work2 = H2S(i,j,bid)*S22(i,j,4) + H2N(i,j,bid)*S22(i,j,3)
work3 = (H2S(i+1,j,bid)*S22(i+1,j,1) + &
H2N(i+1,j,bid)*S22(i+1,j,2))*GE(i,j)
work4 = (H2S(i-1,j,bid)*S22(i-1,j,4) + &
H2N(i-1,j,bid)*S22(i-1,j,3))*GW(i,j)
!*** -<h2*S22>_e*h1e*k1e/2
!*** -<h2*S22>_w*h1w*k1w/2
FX = FX - p125*((work2 + work3)*H1E(i,j,bid)*K1E(i,j,bid) &
+ (work1 + work4)*H1W(i,j,bid)*K1W(i,j,bid))
!-----------------------------------------------------------------------
!
! y-component
!
!-----------------------------------------------------------------------
work1 = H1W(i,j,bid)*S22(i,j,1) + H1E(i,j,bid)*S22(i,j,4)
work2 = H1W(i,j,bid)*S22(i,j,2) + H1E(i,j,bid)*S22(i,j,3)
work3 = (H1W(i,j+1,bid)*S22(i,j+1,1) + &
H1E(i,j+1,bid)*S22(i,j+1,4))*GN(i,j)
work4 = (H1W(i,j-1,bid)*S22(i,j-1,2) + &
H1E(i,j-1,bid)*S22(i,j-1,3))*GS(i,j)
!*** <h1*S22>_n - <h1*S22>_s
FY = p25*(work2 + work3 - work1 - work4)
work1 = H2S(i,j,bid)*S12(i,j,1) + H2N(i,j,bid)*S12(i,j,2)
work2 = H2S(i,j,bid)*S12(i,j,4) + H2N(i,j,bid)*S12(i,j,3)
work3 = (H2S(i+1,j,bid)*S12(i+1,j,1) + &
H2N(i+1,j,bid)*S12(i+1,j,2))*GE(i,j)
work4 = (H2S(i-1,j,bid)*S12(i-1,j,4) + &
H2N(i-1,j,bid)*S12(i-1,j,3))*GW(i,j)
!*** <h2*S12>_e*(1+h1e*k1e/2)
!*** -<h2*S12>_w*(1-h1w*k1w/2)
FY = FY + p25*((work2 + work3)* &
(c1 + p5*H1E(i,j,bid)*K1E(i,j,bid)) &
- (work1 + work4)* &
(c1 - p5*H1W(i,j,bid)*K1W(i,j,bid)))
work1 = H1W(i,j,bid)*S11(i,j,1) + H1E(i,j,bid)*S11(i,j,4)
work2 = H1W(i,j,bid)*S11(i,j,2) + H1E(i,j,bid)*S11(i,j,3)
work3 = (H1W(i,j+1,bid)*S11(i,j+1,1) + &
H1E(i,j+1,bid)*S11(i,j+1,4))*GN(i,j)
work4 = (H1W(i,j-1,bid)*S11(i,j-1,2) + &
H1E(i,j-1,bid)*S11(i,j-1,3))*GS(i,j)
!*** -<h1*S11>_n*h2n*k2n/2 - <h1*S11>_s*h2s*k2s/2
FY = FY - p125*((work2 + work3)*H2N(i,j,bid)*K2N(i,j,bid) &
+ (work1 + work4)*H2S(i,j,bid)*K2S(i,j,bid))
!-----------------------------------------------------------------------
!
! divide by U-cell area
!
!-----------------------------------------------------------------------
if (KMU(i,j,bid) >= k) then
HDUK(i,j) = FX/UAREA(i,j,bid)
HDVK(i,j) = FY/UAREA(i,j,bid)
else
HDUK(i,j) = c0
HDVK(i,j) = c0
endif
enddo ! (i-loop)
enddo ! (j-loop)
!-----------------------------------------------------------------------
!
! compute horiz diffusion cfl diagnostics if required
!
!-----------------------------------------------------------------------
if (ldiag_cfl .and. .not. lsmag_aniso) then
if (lvariable_hmix_aniso) then
HDIFFCFL = merge(c4*F_PARA(:,:,k,bid)* &
(DXUR(:,:,bid)**2 + DYUR(:,:,bid)**2), &
c0, KMU(:,:,bid) > k)
else
HDIFFCFL = merge(c4*visc_para* &
(DXUR(:,:,bid)**2 + DYUR(:,:,bid)**2), &
c0, KMU(:,:,bid) > k)
endif
HDIFFCFL = abs(HDIFFCFL)
call cfl_hdiff(k,bid,HDIFFCFL,2,this_block)
endif
!-----------------------------------------------------------------------
!EOC
end subroutine hdiffu_aniso
!***********************************************************************
!BOP
! !IROUTINE: compute_ccsm_var_viscosity
! !INTERFACE:
subroutine compute_ccsm_var_viscosity 1,3
! !DESCRIPTION:
! This routine computes spatially-varying anisotropic viscosity
! coefficients similar to NCOM.
!
!gokhan THIS SECTION WILL BE UPDATED LATER !!!!!!!!
!
! F_PARA = min(F_PARA, AMAX_CFL),
! F_PERP = min(F_PERP, AMAX_CFL)
!
! are enforced in init_aniso and hdiffu_aniso for the lvariable_hmix_aniso
! and lsmag_aniso choices, respectively.
!
! beta_f (x,y) = 2 * omega * cos(ULAT(i,j)) / radius
! distance (x,y,z) = actual distance to vconst_5 points
! west of the nearest western boundary
! dx (x,y) = DXU(i,j)
! dy (x,y) = DYU(i,j)
! y (x,y) = ULAT(i,j), latitude of u/v grid pts in radians
!
! Also, vconst_x are input parameters defined in namelist hmix_aniso_nml.
! note that vconst_1, vconst_6, vconst_4, and vconst_7 have dimensions of cm^2/s,
! cm^2/s, 1/cm, and degrees (of latitude) respectively. vconst_5 is an INTEGER.
!
! NOTE: The nearest western boundary computations are done along the
! model longitudinal grid lines. Therefore, the viscosity
! coefficients based on these are only approximate in the high
! Northern Hemisphere when used with non-grid-aligned options.
!
!-------------------------------------------------------------------------
!
!
! !REVISION HISTORY:
! written by: Stephen Yeager 3/2000
! modified by: GD (08/2001)
! !INPUT PARAMETERS:
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j,k,ie,is,n,ii, &! dummy loop counters
index, &! index of nearest western boundary point
indexo, &! index of nwbp + buffer
ig,jg,igp1, &! dummy loop counters (global)
iblock, &! block counter
ncount ! number of western boundaries
real (r8) :: &
visc_vel_scale, &! viscous velocity scale
bu, bv ! B_viscosity terms
integer (int_kind), dimension(nx_global) :: &
iwp
integer (int_kind), dimension(nx_global,ny_global) :: &
NWBP_G, &! nearest western boundary point (global)
KMU_G ! kmt array at U points
real (r8), parameter :: &
vvsl = 1500.e2_r8, &! visc_vel_scale_length (cm)
dist_max = 1.e10_r8 ! distance for ACC region
real (r8), dimension(nx_global,ny_global) :: &
HTN_G, &! zonal distance between U points
DIST_G ! distance to nwbp (cm)
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
BETA_F, &! Coriolis parameter
DIST ! distance to nearest western boundary
!-----------------------------------------------------------------------
!
! initialize variables and gather global versions of HTN, KMU
!
!-----------------------------------------------------------------------
do iblock=1,nblocks_clinic
BETA_F(:,:,iblock) = c2*omega*cos(ULAT(:,:,iblock))/radius
F_PARA(:,:,:,iblock) = c0
F_PERP(:,:,:,iblock) = c0
end do
call gather_global(HTN_G, HTN, master_task, distrb_clinic)
call gather_global(KMU_G, KMU, master_task, distrb_clinic)
do k=1,km
!***
!*** determine nearest western boundary
!*** (master loops over global domain)
!***
if (my_task == master_task) then
NWBP_G = 0
do jg=1,ny_global
ncount = 0
iwp = 0
do ig=1,nx_global
igp1 = ig + 1
if ( ig == nx_global ) igp1 = 1
if ( (KMU_G(ig ,jg) < k) .and. &
(KMU_G(igp1,jg) >= k) ) then
ncount = ncount + 1
iwp(ncount) = ig
endif
enddo
if ( ncount > 0 ) then
do n=1,ncount-1
is = iwp(n)
ie = iwp(n+1) - 1
do ig = is, ie
NWBP_G(ig,jg) = is
enddo
enddo
do ig=1,nx_global
if (NWBP_G(ig,jg) == 0 ) &
NWBP_G(ig,jg) = iwp(ncount)
enddo
endif
enddo
!***
!*** determine distance to nearest western boundary
!*** (master loops over global domain)
!***
do jg=1,ny_global
do ig=1,nx_global
index = NWBP_G(ig,jg)
indexo = index + vconst_5
if ( index == 0) then
DIST_G(ig,jg) = dist_max
elseif ( ig >= index .and. ig <= indexo ) then
DIST_G(ig,jg) = c0
elseif ( (ig > indexo) ) then
DIST_G(ig,jg) = HTN_G(ig,jg) + DIST_G(ig-1,jg)
elseif ( ig < index ) then
if (indexo <= nx_global) then
if (ig == 1) then
DIST_G(ig,jg) = c0
do ii=indexo+1,nx_global
DIST_G(ig,jg) = HTN_G(ii,jg) + DIST_G(ig,jg)
enddo
DIST_G(ig,jg) = HTN_G(ig,jg) + DIST_G(ig,jg)
else
DIST_G(ig,jg) = HTN_G(ig,jg) + DIST_G(ig-1,jg)
endif
else
if (ig <= (indexo - nx_global)) then
DIST_G(ig,jg) = c0
else
DIST_G(ig,jg) = HTN_G(ig,jg) + DIST_G(ig-1,jg)
endif
endif
endif
enddo
enddo
endif ! (my_task == master_task)
call scatter_global(DIST, DIST_G, master_task, distrb_clinic, &
field_loc_NEcorner, field_type_scalar)
!***
!*** calculate viscosity coefficients
!*** (all processors loop over local subdomain)
!***
do iblock=1,nblocks_clinic
do j=1,ny_block
do i=1,nx_block
!gokhan
!**** visc_vel_scale = 100.0_r8 * exp(-zt(k)/vvsl)
!
! ! use bu as temp
!
!
!**** bu = 0.15_r8
!**** if ( abs(ULAT(i,j,iblock)*radian) < 30._r8 ) &
!**** bu = 0.425_r8*cos(pi*ULAT(i,j,iblock)*radian/30._r8) + &
!**** 0.575_r8
!
!**** F_PARA(i,j,k,iblock) = max(p5*visc_vel_scale*bu* &
!**** max(DXU(i,j,iblock), &
!**** DYU(i,j,iblock)), vconst_6 )
!
!**** bu = vconst_1*(c1 + vconst_2*(c1 + &
!**** cos((c2*ULAT(i,j,iblock))+pi)))
bv = (min(abs(ULAT(i,j,iblock)*radian),vconst_7) &
* 90._r8 / vconst_7) / radian
bu = vconst_1 * ( c1 + vconst_2*(c1 - cos(c2*bv)) )
bv = vconst_3*BETA_F(i,j,iblock)*(DXU(i,j,iblock)**3)
bv = bv*exp(-(vconst_4*DIST(i,j,iblock))**2)
F_PERP(i,j,k,iblock) = max(bu,bv)
F_PARA(i,j,k,iblock) = max(bv,vconst_6)
! the diffusive CFL criteria will be enforced in subroutine
! init_aniso
enddo ! i
enddo ! j
enddo ! iblock
enddo ! k
!-----------------------------------------------------------------------
!EOC
end subroutine compute_ccsm_var_viscosity
!***********************************************************************
!BOP
! !IROUTINE: read_var_viscosity
! !INTERFACE:
subroutine read_var_viscosity(infile, infile_fmt) 1,15
! !DESCRIPTION:
! This routine reads the spatially-varying anisotropic viscosity
! factors F_PARA, F_PERP from an input file.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
character (*), intent(in) :: &
infile, &! input file name (with path)
infile_fmt ! input file format (bin or nc)
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
type (datafile) :: &
visc_file ! io file type for viscosity file
type (io_field_desc) :: &
F_PARA_d, F_PERP_d ! descriptors for para/perp visc factors
type (io_dim) :: &
i_dim, j_dim, &! dimension descriptors for horiz dims
k_dim ! dimension descriptor for vertical levels
!-----------------------------------------------------------------------
!
! construct io file type and open for reading
!
!-----------------------------------------------------------------------
visc_file = construct_file(infile_fmt, full_name=infile, &
record_length=rec_type_dbl, &
recl_words=nx_global*ny_global)
call data_set(visc_file, 'open_read')
!-----------------------------------------------------------------------
!
! define variables to be read
!
!-----------------------------------------------------------------------
!*** define dimensions
i_dim = construct_io_dim('i', nx_global)
j_dim = construct_io_dim('j', ny_global)
k_dim = construct_io_dim('k', km)
F_PARA_d = construct_io_field('F_PARA', dim1=i_dim, dim2=j_dim, dim3=k_dim, &
long_name='parallel variable viscosity ', &
units ='cm2/s', &
grid_loc ='3221', &
field_loc = field_loc_NEcorner, &
field_type = field_type_scalar, &
d3d_array = F_PARA)
F_PERP_d = construct_io_field('F_PERP', dim1=i_dim, dim2=j_dim, dim3=k_dim, &
long_name='perpendicular variable viscosity', &
units ='cm2/s', &
grid_loc ='3221', &
field_loc = field_loc_NEcorner, &
field_type = field_type_scalar, &
d3d_array = F_PERP)
call data_set (visc_file, 'define', F_PARA_d)
call data_set (visc_file, 'define', F_PERP_d)
!-----------------------------------------------------------------------
!
! read arrays then clean up
!
!-----------------------------------------------------------------------
call data_set (visc_file, 'read', F_PARA_d)
call data_set (visc_file, 'read', F_PERP_d)
call destroy_io_field (F_PARA_d)
call destroy_io_field (F_PERP_d)
call data_set (visc_file, 'close')
call destroy_file(visc_file)
!-----------------------------------------------------------------------
!EOC
end subroutine read_var_viscosity
!***********************************************************************
!BOP
! !IROUTINE: write_var_viscosity
! !INTERFACE:
subroutine write_var_viscosity(outfile, outfile_fmt) 1,15
! !DESCRIPTION:
! This routine writes the spatially-varying anisotropic viscosity
! factors F_PARA, F_PERP which have been computed internally.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
character (*), intent(in) :: &
outfile, &! input file name (with path)
outfile_fmt ! input file format (bin or nc)
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
type (datafile) :: &
visc_file ! io file type for viscosity file
type (io_field_desc) :: &
F_PARA_d, F_PERP_d ! descriptors for para/perp visc factors
type (io_dim) :: &
i_dim, j_dim, &! dimension descriptors for horiz dims
k_dim ! dimension descriptor for vertical levels
!-----------------------------------------------------------------------
!
! construct io file type and open for writing
!
!-----------------------------------------------------------------------
visc_file = construct_file(outfile_fmt, root_name=outfile, &
record_length=rec_type_dbl, &
recl_words=nx_global*ny_global)
call data_set(visc_file, 'open')
!-----------------------------------------------------------------------
!
! define variables to be written
!
!-----------------------------------------------------------------------
!*** define dimensions
i_dim = construct_io_dim('i', nx_global)
j_dim = construct_io_dim('j', ny_global)
k_dim = construct_io_dim('k', km)
F_PARA_d = construct_io_field('F_PARA', dim1=i_dim, dim2=j_dim, dim3=k_dim, &
long_name='parallel variable viscosity', &
units ='cm2/s', &
grid_loc ='3221', &
field_loc = field_loc_NEcorner, &
field_type = field_type_scalar, &
d3d_array = F_PARA)
F_PERP_d = construct_io_field('F_PERP', dim1=i_dim, dim2=j_dim, dim3=k_dim, &
long_name='perpendicular variable viscosity', &
units ='cm2/s', &
grid_loc ='3221', &
field_loc = field_loc_NEcorner, &
field_type = field_type_scalar, &
d3d_array = F_PERP)
call data_set (visc_file, 'define', F_PARA_d)
call data_set (visc_file, 'define', F_PERP_d)
!-----------------------------------------------------------------------
!
! write arrays then clean up
!
!-----------------------------------------------------------------------
call data_set (visc_file, 'write', F_PARA_d)
call data_set (visc_file, 'write', F_PERP_d)
call destroy_io_field (F_PARA_d)
call destroy_io_field (F_PERP_d)
call data_set (visc_file, 'close')
call destroy_file(visc_file)
!-----------------------------------------------------------------------
!EOC
end subroutine write_var_viscosity
!***********************************************************************
end module hmix_aniso
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||