salinity.F90 14.6 KB
Newer Older
gotm's avatar
gotm committed
1 2 3 4 5 6 7 8 9 10
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE:  Salinity
!
! !INTERFACE:
   module salinity
!
! !DESCRIPTION:
11 12 13 14 15
!
! In this module, the salinity equation is processed by
! reading in the namelist {\tt salt} and initialising the salinity field
! (this is done in {\tt init\_salinity}),
! and calculating the advection-diffusion-equation (see {\tt do\_salinity}).
gotm's avatar
gotm committed
16 17
!
! !USES:
18
   use exceptions
19
   use domain, only: imin,jmin,imax,jmax,kmax,ioff,joff
20
   use domain, only: H,az
21
!KB   use get_field, only: get_3d_field
22
   use variables_2d, only: fwf_int
Knut's avatar
Knut committed
23
   use variables_3d, only: rk,S,hn,kmin
24
   use halo_zones, only: update_3d_halo,wait_halo,D_TAG,H_TAG
gotm's avatar
gotm committed
25 26 27 28 29
   IMPLICIT NONE
!
   private
!
! !PUBLIC DATA MEMBERS:
30
   public init_salinity, do_salinity, init_salinity_field
gotm's avatar
gotm committed
31 32
!
! !PRIVATE DATA MEMBERS:
kbk's avatar
kbk committed
33 34
   integer                   :: salt_method=1,salt_format=2
   character(len=PATH_MAX)   :: salt_file="t_and_s.nc"
35
   integer                   :: salt_field_no=1
kbk's avatar
kbk committed
36
   character(len=32)         :: salt_name='salt'
37
   REALTYPE                  :: salt_const=35*_ONE_
38
   integer                   :: salt_adv_split=0
39 40
   integer                   :: salt_adv_hor=1
   integer                   :: salt_adv_ver=1
kbk's avatar
kbk committed
41
   REALTYPE                  :: salt_AH=-_ONE_
42
   integer                   :: salt_check=0
43
   REALTYPE                  :: min_salt=_ZERO_,max_salt=40*_ONE_
gotm's avatar
gotm committed
44
!
kbk's avatar
kbk committed
45 46 47
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
gotm's avatar
gotm committed
48 49 50 51 52 53 54 55
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
56 57
! !IROUTINE: init_salinity - initialisation of salinity
! \label{sec-init-salinity}
gotm's avatar
gotm committed
58 59
!
! !INTERFACE:
60
   subroutine init_salinity()
gotm's avatar
gotm committed
61 62
!
! !DESCRIPTION:
63 64 65 66 67 68
!
! Here, the salinity equation is initialised. First, the namelist
! {\tt salt} is read from {\tt getm.inp}. Then, depending on the
! {\tt salt\_method}, the salinity field is read from a
! hotstart file ({\tt salt\_method}=0), initialised with a constant value
! ({\tt salt\_method}=1), initialised and interpolated
69
! with horizontally homogeneous
70 71 72 73 74 75 76
! salinity from a given salinity profile ({\tt salt\_method}=2),
! or read in and interpolated from a 3D netCDF field ({\tt salt\_method}=3).
! Finally, a number of sanity checks are performed for the chosen
! salinity advection schemes.
!
! Apart from this, there are various options for specific initial
! conditions which are selected by means of compiler options.
gotm's avatar
gotm committed
77 78
!
! !USES:
79 80
   use advection, only: J7
   use advection_3d, only: print_adv_settings_3d
gotm's avatar
gotm committed
81 82 83 84 85 86 87
   IMPLICIT NONE
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
88 89
   integer                   :: i,j,k,n
   integer                   :: status
90 91 92
   NAMELIST /salt/                                            &
            salt_method,salt_const,salt_file,                 &
            salt_format,salt_name,salt_field_no,              &
93
            salt_adv_split,salt_adv_hor,salt_adv_ver,salt_AH, &
