16 use mom_eos,
only : int_specific_vol_dp, query_compressible
18 implicit none ;
private
20 #include <MOM_memory.h>
22 public pressureforce_mont_bouss, pressureforce_mont_nonbouss, set_pbce_bouss
23 public set_pbce_nonbouss, pressureforce_mont_init, pressureforce_mont_end
40 type(time_type),
pointer :: time => null()
43 real,
pointer :: pfu_bc(:,:,:) => null()
44 real,
pointer :: pfv_bc(:,:,:) => null()
47 integer :: id_pfu_bc = -1, id_pfv_bc = -1, id_e_tidal = -1
63 subroutine pressureforce_mont_nonbouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
67 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
69 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: pfu
71 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: pfv
74 real,
dimension(:,:),
optional,
pointer :: p_atm
76 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
77 optional,
intent(out) :: pbce
80 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: eta
83 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
87 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p
91 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
target :: &
97 real,
dimension(SZI_(G)) :: rho_cv_bl
100 real,
dimension(SZI_(G),SZJ_(G)) :: &
110 real :: p_ref(szi_(g))
112 real :: rho_in_situ(szi_(g))
113 real :: pfu_bc, pfv_bc
128 real :: alpha_lay(szk_(g))
129 real :: dalpha_int(szk_(g)+1)
131 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
133 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
134 nkmb=gv%nk_rho_varies
135 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
138 if (
present(p_atm))
then ;
if (
associated(p_atm)) use_p_atm = .true. ;
endif
139 is_split = .false. ;
if (
present(pbce)) is_split = .true.
140 use_eos =
associated(tv%eqn_of_state)
142 if (.not.
associated(cs))
call mom_error(fatal, &
143 "MOM_PressureForce_Mont: Module must be initialized before it is used.")
145 if (query_compressible(tv%eqn_of_state))
call mom_error(fatal, &
146 "PressureForce_Mont_nonBouss: The Montgomery form of the pressure force "//&
147 "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.")
150 i_gearth = 1.0 / (us%L_T_to_m_s**2 * gv%g_Earth)
151 dp_neglect = gv%H_to_Pa * gv%H_subroundoff
152 do k=1,nz ; alpha_lay(k) = 1.0 / gv%Rlay(k) ;
enddo
153 do k=2,nz ; dalpha_int(k) = alpha_lay(k-1) - alpha_lay(k) ;
enddo
157 do j=jsq,jeq+1 ;
do i=isq,ieq+1 ; p(i,j,1) = p_atm(i,j) ;
enddo ;
enddo
160 do j=jsq,jeq+1 ;
do i=isq,ieq+1 ; p(i,j,1) = 0.0 ;
enddo ;
enddo
163 do j=jsq,jeq+1 ;
do k=1,nz ;
do i=isq,ieq+1
164 p(i,j,k+1) = p(i,j,k) + gv%H_to_Pa * h(i,j,k)
165 enddo ;
enddo ;
enddo
167 if (
present(eta))
then
168 pa_to_h = 1.0 / gv%H_to_Pa
171 do j=jsq,jeq+1 ;
do i=isq,ieq+1
172 eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*pa_to_h
176 do j=jsq,jeq+1 ;
do i=isq,ieq+1
177 eta(i,j) = p(i,j,nz+1)*pa_to_h
186 do j=jsq,jeq+1 ;
do i=isq,ieq+1
187 ssh(i,j) = -g%bathyT(i,j)
192 call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), &
193 0.0, g%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1)
196 do j=jsq,jeq+1 ;
do k=1,nz;
do i=isq,ieq+1
197 ssh(i,j) = ssh(i,j) + i_gearth * dz_geo(i,j,k)
198 enddo ;
enddo ;
enddo
201 do j=jsq,jeq+1 ;
do k=1,nz ;
do i=isq,ieq+1
202 ssh(i,j) = ssh(i,j) + (us%m_to_Z*gv%H_to_kg_m2)*h(i,j,k)*alpha_lay(k)
203 enddo ;
enddo ;
enddo
206 call calc_tidal_forcing(cs%Time, ssh, e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
208 do j=jsq,jeq+1 ;
do i=isq,ieq+1
209 geopot_bot(i,j) = -us%L_T_to_m_s**2 * gv%g_Earth*(e_tidal(i,j) + g%bathyT(i,j))
213 do j=jsq,jeq+1 ;
do i=isq,ieq+1
214 geopot_bot(i,j) = -us%L_T_to_m_s**2 * gv%g_Earth*g%bathyT(i,j)
226 tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
227 tv_tmp%eqn_of_state => tv%eqn_of_state
228 do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ;
enddo
231 do k=1,nkmb ;
do i=isq,ieq+1
232 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
235 rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state)
236 do k=nkmb+1,nz ;
do i=isq,ieq+1
237 if (gv%Rlay(k) < rho_cv_bl(i))
then
238 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
240 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
245 tv_tmp%T => tv%T ; tv_tmp%S => tv%S
246 tv_tmp%eqn_of_state => tv%eqn_of_state
247 do i=isq,ieq+1 ; p_ref(i) = 0 ;
enddo
250 do k=1,nz ;
do j=jsq,jeq+1
252 rho_in_situ,isq,ieq-isq+2,tv%eqn_of_state)
253 do i=isq,ieq+1 ; alpha_star(i,j,k) = 1.0 / rho_in_situ(i) ;
enddo
261 m(i,j,nz) = geopot_bot(i,j) + p(i,j,nz+1) * alpha_star(i,j,nz)
263 do k=nz-1,1,-1 ;
do i=isq,ieq+1
264 m(i,j,k) = m(i,j,k+1) + p(i,j,k+1) * (alpha_star(i,j,k) - alpha_star(i,j,k+1))
271 m(i,j,nz) = geopot_bot(i,j) + p(i,j,nz+1) * alpha_lay(nz)
273 do k=nz-1,1,-1 ;
do i=isq,ieq+1
274 m(i,j,k) = m(i,j,k+1) + p(i,j,k+1) * dalpha_int(k+1)
279 if (cs%GFS_scale < 1.0)
then
282 do j=jsq,jeq+1 ;
do i=isq,ieq+1
283 dm(i,j) = (cs%GFS_scale - 1.0) * m(i,j,1)
286 do k=1,nz ;
do j=jsq,jeq+1 ;
do i=isq,ieq+1
287 m(i,j,k) = m(i,j,k) + dm(i,j)
288 enddo ;
enddo ;
enddo
308 if (
present(pbce))
then
309 call set_pbce_nonbouss(p, tv_tmp, g, gv, us, cs%GFS_scale, pbce, alpha_star)
317 do j=jsq,jeq+1 ;
do i=isq,ieq+1
318 dp_star(i,j) = (p(i,j,k+1) - p(i,j,k)) + dp_neglect
320 do j=js,je ;
do i=isq,ieq
322 pfu_bc = (alpha_star(i+1,j,k) - alpha_star(i,j,k)) * (g%IdxCu(i,j) * &
323 ((dp_star(i,j) * dp_star(i+1,j) + (p(i,j,k) * dp_star(i+1,j) + &
324 p(i+1,j,k) * dp_star(i,j))) / (dp_star(i,j) + dp_star(i+1,j))))
325 pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j) + pfu_bc
326 if (
associated(cs%PFu_bc)) cs%PFu_bc(i,j,k) = pfu_bc
328 do j=jsq,jeq ;
do i=is,ie
329 pfv_bc = (alpha_star(i,j+1,k) - alpha_star(i,j,k)) * (g%IdyCv(i,j) * &
330 ((dp_star(i,j) * dp_star(i,j+1) + (p(i,j,k) * dp_star(i,j+1) + &
331 p(i,j+1,k) * dp_star(i,j))) / (dp_star(i,j) + dp_star(i,j+1))))
332 pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j) + pfv_bc
333 if (
associated(cs%PFv_bc)) cs%PFv_bc(i,j,k) = pfv_bc
339 do j=js,je ;
do i=isq,ieq
340 pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j)
342 do j=jsq,jeq ;
do i=is,ie
343 pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j)
348 if (cs%id_PFu_bc>0)
call post_data(cs%id_PFu_bc, cs%PFu_bc, cs%diag)
349 if (cs%id_PFv_bc>0)
call post_data(cs%id_PFv_bc, cs%PFv_bc, cs%diag)
350 if (cs%id_e_tidal>0)
call post_data(cs%id_e_tidal, e_tidal, cs%diag)
352 end subroutine pressureforce_mont_nonbouss
361 subroutine pressureforce_mont_bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
365 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
367 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: pfu
369 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: pfv
372 real,
dimension(:,:),
optional,
pointer :: p_atm
374 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(out) :: pbce
377 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: eta
379 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
383 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e
387 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
target :: &
393 real :: rho_cv_bl(szi_(g))
395 real :: h_star(szi_(g),szj_(g))
397 real :: e_tidal(szi_(g),szj_(g))
400 real :: p_ref(szi_(g))
404 real :: pfu_bc, pfv_bc
416 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
419 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
420 nkmb=gv%nk_rho_varies
421 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
424 if (
present(p_atm))
then ;
if (
associated(p_atm)) use_p_atm = .true. ;
endif
425 is_split = .false. ;
if (
present(pbce)) is_split = .true.
426 use_eos =
associated(tv%eqn_of_state)
428 if (.not.
associated(cs))
call mom_error(fatal, &
429 "MOM_PressureForce_Mont: Module must be initialized before it is used.")
431 if (query_compressible(tv%eqn_of_state))
call mom_error(fatal, &
432 "PressureForce_Mont_Bouss: The Montgomery form of the pressure force "//&
433 "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.")
436 h_neglect = gv%H_subroundoff * gv%H_to_Z
438 g_rho0 = us%L_T_to_m_s**2 * gv%g_Earth/gv%Rho0
447 do i=isq,ieq+1 ; e(i,j,1) = -g%bathyT(i,j) ;
enddo
448 do k=1,nz ;
do i=isq,ieq+1
449 e(i,j,1) = e(i,j,1) + h(i,j,k)*gv%H_to_Z
452 call calc_tidal_forcing(cs%Time, e(:,:,1), e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
458 do j=jsq,jeq+1 ;
do i=isq,ieq+1
459 e(i,j,nz+1) = -(g%bathyT(i,j) + e_tidal(i,j))
463 do j=jsq,jeq+1 ;
do i=isq,ieq+1
464 e(i,j,nz+1) = -g%bathyT(i,j)
468 do j=jsq,jeq+1 ;
do k=nz,1,-1 ;
do i=isq,ieq+1
469 e(i,j,k) = e(i,j,k+1) + h(i,j,k)*gv%H_to_Z
470 enddo ;
enddo ;
enddo
481 tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
482 tv_tmp%eqn_of_state => tv%eqn_of_state
484 do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ;
enddo
487 do k=1,nkmb ;
do i=isq,ieq+1
488 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
491 rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state)
493 do k=nkmb+1,nz ;
do i=isq,ieq+1
494 if (gv%Rlay(k) < rho_cv_bl(i))
then
495 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
497 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
502 tv_tmp%T => tv%T ; tv_tmp%S => tv%S
503 tv_tmp%eqn_of_state => tv%eqn_of_state
504 do i=isq,ieq+1 ; p_ref(i) = 0.0 ;
enddo
510 do k=1,nz+1 ;
do j=jsq,jeq+1
511 call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), &
512 isq,ieq-isq+2,tv%eqn_of_state)
513 do i=isq,ieq+1 ; rho_star(i,j,k) = g_rho0*rho_star(i,j,k) ;
enddo
522 m(i,j,1) = cs%GFS_scale * (rho_star(i,j,1) * e(i,j,1))
523 if (use_p_atm) m(i,j,1) = m(i,j,1) + p_atm(i,j) * i_rho0
525 do k=2,nz ;
do i=isq,ieq+1
526 m(i,j,k) = m(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) * e(i,j,k)
533 m(i,j,1) = us%L_to_m**2*us%s_to_T**2*gv%g_prime(1) * e(i,j,1)
534 if (use_p_atm) m(i,j,1) = m(i,j,1) + p_atm(i,j) * i_rho0
536 do k=2,nz ;
do i=isq,ieq+1
537 m(i,j,k) = m(i,j,k-1) + us%L_to_m**2*us%s_to_T**2*gv%g_prime(k) * e(i,j,k)
542 if (
present(pbce))
then
543 call set_pbce_bouss(e, tv_tmp, g, gv, us, cs%Rho0, cs%GFS_scale, pbce, rho_star)
551 do j=jsq,jeq+1 ;
do i=isq,ieq+1
552 h_star(i,j) = (e(i,j,k) - e(i,j,k+1)) + h_neglect
554 do j=js,je ;
do i=isq,ieq
555 pfu_bc = -1.0*(rho_star(i+1,j,k) - rho_star(i,j,k)) * (g%IdxCu(i,j) * &
556 ((h_star(i,j) * h_star(i+1,j) - (e(i,j,k) * h_star(i+1,j) + &
557 e(i+1,j,k) * h_star(i,j))) / (h_star(i,j) + h_star(i+1,j))))
558 pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j) + pfu_bc
559 if (
associated(cs%PFu_bc)) cs%PFu_bc(i,j,k) = pfu_bc
561 do j=jsq,jeq ;
do i=is,ie
562 pfv_bc = -1.0*(rho_star(i,j+1,k) - rho_star(i,j,k)) * (g%IdyCv(i,j) * &
563 ((h_star(i,j) * h_star(i,j+1) - (e(i,j,k) * h_star(i,j+1) + &
564 e(i,j+1,k) * h_star(i,j))) / (h_star(i,j) + h_star(i,j+1))))
565 pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j) + pfv_bc
566 if (
associated(cs%PFv_bc)) cs%PFv_bc(i,j,k) = pfv_bc
572 do j=js,je ;
do i=isq,ieq
573 pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j)
575 do j=jsq,jeq ;
do i=is,ie
576 pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j)
581 if (
present(eta))
then
587 do j=jsq,jeq+1 ;
do i=isq,ieq+1
588 eta(i,j) = e(i,j,1)*gv%Z_to_H + e_tidal(i,j)*gv%Z_to_H
592 do j=jsq,jeq+1 ;
do i=isq,ieq+1
593 eta(i,j) = e(i,j,1)*gv%Z_to_H
598 if (cs%id_PFu_bc>0)
call post_data(cs%id_PFu_bc, cs%PFu_bc, cs%diag)
599 if (cs%id_PFv_bc>0)
call post_data(cs%id_PFv_bc, cs%PFv_bc, cs%diag)
600 if (cs%id_e_tidal>0)
call post_data(cs%id_e_tidal, e_tidal, cs%diag)
602 end subroutine pressureforce_mont_bouss
606 subroutine set_pbce_bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
609 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: e
612 real,
intent(in) :: rho0
613 real,
intent(in) :: gfs_scale
616 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
620 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
621 optional,
intent(in) :: rho_star
625 real :: ihtot(szi_(g))
626 real :: press(szi_(g))
627 real :: t_int(szi_(g))
628 real :: s_int(szi_(g))
629 real :: dr_dt(szi_(g))
630 real :: dr_ds(szi_(g))
631 real :: rho_in_situ(szi_(g))
638 integer :: isq, ieq, jsq, jeq, nz, i, j, k
640 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
642 rho0xg = rho0*us%L_T_to_m_s**2 * gv%g_Earth
643 g_rho0 = gv%g_Earth / gv%Rho0
644 use_eos =
associated(tv%eqn_of_state)
645 z_neglect = gv%H_subroundoff*gv%H_to_Z
648 if (
present(rho_star))
then
652 ihtot(i) = gv%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
653 pbce(i,j,1) = gfs_scale * us%m_s_to_L_T**2*rho_star(i,j,1) * gv%H_to_Z
655 do k=2,nz ;
do i=isq,ieq+1
656 pbce(i,j,k) = pbce(i,j,k-1) + us%m_s_to_L_T**2*(rho_star(i,j,k)-rho_star(i,j,k-1)) * &
657 ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i))
664 ihtot(i) = gv%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
665 press(i) = -rho0xg*e(i,j,1)
668 isq, ieq-isq+2, tv%eqn_of_state)
670 pbce(i,j,1) = g_rho0*(gfs_scale * rho_in_situ(i)) * gv%H_to_Z
674 press(i) = -rho0xg*e(i,j,k)
675 t_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
676 s_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
679 isq, ieq-isq+2, tv%eqn_of_state)
681 pbce(i,j,k) = pbce(i,j,k-1) + g_rho0 * &
682 ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i)) * &
683 (dr_dt(i)*(tv%T(i,j,k)-tv%T(i,j,k-1)) + &
684 dr_ds(i)*(tv%S(i,j,k)-tv%S(i,j,k-1)))
693 ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
694 pbce(i,j,1) = gv%g_prime(1) * gv%H_to_Z
696 do k=2,nz ;
do i=isq,ieq+1
697 pbce(i,j,k) = pbce(i,j,k-1) + &
698 (gv%g_prime(k)*gv%H_to_Z) * ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i))
703 end subroutine set_pbce_bouss
707 subroutine set_pbce_nonbouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
710 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: p
713 real,
intent(in) :: gfs_scale
716 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: pbce
719 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(in) :: alpha_star
722 real,
dimension(SZI_(G),SZJ_(G)) :: &
726 real :: t_int(szi_(g))
727 real :: s_int(szi_(g))
728 real :: dr_dt(szi_(g))
729 real :: dr_ds(szi_(g))
730 real :: rho_in_situ(szi_(g))
731 real :: alpha_lay(szk_(g))
732 real :: dalpha_int(szk_(g)+1)
739 integer :: isq, ieq, jsq, jeq, nz, i, j, k
741 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
743 use_eos =
associated(tv%eqn_of_state)
745 dp_dh = us%m_s_to_L_T**2*gv%H_to_Pa
746 dp_neglect = gv%H_to_Pa * gv%H_subroundoff
748 do k=1,nz ; alpha_lay(k) = 1.0 / gv%Rlay(k) ;
enddo
749 do k=2,nz ; dalpha_int(k) = alpha_lay(k-1) - alpha_lay(k) ;
enddo
752 if (
present(alpha_star))
then
756 c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
757 pbce(i,j,nz) = dp_dh * alpha_star(i,j,nz)
759 do k=nz-1,1,-1 ;
do i=isq,ieq+1
760 pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1)) * c_htot(i,j)) * &
761 (alpha_star(i,j,k) - alpha_star(i,j,k+1))
768 rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state)
770 c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
771 pbce(i,j,nz) = dp_dh / rho_in_situ(i)
775 t_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1))
776 s_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1))
779 isq, ieq-isq+2, tv%eqn_of_state)
781 isq, ieq-isq+2, tv%eqn_of_state)
783 pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1))*c_htot(i,j)) * &
784 ((dr_dt(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + &
785 dr_ds(i)*(tv%S(i,j,k+1)-tv%S(i,j,k))) / rho_in_situ(i)**2)
794 c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
795 pbce(i,j,nz) = dp_dh * alpha_lay(nz)
797 do k=nz-1,1,-1 ;
do i=isq,ieq+1
798 pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1))*c_htot(i,j)) * &
804 if (gfs_scale < 1.0)
then
807 do j=jsq,jeq+1 ;
do i=isq,ieq+1
808 dpbce(i,j) = (gfs_scale - 1.0) * pbce(i,j,1)
811 do k=1,nz ;
do j=jsq,jeq+1 ;
do i=isq,ieq+1
812 pbce(i,j,k) = pbce(i,j,k) + dpbce(i,j)
813 enddo ;
enddo ;
enddo
816 end subroutine set_pbce_nonbouss
819 subroutine pressureforce_mont_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)
820 type(time_type),
target,
intent(in) :: time
825 type(
diag_ctrl),
target,
intent(inout) :: diag
829 logical :: use_temperature, use_eos
831 #include "version_variable.h"
832 character(len=40) :: mdl
834 if (
associated(cs))
then
835 call mom_error(warning,
"PressureForce_init called with an associated "// &
836 "control structure.")
838 else ;
allocate(cs) ;
endif
840 cs%diag => diag ; cs%Time => time
841 if (
present(tides_csp))
then
842 if (
associated(tides_csp)) cs%tides_CSp => tides_csp
845 mdl =
"MOM_PressureForce_Mont"
847 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
848 "The mean ocean density used with BOUSSINESQ true to "//&
849 "calculate accelerations and the mass for conservation "//&
850 "properties, or with BOUSSINSEQ false to convert some "//&
851 "parameters from vertical units of m to kg m-2.", &
852 units=
"kg m-3", default=1035.0)
853 call get_param(param_file, mdl,
"TIDES", cs%tides, &
854 "If true, apply tidal momentum forcing.", default=.false.)
855 call get_param(param_file, mdl,
"USE_EOS", use_eos, default=.true., &
859 cs%id_PFu_bc = register_diag_field(
'ocean_model',
'PFu_bc', diag%axesCuL, time, &
860 'Density Gradient Zonal Pressure Force Accel.',
"meter second-2")
861 cs%id_PFv_bc = register_diag_field(
'ocean_model',
'PFv_bc', diag%axesCvL, time, &
862 'Density Gradient Meridional Pressure Force Accel.',
"meter second-2")
863 if (cs%id_PFu_bc > 0)
then
864 call safe_alloc_ptr(cs%PFu_bc,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
865 cs%PFu_bc(:,:,:) = 0.0
867 if (cs%id_PFv_bc > 0)
then
868 call safe_alloc_ptr(cs%PFv_bc,g%isd,g%ied,g%JsdB,g%JedB,g%ke)
869 cs%PFv_bc(:,:,:) = 0.0
874 cs%id_e_tidal = register_diag_field(
'ocean_model',
'e_tidal', diag%axesT1, &
875 time,
'Tidal Forcing Astronomical and SAL Height Anomaly',
'meter', conversion=us%Z_to_m)
879 if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
881 call log_param(param_file, mdl,
"GFS / G_EARTH", cs%GFS_scale)
883 end subroutine pressureforce_mont_init
886 subroutine pressureforce_mont_end(CS)
888 if (
associated(cs))
deallocate(cs)
889 end subroutine pressureforce_mont_end