Changeset 12495 for NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/tests/ISOMIP+/MY_SRC/eosbn2.F90
- Timestamp:
- 2020-03-02T09:10:34+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/tests/ISOMIP+/MY_SRC/eosbn2.F90
r12353 r12495 191 191 !! *** ROUTINE eos_insitu *** 192 192 !! 193 !! ** Purpose : Compute the in situ density (ratio rho/r au0) from193 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 194 194 !! potential temperature and salinity using an equation of state 195 195 !! selected in the nameos namelist 196 196 !! 197 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - r au0 ) / rau0197 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 198 198 !! with prd in situ density anomaly no units 199 199 !! t TEOS10: CT or EOS80: PT Celsius … … 201 201 !! z depth meters 202 202 !! rho in situ density kg/m^3 203 !! r au0 reference density kg/m^3203 !! rho0 reference density kg/m^3 204 204 !! 205 205 !! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). … … 210 210 !! 211 211 !! ln_seos : simplified equation of state 212 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / r au0212 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0 213 213 !! linear case function of T only: rn_alpha<>0, other coefficients = 0 214 214 !! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 … … 216 216 !! 217 217 !! ln_leos : linear ISOMIP equation of state 218 !! prd(t,s,z) = ( -a0*(T-T0) + b0*(S-S0) ) / r au0218 !! prd(t,s,z) = ( -a0*(T-T0) + b0*(S-S0) ) / rho0 219 219 !! setup for ISOMIP linear eos 220 220 !! … … 273 273 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 274 274 ! 275 prd(ji,jj,jk) = ( zn * r1_r au0 - 1._wp ) * ztm ! density anomaly (masked)275 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 276 276 ! 277 277 END DO … … 293 293 & - rn_nu * zt * zs 294 294 ! 295 prd(ji,jj,jk) = zn * r1_r au0 * ztm ! density anomaly (masked)295 prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked) 296 296 END DO 297 297 END DO … … 308 308 ztm = tmask(ji,jj,jk) 309 309 ! 310 zn = r au0 * ( - rn_a0 * zt + rn_b0 * zs )310 zn = rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 311 311 ! 312 prd(ji,jj,jk) = zn * r1_r au0 * ztm ! density anomaly (masked)312 prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked) 313 313 END DO 314 314 END DO … … 328 328 !! *** ROUTINE eos_insitu_pot *** 329 329 !! 330 !! ** Purpose : Compute the in situ density (ratio rho/r au0) and the330 !! ** Purpose : Compute the in situ density (ratio rho/rho0) and the 331 331 !! potential volumic mass (Kg/m3) from potential temperature and 332 332 !! salinity fields using an equation of state selected in the … … 410 410 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 411 411 ! 412 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_r au0 - 1._wp ) ! density anomaly (masked)412 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rho0 - 1._wp ) ! density anomaly (masked) 413 413 END DO 414 414 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos … … 454 454 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 455 455 ! 456 prd(ji,jj,jk) = ( zn * r1_r au0 - 1._wp ) * ztm ! density anomaly (masked)456 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 457 457 END DO 458 458 END DO … … 473 473 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs & 474 474 & - rn_nu * zt * zs 475 prhop(ji,jj,jk) = ( r au0 + zn ) * ztm475 prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 476 476 ! ! density anomaly (masked) 477 477 zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 478 prd(ji,jj,jk) = zn * r1_r au0 * ztm478 prd(ji,jj,jk) = zn * r1_rho0 * ztm 479 479 ! 480 480 END DO … … 492 492 ztm = tmask(ji,jj,jk) 493 493 ! ! potential density referenced at the surface 494 zn = r au0 * ( - rn_a0 * zt + rn_b0 * zs )495 prhop(ji,jj,jk) = ( r au0 + zn ) * ztm494 zn = rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 495 prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 496 496 ! ! density anomaly (masked) 497 prd(ji,jj,jk) = zn * r1_r au0 * ztm497 prd(ji,jj,jk) = zn * r1_rho0 * ztm 498 498 ! 499 499 END DO … … 514 514 !! *** ROUTINE eos_insitu_2d *** 515 515 !! 516 !! ** Purpose : Compute the in situ density (ratio rho/r au0) from516 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 517 517 !! potential temperature and salinity using an equation of state 518 518 !! selected in the nameos namelist. * 2D field case … … 569 569 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 570 570 ! 571 prd(ji,jj) = zn * r1_r au0 - 1._wp ! unmasked in situ density anomaly571 prd(ji,jj) = zn * r1_rho0 - 1._wp ! unmasked in situ density anomaly 572 572 ! 573 573 END DO … … 589 589 & - rn_nu * zt * zs 590 590 ! 591 prd(ji,jj) = zn * r1_r au0 ! unmasked in situ density anomaly591 prd(ji,jj) = zn * r1_rho0 ! unmasked in situ density anomaly 592 592 ! 593 593 END DO … … 605 605 zh = pdep (ji,jj) ! depth at the partial step level 606 606 ! 607 zn = r au0 * ( - rn_a0 * zt + rn_b0 * zs )608 ! 609 prd(ji,jj) = zn * r1_r au0 ! unmasked in situ density anomaly607 zn = rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 608 ! 609 prd(ji,jj) = zn * r1_rho0 ! unmasked in situ density anomaly 610 610 ! 611 611 END DO … … 676 676 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 677 677 ! 678 pab(ji,jj,jk,jp_tem) = zn * r1_r au0 * ztm678 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 679 679 ! 680 680 ! beta … … 697 697 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 698 698 ! 699 pab(ji,jj,jk,jp_sal) = zn / zs * r1_r au0 * ztm699 pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 700 700 ! 701 701 END DO … … 714 714 ! 715 715 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 716 pab(ji,jj,jk,jp_tem) = zn * r1_r au0 * ztm ! alpha716 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm ! alpha 717 717 ! 718 718 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 719 pab(ji,jj,jk,jp_sal) = zn * r1_r au0 * ztm ! beta719 pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm ! beta 720 720 ! 721 721 END DO … … 733 733 ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask 734 734 ! 735 zn = rn_a0 * r au0736 pab(ji,jj,jk,jp_tem) = zn * r1_r au0 * ztm ! alpha737 ! 738 zn = rn_b0 * r au0739 pab(ji,jj,jk,jp_sal) = zn * r1_r au0 * ztm ! beta735 zn = rn_a0 * rho0 736 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm ! alpha 737 ! 738 zn = rn_b0 * rho0 739 pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm ! beta 740 740 ! 741 741 END DO … … 809 809 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 810 810 ! 811 pab(ji,jj,jp_tem) = zn * r1_r au0811 pab(ji,jj,jp_tem) = zn * r1_rho0 812 812 ! 813 813 ! beta … … 830 830 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 831 831 ! 832 pab(ji,jj,jp_sal) = zn / zs * r1_r au0832 pab(ji,jj,jp_sal) = zn / zs * r1_rho0 833 833 ! 834 834 ! … … 848 848 ! 849 849 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 850 pab(ji,jj,jp_tem) = zn * r1_r au0 ! alpha850 pab(ji,jj,jp_tem) = zn * r1_rho0 ! alpha 851 851 ! 852 852 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 853 pab(ji,jj,jp_sal) = zn * r1_r au0 ! beta853 pab(ji,jj,jp_sal) = zn * r1_rho0 ! beta 854 854 ! 855 855 END DO … … 867 867 zh = pdep (ji,jj) ! depth at the partial step level 868 868 ! 869 zn = rn_a0 * r au0870 pab(ji,jj,jp_tem) = zn * r1_r au0 ! alpha871 ! 872 zn = rn_b0 * r au0873 pab(ji,jj,jp_sal) = zn * r1_r au0 ! beta869 zn = rn_a0 * rho0 870 pab(ji,jj,jp_tem) = zn * r1_rho0 ! alpha 871 ! 872 zn = rn_b0 * rho0 873 pab(ji,jj,jp_sal) = zn * r1_rho0 ! beta 874 874 ! 875 875 END DO … … 941 941 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 942 942 ! 943 pab(jp_tem) = zn * r1_r au0943 pab(jp_tem) = zn * r1_rho0 944 944 ! 945 945 ! beta … … 962 962 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 963 963 ! 964 pab(jp_sal) = zn / zs * r1_r au0964 pab(jp_sal) = zn / zs * r1_rho0 965 965 ! 966 966 ! … … 973 973 ! 974 974 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 975 pab(jp_tem) = zn * r1_r au0 ! alpha975 pab(jp_tem) = zn * r1_rho0 ! alpha 976 976 ! 977 977 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 978 pab(jp_sal) = zn * r1_r au0 ! beta978 pab(jp_sal) = zn * r1_rho0 ! beta 979 979 ! 980 980 CASE( np_leos ) !== linear ISOMIP EOS ==! … … 984 984 zh = pdep ! depth at the partial step level 985 985 ! 986 zn = rn_a0 * r au0987 pab(jp_tem) = zn * r1_r au0 ! alpha988 ! 989 zn = rn_b0 * r au0990 pab(jp_sal) = zn * r1_r au0 ! beta986 zn = rn_a0 * rho0 987 pab(jp_tem) = zn * r1_rho0 ! alpha 988 ! 989 zn = rn_b0 * rho0 990 pab(jp_sal) = zn * r1_rho0 ! beta 991 991 ! 992 992 CASE DEFAULT … … 1214 1214 !! ** Method : PE is defined analytically as the vertical 1215 1215 !! primitive of EOS times -g integrated between 0 and z>0. 1216 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - r au0 gz ) / rau0 gz - rd1216 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd 1217 1217 !! = 1/z * /int_0^z rd dz - rd 1218 1218 !! where rd is the density anomaly (see eos_rhd function) 1219 1219 !! ab_pe are partial derivatives of PE anomaly with respect to T and S: 1220 !! ab_pe(1) = - 1/(r au0 gz) * dPE/dT + drd/dT = - d(pen)/dT1221 !! ab_pe(2) = 1/(r au0 gz) * dPE/dS + drd/dS = d(pen)/dS1220 !! ab_pe(1) = - 1/(rho0 gz) * dPE/dT + drd/dT = - d(pen)/dT 1221 !! ab_pe(2) = 1/(rho0 gz) * dPE/dS + drd/dS = d(pen)/dS 1222 1222 !! 1223 1223 !! ** Action : - pen : PE anomaly given at T-points … … 1267 1267 zn = ( zn2 * zh + zn1 ) * zh + zn0 1268 1268 ! 1269 ppen(ji,jj,jk) = zn * zh * r1_r au0 * ztm1269 ppen(ji,jj,jk) = zn * zh * r1_rho0 * ztm 1270 1270 ! 1271 1271 ! alphaPE non-linear anomaly … … 1282 1282 zn = ( zn2 * zh + zn1 ) * zh + zn0 1283 1283 ! 1284 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_r au0 * ztm1284 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 1285 1285 ! 1286 1286 ! betaPE non-linear anomaly … … 1297 1297 zn = ( zn2 * zh + zn1 ) * zh + zn0 1298 1298 ! 1299 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_r au0 * ztm1299 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 1300 1300 ! 1301 1301 END DO … … 1312 1312 zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point 1313 1313 ztm = tmask(ji,jj,jk) ! tmask 1314 zn = 0.5_wp * zh * r1_r au0 * ztm1314 zn = 0.5_wp * zh * r1_rho0 * ztm 1315 1315 ! ! Potential Energy 1316 1316 ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn … … 1332 1332 zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point 1333 1333 ztm = tmask(ji,jj,jk) ! tmask 1334 zn = 0.5_wp * zh * r1_r au0 * ztm1334 zn = 0.5_wp * zh * r1_rho0 * ztm 1335 1335 ! ! Potential Energy 1336 1336 ppen(ji,jj,jk) = 0. … … 1379 1379 IF(lwm) WRITE( numond, nameos ) 1380 1380 ! 1381 r au0 = 1027.51_wp !: volumic mass of reference [kg/m3]1381 rho0 = 1027.51_wp !: volumic mass of reference [kg/m3] 1382 1382 rcp = 3974.00_wp !: heat capacity [J/K] 1383 1383 ! … … 1793 1793 WRITE(numout,*) ' ==>>> use of simplified eos: ' 1794 1794 WRITE(numout,*) ' rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT ' 1795 WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / r au0'1795 WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0' 1796 1796 WRITE(numout,*) ' with the following coefficients :' 1797 1797 WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0 … … 1810 1810 WRITE(numout,*) 1811 1811 WRITE(numout,*) ' use of linear ISOMIP eos: rhd(dT=T-(-1),dS=S-(34.2),Z) = ' 1812 WRITE(numout,*) ' [ -a0*dT + b0*dS ]/r au0'1812 WRITE(numout,*) ' [ -a0*dT + b0*dS ]/rho0' 1813 1813 WRITE(numout,*) 1814 1814 WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0 … … 1822 1822 END SELECT 1823 1823 ! 1824 r au0_rcp = rau0 * rcp1825 r1_r au0 = 1._wp / rau01824 rho0_rcp = rho0 * rcp 1825 r1_rho0 = 1._wp / rho0 1826 1826 r1_rcp = 1._wp / rcp 1827 r1_r au0_rcp = 1._wp / rau0_rcp1827 r1_rho0_rcp = 1._wp / rho0_rcp 1828 1828 ! 1829 1829 IF(lwp) THEN … … 1840 1840 IF(lwp) WRITE(numout,*) 1841 1841 IF(lwp) WRITE(numout,*) ' Associated physical constant' 1842 IF(lwp) WRITE(numout,*) ' volumic mass of reference r au0 = ', rau0 , ' kg/m^3'1843 IF(lwp) WRITE(numout,*) ' 1. / r au0 r1_rau0 = ', r1_rau0, ' m^3/kg'1842 IF(lwp) WRITE(numout,*) ' volumic mass of reference rho0 = ', rho0 , ' kg/m^3' 1843 IF(lwp) WRITE(numout,*) ' 1. / rho0 r1_rho0 = ', r1_rho0, ' m^3/kg' 1844 1844 IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin' 1845 IF(lwp) WRITE(numout,*) ' r au0 * rcp rau0_rcp = ', rau0_rcp1846 IF(lwp) WRITE(numout,*) ' 1. / ( r au0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp1845 IF(lwp) WRITE(numout,*) ' rho0 * rcp rho0_rcp = ', rho0_rcp 1846 IF(lwp) WRITE(numout,*) ' 1. / ( rho0 * rcp ) r1_rho0_rcp = ', r1_rho0_rcp 1847 1847 ! 1848 1848 END SUBROUTINE eos_init
Note: See TracChangeset
for help on using the changeset viewer.
