ROMS
Loading...
Searching...
No Matches
rp_diag.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE rp_diag_mod
3#ifdef TL_IOMS
4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md !
10!=======================================================================
11! !
12! This routine computes various representer tangent linear diagnostic !
13! fields. !
14! !
15!=======================================================================
16!
17 implicit none
18!
19 PRIVATE
20 PUBLIC :: rp_diag
21!
22 CONTAINS
23!
24!***********************************************************************
25 SUBROUTINE rp_diag (ng, tile)
26!***********************************************************************
27!
28 USE mod_param
29 USE mod_grid
30 USE mod_ocean
31 USE mod_stepping
32!
33! Imported variable declarations.
34!
35 integer, intent(in) :: ng, tile
36!
37! Local variable declarations.
38!
39 character (len=*), parameter :: myfile = &
40 & __FILE__
41!
42# include "tile.h"
43!
44# ifdef PROFILE
45 CALL wclock_on (ng, irpm, 7, __line__, myfile)
46# endif
47 CALL rp_diag_tile (ng, tile, &
48 & lbi, ubi, lbj, ubj, &
49 & imins, imaxs, jmins, jmaxs, &
50 & nstp(ng), kstp(ng), &
51 & grid(ng) % h, &
52 & grid(ng) % omn, &
53# ifdef SOLVE3D
54 & grid(ng) % tl_Hz, &
55 & grid(ng) % tl_z_r, &
56 & grid(ng) % tl_z_w, &
57 & ocean(ng) % tl_rho, &
58 & ocean(ng) % tl_u, &
59 & ocean(ng) % tl_v, &
60# endif
61 & ocean(ng) % tl_ubar, &
62 & ocean(ng) % tl_vbar, &
63 & ocean(ng) % tl_zeta)
64# ifdef PROFILE
65 CALL wclock_off (ng, irpm, 7, __line__, myfile)
66# endif
67!
68 RETURN
69 END SUBROUTINE rp_diag
70!
71!***********************************************************************
72 SUBROUTINE rp_diag_tile (ng, tile, &
73 & LBi, UBi, LBj, UBj, &
74 & IminS, ImaxS, JminS, JmaxS, &
75 & nstp, kstp, &
76 & h, omn, &
77# ifdef SOLVE3D
78 & tl_Hz, tl_z_r, tl_z_w, &
79 & tl_rho, tl_u, tl_v, &
80# endif
81 & tl_ubar, tl_vbar, tl_zeta)
82!***********************************************************************
83!
84 USE mod_param
85 USE mod_parallel
86 USE mod_iounits
87 USE mod_scalars
88
89# ifdef DISTRIBUTE
90!
91 USE distribute_mod, ONLY : mp_reduce
92# endif
93!
94 implicit none
95!
96! Imported variable declarations.
97!
98 integer, intent(in) :: ng, tile
99 integer, intent(in) :: LBi, UBi, LBj, UBj
100 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
101 integer, intent(in) :: nstp, kstp
102!
103# ifdef ASSUMED_SHAPE
104 real(r8), intent(in) :: h(LBi:,LBj:)
105 real(r8), intent(in) :: omn(LBi:,LBj:)
106# ifdef SOLVE3D
107 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
108 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
109 real(r8), intent(in) :: tl_z_w(LBi:,LBj:,0:)
110 real(r8), intent(in) :: tl_rho(LBi:,LBj:,:)
111 real(r8), intent(in) :: tl_u(LBi:,LBj:,:,:)
112 real(r8), intent(in) :: tl_v(LBi:,LBj:,:,:)
113# endif
114 real(r8), intent(in) :: tl_ubar(LBi:,LBj:,:)
115 real(r8), intent(in) :: tl_vbar(LBi:,LBj:,:)
116 real(r8), intent(in) :: tl_zeta(LBi:,LBj:,:)
117# else
118 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
119 real(r8), intent(in) :: omn(LBi:UBi,LBj:UBj)
120# ifdef SOLVE3D
121 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
122 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,N(ng))
123 real(r8), intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:N(ng))
124 real(r8), intent(in) :: tl_rho(LBi:UBi,LBj:UBj,N(ng))
125 real(r8), intent(in) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
126 real(r8), intent(in) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
127# endif
128 real(r8), intent(in) :: tl_ubar(LBi:UBi,LBj:UBj,:)
129 real(r8), intent(in) :: tl_vbar(LBi:UBi,LBj:UBj,:)
130 real(r8), intent(in) :: tl_zeta(LBi:UBi,LBj:UBj,:)
131# endif
132!
133! Local variable declarations.
134!
135 integer :: NSUB, i, j, k, trd
136
137# ifdef DISTRIBUTE
138 real(r8), dimension(3) :: rbuffer
139 character (len=3), dimension(3) :: op_handle
140# else
141 integer :: my_threadnum
142# endif
143
144 real(r8) :: cff, my_avgke, my_avgpe, my_volume
145
146 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ke2d
147 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: pe2d
148
149 character (len=8) :: kechar, pechar
150
151# include "set_bounds.h"
152!
153!-----------------------------------------------------------------------
154! Compute and report out volume averaged kinetic, potential total
155! energy, and volume of tangent linear fields. Notice that the
156! proper adjoint of these quantities are not coded.
157!-----------------------------------------------------------------------
158!
159 IF (mod(iic(ng)-1,ninfo(ng)).eq.0) THEN
160 DO j=jstr,jend
161# ifdef SOLVE3D
162 DO i=istr,iend
163 ke2d(i,j)=0.0_r8
164 pe2d(i,j)=0.5_r8*g*tl_z_w(i,j,n(ng))*tl_z_w(i,j,n(ng))
165 END DO
166 cff=g/rho0
167 DO k=n(ng),1,-1
168 DO i=istr,iend
169 ke2d(i,j)=ke2d(i,j)+ &
170 & tl_hz(i,j,k)* &
171 & 0.25_r8*(tl_u(i ,j,k,nstp)*tl_u(i ,j,k,nstp)+ &
172 & tl_u(i+1,j,k,nstp)*tl_u(i+1,j,k,nstp)+ &
173 & tl_v(i,j ,k,nstp)*tl_v(i,j ,k,nstp)+ &
174 & tl_v(i,j+1,k,nstp)*tl_v(i,j+1,k,nstp))
175 pe2d(i,j)=pe2d(i,j)+ &
176 & cff*tl_hz(i,j,k)*(tl_rho(i,j,k)+1000.0_r8)* &
177 & (tl_z_r(i,j,k)-tl_z_w(i,j,0))
178 END DO
179 END DO
180# else
181 cff=0.5_r8*g
182 DO i=istr,iend
183 ke2d(i,j)=(tl_zeta(i,j,kstp)+h(i,j))* &
184 & 0.25_r8*(tl_ubar(i ,j,kstp)*tl_ubar(i ,j,kstp)+ &
185 & tl_ubar(i+1,j,kstp)*tl_ubar(i+1,j,kstp)+ &
186 & tl_vbar(i,j ,kstp)*tl_vbar(i,j ,kstp)+ &
187 & tl_vbar(i,j+1,kstp)*tl_vbar(i,j+1,kstp))
188 pe2d(i,j)=cff*tl_zeta(i,j,kstp)*tl_zeta(i,j,kstp)
189 END DO
190# endif
191 END DO
192!
193! Integrate horizontally within one tile. In order to reduce the
194! round-off errors, the summation is performed in two stages. First,
195! the index j is collapsed and then the accumulation is carried out
196! along index i. In this order, the partial sums consist on much
197! fewer number of terms than in a straightforward two-dimensional
198! summation. Thus, adding numbers which are orders of magnitude
199! apart is avoided.
200!
201 DO i=istr,iend
202 pe2d(i,jend+1)=0.0_r8
203 pe2d(i,jstr-1)=0.0_r8
204 ke2d(i,jstr-1)=0.0_r8
205 END DO
206 DO j=jstr,jend
207 DO i=istr,iend
208# ifdef SOLVE3D
209!! pe2d(i,Jend+1)=pe2d(i,Jend+1)+ &
210!! & omn(i,j)*(z_w(i,j,N(ng))-z_w(i,j,0))
211# else
212!! pe2d(i,Jend+1)=pe2d(i,Jend+1)+ &
213!! & omn(i,j)*(zeta(i,j,kstp)+h(i,j))
214# endif
215 pe2d(i,jend+1)=pe2d(i,jend+1)+omn(i,j)*h(i,j)
216 pe2d(i,jstr-1)=pe2d(i,jstr-1)+omn(i,j)*pe2d(i,j)
217 ke2d(i,jstr-1)=ke2d(i,jstr-1)+omn(i,j)*ke2d(i,j)
218 END DO
219 END DO
220 my_volume=0.0_r8
221 my_avgpe=0.0_r8
222 my_avgke=0.0_r8
223 DO i=istr,iend
224 my_volume=my_volume+pe2d(i,jend+1)
225 my_avgpe =my_avgpe +pe2d(i,jstr-1)
226 my_avgke =my_avgke +ke2d(i,jstr-1)
227 END DO
228!
229! Perform global summation: whoever gets first to the critical region
230! resets global sums before global summation starts; after the global
231! summation is completed, thread, which is the last one to enter the
232! critical region, finalizes the computation of diagnostics and prints
233! them out.
234!
235# ifdef DISTRIBUTE
236 nsub=1 ! distributed-memory
237# else
238 IF (domain(ng)%SouthWest_Corner(tile).and. &
239 & domain(ng)%NorthEast_Corner(tile)) THEN
240 nsub=1 ! non-tiled application
241 ELSE
242 nsub=ntilex(ng)*ntilee(ng) ! tiled application
243 END IF
244# endif
245!$OMP CRITICAL (RP_DIAGNOSTICS)
246 IF (tile_count.eq.0) THEN
247 volume=0.0_r8
248 avgke=0.0_r8
249 avgpe=0.0_r8
250 END IF
251 volume=volume+my_volume
252 avgke=avgke+my_avgke
253 avgpe=avgpe+my_avgpe
254 tile_count=tile_count+1
255 IF (tile_count.eq.nsub) THEN
256 tile_count=0
257# ifdef DISTRIBUTE
258 rbuffer(1)=volume
259 rbuffer(2)=avgke
260 rbuffer(3)=avgpe
261 op_handle(1)='SUM'
262 op_handle(2)='SUM'
263 op_handle(3)='SUM'
264 CALL mp_reduce (ng, irpm, 3, rbuffer, op_handle)
265 volume=rbuffer(1)
266 avgke=rbuffer(2)
267 avgpe=rbuffer(3)
268 trd=mymaster
269# else
270 trd=my_threadnum()
271# endif
272 avgke=avgke/volume
273 avgpe=avgpe/volume
274 avgkp=avgke+avgpe
275 IF (first_time(ng).eq.0) THEN
276 first_time(ng)=1
277 IF (master.and.(ng.eq.1)) THEN
278 WRITE (stdout,10) 'TIME-STEP', 'YYYY-MM-DD hh:mm:ss.ss', &
279 & 'KINETIC_ENRG', 'POTEN_ENRG', &
280# ifdef NESTING
281 & 'TOTAL_ENRG', 'NET_VOLUME', 'Grid'
282# else
283 & 'TOTAL_ENRG', 'NET_VOLUME'
284# endif
285# ifdef NESTING
286 10 FORMAT (/,1x,a,1x,a,2x,a,3x,a,4x,a,4x,a,2x,a)
287# else
288 10 FORMAT (/,1x,a,1x,a,2x,a,3x,a,4x,a,4x,a)
289# endif
290 END IF
291 END IF
292 IF (master) THEN ! restart counter after 10 billion steps
293 WRITE (stdout,20) int(mod(real(iic(ng)-1,r8),1.0e+10_r8)), &
294 & time_code(ng), &
295# ifdef NESTING
296 & avgke, avgpe, avgkp, volume, ng
297# else
298 & avgke, avgpe, avgkp, volume
299# endif
300# ifdef NESTING
301 20 FORMAT (i10,1x,a,4(1pe14.6),2x,i2.2)
302# else
303 20 FORMAT (i10,1x,a,4(1pe14.6))
304# endif
305 FLUSH (stdout)
306 END IF
307!
308! If blowing-up, set exit_flag to stop computations.
309!
310 WRITE (kechar,'(1pe8.1)') avgke
311 WRITE (pechar,'(1pe8.1)') avgpe
312 DO i=1,8
313 IF ((kechar(i:i).eq.'N').or.(pechar(i:i).eq.'N').or. &
314 & (kechar(i:i).eq.'n').or.(pechar(i:i).eq.'n').or. &
315 & (kechar(i:i).eq.'*').or.(pechar(i:i).eq.'*')) THEN
316 exit_flag=1
317 END IF
318 END DO
319 END IF
320!$OMP END CRITICAL (RP_DIAGNOSTICS)
321 END IF
322!
323 RETURN
324 END SUBROUTINE rp_diag_tile
325#endif
326 END MODULE rp_diag_mod
distribute_mod::mp_reduce
Definition distribute.F:113
distribute_mod
Definition distribute.F:3
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_iounits
Definition mod_iounits.F:2
mod_iounits::stdout
integer stdout
Definition mod_iounits.F:195
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_parallel
Definition mod_parallel.F:2
mod_parallel::master
logical master
Definition mod_parallel.F:40
mod_parallel::mymaster
integer mymaster
Definition mod_parallel.F:124
mod_parallel::tile_count
integer tile_count
Definition mod_parallel.F:89
mod_param
Definition mod_param.F:2
mod_param::irpm
integer, parameter irpm
Definition mod_param.F:664
mod_param::ntilex
integer, dimension(:), allocatable ntilex
Definition mod_param.F:685
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_param::ntilee
integer, dimension(:), allocatable ntilee
Definition mod_param.F:686
mod_scalars
Definition mod_scalars.F:2
mod_scalars::ninfo
integer, dimension(:), allocatable ninfo
Definition mod_scalars.F:1103
mod_scalars::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
mod_scalars::avgke
real(dp) avgke
Definition mod_scalars.F:662
mod_scalars::avgpe
real(dp) avgpe
Definition mod_scalars.F:663
mod_scalars::volume
real(dp) volume
Definition mod_scalars.F:665
mod_scalars::first_time
integer, dimension(:), allocatable first_time
Definition mod_scalars.F:660
mod_scalars::time_code
character(len=22), dimension(:), allocatable time_code
Definition mod_scalars.F:292
mod_scalars::exit_flag
integer exit_flag
Definition mod_scalars.F:559
mod_scalars::avgkp
real(dp) avgkp
Definition mod_scalars.F:664
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
mod_stepping
Definition mod_stepping.F:2
mod_stepping::kstp
integer, dimension(:), allocatable kstp
Definition mod_stepping.F:66
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
rp_diag_mod
Definition rp_diag.F:2
rp_diag_mod::rp_diag
subroutine, public rp_diag(ng, tile)
Definition rp_diag.F:26
rp_diag_mod::rp_diag_tile
subroutine rp_diag_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nstp, kstp, h, omn, tl_hz, tl_z_r, tl_z_w, tl_rho, tl_u, tl_v, tl_ubar, tl_vbar, tl_zeta)
Definition rp_diag.F:82
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3