les_smagorinsky_sym.F90 6.78 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
#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))
               end if
#endif
               AmX(i,j) =  (                                           &
                              dudx                                     &
#ifndef SLICE_MODEL
                            - dvdy                                     &
#endif
                           )**2
                         + shearX(i,j)**2
               AmX(i,j) = SmagX2_2d(i,j)*DXX*DYX*sqrt(AmX(i,j))
            end if
         end do
#ifndef SLICE_MODEL
      end do
#endif
!$OMP END DO NOWAIT

#ifdef SLICE_MODEL
!$OMP BARRIER
!$OMP SINGLE
      AmX(imin-1:imax,j-1) = AmX(imin-1:imax,j)
      AmX(imin-1:imax,j+1) = AmX(imin-1:imax,j)
!$OMP END SINGLE
#endif

   end if

   if (present(AmU)) then
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
      do j=jmin-1,jmax+1
#endif
         do i=imin-1,imax
!           Note (KK): shearU(au=3) not available
!                      (however we only need AmU(au=[1|2]) in tracer_diffusion)
            if(au(i,j).eq.1 .or. au(i,j).eq.2) then
!              interpolate dudxU (see deformation_rates)
               if (au(i,j) .eq. 1) then
                  dudx = _HALF_*(dudxC(i,j) + dudxC(i+1,j))
               else
                  dudx = _ZERO_
               end if
               AmU(i,j) =  (                                           &
                              dudx                                     &
#ifndef SLICE_MODEL
                            - dvdyU(i,j)                               &
#endif
                           )**2                                        &
                         + shearU(i,j)**2
               AmU(i,j) = SmagU2_2d(i,j)*DXU*DYU*sqrt(AmU(i,j))
            end if
         end do
#ifndef SLICE_MODEL
      end do
#endif
!$OMP END DO NOWAIT

#ifdef SLICE_MODEL
!$OMP BARRIER
!$OMP SINGLE
      AmU(imin-1:imax,j+1) = AmU(imin-1:imax,j)
!$OMP END SINGLE
#endif

   end if

#ifndef SLICE_MODEL
   if (present(AmV)) then
!$OMP DO SCHEDULE(RUNTIME)
      do j=jmin-1,jmax
         do i=imin-1,imax+1
!           Note (KK): we only need AmV(av=[1|2]) in tracer_diffusion
!                      (shearV(av=3) cannot be calculated anyway)
            if(av(i,j).eq.1 .or. av(i,j).eq.2) then
!              interpolate dvdyV and shearV (see deformation_rates)
               if (av(i,j) .eq. 1) then
                  dvdy = _HALF_*(dvdyC(i,j) + dvdyC(i,j+1))
               else
                  dvdy = _ZERO_
               end if
               AmV(i,j) =  (                                           &
                              dudxV(i,j)                               &
                            - dvdy                                     &
                           )**2                                        &
                         + (_HALF_*(shearX(i-1,j) + shearX(i,j)))**2
               AmV(i,j) = SmagV2_2d(i,j)*DXV*DYV*sqrt(AmV(i,j))
            end if
         end do
      end do
!$OMP END DO NOWAIT
   end if
#endif

!$OMP END PARALLEL

   call toc(TIM_SMAG2D)
#ifdef DEBUG
   write(debug,*) 'Leaving les_smagorinsky_sym()'
   write(debug,*)
#endif
   return
   end subroutine les_smagorinsky_sym
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2011 - Hans Burchard and Karsten Bolding               !
!-----------------------------------------------------------------------