salinity.F90 18.4 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
Knut's avatar
Knut committed
19
   use domain, only: imin,jmin,imax,jmax,kmax,H,az,dry_z
20
!KB   use get_field, only: get_3d_field
21
   use variables_2d, only: fwf_int
Knut's avatar
Knut committed
22
   use variables_3d, only: rk,S,hn,kmin
23
   use meteo, only: metforcing,met_method,nudge_sss,sss
Knut's avatar
Knut committed
24
   use meteo, only: METEO_CONST,METEO_FROMFILE,METEO_FROMEXT
25
   use halo_zones, only: update_3d_halo,wait_halo,D_TAG,H_TAG
gotm's avatar
gotm committed
26
   IMPLICIT NONE
27 28 29 30 31 32 33 34 35 36 37 38

   interface
      subroutine tracer_diffusion(f,hn,AH_method,AH_const,AH_Prt,AH_stirr_const, &
                                  phymix)
         use domain, only: imin,imax,jmin,jmax,kmax
         REALTYPE,intent(in)           :: hn(I3DFIELD)
         integer,intent(in)            :: AH_method
         REALTYPE,intent(in)           :: AH_const,AH_Prt,AH_stirr_const
         REALTYPE,intent(inout)        :: f(I3DFIELD)
         REALTYPE,dimension(:,:,:),pointer,intent(out),optional :: phymix
      end subroutine tracer_diffusion
   end interface
gotm's avatar
gotm committed
39 40 41 42
!
   private
!
! !PUBLIC DATA MEMBERS:
43
   public init_salinity, do_salinity, init_salinity_field
gotm's avatar
gotm committed
44 45
!
! !PRIVATE DATA MEMBERS:
kbk's avatar
kbk committed
46 47
   integer                   :: salt_method=1,salt_format=2
   character(len=PATH_MAX)   :: salt_file="t_and_s.nc"
48
   integer                   :: salt_field_no=1
kbk's avatar
kbk committed
49
   character(len=32)         :: salt_name='salt'
50
   REALTYPE                  :: salt_const=35*_ONE_
51
   integer                   :: salt_adv_split=0
52 53
   integer                   :: salt_adv_hor=1
   integer                   :: salt_adv_ver=1
Knut's avatar
Knut committed
54
   integer                   :: salt_AH_method=0
Knut's avatar
Knut committed
55 56
   REALTYPE                  :: salt_AH_const=1.1d-9
   REALTYPE                  :: salt_AH_Prt=_TWO_
Knut's avatar
Knut committed
57
   REALTYPE                  :: salt_AH_stirr_const=_ONE_
Knut's avatar
Knut committed
58
   REALTYPE                  :: sss_nudging_time=-_ONE_
59
   integer                   :: salt_check=0
60
   REALTYPE                  :: min_salt=_ZERO_,max_salt=40*_ONE_
gotm's avatar
gotm committed
61
!
kbk's avatar
kbk committed
62 63 64
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
gotm's avatar
gotm committed
65 66 67 68 69 70 71 72
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
73 74
! !IROUTINE: init_salinity - initialisation of salinity
! \label{sec-init-salinity}
gotm's avatar
gotm committed
75 76
!
! !INTERFACE:
Knut's avatar
Knut committed
77
   subroutine init_salinity()
gotm's avatar
gotm committed
78 79
!
! !DESCRIPTION:
80 81 82 83 84 85
!
! 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
86
! with horizontally homogeneous
87 88 89 90 91 92 93
! 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
94 95
!
! !USES:
Knut's avatar
Knut committed
96
   use advection, only: J7
Knut's avatar
Knut committed
97
   use advection_3d, only: print_adv_settings_3d
98 99 100
   use variables_3d, only: deformC_3d,deformX_3d,deformUV_3d,calc_stirr
   use m2d, only: Am_method,AM_LES
   use les, only: les_mode,LES_TRACER,LES_BOTH
101
   use meteo, only: sss_const
gotm's avatar
gotm committed
102 103 104 105 106 107 108
   IMPLICIT NONE
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
Knut's avatar
Knut committed
109
   integer                   :: i,j,k,n,rc
110
   integer                   :: status
111 112 113
   NAMELIST /salt/                                            &
            salt_method,salt_const,salt_file,                 &
            salt_format,salt_name,salt_field_no,              &
Knut's avatar
Knut committed
114
            salt_adv_split,salt_adv_hor,salt_adv_ver,         &
Knut's avatar
Knut committed
115
            salt_AH_method,salt_AH_const,salt_AH_Prt,         &
Knut's avatar
Knut committed
116 117
            salt_AH_stirr_const,sss_nudging_time,             &
            salt_check,min_salt,max_salt
