Merge branch 'adv' into nonhyd

src/2d/bdy_2d.F90

@@ -388,16 +388,18 @@
             select case (bdy_2d_type(l))
                case (ZERO_GRADIENT,CLAMPED_VEL,FLATHER_VEL)
                   do j = wfj(n),wlj(n)
-                     z(i,j) = z(i+1,j)
+                     if ( az(i+1,j) .ne. 0 ) z(i,j) = z(i+1,j)
                   end do
                case (SOMMERFELD)
                   do j = wfj(n),wlj(n)
+                     if ( az(i+1,j) .ne. 0 ) then
                      cfl = sqrt(g*_HALF_*(D(i,j)+D(i+1,j)))*dtm/DXU
                      z(i,j) = (                                             &
                                 (_ONE_ - _TWO_*cfl*(_ONE_-theta))*z (i  ,j) &
                                +(_ONE_ + _TWO_*cfl*(_ONE_-theta))*zo(i+1,j) &
                                -(_ONE_ - _TWO_*cfl*theta        )*z (i+1,j) &
                               )/(_ONE_ + _TWO_*cfl*theta        )
+                     end if
                   end do
                case (CLAMPED_ELEV,CLAMPED)
                   do j = wfj(n),wlj(n)
@@ -407,13 +409,16 @@
                   end do
                case (FLATHER_ELEV)
                   do j = wfj(n),wlj(n)
+                     if ( az(i+1,j) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i,j)+D(i+1,j))
 !                    Note (KK): note approximation of sse at vel-time stage
-                     a = ramp*bdy_data(kl) &
-                         - _TWO_/sqrt(g*depth)*(U(i,j)-ramp*bdy_data_u(kl)*depth)
-                     z(i,j) = max(a,-H(i,j)+min_depth)
+                     a = _TWO_/sqrt(g*depth)*(U(i,j)-ramp*bdy_data_u(kl)*depth)
+                     else
+                     a = _ZERO_
+                     end if
+                     z(i,j) = max(ramp*bdy_data(kl)-a,-H(i,j)+min_depth)
                      k = k+1
                      kl = kl + 1
                   end do
@@ -427,16 +432,18 @@
             select case (bdy_2d_type(l))
                case (ZERO_GRADIENT,CLAMPED_VEL,FLATHER_VEL)
                   do i = nfi(n),nli(n)
-                     z(i,j) = z(i,j-1)
+                     if ( az(i,j-1) .ne. 0 ) z(i,j) = z(i,j-1)
                   end do
                case (SOMMERFELD)
                   do i = nfi(n),nli(n)
+                     if ( az(i,j-1) .ne. 0 ) then
                      cfl = sqrt(g*_HALF_*(D(i,j-1)+D(i,j)))*dtm/DYVJM1
                      z(i,j) = (                                             &
                                 (_ONE_ - _TWO_*cfl*(_ONE_-theta))*z (i,j  ) &
                                +(_ONE_ + _TWO_*cfl*(_ONE_-theta))*zo(i,j-1) &
                                -(_ONE_ - _TWO_*cfl*theta        )*z (i,j-1) &
                               )/(_ONE_ + _TWO_*cfl*theta        )
+                     end if
                   end do
                case (CLAMPED_ELEV,CLAMPED)
                   do i = nfi(n),nli(n)
@@ -446,13 +453,16 @@
                   end do
                case (FLATHER_ELEV)
                   do i = nfi(n),nli(n)
+                     if ( az(i,j-1) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i,j-1)+D(i,j))
 !                    Note (KK): note approximation of sse at vel-time stage
-                     a = ramp*bdy_data(kl) &
-                         + _TWO_/sqrt(g*depth)*(V(i,j-1)-ramp*bdy_data_v(kl)*depth)
-                     z(i,j) = max(a,-H(i,j)+min_depth)
+                     a = _TWO_/sqrt(g*depth)*(V(i,j-1)-ramp*bdy_data_v(kl)*depth)
+                     else
+                     a = _ZERO_
+                     end if
+                     z(i,j) = max(ramp*bdy_data(kl)+a,-H(i,j)+min_depth)
                      k = k+1
                      kl = kl + 1
                   end do
