How to combine visc2 & smagorinsky visc4?

Report or discuss software problems and other woes

Moderators: arango, robertson

Post Reply
Message
Author
stone

How to combine visc2 & smagorinsky visc4?

#1 Unread post by stone »

Hi, All
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

eg. In "step2d_LF_AM3.h"

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
…
This subroutine is for barotropic mode, and use visc2_r or visc4_r, not related to smagorinsky. Is this the correct way to forward step2d?

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/s
I don’t understand why it blows up since smagorinsky visc4 are the same, and setting traditional visc2=0.
Another 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.

Post Reply