23 implicit none ;
private
25 #include <MOM_memory.h>
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)
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, scale=us%kg_m3_to_R)
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/(us%R_to_kg_m3*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
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 us%m_s_to_L_T**2*(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 us%m_s_to_L_T**2*(p(i,j,1)*(1.0/(us%R_to_kg_m3*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 (us%m_s_to_L_T**2 * 2.0*g%IdxCu(i,j) / &
369 ((dp_bk(ib,jb) + dp_bk(ib+1,jb)) + dp_neglect))
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 (us%m_s_to_L_T**2 * 2.0*g%IdyCv(i,j) / &
379 ((dp_bk(ib,jb) + dp_bk(ib,jb+1)) + dp_neglect))
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
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)
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))) :: &
463 real,
dimension(SZDIB_(G%Block(1)),SZDJ_(G%Block(1))) :: &
467 real,
dimension(SZDI_(G%Block(1)),SZDJB_(G%Block(1))) :: &
472 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
target :: &
477 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
480 real :: rho_in_situ(szi_(g))
481 real :: p_ref(szi_(g))
487 real :: g_earth_mks_z
488 real :: g_earth_z_geo
499 real,
parameter :: c1_6 = 1.0/6.0
500 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
501 integer :: is_bk, ie_bk, js_bk, je_bk, isq_bk, ieq_bk, jsq_bk, jeq_bk
502 integer :: ioff_bk, joff_bk
503 integer :: i, j, k, n, ib, jb
505 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
506 nkmb=gv%nk_rho_varies
507 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
509 if (.not.
associated(cs))
call mom_error(fatal, &
510 "MOM_PressureForce_AFV_Bouss: Module must be initialized before it is used.")
513 if (
present(p_atm))
then ;
if (
associated(p_atm)) use_p_atm = .true. ;
endif
514 use_eos =
associated(tv%eqn_of_state)
515 do i=isq,ieq+1 ; p0(i) = 0.0 ;
enddo
517 if (
associated(ale_csp)) use_ale = cs%reconstruct .and. use_eos
519 h_neglect = gv%H_subroundoff
520 dz_neglect = gv%H_subroundoff * gv%H_to_Z
521 i_rho0 = us%m_s_to_L_T**2 / (us%R_to_kg_m3*gv%Rho0)
522 g_earth_mks_z = us%L_T_to_m_s**2 * gv%g_Earth
523 g_earth_z_geo = us%R_to_kg_m3*us%L_T_to_m_s**2 * gv%g_Earth
524 g_rho0 = gv%g_Earth / gv%Rho0
525 rho_ref_mks = cs%Rho0
526 rho_ref = rho_ref_mks*us%kg_m3_to_R
536 e(i,j,1) = -g%bathyT(i,j)
538 do k=1,nz ;
do i=isq,ieq+1
539 e(i,j,1) = e(i,j,1) + h(i,j,k)*gv%H_to_Z
542 call calc_tidal_forcing(cs%Time, e(:,:,1), e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
548 do j=jsq,jeq+1 ;
do i=isq,ieq+1
549 e(i,j,nz+1) = -(g%bathyT(i,j) + e_tidal(i,j))
553 do j=jsq,jeq+1 ;
do i=isq,ieq+1
554 e(i,j,nz+1) = -g%bathyT(i,j)
558 do j=jsq,jeq+1;
do k=nz,1,-1 ;
do i=isq,ieq+1
559 e(i,j,k) = e(i,j,k+1) + h(i,j,k)*gv%H_to_Z
560 enddo ;
enddo ;
enddo
569 tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
570 tv_tmp%eqn_of_state => tv%eqn_of_state
572 do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ;
enddo
575 do k=1,nkmb ;
do i=isq,ieq+1
576 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
579 rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
581 do k=nkmb+1,nz ;
do i=isq,ieq+1
582 if (gv%Rlay(k) < rho_cv_bl(i))
then
583 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
585 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
590 tv_tmp%T => tv%T ; tv_tmp%S => tv%S
591 tv_tmp%eqn_of_state => tv%eqn_of_state
595 if (cs%GFS_scale < 1.0)
then
602 rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
605 rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
608 dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * rho_in_situ(i)) * e(i,j,1)
613 do j=jsq,jeq+1 ;
do i=isq,ieq+1
614 dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * gv%Rlay(1)) * e(i,j,1)
625 if ( cs%Recon_Scheme == 1 )
then
627 elseif ( cs%Recon_Scheme == 2 )
then
640 is_bk=g%Block(n)%isc ; ie_bk=g%Block(n)%iec
641 js_bk=g%Block(n)%jsc ; je_bk=g%Block(n)%jec
642 isq_bk=g%Block(n)%IscB ; ieq_bk=g%Block(n)%IecB
643 jsq_bk=g%Block(n)%JscB ; jeq_bk=g%Block(n)%JecB
644 ioff_bk = g%Block(n)%idg_offset - g%HI%idg_offset
645 joff_bk = g%Block(n)%jdg_offset - g%HI%jdg_offset
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_geo)*e(i,j,1) + p_atm(i,j)
656 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
657 i = ib+ioff_bk ; j = jb+joff_bk
658 pa_bk(ib,jb) = (rho_ref*g_earth_z_geo)*e(i,j,1)
661 do jb=js_bk,je_bk ;
do ib=isq_bk,ieq_bk
662 intx_pa_bk(ib,jb) = 0.5*(pa_bk(ib,jb) + pa_bk(ib+1,jb))
664 do jb=jsq_bk,jeq_bk ;
do ib=is_bk,ie_bk
665 inty_pa_bk(ib,jb) = 0.5*(pa_bk(ib,jb) + pa_bk(ib,jb+1))
679 if ( cs%Recon_Scheme == 1 )
then
680 call int_density_dz_generic_plm( t_t(:,:,k), t_b(:,:,k), &
681 s_t(:,:,k), s_b(:,:,k), e(:,:,k), e(:,:,k+1), &
682 rho_ref_mks, cs%Rho0, g_earth_mks_z, &
683 dz_neglect, g%bathyT, g%HI, g%Block(n), &
684 tv%eqn_of_state, dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, &
685 usemasswghtinterp = cs%useMassWghtInterp)
686 elseif ( cs%Recon_Scheme == 2 )
then
687 call int_density_dz_generic_ppm( tv%T(:,:,k), t_t(:,:,k), t_b(:,:,k), &
688 tv%S(:,:,k), s_t(:,:,k), s_b(:,:,k), e(:,:,k), e(:,:,k+1), &
689 rho_ref_mks, cs%Rho0, g_earth_mks_z, &
690 g%HI, g%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, &
691 intx_dpa_bk, inty_dpa_bk)
694 call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,k), e(:,:,k+1), &
695 rho_ref_mks, cs%Rho0, g_earth_mks_z, g%HI, g%Block(n), tv%eqn_of_state, &
696 dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, &
697 g%bathyT, dz_neglect, cs%useMassWghtInterp)
699 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
700 intz_dpa_bk(ib,jb) = intz_dpa_bk(ib,jb)*gv%Z_to_H
703 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
704 i = ib+ioff_bk ; j = jb+joff_bk
705 dz_bk(ib,jb) = g_earth_z_geo*gv%H_to_Z*h(i,j,k)
706 dpa_bk(ib,jb) = (gv%Rlay(k) - rho_ref)*dz_bk(ib,jb)
707 intz_dpa_bk(ib,jb) = 0.5*(gv%Rlay(k) - rho_ref) * dz_bk(ib,jb)*h(i,j,k)
709 do jb=js_bk,je_bk ;
do ib=isq_bk,ieq_bk
710 intx_dpa_bk(ib,jb) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_bk(ib,jb)+dz_bk(ib+1,jb))
712 do jb=jsq_bk,jeq_bk ;
do ib=is_bk,ie_bk
713 inty_dpa_bk(ib,jb) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_bk(ib,jb)+dz_bk(ib,jb+1))
718 do jb=js_bk,je_bk ;
do ib=isq_bk,ieq_bk
719 i = ib+ioff_bk ; j = jb+joff_bk
720 pfu(i,j,k) = (((pa_bk(ib,jb)*h(i,j,k) + intz_dpa_bk(ib,jb)) - &
721 (pa_bk(ib+1,jb)*h(i+1,j,k) + intz_dpa_bk(ib+1,jb))) + &
722 ((h(i+1,j,k) - h(i,j,k)) * intx_pa_bk(ib,jb) - &
723 (e(i+1,j,k+1) - e(i,j,k+1)) * intx_dpa_bk(ib,jb) * gv%Z_to_H)) * &
724 ((2.0*i_rho0*g%IdxCu(i,j)) / &
725 ((h(i,j,k) + h(i+1,j,k)) + h_neglect))
726 intx_pa_bk(ib,jb) = intx_pa_bk(ib,jb) + intx_dpa_bk(ib,jb)
729 do jb=jsq_bk,jeq_bk ;
do ib=is_bk,ie_bk
730 i = ib+ioff_bk ; j = jb+joff_bk
731 pfv(i,j,k) = (((pa_bk(ib,jb)*h(i,j,k) + intz_dpa_bk(ib,jb)) - &
732 (pa_bk(ib,jb+1)*h(i,j+1,k) + intz_dpa_bk(ib,jb+1))) + &
733 ((h(i,j+1,k) - h(i,j,k)) * inty_pa_bk(ib,jb) - &
734 (e(i,j+1,k+1) - e(i,j,k+1)) * inty_dpa_bk(ib,jb) * gv%Z_to_H)) * &
735 ((2.0*i_rho0*g%IdyCv(i,j)) / &
736 ((h(i,j,k) + h(i,j+1,k)) + h_neglect))
737 inty_pa_bk(ib,jb) = inty_pa_bk(ib,jb) + inty_dpa_bk(ib,jb)
739 do jb=jsq_bk,jeq_bk+1 ;
do ib=isq_bk,ieq_bk+1
740 pa_bk(ib,jb) = pa_bk(ib,jb) + dpa_bk(ib,jb)
744 if (cs%GFS_scale < 1.0)
then
746 do j=js_bk+joff_bk,je_bk+joff_bk ;
do i=isq_bk+ioff_bk,ieq_bk+ioff_bk
747 pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
749 do j=jsq_bk+joff_bk,jeq_bk+joff_bk ;
do i=is_bk+ioff_bk,ie_bk+ioff_bk
750 pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
756 if (
present(pbce))
then
757 call set_pbce_bouss(e, tv_tmp, g, gv, us, cs%Rho0, cs%GFS_scale, pbce)
760 if (
present(eta))
then
766 do j=jsq,jeq+1 ;
do i=isq,ieq+1
767 eta(i,j) = e(i,j,1)*gv%Z_to_H + e_tidal(i,j)*gv%Z_to_H
771 do j=jsq,jeq+1 ;
do i=isq,ieq+1
772 eta(i,j) = e(i,j,1)*gv%Z_to_H
777 if (cs%id_e_tidal>0)
call post_data(cs%id_e_tidal, e_tidal, cs%diag)
783 type(time_type),
target,
intent(in) :: time
788 type(
diag_ctrl),
target,
intent(inout) :: diag
792 #include "version_variable.h"
793 character(len=40) :: mdl
796 if (
associated(cs))
then
797 call mom_error(warning,
"PressureForce_init called with an associated "// &
798 "control structure.")
800 else ;
allocate(cs) ;
endif
802 cs%diag => diag ; cs%Time => time
803 if (
present(tides_csp))
then
804 if (
associated(tides_csp)) cs%tides_CSp => tides_csp
807 mdl =
"MOM_PressureForce_blk_AFV"
809 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
810 "The mean ocean density used with BOUSSINESQ true to "//&
811 "calculate accelerations and the mass for conservation "//&
812 "properties, or with BOUSSINSEQ false to convert some "//&
813 "parameters from vertical units of m to kg m-2.", &
814 units=
"kg m-3", default=1035.0)
815 call get_param(param_file, mdl,
"TIDES", cs%tides, &
816 "If true, apply tidal momentum forcing.", default=.false.)
817 call get_param(param_file,
"MOM",
"USE_REGRIDDING", use_ale, &
818 "If True, use the ALE algorithm (regridding/remapping). "//&
819 "If False, use the layered isopycnal algorithm.", default=.false. )
820 call get_param(param_file, mdl,
"MASS_WEIGHT_IN_PRESSURE_GRADIENT", cs%useMassWghtInterp, &
821 "If true, use mass weighting when interpolating T/S for "//&
822 "integrals near the bathymetry in AFV pressure gradient "//&
823 "calculations.", default=.false.)
824 call get_param(param_file, mdl,
"RECONSTRUCT_FOR_PRESSURE", cs%reconstruct, &
825 "If True, use vertical reconstruction of T & S within "//&
826 "the integrals of the FV pressure gradient calculation. "//&
827 "If False, use the constant-by-layer algorithm. "//&
828 "The default is set by USE_REGRIDDING.", &
830 call get_param(param_file, mdl,
"PRESSURE_RECONSTRUCTION_SCHEME", cs%Recon_Scheme, &
831 "Order of vertical reconstruction of T/S to use in the "//&
832 "integrals within the FV pressure gradient calculation.\n"//&
833 " 0: PCM or no reconstruction.\n"//&
834 " 1: PLM reconstruction.\n"//&
835 " 2: PPM reconstruction.", default=1)
836 call get_param(param_file, mdl,
"BOUNDARY_EXTRAPOLATION_PRESSURE", cs%boundary_extrap, &
837 "If true, the reconstruction of T & S for pressure in "//&
838 "boundary cells is extrapolated, rather than using PCM "//&
839 "in these cells. If true, the same order polynomial is "//&
840 "used as is used for the interior cells.", default=.true.)
843 cs%id_e_tidal = register_diag_field(
'ocean_model',
'e_tidal', diag%axesT1, &
844 time,
'Tidal Forcing Astronomical and SAL Height Anomaly',
'meter', conversion=us%Z_to_m)
848 if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
850 call log_param(param_file, mdl,
"GFS / G_EARTH", cs%GFS_scale)
858 if (
associated(cs))
deallocate(cs)