@@ -466,16 +476,18 @@
             select case (bdy_2d_type(l))
                case (ZERO_GRADIENT,CLAMPED_VEL,FLATHER_VEL)
                   do j = efj(n),elj(n)
-                     z(i,j) = z(i-1,j)
+                     if ( az(i-1,j) .ne. 0 ) z(i,j) = z(i-1,j)
                   end do
                case (SOMMERFELD)
                   do j = efj(n),elj(n)
+                     if ( az(i-1,j) .ne. 0 ) then
                      cfl = sqrt(g*_HALF_*(D(i-1,j)+D(i,j)))*dtm/DXUIM1
                      z(i,j) = (                                             &
                                 (_ONE_ - _TWO_*cfl*(_ONE_-theta))*z (i  ,j) &
                                +(_ONE_ + _TWO_*cfl*(_ONE_-theta))*zo(i-1,j) &
                                -(_ONE_ - _TWO_*cfl*theta        )*z (i-1,j) &
                               )/(_ONE_ + _TWO_*cfl*theta        )
+                     end if
                   end do
                case (CLAMPED_ELEV,CLAMPED)
                   do j = efj(n),elj(n)
@@ -485,13 +497,16 @@
                   end do
                case (FLATHER_ELEV)
                   do j = efj(n),elj(n)
+                     if ( az(i-1,j) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i-1,j)+D(i,j))
 !                    Note (KK): note approximation of sse at vel-time stage
-                     a = ramp*bdy_data(kl) &
-                         + _TWO_/sqrt(g*depth)*(U(i-1,j)-ramp*bdy_data_u(kl)*depth)
-                     z(i,j) = max(a,-H(i,j)+min_depth)
+                     a = _TWO_/sqrt(g*depth)*(U(i-1,j)-ramp*bdy_data_u(kl)*depth)
+                     else
+                     a = _ZERO_
+                     end if
+                     z(i,j) = max(ramp*bdy_data(kl)+a,-H(i,j)+min_depth)
                      k = k+1
                      kl = kl + 1
                   end do
@@ -505,16 +520,18 @@
             select case (bdy_2d_type(l))
                case (ZERO_GRADIENT,CLAMPED_VEL,FLATHER_VEL)
                   do i = sfi(n),sli(n)
-                     z(i,j) = z(i,j+1)
+                     if ( az(i,j+1) .ne. 0 ) z(i,j) = z(i,j+1)
                   end do
                case (SOMMERFELD)
                   do i = sfi(n),sli(n)
+                     if ( az(i,j+1) .ne. 0 ) then
                      cfl = sqrt(g*_HALF_*(D(i,j)+D(i,j+1)))*dtm/DYV
                      z(i,j) = (                                             &
                                 (_ONE_ - _TWO_*cfl*(_ONE_-theta))*z (i,j  ) &
                                +(_ONE_ + _TWO_*cfl*(_ONE_-theta))*zo(i,j+1) &
                                -(_ONE_ - _TWO_*cfl*theta        )*z (i,j+1) &
                               )/(_ONE_ + _TWO_*cfl*theta        )
+                     end if
                   end do
                case (CLAMPED_ELEV,CLAMPED)
                   do i = sfi(n),sli(n)
@@ -524,13 +541,16 @@
                   end do
                case (FLATHER_ELEV)
                   do i = sfi(n),sli(n)
+                     if ( az(i,j+1) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i,j)+D(i,j+1))
 !                    Note (KK): note approximation of sse at vel-time stage
-                     a = ramp*bdy_data(kl) &
-                         - _TWO_/sqrt(g*depth)*(V(i,j)-ramp*bdy_data_v(kl)*depth)
-                     z(i,j) = max(a,-H(i,j)+min_depth)
+                     a = _TWO_/sqrt(g*depth)*(V(i,j)-ramp*bdy_data_v(kl)*depth)
+                     else
+                     a = _ZERO_
+                     end if
+                     z(i,j) = max(ramp*bdy_data(kl)-a,-H(i,j)+min_depth)
                      k = k+1
                      kl = kl + 1
                   end do
