salinity.F90 13.9 KB
Newer Older
kbk's avatar
kbk committed
1
!$Id: salinity.F90,v 1.20 2006-03-01 15:54:08 kbk Exp $
gotm's avatar
gotm committed
2 3 4 5 6 7 8 9 10 11
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE:  Salinity
!
! !INTERFACE:
   module salinity
!
! !DESCRIPTION:
12 13 14 15 16
!
! 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
17 18
!
! !USES:
19
   use exceptions
20 21 22 23
   use domain, only: imin,jmin,imax,jmax,ioff,joff
#ifdef HAIDVOGEL_TEST
   use domain, only: iextr,jextr
#endif
gotm's avatar
gotm committed
24
   use domain, only: iimin,jjmin,iimax,jjmax,kmax
25
   use domain, only: H,az
26
   use variables_3d, only: S,hn,adv_schemes
27
   use halo_zones, only: update_3d_halo,wait_halo,D_TAG
gotm's avatar
gotm committed
28 29 30 31 32 33 34 35
   IMPLICIT NONE
!
   private
!
! !PUBLIC DATA MEMBERS:
   public init_salinity, do_salinity
!
! !PRIVATE DATA MEMBERS:
kbk's avatar
kbk committed
36 37 38 39
   integer                   :: salt_method=1,salt_format=2
   character(len=PATH_MAX)   :: salt_file="t_and_s.nc"
   character(len=32)         :: salt_name='salt'
   REALTYPE                  :: salt_const=35.
40 41
   integer                   :: salt_hor_adv=1,salt_ver_adv=1
   integer                   :: salt_adv_split=0
kbk's avatar
kbk committed
42
   REALTYPE                  :: salt_AH=-_ONE_
gotm's avatar
gotm committed
43
!
kbk's avatar
kbk committed
44 45 46
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
gotm's avatar
gotm committed
47 48 49 50 51 52 53 54 55
! !LOCAL VARIABLES:
!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 60 61 62
!
! !INTERFACE:
   subroutine init_salinity(adv_method)
!
! !DESCRIPTION:
63 64 65 66 67 68 69 70 71 72 73 74 75 76
!
! 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
! with horizontally homogeneous 
! 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:
kbk's avatar
kbk committed
79
#ifdef FRESHWATER_LENSE_TEST
80 81
   use domain, only: dx
#endif
gotm's avatar
gotm committed
82 83 84
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
85
   integer, intent(in)                 :: adv_method
gotm's avatar
gotm committed
86 87 88 89 90 91
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
92
  integer                    :: i,j,k,n
gotm's avatar
gotm committed
93
#ifdef PECS_TEST
kbk's avatar
kbk committed
94
   integer                   :: cc(1:30)
gotm's avatar
gotm committed
95
#endif
kbk's avatar
kbk committed
96
#ifdef FRESHWATER_LENSE_TEST
kbk's avatar
kbk committed
97 98
   INTEGER                   :: ic,jc
   REALTYPE                  :: dist
gotm's avatar
gotm committed
99
#endif
kbk's avatar
kbk committed
100 101
   integer, parameter        :: nmax=100
   REALTYPE                  :: zlev(nmax),prof(nmax)
102
   integer                   :: salt_field_no=1
kbk's avatar
kbk committed
103 104
   NAMELIST /salt/                                          &
            salt_method,salt_const,salt_file,               &
105
            salt_format,salt_name,salt_field_no,            &
106
            salt_hor_adv,salt_ver_adv,salt_adv_split,salt_AH
gotm's avatar
gotm committed
107 108 109 110 111 112 113 114 115
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'init_salinity() # ',Ncall
#endif

kbk's avatar
kbk committed
116
#ifdef NS_FRESHWATER_LENSE_TEST
117
salt_field_no=1
gotm's avatar
gotm committed
118 119
#endif
#ifdef MED_15X15MINS_TEST
120
salt_field_no=1
gotm's avatar
gotm committed
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
#endif

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

   select case (salt_method)
      case(0)
         LEVEL3 'getting initial fields from hotstart'
      case(1)
         LEVEL3 'setting to constant value'
         forall(i=iimin:iimax,j=jjmin:jjmax, az(i,j) .ne. 0) &
                S(i,j,:) = salt_const
      case(2)
         LEVEL3 'using profile'
         call read_profile(salt_file,nmax,zlev,prof,n)
