17 use mom_eos,
only : int_density_dz, int_specific_vol_dp
18 use mom_eos,
only : int_density_dz_generic_plm, int_density_dz_generic_ppm
19 use mom_eos,
only : int_spec_vol_dp_generic_plm
20 use mom_eos,
only : int_density_dz_generic, int_spec_vol_dp_generic
21 use mom_ale,
only : pressure_gradient_plm, pressure_gradient_ppm,
ale_cs
23 implicit none ;
private
25 #include <MOM_memory.h>
27 public pressureforce_blk_afv, pressureforce_blk_afv_init, pressureforce_blk_afv_end
28 public pressureforce_blk_afv_bouss, pressureforce_blk_afv_nonbouss
42 type(time_type),
pointer :: time
45 logical :: usemasswghtinterp
46 logical :: boundary_extrap
49 logical :: reconstruct
54 integer :: recon_scheme
58 integer :: id_e_tidal = -1
66 subroutine pressureforce_blk_afv(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
70 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
72 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: pfu
73 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: pfv
75 type(
ale_cs),
pointer :: ale_csp
76 real,
dimension(:,:),
optional,
pointer :: p_atm
78 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(out) :: pbce
81 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: eta
85 if (gv%Boussinesq)
then
86 call pressureforce_blk_afv_bouss(h, tv, pfu, pfv, g, gv, us, cs, ale_csp, p_atm, pbce, eta)
88 call pressureforce_blk_afv_nonbouss(h, tv, pfu, pfv, g, gv, us, cs, p_atm, pbce, eta)
91 end subroutine pressureforce_blk_afv
102 subroutine pressureforce_blk_afv_nonbouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
106 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
108 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: pfu
109 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: pfv
111 real,
dimension(:,:),
optional,
pointer :: p_atm
113 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(out) :: pbce
116 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: eta
120 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p
121 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
target :: &
126 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
131 real,
dimension(SZI_(G),SZJ_(G)) :: &
140 real,
dimension(SZDI_(G%Block(1)),SZDJ_(G%Block(1))) :: &
145 real,
dimension(SZI_(G)) :: rho_cv_bl
147 real,
dimension(SZDIB_(G%Block(1)),SZDJ_(G%Block(1))) :: &
150 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
152 real,
dimension(SZDI_(G%Block(1)),SZDJB_(G%Block(1))) :: &
155 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
157 real :: p_ref(szi_(g))
173 real :: rho_in_situ(szi_(g))
176 real,
parameter :: c1_6 = 1.0/6.0
177 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
178 integer :: is_bk, ie_bk, js_bk, je_bk, isq_bk, ieq_bk, jsq_bk, jeq_bk
179 integer :: i, j, k, n, ib, jb, ioff_bk, joff_bk
181 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
182 nkmb=gv%nk_rho_varies
183 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
185 if (.not.
associated(cs))
call mom_error(fatal, &
186 "MOM_PressureForce_AFV_nonBouss: Module must be initialized before it is used.")
189 if (
present(p_atm))
then ;
if (
associated(p_atm)) use_p_atm = .true. ;
endif
190 use_eos =
associated(tv%eqn_of_state)
192 dp_neglect = gv%H_to_Pa * gv%H_subroundoff
193 alpha_ref = 1.0/cs%Rho0
194 g_earth_z = us%L_T_to_m_s**2 * gv%g_Earth
195 i_gearth = 1.0 / g_earth_z
199 do j=jsq,jeq+1 ;
do i=isq,ieq+1
200 p(i,j,1) = p_atm(i,j)
204 do j=jsq,jeq+1 ;
do i=isq,ieq+1
209 do j=jsq,jeq+1 ;
do k=2,nz+1 ;
do i=isq,ieq+1
210 p(i,j,k) = p(i,j,k-1) + gv%H_to_Pa * h(i,j,k-1)
211 enddo ;
enddo ;
enddo
219 tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
220 tv_tmp%eqn_of_state => tv%eqn_of_state
221 do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ;
enddo
224 do k=1,nkmb ;
do i=isq,ieq+1
225 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
228 rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state)
229 do k=nkmb+1,nz ;
do i=isq,ieq+1
230 if (gv%Rlay(k) < rho_cv_bl(i))
then
231 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
233 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
238 tv_tmp%T => tv%T ; tv_tmp%S => tv%S
239 tv_tmp%eqn_of_state => tv%eqn_of_state
248 call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,k), &
249 p(:,:,k+1), alpha_ref, g%HI, tv%eqn_of_state, &
250 dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), &
251 inty_dza(:,:,k), bathyp=p(:,:,nz+1), dp_tiny=dp_neglect, &
252 usemasswghtinterp = cs%useMassWghtInterp)
254 alpha_anom = 1.0/gv%Rlay(k) - alpha_ref
255 do j=jsq,jeq+1 ;
do i=isq,ieq+1
256 dp(i,j) = gv%H_to_Pa * h(i,j,k)
257 dza(i,j,k) = alpha_anom * dp(i,j)
258 intp_dza(i,j,k) = 0.5 * alpha_anom * dp(i,j)**2
260 do j=js,je ;
do i=isq,ieq
261 intx_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i+1,j))
263 do j=jsq,jeq ;
do i=is,ie
264 inty_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i,j+1))
280 za(i,j) = alpha_ref*p(i,j,nz+1) - g_earth_z*g%bathyT(i,j)
282 do k=nz,1,-1 ;
do i=isq,ieq+1
283 za(i,j) = za(i,j) + dza(i,j,k)
290 do j=jsq,jeq+1 ;
do i=isq,ieq+1
291 ssh(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * i_gearth
293 call calc_tidal_forcing(cs%Time, ssh, e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
295 do j=jsq,jeq+1 ;
do i=isq,ieq+1
296 za(i,j) = za(i,j) - g_earth_z * e_tidal(i,j)
300 if (cs%GFS_scale < 1.0)
then
306 rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state)
309 dm(i,j) = (cs%GFS_scale - 1.0) * &
310 (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j))
315 do j=jsq,jeq+1 ;
do i=isq,ieq+1
316 dm(i,j) = (cs%GFS_scale - 1.0) * &
317 (p(i,j,1)*(1.0/gv%Rlay(1) - alpha_ref) + za(i,j))
335 is_bk=g%block(n)%isc ; ie_bk=g%block(n)%iec
336 js_bk=g%block(n)%jsc ; je_bk=g%block(n)%jec
337 isq_bk=g%block(n)%IscB ; ieq_bk=g%block(n)%IecB
338 jsq_bk=g%block(n)%JscB ; jeq_bk=g%block(n)%JecB
339 ioff_bk = g%Block(n)%idg_offset - g%HI%idg_offset
340 joff_bk = g%Block(n)%jdg_offset - g%HI%jdg_offset
341 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
342 i = ib+ioff_bk ; j = jb+joff_bk
343 za_bk(ib,jb) = za(i,j)
345 do jb=js_bk,je_bk ;
do ib=isq_bk,ieq_bk
346 i = ib+ioff_bk ; j = jb+joff_bk
347 intx_za_bk(ib,jb) = 0.5*(za_bk(ib,jb) + za_bk(ib+1,jb))
349 do jb=jsq_bk,jeq_bk ;
do ib=is_bk,ie_bk
350 i = ib+ioff_bk ; j = jb+joff_bk
351 inty_za_bk(ib,jb) = 0.5*(za_bk(ib,jb) + za_bk(ib,jb+1))
356 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
357 i = ib+ioff_bk ; j = jb+joff_bk
358 dp_bk(ib,jb) = gv%H_to_Pa*h(i,j,k)
359 za_bk(ib,jb) = za_bk(ib,jb) - dza(i,j,k)
361 do jb=js_bk,je_bk ;
do ib=isq_bk,ieq_bk
362 i = ib+ioff_bk ; j = jb+joff_bk
363 intx_za_bk(ib,jb) = intx_za_bk(ib,jb) - intx_dza(i,j,k)
364 pfu(i,j,k) = (((za_bk(ib,jb)*dp_bk(ib,jb) + intp_dza(i,j,k)) - &
365 (za_bk(ib+1,jb)*dp_bk(ib+1,jb) + intp_dza(i+1,j,k))) + &
366 ((dp_bk(ib+1,jb) - dp_bk(ib,jb)) * intx_za_bk(ib,jb) - &
367 (p(i+1,j,k) - p(i,j,k)) * intx_dza(i,j,k))) * &
368 (2.0*g%IdxCu(i,j) / ((dp_bk(ib,jb) + dp_bk(ib+1,jb)) + &
371 do jb=jsq_bk,jeq_bk ;
do ib=is_bk,ie_bk
372 i = ib+ioff_bk ; j = jb+joff_bk
373 inty_za_bk(ib,jb) = inty_za_bk(ib,jb) - inty_dza(i,j,k)
374 pfv(i,j,k) = (((za_bk(ib,jb)*dp_bk(ib,jb) + intp_dza(i,j,k)) - &
375 (za_bk(ib,jb+1)*dp_bk(ib,jb+1) + intp_dza(i,j+1,k))) + &
376 ((dp_bk(ib,jb+1) - dp_bk(ib,jb)) * inty_za_bk(ib,jb) - &
377 (p(i,j+1,k) - p(i,j,k)) * inty_dza(i,j,k))) * &
378 (2.0*g%IdyCv(i,j) / ((dp_bk(ib,jb) + dp_bk(ib,jb+1)) + &
382 if (cs%GFS_scale < 1.0)
then
384 do j=js_bk+joff_bk,je_bk+joff_bk ;
do i=isq_bk+ioff_bk,ieq_bk+ioff_bk
385 pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
387 do j=jsq_bk+joff_bk,jeq_bk+joff_bk ;
do i=is_bk+ioff_bk,ie_bk+ioff_bk
388 pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
394 if (
present(pbce))
then
395 call set_pbce_nonbouss(p, tv_tmp, g, gv, us, cs%GFS_scale, pbce)
398 if (
present(eta))
then
399 pa_to_h = 1.0 / gv%H_to_Pa
402 do j=jsq,jeq+1 ;
do i=isq,ieq+1
403 eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*pa_to_h
407 do j=jsq,jeq+1 ;
do i=isq,ieq+1
408 eta(i,j) = p(i,j,nz+1)*pa_to_h
413 if (cs%id_e_tidal>0)
call post_data(cs%id_e_tidal, e_tidal, cs%diag)
415 end subroutine pressureforce_blk_afv_nonbouss
426 subroutine pressureforce_blk_afv_bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
430 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
432 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: pfu
433 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: pfv
435 type(
ale_cs),
pointer :: ale_csp
436 real,
dimension(:,:),
optional,
pointer :: p_atm
438 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(out) :: pbce
441 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: eta
445 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e
446 real,
dimension(SZI_(G),SZJ_(G)) :: &
451 real,
dimension(SZI_(G)) :: &
454 real,
dimension(SZDI_(G%Block(1)),SZDJ_(G%Block(1))) :: &
462 real,
dimension(SZDIB_(G%Block(1)),SZDJ_(G%Block(1))) :: &
466 real,
dimension(SZDI_(G%Block(1)),SZDJB_(G%Block(1))) :: &
471 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
target :: &
476 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
479 real :: rho_in_situ(szi_(g))
480 real :: p_ref(szi_(g))
496 real,
parameter :: c1_6 = 1.0/6.0
497 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
498 integer :: is_bk, ie_bk, js_bk, je_bk, isq_bk, ieq_bk, jsq_bk, jeq_bk
499 integer :: ioff_bk, joff_bk
500 integer :: i, j, k, n, ib, jb
502 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
503 nkmb=gv%nk_rho_varies
504 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
506 if (.not.
associated(cs))
call mom_error(fatal, &
507 "MOM_PressureForce_AFV_Bouss: Module must be initialized before it is used.")
510 if (
present(p_atm))
then ;
if (
associated(p_atm)) use_p_atm = .true. ;
endif
511 use_eos =
associated(tv%eqn_of_state)
512 do i=isq,ieq+1 ; p0(i) = 0.0 ;
enddo
514 if (
associated(ale_csp)) use_ale = cs%reconstruct .and. use_eos
516 h_neglect = gv%H_subroundoff
517 dz_neglect = gv%H_subroundoff * gv%H_to_Z
519 g_earth_z = us%L_T_to_m_s**2 * gv%g_Earth
520 g_rho0 = g_earth_z / gv%Rho0
531 e(i,j,1) = -g%bathyT(i,j)
533 do k=1,nz ;
do i=isq,ieq+1
534 e(i,j,1) = e(i,j,1) + h(i,j,k)*gv%H_to_Z
537 call calc_tidal_forcing(cs%Time, e(:,:,1), e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
543 do j=jsq,jeq+1 ;
do i=isq,ieq+1
544 e(i,j,nz+1) = -(g%bathyT(i,j) + e_tidal(i,j))
548 do j=jsq,jeq+1 ;
do i=isq,ieq+1
549 e(i,j,nz+1) = -g%bathyT(i,j)
553 do j=jsq,jeq+1;
do k=nz,1,-1 ;
do i=isq,ieq+1
554 e(i,j,k) = e(i,j,k+1) + h(i,j,k)*gv%H_to_Z
555 enddo ;
enddo ;
enddo
564 tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
565 tv_tmp%eqn_of_state => tv%eqn_of_state
567 do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ;
enddo
570 do k=1,nkmb ;
do i=isq,ieq+1
571 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
574 rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state)
576 do k=nkmb+1,nz ;
do i=isq,ieq+1
577 if (gv%Rlay(k) < rho_cv_bl(i))
then
578 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
580 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
585 tv_tmp%T => tv%T ; tv_tmp%S => tv%S
586 tv_tmp%eqn_of_state => tv%eqn_of_state
590 if (cs%GFS_scale < 1.0)
then
597 rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state)
600 rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state)
603 dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * rho_in_situ(i)) * e(i,j,1)
608 do j=jsq,jeq+1 ;
do i=isq,ieq+1
609 dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * gv%Rlay(1)) * e(i,j,1)
620 if ( cs%Recon_Scheme == 1 )
then
621 call pressure_gradient_plm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
622 elseif ( cs%Recon_Scheme == 2 )
then
623 call pressure_gradient_ppm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
635 is_bk=g%Block(n)%isc ; ie_bk=g%Block(n)%iec
636 js_bk=g%Block(n)%jsc ; je_bk=g%Block(n)%jec
637 isq_bk=g%Block(n)%IscB ; ieq_bk=g%Block(n)%IecB
638 jsq_bk=g%Block(n)%JscB ; jeq_bk=g%Block(n)%JecB
639 ioff_bk = g%Block(n)%idg_offset - g%HI%idg_offset
640 joff_bk = g%Block(n)%jdg_offset - g%HI%jdg_offset
646 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
647 i = ib+ioff_bk ; j = jb+joff_bk
648 pa_bk(ib,jb) = (rho_ref*g_earth_z)*e(i,j,1) + p_atm(i,j)
651 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
652 i = ib+ioff_bk ; j = jb+joff_bk
653 pa_bk(ib,jb) = (rho_ref*g_earth_z)*e(i,j,1)
656 do jb=js_bk,je_bk ;
do ib=isq_bk,ieq_bk
657 intx_pa_bk(ib,jb) = 0.5*(pa_bk(ib,jb) + pa_bk(ib+1,jb))
659 do jb=jsq_bk,jeq_bk ;
do ib=is_bk,ie_bk
660 inty_pa_bk(ib,jb) = 0.5*(pa_bk(ib,jb) + pa_bk(ib,jb+1))
674 if ( cs%Recon_Scheme == 1 )
then
675 call int_density_dz_generic_plm( t_t(:,:,k), t_b(:,:,k), &
676 s_t(:,:,k), s_b(:,:,k), e(:,:,k), e(:,:,k+1), &
677 rho_ref, cs%Rho0, g_earth_z, &
678 dz_neglect, g%bathyT, g%HI, g%Block(n), &
679 tv%eqn_of_state, dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, &
680 usemasswghtinterp = cs%useMassWghtInterp)
681 elseif ( cs%Recon_Scheme == 2 )
then
682 call int_density_dz_generic_ppm( tv%T(:,:,k), t_t(:,:,k), t_b(:,:,k), &
683 tv%S(:,:,k), s_t(:,:,k), s_b(:,:,k), e(:,:,k), e(:,:,k+1), &
684 rho_ref, cs%Rho0, g_earth_z, &
685 g%HI, g%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, &
686 intx_dpa_bk, inty_dpa_bk)
689 call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,k), e(:,:,k+1), &
690 rho_ref, cs%Rho0, g_earth_z, g%HI, g%Block(n), tv%eqn_of_state, &
691 dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, &
692 g%bathyT, dz_neglect, cs%useMassWghtInterp)
694 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
695 intz_dpa_bk(ib,jb) = intz_dpa_bk(ib,jb)*gv%Z_to_H
698 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
699 i = ib+ioff_bk ; j = jb+joff_bk
700 dz_bk(ib,jb) = g_earth_z*gv%H_to_Z*h(i,j,k)
701 dpa_bk(ib,jb) = (gv%Rlay(k) - rho_ref)*dz_bk(ib,jb)
702 intz_dpa_bk(ib,jb) = 0.5*(gv%Rlay(k) - rho_ref)*dz_bk(ib,jb)*h(i,j,k)
704 do jb=js_bk,je_bk ;
do ib=isq_bk,ieq_bk
705 intx_dpa_bk(ib,jb) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_bk(ib,jb)+dz_bk(ib+1,jb))
707 do jb=jsq_bk,jeq_bk ;
do ib=is_bk,ie_bk
708 inty_dpa_bk(ib,jb) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_bk(ib,jb)+dz_bk(ib,jb+1))
713 do jb=js_bk,je_bk ;
do ib=isq_bk,ieq_bk
714 i = ib+ioff_bk ; j = jb+joff_bk
715 pfu(i,j,k) = (((pa_bk(ib,jb)*h(i,j,k) + intz_dpa_bk(ib,jb)) - &
716 (pa_bk(ib+1,jb)*h(i+1,j,k) + intz_dpa_bk(ib+1,jb))) + &
717 ((h(i+1,j,k) - h(i,j,k)) * intx_pa_bk(ib,jb) - &
718 (e(i+1,j,k+1) - e(i,j,k+1)) * intx_dpa_bk(ib,jb) * gv%Z_to_H)) * &
719 ((2.0*i_rho0*g%IdxCu(i,j)) / &
720 ((h(i,j,k) + h(i+1,j,k)) + h_neglect))
721 intx_pa_bk(ib,jb) = intx_pa_bk(ib,jb) + intx_dpa_bk(ib,jb)
724 do jb=jsq_bk,jeq_bk ;
do ib=is_bk,ie_bk
725 i = ib+ioff_bk ; j = jb+joff_bk
726 pfv(i,j,k) = (((pa_bk(ib,jb)*h(i,j,k) + intz_dpa_bk(ib,jb)) - &
727 (pa_bk(ib,jb+1)*h(i,j+1,k) + intz_dpa_bk(ib,jb+1))) + &
728 ((h(i,j+1,k) - h(i,j,k)) * inty_pa_bk(ib,jb) - &
729 (e(i,j+1,k+1) - e(i,j,k+1)) * inty_dpa_bk(ib,jb) * gv%Z_to_H)) * &
730 ((2.0*i_rho0*g%IdyCv(i,j)) / &
731 ((h(i,j,k) + h(i,j+1,k)) + h_neglect))
732 inty_pa_bk(ib,jb) = inty_pa_bk(ib,jb) + inty_dpa_bk(ib,jb)
734 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
735 pa_bk(ib,jb) = pa_bk(ib,jb) + dpa_bk(ib,jb)
739 if (cs%GFS_scale < 1.0)
then
741 do j=js_bk+joff_bk,je_bk+joff_bk ;
do i=isq_bk+ioff_bk,ieq_bk+ioff_bk
742 pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
744 do j=jsq_bk+joff_bk,jeq_bk+joff_bk ;
do i=is_bk+ioff_bk,ie_bk+ioff_bk
745 pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
751 if (
present(pbce))
then
752 call set_pbce_bouss(e, tv_tmp, g, gv, us, cs%Rho0, cs%GFS_scale, pbce)
755 if (
present(eta))
then
761 do j=jsq,jeq+1 ;
do i=isq,ieq+1
762 eta(i,j) = e(i,j,1)*gv%Z_to_H + e_tidal(i,j)*gv%Z_to_H
766 do j=jsq,jeq+1 ;
do i=isq,ieq+1
767 eta(i,j) = e(i,j,1)*gv%Z_to_H
772 if (cs%id_e_tidal>0)
call post_data(cs%id_e_tidal, e_tidal, cs%diag)
774 end subroutine pressureforce_blk_afv_bouss
777 subroutine pressureforce_blk_afv_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)
778 type(time_type),
target,
intent(in) :: time
783 type(
diag_ctrl),
target,
intent(inout) :: diag
787 #include "version_variable.h"
788 character(len=40) :: mdl
791 if (
associated(cs))
then
792 call mom_error(warning,
"PressureForce_init called with an associated "// &
793 "control structure.")
795 else ;
allocate(cs) ;
endif
797 cs%diag => diag ; cs%Time => time
798 if (
present(tides_csp))
then
799 if (
associated(tides_csp)) cs%tides_CSp => tides_csp
802 mdl =
"MOM_PressureForce_blk_AFV"
804 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
805 "The mean ocean density used with BOUSSINESQ true to "//&
806 "calculate accelerations and the mass for conservation "//&
807 "properties, or with BOUSSINSEQ false to convert some "//&
808 "parameters from vertical units of m to kg m-2.", &
809 units=
"kg m-3", default=1035.0)
810 call get_param(param_file, mdl,
"TIDES", cs%tides, &
811 "If true, apply tidal momentum forcing.", default=.false.)
812 call get_param(param_file,
"MOM",
"USE_REGRIDDING", use_ale, &
813 "If True, use the ALE algorithm (regridding/remapping). "//&
814 "If False, use the layered isopycnal algorithm.", default=.false. )
815 call get_param(param_file, mdl,
"MASS_WEIGHT_IN_PRESSURE_GRADIENT", cs%useMassWghtInterp, &
816 "If true, use mass weighting when interpolating T/S for "//&
817 "integrals near the bathymetry in AFV pressure gradient "//&
818 "calculations.", default=.false.)
819 call get_param(param_file, mdl,
"RECONSTRUCT_FOR_PRESSURE", cs%reconstruct, &
820 "If True, use vertical reconstruction of T & S within "//&
821 "the integrals of the FV pressure gradient calculation. "//&
822 "If False, use the constant-by-layer algorithm. "//&
823 "The default is set by USE_REGRIDDING.", &
825 call get_param(param_file, mdl,
"PRESSURE_RECONSTRUCTION_SCHEME", cs%Recon_Scheme, &
826 "Order of vertical reconstruction of T/S to use in the "//&
827 "integrals within the FV pressure gradient calculation.\n"//&
828 " 0: PCM or no reconstruction.\n"//&
829 " 1: PLM reconstruction.\n"//&
830 " 2: PPM reconstruction.", default=1)
831 call get_param(param_file, mdl,
"BOUNDARY_EXTRAPOLATION_PRESSURE", cs%boundary_extrap, &
832 "If true, the reconstruction of T & S for pressure in "//&
833 "boundary cells is extrapolated, rather than using PCM "//&
834 "in these cells. If true, the same order polynomial is "//&
835 "used as is used for the interior cells.", default=.true.)
838 cs%id_e_tidal = register_diag_field(
'ocean_model',
'e_tidal', diag%axesT1, &
839 time,
'Tidal Forcing Astronomical and SAL Height Anomaly',
'meter', conversion=us%Z_to_m)
843 if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
845 call log_param(param_file, mdl,
"GFS / G_EARTH", cs%GFS_scale)
847 end subroutine pressureforce_blk_afv_init
850 subroutine pressureforce_blk_afv_end(CS)
853 if (
associated(cs))
deallocate(cs)
854 end subroutine pressureforce_blk_afv_end