@@ -625,21 +645,25 @@
             select case (bdy_2d_type(l))
                case (FLATHER_VEL)
                   do j = wfj(n),wlj(n)
+                     if ( az(i+1,j) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i,j)+D(i+1,j))
 !                    Note (KK): note approximation of sse at vel-time stage
                      U(i,j) = ramp*bdy_data_u(kl)*depth &
                               - _HALF_*sqrt(g*depth)*(z(i,j)-ramp*bdy_data(kl))
+                     end if
                      k = k+1
                      kl = kl + 1
                   end do
                case (CLAMPED_VEL,CLAMPED)
                   do j = wfj(n),wlj(n)
+                     if ( az(i+1,j) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i,j)+D(i+1,j))
                      U(i,j) = ramp*bdy_data_u(kl)*depth
+                     end if
                      k = k+1
                      kl = kl + 1
                   end do
@@ -654,21 +678,25 @@
             select case (bdy_2d_type(l))
                case (FLATHER_VEL)
                   do j = efj(n),elj(n)
+                     if ( az(i-1,j) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i-1,j)+D(i,j))
 !                    Note (KK): note approximation of sse at vel-time stage
                      U(i-1,j) = ramp*bdy_data_u(kl)*depth &
                                 + _HALF_*sqrt(g*depth)*(z(i,j)-ramp*bdy_data(kl))
+                     end if
                      k = k+1
                      kl = kl + 1
                   end do
                case (CLAMPED_VEL,CLAMPED)
                   do j = efj(n),elj(n)
+                     if ( az(i-1,j) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i-1,j)+D(i,j))
                      U(i-1,j) = ramp*bdy_data_u(kl)*depth
+                     end if
                      k = k+1
                      kl = kl + 1
                   end do
@@ -686,21 +714,25 @@
             select case (bdy_2d_type(l))
                case (FLATHER_VEL)
                   do i = nfi(n),nli(n)
+                     if ( az(i,j-1) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i,j-1)+D(i,j))
 !                    Note (KK): note approximation of sse at vel-time stage
                      V(i,j-1) = ramp*bdy_data_v(kl)*depth &
                                 + _HALF_*sqrt(g*depth)*(z(i,j)-ramp*bdy_data(kl))
+                     end if
                      k = k+1
                      kl = kl + 1
                   end do
                case (CLAMPED_VEL,CLAMPED)
                   do i = nfi(n),nli(n)
+                     if ( az(i,j-1) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i,j-1)+D(i,j))
                      V(i,j-1) = ramp*bdy_data_v(kl)*depth
+                     end if
                      k = k+1
                      kl = kl + 1
                   end do
@@ -715,21 +747,25 @@
             select case (bdy_2d_type(l))
                case (FLATHER_VEL)
                   do i = sfi(n),sli(n)
+                     if ( az(i,j+1) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i,j)+D(i,j+1))
 !                    Note (KK): note approximation of sse at vel-time stage
                      V(i,j) = ramp*bdy_data_v(kl)*depth &
                               - _HALF_*sqrt(g*depth)*(z(i,j)-ramp*bdy_data(kl))
+                     end if
                      k = k+1
                      kl = kl + 1
                   end do
                case (CLAMPED_VEL,CLAMPED)
                   do i = sfi(n),sli(n)
+                     if ( az(i,j+1) .ne. 0 ) then
 !                    Note (KK): approximate interface depths at vel-time stage
 !                               by spatial mean at last sse-time stage
                      depth = _HALF_*(D(i,j)+D(i,j+1))
                      V(i,j) = ramp*bdy_data_v(kl)*depth
+                     end if
                      k = k+1
                      kl = kl + 1
                   end do

src/2d/m2d.F90

@@ -180,7 +180,6 @@
    integer                   :: i,j
    integer                   :: elev_method=1
    REALTYPE                  :: elev_const=_ZERO_
-   integer,parameter :: rk = kind(_ONE_)
    character(LEN = PATH_MAX) :: elev_file='elev.nc'
    namelist /m2d/ &
           elev_method,elev_const,elev_file,                       &
@@ -449,8 +448,8 @@
 
 !     This is only needed for proper flexible output
       where (az .eq. 0)
