Subroutine waterbalance

  1 : !***********************************************************************
  2 : !*  Layered water balance
  3 : !*----------------------------------------------------------------------
  4 : !*  Author: Stan Schymanski, Max Planck Institute for Biogeochemistry
  5 : !*  Email: sschym@bgc-jena.mpg.de
  6 : !*  06/2008
  7 : !*  Version: REW with layered unsaturated zone, no routing
  8 : !*----------------------------------------------------------------------
  9 : !*
 10 : !* Numbers in the commented parentheses refer to the equation numeration
 11 : !* in the document 'equations.pdf' that comes with the documentation.
 12 : !*
 13 : !*----------------------------------------------------------------------
 14 : !*
 15 : !* This subroutine ('waterbalance') uses the variables defined in
 16 : !* modules.f90, some of which have to be allocated prior to calling
 17 : !* 'waterbalance'.
 18 : !* When calling the 'waterbalance', the 'init' variable must be given
 19 : !* (1 to generate initial conditions, otherwise 0).
 20 : !* This subroutine calculates the water balance for a time step <=dtmax.
 21 : !* The calling program must transfer the new variables ('new' in their names)
 22 : !* to the old variables for the next time step, before calling 'waterbalance'
 23 : !* again. Example: call waterbalance(init) , time=time+dt, zw=zwnew...
 24 : !*
 25 : !*----------------------------------------------------------------------
 26 : !*  Copyright (C) 2008  Stan Schymanski
 27 : !*
 28 : !*    This program is free software: you can redistribute it and/or modify
 29 : !*    it under the terms of the GNU General Public License as published by
 30 : !*    the Free Software Foundation, either version 3 of the License, or
 31 : !*    (at your option) any later version.
 32 : !*
 33 : !*    This program is distributed in the hope that it will be useful,
 34 : !*    but WITHOUT ANY WARRANTY; without even the implied warranty of
 35 : !*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 36 : !*    GNU General Public License for more details.
 37 : !*
 38 : !*    You should have received a copy of the GNU General Public License
 39 : !*    along with this program. If not, see http://www.gnu.org/licenses.
 40 : !*
 41 : !***********************************************************************
 42 : 
 43 :       subroutine waterbalance()
 44 :       use vom_vegwat_mod
 45 :       implicit none
 46 : 
 47 :       INTEGER :: jj
 48 : 
 49 : !     * CURRENT SOIL WATER CONTENT
 50 : 
 51 :       cH2Ol_s(:) = (-su__(:) * s_thetar(:) + su__(:) * s_thetas(:)     &
 52 :      &           + s_thetar(:)) * s_delz(:)
 53 :       wc_ = SUM(cH2Ol_s)
 54 : 
 55 : !     * FLUXES (inf, infx, qbl, esoil__, spgfcf__)
 56 : 
 57 :       call waterbalance_fluxes()
 58 : 
 59 : !     * changes in water storage (waterbalance)
 60 : 
 61 :       if( i_no_veg .eq. 0) then
 62 :          io__     = inf__ - esoil__ - spgfcf__ - SUM(ruptkt__(:)) - SUM(ruptkg__(:))  ! (3.19)
 63 :       else
 64 :          io__     = inf__ - esoil__ - spgfcf__  ! (no vegetation, WB still needs to close)
 65 :       end if
 66 : 
 67 :       iovec(:) = 0.d0
 68 :       iovec(1) = qbl(1) + inf__ - esoil__ - ruptkt__(1) - ruptkg__(1)
 69 :       if (wlayer_ .eq. 1) then
 70 :         iovec(1) = iovec(1) - spgfcf__
 71 :       endif
 72 :       if (wlayer_ .gt. 2) then
 73 :         do jj = 2, wlayer_ - 1
 74 :           iovec(jj) = qbl(jj) - qbl(jj-1) - ruptkt__(jj) - ruptkg__(jj)
 75 :         enddo
 76 :       endif
 77 :       if (wlayer_ .gt. 1) then
 78 :         iovec(wlayer_) = qbl(wlayer_) - qbl(wlayer_-1)                 &
 79 :      &                 - ruptkt__(wlayer_) - ruptkg__(wlayer_) - spgfcf__
 80 :       endif
 81 : 
 82 : !     * change in saturation degree
 83 : 
 84 :       dsu(:) = -iovec(:) / ((s_thetar(:) - s_thetas(:)) * s_delz(:))
 85 : 
 86 : !     * calculation of maximal time step size
 87 : 
 88 :       call waterbalance_timestep()
 89 : 
 90 : !     Calculating state variables at next time step
 91 : !     (sunew, wlayernew, pcapnew, kunsatnew, zwnew)
 92 : 
 93 :       call waterbalance_update_state()
 94 : 
 95 :       call waterbalance_diag()
 96 : 
 97 :       return
 98 :       end subroutine waterbalance