94
            salt_check,min_salt,max_salt
gotm's avatar
gotm committed
95 96 97 98
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
99 100 101
!   integer, save :: Ncall = 0
!   Ncall = Ncall+1
!   write(debug,*) 'init_salinity() # ',Ncall
gotm's avatar
gotm committed
102 103 104 105 106
#endif

   LEVEL2 'init_salinity()'
   read(NAMLST,salt)

107
   call init_salinity_field()
108

109
!  Sanity checks for advection specifications
110 111 112
   LEVEL3 'Advection of salinity'
   if (salt_adv_hor .eq. J7) stop 'init_salinity: J7 not implemented yet'
   call print_adv_settings_3d(salt_adv_split,salt_adv_hor,salt_adv_ver,salt_AH)
113

114 115 116 117 118 119 120 121 122 123 124 125
   LEVEL3 'salt_check=',salt_check
   if (salt_check .ne. 0) then
      LEVEL4 'doing sanity check on salinity'
      LEVEL4 'min_salt=',min_salt
      LEVEL4 'max_salt=',max_salt
      if (salt_check .gt. 0) then
         LEVEL4 'out-of-bound values result in termination of program'
      end if
      if (salt_check .lt. 0) then
         LEVEL4 'out-of-bound values result in warnings only'
      end if

126
      if (salt_method .ne. 0) then
127 128 129 130
      call check_3d_fields(imin,jmin,imax,jmax,kmin,kmax,az, &
                           S,min_salt,max_salt,status)
      if (status .gt. 0) then
         if (salt_check .gt. 0) then
131
            call getm_error("init_salinity()", &
132 133 134
                            "out-of-bound values encountered")
         end if
         if (salt_check .lt. 0) then
135
            LEVEL1 'init_salinity(): ',status, &
136 137 138
                   ' out-of-bound values encountered'
         end if
      end if
139
      end if
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
   end if


#ifdef DEBUG
   write(debug,*) 'Leaving init_salinity()'
   write(debug,*)
#endif
   return
   end subroutine init_salinity
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_salinity_field - initialisation of the salinity field
155
! \label{sec-init-salinity-field}
156 157 158 159 160
!
! !INTERFACE:
   subroutine init_salinity_field()
!
! !DESCRIPTION:
161
! Initialisation of the salinity field as specified by salt\_method
162 163 164 165 166 167 168 169 170 171 172 173 174
! and exchange of the HALO zones
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
   integer                   :: i,j,k,n
175
   integer, parameter        :: nmax=10000
176 177 178 179 180 181 182 183 184 185
   REALTYPE                  :: zlev(nmax),prof(nmax)
   integer                   :: status
!EOP
!-------------------------------------------------------------------------
!BOC

   select case (salt_method)
      case(0)
         LEVEL3 'getting initial fields from hotstart'
      case(1)
Knut's avatar
Knut committed
186
         LEVEL3 'setting to constant value ',real(salt_const)
187 188 189 190
         forall(i=imin:imax,j=jmin:jmax, az(i,j) .ne. 0) &
                S(i,j,:) = salt_const
      case(2)
         LEVEL3 'using profile'
Knut's avatar
Knut committed
191
         LEVEL4 trim(salt_file)
192 193 194 195 196
         call read_profile(salt_file,nmax,zlev,prof,n)
         call ver_interpol(n,zlev,prof,imin,jmin,imax,jmax,kmax, &
                           az,H,hn,S)
      case(3)
         LEVEL3 'interpolating from 3D field'
Knut's avatar
Knut committed
197
         LEVEL4 trim(salt_file)
198
         call get_3d_field(salt_file,salt_name,salt_field_no,.true.,S)
199 200 201 202 203
      case default
         FATAL 'Not valid salt_method specified'
         stop 'init_salinity'
   end select

Knut's avatar
Knut committed
204 205
   S(:,:,0) = -9999._rk
   forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) S(i,j,:) = -9999._rk
