fv_zones.F90

Go to the documentation of this file.

!###############################################################################
!#                                                                             #
!# fv_zones.F90                                                                #
!#                                                                             #
!# Interface for FV (Finite Volume) Model to AED2 modules.                     #
!#   Designed for TUFLOW-FV, released by BMT-WBM:                              #
!#   http://www.tuflow.com/Tuflow%20FV.aspx                                    #
!#                                                                             #
!# This is a support module to allow ability for benthic/sediment zones in     #
!# AED2, including zone-averaging                                              #
!#                                                                             #
!#   -----------------------------------------------------------------------   #
!#                                                                             #
!# Developed by :                                                              #
!#     AquaticEcoDynamics (AED) Group                                          #
!#     School of Earth & Environment                                           #
!# (C) The University of Western Australia                                     #
!#                                                                             #
!# Copyright by the AED-team @ UWA under the GNU Public License - www.gnu.org  #
!#                                                                             #
!#   -----------------------------------------------------------------------   #
!#                                                                             #
!# Created Apr 2015                                                            #
!#                                                                             #
!###############################################################################

#if defined HAVE_AED2
#include "aed2.h"
!#include "aed2_common.h"
#endif

#ifndef DEBUG
#define DEBUG      0
#endif

!###############################################################################
MODULE fv_zones
!-------------------------------------------------------------------------------
#if defined HAVE_AED2
   USE aed2_common
   USE declarations_special

   IMPLICIT NONE

   PUBLIC init_zones, calc_zone_areas, copy_to_zone, copy_from_zone
   PUBLIC compute_zone_benthic_fluxes !, distribute_flux_from_zone
   PUBLIC zone

   !#--------------------------------------------------------------------------#
   !# Module Data

   !# Arrays for environmental variables not supplied externally.
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone
   aed_real,DIMENSION(:,:),ALLOCATABLE,TARGET :: zone_cc
   aed_real,DIMENSION(:,:),ALLOCATABLE,TARGET :: zone_cc_diag
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_area
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_temp
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_salt
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_rho
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_height
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_extc
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_tss
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_par
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_wind
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_rain ! JC not sure it makes sense
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_i_0
   aed_real,DIMENSION(:),  ALLOCATABLE,TARGET :: zone_taub
   INTEGER, DIMENSION(:),  ALLOCATABLE        :: zone_count, zm

   aed_real,DIMENSION(:,:),ALLOCATABLE,TARGET :: flux_pelz
   aed_real,DIMENSION(:,:),ALLOCATABLE,TARGET :: flux_benz

   INTEGER :: nzones, nwq_var, nben_var

!#####################################################

CONTAINS
!===============================================================================


!###############################################################################
SUBROUTINE init_zones(nCols, mat_id, n_aed2_vars, n_vars, n_vars_ben)
!-------------------------------------------------------------------------------
!ARGUMENTS
   INTEGER,INTENT(in) :: nCols
   INTEGER,DIMENSION(:,:),INTENT(in) :: mat_id
   INTEGER,INTENT(in) :: n_aed2_vars, n_vars, n_vars_ben
!
!LOCALS
   INTEGER :: cType, nTypes
   INTEGER :: col, zon
   INTEGER,DIMENSION(:),ALLOCATABLE :: mat_t
!
!-------------------------------------------------------------------------------
!BEGIN
   ALLOCATE(mat_t(ncols)) ; ALLOCATE(zm(ncols))
   !# The new form of zones
   ctype = mat_id(1,1) ; ntypes = 1 ; mat_t(ntypes) = mat_id(1,1)
   DO col=1, ubound(mat_id,2)
      !# use the bottom index to fill the array
      IF ( ctype /= mat_id(1,col) ) THEN
         DO zon=1,ntypes
            IF ( mat_t(zon) .eq. mat_id(1,col) ) THEN
               ctype = mat_id(1,col)
               EXIT
            ENDIF
         ENDDO
      ENDIF
      IF ( ctype /= mat_id(1,col) ) THEN
         ntypes = ntypes + 1
         mat_t(ntypes) = mat_id(1,col)
         ctype = mat_id(1,col)
         zon = ntypes
      ENDIF
      zm(col) = zon
   ENDDO
   WRITE(lu, *) "Number of mats = ", ntypes, " = ", mat_t(1:ntypes)
   DEALLOCATE(mat_t)
   nzones = ntypes

   ALLOCATE(zone(nzones))
   ALLOCATE(zone_area(nzones))
   ALLOCATE(zone_temp(nzones))
   ALLOCATE(zone_salt(nzones))
   ALLOCATE(zone_rho(nzones))
   ALLOCATE(zone_height(nzones))

   ALLOCATE(zone_extc(nzones))
   ALLOCATE(zone_tss(nzones))
   ALLOCATE(zone_par(nzones))
   ALLOCATE(zone_wind(nzones))
   ALLOCATE(zone_rain(nzones)) ! JC
   ALLOCATE(zone_i_0(nzones))
   ALLOCATE(zone_taub(nzones))

   ALLOCATE(zone_count(nzones))

   ALLOCATE(zone_cc(n_vars, nzones))
   ALLOCATE(zone_cc_diag(n_aed2_vars-n_vars, nzones))

   ALLOCATE(flux_pelz(n_vars, nzones))
   ALLOCATE(flux_benz(n_vars, nzones))

   nwq_var = n_vars
   nben_var = n_vars_ben