gotm's avatar
gotm committed
118 119 120 121
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
122 123 124
!   integer, save :: Ncall = 0
!   Ncall = Ncall+1
!   write(debug,*) 'init_salinity() # ',Ncall
gotm's avatar
gotm committed
125 126 127 128 129
#endif

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

130
   call init_salinity_field()
131

132
!  Sanity checks for advection specifications
Knut's avatar
Knut committed
133
   LEVEL3 'Advection of salinity'
134
   if (salt_adv_hor .eq. J7) stop 'init_salinity: J7 not implemented yet'
Knut's avatar
Knut committed
135
   call print_adv_settings_3d(salt_adv_split,salt_adv_hor,salt_adv_ver,_ZERO_)
136

Knut's avatar
Knut committed
137 138 139 140 141 142 143 144 145
   select case (salt_AH_method)
      case(0)
         LEVEL3 'salt_AH_method=0 -> horizontal diffusion of salt disabled'
      case(1)
         LEVEL3 'salt_AH_method=1 -> Using constant horizontal diffusivity of salt'
         if (salt_AH_const .lt. _ZERO_) then
              call getm_error("init_salinity()", &
                         "Constant horizontal diffusivity of salt <0");
         end if
Knut's avatar
Knut committed
146
         LEVEL4 real(salt_AH_const)
Knut's avatar
Knut committed
147 148
      case(2)
         LEVEL3 'salt_AH_method=2 -> using LES parameterisation'
Knut's avatar
Knut committed
149
         LEVEL4 'turbulent Prandtl number: ',real(salt_AH_Prt)
150 151 152 153 154 155 156 157
         deformC_3d =.true.
         deformX_3d =.true.
         deformUV_3d=.true.
         if (Am_method .eq. AM_LES) then
            les_mode = LES_BOTH
         else
            les_mode = LES_TRACER
         end if
Knut's avatar
Knut committed
158 159 160 161 162 163
      case(3)
         LEVEL3 'salt_AH_method=3 -> SGS stirring parameterisation'
         if (salt_AH_stirr_const .lt. _ZERO_) then
              call getm_error("init_salinity()", &
                         "salt_AH_stirr_const <0");
         end if
Knut's avatar
Knut committed
164
         LEVEL4 'stirring constant: ',real(salt_AH_stirr_const)
165 166 167 168
         deformC_3d =.true.
         deformX_3d =.true.
         deformUV_3d=.true.
         calc_stirr=.true.
Knut's avatar
Knut committed
169 170 171 172 173
      case default
         call getm_error("init_salinity()", &
                         "A non valid salt_AH_method has been chosen");
   end select

