stock_constants.F90

!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
!* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
!* for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************
!> @defgroup stock_constants_mod stock_constants_mod
!> @ingroup exchange
!> @brief Parameters, routines, and types for computing stocks in @ref xgrid_mod
 
!> @addtogroup stock_constants_mod
!> @{
module stock_constants_mod
 
  use mpp_mod, only : mpp_pe, mpp_root_pe, mpp_sum
  use fms_mod, only : write_version_number
  use time_manager_mod, only : time_type, get_time
  use time_manager_mod, only : operator(+), operator(-)
  use diag_manager_mod, only : register_diag_field,send_data
  use platform_mod,     only : r8_kind
 
  implicit none
 
  ! Include variable "version" to be written to log file.
#include<file_version.h>
 
  integer,public,    parameter                :: NELEMS=3
  integer,           parameter                :: NELEMS_report=3
  integer,public,    parameter                :: ISTOCK_WATER=1, istock_heat=2, istock_salt=3
  integer,public,    parameter                :: ISTOCK_TOP=1, istock_bottom=2, istock_side=3
  integer,public                              :: stocks_file
  ! Stock related stuff
  ! Shallow (no constructor) data structures holding the starting stock values (per PE) and
  ! flux integrated increments at present time.
 
  integer, parameter :: nsides  = 3         !< top, bottom, side
  !> @}
 
  !> @brief Holds stocks amounts per PE values
  !> @ingroup stock_constants_mod
  type stock_type
     real(r8_kind)  :: q_start = 0.0_r8_kind    !< total stocks at start time
     real(r8_kind)  :: q_now   = 0.0_r8_kind    !< total stocks at time t
 
     ! The dq's below are the stocks increments at the present time
     ! delta_t * surf integr of flux
     ! one for each side (ISTOCK_TOP, ISTOCK_BOTTOM, ISTOCK_SIDE)
     real(r8_kind)  :: dq(nsides)    = 0.0_r8_kind    !< stock increments at present time on the Ice   grid
     real(r8_kind)  :: dq_in(nsides) = 0.0_r8_kind    !< stock increments at present time on the Ocean grid
  end type stock_type
  !> @addtogroup stock_constants_mod
  !> @{
 
  type(stock_type), save, public, dimension(NELEMS) :: atm_stock, ocn_stock, lnd_stock, ice_stock
  type(time_type), save :: init_time
 
  public stocks_report
  public stocks_report_init
  public stocks_set_init_time
 
  integer,private,  parameter                     :: ncomps=4
  integer,private,  parameter                     :: istock_atm=1,istock_lnd=2,istock_ice=3,istock_ocn=4
  character(len=3), parameter, dimension(NCOMPS)  :: comp_names=(/'ATM', 'LND', 'ICE', 'OCN'/)
 
 
  character(len=5) , parameter, dimension(NELEMS)  :: stock_names=(/'water', 'heat ', 'salt '/)
  character(len=12), parameter, dimension(NELEMS)  :: stock_units=(/'[Kg]    ','[Joules]','[Kg]    '/)
 
 
contains
 
    !> Starts a stock report
    subroutine stocks_report_init(Time)
    type(time_type), intent(in) :: time !< Model time
 
    character(len=80)                :: formatstring,space
    integer                          :: i,s
    real(r8_kind), dimension(NELEMS) :: val_atm, val_lnd, val_ice, val_ocn
 
    ! Write the version of this file to the log file
    call write_version_number('STOCK_CONSTANTS_MOD', version)
 
    do i = 1, nelems_report
       val_atm(i) = atm_stock(i)%q_start
       val_lnd(i) = lnd_stock(i)%q_start
       val_ice(i) = ice_stock(i)%q_start
       val_ocn(i) = ocn_stock(i)%q_start
       call mpp_sum(val_atm(i))
       call mpp_sum(val_lnd(i))
       call mpp_sum(val_ice(i))
       call mpp_sum(val_ocn(i))
    enddo
 
 
 
    if(mpp_pe() == mpp_root_pe()) then
!       earth_area = 4.*PI*Radius**2
 
       write(stocks_file,*) '================Stocks Report Guide====================================='
       write(stocks_file,*) ' '
       write(stocks_file,*) 'S(t) = Total amount     of a tracer in the component model at time t.'
       write(stocks_file,*) '       Calculated via the component model itself.'
       write(stocks_file,*) ' '
       write(stocks_file,*) 'F(t) = Cumulative input of a tracer to the component model at time t.'
       write(stocks_file,*) '       Calculated via interchange of fluxes with other component models.'
       write(stocks_file,*) ' '
       write(stocks_file,*) 'S(t) - S(0) = Cumulative increase of the component stocks at time t'
       write(stocks_file,*) '              Calculated by the component itself.'
       write(stocks_file,*) ' '
       write(stocks_file,*) 'In a conserving component F(t)=S(t)-S(0) to within numerical accuracy.'
       write(stocks_file,*) ' '
       write(stocks_file,*) 'Component Model refers to one of OCN, ATM, LND or ICE'
       write(stocks_file,*) ''
       write(stocks_file,*) 'NOTE: When use_lag_fluxes=.true. is used in coupler, the ocean stocks '
       write(stocks_file,*) '      calculations are in error by an order which scales as the inverse'
       write(stocks_file,*) '      of the number of time steps.'
       write(stocks_file,*) ' '
       write(stocks_file,*) '======================================================================='
 
 
       write(stocks_file,*) '======================Initial Stock S(0)==============================='
!The following produces  formatString='(5x,a,a,12x,a,a, 9x)' but is general to handle more elements
       formatstring= '(5x'
       do i=1,nelems_report
          s = 25-len_trim(stock_names(i))-len_trim(stock_units(i))
          write(space,'(i2)') s
          formatstring= trim(formatstring)//',a,a,'//trim(space)
          formatstring= trim(formatstring)//trim('x')
       enddo
       formatstring= trim(formatstring)//')'
 
       write(stocks_file,formatstring) (trim(stock_names(i)),trim(stock_units(i)), i=1,nelems_report)
 
!The following produces  formatString=' (a,x,es22.15,3x,es22.15,3x)' but is general to handle more elements
       formatstring= '(a,x'
       do i=1,nelems_report
          write(space,'(i2)') s
          formatstring= trim(formatstring)//',es22.15,3x'
       enddo
       formatstring= trim(formatstring)//')'
 
 
       write(stocks_file,formatstring) 'ATM', (val_atm(i), i=1,nelems_report)
       write(stocks_file,formatstring) 'LND', (val_lnd(i), i=1,nelems_report)
       write(stocks_file,formatstring) 'ICE', (val_ice(i), i=1,nelems_report)
       write(stocks_file,formatstring) 'OCN', (val_ocn(i), i=1,nelems_report)
 
       write(stocks_file,*) '========================================================================'
       write(stocks_file,'(a)'  ) ' '!blank line
 
    end if
 
    call stocks_set_init_time(time)
 
  end subroutine stocks_report_init
 
  !> Writes update to stock report
  subroutine stocks_report(Time)
    type(time_type)               , intent(in) :: time !< Model time
 
    type(stock_type)                 :: stck
    real(r8_kind), dimension(NCOMPS) :: f_value, f_ice_grid, f_ocn_grid, f_ocn_btf, q_start, q_now,c_value
    character(len=80)                :: formatstring
    integer                          :: iday0, isec0, iday, isec, hours
    real(r8_kind)                    :: days
    integer                          :: diagid , comp,elem,i
    integer, parameter :: initid = -2 !< initial value for diag IDs. Must not be equal to the value
                                      !! that register_diag_field returns when it can't register
                                      !! the filed -- otherwise the registration
                                      !! is attempted every time this subroutine is called
 
    integer, dimension(NCOMPS,NELEMS), save :: f_valuediagid = initid
    integer, dimension(NCOMPS,NELEMS), save :: c_valuediagid = initid
    integer, dimension(NCOMPS,NELEMS), save :: fmc_valuediagid = initid
    integer, dimension(NCOMPS,NELEMS), save :: f_lostdiagid = initid
 
    real(r8_kind)     :: diagfield
    logical           :: used
    character(len=30) :: field_name, units
 
    if(mpp_pe()==mpp_root_pe()) then
       call get_time(init_time, isec0, iday0)
       call get_time(time, isec, iday)
 
       hours = iday*24 + isec/3600 - iday0*24 - isec0/3600
       days  = real(hours, r8_kind)/24.0_r8_kind
       write(stocks_file,*) '==============================================='
       write(stocks_file,'(a,f12.3)') 't = TimeSinceStart[days]= ',days
       write(stocks_file,*) '==============================================='
    endif
 
    do elem = 1,nelems_report
 
       do comp = 1,ncomps
 
          if(comp == istock_atm) stck = atm_stock(elem)
          if(comp == istock_lnd) stck = lnd_stock(elem)
          if(comp == istock_ice) stck = ice_stock(elem)
          if(comp == istock_ocn) stck = ocn_stock(elem)
 
 
          f_ice_grid(comp) = sum(stck%dq)
          f_ocn_grid(comp) = sum(stck%dq_IN)
          f_ocn_btf(comp)  = stck%dq_IN( istock_bottom )
 
          q_start(comp) = stck%q_start
          q_now(comp)   = stck%q_now
 
          call mpp_sum(f_ice_grid(comp))
          call mpp_sum(f_ocn_grid(comp))
          call mpp_sum(f_ocn_btf(comp))
          call mpp_sum(q_start(comp))
          call mpp_sum(q_now(comp))
 
          c_value(comp) = q_now(comp) - q_start(comp)
 
          if(mpp_pe() == mpp_root_pe()) then
 
             if(f_valuediagid(comp,elem) == initid) then
                field_name = trim(comp_names(comp)) // trim(stock_names(elem))
                field_name  = trim(field_name) // 'StocksChange_Flux'
                units = trim(stock_units(elem))
                f_valuediagid(comp,elem) = register_diag_field('stock_print', field_name, time, &
                     units=units)
             endif
 
             if(c_valuediagid(comp,elem) == initid) then
                field_name = trim(comp_names(comp)) // trim(stock_names(elem))
                field_name = trim(field_name) // 'StocksChange_Comp'
                units = trim(stock_units(elem))
                c_valuediagid(comp,elem) = register_diag_field('stock_print', field_name, time, &
                     units=units)
             endif
 
             if(fmc_valuediagid(comp,elem) == initid) then
                field_name = trim(comp_names(comp)) // trim(stock_names(elem))
                field_name = trim(field_name) // 'StocksChange_Diff'
                units = trim(stock_units(elem))
                fmc_valuediagid(comp,elem) = register_diag_field('stock_print', field_name, time, &
                     units=units)
             endif
 
             f_value(comp) = f_ice_grid(comp)
 
             if(comp == istock_ocn) then
 
                f_value(comp) = f_ocn_grid(comp)
 
                if(f_lostdiagid(comp,elem) == initid) then
                   field_name = trim(comp_names(comp)) // trim(stock_names(elem))
                   field_name = trim(field_name) // 'StocksExchangeLost'
                   units = trim(stock_units(elem))
                   f_lostdiagid(comp,elem) = register_diag_field('stock_print', field_name, time, &
                        units=units)
                endif
 
                diagid=f_lostdiagid(comp,elem)
                diagfield = f_ice_grid(comp) - f_ocn_grid(comp)
                if (diagid > 0)  used = send_data(diagid, diagfield, time)
 
             endif
 
 
             diagid=f_valuediagid(comp,elem)
             diagfield = f_value(comp)
             if (diagid > 0)  used = send_data(diagid, diagfield, time)
             diagid=c_valuediagid(comp,elem)
             diagfield = c_value(comp)
             if (diagid > 0)  used = send_data(diagid, diagfield, time)
             diagid=fmc_valuediagid(comp,elem)
             diagfield = f_value(comp)-c_value(comp)
             if (diagid > 0)  used = send_data(diagid, diagfield, time)
 
 
             ! formatString = '(a,a,a,i16,2x,es22.15,2x,es22.15,2x,es22.15,2x,es22.15,2x,es22.15,2x,es22.15)'
             !
             !             write(stocks_file,formatString) trim(COMP_NAMES(comp)),STOCK_NAMES(elem),STOCK_UNITS(elem) &
             !                  ,hours, q_now, q_now-q_start, f_value, f_value - (q_now - q_start),
             !                  (f_value - (q_now - q_start))/q_start
 
 
          endif
       enddo
 
 
       if(mpp_pe()==mpp_root_pe()) then
!          write(stocks_file,'(a)'  ) ' '!blank line
!          write(stocks_file,'(a,f12.3)') 't = TimeSinceStart[days]= ',days
!          write(stocks_file,'(a)'  )   ' '!blank line
!          write(stocks_file,'(a,30x,a,20x,a,20x,a,20x,a)') 'Component ','ATM','LND','ICE','OCN'
!          write(stocks_file,'(55x,a,20x,a,20x,a,20x,a)')  'ATM','LND','ICE','OCN'
!          write(stocks_file,'(a,f12.3,12x,a,20x,a,20x,a,20x,a)') 't = TimeSinceStart[days]=
!          ',days,'ATM','LND','ICE','OCN'
 
          write(stocks_file,'(a,a,40x,a,20x,a,20x,a,20x,a)') 'Stocks of ',trim(stock_names(elem)),'ATM','LND', &
               & 'ICE','OCN'
          formatstring = '(a,a,2x,es22.15,2x,es22.15,2x,es22.15,2x,es22.15)'
 
          write(stocks_file,formatstring) 'Total =S(t)               ',stock_units(elem),&
               ( q_now(i), i=1,ncomps)
          write(stocks_file,formatstring) 'Change=S(t)-S(0)          ',stock_units(elem),&
               ( q_now(i)-q_start(i), i=1,ncomps)
          write(stocks_file,formatstring) 'Input =F(t)               ',stock_units(elem),&
               ( f_value(i), i=1,ncomps)
          write(stocks_file,formatstring) 'Diff  =F(t) - (S(t)-S(0)) ',stock_units(elem),&
               ( f_value(i) - c_value(i), i=1,ncomps)
          write(stocks_file,formatstring) 'Error =Diff/S(0)          ','[NonDim]    ', &
               ((f_value(i) - c_value(i))/(1+q_start(i)), i=1,ncomps)  !added 1 to avoid div by zero.
                                                                       !! Assuming q_start large
 
          write(stocks_file,'(a)'  ) ' '!blank line
          formatstring = '(a,a,a,6x,es22.15)'
          write(stocks_file,formatstring) 'Lost Stocks in the exchange between Ice and Ocean ',trim(stock_names(elem))&
                     &,trim(stock_units(elem)), f_ice_grid(istock_ocn) - f_ocn_grid(istock_ocn) + f_ocn_btf(istock_ocn)
 
          write(stocks_file,'(a)') ' ' !blank line
          write(stocks_file,'(a)') ' ' !blank line
 
       endif
    enddo
 
  end subroutine stocks_report
 
  subroutine stocks_set_init_time(Time)
    type(time_type)     , intent(in) :: time !< init time to set for stock report
    init_time = time
 
  end subroutine stocks_set_init_time
 
end module stock_constants_mod
!> @}
! close documentation grouping