Knut's avatar
Knut committed
206

207
   call update_3d_halo(S,S,az,imin,jmin,imax,jmax,kmax,D_TAG)
gotm's avatar
gotm committed
208
   call wait_halo(D_TAG)
kbk's avatar
kbk committed
209
   call mirror_bdy_3d(S,D_TAG)
gotm's avatar
gotm committed
210 211

#ifdef DEBUG
212
   write(debug,*) 'Leaving init_salinity_field()'
gotm's avatar
gotm committed
213 214 215
   write(debug,*)
#endif
   return
216
   end subroutine init_salinity_field
gotm's avatar
gotm committed
217 218 219 220 221
!EOC

!-----------------------------------------------------------------------
!BOP
!
222
! !IROUTINE:  do_salinity - salinity equation \label{sec-do-salinity}
gotm's avatar
gotm committed
223 224 225 226 227
!
! !INTERFACE:
   subroutine do_salinity(n)
!
! !DESCRIPTION:
228 229 230 231 232 233 234 235
!
! Here, one time step for the salinity equation is performed.
! First, preparations for the call to the advection schemes are
! made, i.e.\ calculating the necessary metric coefficients.
! After the call to the advection schemes, which actually perform
! the advection (and horizontal diffusion) step as an operational
! split step, the tri-diagonal matrix for calculating the
! new salinity by means of a semi-implicit central scheme for the
236
! vertical diffusion is set up.
237 238 239 240 241 242
! There are no source terms on the right hand sides.
! The subroutine is completed by solving the tri-diagonal linear
! equation by means of a tri-diagonal solver.
!
! Also here, there are some specific options for single test cases
! selected by compiler options.
gotm's avatar
gotm committed
243 244 245 246
!
! !USES:
   use advection_3d, only: do_advection_3d
   use variables_3d, only: dt,cnpar,hn,ho,nuh,uu,vv,ww,hun,hvn
247
   use domain,       only: imin,imax,jmin,jmax,kmax,az
248
   use parameters, only: avmols
249
   use getm_timers, only: tic, toc, TIM_SALT, TIM_MIXANALYSIS
250
   use variables_3d, only: do_numerical_analyses
251 252
   use variables_3d, only: nummix3d_S,nummix2d_S
   use variables_3d, only: phymix3d_S,phymix2d_S
bjb's avatar
bjb committed
253
!$ use omp_lib
gotm's avatar
gotm committed
254 255 256 257 258 259
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer, intent(in) :: n
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
260
   integer                   :: i,j,k,rc
bjb's avatar
bjb committed
261 262 263
   REALTYPE, POINTER         :: Res(:)
   REALTYPE, POINTER         :: auxn(:),auxo(:)
   REALTYPE, POINTER         :: a1(:),a2(:),a3(:),a4(:)
264
  REALTYPE                   :: S2(I3DFIELD)
265
  integer                    :: status
gotm's avatar
gotm committed
266 267 268 269 270 271 272 273
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'do_salinity() # ',Ncall
#endif
bjb's avatar
bjb committed
274
   call tic(TIM_SALT)
gotm's avatar
gotm committed
275

276 277
   do j=jmin-HALO,jmax+HALO
      do i=imin-HALO,imax+HALO
278
         if (az(i,j) .eq. 1) then
279 280 281 282
!           Note (KK): This would be the correct dilution if hn was
!                      corrected in start_macro. This also requires,
!                      that ho=hn is done in coordinates!
!            S(i,j,kmax) = S(i,j,kmax)*(_ONE_-fwf_int(i,j)/ho(i,j,kmax))
283
! Developers note:
284 285
!  The parentheses are set to minimize truncation errors for fwf_int=0
            S(i,j,kmax) = S(i,j,kmax)*            &
286 287
                          ( ho(i,j,kmax) / (ho(i,j,kmax)+fwf_int(i,j)) )

288 289 290 291
         end if
      end do
   end do