-         z  = -9999.0d0
-         zo = -9999.0d0
+         z  = -9999._rk
+         zo = -9999._rk
       end where
 
       call depth_update(zo,z,D,Dvel,DU,DV)

src/2d/variables_2d.F90

@@ -22,6 +22,7 @@
    IMPLICIT NONE
 !
 ! !PUBLIC DATA MEMBERS:
+   integer, parameter                  :: rk = kind(_ONE_)
    REALTYPE                            :: dtm
    REALTYPE,dimension(:,:),pointer     :: zo,z
    logical                             :: deformC=.false.
@@ -115,7 +116,6 @@
 !
 ! !LOCAL VARIABLES:
    integer                   :: rc
-   integer,parameter         :: rk = kind(_ONE_)
 !EOP
 !-------------------------------------------------------------------------
 !BOC
@@ -140,7 +140,7 @@
    break_stat = 0
 #endif
 
-   z  = -9999*_ONE_; zo =_ZERO_
+   z  = -9999._rk ; zo =_ZERO_
 
    U = _ZERO_; DU = _ZERO_; Uint = _ZERO_; UEx = _ZERO_
    V = _ZERO_; DV = _ZERO_; Vint = _ZERO_; VEx = _ZERO_
@@ -291,7 +291,6 @@
 !
 ! !LOCAL VARIABLES:
    logical :: used
-   integer,parameter :: rk = kind(_ONE_)
 !EOP
 !-----------------------------------------------------------------------
 !BOC

src/3d/m3d.F90

@@ -397,8 +397,8 @@
                num(i,j,:) = 1.e-15
                nuh(i,j,:) = 1.e-15
 #ifndef NO_BAROCLINIC
-               S(i,j,:)   = -9999.0
-               T(i,j,:)   = -9999.0
+               S(i,j,:)   = -9999._rk
+               T(i,j,:)   = -9999._rk
 #endif
             end if
          end do

src/3d/salinity.F90

@@ -19,7 +19,7 @@
    use domain, only: imin,jmin,imax,jmax,kmax,H,az,dry_z
 !KB   use get_field, only: get_3d_field
    use variables_2d, only: fwf_int
-   use variables_3d, only: S,hn,kmin
+   use variables_3d, only: rk,S,hn,kmin
    use meteo, only: metforcing,met_method,nudge_sss,sss
    use meteo, only: METEO_CONST,METEO_FROMFILE,METEO_FROMEXT
    use halo_zones, only: update_3d_halo,wait_halo,D_TAG,H_TAG
@@ -281,8 +281,8 @@
          stop 'init_salinity'
    end select
 
-   S(:,:,0) = -9999*_ONE_
-   forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) S(i,j,:) = -9999*_ONE_
+   S(:,:,0) = -9999._rk
+   forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) S(i,j,:) = -9999._rk
 
    call update_3d_halo(S,S,az,imin,jmin,imax,jmax,kmax,D_TAG)
    call wait_halo(D_TAG)

src/3d/temperature.F90

@@ -20,7 +20,7 @@
    use domain, only: imin,jmin,imax,kmax,jmax,H,az,dry_z
    use domain, only: ill,ihl,jll,jhl
    use domain, only: ilg,ihg,jlg,jhg
-   use variables_3d, only: T,rad,hn,kmin,A,g1,g2,heatflux_net
+   use variables_3d, only: rk,T,rad,hn,kmin,A,g1,g2,heatflux_net
    use meteo, only: metforcing,met_method,nudge_sst,sst
    use meteo, only: METEO_CONST,METEO_FROMFILE,METEO_FROMEXT
 !KB   use get_field, only: get_3d_field
@@ -359,8 +359,8 @@ end interface
          stop 'init_temperature'
    end select
 
-   T(:,:,0) = -9999*_ONE_
-   forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) T(i,j,:) = -9999*_ONE_
+   T(:,:,0) = -9999._rk
+   forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) T(i,j,:) = -9999._rk
 
    call update_3d_halo(T,T,az,imin,jmin,imax,jmax,kmax,D_TAG)
    call wait_halo(D_TAG)