kbk's avatar
kbk committed
136 137
         call ver_interpol(n,zlev,prof,imin,jmin,imax,jmax,az,H,  &
                           iimin,jjmin,iimax,jjmax,kmax,hn,S)
gotm's avatar
gotm committed
138 139
      case(3)
         LEVEL3 'interpolating from 3D field'
140
         call get_field(salt_file,salt_name,salt_field_no,S)
gotm's avatar
gotm committed
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
#ifdef SALTWEDGE_TEST
      case(4)
      !I need this here for hotstart of salinity!!!
         LEVEL3 'initializing with #ifdef SALTWEDGE'
         S =  _ZERO_
         do i=1,100
            do k=1,kmax
               S(i,2,k)=30.*(1.- tanh(float(i-1)*0.05))
            end do      
         end do   
#endif
      case default
         FATAL 'Not valid salt_method specified'
         stop 'init_salinity'
   end select

157 158 159 160 161 162 163 164 165 166
!  Sanity checks for advection specifications
   LEVEL3 'salt_hor_adv=   ',salt_hor_adv
   LEVEL3 'salt_ver_adv=   ',salt_ver_adv
   LEVEL3 'salt_adv_split= ',salt_adv_split
   if(salt_hor_adv .eq. 1) then
      salt_adv_split=-1
      if(salt_ver_adv .ne. 1) then
         LEVEL3 "setting salt_ver_adv to 1 - since salt_hor_adv is 1"
         salt_ver_adv=1
      end if
167
   end if
kbk's avatar
kbk committed
168 169
   LEVEL3 "horizontal: ",trim(adv_schemes(salt_hor_adv))
   LEVEL3 "vertical:   ",trim(adv_schemes(salt_ver_adv))
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

   select case (salt_adv_split)
      case (-1)
      case (0)
         select case (salt_hor_adv)
            case (2,3,4,5,6)
            case default
               call getm_error("init_3d()", &
                    "salt_adv_split=0: salt_hor_adv not valid (2-6)")
         end select
         select case (salt_ver_adv)
            case (2,3,4,5,6)
            case default
               call getm_error("init_3d()", &
                    "salt_adv_split=0: salt_ver_adv not valid (2-6)")
         end select
         LEVEL3 "1D split --> full u, full v, full w"
      case (1)
         select case (salt_hor_adv)
            case (2,3,4,5,6)
            case default
               call getm_error("init_3d()", &
                    "salt_adv_split=1: salt_hor_adv not valid (2-6)")
         end select
         select case (salt_ver_adv)
            case (2,3,4,5,6)
            case default
               call getm_error("init_3d()", &
                    "salt_adv_split=1: salt_ver_adv not valid (2-6)")
         end select
         LEVEL3 "1D split --> half u, half v, full w, half v, half u"
      case (2)
         select case (salt_hor_adv)
            case (2,7)
            case default
               call getm_error("init_3d()", &
                    "salt_adv_split=2: salt_hor_adv not valid (2,7)")
         end select
         select case (salt_ver_adv)
            case (2,3,4,5,6)
            case default
               call getm_error("init_3d()", &
                    "salt_adv_split=2: salt_ver_adv not valid (2-6)")
         end select
         LEVEL3 "2D-hor, 1D-vert split --> full uv, full w"
      case default
   end select
217

kbk's avatar
kbk committed
218
#ifdef FRESHWATER_LENSE_TEST
gotm's avatar
gotm committed
219 220 221 222 223 224
   S=34.85
   ic=nint(iimax/2.)
   jc=nint(jjmax/2.)
   do i=1,iimax
      do j=1,jjmax
         do k=kmax/2,kmax
225
            dist=sqrt((float(i)-float(ic))**2+(float(j)-float(jc))**2)*dx
kbk's avatar
kbk committed
226
            if (dist.le.3000.) then
gotm's avatar
gotm committed
227
               S(i,j,k)=1.1*(dist/(3000.))**8+33.75
kbk's avatar
kbk committed
228 229 230 231
            else
               S(i,j,k)=34.85
            end if
         end do
gotm's avatar
gotm committed
232 233 234
      end do
   end do
