rivers.F90 16 KB
Newer Older
kbk's avatar
kbk committed
1
!$Id: rivers.F90,v 1.11 2006-03-01 15:54:08 kbk Exp $
gotm's avatar
gotm committed
2 3 4 5
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
6
! !MODULE:  rivers \label{sec-rivers}
gotm's avatar
gotm committed
7 8 9 10 11
!
! !INTERFACE:
   module rivers
!
! !DESCRIPTION:
12
!
gotm's avatar
gotm committed
13
!  This module includes support for river input. Rivers are treated the same
14 15 16 17 18 19 20 21 22 23 24
!  way as meteorology, i.e.\ as external module to the hydrodynamic model 
!  itself.
!  The module follows the same scheme as all other modules, i.e.\
!  {\tt init\_rivers}
!  sets up necessary information, and {\tt do\_rivers} updates
!  the relevant variables.
!  {\tt do\_river} is called in {\tt getm/integration.F90} 
!  between the {\tt 2d} and {\tt 3d} routines as it only
!  updates the sea surface elevation (in {\tt 2d}) and sea surface elevation,
!  and
!  optionally salinity and temperature (in {\tt 3d}). 
gotm's avatar
gotm committed
25 26 27 28 29
!  At present the momentum of the river water is not include, the model
!  however has a direct response to the river water because of the
!  pressure gradient introduced.
!
! !USES:
30
   use domain, only: iimin,jjmin,iimax,jjmax,ioff,joff
gotm's avatar
gotm committed
31 32 33 34 35 36 37
#if defined(SPHERICAL) || defined(CURVILINEAR)
   use domain, only: H,az,kmax,arcd1
#else
   use domain, only: H,az,kmax,ard1
#endif
   use m2d, only: dtm
   use variables_2d, only: z
38
#ifndef NO_BAROCLINIC
gotm's avatar
gotm committed
39 40
   use m3d, only: calc_salt,calc_temp
   use variables_3d, only: hn,ssen,T,S
kbk's avatar
kbk committed
41 42 43 44 45
#endif
#ifdef GETM_BIO
   use bio, only: bio_calc
   use bio_var, only: numc
   use variables_3d, only: cc3d
46
#endif
gotm's avatar
gotm committed
47 48 49 50 51 52
   IMPLICIT NONE
!
   private
!
! !PUBLIC DATA MEMBERS:
   public init_rivers, do_rivers, clean_rivers
kbk's avatar
kbk committed
53 54 55
#ifdef GETM_BIO
   public init_rivers_bio
#endif
56 57 58
   integer, public                     :: river_method=0,nriver=0,rriver=0
   logical,public                      :: use_river_temp = .false.
   logical,public                      :: use_river_salt = .false.
kbk's avatar
kbk committed
59 60
   character(len=64), public           :: river_data="rivers.nc"
   character(len=64), public, allocatable  :: river_name(:)
61
   character(len=64), public, allocatable  :: real_river_name(:)
kbk's avatar
kbk committed
62
   integer, public, allocatable        :: ok(:)
63 64 65
   REALTYPE, public, allocatable       :: river_flow(:)
   REALTYPE, public, allocatable       :: river_salt(:)
   REALTYPE, public, allocatable       :: river_temp(:)
kbk's avatar
kbk committed
66
   REALTYPE, public                    :: river_factor= _ONE_
67 68 69
   REALTYPE, public,parameter          :: temp_missing=-9999.0
   REALTYPE, public,parameter          :: salt_missing=-9999.0
   integer,  public, allocatable       :: river_split(:)
kbk's avatar
kbk committed
70 71 72 73
#ifdef GETM_BIO
   REALTYPE, public, allocatable       :: river_bio(:,:)
   REALTYPE, public, parameter         :: bio_missing=-9999.0
#endif
gotm's avatar
gotm committed
74 75
!
! !PRIVATE DATA MEMBERS:
kbk's avatar
kbk committed
76 77 78 79
   integer                   :: river_format=2
   character(len=64)         :: river_info="riverinfo.dat"
   integer, allocatable      :: ir(:),jr(:)
   REALTYPE, allocatable     :: irr(:)
80 81
   REALTYPE, allocatable     :: macro_height(:)
   REALTYPE, allocatable     :: flow_fraction(:)
gotm's avatar
gotm committed
82
!
kbk's avatar
kbk committed
83 84 85
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
gotm's avatar
gotm committed
86 87 88 89 90 91 92 93 94 95 96 97
! !LOCAL VARIABLES:
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_rivers
!
! !INTERFACE:
98
   subroutine init_rivers