src/3d/variables_3d.F90

@@ -115,6 +115,7 @@
    IMPLICIT NONE
 !
 ! !PUBLIC DATA MEMBERS:
+   integer, parameter                  :: rk = kind(_ONE_)
    REALTYPE                            :: dt,cnpar=0.9
    REALTYPE                            :: avmback=_ZERO_,avhback=_ZERO_
    logical                             :: do_numerical_analyses_3d=.false.
@@ -201,7 +202,6 @@
 !
 ! !LOCAL VARIABLES:
    integer                   :: rc
-   integer,parameter         :: rk = kind(_ONE_)
 !EOP
 !-------------------------------------------------------------------------
 !BOC
@@ -260,8 +260,8 @@
 #endif
 
 !  must be nonzero for gotm_fabm in case of calc_temp=F
-   g1 = -9999*_ONE_
-   g2 = -9999*_ONE_
+   g1 = -9999._rk
+   g2 = -9999._rk
 
 #ifdef STRUCTURE_FRICTION
    sf = _ZERO_
@@ -458,7 +458,7 @@
 !  Original author(s): Karsten Bolding & Jorn Bruggeman
 !
 ! !LOCAL VARIABLES:
-   integer,parameter :: rk = kind(_ONE_)
+!
 !EOP
 !-----------------------------------------------------------------------
 !BOC

src/3d/vertical_coordinates.F90

@@ -171,9 +171,9 @@ stop
    ! calculate the z-coordinate of the cell centers
    ! references to mean sea level
    zc(:,:,0)=-H(:,:)
-   zc(:,:,1)=-H(:,:) + 0.5*hn(:,:,1)
+   zc(:,:,1)=-H(:,:) + _HALF_*hn(:,:,1)
    do k=2,kmax
-      zc(:,:,k)=zc(:,:,k-1)+0.5*(hn(:,:,k-1)+hn(:,:,k))
+      zc(:,:,k)=zc(:,:,k-1)+_HALF_*(hn(:,:,k-1)+hn(:,:,k))
    end do
 
 #ifdef SLICE_MODEL

src/getm/initialise.F90

@@ -275,10 +275,11 @@
 #endif
    end if
 #endif
-   call init_output_processing()
 
    call init_les(runtype)
 
+   call init_output_processing()
+
    call init_register_all_variables(runtype)
 
    allocate(type_getm_host::output_manager_host)
@@ -375,6 +376,8 @@
       end if
    end if
 
+   call finalize_register_all_variables(runtype)
+
    if (.not. dryrun) then
       if (save_initial) then
       call output_manager_prepare_save(julianday, int(secondsofday), 0, int(MinN-1))

src/getm/register_all_variables.F90

@@ -11,6 +11,12 @@
 !
 ! !USES:
    use field_manager
+   use variables_2d, only: register_2d_variables
+   use variables_3d, only: register_3d_variables
+#ifdef _FABM_
+   use getm_fabm, only: register_fabm_variables
+#endif
+   use output_processing, only: register_processed_variables, finalize_register_processed_variables
    IMPLICIT NONE
 !
 !  default: all is private.
@@ -227,6 +233,7 @@
 #ifdef _FABM_
    call finalize_register_fabm_variables(fm)
 #endif
+   call finalize_register_processed_variables(fm)
 
    return
    end subroutine finalize_register_all_variables

src/les/les_smagorinsky_sym.F90

