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