END SUBROUTINE init_zones
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!###############################################################################
SUBROUTINE calc_zone_areas(nCols, temp, salt, h, area, wnd, rho, extcoeff, I_0, par, tss, active, rain)!, bathy ,rain JC
!-------------------------------------------------------------------------------
!ARGUMENTS
   INTEGER,INTENT(in) :: nCols
   aed_real,DIMENSION(:),INTENT(in) :: temp, salt, h, area, wnd, rain
   aed_real,DIMENSION(:),INTENT(in) :: rho, extcoeff, i_0, par, tss
   LOGICAL,DIMENSION(:), INTENT(in) :: active
!
!LOCALS
   INTEGER :: col, zon
!
!-------------------------------------------------------------------------------
!BEGIN
   zone_area = zero_
   zone_temp = zero_
   zone_salt = zero_
   zone_height = zero_
   zone_wind = zero_
   zone_rain = zero_ ! JC
   zone_count = 0

   DO col=1, ncols
      IF (.NOT. active(col)) cycle

      zon = zm(col)
      zone_area(zon)   = zone_area(zon) + area(col)
      zone_temp(zon)   = zone_temp(zon) + temp(col)
      zone_salt(zon)   = zone_salt(zon) + salt(col)
      zone_rho(zon)    = zone_rho(zon) + rho(col)
      zone_height(zon) = zone_height(zon) + h(col)
      zone_extc(zon)   = zone_extc(zon) + extcoeff(col)
      zone_tss(zon)    = zone_tss(zon) + tss(col)
      zone_par(zon)    = zone_par(zon) + par(col)
      zone_wind(zon)   = zone_wind(zon) + wnd(col)
      zone_rain(zon)   = zone_rain(zon) + rain(col) !JC
      zone_i_0(zon)    = zone_i_0(zon) + i_0(col)
    ! zone_taub(zon)   = zone_taub(zon) + col_taub

      zone_count(zon)  = zone_count(zon) + 1
   ENDDO

   zone_temp = zone_temp / zone_count
   zone_salt = zone_salt / zone_count

   zone_rho  =  zone_rho / zone_count
   zone_extc = zone_extc / zone_count
   zone_tss  =  zone_tss / zone_count
   zone_par  =  zone_par / zone_count
   zone_wind = zone_wind / zone_count
   zone_rain = zone_rain / zone_count ! JC
   zone_i_0  =  zone_i_0 / zone_count
   zone_taub = zone_taub / zone_count

END SUBROUTINE calc_zone_areas
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!###############################################################################
SUBROUTINE copy_to_zone(nCols, cc, area, active, benth_map)
!-------------------------------------------------------------------------------
!ARGUMENTS
   INTEGER,INTENT(in)  :: nCols
   aed_real,INTENT(in) :: cc(:,:)       !# (n_vars, n_layers)
   aed_real,DIMENSION(:),INTENT(in) :: area
   LOGICAL,DIMENSION(:), INTENT(in) :: active
   INTEGER,DIMENSION(:), INTENT(in) :: benth_map
!
!LOCALS
   INTEGER :: col, zon, bot
!
!-------------------------------------------------------------------------------
!BEGIN
   DO col=1, ncols
      IF (.NOT. active(col)) cycle

      bot = benth_map(col)
      zon = zm(col)

      zone_cc(1:nwq_var,zon) = zone_cc(1:nwq_var,zon) + (cc(1:nwq_var,bot) * area(col) / zone_area(zon))
   ENDDO
END SUBROUTINE copy_to_zone
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!###############################################################################
SUBROUTINE copy_from_zone(nCols, cc, area, active, benth_map)
!-------------------------------------------------------------------------------
!ARGUMENTS
   INTEGER,INTENT(in)     :: nCols
   aed_real,INTENT(inout) :: cc(:,:)       !# (n_vars, n_layers)
   aed_real,DIMENSION(:),INTENT(in) :: area
   LOGICAL,DIMENSION(:), INTENT(in) :: active
   INTEGER,DIMENSION(:), INTENT(in) :: benth_map
!
!LOCALS
   INTEGER :: col, zon, bot
!
!-------------------------------------------------------------------------------
!BEGIN
   DO col=1, ncols
      IF (.NOT. active(col)) cycle

      bot = benth_map(col)
      zon = zm(col)

      cc(nwq_var+1:nwq_var+nben_var,bot) = zone_cc(nwq_var+1:nwq_var+nben_var,zon)
   ENDDO
END SUBROUTINE copy_from_zone
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!###############################################################################
SUBROUTINE stopit(message)
!-------------------------------------------------------------------------------
!ARGUMENTS
   CHARACTER(*) :: message
!-------------------------------------------------------------------------------
   print *,message
   stop