+#include "cppdefs.h"
+!-----------------------------------------------------------------------
+!BOP
+!
+! !ROUTINE: les_smagorinsky_sym - \label{les_smagorinsky_sym}
+!
+! !INTERFACE:
+   subroutine les_smagorinsky_sym(dudxC,dudxV,   &
+#ifndef SLICE_MODEL
+                                  dvdyC,dvdyU,   &
+#endif
+                                  shearX,shearU, &
+                                  AmC,AmX,AmU,AmV)
+
+!  Note (KK): keep in sync with interface in les.F90
+!
+! !DESCRIPTION:
+!
+!
+! !USES:
+   use variables_les, only: SmagC2_2d,SmagX2_2d,SmagU2_2d,SmagV2_2d
+   use domain, only: imin,imax,jmin,jmax,az,ax,au,av
+   use domain, only: dxc,dyc,dxx,dyx,dxu,dyu,dxv,dyv
+   use getm_timers, only: tic,toc,TIM_SMAG2D
+!$ use omp_lib
+   IMPLICIT NONE
+!
+! !INPUT PARAMETERS:
+   REALTYPE,dimension(E2DFIELD),intent(in) :: dudxC,dudxV
+#ifndef SLICE_MODEL
+   REALTYPE,dimension(E2DFIELD),intent(in) :: dvdyC,dvdyU
+#endif
+   REALTYPE,dimension(E2DFIELD),intent(in) :: shearX,shearU
+!
+! !OUTPUT PARAMETERS:
+   REALTYPE,dimension(E2DFIELD),intent(out),optional :: AmC,AmX,AmU,AmV
+!
+! !REVISION HISTORY:
+!  Original author(s): Knut Klingbeil
+!
+! !LOCAL VARIABLES:
+   REALTYPE :: dudx,dvdy
+   integer  :: i,j
+
+!EOP
+!-----------------------------------------------------------------------
+!BOC
+#ifdef DEBUG
+   integer, save :: Ncall = 0
+   Ncall = Ncall+1
+   write(debug,*) 'les_smagorinsky_sym() # ',Ncall
+#endif
+#ifdef SLICE_MODEL
+   j = jmax/2 ! this MUST NOT be changed!!!
+#endif
+   call tic(TIM_SMAG2D)
+
+!$OMP PARALLEL DEFAULT(SHARED)                                         &
+!$OMP          FIRSTPRIVATE(j)                                         &
+!$OMP          PRIVATE(i,dudx,dvdy)
+
+   if (present(AmC)) then
+!$OMP DO SCHEDULE(RUNTIME)
+#ifndef SLICE_MODEL
+      do j=jmin-1,jmax+1
+#endif
+         do i=imin-1,imax+1
+!           Note (KK): AmC(az=1) needed in uv_diffusion
+            if (az(i,j) .eq. 1) then
+!              interpolate shearC
+!              Note (KK): in W/E open boundary cells shearC(az=2) would
+!                         require shearU outside open boundary
+!                         in N/S open boundary cells shearC(az=2) would
+!                         require shearU(au=3)
+!                         however shearC(az=2) not needed
+               AmC(i,j) =  (                                           &
+                              dudxC(i,j)                               &
+#ifndef SLICE_MODEL
+                            - dvdyC(i,j)                               &
+#endif
+                           )**2                                        &
+                         + (_HALF_*(shearU(i-1,j) + shearU(i,j)))**2
+               AmC(i,j) = SmagC2_2d(i,j)*DXC*DYC*sqrt(AmC(i,j))
+            end if
+         end do
+#ifndef SLICE_MODEL
+      end do
+#endif
+!$OMP END DO NOWAIT
+
+#ifdef SLICE_MODEL
+!$OMP BARRIER
+!$OMP SINGLE
+      AmC(imin-1:imax+1,j+1) =  AmC(imin-1:imax+1,j)
+!$OMP END SINGLE
+#endif
+
+   end if
+
+   if (present(AmX)) then
+!$OMP DO SCHEDULE(RUNTIME)
+#ifndef SLICE_MODEL
+      do j=jmin-1,jmax
+#endif
+         do i=imin-1,imax
+            if (ax(i,j) .eq. 1) then
+!              interpolate dudxX and dvdyX
+               if (av(i,j).eq.3 .or. av(i+1,j).eq.3) then
+!                 Note (KK): western/eastern open bdy
+                  dudx = _ZERO_
+               else
+                  dudx = _HALF_*(dudxV(i,j) + dudxV(i+1,j  ))
+               end if
+#ifndef SLICE_MODEL
+               if (au(i,j).eq.3 .or. au(i,j+1).eq.3) then
+!                 Note (KK): northern/southern open bdy
+                  dvdy = _ZERO_
+               else
+                  dvdy = _HALF_*(dvdyU(i,j) + dvdyU(i  ,j+1))