ROMS
Loading...
Searching...
No Matches
ad_omega.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE ad_omega_mod
3#if defined ADJOINT && defined SOLVE3D
4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md !
10!=======================================================================
11! !
12! This routine computes S-coordinate vertical velocity (m^3/s), !
13! !
14! W=[Hz/(m*n)]*omega, !
15! !
16! diagnostically at horizontal RHO-points and vertical W-points. !
17! !
18! BASIC STATE variables needed: Huon, Hvom, z_w. !
19! !
20!=======================================================================
21!
22 implicit none
23!
24 PRIVATE
25 PUBLIC :: ad_omega
26!
27 CONTAINS
28!
29!***********************************************************************
30 SUBROUTINE ad_omega (ng, tile, model)
31!***********************************************************************
32!
33 USE mod_param
34 USE mod_grid
35 USE mod_ocean
36# if defined SEDIMENT && defined SED_MORPH
37 USE mod_sedbed
38 USE mod_stepping
39# endif
40!
41! Imported variable declarations.
42!
43 integer, intent(in) :: ng, tile, model
44!
45! Local variable declarations.
46!
47 character (len=*), parameter :: myfile = &
48 & __FILE__
49!
50# include "tile.h"
51!
52# ifdef PROFILE
53 CALL wclock_on (ng, model, 13, __line__, myfile)
54# endif
55 CALL ad_omega_tile (ng, tile, model, &
56 & lbi, ubi, lbj, ubj, &
57 & imins, imaxs, jmins, jmaxs, &
58# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
59 & nstp(ng), nnew(ng), &
60 & grid(ng) % omn, &
61 & sedbed(ng) % bed_thick, &
62 & sedbed(ng) % ad_bed_thick, &
63# endif
64 & grid(ng) % Huon, &
65 & grid(ng) % Hvom, &
66 & grid(ng) % z_w, &
67 & grid(ng) % ad_Huon, &
68 & grid(ng) % ad_Hvom, &
69 & grid(ng) % ad_z_w, &
70 & ocean(ng) % W, &
71 & ocean(ng) % ad_W_sol, &
72 & ocean(ng) % ad_W)
73# ifdef PROFILE
74 CALL wclock_off (ng, model, 13, __line__, myfile)
75# endif
76!
77 RETURN
78 END SUBROUTINE ad_omega
79!
80!***********************************************************************
81 SUBROUTINE ad_omega_tile (ng, tile, model, &
82 & LBi, UBi, LBj, UBj, &
83 & IminS, ImaxS, JminS, JmaxS, &
84# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
85 & nstp, nnew, &
86 & omn,
87 & bed_thick, ad_bed_thick, &
88# endif
89 & Huon, Hvom, z_w, &
90 & ad_Huon, ad_Hvom, ad_z_w, &
91 & W, ad_W_sol, ad_W)
92!***********************************************************************
93!
94 USE mod_param
95 USE mod_scalars
96 USE mod_sources
97!
98 USE ad_bc_3d_mod, ONLY : ad_bc_w3d_tile
99 USE bc_3d_mod, ONLY : bc_w3d_tile
100# ifdef DISTRIBUTE
101 USE mp_exchange_mod, ONLY : ad_mp_exchange3d
102 USE mp_exchange_mod, ONLY : mp_exchange3d
103# endif
104!
105! Imported variable declarations.
106!
107 integer, intent(in) :: ng, tile, model
108 integer, intent(in) :: LBi, UBi, LBj, UBj
109 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
110# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
111 integer, intent(in) :: nstp, nnew
112# endif
113!
114# ifdef ASSUMED_SHAPE
115 real(r8), intent(in) :: Huon(LBi:,LBj:,:)
116 real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
117 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
118# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
119 real(r8), intent(in) :: omn(LBi:,LBj:)
120 real(r8), intent(in):: bed_thick(LBi:,LBj:,:)
121 real(r8), intent(inout):: ad_bed_thick(LBi:,LBj:,:)
122# endif
123 real(r8), intent(inout) :: ad_Huon(LBi:,LBj:,:)
124 real(r8), intent(inout) :: ad_Hvom(LBi:,LBj:,:)
125 real(r8), intent(inout) :: ad_z_w(LBi:,LBj:,0:)
126 real(r8), intent(inout) :: ad_W(LBi:,LBj:,0:)
127 real(r8), intent(inout) :: ad_W_sol(LBi:,LBj:,0:)
128
129 real(r8), intent(out) :: W(LBi:,LBj:,0:)
130
131# else
132
133 real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
134 real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
135 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
136# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
137 real(r8), intent(in) :: omn(LBi:UBi,LBj:UBj)
138 real(r8), intent(in):: bed_thick(LBi:UBi,LBj:UBj,2)
139 real(r8), intent(inout):: ad_bed_thick(LBi:UBi,LBj:UBj,2)
140# endif
141 real(r8), intent(inout) :: ad_Huon(LBi:UBi,LBj:UBj,N(ng))
142 real(r8), intent(inout) :: ad_Hvom(LBi:UBi,LBj:UBj,N(ng))
143 real(r8), intent(inout) :: ad_z_w(LBi:UBi,LBj:UBj,0:N(ng))
144 real(r8), intent(inout) :: ad_W(LBi:UBi,LBj:UBj,0:N(ng))
145 real(r8), intent(inout) :: ad_W_sol(LBi:UBi,LBj:UBj,0:N(ng))
146
147 real(r8), intent(out) :: W(LBi:UBi,LBj:UBj,0:N(ng))
148# endif
149!
150! Local variable declarations.
151!
152 integer :: i, ii, is, j, jj, k
153 real(r8) :: cff
154 real(r8) :: ad_cff, adfac
155# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
156 real(r8) :: cff1
157# endif
158 real(r8), dimension(IminS:ImaxS) :: wrk
159 real(r8), dimension(IminS:ImaxS) :: ad_wrk
160
161# include "set_bounds.h"
162!
163!-----------------------------------------------------------------------
164! Initialize adjoint private variables.
165!-----------------------------------------------------------------------
166!
167 ad_cff=0.0_r8
168 DO i=imins,imaxs
169 ad_wrk(i)=0.0_r8
170 END DO
171!
172!-----------------------------------------------------------------------
173! Vertically integrage horizontal mass flux divergence.
174!-----------------------------------------------------------------------
175!
176! Save current adjoint solution for IO purposes.
177!
178 DO k=0,n(ng)
179 DO j=jstrr,jendr
180 DO i=istrr,iendr
181 ad_w_sol(i,j,k)=ad_w(i,j,k)
182 END DO
183 END DO
184 END DO
185!
186! Set lateral boundary conditions.
187!
188# ifdef DISTRIBUTE
189!^ CALL mp_exchange3d (ng, tile, model, 1, &
190!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
191!^ & NghostPoints, &
192!^ & EWperiodic(ng), NSperiodic(ng), &
193!^ & tl_W)
194!^
195 CALL ad_mp_exchange3d (ng, tile, model, 1, &
196 & lbi, ubi, lbj, ubj, 0, n(ng), &
197 & nghostpoints, &
198 & ewperiodic(ng), nsperiodic(ng), &
199 & ad_w)
200# endif
201!^ CALL bc_w3d_tile (ng, tile, &
202!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
203!^ & tl_W)
204!^
205 CALL ad_bc_w3d_tile (ng, tile, &
206 & lbi, ubi, lbj, ubj, 0, n(ng), &
207 & ad_w)
208!
209! In order to insure zero vertical velocity at the free-surface,
210! subtract the vertical velocities of the moving S-coordinates
211! isosurfaces. These isosurfaces are proportional to d(zeta)/d(t).
212! The proportionaly coefficients are a linear function of the
213! S-coordinate with zero value at the bottom (k=0) and unity at
214! the free-surface (k=N).
215!
216! Notice that here we need to recompute the intermediate value
217! of W which is needed for wrk.
218!
219# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
220 cff1=1.0_r8/dt(ng)
221# endif
222 DO j=jstr,jend
223 DO i=istr,iend
224# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
225 w(i,j,0)=-cff1*(bed_thick(i,j,nstp)- &
226 & bed_thick(i,j,nnew))*omn(i,j)
227# else
228 w(i,j,0)=0.0_r8
229# endif
230
231 END DO
232 DO k=1,n(ng)
233 DO i=istr,iend
234 w(i,j,k)=w(i,j,k-1)- &
235 & (huon(i+1,j,k)-huon(i,j,k)+ &
236 & hvom(i,j+1,k)-hvom(i,j,k))
237 END DO
238 END DO
239!
240! Apply mass point sources (volume vertical influx) to the BASIC
241! STATE, if any.
242!
243! Overwrite W(Isrc,Jsrc,k) with the same divergence of Huon,Hvom as
244! above but add in point source Qsrc(k) and reaccumulate the vertical
245! sum to obtain the correct net Qbar given in user input - J. Levin
246! (Jupiter Intelligence Inc.) and J. Wilkin
247!
248! Dsrc(is) = 2, flow across grid cell w-face (positive or negative)
249!
250 IF (lwsrc(ng)) THEN
251 DO is=1,nsrc(ng)
252 IF (int(sources(ng)%Dsrc(is)).eq.2) THEN
253 ii=sources(ng)%Isrc(is)
254 jj=sources(ng)%Jsrc(is)
255 IF (((istrr.le.ii).and.(ii.le.iendr)).and. &
256 & ((jstrr.le.jj).and.(jj.le.jendr)).and. &
257 & (j.eq.jj)) THEN
258 DO k=1,n(ng)
259 w(ii,jj,k)=w(ii,jj,k-1)- &
260 & (huon(ii+1,jj,k)-huon(ii,jj,k)+ &
261 & hvom(ii,jj+1,k)-hvom(ii,jj,k))+ &
262 & sources(ng)%Qsrc(is,k)
263 END DO
264 END IF
265 END IF
266 END DO
267 END IF
268!
269 DO i=istr,iend
270 wrk(i)=w(i,j,n(ng))/(z_w(i,j,n(ng))-z_w(i,j,0))
271 END DO
272!
273! Starting with zero vertical velocity at the bottom, integrate
274! from the bottom (k=0) to the free-surface (k=N). The w(:,:,N)
275! contains the vertical velocity at the free-surface, d(zeta)/d(t).
276! Notice that barotropic mass flux divergence is not used directly.
277!
278 DO i=istr,iend
279!^ tl_W(i,j,N(ng))=0.0_r8
280!^
281 ad_w(i,j,n(ng))=0.0_r8
282 END DO
283 DO k=n(ng)-1,1,-1
284 DO i=istr,iend
285!^ tl_W(i,j,k)=tl_W(i,j,k)- &
286!^ & tl_wrk(i)*(z_w(i,j,k)-z_w(i,j,0))- &
287!^ & wrk(i)*(tl_z_w(i,j,k)-tl_z_w(i,j,0))
288!^
289 adfac=wrk(i)*ad_w(i,j,k)
290 ad_wrk(i)=ad_wrk(i)- &
291 & ad_w(i,j,k)*(z_w(i,j,k)-z_w(i,j,0))
292 ad_z_w(i,j,0)=ad_z_w(i,j,0)+adfac
293 ad_z_w(i,j,k)=ad_z_w(i,j,k)-adfac
294 END DO
295 END DO
296 DO i=istr,iend
297 cff=1.0_r8/(z_w(i,j,n(ng))-z_w(i,j,0))
298!^ tl_wrk(i)=tl_cff*W(i,j,N(ng))+cff*tl_W(i,j,N(ng))
299!^
300 ad_w(i,j,n(ng))=ad_w(i,j,n(ng))+cff*ad_wrk(i)
301 ad_cff=ad_cff+w(i,j,n(ng))*ad_wrk(i)
302 ad_wrk(i)=0.0_r8
303!^ tl_cff=-cff*cff*(tl_z_w(i,j,N(ng))-tl_z_w(i,j,0))
304!^
305 adfac=-cff*cff*ad_cff
306 ad_z_w(i,j,0 )=ad_z_w(i,j,0 )-adfac
307 ad_z_w(i,j,n(ng))=ad_z_w(i,j,n(ng))+adfac
308 ad_cff=0.0_r8
309 END DO
310!
311! Adjoint of apply mass point sources (volume vertical influx), if any.
312!
313! Dsrc(is) = 2, flow across grid cell w-face (positive or negative)
314!
315 IF (lwsrc(ng)) THEN
316 DO is=1,nsrc(ng)
317 IF (int(sources(ng)%Dsrc(is)).eq.2) THEN
318 ii=sources(ng)%Isrc(is)
319 jj=sources(ng)%Jsrc(is)
320 IF (((istrr.le.ii).and.(ii.le.iendr)).and. &
321 & ((jstrr.le.jj).and.(jj.le.jendr)).and. &
322 & (j.eq.jj)) THEN
323!^ DO k=1,N(ng)
324!^
325 DO k=n(ng),1,-1
326!^ tl_W(ii,jj,k)=tl_W(ii,jj,k-1)- &
327!^ & (tl_Huon(ii+1,jj,k)-tl_Huon(ii,jj,k)+ &
328!^ & tl_Hvom(ii,jj+1,k)-tl_Hvom(ii,jj,k))+ &
329!^ & SOURCES(ng)%tl_Qsrc(is,k)
330!^
331 ad_w(ii,jj,k-1)=ad_w(ii,jj,k-1)+ad_w(ii,jj,k)
332 ad_huon(ii ,jj,k)=ad_huon(ii ,jj,k)+ad_w(ii,jj,k)
333 ad_huon(ii+1,jj,k)=ad_huon(ii+1,jj,k)-ad_w(ii,jj,k)
334 ad_hvom(ii,jj ,k)=ad_hvom(ii,jj ,k)+ad_w(ii,jj,k)
335 ad_hvom(ii,jj+1,k)=ad_hvom(ii,jj+1,k)-ad_w(ii,jj,k)
336 sources(ng)%ad_Qsrc(is,k)=sources(ng)%ad_Qsrc(is,k)+ &
337 & ad_w(ii,jj,k)
338!
339! We have to do this here so that we do not double count ad_W(ii,jj,k)
340! in the next loop. It is okay here since the loop is counting
341! backwards in k, and ad_W is cleared after the net loop.
342!
343 ad_w(ii,jj,k)=0.0_r8
344 END DO
345 END IF
346 END IF
347 END DO
348 END IF
349!
350! Adjoint of code added to clear tl_W to be consistent with adjoint.
351!
352!^ DO k=1,N(ng)
353!^
354 DO k=n(ng),1,-1
355 DO i=istr,iend
356!^ tl_W(i,j,k)=tl_W(i,j,k-1)- &
357!^ & (tl_Huon(i+1,j,k)-tl_Huon(i,j,k)+ &
358!^ & tl_Hvom(i,j+1,k)-tl_Hvom(i,j,k))
359!^
360 ad_w(i,j,k-1)=ad_w(i,j,k-1)+ad_w(i,j,k)
361 ad_huon(i ,j,k)=ad_huon(i ,j,k)+ad_w(i,j,k)
362 ad_huon(i+1,j,k)=ad_huon(i+1,j,k)-ad_w(i,j,k)
363 ad_hvom(i,j ,k)=ad_hvom(i,j ,k)+ad_w(i,j,k)
364 ad_hvom(i,j+1,k)=ad_hvom(i,j+1,k)-ad_w(i,j,k)
365 END DO
366 END DO
367 DO k=1,n(ng)
368 DO i=istr,iend
369!^ tl_W(i,j,k)=0.0_r8
370!^
371 ad_w(i,j,k)=0.0_r8
372 END DO
373 END DO
374!
375! Adjoint of starting vertical velocity at the bottom.
376!
377 DO i=istr,iend
378# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
379!^ tl_W(i,j,0)=-cff1*(tl_bed_thick(i,j,nstp)- &
380!^ & tl_bed_thick(i,j,nnew))*omn(i,j)
381!^
382 adfac=cff1*omn(i,j)*ad_w(i,j,0)
383 ad_bed_thick(i,j,nnew)=ad_bed_thick(i,j,nnew)+adfac
384 ad_bed_thick(i,j,nstp)=ad_bed_thick(i,j,nstp)+adfac
385 ad_w(i,j,0)=0.0_r8
386# else
387!^ tl_W(i,j,0)=0.0_r8
388!^
389 ad_w(i,j,0)=0.0_r8
390# endif
391 END DO
392!
393! Complete the computation of BASIC STATE W here so that it is correct
394! for the remainder of the code.
395!
396 DO k=n(ng)-1,1,-1
397 DO i=istr,iend
398 w(i,j,k)=w(i,j,k)-wrk(i)*(z_w(i,j,k)-z_w(i,j,0))
399 END DO
400 END DO
401 DO i=istr,iend
402 w(i,j,n(ng))=0.0_r8
403 END DO
404 END DO
405!
406! Set lateral boundary conditions for basic state.
407!
408 CALL bc_w3d_tile (ng, tile, &
409 & lbi, ubi, lbj, ubj, 0, n(ng), &
410 & w)
411# ifdef DISTRIBUTE
412 CALL mp_exchange3d (ng, tile, model, 1, &
413 & lbi, ubi, lbj, ubj, 0, n(ng), &
414 & nghostpoints, &
415 & ewperiodic(ng), nsperiodic(ng), &
416 & w)
417# endif
418!
419 RETURN
420 END SUBROUTINE ad_omega_tile
421#endif
422 END MODULE ad_omega_mod
ad_bc_3d_mod
Definition ad_bc_3d.F:2
ad_bc_3d_mod::ad_bc_w3d_tile
subroutine ad_bc_w3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, ad_a)
Definition ad_bc_3d.F:733
ad_omega_mod
Definition ad_omega.F:2
ad_omega_mod::ad_omega_tile
subroutine ad_omega_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nstp, nnew, omn,
Definition ad_omega.F:87
ad_omega_mod::ad_omega
subroutine, public ad_omega(ng, tile, model)
Definition ad_omega.F:31
bc_3d_mod
Definition bc_3d.F:2
bc_3d_mod::bc_w3d_tile
subroutine bc_w3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition bc_3d.F:591
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_scalars
Definition mod_scalars.F:2
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::lwsrc
logical, dimension(:), allocatable lwsrc
Definition mod_scalars.F:543
mod_sedbed
Definition sedbed_mod.h:1
mod_sedbed::sedbed
type(t_sedbed), dimension(:), allocatable sedbed
Definition sedbed_mod.h:157
mod_sources
Definition mod_sources.F:2
mod_sources::sources
type(t_sources), dimension(:), allocatable sources
Definition mod_sources.F:90
mod_sources::nsrc
integer, dimension(:), allocatable nsrc
Definition mod_sources.F:97
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::ad_mp_exchange3d
subroutine ad_mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, ad_a, ad_b, ad_c, ad_d)
Definition mp_exchange.F:5451
mp_exchange_mod::mp_exchange3d
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:2037
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