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