!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| module hmix_aniso !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 ! !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) ! !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) !*** _e - _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) !*** _n*(1+h2n*k2n/2) !*** -_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) !*** -_e*h1e*k1e/2 !*** -_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) !*** _n - _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) !*** _e*(1+h1e*k1e/2) !*** -_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) !*** -_n*h2n*k2n/2 - _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 ! !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) ! !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) ! !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 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||