18 implicit none ;
private
20 #include <MOM_memory.h>
40 type(time_type),
pointer :: time => null()
43 real,
pointer :: pfu_bc(:,:,:) => null()
45 real,
pointer :: pfv_bc(:,:,:) => null()
48 integer :: id_pfu_bc = -1, id_pfv_bc = -1, id_e_tidal = -1
64 subroutine pressureforce_mont_nonbouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
68 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
70 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: pfu
72 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: pfv
75 real,
dimension(:,:),
optional,
pointer :: p_atm
77 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
78 optional,
intent(out) :: pbce
81 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: eta
84 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
88 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p
92 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
target :: &
98 real,
dimension(SZI_(G)) :: rho_cv_bl
101 real,
dimension(SZI_(G),SZJ_(G)) :: &
111 real :: p_ref(szi_(g))
113 real :: rho_in_situ(szi_(g))
114 real :: pfu_bc, pfv_bc
131 real :: alpha_lay(szk_(g))
132 real :: dalpha_int(szk_(g)+1)
134 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
136 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
137 nkmb=gv%nk_rho_varies
138 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
141 if (
present(p_atm))
then ;
if (
associated(p_atm)) use_p_atm = .true. ;
endif
142 is_split = .false. ;
if (
present(pbce)) is_split = .true.
143 use_eos =
associated(tv%eqn_of_state)
145 if (.not.
associated(cs))
call mom_error(fatal, &
146 "MOM_PressureForce_Mont: Module must be initialized before it is used.")
149 "PressureForce_Mont_nonBouss: The Montgomery form of the pressure force "//&
150 "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.")
153 pa_to_p_dyn = us%kg_m3_to_R * us%m_s_to_L_T**2
154 i_gearth = 1.0 / (us%L_T_to_m_s**2 * gv%g_Earth)
155 dp_neglect = gv%H_to_Pa * gv%H_subroundoff
156 do k=1,nz ; alpha_lay(k) = 1.0 / (gv%Rlay(k)) ;
enddo
157 do k=2,nz ; dalpha_int(k) = alpha_lay(k-1) - alpha_lay(k) ;
enddo
161 do j=jsq,jeq+1 ;
do i=isq,ieq+1 ; p(i,j,1) = p_atm(i,j) ;
enddo ;
enddo
164 do j=jsq,jeq+1 ;
do i=isq,ieq+1 ; p(i,j,1) = 0.0 ;
enddo ;
enddo
167 do j=jsq,jeq+1 ;
do k=1,nz ;
do i=isq,ieq+1
168 p(i,j,k+1) = p(i,j,k) + gv%H_to_Pa * h(i,j,k)
169 enddo ;
enddo ;
enddo
171 if (
present(eta))
then
172 pa_to_h = 1.0 / gv%H_to_Pa
175 do j=jsq,jeq+1 ;
do i=isq,ieq+1
176 eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*pa_to_h
180 do j=jsq,jeq+1 ;
do i=isq,ieq+1
181 eta(i,j) = p(i,j,nz+1)*pa_to_h
190 do j=jsq,jeq+1 ;
do i=isq,ieq+1
191 ssh(i,j) = -g%bathyT(i,j)
196 call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), &
197 0.0, g%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1)
200 do j=jsq,jeq+1 ;
do k=1,nz;
do i=isq,ieq+1
201 ssh(i,j) = ssh(i,j) + i_gearth * dz_geo(i,j,k)
202 enddo ;
enddo ;
enddo
205 do j=jsq,jeq+1 ;
do k=1,nz ;
do i=isq,ieq+1
206 ssh(i,j) = ssh(i,j) + gv%H_to_RZ * h(i,j,k) * alpha_lay(k)
207 enddo ;
enddo ;
enddo
212 do j=jsq,jeq+1 ;
do i=isq,ieq+1
213 geopot_bot(i,j) = -gv%g_Earth*(e_tidal(i,j) + g%bathyT(i,j))
217 do j=jsq,jeq+1 ;
do i=isq,ieq+1
218 geopot_bot(i,j) = -gv%g_Earth*g%bathyT(i,j)
230 tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
231 tv_tmp%eqn_of_state => tv%eqn_of_state
232 do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ;
enddo
235 do k=1,nkmb ;
do i=isq,ieq+1
236 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
239 rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
240 do k=nkmb+1,nz ;
do i=isq,ieq+1
241 if (gv%Rlay(k) < rho_cv_bl(i))
then
242 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
244 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
249 tv_tmp%T => tv%T ; tv_tmp%S => tv%S
250 tv_tmp%eqn_of_state => tv%eqn_of_state
251 do i=isq,ieq+1 ; p_ref(i) = 0 ;
enddo
254 do k=1,nz ;
do j=jsq,jeq+1
256 rho_in_situ,isq,ieq-isq+2,tv%eqn_of_state, scale=us%kg_m3_to_R)
257 do i=isq,ieq+1 ; alpha_star(i,j,k) = 1.0 / rho_in_situ(i) ;
enddo
265 m(i,j,nz) = geopot_bot(i,j) + pa_to_p_dyn*p(i,j,nz+1) * alpha_star(i,j,nz)
267 do k=nz-1,1,-1 ;
do i=isq,ieq+1
268 m(i,j,k) = m(i,j,k+1) + pa_to_p_dyn*p(i,j,k+1) * (alpha_star(i,j,k) - alpha_star(i,j,k+1))
275 m(i,j,nz) = geopot_bot(i,j) + pa_to_p_dyn*p(i,j,nz+1) * alpha_lay(nz)
277 do k=nz-1,1,-1 ;
do i=isq,ieq+1
278 m(i,j,k) = m(i,j,k+1) + pa_to_p_dyn*p(i,j,k+1) * dalpha_int(k+1)
283 if (cs%GFS_scale < 1.0)
then
286 do j=jsq,jeq+1 ;
do i=isq,ieq+1
287 dm(i,j) = (cs%GFS_scale - 1.0) * m(i,j,1)
290 do k=1,nz ;
do j=jsq,jeq+1 ;
do i=isq,ieq+1
291 m(i,j,k) = m(i,j,k) + dm(i,j)
292 enddo ;
enddo ;
enddo
312 if (
present(pbce))
then
321 do j=jsq,jeq+1 ;
do i=isq,ieq+1
322 dp_star(i,j) = (p(i,j,k+1) - p(i,j,k)) + dp_neglect
324 do j=js,je ;
do i=isq,ieq
326 pfu_bc = (alpha_star(i+1,j,k) - alpha_star(i,j,k)) * (g%IdxCu(i,j) * pa_to_p_dyn * &
327 ((dp_star(i,j) * dp_star(i+1,j) + (p(i,j,k) * dp_star(i+1,j) + &
328 p(i+1,j,k) * dp_star(i,j))) / (dp_star(i,j) + dp_star(i+1,j))))
329 pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j) + pfu_bc
330 if (
associated(cs%PFu_bc)) cs%PFu_bc(i,j,k) = pfu_bc
332 do j=jsq,jeq ;
do i=is,ie
333 pfv_bc = (alpha_star(i,j+1,k) - alpha_star(i,j,k)) * (g%IdyCv(i,j) * pa_to_p_dyn * &
334 ((dp_star(i,j) * dp_star(i,j+1) + (p(i,j,k) * dp_star(i,j+1) + &
335 p(i,j+1,k) * dp_star(i,j))) / (dp_star(i,j) + dp_star(i,j+1))))
336 pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j) + pfv_bc
337 if (
associated(cs%PFv_bc)) cs%PFv_bc(i,j,k) = pfv_bc
343 do j=js,je ;
do i=isq,ieq
344 pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j)
346 do j=jsq,jeq ;
do i=is,ie
347 pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j)
352 if (cs%id_PFu_bc>0)
call post_data(cs%id_PFu_bc, cs%PFu_bc, cs%diag)
353 if (cs%id_PFv_bc>0)
call post_data(cs%id_PFv_bc, cs%PFv_bc, cs%diag)
354 if (cs%id_e_tidal>0)
call post_data(cs%id_e_tidal, e_tidal, cs%diag)
365 subroutine pressureforce_mont_bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
369 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
371 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: pfu
373 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: pfv
376 real,
dimension(:,:),
optional,
pointer :: p_atm
378 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(out) :: pbce
381 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: eta
383 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
387 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e
391 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
target :: &
397 real :: rho_cv_bl(szi_(g))
399 real :: h_star(szi_(g),szj_(g))
401 real :: e_tidal(szi_(g),szj_(g))
404 real :: p_ref(szi_(g))
408 real :: pfu_bc, pfv_bc
419 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
422 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
423 nkmb=gv%nk_rho_varies
424 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
427 if (
present(p_atm))
then ;
if (
associated(p_atm)) use_p_atm = .true. ;
endif
428 is_split = .false. ;
if (
present(pbce)) is_split = .true.
429 use_eos =
associated(tv%eqn_of_state)
431 if (.not.
associated(cs))
call mom_error(fatal, &
432 "MOM_PressureForce_Mont: Module must be initialized before it is used.")
435 "PressureForce_Mont_Bouss: The Montgomery form of the pressure force "//&
436 "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.")
439 h_neglect = gv%H_subroundoff * gv%H_to_Z
441 g_rho0 = gv%g_Earth / gv%Rho0
450 do i=isq,ieq+1 ; e(i,j,1) = -g%bathyT(i,j) ;
enddo
451 do k=1,nz ;
do i=isq,ieq+1
452 e(i,j,1) = e(i,j,1) + h(i,j,k)*gv%H_to_Z
455 call calc_tidal_forcing(cs%Time, e(:,:,1), e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
461 do j=jsq,jeq+1 ;
do i=isq,ieq+1
462 e(i,j,nz+1) = -(g%bathyT(i,j) + e_tidal(i,j))
466 do j=jsq,jeq+1 ;
do i=isq,ieq+1
467 e(i,j,nz+1) = -g%bathyT(i,j)
471 do j=jsq,jeq+1 ;
do k=nz,1,-1 ;
do i=isq,ieq+1
472 e(i,j,k) = e(i,j,k+1) + h(i,j,k)*gv%H_to_Z
473 enddo ;
enddo ;
enddo
484 tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
485 tv_tmp%eqn_of_state => tv%eqn_of_state
487 do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ;
enddo
490 do k=1,nkmb ;
do i=isq,ieq+1
491 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
494 rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
496 do k=nkmb+1,nz ;
do i=isq,ieq+1
497 if (gv%Rlay(k) < rho_cv_bl(i))
then
498 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
500 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
505 tv_tmp%T => tv%T ; tv_tmp%S => tv%S
506 tv_tmp%eqn_of_state => tv%eqn_of_state
507 do i=isq,ieq+1 ; p_ref(i) = 0.0 ;
enddo
513 do k=1,nz+1 ;
do j=jsq,jeq+1
514 call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), &
515 isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R*g_rho0)
524 m(i,j,1) = cs%GFS_scale * (rho_star(i,j,1) * e(i,j,1))
525 if (use_p_atm) m(i,j,1) = m(i,j,1) + us%m_s_to_L_T**2*p_atm(i,j) * i_rho0
527 do k=2,nz ;
do i=isq,ieq+1
528 m(i,j,k) = m(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) * e(i,j,k)
535 m(i,j,1) = gv%g_prime(1) * e(i,j,1)
536 if (use_p_atm) m(i,j,1) = m(i,j,1) + us%m_s_to_L_T**2*p_atm(i,j) * i_rho0
538 do k=2,nz ;
do i=isq,ieq+1
539 m(i,j,k) = m(i,j,k-1) + gv%g_prime(k) * e(i,j,k)
544 if (
present(pbce))
then
545 call set_pbce_bouss(e, tv_tmp, g, gv, us, cs%Rho0, cs%GFS_scale, pbce, rho_star)
553 do j=jsq,jeq+1 ;
do i=isq,ieq+1
554 h_star(i,j) = (e(i,j,k) - e(i,j,k+1)) + h_neglect
556 do j=js,je ;
do i=isq,ieq
557 pfu_bc = -1.0*(rho_star(i+1,j,k) - rho_star(i,j,k)) * (g%IdxCu(i,j) * &
558 ((h_star(i,j) * h_star(i+1,j) - (e(i,j,k) * h_star(i+1,j) + &
559 e(i+1,j,k) * h_star(i,j))) / (h_star(i,j) + h_star(i+1,j))))
560 pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j) + pfu_bc
561 if (
associated(cs%PFu_bc)) cs%PFu_bc(i,j,k) = pfu_bc
563 do j=jsq,jeq ;
do i=is,ie
564 pfv_bc = -1.0*(rho_star(i,j+1,k) - rho_star(i,j,k)) * (g%IdyCv(i,j) * &
565 ((h_star(i,j) * h_star(i,j+1) - (e(i,j,k) * h_star(i,j+1) + &
566 e(i,j+1,k) * h_star(i,j))) / (h_star(i,j) + h_star(i,j+1))))
567 pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j) + pfv_bc
568 if (
associated(cs%PFv_bc)) cs%PFv_bc(i,j,k) = pfv_bc
574 do j=js,je ;
do i=isq,ieq
575 pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j)
577 do j=jsq,jeq ;
do i=is,ie
578 pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j)
583 if (
present(eta))
then
589 do j=jsq,jeq+1 ;
do i=isq,ieq+1
590 eta(i,j) = e(i,j,1)*gv%Z_to_H + e_tidal(i,j)*gv%Z_to_H
594 do j=jsq,jeq+1 ;
do i=isq,ieq+1
595 eta(i,j) = e(i,j,1)*gv%Z_to_H
600 if (cs%id_PFu_bc>0)
call post_data(cs%id_PFu_bc, cs%PFu_bc, cs%diag)
601 if (cs%id_PFv_bc>0)
call post_data(cs%id_PFv_bc, cs%PFv_bc, cs%diag)
602 if (cs%id_e_tidal>0)
call post_data(cs%id_e_tidal, e_tidal, cs%diag)
608 subroutine set_pbce_bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
611 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: e
614 real,
intent(in) :: rho0
615 real,
intent(in) :: gfs_scale
618 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
622 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
623 optional,
intent(in) :: rho_star
627 real :: ihtot(szi_(g))
628 real :: press(szi_(g))
629 real :: t_int(szi_(g))
630 real :: s_int(szi_(g))
631 real :: dr_dt(szi_(g))
632 real :: dr_ds(szi_(g))
633 real :: rho_in_situ(szi_(g))
640 integer :: isq, ieq, jsq, jeq, nz, i, j, k
642 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
644 rho0xg = rho0*us%L_T_to_m_s**2 * gv%g_Earth
645 g_rho0 = gv%g_Earth / gv%Rho0
646 use_eos =
associated(tv%eqn_of_state)
647 z_neglect = gv%H_subroundoff*gv%H_to_Z
650 if (
present(rho_star))
then
654 ihtot(i) = gv%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
655 pbce(i,j,1) = gfs_scale * rho_star(i,j,1) * gv%H_to_Z
657 do k=2,nz ;
do i=isq,ieq+1
658 pbce(i,j,k) = pbce(i,j,k-1) + (rho_star(i,j,k)-rho_star(i,j,k-1)) * &
659 ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i))
666 ihtot(i) = gv%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
667 press(i) = -rho0xg*e(i,j,1)
670 isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
672 pbce(i,j,1) = g_rho0*(gfs_scale * rho_in_situ(i)) * gv%H_to_Z
676 press(i) = -rho0xg*e(i,j,k)
677 t_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
678 s_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
681 isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
683 pbce(i,j,k) = pbce(i,j,k-1) + g_rho0 * &
684 ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i)) * &
685 (dr_dt(i)*(tv%T(i,j,k)-tv%T(i,j,k-1)) + &
686 dr_ds(i)*(tv%S(i,j,k)-tv%S(i,j,k-1)))
695 ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
696 pbce(i,j,1) = gv%g_prime(1) * gv%H_to_Z
698 do k=2,nz ;
do i=isq,ieq+1
699 pbce(i,j,k) = pbce(i,j,k-1) + &
700 (gv%g_prime(k)*gv%H_to_Z) * ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i))
712 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: p
715 real,
intent(in) :: gfs_scale
718 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: pbce
721 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(in) :: alpha_star
724 real,
dimension(SZI_(G),SZJ_(G)) :: &
728 real :: t_int(szi_(g))
729 real :: s_int(szi_(g))
730 real :: dr_dt(szi_(g))
731 real :: dr_ds(szi_(g))
732 real :: rho_in_situ(szi_(g))
733 real :: alpha_lay(szk_(g))
734 real :: dalpha_int(szk_(g)+1)
741 integer :: isq, ieq, jsq, jeq, nz, i, j, k
743 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
745 use_eos =
associated(tv%eqn_of_state)
747 dp_dh = gv%g_Earth * gv%H_to_RZ
748 dp_neglect = gv%H_to_Pa * gv%H_subroundoff
751 if (
present(alpha_star))
then
755 c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
756 pbce(i,j,nz) = dp_dh * alpha_star(i,j,nz)
758 do k=nz-1,1,-1 ;
do i=isq,ieq+1
759 pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1)) * c_htot(i,j)) * &
760 (alpha_star(i,j,k) - alpha_star(i,j,k+1))
767 rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
769 c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
770 pbce(i,j,nz) = dp_dh / (rho_in_situ(i))
774 t_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1))
775 s_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1))
778 isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
780 isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
782 pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1))*c_htot(i,j)) * &
783 ((dr_dt(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + &
784 dr_ds(i)*(tv%S(i,j,k+1)-tv%S(i,j,k))) / (rho_in_situ(i)**2))
791 do k=1,nz ; alpha_lay(k) = 1.0 / (gv%Rlay(k)) ;
enddo
792 do k=2,nz ; dalpha_int(k) = alpha_lay(k-1) - alpha_lay(k) ;
enddo
797 c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
798 pbce(i,j,nz) = dp_dh * alpha_lay(nz)
800 do k=nz-1,1,-1 ;
do i=isq,ieq+1
801 pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1))*c_htot(i,j)) * &
807 if (gfs_scale < 1.0)
then
810 do j=jsq,jeq+1 ;
do i=isq,ieq+1
811 dpbce(i,j) = (gfs_scale - 1.0) * pbce(i,j,1)
814 do k=1,nz ;
do j=jsq,jeq+1 ;
do i=isq,ieq+1
815 pbce(i,j,k) = pbce(i,j,k) + dpbce(i,j)
816 enddo ;
enddo ;
enddo
823 type(time_type),
target,
intent(in) :: time
828 type(
diag_ctrl),
target,
intent(inout) :: diag
833 logical :: use_temperature, use_eos
835 # include "version_variable.h"
836 character(len=40) :: mdl
838 if (
associated(cs))
then
839 call mom_error(warning,
"PressureForce_init called with an associated "// &
840 "control structure.")
842 else ;
allocate(cs) ;
endif
844 cs%diag => diag ; cs%Time => time
845 if (
present(tides_csp))
then
846 if (
associated(tides_csp)) cs%tides_CSp => tides_csp
849 mdl =
"MOM_PressureForce_Mont"
851 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
852 "The mean ocean density used with BOUSSINESQ true to "//&
853 "calculate accelerations and the mass for conservation "//&
854 "properties, or with BOUSSINSEQ false to convert some "//&
855 "parameters from vertical units of m to kg m-2.", &
856 units=
"kg m-3", default=1035.0)
857 call get_param(param_file, mdl,
"TIDES", cs%tides, &
858 "If true, apply tidal momentum forcing.", default=.false.)
859 call get_param(param_file, mdl,
"USE_EOS", use_eos, default=.true., &
863 cs%id_PFu_bc = register_diag_field(
'ocean_model',
'PFu_bc', diag%axesCuL, time, &
864 'Density Gradient Zonal Pressure Force Accel.',
"meter second-2", conversion=us%L_T2_to_m_s2)
865 cs%id_PFv_bc = register_diag_field(
'ocean_model',
'PFv_bc', diag%axesCvL, time, &
866 'Density Gradient Meridional Pressure Force Accel.',
"meter second-2", conversion=us%L_T2_to_m_s2)
867 if (cs%id_PFu_bc > 0)
then
868 call safe_alloc_ptr(cs%PFu_bc,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
869 cs%PFu_bc(:,:,:) = 0.0
871 if (cs%id_PFv_bc > 0)
then
872 call safe_alloc_ptr(cs%PFv_bc,g%isd,g%ied,g%JsdB,g%JedB,g%ke)
873 cs%PFv_bc(:,:,:) = 0.0
878 cs%id_e_tidal = register_diag_field(
'ocean_model',
'e_tidal', diag%axesT1, &
879 time,
'Tidal Forcing Astronomical and SAL Height Anomaly',
'meter', conversion=us%Z_to_m)
883 if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
885 call log_param(param_file, mdl,
"GFS / G_EARTH", cs%GFS_scale)
892 if (
associated(cs))
deallocate(cs)