#endif
hb's avatar
hb committed
235 236 237 238 239 240 241 242 243 244 245 246
#ifdef ARKONA_TEST
   do i=100,135
      do j=256,257
         if (az(i,j).ge.1) S(i,j,0:kmax) = 25.
      end do
   end do
   do i=26,27
      do j=77,100
         S(i,j,0:kmax) = 8.
      end do
   end do
#endif
247 248
#ifdef SLOPE_TEST
   do i=81,82
249
      do j=42,43
250 251 252 253
      S(i,j,0:kmax) = 25.
      end do
   end do
#endif
lars's avatar
lars committed
254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271
#ifdef BALTIC_SLICE_TEST
  j=2
   if (iimax.eq.102) then
      do i=2,21
        S(i,j,0:kmax) = 25.
      end do
   end if
   if (iimax.eq.302) then
      do i=2,61
        S(i,j,0:kmax) = 25.
      end do
   end if
   if (iimax.eq.902) then
      do i=2,181
        S(i,j,0:kmax) = 25.
      end do
   end if
#endif
gotm's avatar
gotm committed
272
#ifdef HAIDVOGEL_TEST
273
   STDERR 'HAIDVOGEL: salinity= ',iimin,iimax,i+ioff,iextr/2
274 275 276 277 278 279 280 281
   do i=iimin-1,iimax+1
      if(i+ioff .le. iextr/2) then
         S(i,jjmin-1:jjmax+1,0:kmax) = 6.4102564
         S(i,jjmin-1:jjmax+1,0:kmax) = 5.
      else
         S(i,jjmin-1:jjmax+1,0:kmax) = 0.
      end if
   end do
gotm's avatar
gotm committed
282 283 284 285 286 287 288 289
#endif
!#else
!#ifdef PECS_TEST
!   S = 10.
!   do i=1,160
!      read(98,*) cc(1:30)
!      do j=1,30
!         if (cc(j).eq.1) then
kbk's avatar
kbk committed
290 291
!            S(i,j+1,1)=20.
!         end if
gotm's avatar
gotm committed
292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
!      end do
!   end do
!#endif
!#endif

   call update_3d_halo(S,S,az,iimin,jjmin,iimax,jjmax,kmax,D_TAG)
   call wait_halo(D_TAG)

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

!-----------------------------------------------------------------------
!BOP
!
311
! !IROUTINE:  do_salinity - salinity equation \label{sec-do-salinity}
gotm's avatar
gotm committed
312 313 314 315 316
!
! !INTERFACE:
   subroutine do_salinity(n)
!
! !DESCRIPTION:
317 318 319 320 321 322 323 324 325 326 327 328 329 330 331
!
! 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
! vertical diffusion is set up. 
! 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
332 333 334 335 336 337 338 339 340 341
!
! !USES:
   use advection_3d, only: do_advection_3d
   use variables_3d, only: dt,cnpar,hn,ho,nuh,uu,vv,ww,hun,hvn
   use domain,       only: iimin,iimax,jjmin,jjmax,kmax,az,au,av
#if defined(SPHERICAL) || defined(CURVILINEAR)
   use domain, only: dxu,dxv,dyu,dyv,arcd1
#else
   use domain, only: dx,dy,ard1
#endif
342
   use parameters, only: avmols
gotm's avatar
gotm committed
343 344 345 346 347 348 349 350 351 352
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer, intent(in) :: n
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
353 354 355 356 357
   integer                   :: i,j,k,rc
   REALTYPE                  :: Res(0:kmax)
   REALTYPE                  :: auxn(1:kmax-1),auxo(1:kmax-1)
   REALTYPE                  :: a1(0:kmax),a2(0:kmax)
   REALTYPE                  :: a3(0:kmax),a4(0:kmax)
kbk's avatar
kbk committed
358
#ifdef FRESHWATER_LENSE_TEST
359
   REALTYPE                  :: SRelax
gotm's avatar
gotm committed
360 361
#endif
#ifdef SALTWEDGE_TEST
kbk's avatar
kbk committed
362
   REALTYPE                  :: SRelax,kk
gotm's avatar
gotm committed
363
#endif
kbk's avatar
kbk committed
364 365 366
  REALTYPE                   :: delxu(I2DFIELD),delxv(I2DFIELD)
  REALTYPE                   :: delyu(I2DFIELD),delyv(I2DFIELD)
  REALTYPE                   :: area_inv(I2DFIELD)
