I want to add a new viscosity scheme which combines traditional harmonic (visc2) and smagorinsky biharmonic viscosity (following [Chassignet, Garraffo, 2001: Viscosity parameterization and the Gulf Stream separation]). I make this option a user defined CPP UV_VIS_C2S4 and actives it in .h file when build oceanM (However UV_VIS4/2 will not be defined when defined this new CPP, this is important!). Since this option is a combination, it should be invoked both like traditional UV_VIS2 and smagorinsky UV_VIS4. There are 2 simple configurations in ROMS to to refer to. Case 1 is "defined Smagorinsky and UV_VIS4”, in this case a Smagorinsky-like Biharmonic viscosity is applied only. The other case (Case 2) is “defined both UV_VIS2 and UV_VIS4” without Smagorinsky scheme.
When combining traditional harmonic (C2 part in UV_VIS_C2S4) and smagorinsky biharmonic (S4 part), S4 should be similar to Case1. Since Case1 defined UV_VIS4, everywhere "defined UV_VIS4" occurs, I change it to “defined UV_VIS4 || defined UV_VIS_C2S4”. At the same time, C2 is requested to behave as a traditional harmonic (visc2), so where there is a “defined UV_VIS2” (except simultaneously defined Smagorinsky, which is the case in “hmixing_mod" as shown below ), I change it to “defined UV_VIS2 || defined UV_VIS_C2S4” to behave as a traditional visc2.
UV_VIS_C2S4 should not be conflicted with itself in ROMS, just like Case 2 ( “defined UV_VIS2 and UV_VIS4” simultaneously without Smagorinsky scheme).
examples I modified are as follows:
eg. “hmixing_mod”,
(it shows we can not simultaneously defined smagorinsky 2/4 order viscosity, since only one “visc3d_r" is defined, here it is used by UV_VIS_C2S4 )
Code: Select all
# ifdef UV_SMAGORINSKY
!
! Smagorinsky viscosity.
!
# if defined UV_VIS2 || defined UV_VIS_C4S2
visc3d_r(i,j,k)=Hviscosity(i,j)+ &
& SmagorCoef*omn(i,j)*DefRate
# [b]elif defined UV_VIS4 || defined UV_VIS_C2S4[/b]
visc3d_r(i,j,k)=Hviscosity(i,j)+ &
& PecletCoef*(omn(i,j)**2)*DefRate
eg. In “rhs3d_mod" the visc2 and visc4 can be called sequentially, after modifying the corresponding subroutines such as uv3dmix_mod etc.
Code: Select all
# ifdef UV_VIS4
!
!-----------------------------------------------------------------------
! Compute horizontal, biharmonic mixing of momentum.
!-----------------------------------------------------------------------
!
CALL uv3dmix4 (ng, tile)
# endif
[b]ifdef UV_VIS_C2S4[/b]
!
!-----------------------------------------------------------------------
! Compute horizontal, harmonic&biharmonic mixing of momentum.
!-----------------------------------------------------------------------
!
CALL uv3dmix2 (ng, tile)
CALL uv3dmix4 (ng, tile)
# endif
Code: Select all
[b]# if defined UV_VIS2 || defined UV_VIS_C2S4[/b]
!
!-----------------------------------------------------------------------
! Add in horizontal harmonic viscosity.
!-----------------------------------------------------------------------
!
! Compute flux-components of the horizontal divergence of the stress
! tensor (m5/s2) in XI- and ETA-directions.
!
DO j=JstrV-1,Jend
DO i=IstrU-1,Iend
….
# if defined UV_VIS4 || defined UV_VIS_C2S4
!
!-----------------------------------------------------------------------
! Add in horizontal biharmonic viscosity. The biharmonic operator
! is computed by applying the harmonic operator twice.
!-----------------------------------------------------------------------
!
! Compute flux-components of the horizontal divergence of the stress
! tensor (m4 s^-3/2) in XI- and ETA-directions. It is assumed here
! that "visc4_r" and "visc4_p" are the squared root of the biharmonic
! viscosity coefficient. For momentum balance purposes, the total
! thickness "D" appears only when computing the second harmonic
! operator.
!
DO j=-1+JV_RANGE
DO i=-1+IU_RANGE
cff=visc4_r(i,j)*0.5_r8* &
& (pmon_r(i,j)* &
& ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,krhs)- &
& (pn(i-1,j)+pn(i ,j))*ubar(i ,j,krhs))- &
& pnom_r(i,j)* &
& ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,krhs)- &
& (pm(i,j-1)+pm(i,j ))*vbar(i,j ,krhs)))
UFx(i,j)=on_r(i,j)*on_r(i,j)*cff
VFe(i,j)=om_r(i,j)*om_r(i,j)*cff
END DO
END DO
DO j=JU_RANGE+1
DO i=IV_RANGE+1
cff=visc4_p(i,j)*0.5_r8* &
& (pmon_p(i,j)* &
& ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,krhs)- &
& (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,krhs))+ &
& pnom_p(i,j)* &
& ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,krhs)- &
& (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,krhs)))
# ifdef MASKING
cff=cff*pmask(i,j)
# endif
UFe(i,j)=om_p(i,j)*om_p(i,j)*cff
VFx(i,j)=on_p(i,j)*on_p(i,j)*cff
…When completed the all the preprocess flag modifications, I ran the model with UV_VIS_C2S4. It blows up at 10th step (10x600sec). However, when I run Case 1("defined Smagorinsky and UV_VIS4”), it runs smoothly. I try to run UV_VIS_C2S4 but setting visc2=0, to imitate Case 1, it blows up as 10th step. I check the smagorinsky coefs, both Case 1 and the blowed UV_VIS_C2S4 has the same values.
Code: Select all
Minimum horizontal diffusion coefficient = 2.46480815E+09 m4/s
Maximum horizontal diffusion coefficient = 1.75788579E+10 m4/s
Minimum horizontal viscosity coefficient = 2.46480815E+09 m4/s
Maximum horizontal viscosity coefficient = 1.75788579E+10 m4/sAnother thing as mentioned before, is using "visc2_r or visc4_r” (not related to smagorinsky) the correct way to forward step2d?
Looking forward to suggestions, especially who has experience to combine traditional harmonic (visc2) and smagorinsky biharmonic viscosity.