gotm's avatar
gotm committed
99 100
!
! !DESCRIPTION:
101 102 103 104 105 106 107 108 109
!
! First of all, the namelist {\tt rivers} is read from getm.F90 and
! a number of vectors with the length of {\tt nriver} (number of
! rivers) is allocated. Then, by looping over all rivers, the 
! ascii file {\tt river\_info} is read, and checked for consistency.
! The number of used rivers {\tt rriver} is calculated and it is checked
! whether they are on land (which gives a warning) or not. When a river name
! occurs more than once in {\tt river\_info}, it means that its runoff
! is split among several grid boxed (for wide river mouths).
gotm's avatar
gotm committed
110 111 112 113 114 115 116 117 118 119 120
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
121
   integer                   :: i,j,n,nn,ni,rc,m
kbk's avatar
kbk committed
122
   integer                   :: unit = 25 ! kbk
123
   logical                   :: outside
124
   REALTYPE                  :: area
kbk's avatar
kbk committed
125
   NAMELIST /rivers/ &
126 127
            river_method,river_info,river_format,river_data,river_factor, &
            use_river_salt,use_river_temp
gotm's avatar
gotm committed
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'init_rivers() # ',Ncall
#endif

   LEVEL1 'init_rivers()'
   read(NAMLST,rivers)

   select case (river_method)
      case (0)
         LEVEL3 'River runoff not included.'
      case (1,2)
         LEVEL2 'river_method= ',river_method
         LEVEL2 'river_data=   ',trim(river_data)
         LEVEL2 'river_format= ',river_format
147 148
         LEVEL2 'use_river_temp= ',use_river_temp
         LEVEL2 'use_river_salt= ',use_river_salt
gotm's avatar
gotm committed
149
         open(unit,file=river_info,action='read',status='old',err=90)
kbk's avatar
kbk committed
150
         read(unit,*) nriver
gotm's avatar
gotm committed
151
         allocate(ir(nriver),stat=rc) ! i index of rivers
152
         if (rc /= 0) stop 'rivers: Error allocating memory (ir)'
gotm's avatar
gotm committed
153
         allocate(jr(nriver),stat=rc) ! j index of rivers
154
         if (rc /= 0) stop 'rivers: Error allocating memory (jr)'
kbk's avatar
kbk committed
155
         allocate(ok(nriver),stat=rc) ! valid river spec.
156 157 158
         if (rc /= 0) stop 'rivers: Error allocating memory (ok)'
         allocate(river_name(nriver),stat=rc) ! NetCDF name of river.
         if (rc /= 0) stop 'rivers: Error allocating memory (river_name)'
159
         allocate(river_flow(nriver),stat=rc) ! river flux
160
         if (rc /= 0) stop 'rivers: Error allocating memory (river_flow)'
161
         allocate(macro_height(nriver),stat=rc) ! height over a macro tims-step
162 163 164 165 166 167 168 169 170
         if (rc /= 0) stop 'rivers: Error allocating memory (macro_height)'
         allocate(river_temp(nriver),stat=rc) ! temperature of river water
         if (rc /= 0) stop 'rivers: Error allocating memory (river_temp)'
         allocate(river_salt(nriver),stat=rc) ! salinity of river water
         if (rc /= 0) stop 'rivers: Error allocating memory (river_salt)'
         allocate(river_split(nriver),stat=rc) ! split factor for river water
         if (rc /= 0) stop 'rivers: Error allocating memory (river_split)'
         allocate(flow_fraction(nriver),stat=rc) ! areafactor of data for river 
         if (rc /= 0) stop 'rivers: Error allocating memory (flow_fraction)'
gotm's avatar
gotm committed
171
         allocate(irr(nriver),stat=rc) ! integrated river runoff
172 173
         if (rc /= 0) stop 'rivers: Error allocating memory (irr)'

kbk's avatar
kbk committed
174
         ok = 0
175 176
         rriver = 0 ! number of real existing rivers...
         flow_fraction = _ZERO_
gotm's avatar
gotm committed
177 178
         do n=1,nriver
            read(unit,*) ir(n),jr(n),river_name(n)
179
            river_name(n) = trim(river_name(n))
gotm's avatar
gotm committed
180
            LEVEL3 trim(river_name(n)),':',ir(n),jr(n)
181 182
            i = ir(n)-ioff
            j = jr(n)-joff
183 184 185
            river_temp(n) = temp_missing
            river_salt(n) = salt_missing
            river_flow(n) = _ZERO_
gotm's avatar
gotm committed
186
            irr(n) = _ZERO_
187
            macro_height(n) = _ZERO_
