14 implicit none ;
private
16 #include <MOM_memory.h>
67 real,
parameter ::
a0 = 7.057924e-4,
a1 = 3.480336e-7,
a2 = -1.112733e-7
68 real,
parameter ::
b0 = 5.790749e8,
b1 = 3.516535e6,
b2 = -4.002714e4
69 real,
parameter ::
b3 = 2.084372e2,
b4 = 5.944068e5,
b5 = -9.643486e3
70 real,
parameter ::
c0 = 1.704853e5,
c1 = 7.904722e2,
c2 = -7.984422
71 real,
parameter ::
c3 = 5.140652e-2,
c4 = -2.302158e2,
c5 = -3.079464
83 real,
intent(in) :: pressure
84 real,
intent(out) :: rho
85 real,
optional,
intent(in) :: rho_ref
95 real,
dimension(1) :: T0, S0, pressure0, rho0
99 pressure0(1) = pressure
111 real,
dimension(:),
intent(in) :: T
112 real,
dimension(:),
intent(in) :: S
113 real,
dimension(:),
intent(in) :: pressure
114 real,
dimension(:),
intent(out) :: rho
115 integer,
intent(in) :: start
116 integer,
intent(in) :: npts
117 real,
optional,
intent(in) :: rho_ref
121 real :: al0, p0, lambda
122 real :: al_TS, p_TSp, lam_TS, pa_000
125 if (
present(rho_ref)) pa_000 = (
b0*(1.0 -
a0*rho_ref) - rho_ref*
c0)
126 if (
present(rho_ref))
then ;
do j=start,start+npts-1
127 al_ts =
a1*t(j) +
a2*s(j)
129 p_tsp = pressure(j) + (
b4*s(j) + t(j) * (
b1 + (t(j)*(
b2 +
b3*t(j)) +
b5*s(j))))
130 lam_ts =
c4*s(j) + t(j) * (
c1 + (t(j)*(
c2 +
c3*t(j)) +
c5*s(j)))
134 rho(j) = (pa_000 + (p_tsp - rho_ref*(p_tsp*al0 + (
b0*al_ts + lam_ts)))) / &
135 ( (
c0 + lam_ts) + al0*(
b0 + p_tsp) )
136 enddo ;
else ;
do j=start,start+npts-1
137 al0 = (
a0 +
a1*t(j)) +
a2*s(j)
138 p0 = (
b0 +
b4*s(j)) + t(j) * (
b1 + t(j)*(
b2 +
b3*t(j)) +
b5*s(j))
139 lambda = (
c0 +
c4*s(j)) + t(j) * (
c1 + t(j)*(
c2 +
c3*t(j)) +
c5*s(j))
140 rho(j) = (pressure(j) + p0) / (lambda + al0*(pressure(j) + p0))
151 real,
intent(in) :: T
152 real,
intent(in) :: S
153 real,
intent(in) :: pressure
154 real,
intent(out) :: specvol
155 real,
optional,
intent(in) :: spv_ref
158 real,
dimension(1) :: T0, S0, pressure0, spv0
160 t0(1) = t ; s0(1) = s ; pressure0(1) = pressure
172 real,
dimension(:),
intent(in) :: T
174 real,
dimension(:),
intent(in) :: S
175 real,
dimension(:),
intent(in) :: pressure
176 real,
dimension(:),
intent(out) :: specvol
177 integer,
intent(in) :: start
178 integer,
intent(in) :: npts
179 real,
optional,
intent(in) :: spv_ref
182 real :: al0, p0, lambda
185 do j=start,start+npts-1
186 al0 = (
a0 +
a1*t(j)) +
a2*s(j)
187 p0 = (
b0 +
b4*s(j)) + t(j) * (
b1 + t(j)*((
b2 +
b3*t(j))) +
b5*s(j))
188 lambda = (
c0 +
c4*s(j)) + t(j) * (
c1 + t(j)*((
c2 +
c3*t(j))) +
c5*s(j))
190 if (
present(spv_ref))
then
191 specvol(j) = (lambda + (al0 - spv_ref)*(pressure(j) + p0)) / (pressure(j) + p0)
193 specvol(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0)
200 real,
intent(in),
dimension(:) :: T
202 real,
intent(in),
dimension(:) :: S
203 real,
intent(in),
dimension(:) :: pressure
204 real,
intent(out),
dimension(:) :: drho_dT
206 real,
intent(out),
dimension(:) :: drho_dS
208 integer,
intent(in) :: start
209 integer,
intent(in) :: npts
212 real :: al0, p0, lambda, I_denom2
215 do j=start,start+npts-1
216 al0 = (
a0 +
a1*t(j)) +
a2*s(j)
217 p0 = (
b0 +
b4*s(j)) + t(j) * (
b1 + t(j)*((
b2 +
b3*t(j))) +
b5*s(j))
218 lambda = (
c0 +
c4*s(j)) + t(j) * (
c1 + t(j)*((
c2 +
c3*t(j))) +
c5*s(j))
220 i_denom2 = 1.0 / (lambda + al0*(pressure(j) + p0))
221 i_denom2 = i_denom2 *i_denom2
222 drho_dt(j) = i_denom2 * &
223 (lambda* (
b1 + t(j)*(2.0*
b2 + 3.0*
b3*t(j)) +
b5*s(j)) - &
224 (pressure(j)+p0) * ( (pressure(j)+p0)*
a1 + &
225 (
c1 + t(j)*(
c2*2.0 +
c3*3.0*t(j)) +
c5*s(j)) ))
226 drho_ds(j) = i_denom2 * (lambda* (
b4 +
b5*t(j)) - &
227 (pressure(j)+p0) * ( (pressure(j)+p0)*
a2 + (
c4 +
c5*t(j)) ))
235 real,
intent(in) :: T
236 real,
intent(in) :: S
237 real,
intent(in) :: pressure
238 real,
intent(out) :: drho_dT
240 real,
intent(out) :: drho_dS
244 real,
dimension(1) :: T0, S0, P0
245 real,
dimension(1) :: drdt0, drds0
258 drho_ds_dp, drho_dt_dp, start, npts)
259 real,
dimension(:),
intent(in ) :: T
260 real,
dimension(:),
intent(in ) :: S
261 real,
dimension(:),
intent(in ) :: P
262 real,
dimension(:),
intent( out) :: drho_ds_ds
264 real,
dimension(:),
intent( out) :: drho_ds_dt
266 real,
dimension(:),
intent( out) :: drho_dt_dt
268 real,
dimension(:),
intent( out) :: drho_ds_dp
270 real,
dimension(:),
intent( out) :: drho_dt_dp
272 integer,
intent(in ) :: start
273 integer,
intent(in ) :: npts
276 real :: z0, z1, z2, z3, z4, z5, z6 ,z7, z8, z9, z10, z11, z2_2, z2_3
281 do j = start,start+npts-1
282 z0 = t(j)*(
b1 +
b5*s(j) + t(j)*(
b2 +
b3*t(j)))
283 z1 = (
b0 + p(j) +
b4*s(j) + z0)
284 z3 = (
b1 +
b5*s(j) + t(j)*(2.*
b2 + 2.*
b3*t(j)))
285 z4 = (
c0 +
c4*s(j) + t(j)*(
c1 +
c5*s(j) + t(j)*(
c2 +
c3*t(j))))
286 z5 = (
b1 +
b5*s(j) + t(j)*(
b2 +
b3*t(j)) + t(j)*(
b2 + 2.*
b3*t(j)))
287 z6 =
c1 +
c5*s(j) + t(j)*(
c2 +
c3*t(j)) + t(j)*(
c2 + 2.*
c3*t(j))
288 z7 = (
c4 +
c5*t(j) +
a2*z1)
289 z8 = (
c1 +
c5*s(j) + t(j)*(2.*
c2 + 3.*
c3*t(j)) +
a1*z1)
290 z9 = (
a0 +
a2*s(j) +
a1*t(j))
292 z11 = (z10*z4 - z1*z7)
293 z2 = (
c0 +
c4*s(j) + t(j)*(
c1 +
c5*s(j) + t(j)*(
c2 +
c3*t(j))) + z9*z1)
297 drho_ds_ds(j) = (z10*(
c4 +
c5*t(j)) -
a2*z10*z1 - z10*z7)/z2_2 - (2.*(
c4 +
c5*t(j) + z9*z10 +
a2*z1)*z11)/z2_3
298 drho_ds_dt(j) = (z10*z6 - z1*(
c5 +
a2*z5) +
b5*z4 - z5*z7)/z2_2 - (2.*(z6 + z9*z5 +
a1*z1)*z11)/z2_3
299 drho_dt_dt(j) = (z3*z6 - z1*(2.*
c2 + 6.*
c3*t(j) +
a1*z5) + (2.*
b2 + 4.*
b3*t(j))*z4 - z5*z8)/z2_2 - &
300 (2.*(z6 + z9*z5 +
a1*z1)*(z3*z4 - z1*z8))/z2_3
301 drho_ds_dp(j) = (-
c4 -
c5*t(j) - 2.*
a2*z1)/z2_2 - (2.*z9*z11)/z2_3
302 drho_dt_dp(j) = (-
c1 -
c5*s(j) - t(j)*(2.*
c2 + 3.*
c3*t(j)) - 2.*
a1*z1)/z2_2 - (2.*z9*(z3*z4 - z1*z8))/z2_3
310 drho_ds_dp, drho_dt_dp)
311 real,
intent(in ) :: T
312 real,
intent(in ) :: S
313 real,
intent(in ) :: P
314 real,
intent( out) :: drho_ds_ds
316 real,
intent( out) :: drho_ds_dt
318 real,
intent( out) :: drho_dt_dt
320 real,
intent( out) :: drho_ds_dp
322 real,
intent( out) :: drho_dt_dp
325 real,
dimension(1) :: T0, S0, P0
326 real,
dimension(1) :: drdsds, drdsdt, drdtdt, drdsdp, drdtdp
332 drho_ds_ds = drdsds(1)
333 drho_ds_dt = drdsdt(1)
334 drho_dt_dt = drdtdt(1)
335 drho_ds_dp = drdsdp(1)
336 drho_dt_dp = drdtdp(1)
343 real,
intent(in),
dimension(:) :: t
344 real,
intent(in),
dimension(:) :: s
345 real,
intent(in),
dimension(:) :: pressure
346 real,
intent(out),
dimension(:) :: dsv_dt
348 real,
intent(out),
dimension(:) :: dsv_ds
350 integer,
intent(in) :: start
351 integer,
intent(in) :: npts
354 real :: al0, p0, lambda, i_denom
357 do j=start,start+npts-1
359 p0 = (
b0 +
b4*s(j)) + t(j) * (
b1 + t(j)*((
b2 +
b3*t(j))) +
b5*s(j))
360 lambda = (
c0 +
c4*s(j)) + t(j) * (
c1 + t(j)*((
c2 +
c3*t(j))) +
c5*s(j))
364 i_denom = 1.0 / (pressure(j) + p0)
365 dsv_dt(j) = (
a1 + i_denom * (
c1 + t(j)*((2.0*
c2 + 3.0*
c3*t(j))) +
c5*s(j))) - &
366 (i_denom**2 * lambda) * (
b1 + t(j)*((2.0*
b2 + 3.0*
b3*t(j))) +
b5*s(j))
367 dsv_ds(j) = (
a2 + i_denom * (
c4 +
c5*t(j))) - &
368 (i_denom**2 * lambda) * (
b4 +
b5*t(j))
380 real,
intent(in),
dimension(:) :: t
381 real,
intent(in),
dimension(:) :: s
382 real,
intent(in),
dimension(:) :: pressure
383 real,
intent(out),
dimension(:) :: rho
384 real,
intent(out),
dimension(:) :: drho_dp
387 integer,
intent(in) :: start
388 integer,
intent(in) :: npts
392 real :: al0, p0, lambda, i_denom
395 do j=start,start+npts-1
396 al0 = (
a0 +
a1*t(j)) +
a2*s(j)
397 p0 = (
b0 +
b4*s(j)) + t(j) * (
b1 + t(j)*((
b2 +
b3*t(j))) +
b5*s(j))
398 lambda = (
c0 +
c4*s(j)) + t(j) * (
c1 + t(j)*((
c2 +
c3*t(j))) +
c5*s(j))
400 i_denom = 1.0 / (lambda + al0*(pressure(j) + p0))
401 rho(j) = (pressure(j) + p0) * i_denom
402 drho_dp(j) = lambda * i_denom * i_denom
410 dpa, intz_dpa, intx_dpa, inty_dpa, &
411 bathyT, dz_neglect, useMassWghtInterp)
414 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
417 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
419 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
421 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
423 real,
intent(in) :: rho_ref
426 real,
intent(in) :: rho_0
429 real,
intent(in) :: g_e
430 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
433 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
434 optional,
intent(out) :: intz_dpa
437 real,
dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
438 optional,
intent(out) :: intx_dpa
441 real,
dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
442 optional,
intent(out) :: inty_dpa
445 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
446 optional,
intent(in) :: bathyt
447 real,
optional,
intent(in) :: dz_neglect
448 logical,
optional,
intent(in) :: usemasswghtinterp
452 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed) :: al0_2d, p0_2d, lambda_2d
453 real :: al0, p0, lambda
455 real :: eps, eps2, rem
457 real :: p_ave, i_al0, i_lzz
462 real :: hwt_ll, hwt_lr
463 real :: hwt_rl, hwt_rr
468 logical :: do_massweight
469 real,
parameter :: c1_3 = 1.0/3.0, c1_7 = 1.0/7.0
470 real,
parameter :: c1_9 = 1.0/9.0, c1_90 = 1.0/90.0
471 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, ioff, joff, m
473 ioff = hio%idg_offset - hii%idg_offset
474 joff = hio%jdg_offset - hii%jdg_offset
478 isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
479 jsq = hio%JscB + joff ; jeq = hio%JecB + joff
480 is = hio%isc + ioff ; ie = hio%iec + ioff
481 js = hio%jsc + joff ; je = hio%jec + joff
486 do_massweight = .false.
487 if (
present(usemasswghtinterp))
then ;
if (usemasswghtinterp)
then
488 do_massweight = .true.
495 do j=jsq,jeq+1 ;
do i=isq,ieq+1
496 al0_2d(i,j) = (
a0 +
a1*t(i,j)) +
a2*s(i,j)
497 p0_2d(i,j) = (
b0 +
b4*s(i,j)) + t(i,j) * (
b1 + t(i,j)*((
b2 +
b3*t(i,j))) +
b5*s(i,j))
498 lambda_2d(i,j) = (
c0 +
c4*s(i,j)) + t(i,j) * (
c1 + t(i,j)*((
c2 +
c3*t(i,j))) +
c5*s(i,j))
500 al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j)
502 dz = z_t(i,j) - z_b(i,j)
503 p_ave = -0.5*gxrho*(z_t(i,j)+z_b(i,j))
506 i_lzz = 1.0 / (p0 + (lambda * i_al0) + p_ave)
507 eps = 0.5*gxrho*dz*i_lzz ; eps2 = eps*eps
511 rho_anom = (p0 + p_ave)*(i_lzz*i_al0) - rho_ref
512 rem = i_rho * (lambda * i_al0**2) * eps2 * &
513 (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
514 dpa(i-ioff,j-joff) = g_e*rho_anom*dz - 2.0*eps*rem
515 if (
present(intz_dpa)) &
516 intz_dpa(i-ioff,j-joff) = 0.5*g_e*rho_anom*dz**2 - dz*(1.0+eps)*rem
519 if (
present(intx_dpa))
then ;
do j=js,je ;
do i=isq,ieq
525 hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
527 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
528 hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
529 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
530 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
531 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
532 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
534 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
537 intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
539 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
540 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
542 al0 = wtt_l*al0_2d(i,j) + wtt_r*al0_2d(i+1,j)
543 p0 = wtt_l*p0_2d(i,j) + wtt_r*p0_2d(i+1,j)
544 lambda = wtt_l*lambda_2d(i,j) + wtt_r*lambda_2d(i+1,j)
546 dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i+1,j) - z_b(i+1,j))
547 p_ave = -0.5*gxrho*(wt_l*(z_t(i,j)+z_b(i,j)) + &
548 wt_r*(z_t(i+1,j)+z_b(i+1,j)))
551 i_lzz = 1.0 / (p0 + (lambda * i_al0) + p_ave)
552 eps = 0.5*gxrho*dz*i_lzz ; eps2 = eps*eps
554 intz(m) = g_e*dz*((p0 + p_ave)*(i_lzz*i_al0) - rho_ref) - 2.0*eps * &
555 i_rho * (lambda * i_al0**2) * eps2 * &
556 (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
559 intx_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
561 enddo ;
enddo ;
endif
563 if (
present(inty_dpa))
then ;
do j=jsq,jeq ;
do i=is,ie
569 hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
571 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
572 hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
573 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
574 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
575 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
576 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
578 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
581 intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i-ioff,j+1-joff)
583 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
584 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
586 al0 = wtt_l*al0_2d(i,j) + wtt_r*al0_2d(i,j+1)
587 p0 = wtt_l*p0_2d(i,j) + wtt_r*p0_2d(i,j+1)
588 lambda = wtt_l*lambda_2d(i,j) + wtt_r*lambda_2d(i,j+1)
590 dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i,j+1) - z_b(i,j+1))
591 p_ave = -0.5*gxrho*(wt_l*(z_t(i,j)+z_b(i,j)) + &
592 wt_r*(z_t(i,j+1)+z_b(i,j+1)))
595 i_lzz = 1.0 / (p0 + (lambda * i_al0) + p_ave)
596 eps = 0.5*gxrho*dz*i_lzz ; eps2 = eps*eps
598 intz(m) = g_e*dz*((p0 + p_ave)*(i_lzz*i_al0) - rho_ref) - 2.0*eps * &
599 i_rho * (lambda * i_al0**2) * eps2 * &
600 (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
603 inty_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
605 enddo ;
enddo ;
endif
615 intp_dza, intx_dza, inty_dza, halo_size, &
616 bathyP, dP_neglect, useMassWghtInterp)
618 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
621 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
623 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
625 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
627 real,
intent(in) :: spv_ref
631 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
634 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
635 optional,
intent(out) :: intp_dza
638 real,
dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
639 optional,
intent(out) :: intx_dza
642 real,
dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
643 optional,
intent(out) :: inty_dza
646 integer,
optional,
intent(in) :: halo_size
648 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
649 optional,
intent(in) :: bathyp
650 real,
optional,
intent(in) :: dp_neglect
652 logical,
optional,
intent(in) :: usemasswghtinterp
656 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d, p0_2d, lambda_2d
657 real :: al0, p0, lambda
659 real :: rem, eps, eps2
665 real :: hwt_ll, hwt_lr
666 real :: hwt_rl, hwt_rr
671 logical :: do_massweight
672 real,
parameter :: c1_3 = 1.0/3.0, c1_7 = 1.0/7.0
673 real,
parameter :: c1_9 = 1.0/9.0, c1_90 = 1.0/90.0
674 integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, halo
676 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
677 halo = 0 ;
if (
present(halo_size)) halo = max(halo_size,0)
678 ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
679 if (
present(intx_dza))
then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh);
endif
680 if (
present(inty_dza))
then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh);
endif
682 do_massweight = .false.
683 if (
present(usemasswghtinterp))
then ;
if (usemasswghtinterp)
then
684 do_massweight = .true.
691 do j=jsh,jeh ;
do i=ish,ieh
692 al0_2d(i,j) = (
a0 +
a1*t(i,j)) +
a2*s(i,j)
693 p0_2d(i,j) = (
b0 +
b4*s(i,j)) + t(i,j) * (
b1 + t(i,j)*((
b2 +
b3*t(i,j))) +
b5*s(i,j))
694 lambda_2d(i,j) = (
c0 +
c4*s(i,j)) + t(i,j) * (
c1 + t(i,j)*((
c2 +
c3*t(i,j))) +
c5*s(i,j))
696 al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j)
697 dp = p_b(i,j) - p_t(i,j)
698 p_ave = 0.5*(p_t(i,j)+p_b(i,j))
700 eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
701 alpha_anom = al0 + lambda / (p0 + p_ave) - spv_ref
702 rem = lambda * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
703 dza(i,j) = alpha_anom*dp + 2.0*eps*rem
704 if (
present(intp_dza)) &
705 intp_dza(i,j) = 0.5*alpha_anom*dp**2 - dp*(1.0-eps)*rem
708 if (
present(intx_dza))
then ;
do j=hi%jsc,hi%jec ;
do i=isq,ieq
714 hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
716 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
717 hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
718 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
719 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
720 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
721 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
723 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
726 intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
728 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
729 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
733 al0 = wtt_l*al0_2d(i,j) + wtt_r*al0_2d(i+1,j)
734 p0 = wtt_l*p0_2d(i,j) + wtt_r*p0_2d(i+1,j)
735 lambda = wtt_l*lambda_2d(i,j) + wtt_r*lambda_2d(i+1,j)
737 dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i+1,j) - p_t(i+1,j))
738 p_ave = 0.5*(wt_l*(p_t(i,j)+p_b(i,j)) + wt_r*(p_t(i+1,j)+p_b(i+1,j)))
740 eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
741 intp(m) = (al0 + lambda / (p0 + p_ave) - spv_ref)*dp + 2.0*eps* &
742 lambda * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
745 intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
747 enddo ;
enddo ;
endif
749 if (
present(inty_dza))
then ;
do j=jsq,jeq ;
do i=hi%isc,hi%iec
755 hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
757 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
758 hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
759 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
760 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
761 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
762 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
764 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
767 intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
769 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
770 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
774 al0 = wt_l*al0_2d(i,j) + wt_r*al0_2d(i,j+1)
775 p0 = wt_l*p0_2d(i,j) + wt_r*p0_2d(i,j+1)
776 lambda = wt_l*lambda_2d(i,j) + wt_r*lambda_2d(i,j+1)
778 dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i,j+1) - p_t(i,j+1))
779 p_ave = 0.5*(wt_l*(p_t(i,j)+p_b(i,j)) + wt_r*(p_t(i,j+1)+p_b(i,j+1)))
781 eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
782 intp(m) = (al0 + lambda / (p0 + p_ave) - spv_ref)*dp + 2.0*eps* &
783 lambda * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
786 inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
788 enddo ;
enddo ;
endif