13 implicit none ;
private
15 #include <MOM_memory.h>
26 sfc_flux, btm_flux, btm_reservoir, sink_rate, convert_flux_in)
29 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h_old
31 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: ea
33 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: eb
35 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: tr
36 real,
intent(in) :: dt
37 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(in) :: sfc_flux
41 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(in) :: btm_flux
44 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(inout) :: btm_reservoir
46 real,
optional,
intent(in) :: sink_rate
48 logical,
optional,
intent(in) :: convert_flux_in
53 real,
dimension(SZI_(G),SZJ_(G)) :: &
54 sfc_src, & !< The time-integrated surface source of the tracer [CU H ~> CU m or CU kg m-2].
56 real,
dimension(SZI_(G)) :: &
57 b1, & !< b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
59 real :: c1(szi_(g),szk_(gv))
60 real :: h_minus_dsink(szi_(g),szk_(gv))
63 real :: sink(szi_(g),szk_(gv)+1)
71 logical :: convert_flux = .true.
74 integer :: i, j, k, is, ie, js, je, nz
75 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
78 call mom_error(warning,
"MOM_tracer_diabatic.F90, tracer_vertdiff called "//&
79 "with only one vertical level")
83 if (
present(convert_flux_in)) convert_flux = convert_flux_in
84 h_neglect = gv%H_subroundoff
86 if (
present(sink_rate)) sink_dist = (dt*sink_rate) * gv%m_to_H
92 do j=js,je;
do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ;
enddo ;
enddo
93 if (
present(sfc_flux))
then
94 if (convert_flux)
then
96 do j = js, je;
do i = is,ie
97 sfc_src(i,j) = (sfc_flux(i,j)*dt) * gv%kg_m2_to_H
101 do j = js, je;
do i = is,ie
102 sfc_src(i,j) = sfc_flux(i,j)
106 if (
present(btm_flux))
then
107 if (convert_flux)
then
109 do j = js, je;
do i = is,ie
110 btm_src(i,j) = (btm_flux(i,j)*dt) * gv%kg_m2_to_H
114 do j = js, je;
do i = is,ie
115 btm_src(i,j) = btm_flux(i,j)
120 if (
present(sink_rate))
then
127 if (
present(btm_reservoir))
then
128 do i=is,ie ; sink(i,nz+1) = sink_dist ;
enddo
129 do k=2,nz ;
do i=is,ie
130 sink(i,k) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k)
133 do i=is,ie ; sink(i,nz+1) = 0.0 ;
enddo
135 do k=nz,2,-1 ;
do i=is,ie
136 if (sink(i,k+1) >= sink_dist)
then
137 sink(i,k) = sink_dist
138 h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,k+1) - sink(i,k))
139 elseif (sink(i,k+1) + h_old(i,j,k) < sink_dist)
then
140 sink(i,k) = sink(i,k+1) + h_old(i,j,k)
141 h_minus_dsink(i,k) = 0.0
143 sink(i,k) = sink_dist
144 h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,k+1)) - sink(i,k)
149 sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2))
153 do i=is,ie ;
if (g%mask2dT(i,j) > 0.5)
then
154 b_denom_1 = h_minus_dsink(i,1) + ea(i,j,1) + h_neglect
155 b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
156 d1(i) = b_denom_1 * b1(i)
157 h_tr = h_old(i,j,1) + h_neglect
158 tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j)
160 do k=2,nz-1 ;
do i=is,ie ;
if (g%mask2dT(i,j) > 0.5)
then
161 c1(i,k) = eb(i,j,k-1) * b1(i)
162 b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ea(i,j,k) + sink(i,k)) + &
164 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
165 d1(i) = b_denom_1 * b1(i)
166 h_tr = h_old(i,j,k) + h_neglect
167 tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + &
168 (ea(i,j,k) + sink(i,k)) * tr(i,j,k-1))
169 endif ;
enddo ;
enddo
170 do i=is,ie ;
if (g%mask2dT(i,j) > 0.5)
then
171 c1(i,nz) = eb(i,j,nz-1) * b1(i)
172 b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ea(i,j,nz) + sink(i,nz)) + &
174 b1(i) = 1.0 / (b_denom_1 + eb(i,j,nz))
175 h_tr = h_old(i,j,nz) + h_neglect
176 tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + &
177 (ea(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1))
179 if (
present(btm_reservoir))
then ;
do i=is,ie ;
if (g%mask2dT(i,j)>0.5)
then
180 btm_reservoir(i,j) = btm_reservoir(i,j) + &
181 (sink(i,nz+1)*tr(i,j,nz)) * gv%H_to_kg_m2
182 endif ;
enddo ;
endif
184 do k=nz-1,1,-1 ;
do i=is,ie ;
if (g%mask2dT(i,j) > 0.5)
then
185 tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
186 endif ;
enddo ;
enddo
191 do i=is,ie ;
if (g%mask2dT(i,j) > -0.5)
then
192 h_tr = h_old(i,j,1) + h_neglect
193 b_denom_1 = h_tr + ea(i,j,1)
194 b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
196 tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j)
199 do k=2,nz-1 ;
do i=is,ie ;
if (g%mask2dT(i,j) > -0.5)
then
200 c1(i,k) = eb(i,j,k-1) * b1(i)
201 h_tr = h_old(i,j,k) + h_neglect
202 b_denom_1 = h_tr + d1(i) * ea(i,j,k)
203 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
204 d1(i) = b_denom_1 * b1(i)
205 tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ea(i,j,k) * tr(i,j,k-1))
206 endif ;
enddo ;
enddo
207 do i=is,ie ;
if (g%mask2dT(i,j) > -0.5)
then
208 c1(i,nz) = eb(i,j,nz-1) * b1(i)
209 h_tr = h_old(i,j,nz) + h_neglect
210 b_denom_1 = h_tr + d1(i)*ea(i,j,nz)
211 b1(i) = 1.0 / ( b_denom_1 + eb(i,j,nz))
212 tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + &
213 ea(i,j,nz) * tr(i,j,nz-1))
215 do k=nz-1,1,-1 ;
do i=is,ie ;
if (g%mask2dT(i,j) > -0.5)
then
216 tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
217 endif ;
enddo ;
enddo
229 in_flux_optional, out_flux_optional, update_h_opt)
233 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: tr
234 real,
intent(in ) :: dt
235 type(
forcing),
intent(in ) :: fluxes
236 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
237 real,
intent(in ) :: evap_cfl_limit
240 real,
intent(in ) :: minimum_forcing_depth
242 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(in ) :: in_flux_optional
244 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(in) :: out_flux_optional
246 logical,
optional,
intent(in) :: update_h_opt
249 integer,
parameter :: maxgroundings = 5
250 integer :: numberofgroundings, iground(maxgroundings), jground(maxgroundings)
251 real :: h_limit_fluxes, iforcingdepthscale
252 real :: dthickness, dtracer
253 real :: fractionofforcing, hold, ithickness
254 real :: rivermixconst
255 real,
dimension(SZI_(G)) :: &
260 real,
dimension(SZI_(G), SZK_(G)) :: h2d, tr2d
261 real,
dimension(SZI_(G),SZJ_(G)) :: in_flux
263 real,
dimension(SZI_(G),SZJ_(G)) :: out_flux
265 real,
dimension(SZI_(G)) :: in_flux_1d, out_flux_1d
266 real :: hgrounding(maxgroundings)
269 integer :: i, j, is, ie, js, je, k, nz, n, nsw
270 character(len=45) :: mesg
272 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
275 if ( (.not.
associated(fluxes%netMassIn)) .or. (.not.
associated(fluxes%netMassOut)) )
return
277 in_flux(:,:) = 0.0 ; out_flux(:,:) = 0.0
278 if (
present(in_flux_optional))
then
279 do j=js,je ;
do i=is,ie
280 in_flux(i,j) = in_flux_optional(i,j)
283 if (
present(out_flux_optional))
then
284 do j=js,je ;
do i=is,ie
285 out_flux(i,j) = out_flux_optional(i,j)
289 if (
present(update_h_opt))
then
290 update_h = update_h_opt
295 numberofgroundings = 0
310 do k=1,nz ;
do i=is,ie
312 tr2d(i,k) = tr(i,j,k)
316 in_flux_1d(i) = in_flux(i,j)
317 out_flux_1d(i) = out_flux(i,j)
329 netmassout(i) = fluxes%netMassOut(i,j)
330 netmassin(i) = fluxes%netMassIn(i,j)
337 if (g%mask2dT(i,j)>0.)
then
343 dthickness = netmassin(i)
348 netmassin(i) = netmassin(i) - dthickness
349 dtracer = dtracer + in_flux_1d(i)
354 h2d(i,k) = h2d(i,k) + dthickness
355 if (h2d(i,k) > 0.0)
then
356 ithickness = 1.0/h2d(i,k)
358 if (dthickness /= 0. .or. dtracer /= 0.) tr2d(i,k) = (hold*tr2d(i,k)+ dtracer)*ithickness
369 iforcingdepthscale = 1. / max(gv%H_subroundoff, minimum_forcing_depth - netmassout(i) )
371 fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale)
376 if (-fractionofforcing*netmassout(i) > evap_cfl_limit*h2d(i,k))
then
377 fractionofforcing = -evap_cfl_limit*h2d(i,k)/netmassout(i)
381 dthickness = max( fractionofforcing*netmassout(i), -h2d(i,k) )
383 dtracer = fractionofforcing*out_flux_1d(i)
387 netmassout(i) = netmassout(i) - dthickness
388 out_flux_1d(i) = out_flux_1d(i) - dtracer
392 h2d(i,k) = h2d(i,k) + dthickness
393 if (h2d(i,k) > 0.)
then
394 ithickness = 1.0/h2d(i,k)
395 tr2d(i,k) = (hold*tr2d(i,k) + dtracer)*ithickness
403 if (abs(in_flux_1d(i))+abs(out_flux_1d(i)) /= 0.0)
then
405 numberofgroundings = numberofgroundings +1
406 if (numberofgroundings<=maxgroundings)
then
407 iground(numberofgroundings) = i
408 jground(numberofgroundings) = j
409 hgrounding(numberofgroundings) = abs(in_flux_1d(i))+abs(out_flux_1d(i))
417 do k=1,nz ;
do i=is,ie
418 tr(i,j,k) = tr2d(i,k)
422 do k=1,nz ;
do i=is,ie
429 if (numberofgroundings>0)
then
430 do i = 1, min(numberofgroundings, maxgroundings)
431 write(mesg(1:45),
'(3es15.3)') g%geoLonT( iground(i), jground(i) ), &
432 g%geoLatT( iground(i), jground(i)) , hgrounding(i)
433 call mom_error(warning,
"MOM_tracer_diabatic.F90, applyTracerBoundaryFluxesInOut(): "//&
434 "Tracer created. x,y,dh= "//trim(mesg), all_print=.true.)
437 if (numberofgroundings - maxgroundings > 0)
then
438 write(mesg,
'(i4)') numberofgroundings - maxgroundings
439 call mom_error(warning,
"MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//&
440 trim(mesg) //
" groundings remaining")