188 189 190 191 192 193 194 195
!           calculate the number of used rivers, they must be 
!           in sequence !
            rriver = rriver +1
            if ( n .gt. 1 ) then
               if (river_name(n) .eq. river_name(n-1)) then 
                  rriver = rriver-1
               end if
            end if
196 197 198 199 200 201 202 203 204 205
            outside= &
                    i .lt. iimin .or. i .gt. iimax .or.  &
                    j .lt. jjmin .or. j .gt. jjmax
            if( .not. outside) then
               if(az(i,j) .eq. 0) then
                  LEVEL3 'Warning:  river# ',n,' at (',i,j,') is on land'
                  LEVEL3 '          setting ok to 0'
                  ok(n) = 0
               else
                  ok(n) = 1
206
                  flow_fraction(n) = _ONE_/ARCD1
207
               end if
gotm's avatar
gotm committed
208
            else
kbk's avatar
kbk committed
209
              LEVEL3 'Outside: river# ',n
gotm's avatar
gotm committed
210 211
            end if
         end do
212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264

!        calculate the number of used gridboxes, they must be 
!        in sequence !
         LEVEL3 'Number of unique rivers: ',rriver
         allocate(real_river_name(rriver),stat=rc) ! NetCDF name of river.
         if (rc /= 0) stop 'rivers: Error allocating memory (rivers)'
         river_split = 1    ! normal case
         do n=2,nriver
               if (river_name(n) .eq. river_name(n-1)) then 
                  river_split(n)=river_split(n-1)+1
            end if
         end do
         ni= nriver
         do n=1,nriver
            if (ni .ge. 1) then
               if ( river_split(ni) .gt. 1 ) then  
                  do m=1,river_split(ni)
                     river_split(ni-m+1) =  river_split(ni)
                  end do
               end if
               ni = ni - river_split(ni)
            end if
         end do
         LEVEL3 'split:',river_split
!        now river_split contains the number of gridboxes used 
!        for a single river
         nn = 1
         ni = 1
         do n=1,nriver
            if (ni .le. nriver) then
               real_river_name(nn) = river_name(ni)
               if ( river_split(ni) .gt. 1 ) then
                  area = _ZERO_
                  do m=1,river_split(ni) 
                     area = area +  flow_fraction(ni+m-1)
                  end do
                  do m=1,river_split(ni)
                     if ( area .gt. _ZERO_ ) then
                        flow_fraction(ni+m-1) = flow_fraction(ni+m-1)/area
                     else
                        flow_fraction(ni+m-1) = _ZERO_
                     end if
                  end do
               else
                  flow_fraction(ni) = _ONE_
               end if
               nn = nn + 1  
               ni = ni + river_split(ni)
            end if 
            if (ok(n) .eq. 0) then
               flow_fraction(n) = _ZERO_
            end if
         end do
gotm's avatar
gotm committed
265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281
      case default
         FATAL 'A non valid river_method has been selected'
         stop 'init_rivers'
   end select
   return

90 LEVEL2 'could not open ',trim(river_info),' for reading info on rivers'
   stop 'init_rivers()'

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

kbk's avatar
kbk committed
282 283 284 285 286 287 288 289 290 291
#ifdef GETM_BIO
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_rivers_bio
!
! !INTERFACE:
   subroutine init_rivers_bio()
!
! !DESCRIPTION:
292 293 294
! First, memory for storing the biological loads from rivers is 
! allocated.
! The variable - {\tt river\_bio} - is initialised to  - {\tt bio\_missing}.
kbk's avatar
kbk committed
295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
   integer                   :: rc
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'init_rivers_bio() # ',Ncall
#endif

   LEVEL1 'init_rivers_bio()'

   allocate(river_bio(nriver,numc),stat=rc)
   if (rc /= 0) stop 'rivers: Error allocating memory (river_bio)'

   river_bio = bio_missing


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

gotm's avatar
gotm committed
333 334 335
!-----------------------------------------------------------------------
!BOP
!
336
! !IROUTINE:  do_rivers - updating river points \label{sec-do-rivers}
gotm's avatar
gotm committed
337 338 339 340 341
!
! !INTERFACE:
   subroutine do_rivers(do_3d)
!
! !DESCRIPTION:
342 343 344 345 346 347 348
! 
! Here, the temperature, salinity, sea surface elevation and layer heights
! are updated in the river inflow grid boxes. Temperature and salinity
! are mixed with riverine values proportional to the old volume and the
! river inflow volume at that time step, sea surface elevation is simply
! increased by the inflow volume divided by the grid box area, and
! the layer heights are increased proportionally. 
gotm's avatar
gotm committed
349 350 351 352 353
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
354
   logical, intent(in)                 :: do_3d