292
   if (do_numerical_analyses) then
293
      call toc(TIM_SALT)
294
      call tic(TIM_MIXANALYSIS)
295 296
! OMP-note: The following array-based line could be implemented
!    with OMP as a WORKSHARE construct. However, it would require a dedicated
297 298
!    PARALLEL region (as the various advection schemes have their own regions),
!    so the overhead of the contruct would be rather large.
299
      S2 = S**2
300 301
      call do_advection_3d(dt,S2,uu,vv,ww,hun,hvn,ho,hn,                           &
                           salt_adv_split,salt_adv_hor,salt_adv_ver,salt_AH,H_TAG)
302 303
      call toc(TIM_MIXANALYSIS)
      call tic(TIM_SALT)
304 305
   end if

306 307
   call do_advection_3d(dt,S,uu,vv,ww,hun,hvn,ho,hn,                            &
                        salt_adv_split,salt_adv_hor,salt_adv_ver,salt_AH,H_TAG)
gotm's avatar
gotm committed
308

309
   if (do_numerical_analyses) then
310
      call toc(TIM_SALT)
311
      call tic(TIM_MIXANALYSIS)
312

313
      call numerical_mixing(S2,S,nummix3d_S,nummix2d_S)
314 315 316

      call update_3d_halo(S,S,az,imin,jmin,imax,jmax,kmax,D_TAG)
      call wait_halo(D_TAG)
317
      call physical_mixing(S,salt_AH,avmols,phymix3d_S,phymix2d_S)
318

319
      call toc(TIM_MIXANALYSIS)
320
      call tic(TIM_SALT)
321 322
   end if

gotm's avatar
gotm committed
323

324
! OMP-NOTE: Pointers are used to for each thread to use its
bjb's avatar
bjb committed
325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
!           own work storage.
!$OMP PARALLEL DEFAULT(SHARED)                                         &
!$OMP    PRIVATE(i,j,k,rc)                                             &
!$OMP    PRIVATE(Res,auxn,auxo,a1,a2,a3,a4)