Knut's avatar
Knut committed
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
   if (metforcing .and. sss_nudging_time.gt._ZERO_) then
      nudge_sss = .True.
      LEVEL3 'nudging of SSS enabled'
      LEVEL4 'sss_nudging_time=',real(sss_nudging_time)
      allocate(sss(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_salinity: Error allocating memory (sss)'
      select case(met_method)
         case (METEO_CONST)
            if (sss_const .gt. _ZERO_) then
               LEVEL4 'constant sss=',real(sss_const)
               sss = sss_const
            else
               call getm_error("init_salinity()", &
                               "non-positive sss_const")
            end if
         case (METEO_FROMFILE,METEO_FROMEXT)
            LEVEL4 'sss read'
      end select
   end if

194 195 196 197 198 199 200 201 202 203 204 205
   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

206
      if (salt_method .ne. 0) then
207 208 209 210
      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
211
            call getm_error("init_salinity()", &
212 213 214
                            "out-of-bound values encountered")
         end if
         if (salt_check .lt. 0) then
215
            LEVEL1 'init_salinity(): ',status, &
216 217 218
                   ' out-of-bound values encountered'
         end if
      end if
219
      end if
220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
   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
235
! \label{sec-init-salinity-field}
236 237 238 239 240
!
! !INTERFACE:
   subroutine init_salinity_field()
!
! !DESCRIPTION:
241
! Initialisation of the salinity field as specified by salt\_method
242 243 244 245 246 247 248 249 250 251 252 253 254
! and exchange of the HALO zones
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
   integer                   :: i,j,k,n
255
   integer, parameter        :: nmax=10000
256 257 258 259 260 261 262 263 264 265
   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
266
         LEVEL3 'setting to constant value ',real(salt_const)
267 268 269 270
         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
271
         LEVEL4 trim(salt_file)
272 273 274 275 276
         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
277
         LEVEL4 trim(salt_file)
278
         call get_3d_field(salt_file,salt_name,salt_field_no,.true.,S)
279 280 281 282 283
      case default
         FATAL 'Not valid salt_method specified'
         stop 'init_salinity'
   end select

Knut's avatar
Knut committed
284 285
   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
286

287
   call update_3d_halo(S,S,az,imin,jmin,imax,jmax,kmax,D_TAG)
gotm's avatar
gotm committed
288
   call wait_halo(D_TAG)
kbk's avatar
kbk committed
289
   call mirror_bdy_3d(S,D_TAG)
gotm's avatar
gotm committed
290 291

#ifdef DEBUG
292
   write(debug,*) 'Leaving init_salinity_field()'
gotm's avatar
gotm committed
293 294 295
   write(debug,*)
#endif
   return
296
   end subroutine init_salinity_field
gotm's avatar
gotm committed
297 298 299 300 301
!EOC

!-----------------------------------------------------------------------
!BOP
!
302
! !IROUTINE:  do_salinity - salinity equation \label{sec-do-salinity}
gotm's avatar
gotm committed
303 304 305 306 307
!
! !INTERFACE:
   subroutine do_salinity(n)
!
! !DESCRIPTION:
308 309 310 311 312 313 314 315
!
! 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
316
! vertical diffusion is set up.
317 318 319 320 321 322
! 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
323 324 325 326
!
! !USES:
   use advection_3d, only: do_advection_3d
   use variables_3d, only: dt,cnpar,hn,ho,nuh,uu,vv,ww,hun,hvn
Knut's avatar
Knut committed
327
   use domain,       only: imin,imax,jmin,jmax,kmax,az
328
   use parameters, only: avmols
Knut's avatar
Knut committed
329
   use getm_timers, only: tic,toc,TIM_SALT,TIM_SALTH,TIM_MIXANALYSIS
Knut's avatar
Knut committed
330
   use variables_3d, only: do_numerical_analyses_3d
331
   use variables_3d, only: nummix_S,nummix_S_old,nummix_S_int
Knut's avatar
Knut committed
332
   use variables_3d, only: phymix_S,phymix_S_int
bjb's avatar
bjb committed
333
!$ use omp_lib
gotm's avatar
gotm committed
334 335 336 337 338 339
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer, intent(in) :: n
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
340
   integer                   :: i,j,k,rc
bjb's avatar
bjb committed
341 342 343
   REALTYPE, POINTER         :: Res(:)
   REALTYPE, POINTER         :: auxn(:),auxo(:)
   REALTYPE, POINTER         :: a1(:),a2(:),a3(:),a4(:)
344
#ifdef _NUMERICAL_ANALYSES_OLD_
345
  REALTYPE                   :: S2(I3DFIELD)
346
#endif
347
  integer                    :: status
gotm's avatar
gotm committed
348 349 350 351 352 353 354 355
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'do_salinity() # ',Ncall
#endif
bjb's avatar
bjb committed
356
   call tic(TIM_SALT)
gotm's avatar
gotm committed
357

358 359
   do j=jmin-HALO,jmax+HALO
      do i=imin-HALO,imax+HALO
360
         if (az(i,j) .eq. 1) then
361
            S(i,j,kmax) = S(i,j,kmax)*(_ONE_-fwf_int(i,j)/ho(i,j,kmax))
362 363 364 365
         end if
      end do
   end do

366
#ifdef _NUMERICAL_ANALYSES_OLD_
Knut's avatar
Knut committed
367
   if (do_numerical_analyses_3d) then
368
      call toc(TIM_SALT)
369
      call tic(TIM_MIXANALYSIS)
370 371
! OMP-note: The following array-based line could be implemented
!    with OMP as a WORKSHARE construct. However, it would require a dedicated
372 373
!    PARALLEL region (as the various advection schemes have their own regions),
!    so the overhead of the contruct would be rather large.
374
      S2 = S**2
375
      call do_advection_3d(dt,S2,uu,vv,ww,hun,hvn,ho,hn,                           &
Knut's avatar
Knut committed
376
                           salt_adv_split,salt_adv_hor,salt_adv_ver,_ZERO_,H_TAG)
377 378
      call toc(TIM_MIXANALYSIS)
      call tic(TIM_SALT)
379
   end if
380
#endif
381

Knut's avatar
Knut committed
382 383
   call do_advection_3d(dt,S,uu,vv,ww,hun,hvn,ho,hn,                           &
                        salt_adv_split,salt_adv_hor,salt_adv_ver,_ZERO_,H_TAG, &
384
                        nvd=nummix_S)
gotm's avatar
gotm committed
385

Knut's avatar
Knut committed
386
#ifdef _NUMERICAL_ANALYSES_OLD_
Knut's avatar
Knut committed
387
   if (do_numerical_analyses_3d) then
388
      call toc(TIM_SALT)
389
      call tic(TIM_MIXANALYSIS)
390
      call numerical_mixing(S2,S,nummix_S_old,nummix_S_int)
391
      call toc(TIM_MIXANALYSIS)
392
      call tic(TIM_SALT)
393
   end if
394
#endif
395

Knut's avatar
Knut committed
396 397 398 399 400 401 402
   if (salt_AH_method .gt. 0) then
!     S is not halo updated after advection
      call tic(TIM_SALTH)
      call update_3d_halo(S,S,az,imin,jmin,imax,jmax,kmax,D_TAG)
      call wait_halo(D_TAG)
      call toc(TIM_SALTH)

403 404
      call tracer_diffusion(S,hn,salt_AH_method,salt_AH_const,salt_AH_Prt,salt_AH_stirr_const, &
                            phymix_S)
Knut's avatar
Knut committed
405 406
   end if

407
   if (do_numerical_analyses_3d) then
Knut's avatar
Knut committed
408 409
      call toc(TIM_SALT)
      call tic(TIM_MIXANALYSIS)
Knut's avatar
Knut committed
410
      call physical_mixing(S,avmols,phymix_S,phymix_S_int,salt_AH_method)
Knut's avatar
Knut committed
411 412
      call toc(TIM_MIXANALYSIS)
      call tic(TIM_SALT)
Knut's avatar
Knut committed
413 414
   end if

gotm's avatar
gotm committed
415

416
! OMP-NOTE: Pointers are used to for each thread to use its
bjb's avatar
bjb committed
417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438
!           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.
439
!   Tested BJB 2009-09-25.
bjb's avatar
bjb committed
440 441 442 443 444



!  Advection and vertical diffusion and of salinity
!$OMP DO SCHEDULE(RUNTIME)
445 446
   do j=jmin,jmax
      do i=imin,imax
gotm's avatar
gotm committed
447 448 449 450
         if (az(i,j) .eq. 1) then
            if (kmax.gt.1) then
!     Auxilury terms, old and new time level,
               do k=1,kmax-1
451
                  auxo(k)=_TWO_*(1-cnpar)*dt*(nuh(i,j,k)+avmols)/ &
kbk's avatar
kbk committed
452
                             (hn(i,j,k+1)+hn(i,j,k))
453
                  auxn(k)=_TWO_*   cnpar *dt*(nuh(i,j,k)+avmols)/ &
kbk's avatar
kbk committed
454
                             (hn(i,j,k+1)+hn(i,j,k))
kbk's avatar
kbk committed
455
               end do
gotm's avatar
gotm committed
456 457

!        Matrix elements for surface layer
kbk's avatar
kbk committed
458 459 460 461
               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)
Knut's avatar
Knut committed
462 463 464 465 466 467 468
               if (nudge_sss) then
!                 implicit nudging
                  a2(kmax) = a2(kmax) + dry_z(i,j)*hn(i,j,kmax)*dt/sss_nudging_time
                  a4(kmax) = a4(kmax) + dry_z(i,j)*hn(i,j,kmax)*dt/sss_nudging_time*sss(i,j)
!                 explicit nudging
                  !a4(kmax) = a4(kmax) - dry_z(i,j)*dt*hn(i,j,kmax)*(S(i,j,kmax)-sss(i,j))/sss_nudging_time
               end if
gotm's avatar
gotm committed
469 470 471

!        Matrix elements for inner layers
               do k=2,kmax-1
kbk's avatar
kbk committed
472 473 474 475 476 477
                  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
478 479 480 481
               end do

!        Matrix elements for bottom layer
               k=1
kbk's avatar
kbk committed
482 483 484 485
               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
486 487 488

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

kbk's avatar
kbk committed
489 490 491
               do k=1,kmax
                  S(i,j,k)=Res(k)
               end do
gotm's avatar
gotm committed
492 493 494 495 496

            end if
         end if
      end do
   end do
bjb's avatar
bjb committed
497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513
!$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)'
514

bjb's avatar
bjb committed
515
!$OMP END PARALLEL
516

517

518
   if (salt_check .ne. 0 .and. mod(n,abs(salt_check)) .eq. 0) then
519 520
      call check_3d_fields(imin,jmin,imax,jmax,kmin,kmax,az, &
                           S,min_salt,max_salt,status)
521 522 523 524 525 526 527 528 529 530 531 532
      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
533
   call toc(TIM_SALT)
gotm's avatar
gotm committed
534 535 536 537 538 539 540 541 542 543 544 545 546 547 548
#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               !
!-----------------------------------------------------------------------