gotm's avatar
gotm committed
367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'do_salinity() # ',Ncall
#endif

#if defined(SPHERICAL) || defined(CURVILINEAR)
   delxu=dxu
   delxv=dxv
   delyu=dyu
   delyv=dyv
   area_inv=arcd1
#else
   delxu=dx
   delxv=dx
   delyu=dy
   delyv=dy
   area_inv=ard1
#endif
   call do_advection_3d(dt,S,uu,vv,ww,hun,hvn,ho,hn,    &
                        delxu,delxv,delyu,delyv,area_inv,az,au,av,   &
391
                        salt_hor_adv,salt_ver_adv,salt_adv_split,salt_AH)
gotm's avatar
gotm committed
392 393 394 395 396

#ifdef PECS_TEST
   S(iimin:iimin,jjmin:jjmax,1:kmax)=10.
   S(iimax:iimax,jjmin:jjmax,1:kmax)=10.
#endif
kbk's avatar
kbk committed
397

kbk's avatar
kbk committed
398
#ifdef FRESHWATER_LENSE_TEST
gotm's avatar
gotm committed
399
   SRelax=34.85
400 401 402 403
   S(iimin:iimin+3,jjmin:jjmax,1:kmax)=SRelax
   S(iimax-3:iimax,jjmin:jjmax,1:kmax)=SRelax
   S(iimin:iimax,jjmin:jjmin+3,1:kmax)=SRelax
   S(iimin:iimax,jjmax-3:jjmax,1:kmax)=SRelax
gotm's avatar
gotm committed
404
#endif
kbk's avatar
kbk committed
405

gotm's avatar
gotm committed
406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425
#ifdef SALTWEDGE_TEST
   SRelax=30.
   j=2
   do k=1,kmax
      do i=1,100
         kk=  1.- tanh(float(i-1)*0.05)
         S(i,j,k)=(1.-kk)*S(i,j,k)+kk*SRelax
      end do
   end do
   S(imax-1,2,:)=0. !river
#endif 

!  Advection and vertical diffusion and of salinity

   do j=jjmin,jjmax
      do i=iimin,iimax
         if (az(i,j) .eq. 1) then
            if (kmax.gt.1) then
!     Auxilury terms, old and new time level,
               do k=1,kmax-1
426
                  auxo(k)=2.*(1-cnpar)*dt*(nuh(i,j,k)+avmols)/ &
kbk's avatar
kbk committed
427
                             (hn(i,j,k+1)+hn(i,j,k))
428
                  auxn(k)=2.*   cnpar *dt*(nuh(i,j,k)+avmols)/ &
kbk's avatar
kbk committed
429
                             (hn(i,j,k+1)+hn(i,j,k))
kbk's avatar
kbk committed
430
               end do
gotm's avatar
gotm committed
431 432

!        Matrix elements for surface layer
kbk's avatar
kbk committed
433 434 435 436
               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
437 438 439

!        Matrix elements for inner layers
               do k=2,kmax-1
kbk's avatar
kbk committed
440 441 442 443 444 445
                  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
446 447 448 449
               end do

!        Matrix elements for bottom layer
               k=1
kbk's avatar
kbk committed
450 451 452 453
               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
454 455 456

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

kbk's avatar
kbk committed
457 458 459
               do k=1,kmax
                  S(i,j,k)=Res(k)
               end do
gotm's avatar
gotm committed
460 461 462 463 464

            end if
         end if
      end do
   end do
465

hb's avatar
hb committed
466 467 468 469 470 471 472 473 474 475 476 477 478
#ifdef ARKONA_TEST
   do i=100,135
      do j=256,257
         if (az(i,j).ge.1) S(i,j,0:kmax) = 25.
      end do
   end do
   do i=26,27
      do j=77,100
      S(i,j,0:kmax) = 8.
      end do
   end do
#endif

gotm's avatar
gotm committed
479

480 481
#ifdef SLOPE_TEST
   do i=81,82
482
      do j=42,43
483 484 485 486 487
      S(i,j,0:kmax) = 25.
      end do
   end do
#endif

gotm's avatar
gotm committed
488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505
   call update_3d_halo(S,S,az,iimin,jjmin,iimax,jjmax,kmax,D_TAG)
   call wait_halo(D_TAG)

#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               !
!-----------------------------------------------------------------------