exchange/stock_constants.F90

   1 !***********************************************************************
   2 !*                   GNU Lesser General Public License
   3 !*
   4 !* This file is part of the GFDL Flexible Modeling System (FMS).
   5 !*
   6 !* FMS is free software: you can redistribute it and/or modify it under
   7 !* the terms of the GNU Lesser General Public License as published by
   8 !* the Free Software Foundation, either version 3 of the License, or (at
   9 !* your option) any later version.
  10 !*
  11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
  12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 !* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 !* for more details.
  15 !*
  16 !* You should have received a copy of the GNU Lesser General Public
  17 !* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
  18 !***********************************************************************
  19 !> @defgroup stock_constants_mod stock_constants_mod
  20 !> @ingroup exchange
  21 !> @brief Parameters, routines, and types for computing stocks in @ref xgrid_mod
  22
  23 !> @addtogroup stock_constants_mod
  24 !> @{
  25 module stock_constants_mod
  26
  27   use mpp_mod, only : mpp_pe, mpp_root_pe, mpp_sum
  28   use fms_mod, only : write_version_number
  29   use time_manager_mod, only : time_type, get_time
  30   use time_manager_mod, only : operator(+), operator(-)
  31   use diag_manager_mod, only : register_diag_field,send_data
  32   use platform_mod,     only : r8_kind
  33
  34   implicit none
  35
  36   ! Include variable "version" to be written to log file.
  37 #include<file_version.h>
  38
  39   integer,public,    parameter                :: NELEMS=3
  40   integer,           parameter                :: NELEMS_report=3
  41   integer,public,    parameter                :: ISTOCK_WATER=1, ISTOCK_HEAT=2, ISTOCK_SALT=3
  42   integer,public,    parameter                :: ISTOCK_TOP=1, ISTOCK_BOTTOM=2, ISTOCK_SIDE=3
  43   integer,public                              :: stocks_file
  44   ! Stock related stuff
  45   ! Shallow (no constructor) data structures holding the starting stock values (per PE) and
  46   ! flux integrated increments at present time.
  47
  48   integer, parameter :: NSIDES  = 3         !< top, bottom, side
  49   !> @}
  50
  51   !> @brief Holds stocks amounts per PE values
  52   !> @ingroup stock_constants_mod
  53   type stock_type
  54      real(r8_kind)  :: q_start = 0.0_r8_kind    !< total stocks at start time
  55      real(r8_kind)  :: q_now   = 0.0_r8_kind    !< total stocks at time t
  56
  57      ! The dq's below are the stocks increments at the present time
  58      ! delta_t * surf integr of flux
  59      ! one for each side (ISTOCK_TOP, ISTOCK_BOTTOM, ISTOCK_SIDE)
  60      real(r8_kind)  :: dq(NSIDES)    = 0.0_r8_kind    !< stock increments at present time on the Ice   grid
  61      real(r8_kind)  :: dq_IN(NSIDES) = 0.0_r8_kind    !< stock increments at present time on the Ocean grid
  62   end type stock_type
  63   !> @addtogroup stock_constants_mod
  64   !> @{
  65
  66   type(stock_type), save, public, dimension(NELEMS) :: Atm_stock, Ocn_stock, Lnd_stock, Ice_stock
  67   type(time_type), save :: init_time
  68
  69   public stocks_report
  70   public stocks_report_init
  71   public stocks_set_init_time
  72
  73   integer,private,  parameter                     :: NCOMPS=4
  74   integer,private,  parameter                     :: ISTOCK_ATM=1,ISTOCK_LND=2,ISTOCK_ICE=3,ISTOCK_OCN=4
  75   character(len=3), parameter, dimension(NCOMPS)  :: COMP_NAMES=(/'ATM', 'LND', 'ICE', 'OCN'/)
  76
  77
  78   character(len=5) , parameter, dimension(NELEMS)  :: STOCK_NAMES=(/'water', 'heat ', 'salt '/)
  79   character(len=12), parameter, dimension(NELEMS)  :: STOCK_UNITS=(/'[Kg]    ','[Joules]','[Kg]    '/)
  80
  81
  82 contains
  83
  84     !> Starts a stock report
  85     subroutine stocks_report_init(Time)
  86     type(time_type), intent(in) :: Time !< Model time
  87
  88     character(len=80)                :: formatString,space
  89     integer                          :: i,s
  90     real(r8_kind), dimension(NELEMS) :: val_atm, val_lnd, val_ice, val_ocn
  91
  92     ! Write the version of this file to the log file
  93     call write_version_number('STOCK_CONSTANTS_MOD', version)
  94
  95     do i = 1, NELEMS_report
  96        val_atm(i) = Atm_stock(i)%q_start
  97        val_lnd(i) = Lnd_stock(i)%q_start
  98        val_ice(i) = Ice_stock(i)%q_start
  99        val_ocn(i) = Ocn_stock(i)%q_start
 100        call mpp_sum(val_atm(i))
 101        call mpp_sum(val_lnd(i))
 102        call mpp_sum(val_ice(i))
 103        call mpp_sum(val_ocn(i))
 104     enddo
 105
 106
 107
 108     if(mpp_pe() == mpp_root_pe()) then
 109 !       earth_area = 4.*PI*Radius**2
 110
 111        write(stocks_file,*) '================Stocks Report Guide====================================='
 112        write(stocks_file,*) ' '
 113        write(stocks_file,*) 'S(t) = Total amount     of a tracer in the component model at time t.'
 114        write(stocks_file,*) '       Calculated via the component model itself.'
 115        write(stocks_file,*) ' '
 116        write(stocks_file,*) 'F(t) = Cumulative input of a tracer to the component model at time t.'
 117        write(stocks_file,*) '       Calculated via interchange of fluxes with other component models.'
 118        write(stocks_file,*) ' '
 119        write(stocks_file,*) 'S(t) - S(0) = Cumulative increase of the component stocks at time t'
 120        write(stocks_file,*) '              Calculated by the component itself.'
 121        write(stocks_file,*) ' '
 122        write(stocks_file,*) 'In a conserving component F(t)=S(t)-S(0) to within numerical accuracy.'
 123        write(stocks_file,*) ' '
 124        write(stocks_file,*) 'Component Model refers to one of OCN, ATM, LND or ICE'
 125        write(stocks_file,*) ''
 126        write(stocks_file,*) 'NOTE: When use_lag_fluxes=.true. is used in coupler, the ocean stocks '
 127        write(stocks_file,*) '      calculations are in error by an order which scales as the inverse'
 128        write(stocks_file,*) '      of the number of time steps.'
 129        write(stocks_file,*) ' '
 130        write(stocks_file,*) '======================================================================='
 131
 132
 133        write(stocks_file,*) '======================Initial Stock S(0)==============================='
 134 !The following produces  formatString='(5x,a,a,12x,a,a, 9x)' but is general to handle more elements
 135        formatString= '(5x'
 136        do i=1,NELEMS_report
 137           s = 25-len_trim(STOCK_NAMES(i))-len_trim(STOCK_UNITS(i))
 138           write(space,'(i2)') s
 139           formatString= trim(formatString)//',a,a,'//trim(space)
 140           formatString= trim(formatString)//trim('x')
 141        enddo
 142        formatString= trim(formatString)//')'
 143
 144        write(stocks_file,formatString) (trim(STOCK_NAMES(i)),trim(STOCK_UNITS(i)), i=1,NELEMS_report)
 145
 146 !The following produces  formatString=' (a,x,es22.15,3x,es22.15,3x)' but is general to handle more elements
 147        formatString= '(a,x'
 148        do i=1,NELEMS_report
 149           write(space,'(i2)') s
 150           formatString= trim(formatString)//',es22.15,3x'
 151        enddo
 152        formatString= trim(formatString)//')'
 153
 154
 155        write(stocks_file,formatString) 'ATM', (val_atm(i), i=1,NELEMS_report)
 156        write(stocks_file,formatString) 'LND', (val_lnd(i), i=1,NELEMS_report)
 157        write(stocks_file,formatString) 'ICE', (val_ice(i), i=1,NELEMS_report)
 158        write(stocks_file,formatString) 'OCN', (val_ocn(i), i=1,NELEMS_report)
 159
 160        write(stocks_file,*) '========================================================================'
 161        write(stocks_file,'(a)'  ) ' '!blank line
 162
 163     end if
 164
 165     call stocks_set_init_time(Time)
 166
 167   end subroutine stocks_report_init
 168
 169   !> Writes update to stock report
 170   subroutine stocks_report(Time)
 171     type(time_type)               , intent(in) :: Time !< Model time
 172
 173     type(stock_type)                 :: stck
 174     real(r8_kind), dimension(NCOMPS) :: f_value, f_ice_grid, f_ocn_grid, f_ocn_btf, q_start, q_now,c_value
 175     character(len=80)                :: formatString
 176     integer                          :: iday0, isec0, iday, isec, hours
 177     real(r8_kind)                    :: days
 178     integer                          :: diagID , comp,elem,i
 179     integer, parameter :: initID = -2 !< initial value for diag IDs. Must not be equal to the value
 180                                       !! that register_diag_field returns when it can't register
 181                                       !! the filed -- otherwise the registration
 182                                       !! is attempted every time this subroutine is called
 183
 184     integer, dimension(NCOMPS,NELEMS), save :: f_valueDiagID = initID
 185     integer, dimension(NCOMPS,NELEMS), save :: c_valueDiagID = initID
 186     integer, dimension(NCOMPS,NELEMS), save :: fmc_valueDiagID = initID
 187     integer, dimension(NCOMPS,NELEMS), save :: f_lostDiagID = initID
 188
 189     real(r8_kind)     :: diagField
 190     logical           :: used
 191     character(len=30) :: field_name, units
 192
 193     if(mpp_pe()==mpp_root_pe()) then
 194        call get_time(init_time, isec0, iday0)
 195        call get_time(Time, isec, iday)
 196
 197        hours = iday*24 + isec/3600 - iday0*24 - isec0/3600
 198        days  = real(hours, r8_kind)/24.0_r8_kind
 199        write(stocks_file,*) '==============================================='
 200        write(stocks_file,'(a,f12.3)') 't = TimeSinceStart[days]= ',days
 201        write(stocks_file,*) '==============================================='
 202     endif
 203
 204     do elem = 1,NELEMS_report
 205
 206        do comp = 1,NCOMPS
 207
 208           if(comp == ISTOCK_ATM) stck = Atm_stock(elem)
 209           if(comp == ISTOCK_LND) stck = Lnd_stock(elem)
 210           if(comp == ISTOCK_ICE) stck = Ice_stock(elem)
 211           if(comp == ISTOCK_OCN) stck = Ocn_stock(elem)
 212
 213
 214           f_ice_grid(comp) = sum(stck%dq)
 215           f_ocn_grid(comp) = sum(stck%dq_IN)
 216           f_ocn_btf(comp)  = stck%dq_IN( ISTOCK_BOTTOM )
 217
 218           q_start(comp) = stck%q_start
 219           q_now(comp)   = stck%q_now
 220
 221           call mpp_sum(f_ice_grid(comp))
 222           call mpp_sum(f_ocn_grid(comp))
 223           call mpp_sum(f_ocn_btf(comp))
 224           call mpp_sum(q_start(comp))
 225           call mpp_sum(q_now(comp))
 226
 227           c_value(comp) = q_now(comp) - q_start(comp)
 228
 229           if(mpp_pe() == mpp_root_pe()) then
 230
 231              if(f_valueDiagID(comp,elem) == initID) then
 232                 field_name = trim(COMP_NAMES(comp)) // trim(STOCK_NAMES(elem))
 233                 field_name  = trim(field_name) // 'StocksChange_Flux'
 234                 units = trim(STOCK_UNITS(elem))
 235                 f_valueDiagID(comp,elem) = register_diag_field('stock_print', field_name, Time, &
 236                      units=units)
 237              endif
 238
 239              if(c_valueDiagID(comp,elem) == initID) then
 240                 field_name = trim(COMP_NAMES(comp)) // trim(STOCK_NAMES(elem))
 241                 field_name = trim(field_name) // 'StocksChange_Comp'
 242                 units = trim(STOCK_UNITS(elem))
 243                 c_valueDiagID(comp,elem) = register_diag_field('stock_print', field_name, Time, &
 244                      units=units)
 245              endif
 246
 247              if(fmc_valueDiagID(comp,elem) == initID) then
 248                 field_name = trim(COMP_NAMES(comp)) // trim(STOCK_NAMES(elem))
 249                 field_name = trim(field_name) // 'StocksChange_Diff'
 250                 units = trim(STOCK_UNITS(elem))
 251                 fmc_valueDiagID(comp,elem) = register_diag_field('stock_print', field_name, Time, &
 252                      units=units)
 253              endif
 254
 255              f_value(comp) = f_ice_grid(comp)
 256
 257              if(comp == ISTOCK_OCN) then
 258
 259                 f_value(comp) = f_ocn_grid(comp)
 260
 261                 if(f_lostDiagID(comp,elem) == initID) then
 262                    field_name = trim(COMP_NAMES(comp)) // trim(STOCK_NAMES(elem))
 263                    field_name = trim(field_name) // 'StocksExchangeLost'
 264                    units = trim(STOCK_UNITS(elem))
 265                    f_lostDiagID(comp,elem) = register_diag_field('stock_print', field_name, Time, &
 266                         units=units)
 267                 endif
 268
 269                 DiagID=f_lostDiagID(comp,elem)
 270                 diagField = f_ice_grid(comp) - f_ocn_grid(comp)
 271                 if (DiagID > 0)  used = send_data(DiagID, diagField, Time)
 272
 273              endif
 274
 275
 276              DiagID=f_valueDiagID(comp,elem)
 277              diagField = f_value(comp)
 278              if (DiagID > 0)  used = send_data(DiagID, diagField, Time)
 279              DiagID=c_valueDiagID(comp,elem)
 280              diagField = c_value(comp)
 281              if (DiagID > 0)  used = send_data(DiagID, diagField, Time)
 282              DiagID=fmc_valueDiagID(comp,elem)
 283              diagField = f_value(comp)-c_value(comp)
 284              if (DiagID > 0)  used = send_data(DiagID, diagField, Time)
 285
 286
 287              ! formatString = '(a,a,a,i16,2x,es22.15,2x,es22.15,2x,es22.15,2x,es22.15,2x,es22.15,2x,es22.15)'
 288              !
 289              !             write(stocks_file,formatString) trim(COMP_NAMES(comp)),STOCK_NAMES(elem),STOCK_UNITS(elem) &
 290              !                  ,hours, q_now, q_now-q_start, f_value, f_value - (q_now - q_start),
 291              !                  (f_value - (q_now - q_start))/q_start
 292
 293
 294           endif
 295        enddo
 296
 297
 298        if(mpp_pe()==mpp_root_pe()) then
 299 !          write(stocks_file,'(a)'  ) ' '!blank line
 300 !          write(stocks_file,'(a,f12.3)') 't = TimeSinceStart[days]= ',days
 301 !          write(stocks_file,'(a)'  )   ' '!blank line
 302 !          write(stocks_file,'(a,30x,a,20x,a,20x,a,20x,a)') 'Component ','ATM','LND','ICE','OCN'
 303 !          write(stocks_file,'(55x,a,20x,a,20x,a,20x,a)')  'ATM','LND','ICE','OCN'
 304 !          write(stocks_file,'(a,f12.3,12x,a,20x,a,20x,a,20x,a)') 't = TimeSinceStart[days]=
 305 !          ',days,'ATM','LND','ICE','OCN'
 306
 307           write(stocks_file,'(a,a,40x,a,20x,a,20x,a,20x,a)') 'Stocks of ',trim(STOCK_NAMES(elem)),'ATM','LND', &
 308                & 'ICE','OCN'
 309           formatString = '(a,a,2x,es22.15,2x,es22.15,2x,es22.15,2x,es22.15)'
 310
 311           write(stocks_file,formatString) 'Total =S(t)               ',STOCK_UNITS(elem),&
 312                ( q_now(i), i=1,NCOMPS)
 313           write(stocks_file,formatString) 'Change=S(t)-S(0)          ',STOCK_UNITS(elem),&
 314                ( q_now(i)-q_start(i), i=1,NCOMPS)
 315           write(stocks_file,formatString) 'Input =F(t)               ',STOCK_UNITS(elem),&
 316                ( f_value(i), i=1,NCOMPS)
 317           write(stocks_file,formatString) 'Diff  =F(t) - (S(t)-S(0)) ',STOCK_UNITS(elem),&
 318                ( f_value(i) - c_value(i), i=1,NCOMPS)
 319           write(stocks_file,formatString) 'Error =Diff/S(0)          ','[NonDim]    ', &
 320                ((f_value(i) - c_value(i))/(1+q_start(i)), i=1,NCOMPS)  !added 1 to avoid div by zero.
 321                                                                        !! Assuming q_start large
 322
 323           write(stocks_file,'(a)'  ) ' '!blank line
 324           formatString = '(a,a,a,6x,es22.15)'
 325           write(stocks_file,formatString) 'Lost Stocks in the exchange between Ice and Ocean ',trim(STOCK_NAMES(elem))&
 326                      &,trim(STOCK_UNITS(elem)), f_ice_grid(ISTOCK_OCN) - f_ocn_grid(ISTOCK_OCN) + f_ocn_btf(ISTOCK_OCN)
 327
 328           write(stocks_file,'(a)') ' ' !blank line
 329           write(stocks_file,'(a)') ' ' !blank line
 330
 331        endif
 332     enddo
 333
 334   end subroutine stocks_report
 335
 336   subroutine stocks_set_init_time(Time)
 337     type(time_type)     , intent(in) :: Time !< init time to set for stock report
 338     init_time = Time
 339
 340   end subroutine stocks_set_init_time
 341
 342 end module stock_constants_mod
 343 !> @}
 344 ! close documentation grouping