END SUBROUTINE stopit
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!###############################################################################
SUBROUTINE define_column_zone(column, zon, n_aed2_vars)
!-------------------------------------------------------------------------------
!ARGUMENTS
   type(aed2_column_t), INTENT(inout) :: column(:)
   INTEGER, INTENT(in) :: zon, n_aed2_vars
!
!LOCALS
   INTEGER :: av
   INTEGER :: v, d, sv, sd, ev
   TYPE(aed2_variable_t),POINTER :: tvar
!
!-------------------------------------------------------------------------------
!BEGIN
   v = 0 ; d = 0; sv = 0; sd = 0 ; ev = 0
   DO av=1,n_aed2_vars

      IF ( .NOT. aed2_get_var(av, tvar) ) stop "Error getting variable info"

      IF ( tvar%extern ) THEN !# global variable
         ev = ev + 1
         SELECT CASE (tvar%name)
            CASE ( 'temperature' ) ; column(av)%cell => zone_temp
            CASE ( 'salinity' )    ; column(av)%cell => zone_salt
            CASE ( 'density' )     ; column(av)%cell => zone_rho
            CASE ( 'layer_ht' )    ; column(av)%cell => zone_height
            CASE ( 'layer_area' )  ; column(av)%cell_sheet => zone_area(zon)
            CASE ( 'rain' )        ; column(av)%cell_sheet => zone_rain(zon) ! JC
       !    CASE ( 'rainloss' )    ; column(av)%cell_sheet => zone_rainloss(zon) ! JC
            CASE ( 'material' )    ; column(av)%cell_sheet => zone(zon)
       !    CASE ( 'bathy' )       ; column(av)%cell_sheet => zone_bathy(zon)! Does it mean this var is already pointed to the correct variable of TFFV? JC
            CASE ( 'extc_coef' )   ; column(av)%cell => zone_extc
            CASE ( 'tss' )         ; column(av)%cell => zone_tss
       !    CASE ( 'par' )         ; column(av)%cell => zone_par
       !    CASE ( 'par' )         ; IF (link_ext_par) THEN
       !                                column(av)%cell => lpar(top:bot)
       !                             ELSE
       !                                column(av)%cell => par(top:bot)
       !                             ENDIF
       !    CASE ( 'nir' )         ; column(av)%cell => zone_nir
       !    CASE ( 'uva' )         ; column(av)%cell => zone_uva
       !    CASE ( 'uvb' )         ; column(av)%cell => zone_uvb
            CASE ( 'sed_zone' )    ; column(av)%cell_sheet => zone(zon)
            CASE ( 'wind_speed' )  ; column(av)%cell_sheet => zone_wind(zon)
            CASE ( 'par_sf' )      ; column(av)%cell_sheet => zone_i_0(zon)
       !    CASE ( 'taub' )        ; column(av)%cell_sheet => zone_taub
       !    CASE ( 'air_temp' )    ; column(av)%cell_sheet => zone_air_temp(zon)
            CASE DEFAULT ; CALL stopit("ERROR: external variable "//trim(tvar%name)//" not found.")
         END SELECT
      ELSEIF ( tvar%diag ) THEN  !# Diagnostic variable
         d = d + 1
         IF ( tvar%sheet ) THEN
            column(av)%cell_sheet => zone_cc_diag(d,zon)
         ELSE
            column(av)%cell => zone_cc_diag(d,zon:zon)
         ENDIF
      ELSE    !# state variable
         IF ( tvar%sheet ) THEN
            sv = sv + 1
            IF ( tvar%bot ) THEN
               column(av)%cell_sheet => zone_cc(nwq_var+sv, zon)
            ELSEIF ( tvar%top ) THEN
    !          column(av)%cell_sheet => zone_cc(nwq_var+sv, top)
            ENDIF
            column(av)%flux_ben => flux_benz(nwq_var+sv, zon)
    !       column(av)%flux_atm => flux_atm(nwq_var+sv)
         ELSE
            v = v + 1
            column(av)%cell => zone_cc(v, 1:1)
            column(av)%flux_pel => flux_pelz(v, 1:1)
            column(av)%flux_ben => flux_benz(v, zon)
    !       column(av)%flux_atm => flux_atm(v)
         ENDIF
      ENDIF
   ENDDO
END SUBROUTINE define_column_zone
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!###############################################################################
SUBROUTINE compute_zone_benthic_fluxes(n_aed2_vars, dt)
!-------------------------------------------------------------------------------
!ARGUMENTS
   INTEGER,INTENT(in) :: n_aed2_vars
   aed_real,INTENT(in):: dt
!
!LOCALS
   INTEGER :: zon, v
   type(aed2_column_t) :: column(n_aed2_vars)
!
!-------------------------------------------------------------------------------
!BEGIN
   DO zon=1, nzones
      CALL define_column_zone(column, zon, n_aed2_vars)!, nwq_var)

      CALL aed2_calculate_benthic(column, 1)
      DO v=nwq_var+1,nwq_var+nben_var
         zone_cc(v, 1) = zone_cc(v, 1) + dt*flux_benz(v, zon);
      ENDDO
   ENDDO
END SUBROUTINE compute_zone_benthic_fluxes
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


#endif
!===============================================================================
END MODULE fv_zones