!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| module state_mod !BOP ! !MODULE: state_mod ! ! !DESCRIPTION: ! This module contains routines necessary for computing the density ! from model temperature and salinity using an equation of state. ! ! The module supports four forms of EOS: ! \begin{enumerate} ! \item The UNESCO equation of state computed using the ! potential-temperature-based bulk modulus from Jackett and ! McDougall, JTECH, Vol.12, pp 381-389, April, 1995. ! \item The faster and more accurate alternative to the UNESCO eos ! of McDougall, Wright, Jackett and Feistel (hereafter ! MWJF, 2001 submission to JTECH). ! \item a polynomial fit to the full UNESCO EOS ! \item a simple linear EOS based on constant expansion coeffs ! \end{enumerate} ! ! !REVISION HISTORY: ! SVN:$Id: state_mod.F90 26114 2010-12-17 20:29:34Z njn01 $ ! !USES: use kinds_mod use blocks use distribution use domain use constants use grid use io use broadcast use time_management use exit_mod use shr_sys_mod #ifdef CCSMCOUPLED !*** ccsm use shr_vmath_mod #endif implicit none private save ! !PUBLIC MEMBER FUNCTIONS: public :: state, & init_state, & ref_pressure, & state_singlept !EOP !BOC !----------------------------------------------------------------------- ! ! valid ranges and pressure as function of depth ! !----------------------------------------------------------------------- real (r8), dimension(km) :: & tmin, tmax, &! valid temperature range for level k smin, smax, &! valid salinity range for level k pressz ! ref pressure (bars) at each level !----------------------------------------------------------------------- ! ! choices for eos type and valid range checks ! !----------------------------------------------------------------------- integer (int_kind), public, parameter :: & state_type_jmcd = 1, &! integer ids for state choice state_type_mwjf = 2, & state_type_polynomial = 3, & state_type_linear = 4 integer (int_kind), public :: & state_itype ! input state type chosen integer (int_kind), public :: & state_range_iopt, &! option for checking valid T,S range state_range_freq ! freq (in steps) for checking T,S range integer (int_kind), public, parameter :: & state_range_ignore = 1, &! do not check T,S range state_range_check = 2, &! check T,S range and report invalid state_range_enforce = 3 ! force polynomial eval within range !----------------------------------------------------------------------- ! ! UNESCO EOS constants and JMcD bulk modulus constants ! !----------------------------------------------------------------------- !*** for density of fresh water (standard UNESCO) real (r8), parameter :: & unt0 = 999.842594_r8, & unt1 = 6.793952e-2_r8, & unt2 = -9.095290e-3_r8, & unt3 = 1.001685e-4_r8, & unt4 = -1.120083e-6_r8, & unt5 = 6.536332e-9_r8 !*** for dependence of surface density on salinity (UNESCO) real (r8), parameter :: & uns1t0 = 0.824493_r8 , & uns1t1 = -4.0899e-3_r8, & uns1t2 = 7.6438e-5_r8, & uns1t3 = -8.2467e-7_r8, & uns1t4 = 5.3875e-9_r8, & unsqt0 = -5.72466e-3_r8, & unsqt1 = 1.0227e-4_r8, & unsqt2 = -1.6546e-6_r8, & uns2t0 = 4.8314e-4_r8 !*** from Table A1 of Jackett and McDougall real (r8), parameter :: & bup0s0t0 = 1.965933e+4_r8, & bup0s0t1 = 1.444304e+2_r8, & bup0s0t2 = -1.706103_r8 , & bup0s0t3 = 9.648704e-3_r8, & bup0s0t4 = -4.190253e-5_r8 real (r8), parameter :: & bup0s1t0 = 5.284855e+1_r8, & bup0s1t1 = -3.101089e-1_r8, & bup0s1t2 = 6.283263e-3_r8, & bup0s1t3 = -5.084188e-5_r8 real (r8), parameter :: & bup0sqt0 = 3.886640e-1_r8, & bup0sqt1 = 9.085835e-3_r8, & bup0sqt2 = -4.619924e-4_r8 real (r8), parameter :: & bup1s0t0 = 3.186519_r8 , & bup1s0t1 = 2.212276e-2_r8, & bup1s0t2 = -2.984642e-4_r8, & bup1s0t3 = 1.956415e-6_r8 real (r8), parameter :: & bup1s1t0 = 6.704388e-3_r8, & bup1s1t1 = -1.847318e-4_r8, & bup1s1t2 = 2.059331e-7_r8, & bup1sqt0 = 1.480266e-4_r8 real (r8), parameter :: & bup2s0t0 = 2.102898e-4_r8, & bup2s0t1 = -1.202016e-5_r8, & bup2s0t2 = 1.394680e-7_r8, & bup2s1t0 = -2.040237e-6_r8, & bup2s1t1 = 6.128773e-8_r8, & bup2s1t2 = 6.207323e-10_r8 !----------------------------------------------------------------------- ! ! MWJF EOS coefficients ! !----------------------------------------------------------------------- !*** these constants will be used to construct the numerator !*** factor unit change (kg/m^3 -> g/cm^3) into numerator terms real (r8), parameter :: & mwjfnp0s0t0 = 9.99843699e+2_r8 * p001, & mwjfnp0s0t1 = 7.35212840e+0_r8 * p001, & mwjfnp0s0t2 = -5.45928211e-2_r8 * p001, & mwjfnp0s0t3 = 3.98476704e-4_r8 * p001, & mwjfnp0s1t0 = 2.96938239e+0_r8 * p001, & mwjfnp0s1t1 = -7.23268813e-3_r8 * p001, & mwjfnp0s2t0 = 2.12382341e-3_r8 * p001, & mwjfnp1s0t0 = 1.04004591e-2_r8 * p001, & mwjfnp1s0t2 = 1.03970529e-7_r8 * p001, & mwjfnp1s1t0 = 5.18761880e-6_r8 * p001, & mwjfnp2s0t0 = -3.24041825e-8_r8 * p001, & mwjfnp2s0t2 = -1.23869360e-11_r8* p001 !*** these constants will be used to construct the denominator real (kind=r8), parameter :: & mwjfdp0s0t0 = 1.0e+0_r8, & mwjfdp0s0t1 = 7.28606739e-3_r8, & mwjfdp0s0t2 = -4.60835542e-5_r8, & mwjfdp0s0t3 = 3.68390573e-7_r8, & mwjfdp0s0t4 = 1.80809186e-10_r8, & mwjfdp0s1t0 = 2.14691708e-3_r8, & mwjfdp0s1t1 = -9.27062484e-6_r8, & mwjfdp0s1t3 = -1.78343643e-10_r8, & mwjfdp0sqt0 = 4.76534122e-6_r8, & mwjfdp0sqt2 = 1.63410736e-9_r8, & mwjfdp1s0t0 = 5.30848875e-6_r8, & mwjfdp2s0t3 = -3.03175128e-16_r8, & mwjfdp3s0t1 = -1.27934137e-17_r8 !----------------------------------------------------------------------- ! ! coeffs and reference values for polynomial eos ! !----------------------------------------------------------------------- real (r8), dimension(:), allocatable :: & to, &! reference temperature for level k so, &! reference salinity for level k sigo ! reference density for level k real (r8), dimension(:,:), allocatable :: & state_coeffs ! coefficients for polynomial eos !----------------------------------------------------------------------- ! ! parameters for linear eos ! !----------------------------------------------------------------------- real (r8), parameter :: & T_leos_ref = 19.0_r8, &! reference T for linear eos (deg C) S_leos_ref = 0.035_r8, &! reference S for linear eos (msu) rho_leos_ref = 1.025022_r8, &! ref dens (g/cm3) at ref T,S and 0 bar alf = 2.55e-4_r8, &! expansion coeff -(drho/dT) (gr/cm^3/K) bet = 7.64e-1_r8 ! expansion coeff (drho/dS) (gr/cm^3/msu) !----------------------------------------------------------------------- ! ! logical control parameters ! !----------------------------------------------------------------------- logical (log_kind) :: & lccsm_control_compatible !EOC !*********************************************************************** contains !*********************************************************************** !BOP ! !IROUTINE: state ! !INTERFACE: subroutine state(k, kk, TEMPK, SALTK, this_block, & RHOOUT, RHOFULL, DRHODT, DRHODS) ! !DESCRIPTION: ! Returns the density of water at level k from equation of state ! $\rho = \rho(d,\theta,S)$ where $d$ is depth, $\theta$ is ! potential temperature, and $S$ is salinity. the density can be ! returned as a perturbation (RHOOUT) or as the full density ! (RHOFULL). Note that only the polynomial EOS choice will return ! a perturbation density; in other cases the full density is returned ! regardless of which argument is requested. ! ! This routine also computes derivatives of density with respect ! to temperature and salinity at level k from equation of state ! if requested (ie the optional arguments are present). ! ! If $k = kk$ are equal the density for level k is returned. ! If $k \neq kk$ the density returned is that for a parcel ! adiabatically displaced from level k to level kk. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: integer (int_kind), intent(in) :: & k, &! depth level index kk ! level to which water is adiabatically ! displaced real (r8), dimension(nx_block,ny_block), intent(in) :: & TEMPK, &! temperature at level k SALTK ! salinity at level k type (block), intent(in) :: & this_block ! block information for current block ! !OUTPUT PARAMETERS: real (r8), dimension(nx_block,ny_block), optional, intent(out) :: & RHOOUT, &! perturbation density of water RHOFULL, &! full density of water DRHODT, &! derivative of density with respect to temperature DRHODS ! derivative of density with respect to salinity !EOP !BOC !----------------------------------------------------------------------- ! ! local variables: ! !----------------------------------------------------------------------- integer (int_kind) :: & ib,ie,jb,je, &! extent of physical domain bid, &! local block index out_of_range ! counter for out-of-range T,S values real (r8), dimension(nx_block,ny_block) :: & TQ,SQ, &! adjusted T,S BULK_MOD, &! Bulk modulus RHO_S, &! density at the surface DRDT0, &! d(density)/d(temperature), for surface DRDS0, &! d(density)/d(salinity ), for surface DKDT, &! d(bulk modulus)/d(pot. temp.) DKDS, &! d(bulk modulus)/d(salinity ) SQR,DENOMK, &! work arrays WORK1, WORK2, WORK3, WORK4, T2 real (r8) :: p, p2 ! temporary pressure scalars !*** MWJF numerator coefficients including pressure real (r8) :: & mwjfnums0t0, mwjfnums0t1, mwjfnums0t2, mwjfnums0t3, & mwjfnums1t0, mwjfnums1t1, mwjfnums2t0, & mwjfdens0t0, mwjfdens0t1, mwjfdens0t2, mwjfdens0t3, mwjfdens0t4, & mwjfdens1t0, mwjfdens1t1, mwjfdens1t3, & mwjfdensqt0, mwjfdensqt2 !----------------------------------------------------------------------- ! ! first check for valid range if requested ! !----------------------------------------------------------------------- bid = this_block%local_id select case (state_range_iopt) case (state_range_ignore) !*** prevent problems with garbage on land points or ghost cells TQ = min(TEMPK, c1000) TQ = max(TQ, -c1000) SQ = min(SALTK, c1000) SQ = max(SALTK, c0) case (state_range_check) call exit_POP(sigAbort, & '(state) ERROR unsupported option -- must define time flag and use check_time_flag') !*** if (time_to_do(freq_opt_nstep, state_range_freq)) then ib = this_block%ib ie = this_block%ie jb = this_block%jb je = this_block%je out_of_range = count((TEMPK(ib:ie,jb:je) < tmin(kk) .or. & TEMPK(ib:ie,jb:je) > tmax(kk)) .and. & k <= KMT(ib:ie,jb:je,bid)) if (out_of_range /= 0) & write(stdout,'(a9,i6,a44,i3)') 'WARNING: ',out_of_range, & 'points outside of valid temp range at level ',kk out_of_range = count((SALTK(ib:ie,jb:je) < smin(kk) .or. & SALTK(ib:ie,jb:je) > smax(kk)) .and. & k <= KMT(ib:ie,jb:je,bid)) if (out_of_range /= 0) & write(stdout,'(a9,i6,a44,i3)') 'WARNING: ',out_of_range, & 'points outside of valid salt range at level ',kk !*** endif TQ = TEMPK SQ = SALTK case (state_range_enforce) TQ = min(TEMPK,tmax(kk)) TQ = max(TQ,tmin(kk)) SQ = min(SALTK,smax(kk)) SQ = max(SQ,smin(kk)) end select !----------------------------------------------------------------------- ! ! now compute density or expansion coefficients ! !----------------------------------------------------------------------- select case (state_itype) !----------------------------------------------------------------------- ! ! McDougall, Wright, Jackett, and Feistel EOS ! test value : rho = 1.033213242 for ! S = 35.0 PSU, theta = 20.0, pressz = 200.0 ! !----------------------------------------------------------------------- case (state_type_mwjf) p = c10*pressz(kk) SQ = c1000*SQ #ifdef CCSMCOUPLED call shr_vmath_sqrt(SQ, SQR, nx_block*ny_block) #else SQR = sqrt(SQ) #endif !*** !*** first calculate numerator of MWJF density [P_1(S,T,p)] !*** mwjfnums0t0 = mwjfnp0s0t0 + p*(mwjfnp1s0t0 + p*mwjfnp2s0t0) mwjfnums0t1 = mwjfnp0s0t1 mwjfnums0t2 = mwjfnp0s0t2 + p*(mwjfnp1s0t2 + p*mwjfnp2s0t2) mwjfnums0t3 = mwjfnp0s0t3 mwjfnums1t0 = mwjfnp0s1t0 + p*mwjfnp1s1t0 mwjfnums1t1 = mwjfnp0s1t1 mwjfnums2t0 = mwjfnp0s2t0 WORK1 = mwjfnums0t0 + TQ * (mwjfnums0t1 + TQ * (mwjfnums0t2 + & mwjfnums0t3 * TQ)) + SQ * (mwjfnums1t0 + & mwjfnums1t1 * TQ + mwjfnums2t0 * SQ) !*** !*** now calculate denominator of MWJF density [P_2(S,T,p)] !*** mwjfdens0t0 = mwjfdp0s0t0 + p*mwjfdp1s0t0 mwjfdens0t1 = mwjfdp0s0t1 + p**3 * mwjfdp3s0t1 mwjfdens0t2 = mwjfdp0s0t2 mwjfdens0t3 = mwjfdp0s0t3 + p**2 * mwjfdp2s0t3 mwjfdens0t4 = mwjfdp0s0t4 mwjfdens1t0 = mwjfdp0s1t0 mwjfdens1t1 = mwjfdp0s1t1 mwjfdens1t3 = mwjfdp0s1t3 mwjfdensqt0 = mwjfdp0sqt0 mwjfdensqt2 = mwjfdp0sqt2 WORK2 = mwjfdens0t0 + TQ * (mwjfdens0t1 + TQ * (mwjfdens0t2 + & TQ * (mwjfdens0t3 + mwjfdens0t4 * TQ))) + & SQ * (mwjfdens1t0 + TQ * (mwjfdens1t1 + TQ*TQ*mwjfdens1t3)+ & SQR * (mwjfdensqt0 + TQ*TQ*mwjfdensqt2)) DENOMK = c1/WORK2 if (present(RHOOUT)) then RHOOUT = WORK1*DENOMK endif if (present(RHOFULL)) then RHOFULL = WORK1*DENOMK endif if (present(DRHODT)) then WORK3 = &! dP_1/dT mwjfnums0t1 + TQ * (c2*mwjfnums0t2 + & c3*mwjfnums0t3 * TQ) + mwjfnums1t1 * SQ WORK4 = &! dP_2/dT mwjfdens0t1 + SQ * mwjfdens1t1 + & TQ * (c2*(mwjfdens0t2 + SQ*SQR*mwjfdensqt2) + & TQ * (c3*(mwjfdens0t3 + SQ * mwjfdens1t3) + & TQ * c4*mwjfdens0t4)) DRHODT = (WORK3 - WORK1*DENOMK*WORK4)*DENOMK endif if (present(DRHODS)) then WORK3 = &! dP_1/dS mwjfnums1t0 + mwjfnums1t1 * TQ + c2*mwjfnums2t0 * SQ WORK4 = mwjfdens1t0 + &! dP_2/dS TQ * (mwjfdens1t1 + TQ*TQ*mwjfdens1t3) + & c1p5*SQR*(mwjfdensqt0 + TQ*TQ*mwjfdensqt2) DRHODS = (WORK3 - WORK1*DENOMK*WORK4)*DENOMK * c1000 endif !----------------------------------------------------------------------- ! ! Jackett and McDougall EOS ! !----------------------------------------------------------------------- case (state_type_jmcd) p = pressz(kk) p2 = p*p SQ = c1000*SQ SQR = sqrt(SQ) T2 = TQ*TQ !*** !*** first calculate surface (p=0) values from UNESCO eqns. !*** WORK1 = uns1t0 + uns1t1*TQ + & (uns1t2 + uns1t3*TQ + uns1t4*T2)*T2 WORK2 = SQR*(unsqt0 + unsqt1*TQ + unsqt2*T2) RHO_S = unt1*TQ + (unt2 + unt3*TQ + (unt4 + unt5*TQ)*T2)*T2 & + (uns2t0*SQ + WORK1 + WORK2)*SQ !*** !*** now calculate bulk modulus at pressure p from !*** Jackett and McDougall formula !*** WORK3 = bup0s1t0 + bup0s1t1*TQ + & (bup0s1t2 + bup0s1t3*TQ)*T2 + & p *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + & p2*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2) WORK4 = SQR*(bup0sqt0 + bup0sqt1*TQ + bup0sqt2*T2 + & bup1sqt0*p) BULK_MOD = bup0s0t0 + bup0s0t1*TQ + & (bup0s0t2 + bup0s0t3*TQ + bup0s0t4*T2)*T2 + & p *(bup1s0t0 + bup1s0t1*TQ + & (bup1s0t2 + bup1s0t3*TQ)*T2) + & p2*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + & SQ*(WORK3 + WORK4) DENOMK = c1/(BULK_MOD - p) !*** !*** now calculate required fields !*** if (present(RHOOUT)) then RHOOUT = merge(((unt0 + RHO_S)*BULK_MOD*DENOMK)*p001, & c0, KMT(:,:,bid) >= k) endif if (present(RHOFULL)) then RHOFULL = merge(((unt0 + RHO_S)*BULK_MOD*DENOMK)*p001, & c0, KMT(:,:,bid) >= k) endif if (present(DRHODT)) then DRDT0 = unt1 + c2*unt2*TQ + & (c3*unt3 + c4*unt4*TQ + c5*unt5*T2)*T2 + & (uns1t1 + c2*uns1t2*TQ + & (c3*uns1t3 + c4*uns1t4*TQ)*T2 + & (unsqt1 + c2*unsqt2*TQ)*SQR )*SQ DKDT = bup0s0t1 + c2*bup0s0t2*TQ + & (c3*bup0s0t3 + c4*bup0s0t4*TQ)*T2 + & p *(bup1s0t1 + c2*bup1s0t2*TQ + c3*bup1s0t3*T2) + & p2*(bup2s0t1 + c2*bup2s0t2*TQ) + & SQ*(bup0s1t1 + c2*bup0s1t2*TQ + c3*bup0s1t3*T2 + & p *(bup1s1t1 + c2*bup1s1t2*TQ) + & p2 *(bup2s1t1 + c2*bup2s1t2*TQ) + & SQR*(bup0sqt1 + c2*bup0sqt2*TQ)) DRHODT = merge ((DENOMK*(DRDT0*BULK_MOD - & pressz(kk)*(unt0+RHO_S)*DKDT*DENOMK))*p001, & c0, KMT(:,:,bid) >= k) endif if (present(DRHODS)) then DRDS0 = c2*uns2t0*SQ + WORK1 + c1p5*WORK2 DKDS = WORK3 + c1p5*WORK4 DRHODS = merge (DENOMK*(DRDS0*BULK_MOD - & pressz(kk)*(unt0+RHO_S)*DKDS*DENOMK), & c0, KMT(:,:,bid) >= k) endif !----------------------------------------------------------------------- ! ! polynomial EOS ! !----------------------------------------------------------------------- case (state_type_polynomial) TQ = TQ - to(kk) SQ = SQ - so(kk) - 0.035_r8 if (present(RHOOUT)) then !*** density as perturbation from a k-level reference RHOOUT = merge((state_coeffs(1,kk) + & (state_coeffs(4,kk) + & state_coeffs(7,kk)*SQ)*SQ + & (state_coeffs(3,kk) + & state_coeffs(8,kk)*SQ + & state_coeffs(6,kk)*TQ)*TQ)*TQ + & (state_coeffs(2,kk) + & (state_coeffs(5,kk) + & state_coeffs(9,kk)*SQ)*SQ)*SQ & ,c0, KMT(:,:,bid) >= k) endif if (present(RHOFULL)) then !*** full density (adding the reference density) RHOFULL=merge(((state_coeffs(1,kk) + & (state_coeffs(4,kk) + & state_coeffs(7,kk)*SQ)*SQ + & (state_coeffs(3,kk) + & state_coeffs(8,kk)*SQ + & state_coeffs(6,kk)*TQ)*TQ)*TQ + & (state_coeffs(2,kk) + & (state_coeffs(5,kk) + & state_coeffs(9,kk)*SQ)*SQ)*SQ) + & sigo(kk)*p001 + c1 & ,c1, KMT(:,:,bid) >= k) endif if (present(DRHODT)) then DRHODT = merge (state_coeffs(1,kk) + & (state_coeffs(4,kk) + & state_coeffs(7,kk)*SQ)*SQ + & (c2*state_coeffs(3,kk) + & c2*state_coeffs(8,kk)*SQ + & c3*state_coeffs(6,kk)*TQ)*TQ, & c0, KMT(:,:,bid) >= k) endif if (present(DRHODS)) then DRHODS = merge((state_coeffs(4,kk) + & c2*state_coeffs(7,kk)*SQ + & state_coeffs(8,kk)*TQ)*TQ + & state_coeffs(2,kk) + & (c2*state_coeffs(5,kk) + & c3*state_coeffs(9,kk)*SQ)*SQ, & c0, KMT(:,:,bid) >= k) endif !----------------------------------------------------------------------- ! ! linear EOS ! !----------------------------------------------------------------------- case (state_type_linear) if (present(RHOOUT )) RHOOUT = bet*(SALTK - S_leos_ref) - & alf*(TEMPK - T_leos_ref) if (present(RHOFULL)) RHOFULL = rho_leos_ref + & bet*(SALTK - S_leos_ref) - & alf*(TEMPK - T_leos_ref) if (present(DRHODT )) DRHODT = -alf if (present(DRHODS )) DRHODS = bet !----------------------------------------------------------------------- end select !----------------------------------------------------------------------- !EOC end subroutine state subroutine state_singlept(TEMPK,SALTK,DEPTHM,RHOFULL) ! !DESCRIPTION: ! Evaluate state for a single point- for given T,S, and depth, compute and return density !----------------------------------------------------------------------- ! This routine was originally called ovf_state and was located in the ! overflows.F90 module. ! ! calculate density of water at level k from equation of state ! at a single point. ! rho = rho(depth,pot.temperature,salinity). the density can be ! returned as a perturbation (RHOOUT) or as the full density ! (RHOFULL). ! ! This routine also computes derivatives of density with respect ! to temperature and salinity at level k from equation of state ! if requested (ie the optional arguments are present) ! ! * normally kk = k, but kk may = k+1 for adiabatic displacement ! to next lower level ! ! this routine supports three forms of the EOS (see module comments ! at the top) ! !----------------------------------------------------------------------- ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! input variables: !----------------------------------------------------------------------- real (r8) :: & TEMPK, & ! temperature at depthm SALTK ! salinity at depthm real (r8) :: & DEPTHM ! depth in centimeters !----------------------------------------------------------------------- ! output variables: !----------------------------------------------------------------------- real (r8) :: & RHOFULL ! full water density (g/cm3) at level k after adiab disp !----------------------------------------------------------------------- ! local variables: !----------------------------------------------------------------------- save real (r8) :: & TQ,SQ, & ! adjusted T,S BULK_MOD, & ! Bulk modulus RHO0, & ! density at the surface SQR,DENOMK, & ! work arrays WORK1, WORK2, & ! work arrays WORK3, WORK4, T2 ! MWJF numerator coefficients including pressure real (r8) :: & mwjfnums0t0, mwjfnums0t1, mwjfnums0t2, mwjfnums0t3, & mwjfnums1t0, mwjfnums1t1, mwjfnums2t0, & mwjfdens0t0, mwjfdens0t1, mwjfdens0t2, mwjfdens0t3, mwjfdens0t4, & mwjfdens1t0, mwjfdens1t1, mwjfdens1t3, & mwjfdensqt0, mwjfdensqt2 real (r8) :: p, p2 ! temporary pressure scalars logical (log_kind), parameter :: prnt = .false. !----------------------------------------------------------------------- if(prnt) then write(stdout,1234) TEMPK,SALTK*c1000,DEPTHM 1234 format(' Density test input TEMPK,SALTK,DEPTHM = ', & 2(f8.4,2x),f10.1) call shr_sys_flush(stdout) endif if (lccsm_control_compatible) then TQ = TEMPK SQ = SALTK*c1000 ! input salt in fraction; convert to ppt else TQ = max(-2.0_r8,TEMPK) SQ = max(c0,SALTK*c1000) ! input salt in fraction; convert to ppt endif !----------------------------------------------------------------------- ! now compute density or expansion coefficients !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! McDougall, Wright, Jackett, and Feistel EOS ! test value : rho = 1.033213387 for ! S = 35.0 PSU, theta = 20.0, pressz = 200.0 !----------------------------------------------------------------------- ! depth input as cm; m internal, so convert p = depthm/100.0_r8 ! first calculate numerator of MWJF density [P_1(S,T,p)] mwjfnums0t0 = mwjfnp0s0t0 + p*(mwjfnp1s0t0 + p*mwjfnp2s0t0) mwjfnums0t1 = mwjfnp0s0t1 mwjfnums0t2 = mwjfnp0s0t2 + p*(mwjfnp1s0t2 + p*mwjfnp2s0t2) mwjfnums0t3 = mwjfnp0s0t3 SQR = sqrt(SQ) ! first calculate numerator of MWJF density [P_1(S,T,p)] mwjfnums0t0 = mwjfnp0s0t0 + p*(mwjfnp1s0t0 + p*mwjfnp2s0t0) mwjfnums0t1 = mwjfnp0s0t1 mwjfnums0t2 = mwjfnp0s0t2 + p*(mwjfnp1s0t2 + p*mwjfnp2s0t2) mwjfnums0t3 = mwjfnp0s0t3 mwjfnums1t0 = mwjfnp0s1t0 + p*mwjfnp1s1t0 mwjfnums1t1 = mwjfnp0s1t1 mwjfnums2t0 = mwjfnp0s2t0 WORK1 = mwjfnums0t0 + TQ * (mwjfnums0t1 + TQ * (mwjfnums0t2 + & mwjfnums0t3 * TQ)) + SQ * (mwjfnums1t0 + & mwjfnums1t1 * TQ + mwjfnums2t0 * SQ) ! now calculate denominator of MWJF density [P_2(S,T,p)] mwjfdens0t0 = mwjfdp0s0t0 + p*mwjfdp1s0t0 mwjfdens0t1 = mwjfdp0s0t1 + p**3 * mwjfdp3s0t1 mwjfdens0t2 = mwjfdp0s0t2 mwjfdens0t3 = mwjfdp0s0t3 + p**2 * mwjfdp2s0t3 mwjfdens0t4 = mwjfdp0s0t4 mwjfdens1t0 = mwjfdp0s1t0 mwjfdens1t1 = mwjfdp0s1t1 mwjfdens1t3 = mwjfdp0s1t3 mwjfdensqt0 = mwjfdp0sqt0 mwjfdensqt2 = mwjfdp0sqt2 WORK2 = mwjfdens0t0 + TQ * (mwjfdens0t1 + TQ * (mwjfdens0t2 + & TQ * (mwjfdens0t3 + mwjfdens0t4 * TQ))) + & SQ * (mwjfdens1t0 + TQ * (mwjfdens1t1 + TQ*TQ*mwjfdens1t3)+ & SQR * (mwjfdensqt0 + TQ*TQ*mwjfdensqt2)) DENOMK = c1/WORK2 ! output density in g/cm3 RHOFULL = WORK1*DENOMK if(prnt) then write(stdout,2345) RHOFULL 2345 format(' state_singlept: RHOFULL output = ',f10.5) call shr_sys_flush(stdout) endif !---------------------------------------------------------------------- !EOC end subroutine state_singlept !*********************************************************************** !BOP ! !IROUTINE: ref_pressure ! !INTERFACE: function ref_pressure(k) ! !DESCRIPTION: ! This function returns a reference pressure at level k. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: integer (int_kind), intent(in) :: & k ! vertical level index ! !OUTPUT PARAMETERS: real (r8) :: & ref_pressure ! reference pressure at level k !EOP !BOC !----------------------------------------------------------------------- ! ! return pre-computed reference pressure at level k ! !----------------------------------------------------------------------- ref_pressure = pressz(k) !----------------------------------------------------------------------- !EOC end function ref_pressure !*********************************************************************** !BOP ! !IROUTINE: init_state ! !INTERFACE: subroutine init_state ! !DESCRIPTION: ! Initializes eos choice and initializes coefficients for ! equation of state. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! namelist for input choice and coefficient file ! !----------------------------------------------------------------------- integer (int_kind) :: & k, &! vertical loop index nu, &! unit number for coeff input file m, &! dummy loop index rec_length, &! record length for input polynomial file nml_error ! namelist i/o error flag character (char_len) :: & state_choice, &! choice of which EOS to use state_file, &! file containing polynomial eos coeffs state_range_opt ! option for checking valid T,S range namelist /state_nml/ state_choice, state_file, & state_range_opt, state_range_freq !----------------------------------------------------------------------- ! ! read input choice for EOS type ! !----------------------------------------------------------------------- state_itype = state_type_mwjf state_file = 'internal' state_range_iopt = state_range_ignore state_range_freq = 100000 if (my_task == master_task) then open (nml_in, file=nml_filename, status='old',iostat=nml_error) if (nml_error /= 0) then nml_error = -1 else nml_error = 1 endif do while (nml_error > 0) read(nml_in, nml=state_nml,iostat=nml_error) end do if (nml_error == 0) close(nml_in) endif call broadcast_scalar(nml_error, master_task) if (nml_error /= 0) then call exit_POP(sigAbort,'ERROR reading state_nml') endif if (my_task == master_task) then write(stdout,blank_fmt) write(stdout,ndelim_fmt) write(stdout,blank_fmt) write(stdout,'(a25)') 'Equation of state options' write(stdout,blank_fmt) write(stdout,delim_fmt) select case (state_choice(1:4)) case ('jmcd') state_itype = state_type_jmcd write(stdout,'(a29)') 'Using Jackett & McDougall EOS' case ('mwjf') state_itype = state_type_mwjf write(stdout,'(a48)') & 'Using McDougall, Wright ,Jackett and Feistel EOS' case ('poly') state_itype = state_type_polynomial write(stdout,'(a20)') 'Using polynomial EOS' case ('line') state_itype = state_type_linear write(stdout,'(a16)') 'Using linear EOS' case default state_itype = -1000 end select select case (state_range_opt) case ('ignore') state_range_iopt = state_range_ignore case ('check') state_range_iopt = state_range_check case ('enforce') state_range_iopt = state_range_enforce case default state_range_iopt = -1000 end select endif call broadcast_scalar(state_itype, master_task) call broadcast_scalar(state_file, master_task) call broadcast_scalar(state_range_iopt, master_task) call broadcast_scalar(state_range_freq, master_task) if (state_itype == -1000) then call exit_POP(sigAbort, & 'unknown choice for equation of state type') endif if (state_range_iopt == -1000) then call exit_POP(sigAbort,'unknown choice for state range option') endif !----------------------------------------------------------------------- ! ! write eos options to stdout ! !----------------------------------------------------------------------- if (my_task == master_task) then select case (state_range_iopt) case (state_range_ignore) write(stdout,'(a30)') 'No checking of valid T,S range' case (state_range_check) write(stdout,'(a30,i6,a7)') 'Valid T,S range checked every ', & state_range_freq,' steps.' case (state_range_enforce) write(stdout,'(a37)') 'EOS computed as if in valid T,S range' end select if (state_itype == state_type_polynomial) then if (state_file == 'internal') then write(stdout,'(a39)') & 'Calculating EOS coefficients internally' else write(stdout,*) & 'Reading EOS coefficients from file: ', trim(state_file) endif endif endif !----------------------------------------------------------------------- ! ! calculate pressure as function of depth ! !----------------------------------------------------------------------- do k=1,km pressz(k) = pressure(zt(k)*mpercm) end do !----------------------------------------------------------------------- ! ! initialize various arrays, constants for each state type ! !----------------------------------------------------------------------- select case (state_itype) !----------------------------------------------------------------------- ! ! McDougall, Wright, Jackett, and Feistel ! !----------------------------------------------------------------------- case (state_type_mwjf) tmin = -2.0_r8 ! limited on the low end tmax = 999.0_r8 ! unlimited on the high end smin = 0.0_r8 ! limited on the low end smax = 0.999_r8 ! unlimited on the high end !----------------------------------------------------------------------- ! ! Jackett and McDougall ! !----------------------------------------------------------------------- case (state_type_jmcd) tmin = -2.0_r8 ! valid pot. temp. range for level k tmax = 40.0_r8 smin = 0.0_r8 ! valid salinity range for level k smax = 0.042_r8 !----------------------------------------------------------------------- ! ! polynomial - initialize polynomial coefficients and ranges ! !----------------------------------------------------------------------- case (state_type_polynomial) allocate (to(km), so(km), sigo(km)) allocate (state_coeffs(9,km)) if (state_file == 'internal') then call init_state_coeffs else call get_unit(nu) if (my_task == master_task) then inquire(iolength = rec_length) to open(nu,file=state_file,access='direct',form='unformatted',& recl=rec_length,status='unknown') read(nu,rec=1) to read(nu,rec=2) so read(nu,rec=3) sigo do m = 1,9 read(nu,rec=m+3) state_coeffs(m,:) enddo if (state_range_iopt /= state_range_ignore) then read(nu,rec=13) tmin read(nu,rec=14) tmax read(nu,rec=15) smin read(nu,rec=16) smax endif close (nu) endif call release_unit(nu) call broadcast_array(to, master_task) call broadcast_array(so, master_task) call broadcast_array(sigo, master_task) do m = 1,9 call broadcast_array(state_coeffs(m,:), master_task) enddo if (state_range_iopt /= state_range_ignore) then call broadcast_array(tmin, master_task) call broadcast_array(tmax, master_task) call broadcast_array(smin, master_task) call broadcast_array(smax, master_task) endif endif !*** !*** convert smin,smax to model units !*** if (state_range_iopt /= state_range_ignore) then smin = smin*ppt_to_salt smax = smax*ppt_to_salt endif !----------------------------------------------------------------------- ! ! linear EOS - intialize valid range ! !----------------------------------------------------------------------- case (state_type_linear) tmin = -2.0_r8 ! valid pot. temp. range for level k tmax = 40.0_r8 smin = 0.0_r8 ! valid salinity range for level k smax = 0.042_r8 end select !----------------------------------------------------------------------- ! ! backwards compatability flag for certain CMIP5 runs ! !----------------------------------------------------------------------- lccsm_control_compatible = registry_match('lccsm_control_compatible') !----------------------------------------------------------------------- !EOC end subroutine init_state !*********************************************************************** !BOP ! !IROUTINE: init_state_coeffs ! !INTERFACE: subroutine init_state_coeffs ! !DESCRIPTION: ! This routine calculates coefficients for the polynomial equation of ! state option. This routine calculates the 9 coefficients of a third ! order (in temperature and salinity) polynomial approximation to the ! equation of state for sea water. The coefficients are calculated by ! first sampling a range of temperature and salinity at each depth. ! The density is then computed at each of the sampled points using the ! full UNESCO equation of state. A least squares method is used to ! fit the coefficients of the polynomial to the sampled points. More ! specifically, the densities calculated from the polynomial ! formula are in the form of sigma anomalies. The method is ! taken from that described by Bryan and Cox (1972). ! By default, the program uses the equation of state set by the ! Joint Panel on Oceanographic Tables and Standards (UNESCO, 1981) ! an described by Gill (1982). ! ! This was originall adopted from the GFDL MOM model and extensively ! modified to improve readability (February (1999). ! ! Additionally, improvements to the accuracy of the EOS have been ! added (March 1999). There are two improvements: (a) A more ! accurate depth-to-pressure conversion based on the Levitus 1994 ! climatology, and (b) Use of a more recent formula for the ! computation of potential temperature (Bryden, 1973) in place of an ! older algorithm (Fofonoff, 1962). See Dukowicz (2000). ! ! References: ! ! Bryan, K. and M. Cox, An approximate equation of state ! for numerical models of ocean circulation, J. Phys. ! Oceanogr., 2, 510-514, 1972. ! ! Bryden, H.L., New polynomials for thermal expansion, adiabatic ! temperature gradient and potential temperature of sea water, ! Deap-Sea Res., 20, 401-408, 1973. ! ! Dukowicz, J. K., 2000: Reduction of Pressure and Pressure ! Gradient Errors in Ocean Simulations, J. Phys. Oceanogr., ! submitted. ! ! Fofonoff, N., The Sea: Vol 1, (ed. M. Hill). Interscience, ! New York, 1962, pp 3-30. ! ! Gill, A., Atmosphere-Ocean Dynamics: International Geophysical ! Series No. 30. Academic Press, London, 1982, pp 599-600. ! ! Hanson, R., and C. Lawson, Extensions and applications of the ! Householder algorithm for solving linear least squares ! problems. Math. Comput., 23, 787-812, 1969. ! ! UNESCO, 10th report of the joint panel on oceanographic tables ! and standards. UNESCO Tech. Papers in Marine Sci. No. 36, ! Paris, 1981. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind), parameter :: & nsample_salt = 5, &! number of pts to sample in S nsample_temp = 2*nsample_salt, &! number of pts in sample in T nsample_all = nsample_salt*nsample_temp ! total sample points integer (int_kind) :: & i,j,k,n, &! dummy loop counters nsample, &! sample number itmax, &! max iters for least squares nlines, &! number of data lines in km output nwords_last_line ! num data words in last line of km output real (r8) :: & dtemp, dsalt, &! T,S increments for curve fit temp_sample, &! sample temperature in temp range salt_sample, &! sample salinity in salt range density, &! density in MKS units tanom, sanom, &! T,S anomalies (from level average) enorm ! norm of least squares residual real (r8), dimension(:), allocatable :: & avg_theta ! average of sample pot. temps for level k real (r8), dimension(nsample_all) :: & tsamp, &! temperature at each sample point ssamp, &! salinity at each sample point thsamp, &! potential temp at each sample point sigma, &! density at each sample point sigman ! density anomaly at each sample point real (r8), dimension(nsample_all,9) :: & lsqarray ! least squares array to find coeffs real (r8), dimension(9) :: & lsqcoeffs ! polynomial coeffs returned by lsq routine !*** !*** bounds for polynomial fit using a reference model of !*** 33 levels from surface to 8000m at z=(k-1)*250 meters !*** The user should review the appropriateness of the !*** reference values set below, and modify them if the !*** intended modelling application could be expected to yield !*** temperature and salinity values outside of the ranges set !*** by default. !*** real (r8), dimension(33) :: & trefmin = (/ -2.0_r8, -2.0_r8, -2.0_r8, & & -2.0_r8, -1.0_r8, -1.0_r8, & & -1.0_r8, -1.0_r8, -1.0_r8, & & -1.0_r8, -1.0_r8, -1.0_r8, & & -1.0_r8, -1.0_r8, -1.0_r8, & & -1.0_r8, -1.0_r8, -1.0_r8, & & -1.0_r8, 0.0_r8, 0.0_r8, & & 0.0_r8, 0.0_r8, 0.0_r8, & & 0.0_r8, 0.0_r8, 0.0_r8, & & 0.0_r8, 0.0_r8, 0.0_r8, & & 0.0_r8, 0.0_r8, 0.0_r8 /) real (r8), dimension(33) :: & trefmax = (/ 29.0_r8, 19.0_r8, 14.0_r8, & & 11.0_r8, 9.0_r8, 7.0_r8, & & 7.0_r8, 7.0_r8, 7.0_r8, & & 7.0_r8, 7.0_r8, 7.0_r8, & & 7.0_r8, 7.0_r8, 7.0_r8, & & 7.0_r8, 7.0_r8, 7.0_r8, & & 7.0_r8, 7.0_r8, 7.0_r8, & & 7.0_r8, 7.0_r8, 7.0_r8, & & 7.0_r8, 7.0_r8, 7.0_r8, & & 7.0_r8, 7.0_r8, 7.0_r8, & & 7.0_r8, 7.0_r8, 7.0_r8 /) real (r8), dimension(33) :: & srefmin = (/ 28.5_r8, 33.7_r8, 34.0_r8, & & 34.1_r8, 34.2_r8, 34.4_r8, & & 34.5_r8, 34.5_r8, 34.6_r8, & & 34.6_r8, 34.6_r8, 34.6_r8, & & 34.6_r8, 34.6_r8, 34.6_r8, & & 34.6_r8, 34.6_r8, 34.6_r8, & & 34.6_r8, 34.6_r8, 34.6_r8, & & 34.6_r8, 34.6_r8, 34.7_r8, & & 34.7_r8, 34.7_r8, 34.7_r8, & & 34.7_r8, 34.7_r8, 34.7_r8, & & 34.7_r8, 34.7_r8, 34.7_r8 /) real (r8), dimension(33) :: & srefmax = (/ 37.0_r8, 36.6_r8, 35.8_r8, & & 35.7_r8, 35.3_r8, 35.1_r8, & & 35.1_r8, 35.0_r8, 35.0_r8, & & 35.0_r8, 35.0_r8, 35.0_r8, & & 35.0_r8, 35.0_r8, 35.0_r8, & & 35.0_r8, 35.0_r8, 35.0_r8, & & 35.0_r8, 35.0_r8, 35.0_r8, & & 35.0_r8, 35.0_r8, 35.0_r8, & & 35.0_r8, 35.0_r8, 35.0_r8, & & 35.0_r8, 35.0_r8, 35.0_r8, & & 35.0_r8, 35.0_r8, 35.0_r8 /) !----------------------------------------------------------------------- ! ! set the temperature and salinity ranges to be used for each ! model level when performing the polynomial fitting ! !----------------------------------------------------------------------- do k=1,km i = int(zt(k)*mpercm/250.0_r8) + 1 tmin(k) = trefmin(i) tmax(k) = trefmax(i) smin(k) = srefmin(i) smax(k) = srefmax(i) end do allocate (avg_theta(km)) !----------------------------------------------------------------------- ! ! loop over all model levels ! !----------------------------------------------------------------------- do k=1,km !----------------------------------------------------------------------- ! ! sample the temperature range with nsample_temp points ! and the salinity range with nsample_salt points and ! create an array of possible combinations of T,S samples ! !----------------------------------------------------------------------- dtemp = (tmax(k)-tmin(k)) / (nsample_temp-c1) dsalt = (smax(k)-smin(k)) / (nsample_salt-c1) nsample = 0 do i=1,nsample_temp temp_sample = tmin(k) + (i-1)*dtemp do j=1,nsample_salt nsample = nsample + 1 salt_sample = smin(k) + (j-1)*dsalt tsamp(nsample) = temp_sample ssamp(nsample) = salt_sample end do end do !----------------------------------------------------------------------- ! ! loop over the number of samples ! !----------------------------------------------------------------------- !*** !*** initialize averaging sums !*** to (k) = c0 so (k) = c0 sigo(k) = c0 avg_theta(k) = c0 do n=1,nsample_all !----------------------------------------------------------------------- ! ! calculate density (sigma) for each t,s combintion at ! this depth using full unesco equation of state ! unesco returns density (kg per m**3) ! !----------------------------------------------------------------------- call unesco(tsamp(n),ssamp(n),pressz(k),density) sigma(n) = density - 1.0e3_r8 !----------------------------------------------------------------------- ! ! calculate potential temp. from from insitu temperature, ! salinity, and pressure ! !----------------------------------------------------------------------- call potem(tsamp(n),ssamp(n),pressz(k),thsamp(n)) !----------------------------------------------------------------------- ! ! accumulate level averages and end loop over samples ! !----------------------------------------------------------------------- to (k) = to (k) + tsamp(n) so (k) = so (k) + ssamp(n) sigo(k) = sigo(k) + sigma(n) avg_theta(k) = avg_theta(k) + thsamp(n) end do ! loop over samples !----------------------------------------------------------------------- ! ! complete layer averages ! !----------------------------------------------------------------------- to (k) = to (k)/real(nsample_all) so (k) = so (k)/real(nsample_all) sigo(k) = sigo(k)/real(nsample_all) avg_theta(k) = avg_theta(k)/real(nsample_all) !----------------------------------------------------------------------- ! ! recompute average (reference) density based on level average ! values of T, S, and pressure. ! use average potential temperature in place of average temp ! !----------------------------------------------------------------------- call unesco (to(k), so(k), pressz(k), density) sigo(k) = density - 1.0e3_r8 to (k) = avg_theta(k) !----------------------------------------------------------------------- ! ! fill array for least squares routine with anomalies ! and their products (the terms in the desired cubic ! polynomial. ! !----------------------------------------------------------------------- do n=1,nsample_all tsamp(n) = thsamp(n) !*** replace temp with potential temp tanom = tsamp(n) - to(k) sanom = ssamp(n) - so(k) sigman(n) = sigma(n) - sigo(k) lsqarray(n,1) = tanom lsqarray(n,2) = sanom lsqarray(n,3) = tanom * tanom lsqarray(n,4) = tanom * sanom lsqarray(n,5) = sanom * sanom lsqarray(n,6) = tanom * tanom * tanom lsqarray(n,7) = sanom * sanom * tanom lsqarray(n,8) = tanom * tanom * sanom lsqarray(n,9) = sanom * sanom * sanom end do !----------------------------------------------------------------------- ! ! LSQL2 is a Jet Propulsion Laboratory subroutine that ! computes the least squares fit in an iterative manner for ! overdetermined systems. it is called here to iteratively ! determine polynomial coefficients for cubic polynomials ! that will best fit the density at the sampled points ! !----------------------------------------------------------------------- itmax = 4 call lsqsl2 (lsqarray, sigman, lsqcoeffs, itmax, enorm, 1.0e-7_r8) !----------------------------------------------------------------------- ! ! store coefficients ! !----------------------------------------------------------------------- state_coeffs(:,k) = lsqcoeffs !----------------------------------------------------------------------- ! ! end of loop over levels ! !----------------------------------------------------------------------- end do !----------------------------------------------------------------------- ! ! rescale some of the coefficients and reference values for correct ! units ! !----------------------------------------------------------------------- do k=1,km sigo(k) = p001 * sigo(k) so (k) = ppt_to_salt * so (k) - 0.035_r8 state_coeffs(1,k) = 1.e-3_r8 * state_coeffs(1,k) state_coeffs(3,k) = 1.e-3_r8 * state_coeffs(3,k) state_coeffs(5,k) = 1.e+3_r8 * state_coeffs(5,k) state_coeffs(6,k) = 1.e-3_r8 * state_coeffs(6,k) state_coeffs(7,k) = 1.e+3_r8 * state_coeffs(7,k) state_coeffs(9,k) = 1.e+6_r8 * state_coeffs(9,k) end do do k=1,km sigo(k) = sigo(k) * c1000 end do deallocate (avg_theta) !----------------------------------------------------------------------- !EOC end subroutine init_state_coeffs !*********************************************************************** !BOP ! !IROUTINE: potem ! !INTERFACE: subroutine potem (temp, salt, pbars, theta) ! !DESCRIPTION: ! This subroutine calculates potential temperature as a function ! of in-situ temperature, salinity, and pressure. ! ! Reference: ! Bryden, H.L., New polynomials for thermal expansion, adiabatic ! temperature gradient and potential temperature of sea water, ! Deap-Sea Res., 20, 401-408, 1973. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: real (r8), intent(in) :: & temp, &! in-situ temperature [degrees C] salt, &! salinity [per mil] pbars ! pressure [bars] ! !OUTPUT PARAMETERS: real (r8), intent(out) :: & theta ! potential temperature [degrees C] !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- real (r8) :: prss2, prss3, potmp !----------------------------------------------------------------------- ! ! compute powers of several variables ! !----------------------------------------------------------------------- prss2 = pbars*pbars prss3 = prss2*pbars !----------------------------------------------------------------------- ! ! compute potential temperature from polynomial ! !----------------------------------------------------------------------- potmp = pbars*(3.6504e-4_r8 + temp*(8.3198e-5_r8 + & temp*(-5.4065e-7_r8 + temp*4.0274e-9_r8))) + & pbars*(salt - 35.0_r8)*(1.7439e-5_r8 - & temp*2.9778e-7_r8) + prss2*(8.9309e-7_r8 + & temp*(-3.1628e-8_r8 + temp*2.1987e-10_r8)) - & 4.1057e-9_r8*prss2*(salt - 35.0_r8) + & prss3*(-1.6056e-10_r8 + temp*5.0484e-12_r8) theta = temp - potmp !----------------------------------------------------------------------- !EOC end subroutine potem !*********************************************************************** !BOP ! !IROUTINE: unesco ! !INTERFACE: subroutine unesco (temp, salt, pbars, rho) ! !DESCRIPTION: ! This subroutine calculates the density of seawater using the ! standard equation of state recommended by unesco (1981). ! ! References: ! ! Gill, A., Atmosphere-Ocean Dynamics: International Geophysical ! Series No. 30. Academic Press, London, 1982, pp 599-600. ! ! UNESCO, 10th report of the joint panel on oceanographic tables ! and standards. UNESCO Tech. Papers in Marine Sci. No. 36, ! Paris, 1981. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: real (r8), intent(in) :: & temp, &! in-situ temperature [degrees C] salt, &! salinity [practical salinity units] pbars ! pressure [bars] ! !OUTPUT PARAMETERS: real (r8), intent(out) :: & rho ! density in kilograms per cubic meter !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- real (r8) :: rw, rsto, xkw, xksto, xkstp, & tem2, tem3, tem4, tem5, slt2, st15, pbar2 !----------------------------------------------------------------------- ! ! compute powers of several variables ! !----------------------------------------------------------------------- tem2 = temp**2 tem3 = temp**3 tem4 = temp**4 tem5 = temp**5 slt2 = salt**2 st15 = salt**(1.5_r8) pbar2 = pbars**2 !----------------------------------------------------------------------- ! ! compute density ! !----------------------------------------------------------------------- rw = unt0 + unt1*temp + unt2*tem2 + unt3*tem3 + & unt4*tem4 + unt5*tem5 rsto = rw + & (uns1t0 + uns1t1*temp + uns1t2*tem2 + & uns1t3*tem3 + uns1t4*tem4)*salt & + (unsqt0 + unsqt1*temp + unsqt2*tem2)*st15 + uns2t0*slt2 xkw = 1.965221e+4_r8 + & 1.484206e+2_r8*temp - 2.327105e+0_r8*tem2 + & 1.360477e-2_r8*tem3 - 5.155288e-5_r8*tem4 xksto = xkw + & (5.46746e+1_r8 - 6.03459e-1_r8*temp + & 1.09987e-2_r8*tem2 - 6.1670e-5_r8*tem3)*salt & + (7.944e-2_r8 + & 1.6483e-2_r8 *temp - 5.3009e-4_r8*tem2)*st15 xkstp = xksto + & (3.239908e+0_r8 + 1.43713e-3_r8*temp + & 1.16092e-4_r8*tem2 - 5.77905e-7_r8*tem3)*pbars & + (2.2838e-3_r8 - 1.0981e-5_r8 *temp - & 1.6078e-6_r8 *tem2)*pbars*salt & + 1.91075e-4_r8 *pbars*st15 & + (8.50935e-5_r8 - 6.12293e-6_r8*temp + & 5.2787e-8_r8 *tem2)*pbar2 & + (-9.9348e-7_r8 + 2.0816e-8_r8 *temp + & 9.1697e-10_r8*tem2)*pbar2*salt rho = rsto / (c1 - pbars/xkstp) !----------------------------------------------------------------------- !EOC end subroutine unesco !*********************************************************************** !BOP ! !IROUTINE: pressure ! !INTERFACE: function pressure(depth) ! !DESCRIPTION: ! This function computes pressure in bars from depth in meters ! using a mean density derived from depth-dependent global ! average temperatures and salinities from Levitus 1994, and ! integrating using hydrostatic balance. ! ! References: ! ! Levitus, S., R. Burgett, and T.P. Boyer, World Ocean Atlas ! 1994, Volume 3: Salinity, NOAA Atlas NESDIS 3, US Dept. of ! Commerce, 1994. ! ! Levitus, S. and T.P. Boyer, World Ocean Atlas 1994, ! Volume 4: Temperature, NOAA Atlas NESDIS 4, US Dept. of ! Commerce, 1994. ! ! Dukowicz, J. K., 2000: Reduction of Pressure and Pressure ! Gradient Errors in Ocean Simulations, J. Phys. Oceanogr., ! submitted. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: real (r8), intent(in) :: depth ! depth in meters ! !OUTPUT PARAMETERS: real (r8) :: pressure ! pressure [bars] !EOP !BOC !----------------------------------------------------------------------- ! ! convert depth in meters to pressure in bars ! !----------------------------------------------------------------------- pressure = 0.059808_r8*(exp(-0.025_r8*depth) - c1) & + 0.100766_r8*depth + 2.28405e-7_r8*depth**2 !----------------------------------------------------------------------- !EOC end function pressure !*********************************************************************** !BOP ! !IROUTINE: lsqsl2 ! !INTERFACE: subroutine lsqsl2 (a,b,x,itmax,enorm,eps1) ! !DESCRIPTION: ! This routine is a modification of a very old ugly bit of code ! from way back in March,1968 and originally written by a ! gentleman (presumably) named R. Hanson of unknown descent. ! The routine finds a linear least squares minimization of Ax-b. ! ! !REVISION HISTORY: ! same as module ! !INPUT/OUTPUT PARAMETERS: real (r8), dimension(:,:), intent(inout) :: & a ! input matrix A with the value of each of ! the fitting functions at each sample point ! !INPUT PARAMETERS: real (r8), dimension(size(a,dim=1)), intent(in) :: & b ! data points to fit integer (int_kind), intent(in) :: & itmax ! max number of iterations for least-square real (r8), intent(in) :: & eps1 ! some sort of tolerance ! !OUTPUT PARAMETERS: real (r8), dimension(size(a,dim=2)), intent(out) :: & x ! coefficient of each fitting function ! that provides optimal fit of sampled pts real (r8), intent(out) :: & enorm ! norm of least squares residuals !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: irank,it real (r8) :: sj,dp,up,bp,aj real (r8), dimension(size(a,dim=1),size(a,dim=2)) :: & aa real (r8), dimension(size(a,dim=1) + 8*size(a,dim=2)) :: & s real (r8), dimension(size(a,dim=1) + 4*size(a,dim=2)) :: & r integer (int_kind) :: i,j,k,l,m,n,isw,ip,lm,irp1,irm1,n1,ns,k1,k2, & j1, j2, j3, j4, j5, j6, j7, j8, j9 real (r8) :: epsloc,am,a1,sp,top,top1,top2,enm1,a2 !----------------------------------------------------------------------- ! ! initialize various useful numbers ! !----------------------------------------------------------------------- k=size(a,dim=1) n=size(a,dim=2) it = 0 irank = 0 epsloc = 1.e-16_r8 isw = 1 l = min(k,n) m = max(k,n) j1 = m j2 = j1 + n j3 = j2 + n j4 = j3 + l j5 = j4 + l j6 = j5 + l j7 = j6 + l j8 = j7 + n j9 = j8 + n lm = l !----------------------------------------------------------------------- ! ! equilibrate columns of a (1)-(2). ! !----------------------------------------------------------------------- aa = a do j=1,n am = c0 do i=1,k am = max(am,abs(a(i,j))) end do !*** s(m+n+1)-s(m+2n) contains scaling for output variables. s(j2+j) = c1/am do i=1,k a(i,j) = a(i,j)*s(j2+j) end do end do !----------------------------------------------------------------------- ! ! compute column lengths ! !----------------------------------------------------------------------- do j=1,n s(j7+j) = s(j2+j) end do !----------------------------------------------------------------------- ! ! s(m+1)-s(m+ n) contains variable permutations. ! set permutation to identity. ! !----------------------------------------------------------------------- do j=1,n s(j1+j) = j end do !----------------------------------------------------------------------- ! ! begin elimination on the matrix a with orthogonal matrices . ! ip=pivot row ! !----------------------------------------------------------------------- do ip=1,lm dp = c0 k2 = ip do j=ip,n sj = c0 do i=ip,k sj = sj + a(i,j)**2 end do if (sj >= dp) then dp = sj k2 = j endif end do !----------------------------------------------------------------------- ! ! maximize (sigma)**2 by column interchange. ! exchange columns if necessary. ! !----------------------------------------------------------------------- if (k2 /= ip) then do i=1,k a1 = a(i,ip) a(i,ip) = a(i,k2) a(i,k2) = a1 end do a1 = s(j1+k2) s(j1+k2) = s(j1+ip) s(j1+ip) = a1 a1 = s(j7+k2) s(j7+k2) = s(j7+ip) s(j7+ip) = a1 endif if (ip /= 1) then a1 = c0 do i=1,ip-1 a1 = a1 + a(i,ip)**2 end do if (a1 > c0) go to 190 end if if (dp > c0) go to 200 !----------------------------------------------------------------------- ! ! test for rank deficiency. ! !----------------------------------------------------------------------- 190 if (sqrt(dp/a1) > eps1) go to 200 irank = ip-1 go to 260 !----------------------------------------------------------------------- ! ! (eps1) rank is deficient. ! !----------------------------------------------------------------------- 200 sp = sqrt(dp) !----------------------------------------------------------------------- ! ! begin front elimination on column ip. ! sp=sqroot(sigma**2). ! !----------------------------------------------------------------------- bp = c1/(dp+sp*abs(a(ip,ip))) !----------------------------------------------------------------------- ! ! store beta in s(3n+1)-s(3n+l). ! !----------------------------------------------------------------------- if (ip == k) bp = c0 r(k+2*n+ip) = bp up = sign(sp + abs(a(ip,ip)), a(ip,ip)) if (ip < k) then if (ip < n) then do j=ip+1,n sj = c0 do i=ip+1,k sj=sj+a(i,j)*a(i,ip) end do sj = sj + up*a(ip,j) sj = bp*sj !*** sj=yj now do i=ip+1,k a(i,j) = a(i,j) - a(i,ip)*sj end do a(ip,j) = a(ip,j) - sj*up end do endif a(ip,ip) = -sign(sp,a(ip,ip)) r(k+3*n+ip) = up endif end do irank = lm 260 irp1 = irank+1 irm1 = irank-1 if (irank == 0 .or. irank == n) go to 360 !----------------------------------------------------------------------- ! ! begin back processing for rank deficiency case ! if irank is less than n. ! !----------------------------------------------------------------------- do j=1,n l = min(j,irank) !*** unscale columns for rank deficient matrices do i=1,l a(i,j) = a(i,j)/s(j7+j) end do s(j7+j) = c1 s(j2+j) = c1 end do ip = irank 300 sj = c0 do j=irp1,n sj = sj + a(ip,j)**2 end do sj = sj + a(ip,ip)**2 aj = sqrt(sj) up = sign(aj+abs(a(ip,ip)), a(ip,ip)) !----------------------------------------------------------------------- ! ! ip th element of u vector calculated. ! bp = 2/length of u squared. ! !----------------------------------------------------------------------- bp = c1/(sj+abs(a(ip,ip))*aj) if (ip-1 > 0) then do i=1,ip-1 dp = a(i,ip)*up do j=irp1,n dp = dp + a(i,j)*a(ip,j) end do dp = dp/(sj+abs(a(ip,ip))*aj) !*** calc. (aj,u), where aj=jth row of a a(i,ip) = a(i,ip) - up*dp !*** modify array a. do j=irp1,n a(i,j) = a(i,j) - a(ip,j)*dp end do end do endif a(ip,ip) = -sign(aj,a(ip,ip)) !----------------------------------------------------------------------- ! ! calc. modified pivot. ! save beta and ip th element of u vector in r array. ! !----------------------------------------------------------------------- r(k+ip) = bp r(k+n+ip) = up !----------------------------------------------------------------------- ! ! test for end of back processing. ! !----------------------------------------------------------------------- if (ip-1 > 0) then ip=ip-1 go to 300 endif 360 continue do j=1,k r(j) = b(j) end do !----------------------------------------------------------------------- ! ! set initial x vector to zero. ! !----------------------------------------------------------------------- do j=1,n x(j) = c0 end do if (irank == 0) go to 690 !----------------------------------------------------------------------- ! ! apply q to rt. hand side. ! !----------------------------------------------------------------------- 390 do ip=1,irank sj = r(k+3*n+ip)*r(ip) if (ip+1 <= k) then do i=ip+1,k sj = sj + a(i,ip)*r(i) end do endif bp = r(k+2*n+ip) if (ip+1 <= k) then do i=ip+1,k r(i) = r(i) - bp*a(i,ip)*sj end do endif r(ip) = r(ip) - bp*r(k+3*n+ip)*sj end do do j=1,irank s(j) = r(j) end do enorm = c0 if (irp1.gt.k) go to 510 do j=irp1,k enorm = enorm + r(j)**2 end do enorm = sqrt(enorm) go to 510 460 do j=1,n sj = c0 ip = s(j1+j) do i=1,k sj = sj + r(i)*aa(i,ip) end do !*** apply to rt. hand side. apply scaling. r(k+n+j) = sj*s(j2+ip) end do s(1) = r(k+n+1)/a(1,1) if (n /= 1) then do j=2,n sj = c0 do i=1,j-1 sj = sj + a(i,j)*s(i) end do s(j) = (r(k+j+n)-sj)/a(j,j) end do endif !----------------------------------------------------------------------- ! ! entry to continue iterating. solves tz = c = 1st irank ! components of qb . ! !----------------------------------------------------------------------- 510 s(irank) = s(irank)/a(irank,irank) if (irm1 /= 0) then do j=1,irm1 n1 = irank-j sj = 0. do i=n1+1,irank sj = sj + a(n1,i)*s(i) end do s(n1) = (s(n1)-sj)/a(n1,n1) end do endif !----------------------------------------------------------------------- ! ! z calculated. compute x = sz. ! !----------------------------------------------------------------------- if (irank /= n) then do j=irp1,n s(j) = c0 end do do i=1,irank sj = r(k+n+i)*s(i) do j=irp1,n sj = sj + a(i,j)*s(j) end do do j=irp1,n s(j) = s(j) - a(i,j)*r(k+i)*sj end do s(i) = s(i) - r(k+i)*r(k+n+i)*sj end do endif !----------------------------------------------------------------------- ! ! increment for x of minimal length calculated. ! !----------------------------------------------------------------------- do i=1,n x(i) = x(i) + s(i) end do !----------------------------------------------------------------------- ! ! calc. sup norm of increment and residuals ! !----------------------------------------------------------------------- top1 = c0 do j=1,n top1 = max(top1,abs(s(j))*s(j7+j)) end do do i=1,k sj = c0 do j=1,n ip = s(j1+j) sj = sj + aa(i,ip)*x(j)*s(j2+ip) end do r(i) = b(i) - sj end do !----------------------------------------------------------------------- ! ! calc. sup norm of x. ! !----------------------------------------------------------------------- top = c0 do j=1,n top = max(top,abs(x(j))*s(j7+j)) end do !----------------------------------------------------------------------- ! ! compare relative change in x with tolerance eps . ! !----------------------------------------------------------------------- if (top1-top*epsloc > 0) then if (it-itmax < 0) then it = it+1 if (it /= 1) then if (top1 > p25*top2) go to 690 endif top2 = top1 if (isw == 1) then go to 390 else if (isw ==2) then go to 460 endif endif it=0 endif 690 sj=c0 do j=1,k sj = sj + r(j)**2 end do enorm = sqrt(sj) if (irank == n .and. isw == 1) then enm1 = enorm !*** save x array do j=1,n r(k+j) = x(j) end do isw = 2 it = 0 go to 460 endif !----------------------------------------------------------------------- ! ! choose best solution ! !----------------------------------------------------------------------- if (irank >= n .and. enorm > enm1) then do j=1,n x(j) = r(k+j) end do enorm = enm1 endif !----------------------------------------------------------------------- ! ! norm of ax - b located in the cell enorm . ! rearrange variables. ! !----------------------------------------------------------------------- do j=1,n s(j) = s(j1+j) end do do j=1,n do i=j,n ip = s(i) if (j == ip) go to 780 end do 780 s(i) = s(j) s(j) = j sj = x(j) x(j) = x(i) x(i) = sj end do !----------------------------------------------------------------------- ! ! scale variables. ! !----------------------------------------------------------------------- do j=1,n x(j) = x(j)*s(j2+j) end do !----------------------------------------------------------------------- !EOC end subroutine lsqsl2 !*********************************************************************** end module state_mod !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||