! Each thread allocates its own HEAP storage:
   allocate(Res(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'do_salinity: Error allocating memory (Res)'
   allocate(auxn(1:kmax-1),stat=rc)    ! work array
   if (rc /= 0) stop 'do_salinity: Error allocating memory (auxn)'
   allocate(auxo(1:kmax-1),stat=rc)    ! work array
   if (rc /= 0) stop 'do_salinity: Error allocating memory (auxo)'
   allocate(a1(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'do_salinity: Error allocating memory (a1)'
   allocate(a2(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'do_salinity: Error allocating memory (a2)'
   allocate(a3(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'do_salinity: Error allocating memory (a3)'
   allocate(a4(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'do_salinity: Error allocating memory (auxo)'

! Note: We do not need to initialize these work arrays.
347
!   Tested BJB 2009-09-25.
bjb's avatar
bjb committed
348 349 350 351 352



!  Advection and vertical diffusion and of salinity
!$OMP DO SCHEDULE(RUNTIME)
353 354
   do j=jmin,jmax
      do i=imin,imax
gotm's avatar
gotm committed
355 356 357 358
         if (az(i,j) .eq. 1) then
            if (kmax.gt.1) then
!     Auxilury terms, old and new time level,
               do k=1,kmax-1
359
                  auxo(k)=_TWO_*(1-cnpar)*dt*(nuh(i,j,k)+avmols)/ &
kbk's avatar
kbk committed
360
                             (hn(i,j,k+1)+hn(i,j,k))
361
                  auxn(k)=_TWO_*   cnpar *dt*(nuh(i,j,k)+avmols)/ &
kbk's avatar
kbk committed
362
                             (hn(i,j,k+1)+hn(i,j,k))
kbk's avatar
kbk committed
363
               end do
gotm's avatar
gotm committed
364 365

!        Matrix elements for surface layer
kbk's avatar
kbk committed
366 367 368 369
               k=kmax
               a1(k)=-auxn(k-1)
               a2(k)=hn(i,j,k)+auxn(k-1)
               a4(k)=S(i,j,k)*(hn(i,j,k)-auxo(k-1))+S(i,j,k-1)*auxo(k-1)
gotm's avatar
gotm committed
370 371 372

!        Matrix elements for inner layers
               do k=2,kmax-1
kbk's avatar
kbk committed
373 374 375 376 377 378
                  a3(k)=-auxn(k  )
                  a1(k)=-auxn(k-1)
                  a2(k)=hn(i,j,k)+auxn(k)+auxn(k-1)
                  a4(k)=S(i,j,k+1)*auxo(k)                           &
                       +S(i,j,k  )*(hn(i,j,k)-auxo(k)-auxo(k-1))     &
                       +S(i,j,k-1)*auxo(k-1)
gotm's avatar
gotm committed
379 380 381 382
               end do

!        Matrix elements for bottom layer
               k=1
kbk's avatar
kbk committed
383 384 385 386
               a3(k)=-auxn(k  )
               a2(k)=hn(i,j,k)+auxn(k)
               a4(k)=S(i,j,k+1)*auxo(k)                              &
                    +S(i,j,k  )*(hn(i,j,k)-auxo(k))
gotm's avatar
gotm committed
387 388 389

               call getm_tridiagonal(kmax,1,kmax,a1,a2,a3,a4,Res)

kbk's avatar
kbk committed
390 391 392
               do k=1,kmax
                  S(i,j,k)=Res(k)
               end do
gotm's avatar
gotm committed
393 394 395 396 397

            end if
         end if
      end do
   end do
bjb's avatar
bjb committed
398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414
!$OMP END DO

! Each thread must deallocate its own HEAP storage:
   deallocate(Res,stat=rc)
   if (rc /= 0) stop 'do_salinity: Error deallocating memory (Res)'
   deallocate(auxn,stat=rc)
   if (rc /= 0) stop 'do_salinity: Error deallocating memory (auxn)'
   deallocate(auxo,stat=rc)
   if (rc /= 0) stop 'do_salinity: Error deallocating memory (auxo)'
   deallocate(a1,stat=rc)
   if (rc /= 0) stop 'do_salinity: Error deallocating memory (a1)'
   deallocate(a2,stat=rc)
   if (rc /= 0) stop 'do_salinity: Error deallocating memory (a2)'
   deallocate(a3,stat=rc)
   if (rc /= 0) stop 'do_salinity: Error deallocating memory (a3)'
   deallocate(a4,stat=rc)
   if (rc /= 0) stop 'do_salinity: Error deallocating memory (a4)'
415

bjb's avatar
bjb committed
416
!$OMP END PARALLEL
417

418

419
   if (salt_check .ne. 0 .and. mod(n,abs(salt_check)) .eq. 0) then
420 421
      call check_3d_fields(imin,jmin,imax,jmax,kmin,kmax,az, &
                           S,min_salt,max_salt,status)
422 423 424 425 426 427 428 429 430 431 432 433
      if (status .gt. 0) then
         if (salt_check .gt. 0) then
            call getm_error("do_salinity()", &
                            "out-of-bound values encountered")
         end if
         if (salt_check .lt. 0) then
            LEVEL1 'do_salinity(): ',status, &
                   ' out-of-bound values encountered'
         end if
      end if
   end if

bjb's avatar
bjb committed
434
   call toc(TIM_SALT)
gotm's avatar
gotm committed
435 436 437 438 439 440 441 442 443 444 445 446 447 448 449
#ifdef DEBUG
   write(debug,*) 'Leaving do_salinity()'
   write(debug,*)
#endif
   return
   end subroutine do_salinity
!EOC

!-----------------------------------------------------------------------

   end module salinity

!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding               !
!-----------------------------------------------------------------------