gotm's avatar
gotm committed
355 356 357 358 359 360
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
361
   integer                   :: i,j,k,m,n
kbk's avatar
kbk committed
362
   REALTYPE                  :: rvol,height
363
   REALTYPE                  :: svol,tvol,vol
gotm's avatar
gotm committed
364 365 366 367 368 369 370 371 372 373 374 375 376 377 378
!
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'do_rivers() # ',Ncall
#endif

   select case (river_method)
      case(0)
      case(1,2)
         do n=1,nriver
            if(ok(n) .gt. 0) then
379
               i = ir(n)-ioff; j = jr(n)-joff
380
               rvol = dtm * river_flow(n) * flow_fraction(n)
gotm's avatar
gotm committed
381
               irr(n) = irr(n) + rvol
382
               height = rvol * ARCD1
gotm's avatar
gotm committed
383
               z(i,j) = z(i,j) + height
384
#ifndef NO_BAROCLINIC
385 386
               macro_height(n)=macro_height(n)+height
!              on macrotime step adjust 3d fields
gotm's avatar
gotm committed
387 388
               if (do_3d) then
                  if (calc_salt) then
389 390 391 392 393 394 395 396
                     if ( river_salt(n) .ne. salt_missing ) then
                        S(i,j,1:kmax) = (S(i,j,1:kmax)*(H(i,j)+ssen(i,j))   &
                                      + river_salt(n)*macro_height(n))      &
                                      / (H(i,j)+ssen(i,j)+macro_height(n))
                     else
                        S(i,j,1:kmax) = S(i,j,1:kmax)*(H(i,j)+ssen(i,j))   &
                                      / (H(i,j)+ssen(i,j)+macro_height(n))
                     end if
gotm's avatar
gotm committed
397
                  end if
398 399 400 401 402
                  if (calc_temp .and. river_temp(n) .ne. temp_missing) then
                     T(i,j,1:kmax) = (T(i,j,1:kmax)*(H(i,j)+ssen(i,j))   &
                                      + river_temp(n)*macro_height(n))      &
                                      / (H(i,j)+ssen(i,j)+macro_height(n))
                  end if
kbk's avatar
kbk committed
403 404 405 406 407 408 409 410 411 412 413 414
#ifdef GETM_BIO
                  if (bio_calc) then
                     do m=1,numc
                        if ( river_bio(n,m) .ne. bio_missing ) then
                           cc3d(m,i,j,1:kmax) = &
                                 (cc3d(m,i,j,1:kmax)*(H(i,j)+ssen(i,j)) &
                                 + river_bio(n,m)*macro_height(n))      &
                                 / (H(i,j)+ssen(i,j)+macro_height(n))
                        end if
                     end do
                  end if
#endif
415
!                 Changes of total and layer height due to river inflow:
gotm's avatar
gotm committed
416
                  hn(i,j,1:kmax) = hn(i,j,1:kmax)/(H(i,j)+ssen(i,j)) &
417 418 419
                                  *(H(i,j)+ssen(i,j)+macro_height(n))
                  ssen(i,j) = ssen(i,j)+macro_height(n)
                  macro_height(n) = _ZERO_
gotm's avatar
gotm committed
420
               end if
421
#endif
gotm's avatar
gotm committed
422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439
            end if
         end do
      case default
         FATAL 'Not valid rivers_method specified'
         stop 'init_rivers'
   end select

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

!-----------------------------------------------------------------------
!BOP
!
440
! !IROUTINE:  clean_rivers 
gotm's avatar
gotm committed
441 442
!
! !INTERFACE:
443
   subroutine clean_rivers
gotm's avatar
gotm committed
444 445
!
! !DESCRIPTION:
446 447 448 449
!  
! This routine closes the river handling by writing the integrated
! river run-off for each river to standard output.
! 
gotm's avatar
gotm committed
450 451 452 453 454 455 456 457 458 459 460
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
461 462
   integer                   :: i,j,n
   REALTYPE                  :: tot=_ZERO_
gotm's avatar
gotm committed
463 464 465 466 467 468 469 470 471 472 473 474 475 476 477
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'clean_rivers() # ',Ncall
#endif

   select case (river_method)
      case(0)
      case(1,2)
         do n=1,nriver
            if(ok(n) .gt. 0) then
               i = ir(n); j = jr(n)
478
               LEVEL2 trim(river_name(n)),':  ' ,irr(n)/1.e6, '10^6 m3'
kbk's avatar
kbk committed
479
               tot = tot+irr(n)
gotm's avatar
gotm committed
480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501
            end if
         end do
      case default
         FATAL 'Not valid rivers_method specified'
         stop 'init_rivers'
   end select

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

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

   end module rivers

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