35 implicit none ;
private
37 #include <MOM_memory.h>
51 public gsw_sp_from_sr, gsw_pt_from_ct
92 integer :: form_of_eos = 0
93 integer :: form_of_tfreeze = 0
95 logical :: eos_quadrature
97 logical :: compressible = .true.
139 real,
intent(in) :: T
140 real,
intent(in) :: S
141 real,
intent(in) :: pressure
142 real,
intent(out) :: rho
144 real,
optional,
intent(in) :: rho_ref
145 real,
optional,
intent(in) :: scale
148 if (.not.
associated(eos))
call mom_error(fatal, &
149 "calculate_density_scalar called with an unassociated EOS_type EOS.")
151 select case (eos%form_of_EOS)
154 eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, rho_ref)
165 "calculate_density_scalar: EOS is not valid.")
168 if (
present(scale))
then ;
if (scale /= 1.0)
then
177 real,
dimension(:),
intent(in) :: T
178 real,
dimension(:),
intent(in) :: S
179 real,
dimension(:),
intent(in) :: pressure
180 real,
dimension(:),
intent(out) :: rho
181 integer,
intent(in) :: start
182 integer,
intent(in) :: npts
184 real,
optional,
intent(in) :: rho_ref
185 real,
optional,
intent(in) :: scale
189 if (.not.
associated(eos))
call mom_error(fatal, &
190 "calculate_density_array called with an unassociated EOS_type EOS.")
192 select case (eos%form_of_EOS)
195 eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, rho_ref)
206 "calculate_density_array: EOS%form_of_EOS is not valid.")
209 if (
present(scale))
then ;
if (scale /= 1.0)
then
210 do j=start,start+npts-1 ; rho(j) = scale * rho(j) ;
enddo
218 real,
intent(in) :: T
219 real,
intent(in) :: S
220 real,
intent(in) :: pressure
221 real,
intent(out) :: specvol
223 real,
optional,
intent(in) :: spv_ref
224 real,
optional,
intent(in) :: scale
229 if (.not.
associated(eos))
call mom_error(fatal, &
230 "calculate_spec_vol_scalar called with an unassociated EOS_type EOS.")
232 select case (eos%form_of_EOS)
235 eos%rho_T0_S0, eos%drho_dT, eos%drho_dS, spv_ref)
244 if (
present(spv_ref))
then
245 specvol = 1.0 / rho - spv_ref
251 "calculate_spec_vol_scalar: EOS is not valid.")
254 if (
present(scale))
then ;
if (scale /= 1.0)
then
255 specvol = scale * specvol
264 real,
dimension(:),
intent(in) :: T
266 real,
dimension(:),
intent(in) :: S
267 real,
dimension(:),
intent(in) :: pressure
268 real,
dimension(:),
intent(out) :: specvol
269 integer,
intent(in) :: start
270 integer,
intent(in) :: npts
272 real,
optional,
intent(in) :: spv_ref
273 real,
optional,
intent(in) :: scale
276 real,
dimension(size(specvol)) :: rho
279 if (.not.
associated(eos))
call mom_error(fatal, &
280 "calculate_spec_vol_array called with an unassociated EOS_type EOS.")
282 select case (eos%form_of_EOS)
285 eos%rho_T0_S0, eos%drho_dT, eos%drho_dS, spv_ref)
294 if (
present(spv_ref))
then
295 specvol(:) = 1.0 / rho(:) - spv_ref
297 specvol(:) = 1.0 / rho(:)
301 "calculate_spec_vol_array: EOS%form_of_EOS is not valid.")
304 if (
present(scale))
then ;
if (scale /= 1.0)
then ;
do j=start,start+npts-1
305 specvol(j) = scale * specvol(j)
306 enddo ;
endif ;
endif
313 real,
intent(in) :: S
314 real,
intent(in) :: pressure
315 real,
intent(out) :: T_fr
319 if (.not.
associated(eos))
call mom_error(fatal, &
320 "calculate_TFreeze_scalar called with an unassociated EOS_type EOS.")
322 select case (eos%form_of_TFreeze)
325 eos%dTFr_dS, eos%dTFr_dp)
332 "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
339 real,
dimension(:),
intent(in) :: S
340 real,
dimension(:),
intent(in) :: pressure
341 real,
dimension(:),
intent(out) :: T_fr
343 integer,
intent(in) :: start
344 integer,
intent(in) :: npts
347 if (.not.
associated(eos))
call mom_error(fatal, &
348 "calculate_TFreeze_scalar called with an unassociated EOS_type EOS.")
350 select case (eos%form_of_TFreeze)
353 eos%TFr_S0_P0, eos%dTFr_dS, eos%dTFr_dp)
360 "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
367 real,
dimension(:),
intent(in) :: T
368 real,
dimension(:),
intent(in) :: S
369 real,
dimension(:),
intent(in) :: pressure
370 real,
dimension(:),
intent(out) :: drho_dT
372 real,
dimension(:),
intent(out) :: drho_dS
374 integer,
intent(in) :: start
375 integer,
intent(in) :: npts
377 real,
optional,
intent(in) :: scale
381 if (.not.
associated(eos))
call mom_error(fatal, &
382 "calculate_density_derivs called with an unassociated EOS_type EOS.")
384 select case (eos%form_of_EOS)
387 eos%dRho_dT, eos%dRho_dS, start, npts)
389 call calculate_density_derivs_unesco(t, s, pressure, drho_dt, drho_ds, start, npts)
398 "calculate_density_derivs_array: EOS%form_of_EOS is not valid.")
401 if (
present(scale))
then ;
if (scale /= 1.0)
then ;
do j=start,start+npts-1
402 drho_dt(j) = scale * drho_dt(j)
403 drho_ds(j) = scale * drho_ds(j)
404 enddo ;
endif ;
endif
411 real,
intent(in) :: T
412 real,
intent(in) :: S
413 real,
intent(in) :: pressure
414 real,
intent(out) :: drho_dT
416 real,
intent(out) :: drho_dS
419 real,
optional,
intent(in) :: scale
422 if (.not.
associated(eos))
call mom_error(fatal, &
423 "calculate_density_derivs called with an unassociated EOS_type EOS.")
425 select case (eos%form_of_EOS)
428 eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
435 "calculate_density_derivs_scalar: EOS%form_of_EOS is not valid.")
438 if (
present(scale))
then ;
if (scale /= 1.0)
then
439 drho_dt = scale * drho_dt
440 drho_ds = scale * drho_ds
447 drho_dS_dP, drho_dT_dP, start, npts, EOS, scale)
448 real,
dimension(:),
intent(in) :: T
449 real,
dimension(:),
intent(in) :: S
450 real,
dimension(:),
intent(in) :: pressure
451 real,
dimension(:),
intent(out) :: drho_dS_dS
453 real,
dimension(:),
intent(out) :: drho_dS_dT
455 real,
dimension(:),
intent(out) :: drho_dT_dT
457 real,
dimension(:),
intent(out) :: drho_dS_dP
459 real,
dimension(:),
intent(out) :: drho_dT_dP
461 integer,
intent(in) :: start
462 integer,
intent(in) :: npts
464 real,
optional,
intent(in) :: scale
468 if (.not.
associated(eos))
call mom_error(fatal, &
469 "calculate_density_derivs called with an unassociated EOS_type EOS.")
471 select case (eos%form_of_EOS)
474 drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
477 drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
480 drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
483 "calculate_density_derivs: EOS%form_of_EOS is not valid.")
486 if (
present(scale))
then ;
if (scale /= 1.0)
then ;
do j=start,start+npts-1
487 drho_ds_ds(j) = scale * drho_ds_ds(j)
488 drho_ds_dt(j) = scale * drho_ds_dt(j)
489 drho_dt_dt(j) = scale * drho_dt_dt(j)
490 drho_ds_dp(j) = scale * drho_ds_dp(j)
491 drho_dt_dp(j) = scale * drho_dt_dp(j)
492 enddo ;
endif ;
endif
498 drho_dS_dP, drho_dT_dP, EOS, scale)
499 real,
intent(in) :: T
500 real,
intent(in) :: S
501 real,
intent(in) :: pressure
502 real,
intent(out) :: drho_dS_dS
504 real,
intent(out) :: drho_dS_dT
506 real,
intent(out) :: drho_dT_dT
508 real,
intent(out) :: drho_dS_dP
510 real,
intent(out) :: drho_dT_dP
513 real,
optional,
intent(in) :: scale
516 if (.not.
associated(eos))
call mom_error(fatal, &
517 "calculate_density_derivs called with an unassociated EOS_type EOS.")
519 select case (eos%form_of_EOS)
522 drho_dt_dt, drho_ds_dp, drho_dt_dp)
525 drho_dt_dt, drho_ds_dp, drho_dt_dp)
528 drho_dt_dt, drho_ds_dp, drho_dt_dp)
531 "calculate_density_derivs: EOS%form_of_EOS is not valid.")
534 if (
present(scale))
then ;
if (scale /= 1.0)
then
535 drho_ds_ds = scale * drho_ds_ds
536 drho_ds_dt = scale * drho_ds_dt
537 drho_dt_dt = scale * drho_dt_dt
538 drho_ds_dp = scale * drho_ds_dp
539 drho_dt_dp = scale * drho_dt_dp
546 real,
dimension(:),
intent(in) :: t
547 real,
dimension(:),
intent(in) :: s
548 real,
dimension(:),
intent(in) :: pressure
549 real,
dimension(:),
intent(out) :: dsv_dt
551 real,
dimension(:),
intent(out) :: dsv_ds
553 integer,
intent(in) :: start
554 integer,
intent(in) :: npts
556 real,
optional,
intent(in) :: scale
560 real,
dimension(size(T)) :: drho_dt, drho_ds, rho
563 if (.not.
associated(eos))
call mom_error(fatal, &
564 "calculate_density_derivs called with an unassociated EOS_type EOS.")
566 select case (eos%form_of_EOS)
568 call calculate_specvol_derivs_linear(t, s, pressure, dsv_dt, dsv_ds, start, &
569 npts, eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
572 call calculate_density_derivs_unesco(t, s, pressure, drho_dt, drho_ds, start, npts)
573 do j=start,start+npts-1
574 dsv_dt(j) = -drho_dt(j)/(rho(j)**2)
575 dsv_ds(j) = -drho_ds(j)/(rho(j)**2)
578 call calculate_specvol_derivs_wright(t, s, pressure, dsv_dt, dsv_ds, start, npts)
580 call calculate_specvol_derivs_teos10(t, s, pressure, dsv_dt, dsv_ds, start, npts)
584 do j=start,start+npts-1
585 dsv_dt(j) = -drho_dt(j)/(rho(j)**2)
586 dsv_ds(j) = -drho_ds(j)/(rho(j)**2)
590 "calculate_density_derivs: EOS%form_of_EOS is not valid.")
593 if (
present(scale))
then ;
if (scale /= 1.0)
then ;
do j=start,start+npts-1
594 dsv_dt(j) = scale * dsv_dt(j)
595 dsv_ds(j) = scale * dsv_ds(j)
596 enddo ;
endif ;
endif
603 real,
dimension(:),
intent(in) :: T
604 real,
dimension(:),
intent(in) :: S
605 real,
dimension(:),
intent(in) :: pressure
606 real,
dimension(:),
intent(out) :: rho
607 real,
dimension(:),
intent(out) :: drho_dp
609 integer,
intent(in) :: start
610 integer,
intent(in) :: npts
613 if (.not.
associated(eos))
call mom_error(fatal, &
614 "calculate_compress called with an unassociated EOS_type EOS.")
616 select case (eos%form_of_EOS)
619 eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
623 call calculate_compress_wright(t, s, pressure, rho, drho_dp, start, npts)
625 call calculate_compress_teos10(t, s, pressure, rho, drho_dp, start, npts)
630 "calculate_compress: EOS%form_of_EOS is not valid.")
638 real,
intent(in) :: T
639 real,
intent(in) :: S
640 real,
intent(in) :: pressure
641 real,
intent(out) :: rho
642 real,
intent(out) :: drho_dp
646 real,
dimension(1) :: Ta, Sa, pa, rhoa, drho_dpa
648 if (.not.
associated(eos))
call mom_error(fatal, &
649 "calculate_compress called with an unassociated EOS_type EOS.")
650 ta(1) = t ; sa(1) = s; pa(1) = pressure
653 rho = rhoa(1) ; drho_dp = drho_dpa(1)
663 dza, intp_dza, intx_dza, inty_dza, halo_size, &
664 bathyP, dP_tiny, useMassWghtInterp)
666 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
668 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
670 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
672 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
674 real,
intent(in) :: alpha_ref
679 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
682 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
683 optional,
intent(out) :: intp_dza
686 real,
dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
687 optional,
intent(out) :: intx_dza
690 real,
dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
691 optional,
intent(out) :: inty_dza
694 integer,
optional,
intent(in) :: halo_size
695 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
696 optional,
intent(in) :: bathyp
697 real,
optional,
intent(in) :: dp_tiny
699 logical,
optional,
intent(in) :: usemasswghtinterp
702 if (.not.
associated(eos))
call mom_error(fatal, &
703 "int_specific_vol_dp called with an unassociated EOS_type EOS.")
705 if (eos%EOS_quadrature)
then
707 dza, intp_dza, intx_dza, inty_dza, halo_size, &
708 bathyp, dp_tiny, usemasswghtinterp)
709 else ;
select case (eos%form_of_EOS)
712 eos%dRho_dT, eos%dRho_dS, dza, intp_dza, &
713 intx_dza, inty_dza, halo_size, &
714 bathyp, dp_tiny, usemasswghtinterp)
716 call int_spec_vol_dp_wright(t, s, p_t, p_b, alpha_ref, hi, dza, &
717 intp_dza, intx_dza, inty_dza, halo_size, &
718 bathyp, dp_tiny, usemasswghtinterp)
721 dza, intp_dza, intx_dza, inty_dza, halo_size, &
722 bathyp, dp_tiny, usemasswghtinterp)
730 subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, &
731 dpa, intz_dpa, intx_dpa, inty_dpa, &
732 bathyT, dz_neglect, useMassWghtInterp)
735 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
737 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
739 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
741 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
743 real,
intent(in) :: rho_ref
745 real,
intent(in) :: rho_0
747 real,
intent(in) :: g_e
749 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
751 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
752 optional,
intent(out) :: intz_dpa
755 real,
dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
756 optional,
intent(out) :: intx_dpa
759 real,
dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
760 optional,
intent(out) :: inty_dpa
763 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
764 optional,
intent(in) :: bathyt
765 real,
optional,
intent(in) :: dz_neglect
766 logical,
optional,
intent(in) :: usemasswghtinterp
769 if (.not.
associated(eos))
call mom_error(fatal, &
770 "int_density_dz called with an unassociated EOS_type EOS.")
772 if (eos%EOS_quadrature)
then
774 eos, dpa, intz_dpa, intx_dpa, inty_dpa, &
775 bathyt, dz_neglect, usemasswghtinterp)
776 else ;
select case (eos%form_of_EOS)
778 call int_density_dz_linear(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio, &
779 eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, &
780 dpa, intz_dpa, intx_dpa, inty_dpa, &
781 bathyt, dz_neglect, usemasswghtinterp)
783 call int_density_dz_wright(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio, &
784 dpa, intz_dpa, intx_dpa, inty_dpa, &
785 bathyt, dz_neglect, usemasswghtinterp)
788 eos, dpa, intz_dpa, intx_dpa, inty_dpa, &
789 bathyt, dz_neglect, usemasswghtinterp)
798 if (.not.
associated(eos))
call mom_error(fatal, &
799 "query_compressible called with an unassociated EOS_type EOS.")
805 subroutine eos_init(param_file, EOS)
809 #include "version_variable.h"
810 character(len=40) :: mdl =
"MOM_EOS"
811 character(len=40) :: tmpstr
818 call get_param(param_file, mdl,
"EQN_OF_STATE", tmpstr, &
819 "EQN_OF_STATE determines which ocean equation of state "//&
820 "should be used. Currently, the valid choices are "//&
821 '"LINEAR", "UNESCO", "WRIGHT", "NEMO" and "TEOS10". '//&
822 "This is only used if USE_EOS is true.", default=
eos_default)
835 call mom_error(fatal,
"interpret_eos_selection: EQN_OF_STATE "//&
836 trim(tmpstr) //
"in input file is invalid.")
838 call mom_mesg(
'interpret_eos_selection: equation of state set to "' // &
839 trim(tmpstr)//
'"', 5)
842 eos%Compressible = .false.
843 call get_param(param_file, mdl,
"RHO_T0_S0", eos%Rho_T0_S0, &
845 "this is the density at T=0, S=0.", units=
"kg m-3", &
847 call get_param(param_file, mdl,
"DRHO_DT", eos%dRho_dT, &
849 "this is the partial derivative of density with "//&
850 "temperature.", units=
"kg m-3 K-1", default=-0.2)
851 call get_param(param_file, mdl,
"DRHO_DS", eos%dRho_dS, &
853 "this is the partial derivative of density with "//&
854 "salinity.", units=
"kg m-3 PSU-1", default=0.8)
857 call get_param(param_file, mdl,
"EOS_QUADRATURE", eos%EOS_quadrature, &
858 "If true, always use the generic (quadrature) code "//&
859 "code for the integrals of density.", default=.false.)
861 call get_param(param_file, mdl,
"TFREEZE_FORM", tmpstr, &
862 "TFREEZE_FORM determines which expression should be "//&
863 "used for the freezing point. Currently, the valid "//&
864 'choices are "LINEAR", "MILLERO_78", "TEOS10"', &
874 call mom_error(fatal,
"interpret_eos_selection: TFREEZE_FORM "//&
875 trim(tmpstr) //
"in input file is invalid.")
879 call get_param(param_file, mdl,
"TFREEZE_S0_P0",eos%TFr_S0_P0, &
881 "this is the freezing potential temperature at "//&
882 "S=0, P=0.", units=
"deg C", default=0.0)
883 call get_param(param_file, mdl,
"DTFREEZE_DS",eos%dTFr_dS, &
885 "this is the derivative of the freezing potential "//&
886 "temperature with salinity.", &
887 units=
"deg C PSU-1", default=-0.054)
888 call get_param(param_file, mdl,
"DTFREEZE_DP",eos%dTFr_dP, &
890 "this is the derivative of the freezing potential "//&
891 "temperature with pressure.", &
892 units=
"deg C Pa-1", default=0.0)
897 call mom_error(fatal,
"interpret_eos_selection: EOS_TEOS10 or EOS_NEMO \n" //&
898 "should only be used along with TFREEZE_FORM = TFREEZE_TEOS10 .")
905 subroutine eos_manual_init(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, &
906 Rho_T0_S0, drho_dT, dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp)
908 integer,
optional,
intent(in) :: form_of_eos
909 integer,
optional,
intent(in) :: form_of_tfreeze
911 logical,
optional,
intent(in) :: eos_quadrature
913 logical,
optional,
intent(in) :: compressible
914 real ,
optional,
intent(in) :: rho_t0_s0
915 real ,
optional,
intent(in) :: drho_dt
917 real ,
optional,
intent(in) :: drho_ds
919 real ,
optional,
intent(in) :: tfr_s0_p0
920 real ,
optional,
intent(in) :: dtfr_ds
922 real ,
optional,
intent(in) :: dtfr_dp
925 if (
present(form_of_eos )) eos%form_of_EOS = form_of_eos
926 if (
present(form_of_tfreeze)) eos%form_of_TFreeze = form_of_tfreeze
927 if (
present(eos_quadrature )) eos%EOS_quadrature = eos_quadrature
928 if (
present(compressible )) eos%Compressible = compressible
929 if (
present(rho_t0_s0 )) eos%Rho_T0_S0 = rho_t0_s0
930 if (
present(drho_dt )) eos%drho_dT = drho_dt
931 if (
present(drho_ds )) eos%dRho_dS = drho_ds
932 if (
present(tfr_s0_p0 )) eos%TFr_S0_P0 = tfr_s0_p0
933 if (
present(dtfr_ds )) eos%dTFr_dS = dtfr_ds
934 if (
present(dtfr_dp )) eos%dTFr_dp = dtfr_dp
942 if (.not.
associated(eos))
allocate(eos)
949 if (
associated(eos))
deallocate(eos)
957 subroutine eos_use_linear(Rho_T0_S0, dRho_dT, dRho_dS, EOS, use_quadrature)
958 real,
intent(in) :: rho_t0_s0
959 real,
intent(in) :: drho_dt
960 real,
intent(in) :: drho_ds
961 logical,
optional,
intent(in) :: use_quadrature
965 if (.not.
associated(eos))
call mom_error(fatal, &
966 "MOM_EOS.F90: EOS_use_linear() called with an unassociated EOS_type EOS.")
969 eos%Compressible = .false.
970 eos%Rho_T0_S0 = rho_t0_s0
971 eos%dRho_dT = drho_dt
972 eos%dRho_dS = drho_ds
973 eos%EOS_quadrature = .false.
974 if (
present(use_quadrature)) eos%EOS_quadrature = use_quadrature
982 EOS, dpa, intz_dpa, intx_dpa, inty_dpa, &
983 bathyT, dz_neglect, useMassWghtInterp)
986 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
988 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
990 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
992 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
994 real,
intent(in) :: rho_ref
997 real,
intent(in) :: rho_0
1000 real,
intent(in) :: g_e
1002 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1005 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1006 optional,
intent(out) :: intz_dpa
1009 real,
dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
1010 optional,
intent(out) :: intx_dpa
1013 real,
dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
1014 optional,
intent(out) :: inty_dpa
1017 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1018 optional,
intent(in) :: bathyt
1019 real,
optional,
intent(in) :: dz_neglect
1020 logical,
optional,
intent(in) :: usemasswghtinterp
1022 real :: t5(5), s5(5), p5(5), r5(5)
1024 real :: w_left, w_right
1025 real,
parameter :: c1_90 = 1.0/90.0
1026 real :: gxrho, i_rho
1031 real :: hwt_ll, hwt_lr
1032 real :: hwt_rl, hwt_rr
1034 real :: wtt_l, wtt_r
1037 logical :: do_massweight
1038 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m, n, ioff, joff
1040 ioff = hio%idg_offset - hii%idg_offset
1041 joff = hio%jdg_offset - hii%jdg_offset
1045 isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
1046 jsq = hio%JscB + joff ; jeq = hio%JecB + joff
1047 is = hio%isc + ioff ; ie = hio%iec + ioff
1048 js = hio%jsc + joff ; je = hio%jec + joff
1053 do_massweight = .false.
1054 if (
present(usemasswghtinterp))
then ;
if (usemasswghtinterp)
then
1055 do_massweight = .true.
1056 if (.not.
present(bathyt))
call mom_error(fatal,
"int_density_dz_generic: "//&
1057 "bathyT must be present if useMassWghtInterp is present and true.")
1058 if (.not.
present(dz_neglect))
call mom_error(fatal,
"int_density_dz_generic: "//&
1059 "dz_neglect must be present if useMassWghtInterp is present and true.")
1062 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1063 dz = z_t(i,j) - z_b(i,j)
1065 t5(n) = t(i,j) ; s5(n) = s(i,j)
1066 p5(n) = -gxrho*(z_t(i,j) - 0.25*real(n-1)*dz)
1071 rho_anom = c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))
1072 dpa(i-ioff,j-joff) = g_e*dz*rho_anom
1075 if (
present(intz_dpa)) intz_dpa(i-ioff,j-joff) = 0.5*g_e*dz**2 * &
1076 (rho_anom - c1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) )
1079 if (
present(intx_dpa))
then ;
do j=js,je ;
do i=isq,ieq
1084 if (do_massweight) &
1085 hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
1086 if (hwght > 0.)
then
1087 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
1088 hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
1089 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1090 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
1091 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
1092 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
1094 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
1097 intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
1101 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
1102 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
1103 dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i+1,j) - z_b(i+1,j))
1104 t5(1) = wtt_l*t(i,j) + wtt_r*t(i+1,j)
1105 s5(1) = wtt_l*s(i,j) + wtt_r*s(i+1,j)
1106 p5(1) = -gxrho*(wt_l*z_t(i,j) + wt_r*z_t(i+1,j))
1108 t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) + gxrho*0.25*dz
1113 intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)))
1116 intx_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
1118 enddo ;
enddo ;
endif
1120 if (
present(inty_dpa))
then ;
do j=jsq,jeq ;
do i=is,ie
1125 if (do_massweight) &
1126 hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
1127 if (hwght > 0.)
then
1128 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
1129 hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
1130 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1131 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
1132 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
1133 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
1135 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
1138 intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i-ioff,j-joff+1)
1142 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
1143 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
1144 dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i,j+1) - z_b(i,j+1))
1145 t5(1) = wtt_l*t(i,j) + wtt_r*t(i,j+1)
1146 s5(1) = wtt_l*s(i,j) + wtt_r*s(i,j+1)
1147 p5(1) = -gxrho*(wt_l*z_t(i,j) + wt_r*z_t(i,j+1))
1149 t5(n) = t5(1) ; s5(n) = s5(1)
1150 p5(n) = p5(n-1) + gxrho*0.25*dz
1155 intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)))
1158 inty_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
1160 enddo ;
enddo ;
endif
1168 rho_0, G_e, dz_subroundoff, bathyT, HII, HIO, EOS, dpa, &
1169 intz_dpa, intx_dpa, inty_dpa, &
1173 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1175 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1177 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1179 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1181 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1184 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1186 real,
intent(in) :: rho_ref
1188 real,
intent(in) :: rho_0
1190 real,
intent(in) :: g_e
1191 real,
intent(in) :: dz_subroundoff
1192 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1193 intent(in) :: bathyt
1195 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1197 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1198 optional,
intent(out) :: intz_dpa
1201 real,
dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
1202 optional,
intent(out) :: intx_dpa
1205 real,
dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
1206 optional,
intent(out) :: inty_dpa
1209 logical,
optional,
intent(in) :: usemasswghtinterp
1223 real :: t5((5*hio%iscb+1):(5*(hio%iecb+2)))
1224 real :: s5((5*hio%iscb+1):(5*(hio%iecb+2)))
1225 real :: p5((5*hio%iscb+1):(5*(hio%iecb+2)))
1226 real :: r5((5*hio%iscb+1):(5*(hio%iecb+2)))
1227 real :: t15((15*hio%iscb+1):(15*(hio%iecb+1)))
1228 real :: s15((15*hio%iscb+1):(15*(hio%iecb+1)))
1229 real :: p15((15*hio%iscb+1):(15*(hio%iecb+1)))
1230 real :: r15((15*hio%iscb+1):(15*(hio%iecb+1)))
1231 real :: wt_t(5), wt_b(5)
1233 real :: w_left, w_right
1236 real,
parameter :: c1_90 = 1.0/90.0
1239 real :: dz(hio%iscb:hio%iecb+1)
1240 real :: dz_x(5,hio%iscb:hio%iecb)
1241 real :: dz_y(5,hio%isc:hio%iec)
1242 real :: weight_t, weight_b
1243 real :: massweighttoggle
1244 real :: ttl, tbl, ttr, tbr
1245 real :: stl, sbl, str, sbr
1249 integer :: isq, ieq, jsq, jeq, i, j, m, n
1250 integer :: iin, jin, ioff, joff
1253 ioff = hio%idg_offset - hii%idg_offset
1254 joff = hio%jdg_offset - hii%jdg_offset
1256 isq = hio%IscB ; ieq = hio%IecB ; jsq = hio%JscB ; jeq = hio%JecB
1260 massweighttoggle = 0.
1261 if (
present(usemasswghtinterp))
then
1262 if (usemasswghtinterp) massweighttoggle = 1.
1266 wt_t(n) = 0.25 * real(5-n)
1267 wt_b(n) = 1.0 - wt_t(n)
1275 do i = isq,ieq+1 ; iin = i+ioff
1276 dz(i) = z_t(iin,jin) - z_b(iin,jin)
1278 p5(i*5+n) = -gxrho*(z_t(iin,jin) - 0.25*real(n-1)*dz(i))
1280 s5(i*5+n) = wt_t(n) * s_t(iin,jin) + wt_b(n) * s_b(iin,jin)
1281 t5(i*5+n) = wt_t(n) * t_t(iin,jin) + wt_b(n) * t_b(iin,jin)
1286 do i=isq,ieq+1 ; iin = i+ioff
1288 rho_anom = c1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3))
1289 dpa(i,j) = g_e*dz(i)*rho_anom
1290 if (
present(intz_dpa))
then
1293 intz_dpa(i,j) = 0.5*g_e*dz(i)**2 * &
1294 (rho_anom - c1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) )
1303 if (
present(intx_dpa))
then ;
do j=hio%jsc,hio%jec ; jin = j+joff
1304 do i=isq,ieq ; iin = i+ioff
1312 hwght = massweighttoggle * &
1313 max(0., -bathyt(iin,jin)-z_t(iin+1,jin), -bathyt(iin+1,jin)-z_t(iin,jin))
1314 if (hwght > 0.)
then
1315 hl = (z_t(iin,jin) - z_b(iin,jin)) + dz_subroundoff
1316 hr = (z_t(iin+1,jin) - z_b(iin+1,jin)) + dz_subroundoff
1317 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1318 idenom = 1./( hwght*(hr + hl) + hl*hr )
1319 ttl = ( (hwght*hr)*t_t(iin+1,jin) + (hwght*hl + hr*hl)*t_t(iin,jin) ) * idenom
1320 ttr = ( (hwght*hl)*t_t(iin,jin) + (hwght*hr + hr*hl)*t_t(iin+1,jin) ) * idenom
1321 tbl = ( (hwght*hr)*t_b(iin+1,jin) + (hwght*hl + hr*hl)*t_b(iin,jin) ) * idenom
1322 tbr = ( (hwght*hl)*t_b(iin,jin) + (hwght*hr + hr*hl)*t_b(iin+1,jin) ) * idenom
1323 stl = ( (hwght*hr)*s_t(iin+1,jin) + (hwght*hl + hr*hl)*s_t(iin,jin) ) * idenom
1324 str = ( (hwght*hl)*s_t(iin,jin) + (hwght*hr + hr*hl)*s_t(iin+1,jin) ) * idenom
1325 sbl = ( (hwght*hr)*s_b(iin+1,jin) + (hwght*hl + hr*hl)*s_b(iin,jin) ) * idenom
1326 sbr = ( (hwght*hl)*s_b(iin,jin) + (hwght*hr + hr*hl)*s_b(iin+1,jin) ) * idenom
1328 ttl = t_t(iin,jin); tbl = t_b(iin,jin); ttr = t_t(iin+1,jin); tbr = t_b(iin+1,jin)
1329 stl = s_t(iin,jin); sbl = s_b(iin,jin); str = s_t(iin+1,jin); sbr = s_b(iin+1,jin)
1333 w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
1334 dz_x(m,i) = w_left*(z_t(iin,jin) - z_b(iin,jin)) + w_right*(z_t(iin+1,jin) - z_b(iin+1,jin))
1341 t15(pos+1) = w_left*ttl + w_right*ttr
1342 t15(pos+5) = w_left*tbl + w_right*tbr
1344 s15(pos+1) = w_left*stl + w_right*str
1345 s15(pos+5) = w_left*sbl + w_right*sbr
1347 p15(pos+1) = -gxrho*(w_left*z_t(iin,jin) + w_right*z_t(iin+1,jin))
1351 p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_x(m,i)
1356 weight_t = 0.25 * real(5-n)
1357 weight_b = 1.0 - weight_t
1358 s15(pos+n) = weight_t * s15(pos+1) + weight_b * s15(pos+5)
1359 t15(pos+n) = weight_t * t15(pos+1) + weight_b * t15(pos+5)
1366 do i=isq,ieq ; iin = i+ioff
1367 intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
1372 intz(m) = g_e*dz_x(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + 32.0*(r15(pos+2)+r15(pos+4)) + &
1376 intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
1384 if (
present(inty_dpa))
then ;
do j=jsq,jeq ; jin = j+joff
1385 do i=hio%isc,hio%iec ; iin = i+ioff
1393 hwght = massweighttoggle * &
1394 max(0., -bathyt(i,j)-z_t(iin,jin+1), -bathyt(i,j+1)-z_t(iin,jin))
1395 if (hwght > 0.)
then
1396 hl = (z_t(iin,jin) - z_b(iin,jin)) + dz_subroundoff
1397 hr = (z_t(iin,jin+1) - z_b(iin,jin+1)) + dz_subroundoff
1398 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1399 idenom = 1./( hwght*(hr + hl) + hl*hr )
1400 ttl = ( (hwght*hr)*t_t(iin,jin+1) + (hwght*hl + hr*hl)*t_t(iin,jin) ) * idenom
1401 ttr = ( (hwght*hl)*t_t(iin,jin) + (hwght*hr + hr*hl)*t_t(iin,jin+1) ) * idenom
1402 tbl = ( (hwght*hr)*t_b(iin,jin+1) + (hwght*hl + hr*hl)*t_b(iin,jin) ) * idenom
1403 tbr = ( (hwght*hl)*t_b(iin,jin) + (hwght*hr + hr*hl)*t_b(iin,jin+1) ) * idenom
1404 stl = ( (hwght*hr)*s_t(iin,jin+1) + (hwght*hl + hr*hl)*s_t(iin,jin) ) * idenom
1405 str = ( (hwght*hl)*s_t(iin,jin) + (hwght*hr + hr*hl)*s_t(iin,jin+1) ) * idenom
1406 sbl = ( (hwght*hr)*s_b(iin,jin+1) + (hwght*hl + hr*hl)*s_b(iin,jin) ) * idenom
1407 sbr = ( (hwght*hl)*s_b(iin,jin) + (hwght*hr + hr*hl)*s_b(iin,jin+1) ) * idenom
1409 ttl = t_t(iin,jin); tbl = t_b(iin,jin); ttr = t_t(iin,jin+1); tbr = t_b(iin,jin+1)
1410 stl = s_t(iin,jin); sbl = s_b(iin,jin); str = s_t(iin,jin+1); sbr = s_b(iin,jin+1)
1414 w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
1415 dz_y(m,i) = w_left*(z_t(iin,jin) - z_b(iin,jin)) + w_right*(z_t(iin,jin+1) - z_b(iin,jin+1))
1422 t15(pos+1) = w_left*ttl + w_right*ttr
1423 t15(pos+5) = w_left*tbl + w_right*tbr
1425 s15(pos+1) = w_left*stl + w_right*str
1426 s15(pos+5) = w_left*sbl + w_right*sbr
1428 p15(pos+1) = -gxrho*(w_left*z_t(iin,jin) + w_right*z_t(iin,jin+1))
1431 do n=2,5 ; p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_y(m,i) ;
enddo
1435 weight_t = 0.25 * real(5-n)
1436 weight_b = 1.0 - weight_t
1437 s15(pos+n) = weight_t * s15(pos+1) + weight_b * s15(pos+5)
1438 t15(pos+n) = weight_t * t15(pos+1) + weight_b * t15(pos+5)
1444 r15(15*hio%isc+1:), 1, 15*(hio%iec-hio%isc+1), eos, rho_ref)
1445 do i=hio%isc,hio%iec ; iin = i+ioff
1446 intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
1451 intz(m) = g_e*dz_y(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + &
1452 32.0*(r15(pos+2)+r15(pos+4)) + &
1456 inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
1469 rho_ref, G_e, EOS, P_b, z_out, z_tol)
1470 real,
intent(in) :: t_t
1471 real,
intent(in) :: t_b
1472 real,
intent(in) :: s_t
1473 real,
intent(in) :: s_b
1474 real,
intent(in) :: z_t
1475 real,
intent(in) :: z_b
1476 real,
intent(in) :: p_t
1477 real,
intent(in) :: p_tgt
1478 real,
intent(in) :: rho_ref
1479 real,
intent(in) :: g_e
1481 real,
intent(out) :: p_b
1482 real,
intent(out) :: z_out
1483 real,
optional,
intent(in) :: z_tol
1485 real :: top_weight, bottom_weight, rho_anom, w_left, w_right, gxrho, dz, dp, f_guess, f_l, f_r
1486 real :: pa, pa_left, pa_right, pa_tol
1488 gxrho = g_e * rho_ref
1491 dp =
frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, 1.0, eos)
1495 if (p_tgt <= p_t )
then
1500 if (p_tgt >= p_b)
then
1506 pa_left = p_t - p_tgt
1508 pa_right = p_b - p_tgt
1509 pa_tol = gxrho * 1.e-5
1510 if (
present(z_tol)) pa_tol = gxrho * z_tol
1511 f_guess = f_l - pa_left / ( pa_right -pa_left ) * ( f_r - f_l )
1512 pa = pa_right - pa_left
1513 do while ( abs(pa) > pa_tol )
1515 z_out = z_t + ( z_b - z_t ) * f_guess
1516 pa =
frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, f_guess, eos) - ( p_tgt - p_t )
1518 if (pa<pa_left)
then
1519 write(0,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
1520 stop
'Blurgh! Too negative'
1524 elseif (pa>pa_right)
then
1525 write(0,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
1526 stop
'Blurgh! Too positive'
1533 f_guess = f_l - pa_left / ( pa_right -pa_left ) * ( f_r - f_l )
1541 real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS)
1542 real,
intent(in) :: t_t
1543 real,
intent(in) :: t_b
1544 real,
intent(in) :: s_t
1545 real,
intent(in) :: s_b
1546 real,
intent(in) :: z_t
1547 real,
intent(in) :: z_b
1548 real,
intent(in) :: rho_ref
1550 real,
intent(in) :: g_e
1551 real,
intent(in) :: pos
1554 real,
parameter :: c1_90 = 1.0/90.0
1555 real :: dz, top_weight, bottom_weight, rho_ave
1556 real,
dimension(5) :: t5, s5, p5, rho5
1561 bottom_weight = 0.25*real(n-1) * pos
1562 top_weight = 1.0 - bottom_weight
1564 s5(n) = top_weight * s_t + bottom_weight * s_b
1565 t5(n) = top_weight * t_t + bottom_weight * t_b
1566 p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( g_e * rho_ref )
1572 rho_ave = c1_90*(7.0*(rho5(1)+rho5(5)) + 32.0*(rho5(2)+rho5(4)) + 12.0*rho5(3))
1574 dz = ( z_t - z_b ) * pos
1583 z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, &
1584 EOS, dpa, intz_dpa, intx_dpa, inty_dpa)
1588 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1590 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1592 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1594 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1596 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1598 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1600 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1602 real,
dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1604 real,
intent(in) :: rho_ref
1606 real,
intent(in) :: rho_0
1608 real,
intent(in) :: g_e
1610 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1612 real,
dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1613 optional,
intent(out) :: intz_dpa
1616 real,
dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
1617 optional,
intent(out) :: intx_dpa
1620 real,
dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
1621 optional,
intent(out) :: inty_dpa
1637 real :: t5(5), s5(5), p5(5), r5(5)
1639 real :: w_left, w_right, intz(5)
1640 real,
parameter :: c1_90 = 1.0/90.0
1641 real :: gxrho, i_rho
1643 real :: weight_t, weight_b
1647 real :: t_top, t_mid, t_bot
1648 real :: s_top, s_mid, s_bot
1649 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m, n, ioff, joff
1650 real,
dimension(4) :: x, y
1651 real,
dimension(9) :: s_node, t_node, p_node, r_node
1655 "int_density_dz_generic_ppm: the implementation is not done yet, contact developer")
1657 ioff = hio%idg_offset - hii%idg_offset
1658 joff = hio%jdg_offset - hii%jdg_offset
1662 isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
1663 jsq = hio%JscB + joff ; jeq = hio%JecB + joff
1664 is = hio%isc + ioff ; ie = hio%iec + ioff
1665 js = hio%jsc + joff ; je = hio%jec + joff
1673 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1674 dz = z_t(i,j) - z_b(i,j)
1678 s1 = 6.0 * s(i,j) - 4.0 * s_t(i,j) - 2.0 * s_b(i,j)
1679 s2 = 3.0 * ( s_t(i,j) + s_b(i,j) - 2.0*s(i,j) )
1683 t1 = 6.0 * t(i,j) - 4.0 * t_t(i,j) - 2.0 * t_b(i,j)
1684 t2 = 3.0 * ( t_t(i,j) + t_b(i,j) - 2.0*t(i,j) )
1687 p5(n) = -gxrho*(z_t(i,j) - 0.25*real(n-1)*dz)
1690 xi = 0.25 * ( n - 1 )
1691 s5(n) = s0 + s1 * xi + s2 * xi**2
1692 t5(n) = t0 + t1 * xi + t2 * xi**2
1701 rho_anom = 1000.0 + s(i,j) - rho_ref
1702 dpa(i-ioff,j-joff) = g_e*dz*rho_anom
1710 intz_dpa(i-ioff,j-joff) = 0.5 * g_e * dz**2 * ( 1000.0 - rho_ref + s0 + s1/3.0 + &
1717 if (
present(intx_dpa))
then ;
do j=js,je ;
do i=isq,ieq
1718 intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
1720 w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
1721 dz = w_left*(z_t(i,j) - z_b(i,j)) + w_right*(z_t(i+1,j) - z_b(i+1,j))
1727 t_top = w_left*t_t(i,j) + w_right*t_t(i+1,j)
1728 t_mid = w_left*t(i,j) + w_right*t(i+1,j)
1729 t_bot = w_left*t_b(i,j) + w_right*t_b(i+1,j)
1731 s_top = w_left*s_t(i,j) + w_right*s_t(i+1,j)
1732 s_mid = w_left*s(i,j) + w_right*s(i+1,j)
1733 s_bot = w_left*s_b(i,j) + w_right*s_b(i+1,j)
1735 p5(1) = -gxrho*(w_left*z_t(i,j) + w_right*z_t(i+1,j))
1739 p5(n) = p5(n-1) + gxrho*0.25*dz
1744 s1 = 6.0 * s_mid - 4.0 * s_top - 2.0 * s_bot
1745 s2 = 3.0 * ( s_top + s_bot - 2.0*s_mid )
1749 t1 = 6.0 * t_mid - 4.0 * t_top - 2.0 * t_bot
1750 t2 = 3.0 * ( t_top + t_bot - 2.0*t_mid )
1754 xi = 0.25 * ( n - 1 )
1755 s5(n) = s0 + s1 * xi + s2 * xi**2
1756 t5(n) = t0 + t1 * xi + t2 * xi**2
1762 intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + &
1763 12.0*r5(3)) - rho_ref)
1765 intx_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
1788 s1 = 6.0 * s(i,j) - 4.0 * s_t(i,j) - 2.0 * s_b(i,j)
1789 s2 = 3.0 * ( s_t(i,j) + s_b(i,j) - 2.0 * s(i,j) )
1791 s_node(6) = s0 + 0.5 * s1 + 0.25 * s2
1792 s_node(3) = s0 + s1 + s2
1796 s1 = 6.0 * s(i+1,j) - 4.0 * s_t(i+1,j) - 2.0 * s_b(i+1,j)
1797 s2 = 3.0 * ( s_t(i+1,j) + s_b(i+1,j) - 2.0 * s(i+1,j) )
1799 s_node(8) = s0 + 0.5 * s1 + 0.25 * s2
1800 s_node(4) = s0 + s1 + s2
1802 s_node(5) = 0.5 * ( s_node(2) + s_node(1) )
1803 s_node(9) = 0.5 * ( s_node(6) + s_node(8) )
1804 s_node(7) = 0.5 * ( s_node(3) + s_node(4) )
1807 r_node = r_node - rho_ref
1811 intx_dpa(i-ioff,j-joff) = intx_dpa(i-ioff,j-joff) * g_e
1813 enddo ;
enddo ;
endif
1818 if (
present(inty_dpa))
then
1819 call mom_error(warning,
"int_density_dz_generic_ppm still needs to be written for inty_dpa!")
1820 do j=jsq,jeq ;
do i=is,ie
1822 inty_dpa(i-ioff,j-joff) = 0.0
1834 real,
dimension(4),
intent(in) :: x
1835 real,
dimension(4),
intent(in) :: y
1836 real,
dimension(9),
intent(in) :: f
1837 real,
intent(out) :: integral
1841 real,
dimension(9) :: weight, xi, eta
1843 real :: dxdxi, dxdeta
1844 real :: dydxi, dydeta
1845 real,
dimension(4) :: phiiso, dphiisodxi, dphiisodeta
1846 real,
dimension(9) :: phi, dphidxi, dphideta
1863 weight(1) = 25.0/81.0 ; xi(1) = -t ; eta(1) = t
1864 weight(2) = 40.0/81.0 ; xi(2) = .0 ; eta(2) = t
1865 weight(3) = 25.0/81.0 ; xi(3) = t ; eta(3) = t
1866 weight(4) = 40.0/81.0 ; xi(4) = -t ; eta(4) = .0
1867 weight(5) = 64.0/81.0 ; xi(5) = .0 ; eta(5) = .0
1868 weight(6) = 40.0/81.0 ; xi(6) = t ; eta(6) = .0
1869 weight(7) = 25.0/81.0 ; xi(7) = -t ; eta(7) = -t
1870 weight(8) = 40.0/81.0 ; xi(8) = .0 ; eta(8) = -t
1871 weight(9) = 25.0/81.0 ; xi(9) = t ; eta(9) = -t
1880 dphiisodxi, dphiisodeta )
1889 dxdxi = dxdxi + x(i) * dphiisodxi(i)
1890 dxdeta = dxdeta + x(i) * dphiisodeta(i)
1891 dydxi = dydxi + y(i) * dphiisodxi(i)
1892 dydeta = dydeta + y(i) * dphiisodeta(i)
1896 jacobian_k = dxdxi*dydeta - dydxi*dxdeta
1904 f_k = f_k + f(i) * phi(i)
1907 integral = integral + weight(k) * f_k * jacobian_k
1917 real,
intent(in) :: xi
1918 real,
intent(in) :: eta
1919 real,
dimension(4),
intent(out) :: phi
1920 real,
dimension(4),
intent(out) :: dphidxi
1922 real,
dimension(4),
intent(out) :: dphideta
1935 phi(1) = 0.25 * ( 1 + xi ) * ( 1 + eta )
1936 phi(2) = 0.25 * ( 1 - xi ) * ( 1 + eta )
1937 phi(3) = 0.25 * ( 1 - xi ) * ( 1 - eta )
1938 phi(4) = 0.25 * ( 1 + xi ) * ( 1 - eta )
1940 dphidxi(1) = 0.25 * ( 1 + eta )
1941 dphidxi(2) = - 0.25 * ( 1 + eta )
1942 dphidxi(3) = - 0.25 * ( 1 - eta )
1943 dphidxi(4) = 0.25 * ( 1 - eta )
1945 dphideta(1) = 0.25 * ( 1 + xi )
1946 dphideta(2) = 0.25 * ( 1 - xi )
1947 dphideta(3) = - 0.25 * ( 1 - xi )
1948 dphideta(4) = - 0.25 * ( 1 + xi )
1958 real,
intent(in) :: xi
1959 real,
intent(in) :: eta
1960 real,
dimension(9),
intent(out) :: phi
1962 real,
dimension(9),
intent(out) :: dphidxi
1964 real,
dimension(9),
intent(out) :: dphideta
1984 phi(1) = 0.25 * xi * ( 1 + xi ) * eta * ( 1 + eta )
1985 phi(2) = - 0.25 * xi * ( 1 - xi ) * eta * ( 1 + eta )
1986 phi(3) = 0.25 * xi * ( 1 - xi ) * eta * ( 1 - eta )
1987 phi(4) = - 0.25 * xi * ( 1 + xi ) * eta * ( 1 - eta )
1988 phi(5) = 0.5 * ( 1 + xi ) * ( 1 - xi ) * eta * ( 1 + eta )
1989 phi(6) = - 0.5 * xi * ( 1 - xi ) * ( 1 - eta ) * ( 1 + eta )
1990 phi(7) = - 0.5 * ( 1 - xi ) * ( 1 + xi ) * eta * ( 1 - eta )
1991 phi(8) = 0.5 * xi * ( 1 + xi ) * ( 1 - eta ) * ( 1 + eta )
1992 phi(9) = ( 1 - xi ) * ( 1 + xi ) * ( 1 - eta ) * ( 1 + eta )
2022 dza, intp_dza, intx_dza, inty_dza, halo_size, &
2023 bathyP, dP_neglect, useMassWghtInterp)
2025 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2027 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2029 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2031 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2033 real,
intent(in) :: alpha_ref
2039 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2042 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2043 optional,
intent(out) :: intp_dza
2046 real,
dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
2047 optional,
intent(out) :: intx_dza
2050 real,
dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
2051 optional,
intent(out) :: inty_dza
2054 integer,
optional,
intent(in) :: halo_size
2055 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2056 optional,
intent(in) :: bathyp
2057 real,
optional,
intent(in) :: dp_neglect
2059 logical,
optional,
intent(in) :: usemasswghtinterp
2069 real :: t5(5), s5(5), p5(5), a5(5)
2076 real :: hwt_ll, hwt_lr
2077 real :: hwt_rl, hwt_rr
2079 real :: wtt_l, wtt_r
2082 logical :: do_massweight
2083 real,
parameter :: c1_90 = 1.0/90.0
2084 integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, n, halo
2086 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
2087 halo = 0 ;
if (
present(halo_size)) halo = max(halo_size,0)
2088 ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
2089 if (
present(intx_dza))
then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh);
endif
2090 if (
present(inty_dza))
then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh);
endif
2092 do_massweight = .false.
2093 if (
present(usemasswghtinterp))
then ;
if (usemasswghtinterp)
then
2094 do_massweight = .true.
2095 if (.not.
present(bathyp))
call mom_error(fatal,
"int_spec_vol_dp_generic: "//&
2096 "bathyP must be present if useMassWghtInterp is present and true.")
2097 if (.not.
present(dp_neglect))
call mom_error(fatal,
"int_spec_vol_dp_generic: "//&
2098 "dP_neglect must be present if useMassWghtInterp is present and true.")
2101 do j=jsh,jeh ;
do i=ish,ieh
2102 dp = p_b(i,j) - p_t(i,j)
2104 t5(n) = t(i,j) ; s5(n) = s(i,j)
2105 p5(n) = p_b(i,j) - 0.25*real(n-1)*dp
2110 alpha_anom = c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3))
2111 dza(i,j) = dp*alpha_anom
2114 if (
present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &
2115 (alpha_anom - c1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) )
2118 if (
present(intx_dza))
then ;
do j=hi%jsc,hi%jec ;
do i=isq,ieq
2123 if (do_massweight) &
2124 hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
2125 if (hwght > 0.)
then
2126 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
2127 hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
2128 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
2129 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
2130 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
2131 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
2133 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
2136 intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
2138 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
2139 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
2143 p5(1) = wt_l*p_b(i,j) + wt_r*p_b(i+1,j)
2144 dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i+1,j) - p_t(i+1,j))
2145 t5(1) = wtt_l*t(i,j) + wtt_r*t(i+1,j)
2146 s5(1) = wtt_l*s(i,j) + wtt_r*s(i+1,j)
2149 t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) - 0.25*dp
2154 intp(m) = dp*( c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + &
2158 intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
2160 enddo ;
enddo ;
endif
2162 if (
present(inty_dza))
then ;
do j=jsq,jeq ;
do i=hi%isc,hi%iec
2167 if (do_massweight) &
2168 hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
2169 if (hwght > 0.)
then
2170 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
2171 hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
2172 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
2173 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
2174 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
2175 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
2177 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
2180 intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
2182 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
2183 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
2187 p5(1) = wt_l*p_b(i,j) + wt_r*p_b(i,j+1)
2188 dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i,j+1) - p_t(i,j+1))
2189 t5(1) = wtt_l*t(i,j) + wtt_r*t(i,j+1)
2190 s5(1) = wtt_l*s(i,j) + wtt_r*s(i,j+1)
2192 t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) - 0.25*dp
2197 intp(m) = dp*( c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + &
2201 inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
2203 enddo ;
enddo ;
endif
2212 dP_neglect, bathyP, HI, EOS, dza, &
2213 intp_dza, intx_dza, inty_dza, useMassWghtInterp)
2215 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2217 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2219 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2221 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2223 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2225 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2227 real,
intent(in) :: alpha_ref
2232 real,
intent(in) :: dp_neglect
2234 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2235 intent(in) :: bathyp
2237 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2240 real,
dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2241 optional,
intent(out) :: intp_dza
2244 real,
dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
2245 optional,
intent(out) :: intx_dza
2248 real,
dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
2249 optional,
intent(out) :: inty_dza
2252 logical,
optional,
intent(in) :: usemasswghtinterp
2262 real,
dimension(5) :: t5, s5, p5, a5
2263 real,
dimension(15) :: t15, s15, p15, a15
2264 real :: wt_t(5), wt_b(5)
2265 real :: t_top, t_bot, s_top, s_bot, p_top, p_bot
2273 real :: hwt_ll, hwt_lr
2274 real :: hwt_rl, hwt_rr
2276 real :: wtt_l, wtt_r
2279 real,
parameter :: c1_90 = 1.0/90.0
2280 logical :: do_massweight
2281 integer :: isq, ieq, jsq, jeq, i, j, m, n, pos
2283 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
2285 do_massweight = .false.
2286 if (
present(usemasswghtinterp)) do_massweight = usemasswghtinterp
2289 wt_t(n) = 0.25 * real(n-1)
2290 wt_b(n) = 1.0 - wt_t(n)
2296 do j=jsq,jeq+1;
do i=isq,ieq+1
2297 dp = p_b(i,j) - p_t(i,j)
2299 p5(n) = wt_t(n) * p_t(i,j) + wt_b(n) * p_b(i,j)
2300 s5(n) = wt_t(n) * s_t(i,j) + wt_b(n) * s_b(i,j)
2301 t5(n) = wt_t(n) * t_t(i,j) + wt_b(n) * t_b(i,j)
2306 alpha_anom = c1_90*((7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4))) + 12.0*a5(3))
2307 dza(i,j) = dp*alpha_anom
2310 if (
present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &
2311 (alpha_anom - c1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) )
2317 if (
present(intx_dza))
then ;
do j=hi%jsc,hi%jec ;
do i=isq,ieq
2324 if (do_massweight) &
2325 hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
2326 if (hwght > 0.)
then
2327 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
2328 hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
2329 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
2330 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
2331 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
2332 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
2334 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
2338 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
2339 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
2343 p_top = wt_l*p_t(i,j) + wt_r*p_t(i+1,j)
2344 p_bot = wt_l*p_b(i,j) + wt_r*p_b(i+1,j)
2345 t_top = wtt_l*t_t(i,j) + wtt_r*t_t(i+1,j)
2346 t_bot = wtt_l*t_b(i,j) + wtt_r*t_b(i+1,j)
2347 s_top = wtt_l*s_t(i,j) + wtt_r*s_t(i+1,j)
2348 s_bot = wtt_l*s_b(i,j) + wtt_r*s_b(i+1,j)
2349 dp_90(m) = c1_90*(p_bot - p_top)
2354 p15(pos+n) = wt_t(n) * p_top + wt_b(n) * p_bot
2355 s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot
2356 t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot
2362 intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
2367 intp(m) = dp_90(m)*((7.0*(a15(pos+1)+a15(pos+5)) + &
2368 32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3))
2371 intx_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
2373 enddo ;
enddo ;
endif
2378 if (
present(inty_dza))
then ;
do j=jsq,jeq ;
do i=hi%isc,hi%iec
2383 if (do_massweight) &
2384 hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
2385 if (hwght > 0.)
then
2386 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
2387 hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
2388 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
2389 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
2390 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
2391 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
2393 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
2397 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
2398 wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
2402 p_top = wt_l*p_t(i,j) + wt_r*p_t(i,j+1)
2403 p_bot = wt_l*p_b(i,j) + wt_r*p_b(i,j+1)
2404 t_top = wtt_l*t_t(i,j) + wtt_r*t_t(i,j+1)
2405 t_bot = wtt_l*t_b(i,j) + wtt_r*t_b(i,j+1)
2406 s_top = wtt_l*s_t(i,j) + wtt_r*s_t(i,j+1)
2407 s_bot = wtt_l*s_b(i,j) + wtt_r*s_b(i,j+1)
2408 dp_90(m) = c1_90*(p_bot - p_top)
2413 p15(pos+n) = wt_t(n) * p_top + wt_b(n) * p_bot
2414 s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot
2415 t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot
2421 intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
2426 intp(m) = dp_90(m) * ((7.0*(a15(pos+1)+a15(pos+5)) + &
2427 32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3))
2430 inty_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
2432 enddo ;
enddo ;
endif
2441 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)), &
2443 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)), &
2445 real,
dimension(:),
intent(in) :: press
2447 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)), &
2448 intent(in) :: mask_z
2449 integer,
intent(in) :: kd
2452 real :: gsw_sr_from_sp, gsw_ct_from_pt, gsw_sa_from_sp
2455 if (.not.
associated(eos))
call mom_error(fatal, &
2456 "convert_temp_salt_to_TEOS10 called with an unassociated EOS_type EOS.")
2460 do k=1,kd ;
do j=g%jsc,g%jec ;
do i=g%isc,g%iec
2461 if (mask_z(i,j,k) >= 1.0)
then
2462 s(i,j,k) = gsw_sr_from_sp(s(i,j,k))
2465 t(i,j,k) = gsw_ct_from_pt(s(i,j,k),t(i,j,k))
2467 enddo ;
enddo ;
enddo
2471 subroutine extract_member_eos(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, &
2472 Rho_T0_S0, drho_dT, dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp)
2474 integer,
optional,
intent(out) :: form_of_eos
2475 integer,
optional,
intent(out) :: form_of_tfreeze
2477 logical,
optional,
intent(out) :: eos_quadrature
2479 logical,
optional,
intent(out) :: compressible
2480 real ,
optional,
intent(out) :: rho_t0_s0
2481 real ,
optional,
intent(out) :: drho_dt
2483 real ,
optional,
intent(out) :: drho_ds
2485 real ,
optional,
intent(out) :: tfr_s0_p0
2486 real ,
optional,
intent(out) :: dtfr_ds
2488 real ,
optional,
intent(out) :: dtfr_dp
2491 if (
present(form_of_eos )) form_of_eos = eos%form_of_EOS
2492 if (
present(form_of_tfreeze)) form_of_tfreeze = eos%form_of_TFreeze
2493 if (
present(eos_quadrature )) eos_quadrature = eos%EOS_quadrature
2494 if (
present(compressible )) compressible = eos%Compressible
2495 if (
present(rho_t0_s0 )) rho_t0_s0 = eos%Rho_T0_S0
2496 if (
present(drho_dt )) drho_dt = eos%drho_dT
2497 if (
present(drho_ds )) drho_ds = eos%dRho_dS
2498 if (
present(tfr_s0_p0 )) tfr_s0_p0 = eos%TFr_S0_P0
2499 if (
present(dtfr_ds )) dtfr_ds = eos%dTFr_dS
2500 if (
present(dtfr_dp )) dtfr_dp = eos%dTFr_dp