[2945] | 1 | function eos_nemo, t_in, psal |
---|
| 2 | ;------------------------------------------------------------------------------ |
---|
| 3 | ; |
---|
| 4 | ; Purpose: |
---|
| 5 | ; To calculate the density using the equation of state used in NEMO. |
---|
| 6 | ; This is based on the Jackett and McDougall (1994) equation of state |
---|
| 7 | ; for calculating the in situ density based on potential temperature |
---|
| 8 | ; and salinity. |
---|
| 9 | ; |
---|
| 10 | ; Inputs: |
---|
| 11 | ; temperature => 1d array potential temperature (deg C) |
---|
| 12 | ; salinity => 1d array salinity (PSU) |
---|
| 13 | ; |
---|
| 14 | ; Outputs: |
---|
| 15 | ; rhop => 1d array pot density (kg/m**3) |
---|
| 16 | ; |
---|
| 17 | ; could add use a difference reference pressure |
---|
| 18 | ; observations are in-situ temperature? |
---|
| 19 | ; |
---|
| 20 | ; Author: |
---|
| 21 | ; D. J. Lea. Dec 2006. |
---|
| 22 | ;------------------------------------------------------------------------------ |
---|
| 23 | |
---|
| 24 | zws = SQRT( ABS( psal ) ) |
---|
| 25 | |
---|
| 26 | sz=size(t_in) |
---|
| 27 | ndim=sz(0) |
---|
| 28 | NO_LEVS=sz(ndim) |
---|
| 29 | |
---|
| 30 | prhop=psal*0. |
---|
| 31 | |
---|
| 32 | ; |
---|
| 33 | for jk = 0L, NO_LEVS-1 do begin |
---|
| 34 | |
---|
| 35 | zt = T_IN(jk) |
---|
| 36 | zs = psal(jk) |
---|
| 37 | ; * depth |
---|
| 38 | ; zh = O_DEP_LEVS(jk) ;used in calculating insitu density only |
---|
| 39 | zsr= zws(jk) |
---|
| 40 | ; * compute volumic mass pure water at atm pressure |
---|
| 41 | zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt $ |
---|
| 42 | -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 |
---|
| 43 | ; * seawater volumic mass atm pressure |
---|
| 44 | zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt $ |
---|
| 45 | -4.0899e-3 ) *zt+0.824493 |
---|
| 46 | zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 |
---|
| 47 | zr4= 4.8314e-4 |
---|
| 48 | |
---|
| 49 | ; * potential volumic mass (reference to the surface) |
---|
| 50 | zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 |
---|
| 51 | |
---|
| 52 | ; * save potential volumic mass |
---|
| 53 | prhop(jk) = zrhop |
---|
| 54 | |
---|
| 55 | ; * add the compression terms |
---|
| 56 | ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 |
---|
| 57 | zbw= ( 1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 |
---|
| 58 | zb = zbw + ze * zs |
---|
| 59 | |
---|
| 60 | zd = -2.042967e-2 |
---|
| 61 | zc = (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 |
---|
| 62 | zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 |
---|
| 63 | za = ( zd*zsr + zc ) *zs + zaw |
---|
| 64 | |
---|
| 65 | zb1= (-0.1909078*zt+7.390729 ) *zt-55.87545 |
---|
| 66 | za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 |
---|
| 67 | zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 |
---|
| 68 | zk0= ( zb1*zsr + za1 )*zs + zkw |
---|
| 69 | |
---|
| 70 | ; ; in situ density anomaly |
---|
| 71 | ; prd(jk) = zrhop / ( 1.0 - zh / ( zk0 - zh * ( za - zh * zb ) ) ) |
---|
| 72 | ; |
---|
| 73 | endfor |
---|
| 74 | |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | return, prhop |
---|
| 78 | |
---|
| 79 | END |
---|