8 implicit none ;
private
10 #include <MOM_memory.h>
57 Rho_T0_S0, dRho_dT, dRho_dS, rho_ref)
60 real,
intent(in) :: pressure
61 real,
intent(out) :: rho
62 real,
intent(in) :: rho_t0_s0
63 real,
intent(in) :: drho_dt
65 real,
intent(in) :: drho_ds
67 real,
optional,
intent(in) :: rho_ref
69 if (
present(rho_ref))
then
70 rho = (rho_t0_s0 - rho_ref) + (drho_dt*t + drho_ds*s)
72 rho = rho_t0_s0 + drho_dt*t + drho_ds*s
81 Rho_T0_S0, dRho_dT, dRho_dS, rho_ref)
82 real,
dimension(:),
intent(in) :: t
83 real,
dimension(:),
intent(in) :: s
84 real,
dimension(:),
intent(in) :: pressure
85 real,
dimension(:),
intent(out) :: rho
86 integer,
intent(in) :: start
87 integer,
intent(in) :: npts
88 real,
intent(in) :: rho_t0_s0
89 real,
intent(in) :: drho_dt
91 real,
intent(in) :: drho_ds
93 real,
optional,
intent(in) :: rho_ref
97 if (
present(rho_ref))
then ;
do j=start,start+npts-1
98 rho(j) = (rho_t0_s0 - rho_ref) + (drho_dt*t(j) + drho_ds*s(j))
99 enddo ;
else ;
do j=start,start+npts-1
100 rho(j) = rho_t0_s0 + drho_dt*t(j) + drho_ds*s(j)
110 Rho_T0_S0, dRho_dT, dRho_dS, spv_ref)
111 real,
intent(in) :: T
113 real,
intent(in) :: S
114 real,
intent(in) :: pressure
115 real,
intent(out) :: specvol
116 real,
intent(in) :: Rho_T0_S0
117 real,
intent(in) :: dRho_dT
118 real,
intent(in) :: dRho_dS
119 real,
optional,
intent(in) :: spv_ref
123 if (
present(spv_ref))
then
124 specvol = ((1.0 - rho_t0_s0*spv_ref) + spv_ref*(drho_dt*t + drho_ds*s)) / &
125 ( rho_t0_s0 + (drho_dt*t + drho_ds*s))
127 specvol = 1.0 / ( rho_t0_s0 + (drho_dt*t + drho_ds*s))
137 Rho_T0_S0, dRho_dT, dRho_dS, spv_ref)
138 real,
dimension(:),
intent(in) :: T
140 real,
dimension(:),
intent(in) :: S
141 real,
dimension(:),
intent(in) :: pressure
142 real,
dimension(:),
intent(out) :: specvol
143 integer,
intent(in) :: start
144 integer,
intent(in) :: npts
145 real,
intent(in) :: Rho_T0_S0
146 real,
intent(in) :: dRho_dT
147 real,
intent(in) :: dRho_dS
148 real,
optional,
intent(in) :: spv_ref
152 if (
present(spv_ref))
then ;
do j=start,start+npts-1
153 specvol(j) = ((1.0 - rho_t0_s0*spv_ref) + spv_ref*(drho_dt*t(j) + drho_ds*s(j))) / &
154 ( rho_t0_s0 + (drho_dt*t(j) + drho_ds*s(j)))
155 enddo ;
else ;
do j=start,start+npts-1
156 specvol(j) = 1.0 / ( rho_t0_s0 + (drho_dt*t(j) + drho_ds*s(j)))
164 drho_dS_out, Rho_T0_S0, dRho_dT, dRho_dS, start, npts)
165 real,
intent(in),
dimension(:) :: T
167 real,
intent(in),
dimension(:) :: S
168 real,
intent(in),
dimension(:) :: pressure
169 real,
intent(out),
dimension(:) :: drho_dT_out
171 real,
intent(out),
dimension(:) :: drho_dS_out
173 real,
intent(in) :: Rho_T0_S0
174 real,
intent(in) :: dRho_dT
175 real,
intent(in) :: dRho_dS
176 integer,
intent(in) :: start
177 integer,
intent(in) :: npts
181 do j=start,start+npts-1
182 drho_dt_out(j) = drho_dt
183 drho_ds_out(j) = drho_ds
191 drho_dS_out, Rho_T0_S0, dRho_dT, dRho_dS)
192 real,
intent(in) :: t
194 real,
intent(in) :: s
195 real,
intent(in) :: pressure
196 real,
intent(out) :: drho_dt_out
198 real,
intent(out) :: drho_ds_out
200 real,
intent(in) :: rho_t0_s0
201 real,
intent(in) :: drho_dt
202 real,
intent(in) :: drho_ds
203 drho_dt_out = drho_dt
204 drho_ds_out = drho_ds
211 drho_dT_dT, drho_dS_dP, drho_dT_dP)
212 real,
intent(in) :: T
213 real,
intent(in) :: S
214 real,
intent(in) :: pressure
215 real,
intent(out) :: drho_dS_dS
217 real,
intent(out) :: drho_dS_dT
219 real,
intent(out) :: drho_dT_dT
221 real,
intent(out) :: drho_dS_dP
223 real,
intent(out) :: drho_dT_dP
237 drho_dS_dP, drho_dT_dP, start, npts)
238 real,
dimension(:),
intent(in) :: T
239 real,
dimension(:),
intent(in) :: S
240 real,
dimension(:),
intent(in) :: pressure
241 real,
dimension(:),
intent(out) :: drho_dS_dS
243 real,
dimension(:),
intent(out) :: drho_dS_dT
245 real,
dimension(:),
intent(out) :: drho_dT_dT
247 real,
dimension(:),
intent(out) :: drho_dS_dP
249 real,
dimension(:),
intent(out) :: drho_dT_dP
251 integer,
intent(in) :: start
252 integer,
intent(in) :: npts
255 do j=start,start+npts-1
267 start, npts, Rho_T0_S0, dRho_dT, dRho_dS)
268 real,
intent(in),
dimension(:) :: t
270 real,
intent(in),
dimension(:) :: s
271 real,
intent(in),
dimension(:) :: pressure
272 real,
intent(out),
dimension(:) :: dsv_ds
274 real,
intent(out),
dimension(:) :: dsv_dt
276 integer,
intent(in) :: start
277 integer,
intent(in) :: npts
278 real,
intent(in) :: rho_t0_s0
279 real,
intent(in) :: drho_dt
281 real,
intent(in) :: drho_ds
287 do j=start,start+npts-1
289 i_rho2 = 1.0 / (rho_t0_s0 + (drho_dt*t(j) + drho_ds*s(j)))**2
290 dsv_dt(j) = -drho_dt * i_rho2
291 dsv_ds(j) = -drho_ds * i_rho2
300 Rho_T0_S0, dRho_dT, dRho_dS)
301 real,
intent(in),
dimension(:) :: t
303 real,
intent(in),
dimension(:) :: s
304 real,
intent(in),
dimension(:) :: pressure
305 real,
intent(out),
dimension(:) :: rho
306 real,
intent(out),
dimension(:) :: drho_dp
309 integer,
intent(in) :: start
310 integer,
intent(in) :: npts
311 real,
intent(in) :: rho_t0_s0
312 real,
intent(in) :: drho_dt
314 real,
intent(in) :: drho_ds
319 do j=start,start+npts-1
320 rho(j) = rho_t0_s0 + drho_dt*t(j) + drho_ds*s(j)
329 Rho_T0_S0, dRho_dT, dRho_dS, dpa, intz_dpa, intx_dpa, inty_dpa, &
330 bathyT, dz_neglect, useMassWghtInterp)
333 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
336 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
338 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
340 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
342 real,
intent(in) :: rho_ref
345 real,
intent(in) :: rho_0_pres
349 real,
intent(in) :: g_e
350 real,
intent(in) :: rho_t0_s0
351 real,
intent(in) :: drho_dt
353 real,
intent(in) :: drho_ds
355 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
358 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
359 optional,
intent(out) :: intz_dpa
362 real,
dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
363 optional,
intent(out) :: intx_dpa
366 real,
dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
367 optional,
intent(out) :: inty_dpa
370 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
371 optional,
intent(in) :: bathyt
372 real,
optional,
intent(in) :: dz_neglect
373 logical,
optional,
intent(in) :: usemasswghtinterp
382 real :: hwt_ll, hwt_lr
383 real :: hwt_rl, hwt_rr
388 logical :: do_massweight
389 real,
parameter :: c1_6 = 1.0/6.0, c1_90 = 1.0/90.0
390 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, ioff, joff, m
392 ioff = hio%idg_offset - hii%idg_offset
393 joff = hio%jdg_offset - hii%jdg_offset
397 isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
398 jsq = hio%JscB + joff ; jeq = hio%JecB + joff
399 is = hio%isc + ioff ; ie = hio%iec + ioff
400 js = hio%jsc + joff ; je = hio%jec + joff
402 do_massweight = .false.
403 if (
present(usemasswghtinterp))
then ;
if (usemasswghtinterp)
then
404 do_massweight = .true.
411 do j=jsq,jeq+1 ;
do i=isq,ieq+1
412 dz = z_t(i,j) - z_b(i,j)
413 rho_anom = (rho_t0_s0 - rho_ref) + drho_dt*t(i,j) + drho_ds*s(i,j)
414 dpa(i-ioff,j-joff) = g_e*rho_anom*dz
415 if (
present(intz_dpa)) intz_dpa(i-ioff,j-joff) = 0.5*g_e*rho_anom*dz**2
418 if (
present(intx_dpa))
then ;
do j=js,je ;
do i=isq,ieq
424 hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
426 if (hwght <= 0.0)
then
427 dzl = z_t(i,j) - z_b(i,j) ; dzr = z_t(i+1,j) - z_b(i+1,j)
428 ral = (rho_t0_s0 - rho_ref) + (drho_dt*t(i,j) + drho_ds*s(i,j))
429 rar = (rho_t0_s0 - rho_ref) + (drho_dt*t(i+1,j) + drho_ds*s(i+1,j))
431 intx_dpa(i-ioff,j-joff) = g_e*c1_6 * (dzl*(2.0*ral + rar) + dzr*(2.0*rar + ral))
433 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
434 hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
435 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
436 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
437 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
438 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
440 intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
442 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
443 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
445 dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i+1,j) - z_b(i+1,j))
446 rho_anom = (rho_t0_s0 - rho_ref) + &
447 (drho_dt * (wtt_l*t(i,j) + wtt_r*t(i+1,j)) + &
448 drho_ds * (wtt_l*s(i,j) + wtt_r*s(i+1,j)))
449 intz(m) = g_e*rho_anom*dz
452 intx_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
455 enddo ;
enddo ;
endif
457 if (
present(inty_dpa))
then ;
do j=jsq,jeq ;
do i=is,ie
463 hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
465 if (hwght <= 0.0)
then
466 dzl = z_t(i,j) - z_b(i,j) ; dzr = z_t(i,j+1) - z_b(i,j+1)
467 ral = (rho_t0_s0 - rho_ref) + (drho_dt*t(i,j) + drho_ds*s(i,j))
468 rar = (rho_t0_s0 - rho_ref) + (drho_dt*t(i,j+1) + drho_ds*s(i,j+1))
470 inty_dpa(i-ioff,j-joff) = g_e*c1_6 * (dzl*(2.0*ral + rar) + dzr*(2.0*rar + ral))
472 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
473 hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
474 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
475 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
476 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
477 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
479 intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
481 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
482 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
484 dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i,j+1) - z_b(i,j+1))
485 rho_anom = (rho_t0_s0 - rho_ref) + &
486 (drho_dt * (wtt_l*t(i,j) + wtt_r*t(i,j+1)) + &
487 drho_ds * (wtt_l*s(i,j) + wtt_r*s(i,j+1)))
488 intz(m) = g_e*rho_anom*dz
491 inty_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
495 enddo ;
enddo ;
endif
503 dRho_dT, dRho_dS, dza, intp_dza, intx_dza, inty_dza, halo_size, &
504 bathyP, dP_neglect, useMassWghtInterp)
506 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
509 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
511 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
513 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
515 real,
intent(in) :: alpha_ref
519 real,
intent(in) :: rho_t0_s0
520 real,
intent(in) :: drho_dt
522 real,
intent(in) :: drho_ds
524 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
527 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
528 optional,
intent(out) :: intp_dza
531 real,
dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
532 optional,
intent(out) :: intx_dza
536 real,
dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
537 optional,
intent(out) :: inty_dza
541 integer,
optional,
intent(in) :: halo_size
542 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
543 optional,
intent(in) :: bathyp
544 real,
optional,
intent(in) :: dp_neglect
546 logical,
optional,
intent(in) :: usemasswghtinterp
556 real :: hwt_ll, hwt_lr
557 real :: hwt_rl, hwt_rr
562 logical :: do_massweight
563 real,
parameter :: c1_6 = 1.0/6.0, c1_90 = 1.0/90.0
564 integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, halo
566 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
567 halo = 0 ;
if (
present(halo_size)) halo = max(halo_size,0)
568 ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
569 if (
present(intx_dza))
then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh);
endif
570 if (
present(inty_dza))
then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh);
endif
572 do_massweight = .false.
573 if (
present(usemasswghtinterp))
then ;
if (usemasswghtinterp)
then
574 do_massweight = .true.
581 do j=jsh,jeh ;
do i=ish,ieh
582 dp = p_b(i,j) - p_t(i,j)
583 drho_ts = drho_dt*t(i,j) + drho_ds*s(i,j)
585 alpha_anom = ((1.0-rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
586 dza(i,j) = alpha_anom*dp
587 if (
present(intp_dza)) intp_dza(i,j) = 0.5*alpha_anom*dp**2
590 if (
present(intx_dza))
then ;
do j=hi%jsc,hi%jec ;
do i=isq,ieq
596 hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
598 if (hwght <= 0.0)
then
599 dpl = p_b(i,j) - p_t(i,j) ; dpr = p_b(i+1,j) - p_t(i+1,j)
600 drho_ts = drho_dt*t(i,j) + drho_ds*s(i,j)
601 aal = ((1.0 - rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
602 drho_ts = drho_dt*t(i+1,j) + drho_ds*s(i+1,j)
603 aar = ((1.0 - rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
605 intx_dza(i,j) = c1_6 * (2.0*(dpl*aal + dpr*aar) + (dpl*aar + dpr*aal))
607 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
608 hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
609 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
610 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
611 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
612 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
614 intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
616 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
617 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
621 dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i+1,j) - p_t(i+1,j))
623 drho_ts = drho_dt*(wtt_l*t(i,j) + wtt_r*t(i+1,j)) + &
624 drho_ds*(wtt_l*s(i,j) + wtt_r*s(i+1,j))
626 alpha_anom = ((1.0-rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
627 intp(m) = alpha_anom*dp
630 intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
633 enddo ;
enddo ;
endif
635 if (
present(inty_dza))
then ;
do j=jsq,jeq ;
do i=hi%isc,hi%iec
641 hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
643 if (hwght <= 0.0)
then
644 dpl = p_b(i,j) - p_t(i,j) ; dpr = p_b(i,j+1) - p_t(i,j+1)
645 drho_ts = drho_dt*t(i,j) + drho_ds*s(i,j)
646 aal = ((1.0 - rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
647 drho_ts = drho_dt*t(i,j+1) + drho_ds*s(i,j+1)
648 aar = ((1.0 - rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
650 inty_dza(i,j) = c1_6 * (2.0*(dpl*aal + dpr*aar) + (dpl*aar + dpr*aal))
652 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
653 hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
654 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
655 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
656 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
657 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
659 intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
661 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
662 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
666 dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i,j+1) - p_t(i,j+1))
668 drho_ts = drho_dt*(wtt_l*t(i,j) + wtt_r*t(i,j+1)) + &
669 drho_ds*(wtt_l*s(i,j) + wtt_r*s(i,j+1))
671 alpha_anom = ((1.0-rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
672 intp(m) = alpha_anom*dp
675 inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
678 enddo ;
enddo ;
endif