Chapter 7: Example of Advanced use

7.1 Flux correction

1. Flux correction (I) (ROMS)

How to implement a flux correction procedure ?

\(\frac{\partial T}{\partial t}+\frac{\partial F^x}{\partial x}+\frac{\partial F^y}{\partial y}=0\)

\(T^{n+1}_{i,j}=T^n_{i,j}-\frac{\Delta t}{\Delta x}\left( F^x_{i,j}-F^x_{i-1,j} \right)-\frac{\Delta t}{\Delta y}\left( F^y_{i,j}-F^y_{i,j-1} \right)\)

Store the fluxes on the fine and coarse grid and use the differences to correct the coarse grid points near the interface.

2. Flux correction (II) (ROMS)

  1. Declare flux arrays (\(F_x, F_y\)), compute then and make the summation on the fine grids over a parent time step

b. Declare the profiles corresponding to these fluxes

Call Agrif_Declare_Variable((/1,2/),(/1,1/),(/'x','y'/),(/0,1/),(/nx,ny/),F_xid)

Call Agrif_Declare_Variable((/2,1/),(/1,1/),(/'x','y'/),(/1,0/),(/nx,ny/),F_yid)

c. Set the update type

Call Agrif_Set_UpdateType(F_xid,update1=Agrif_Update_Copy , update2=Agrif_Update_Average )

Call Agrif_Set_UpdateType(F_yid,update1=Agrif_Update_Average , update2= Agrif_Update_Copy )

3. Make the flux correction inside a call to Agrif_Update

Call Agrif_Update_Variable (F_xid, procname = Correct_Flux_x)

Call Agrif_Update_Variable (F_yid, procname = Correct_Flux_y)

4. Flux correction (III) (ROMS)

Example: Correct_Flux_x

Subroutine Correct_Flux_x(tabres,i1,i2,j1,j2,before,nb,ndir)

logical :: western_side, eastern_side

if (before) then

        tabres = Fx(i1:i2,j1:j2)

else

        western_side  = (nb == 1).AND.(ndir == 1)

        eastern_side  = (nb == 1).AND.(ndir == 2)

        if (western_side) then

          do j=j1,j2
                 My_traceur(i1,j)   = My_traceur(i1,j)-(dt/(dx))(tabres(i1,j)-Fx(i1,j))
          enddo

        endif

        if (eastern_side) then

          do j=j1,j2
                My_traceur(i2+1,j) = My_traceur(i2+1,j)+(dt/(dx))(tabres(i2,j)-Fx(i2,j))
          enddo

        endif

endif

End Subroutine Correct_Flux_x

7.2 Wetting and drying

Wetting and drying (I) (MARS)

How to specify if we are in land or water through the SpecialValue ?

Example: Interp_Traceur

Agrif_Use_SpecialValue = .TRUE.

Agrif_SpecialValue = -999.

Call Agrif_Bc_Variable(tabtemp, my_Traceur_id,procname = Interp_my_Traceur)

Subroutine Inter_my_Traceur(tabres,i1,i2,j1,j2,before)

        if (before) then

                where ((h0(i1:i2,j1:j2)+ssh(i1:i2,j1:j2)) < min_depth)
                        tabres = Agrif_SpecialValue
                elsewhere
                        tabres = my_Traceur(i1:i2,j1:j2)
                endwhere

        else

                where (tabres /= Agrif_SpecialValue)
                        My_traceur = tabres
                endwhere

        endif

End Subroutine Inter_my_Traceur

7.3 Vertical refinement

7.3.1. Different vertical grids (I) (NEMO)

The grids have different vertical meshes. The vertical remapping between the two grids are done in the procnames routines.

Example: Update with different vertical grids

Subroutine Update_My_Traceur(tabres,i1,i2,k1,k2,before)

real,dimension(i1:i2,k1:k2) :: tabres

logical before

if (before) then
! on the child grid
        tabres = My_Traceur(i1:i2,k1:k2)
else
! on the parent grid (NB: tabres has a vertical dimension corresponding to the one of fine grid)

        do i=i1,i2
          call remap(tabres(i,:),My_traceur(i,:)) ! remap is a vertical remapping procedure
        enddo

endif

End Subroutine Update_My_Traceur

7.3.2. Different vertical grids (II) (NEMO)

The grids have different vertical mesh. The vertical remapping between the two grids are done in the procnames routines.

Example: Interpolation with different vertical grids

Subroutine Interp_My_Traceur(tabres,i1,i2,k1,k2,before)

real,dimension(i1:i2,k1:k2) :: tabres

logical before

if (before) then

! on the parent grid
        tabres = My_Traceur(i1:i2,k1:k2)

else

! on the child grid (NB: tabres has a vertical dimension corresponding to the one of parent grid)

        do i=i1,i2
          call remap(tabres(i,:),My_traceur(i,:)) ! remap is a vertical remapping procedure
        enddo

endif

End Subroutine Interp_My_Traceur

7.4 Interpolation with fine grid values

7.4.1. Tracer interpolation (I) (NEMO)

We want to use internal fine grid values during the interpolation

_images/interpolation_finegrid.png

7.4.2. Tracer interpolation (II) (NEMO)

Profiles

decal = Agrif_IRhox() + 1

Call Agrif_Declare_Profiles((/2,2/),(/1,1/),(/-decal,-decal/),(/nx+decal,ny+decal/),My_Traceur_id)

Call Agrif_Set_Bc(My_traceur_id,(/decal-1,decal/))

! NB: if not specified interp type is Agrif_constant

Call Agrif_Set_Bcinterp(My_traceur_id,interp12=Agrif_ppm,interp21=Agrif_ppm)

7.4.3. Tracer interpolation (III) (NEMO)

Interpolation procnames

Call Agrif_Bc_Variable(tabtemp, My_traceur_id,procname = Interp_My_Traceur)

Subroutine Interp_My_Traceur(tabres,i1,i2,j1,j2,before,nb,ndir)

if (before) then

        tabres = My_Traceur(i1:i2,j1:j2)

else

        western_side  = (nb == 1).AND.(ndir == 1)

        if (western_edge) then

                do j=j1,j2

                        if (u(1,j) <0.) then
                                My_Traceur(0,j)=a1*tabres(-Agrif_irhox(),j) + b1*My_Traceur(1,j) + (1.-a1-b1)*My_Traceur(2,j)
                        else
                                My_Traceur(0,j)=a2*tabres(-Agrif_irhox()-1,j)+b2*tabres(-Agrif_irhox(),j)+(1.-a2-b2)*My_Traceur(1,j)
                        endif

                enddo

        endif

endif

End Subroutine Interp_My_Traceur