Changeset 5837 for branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
- Timestamp:
- 2015-10-26T15:59:39+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r4624 r5837 15 15 !! - ! 2002-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d 16 16 !! - ! 2003-08 (G. Madec) F90, free form 17 !! 3.0 ! 2006-08 (G. Madec) add tfreez function 17 !! 3.0 ! 2006-08 (G. Madec) add tfreez function (now eos_fzp function) 18 18 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 19 !! - ! 2010-10 (G. Nurser, G. Madec) add eos_alpbet used in ldfslp 19 !! - ! 2010-10 (G. Nurser, G. Madec) add alpha/beta used in ldfslp 20 !! 3.7 ! 2012-03 (F. Roquet, G. Madec) add primitive of alpha and beta used in PE computation 21 !! - ! 2012-05 (F. Roquet) add Vallis and original JM95 equation of state 22 !! - ! 2013-04 (F. Roquet, G. Madec) add eos_rab, change bn2 computation and reorganize the module 23 !! - ! 2014-09 (F. Roquet) add TEOS-10, S-EOS, and modify EOS-80 20 24 !!---------------------------------------------------------------------- 21 25 … … 23 27 !! eos : generic interface of the equation of state 24 28 !! eos_insitu : Compute the in situ density 25 !! eos_insitu_pot : Compute the insitu and surface referenced potential 26 !! volumic mass 29 !! eos_insitu_pot : Compute the insitu and surface referenced potential volumic mass 27 30 !! eos_insitu_2d : Compute the in situ density for 2d fields 28 !! eos_bn2 : Compute the Brunt-Vaisala frequency 29 !! eos_alpbet : calculates the in situ thermal/haline expansion ratio 30 !! tfreez : Compute the surface freezing temperature 31 !! bn2 : Compute the Brunt-Vaisala frequency 32 !! eos_rab : generic interface of in situ thermal/haline expansion ratio 33 !! eos_rab_3d : compute in situ thermal/haline expansion ratio 34 !! eos_rab_2d : compute in situ thermal/haline expansion ratio for 2d fields 35 !! eos_fzp_2d : freezing temperature for 2d fields 36 !! eos_fzp_0d : freezing temperature for scalar 31 37 !! eos_init : set eos parameters (namelist) 32 38 !!---------------------------------------------------------------------- 33 39 USE dom_oce ! ocean space and time domain 34 40 USE phycst ! physical constants 35 USE zdfddm ! vertical physics: double diffusion41 ! 36 42 USE in_out_manager ! I/O manager 37 43 USE lib_mpp ! MPP library 44 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 38 45 USE prtctl ! Print control 39 46 USE wrk_nemo ! Memory Allocation 47 USE lbclnk ! ocean lateral boundary conditions 40 48 USE timing ! Timing 49 USE stopar ! Stochastic T/S fluctuations 50 USE stopts ! Stochastic T/S fluctuations 41 51 42 52 IMPLICIT NONE … … 47 57 MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 48 58 END INTERFACE 49 INTERFACE bn2 50 MODULE PROCEDURE eos_bn2 59 ! 60 INTERFACE eos_rab 61 MODULE PROCEDURE rab_3d, rab_2d, rab_0d 51 62 END INTERFACE 52 53 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules 54 PUBLIC eos_init ! called by istate module 55 PUBLIC bn2 ! called by step module 56 PUBLIC eos_alpbet ! called by ldfslp module 57 PUBLIC tfreez ! called by sbcice_... modules 58 59 ! !!* Namelist (nameos) * 60 INTEGER , PUBLIC :: nn_eos !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 61 REAL(wp), PUBLIC :: rn_alpha !: thermal expension coeff. (linear equation of state) 62 REAL(wp), PUBLIC :: rn_beta !: saline expension coeff. (linear equation of state) 63 64 REAL(wp), PUBLIC :: ralpbet !: alpha / beta ratio 63 ! 64 INTERFACE eos_fzp 65 MODULE PROCEDURE eos_fzp_2d, eos_fzp_0d 66 END INTERFACE 67 ! 68 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules 69 PUBLIC bn2 ! called by step module 70 PUBLIC eos_rab ! called by ldfslp, zdfddm, trabbl 71 PUBLIC eos_pt_from_ct ! called by sbcssm 72 PUBLIC eos_fzp ! called by traadv_cen2 and sbcice_... modules 73 PUBLIC eos_pen ! used for pe diagnostics in trdpen module 74 PUBLIC eos_init ! called by istate module 75 76 ! !!* Namelist (nameos) * 77 INTEGER , PUBLIC :: nn_eos ! = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 78 LOGICAL , PUBLIC :: ln_useCT ! determine if eos_pt_from_ct is used to compute sst_m 79 80 ! !!! simplified eos coefficients 81 ! default value: Vallis 2006 82 REAL(wp) :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff. 83 REAL(wp) :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff. 84 REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2 85 REAL(wp) :: rn_lambda2 = 5.4914e-4_wp ! cabbeling coeff. in S^2 86 REAL(wp) :: rn_mu1 = 1.4970e-4_wp ! thermobaric coeff. in T 87 REAL(wp) :: rn_mu2 = 1.1090e-5_wp ! thermobaric coeff. in S 88 REAL(wp) :: rn_nu = 2.4341e-3_wp ! cabbeling coeff. in theta*salt 89 90 ! TEOS10/EOS80 parameters 91 REAL(wp) :: r1_S0, r1_T0, r1_Z0, rdeltaS 92 93 ! EOS parameters 94 REAL(wp) :: EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600 95 REAL(wp) :: EOS010 , EOS110 , EOS210 , EOS310 , EOS410 , EOS510 96 REAL(wp) :: EOS020 , EOS120 , EOS220 , EOS320 , EOS420 97 REAL(wp) :: EOS030 , EOS130 , EOS230 , EOS330 98 REAL(wp) :: EOS040 , EOS140 , EOS240 99 REAL(wp) :: EOS050 , EOS150 100 REAL(wp) :: EOS060 101 REAL(wp) :: EOS001 , EOS101 , EOS201 , EOS301 , EOS401 102 REAL(wp) :: EOS011 , EOS111 , EOS211 , EOS311 103 REAL(wp) :: EOS021 , EOS121 , EOS221 104 REAL(wp) :: EOS031 , EOS131 105 REAL(wp) :: EOS041 106 REAL(wp) :: EOS002 , EOS102 , EOS202 107 REAL(wp) :: EOS012 , EOS112 108 REAL(wp) :: EOS022 109 REAL(wp) :: EOS003 , EOS103 110 REAL(wp) :: EOS013 111 112 ! ALPHA parameters 113 REAL(wp) :: ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500 114 REAL(wp) :: ALP010 , ALP110 , ALP210 , ALP310 , ALP410 115 REAL(wp) :: ALP020 , ALP120 , ALP220 , ALP320 116 REAL(wp) :: ALP030 , ALP130 , ALP230 117 REAL(wp) :: ALP040 , ALP140 118 REAL(wp) :: ALP050 119 REAL(wp) :: ALP001 , ALP101 , ALP201 , ALP301 120 REAL(wp) :: ALP011 , ALP111 , ALP211 121 REAL(wp) :: ALP021 , ALP121 122 REAL(wp) :: ALP031 123 REAL(wp) :: ALP002 , ALP102 124 REAL(wp) :: ALP012 125 REAL(wp) :: ALP003 126 127 ! BETA parameters 128 REAL(wp) :: BET000 , BET100 , BET200 , BET300 , BET400 , BET500 129 REAL(wp) :: BET010 , BET110 , BET210 , BET310 , BET410 130 REAL(wp) :: BET020 , BET120 , BET220 , BET320 131 REAL(wp) :: BET030 , BET130 , BET230 132 REAL(wp) :: BET040 , BET140 133 REAL(wp) :: BET050 134 REAL(wp) :: BET001 , BET101 , BET201 , BET301 135 REAL(wp) :: BET011 , BET111 , BET211 136 REAL(wp) :: BET021 , BET121 137 REAL(wp) :: BET031 138 REAL(wp) :: BET002 , BET102 139 REAL(wp) :: BET012 140 REAL(wp) :: BET003 141 142 ! PEN parameters 143 REAL(wp) :: PEN000 , PEN100 , PEN200 , PEN300 , PEN400 144 REAL(wp) :: PEN010 , PEN110 , PEN210 , PEN310 145 REAL(wp) :: PEN020 , PEN120 , PEN220 146 REAL(wp) :: PEN030 , PEN130 147 REAL(wp) :: PEN040 148 REAL(wp) :: PEN001 , PEN101 , PEN201 149 REAL(wp) :: PEN011 , PEN111 150 REAL(wp) :: PEN021 151 REAL(wp) :: PEN002 , PEN102 152 REAL(wp) :: PEN012 153 154 ! ALPHA_PEN parameters 155 REAL(wp) :: APE000 , APE100 , APE200 , APE300 156 REAL(wp) :: APE010 , APE110 , APE210 157 REAL(wp) :: APE020 , APE120 158 REAL(wp) :: APE030 159 REAL(wp) :: APE001 , APE101 160 REAL(wp) :: APE011 161 REAL(wp) :: APE002 162 163 ! BETA_PEN parameters 164 REAL(wp) :: BPE000 , BPE100 , BPE200 , BPE300 165 REAL(wp) :: BPE010 , BPE110 , BPE210 166 REAL(wp) :: BPE020 , BPE120 167 REAL(wp) :: BPE030 168 REAL(wp) :: BPE001 , BPE101 169 REAL(wp) :: BPE011 170 REAL(wp) :: BPE002 65 171 66 172 !! * Substitutions … … 68 174 # include "vectopt_loop_substitute.h90" 69 175 !!---------------------------------------------------------------------- 70 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)176 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 71 177 !! $Id$ 72 178 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 82 188 !! defined through the namelist parameter nn_eos. 83 189 !! 84 !! ** Method : 3 cases: 85 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. 86 !! the in situ density is computed directly as a function of 87 !! potential temperature relative to the surface (the opa t 88 !! variable), salt and pressure (assuming no pressure variation 89 !! along geopotential surfaces, i.e. the pressure p in decibars 90 !! is approximated by the depth in meters. 91 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 92 !! with pressure p decibars 93 !! potential temperature t deg celsius 94 !! salinity s psu 95 !! reference volumic mass rau0 kg/m**3 96 !! in situ volumic mass rho kg/m**3 97 !! in situ density anomalie prd no units 98 !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, 99 !! t = 40 deg celcius, s=40 psu 100 !! nn_eos = 1 : linear equation of state function of temperature only 101 !! prd(t) = 0.0285 - rn_alpha * t 102 !! nn_eos = 2 : linear equation of state function of temperature and 103 !! salinity 104 !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 105 !! Note that no boundary condition problem occurs in this routine 106 !! as pts are defined over the whole domain. 190 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0 191 !! with prd in situ density anomaly no units 192 !! t TEOS10: CT or EOS80: PT Celsius 193 !! s TEOS10: SA or EOS80: SP TEOS10: g/kg or EOS80: psu 194 !! z depth meters 195 !! rho in situ density kg/m^3 196 !! rau0 reference density kg/m^3 197 !! 198 !! nn_eos = -1 : polynomial TEOS-10 equation of state is used for rho(t,s,z). 199 !! Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celcius, sa=35.5 g/kg 200 !! 201 !! nn_eos = 0 : polynomial EOS-80 equation of state is used for rho(t,s,z). 202 !! Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celcius, sp=35.5 psu 203 !! 204 !! nn_eos = 1 : simplified equation of state 205 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rau0 206 !! linear case function of T only: rn_alpha<>0, other coefficients = 0 207 !! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 208 !! Vallis like equation: use default values of coefficients 107 209 !! 108 210 !! ** Action : compute prd , the in situ density (no units) 109 211 !! 110 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 111 !!---------------------------------------------------------------------- 112 !! 113 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 114 ! ! 2 : salinity [psu] 115 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-] 116 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m] 117 !! 118 INTEGER :: ji, jj, jk ! dummy loop indices 119 REAL(wp) :: zt , zs , zh , zsr ! local scalars 120 REAL(wp) :: zr1, zr2, zr3, zr4 ! - - 121 REAL(wp) :: zrhop, ze, zbw, zb ! - - 122 REAL(wp) :: zd , zc , zaw, za ! - - 123 REAL(wp) :: zb1, za1, zkw, zk0 ! - - 124 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 125 !!---------------------------------------------------------------------- 126 127 ! 128 IF( nn_timing == 1 ) CALL timing_start('eos') 129 ! 130 CALL wrk_alloc( jpi, jpj, jpk, zws ) 212 !! References : Roquet et al, Ocean Modelling, in preparation (2014) 213 !! Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 214 !! TEOS-10 Manual, 2010 215 !!---------------------------------------------------------------------- 216 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 217 ! ! 2 : salinity [psu] 218 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-] 219 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m] 220 ! 221 INTEGER :: ji, jj, jk ! dummy loop indices 222 REAL(wp) :: zt , zh , zs , ztm ! local scalars 223 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 224 !!---------------------------------------------------------------------- 225 ! 226 IF( nn_timing == 1 ) CALL timing_start('eos-insitu') 131 227 ! 132 228 SELECT CASE( nn_eos ) 133 229 ! 134 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 135 !CDIR NOVERRCHK 136 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 230 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 137 231 ! 138 232 DO jk = 1, jpkm1 139 233 DO jj = 1, jpj 140 234 DO ji = 1, jpi 141 zt = pts (ji,jj,jk,jp_tem) 142 zs = pts (ji,jj,jk,jp_sal) 143 zh = pdep(ji,jj,jk) ! depth 144 zsr= zws (ji,jj,jk) ! square root salinity 145 ! 146 ! compute volumic mass pure water at atm pressure 147 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 148 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 149 ! seawater volumic mass atm pressure 150 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 151 & -4.0899e-3_wp ) *zt+0.824493_wp 152 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 153 zr4= 4.8314e-4_wp 154 ! 155 ! potential volumic mass (reference to the surface) 156 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 157 ! 158 ! add the compression terms 159 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 160 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 161 zb = zbw + ze * zs 162 ! 163 zd = -2.042967e-2_wp 164 zc = (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 165 zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 166 za = ( zd*zsr + zc ) *zs + zaw 167 ! 168 zb1= (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp 169 za1= ( ( 2.326469e-3_wp*zt+1.553190_wp) *zt-65.00517_wp ) *zt+1044.077_wp 170 zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 171 zk0= ( zb1*zsr + za1 )*zs + zkw 172 ! 173 ! masked in situ density anomaly 174 prd(ji,jj,jk) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) & 175 & - rau0 ) * r1_rau0 * tmask(ji,jj,jk) 235 ! 236 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 237 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 238 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 239 ztm = tmask(ji,jj,jk) ! tmask 240 ! 241 zn3 = EOS013*zt & 242 & + EOS103*zs+EOS003 243 ! 244 zn2 = (EOS022*zt & 245 & + EOS112*zs+EOS012)*zt & 246 & + (EOS202*zs+EOS102)*zs+EOS002 247 ! 248 zn1 = (((EOS041*zt & 249 & + EOS131*zs+EOS031)*zt & 250 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 251 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 252 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 253 ! 254 zn0 = (((((EOS060*zt & 255 & + EOS150*zs+EOS050)*zt & 256 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 257 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 258 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 259 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 260 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 261 ! 262 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 263 ! 264 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 265 ! 176 266 END DO 177 267 END DO 178 268 END DO 179 269 ! 180 CASE( 1 ) !== Linear formulation function of temperature only ==! 270 CASE( 1 ) !== simplified EOS ==! 271 ! 181 272 DO jk = 1, jpkm1 182 prd(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 273 DO jj = 1, jpj 274 DO ji = 1, jpi 275 zt = pts (ji,jj,jk,jp_tem) - 10._wp 276 zs = pts (ji,jj,jk,jp_sal) - 35._wp 277 zh = pdep (ji,jj,jk) 278 ztm = tmask(ji,jj,jk) 279 ! 280 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt & 281 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 282 & - rn_nu * zt * zs 283 ! 284 prd(ji,jj,jk) = zn * r1_rau0 * ztm ! density anomaly (masked) 285 END DO 286 END DO 183 287 END DO 184 288 ! 185 CASE( 2 ) !== Linear formulation function of temperature and salinity ==!186 DO jk = 1, jpkm1187 prd(:,:,jk) = ( rn_beta * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)188 END DO189 !190 289 END SELECT 191 290 ! 192 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos : ', ovlap=1, kdim=jpk ) 193 ! 194 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 195 ! 196 IF( nn_timing == 1 ) CALL timing_stop('eos') 291 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu : ', ovlap=1, kdim=jpk ) 292 ! 293 IF( nn_timing == 1 ) CALL timing_stop('eos-insitu') 197 294 ! 198 295 END SUBROUTINE eos_insitu … … 208 305 !! namelist parameter nn_eos. 209 306 !! 210 !! ** Method :211 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state.212 !! the in situ density is computed directly as a function of213 !! potential temperature relative to the surface (the opa t214 !! variable), salt and pressure (assuming no pressure variation215 !! along geopotential surfaces, i.e. the pressure p in decibars216 !! is approximated by the depth in meters.217 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0218 !! rhop(t,s) = rho(t,s,0)219 !! with pressure p decibars220 !! potential temperature t deg celsius221 !! salinity s psu222 !! reference volumic mass rau0 kg/m**3223 !! in situ volumic mass rho kg/m**3224 !! in situ density anomalie prd no units225 !!226 !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,227 !! t = 40 deg celcius, s=40 psu228 !!229 !! nn_eos = 1 : linear equation of state function of temperature only230 !! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t231 !! rhop(t,s) = rho(t,s)232 !!233 !! nn_eos = 2 : linear equation of state function of temperature and234 !! salinity235 !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0236 !! = rn_beta * s - rn_alpha * tn - 1.237 !! rhop(t,s) = rho(t,s)238 !! Note that no boundary condition problem occurs in this routine239 !! as (tn,sn) or (ta,sa) are defined over the whole domain.240 !!241 307 !! ** Action : - prd , the in situ density (no units) 242 308 !! - prhop, the potential volumic mass (Kg/m3) 243 309 !! 244 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 245 !! Brown and Campana, Mon. Weather Rev., 1978 246 !!---------------------------------------------------------------------- 247 !! 310 !!---------------------------------------------------------------------- 248 311 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 249 312 ! ! 2 : salinity [psu] … … 252 315 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m] 253 316 ! 254 INTEGER :: ji, jj, jk ! dummy loop indices 255 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! local scalars 256 REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0 ! - - 257 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 258 !!---------------------------------------------------------------------- 259 ! 260 IF( nn_timing == 1 ) CALL timing_start('eos-p') 261 ! 262 CALL wrk_alloc( jpi, jpj, jpk, zws ) 317 INTEGER :: ji, jj, jk, jsmp ! dummy loop indices 318 INTEGER :: jdof 319 REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars 320 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 321 REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors 322 !!---------------------------------------------------------------------- 323 ! 324 IF( nn_timing == 1 ) CALL timing_start('eos-pot') 263 325 ! 264 326 SELECT CASE ( nn_eos ) 265 327 ! 266 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 267 !CDIR NOVERRCHK 268 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 328 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 329 ! 330 ! Stochastic equation of state 331 IF ( ln_sto_eos ) THEN 332 ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 333 ALLOCATE(zn_sto(1:2*nn_sto_eos)) 334 ALLOCATE(zsign(1:2*nn_sto_eos)) 335 DO jsmp = 1, 2*nn_sto_eos, 2 336 zsign(jsmp) = 1._wp 337 zsign(jsmp+1) = -1._wp 338 END DO 339 ! 340 DO jk = 1, jpkm1 341 DO jj = 1, jpj 342 DO ji = 1, jpi 343 ! 344 ! compute density (2*nn_sto_eos) times: 345 ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 346 ! (2) for t-dt, s-ds (with the opposite fluctuation) 347 DO jsmp = 1, nn_sto_eos*2 348 jdof = (jsmp + 1) / 2 349 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 350 zt = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0 ! temperature 351 zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 352 zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity 353 ztm = tmask(ji,jj,jk) ! tmask 354 ! 355 zn3 = EOS013*zt & 356 & + EOS103*zs+EOS003 357 ! 358 zn2 = (EOS022*zt & 359 & + EOS112*zs+EOS012)*zt & 360 & + (EOS202*zs+EOS102)*zs+EOS002 361 ! 362 zn1 = (((EOS041*zt & 363 & + EOS131*zs+EOS031)*zt & 364 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 365 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 366 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 367 ! 368 zn0_sto(jsmp) = (((((EOS060*zt & 369 & + EOS150*zs+EOS050)*zt & 370 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 371 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 372 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 373 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 374 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 375 ! 376 zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 377 END DO 378 ! 379 ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 380 prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 381 DO jsmp = 1, nn_sto_eos*2 382 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 383 ! 384 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rau0 - 1._wp ) ! density anomaly (masked) 385 END DO 386 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 387 prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos 388 END DO 389 END DO 390 END DO 391 DEALLOCATE(zn0_sto,zn_sto,zsign) 392 ! Non-stochastic equation of state 393 ELSE 394 DO jk = 1, jpkm1 395 DO jj = 1, jpj 396 DO ji = 1, jpi 397 ! 398 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 399 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 400 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 401 ztm = tmask(ji,jj,jk) ! tmask 402 ! 403 zn3 = EOS013*zt & 404 & + EOS103*zs+EOS003 405 ! 406 zn2 = (EOS022*zt & 407 & + EOS112*zs+EOS012)*zt & 408 & + (EOS202*zs+EOS102)*zs+EOS002 409 ! 410 zn1 = (((EOS041*zt & 411 & + EOS131*zs+EOS031)*zt & 412 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 413 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 414 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 415 ! 416 zn0 = (((((EOS060*zt & 417 & + EOS150*zs+EOS050)*zt & 418 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 419 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 420 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 421 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 422 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 423 ! 424 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 425 ! 426 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 427 ! 428 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 429 END DO 430 END DO 431 END DO 432 ENDIF 433 434 CASE( 1 ) !== simplified EOS ==! 269 435 ! 270 436 DO jk = 1, jpkm1 271 437 DO jj = 1, jpj 272 438 DO ji = 1, jpi 273 zt = pts (ji,jj,jk,jp_tem) 274 zs = pts (ji,jj,jk,jp_sal) 275 zh = pdep(ji,jj,jk) ! depth 276 zsr= zws (ji,jj,jk) ! square root salinity 277 ! 278 ! compute volumic mass pure water at atm pressure 279 zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt & 280 & -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 281 ! seawater volumic mass atm pressure 282 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 283 & -4.0899e-3_wp ) *zt+0.824493_wp 284 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 285 zr4= 4.8314e-4_wp 286 ! 287 ! potential volumic mass (reference to the surface) 288 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 289 ! 290 ! save potential volumic mass 291 prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk) 292 ! 293 ! add the compression terms 294 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 295 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 296 zb = zbw + ze * zs 297 ! 298 zd = -2.042967e-2_wp 299 zc = (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 300 zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 301 za = ( zd*zsr + zc ) *zs + zaw 302 ! 303 zb1= ( -0.1909078_wp *zt+7.390729_wp ) *zt-55.87545_wp 304 za1= ( ( 2.326469e-3_wp*zt+1.553190_wp ) *zt-65.00517_wp ) *zt + 1044.077_wp 305 zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 306 zk0= ( zb1*zsr + za1 )*zs + zkw 307 ! 308 ! masked in situ density anomaly 309 prd(ji,jj,jk) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) & 310 & - rau0 ) * r1_rau0 * tmask(ji,jj,jk) 439 zt = pts (ji,jj,jk,jp_tem) - 10._wp 440 zs = pts (ji,jj,jk,jp_sal) - 35._wp 441 zh = pdep (ji,jj,jk) 442 ztm = tmask(ji,jj,jk) 443 ! ! potential density referenced at the surface 444 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt & 445 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs & 446 & - rn_nu * zt * zs 447 prhop(ji,jj,jk) = ( rau0 + zn ) * ztm 448 ! ! density anomaly (masked) 449 zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 450 prd(ji,jj,jk) = zn * r1_rau0 * ztm 451 ! 311 452 END DO 312 453 END DO 313 454 END DO 314 455 ! 315 CASE( 1 ) !== Linear formulation = F( temperature ) ==!316 DO jk = 1, jpkm1317 prd (:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)318 prhop(:,:,jk) = ( 1.e0_wp + prd (:,:,jk) ) * rau0 * tmask(:,:,jk)319 END DO320 !321 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==!322 DO jk = 1, jpkm1323 prd (:,:,jk) = ( rn_beta * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)324 prhop(:,:,jk) = ( 1.e0_wp + prd (:,:,jk) ) * rau0 * tmask(:,:,jk)325 END DO326 !327 456 END SELECT 328 457 ! 329 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk ) 330 ! 331 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 332 ! 333 IF( nn_timing == 1 ) CALL timing_stop('eos-p') 458 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk ) 459 ! 460 IF( nn_timing == 1 ) CALL timing_stop('eos-pot') 334 461 ! 335 462 END SUBROUTINE eos_insitu_pot … … 344 471 !! defined through the namelist parameter nn_eos. * 2D field case 345 472 !! 346 !! ** Method : 347 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. 348 !! the in situ density is computed directly as a function of 349 !! potential temperature relative to the surface (the opa t 350 !! variable), salt and pressure (assuming no pressure variation 351 !! along geopotential surfaces, i.e. the pressure p in decibars 352 !! is approximated by the depth in meters. 353 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 354 !! with pressure p decibars 355 !! potential temperature t deg celsius 356 !! salinity s psu 357 !! reference volumic mass rau0 kg/m**3 358 !! in situ volumic mass rho kg/m**3 359 !! in situ density anomalie prd no units 360 !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, 361 !! t = 40 deg celcius, s=40 psu 362 !! nn_eos = 1 : linear equation of state function of temperature only 363 !! prd(t) = 0.0285 - rn_alpha * t 364 !! nn_eos = 2 : linear equation of state function of temperature and 365 !! salinity 366 !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 367 !! Note that no boundary condition problem occurs in this routine 368 !! as pts are defined over the whole domain. 369 !! 370 !! ** Action : - prd , the in situ density (no units) 371 !! 372 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 373 !!---------------------------------------------------------------------- 374 !! 473 !! ** Action : - prd , the in situ density (no units) (unmasked) 474 !! 475 !!---------------------------------------------------------------------- 375 476 REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 376 477 ! ! 2 : salinity [psu] 377 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m]478 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 378 479 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 379 !! 380 INTEGER :: ji, jj ! dummy loop indices 381 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! temporary scalars 382 REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zmask ! - - 383 REAL(wp), POINTER, DIMENSION(:,:) :: zws 384 !!---------------------------------------------------------------------- 385 ! 386 IF( nn_timing == 1 ) CALL timing_start('eos2d') 387 ! 388 CALL wrk_alloc( jpi, jpj, zws ) 389 ! 390 480 ! 481 INTEGER :: ji, jj, jk ! dummy loop indices 482 REAL(wp) :: zt , zh , zs ! local scalars 483 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 484 !!---------------------------------------------------------------------- 485 ! 486 IF( nn_timing == 1 ) CALL timing_start('eos2d') 487 ! 391 488 prd(:,:) = 0._wp 392 489 ! 393 490 SELECT CASE( nn_eos ) 394 491 ! 395 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 396 ! 397 !CDIR NOVERRCHK 492 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 493 ! 398 494 DO jj = 1, jpjm1 399 !CDIR NOVERRCHK400 495 DO ji = 1, fs_jpim1 ! vector opt. 401 zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) ) 496 ! 497 zh = pdep(ji,jj) * r1_Z0 ! depth 498 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 499 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 500 ! 501 zn3 = EOS013*zt & 502 & + EOS103*zs+EOS003 503 ! 504 zn2 = (EOS022*zt & 505 & + EOS112*zs+EOS012)*zt & 506 & + (EOS202*zs+EOS102)*zs+EOS002 507 ! 508 zn1 = (((EOS041*zt & 509 & + EOS131*zs+EOS031)*zt & 510 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 511 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 512 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 513 ! 514 zn0 = (((((EOS060*zt & 515 & + EOS150*zs+EOS050)*zt & 516 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 517 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 518 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 519 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 520 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 521 ! 522 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 523 ! 524 prd(ji,jj) = zn * r1_rau0 - 1._wp ! unmasked in situ density anomaly 525 ! 402 526 END DO 403 527 END DO 528 ! 529 CALL lbc_lnk( prd, 'T', 1. ) ! Lateral boundary conditions 530 ! 531 CASE( 1 ) !== simplified EOS ==! 532 ! 404 533 DO jj = 1, jpjm1 405 534 DO ji = 1, fs_jpim1 ! vector opt. 406 zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask 407 zt = pts (ji,jj,jp_tem) ! interpolated T 408 zs = pts (ji,jj,jp_sal) ! interpolated S 409 zsr = zws (ji,jj) ! square root of interpolated S 410 zh = pdep (ji,jj) ! depth at the partial step level 411 ! 412 ! compute volumic mass pure water at atm pressure 413 zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt & 414 & -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 415 ! seawater volumic mass atm pressure 416 zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt & 417 & -4.0899e-3_wp ) *zt+0.824493_wp 418 zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 419 zr4 = 4.8314e-4_wp 420 ! 421 ! potential volumic mass (reference to the surface) 422 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 423 ! 424 ! add the compression terms 425 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 426 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 427 zb = zbw + ze * zs 428 ! 429 zd = -2.042967e-2_wp 430 zc = (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 431 zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt -4.721788_wp 432 za = ( zd*zsr + zc ) *zs + zaw 433 ! 434 zb1= (-0.1909078_wp *zt+7.390729_wp ) *zt-55.87545_wp 435 za1= ( ( 2.326469e-3_wp*zt+1.553190_wp ) *zt-65.00517_wp ) *zt+1044.077_wp 436 zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt & 437 & +2098.925_wp ) *zt+190925.6_wp 438 zk0= ( zb1*zsr + za1 )*zs + zkw 439 ! 440 ! masked in situ density anomaly 441 prd(ji,jj) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) - rau0 ) / rau0 * zmask 535 ! 536 zt = pts (ji,jj,jp_tem) - 10._wp 537 zs = pts (ji,jj,jp_sal) - 35._wp 538 zh = pdep (ji,jj) ! depth at the partial step level 539 ! 540 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt & 541 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 542 & - rn_nu * zt * zs 543 ! 544 prd(ji,jj) = zn * r1_rau0 ! unmasked in situ density anomaly 545 ! 442 546 END DO 443 547 END DO 444 548 ! 445 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 446 DO jj = 1, jpjm1 447 DO ji = 1, fs_jpim1 ! vector opt. 448 prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 449 END DO 450 END DO 451 ! 452 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 453 DO jj = 1, jpjm1 454 DO ji = 1, fs_jpim1 ! vector opt. 455 prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 456 END DO 457 END DO 549 CALL lbc_lnk( prd, 'T', 1. ) ! Lateral boundary conditions 458 550 ! 459 551 END SELECT 460 552 ! 461 553 IF(ln_ctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 462 554 ! 463 CALL wrk_dealloc( jpi, jpj, zws ) 464 ! 465 IF( nn_timing == 1 ) CALL timing_stop('eos2d') 555 IF( nn_timing == 1 ) CALL timing_stop('eos2d') 466 556 ! 467 557 END SUBROUTINE eos_insitu_2d 468 558 469 559 470 SUBROUTINE eos_bn2( pts, pn2 ) 471 !!---------------------------------------------------------------------- 472 !! *** ROUTINE eos_bn2 *** 473 !! 474 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the time- 475 !! step of the input arguments 476 !! 477 !! ** Method : 478 !! * nn_eos = 0 : UNESCO sea water properties 479 !! The brunt-vaisala frequency is computed using the polynomial 480 !! polynomial expression of McDougall (1987): 481 !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 482 !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 483 !! computed and used in zdfddm module : 484 !! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 485 !! * nn_eos = 1 : linear equation of state (temperature only) 486 !! N^2 = grav * rn_alpha * dk[ t ]/e3w 487 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 488 !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 489 !! The use of potential density to compute N^2 introduces e r r o r 490 !! in the sign of N^2 at great depths. We recommand the use of 491 !! nn_eos = 0, except for academical studies. 492 !! Macro-tasked on horizontal slab (jk-loop) 493 !! N.B. N^2 is set to zero at the first level (JK=1) in inidtr 494 !! and is never used at this level. 495 !! 496 !! ** Action : - pn2 : the brunt-vaisala frequency 497 !! 498 !! References : McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987. 499 !!---------------------------------------------------------------------- 500 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 501 ! ! 2 : salinity [psu] 502 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pn2 ! Brunt-Vaisala frequency [s-1] 503 !! 504 INTEGER :: ji, jj, jk ! dummy loop indices 505 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 506 #if defined key_zdfddm 507 REAL(wp) :: zds ! local scalars 508 #endif 509 !!---------------------------------------------------------------------- 510 511 ! 512 IF( nn_timing == 1 ) CALL timing_start('bn2') 513 ! 514 ! pn2 : interior points only (2=< jk =< jpkm1 ) 515 ! -------------------------- 516 ! 517 SELECT CASE( nn_eos ) 518 ! 519 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 520 DO jk = 2, jpkm1 560 SUBROUTINE rab_3d( pts, pab ) 561 !!---------------------------------------------------------------------- 562 !! *** ROUTINE rab_3d *** 563 !! 564 !! ** Purpose : Calculates thermal/haline expansion ratio at T-points 565 !! 566 !! ** Method : calculates alpha / beta at T-points 567 !! 568 !! ** Action : - pab : thermal/haline expansion ratio at T-points 569 !!---------------------------------------------------------------------- 570 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 571 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab ! thermal/haline expansion ratio 572 ! 573 INTEGER :: ji, jj, jk ! dummy loop indices 574 REAL(wp) :: zt , zh , zs , ztm ! local scalars 575 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 576 !!---------------------------------------------------------------------- 577 ! 578 IF( nn_timing == 1 ) CALL timing_start('rab_3d') 579 ! 580 SELECT CASE ( nn_eos ) 581 ! 582 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 583 ! 584 DO jk = 1, jpkm1 521 585 DO jj = 1, jpj 522 586 DO ji = 1, jpi 523 zgde3w = grav / fse3w(ji,jj,jk) 524 zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-pt 525 zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0 ! salinity anomaly (s-35) at w-pt 526 zh = fsdepw(ji,jj,jk) ! depth in meters at w-point 527 ! 528 zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & ! ratio alpha/beta 529 & - 0.203814e-03_wp ) * zt & 530 & + 0.170907e-01_wp ) * zt & 531 & + 0.665157e-01_wp & 532 & + ( - 0.678662e-05_wp * zs & 533 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs & 534 & + ( ( - 0.302285e-13_wp * zh & 535 & - 0.251520e-11_wp * zs & 536 & + 0.512857e-12_wp * zt * zt ) * zh & 537 & - 0.164759e-06_wp * zs & 538 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt & 539 & + 0.380374e-04_wp ) * zh 540 ! 541 zbeta = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & ! beta 542 & - 0.301985e-05_wp ) * zt & 543 & + 0.785567e-03_wp & 544 & + ( 0.515032e-08_wp * zs & 545 & + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs & 546 & + ( ( 0.121551e-17_wp * zh & 547 & - 0.602281e-15_wp * zs & 548 & - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh & 549 & + 0.408195e-10_wp * zs & 550 & + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt & 551 & - 0.121555e-07_wp ) * zh 552 ! 553 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 554 & * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 555 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 556 #if defined key_zdfddm 557 ! !!bug **** caution a traiter zds=dk[S]= 0 !!!! 558 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s]) 559 IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 560 rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 561 #endif 587 ! 588 zh = fsdept(ji,jj,jk) * r1_Z0 ! depth 589 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 590 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 591 ztm = tmask(ji,jj,jk) ! tmask 592 ! 593 ! alpha 594 zn3 = ALP003 595 ! 596 zn2 = ALP012*zt + ALP102*zs+ALP002 597 ! 598 zn1 = ((ALP031*zt & 599 & + ALP121*zs+ALP021)*zt & 600 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 601 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 602 ! 603 zn0 = ((((ALP050*zt & 604 & + ALP140*zs+ALP040)*zt & 605 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 606 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 607 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 608 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 609 ! 610 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 611 ! 612 pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm 613 ! 614 ! beta 615 zn3 = BET003 616 ! 617 zn2 = BET012*zt + BET102*zs+BET002 618 ! 619 zn1 = ((BET031*zt & 620 & + BET121*zs+BET021)*zt & 621 & + (BET211*zs+BET111)*zs+BET011)*zt & 622 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 623 ! 624 zn0 = ((((BET050*zt & 625 & + BET140*zs+BET040)*zt & 626 & + (BET230*zs+BET130)*zs+BET030)*zt & 627 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 628 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 629 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 630 ! 631 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 632 ! 633 pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm 634 ! 562 635 END DO 563 636 END DO 564 637 END DO 565 638 ! 566 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 567 DO jk = 2, jpkm1 568 pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 569 END DO 570 ! 571 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 572 DO jk = 2, jpkm1 573 pn2(:,:,jk) = grav * ( rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) & 574 & - rn_beta * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) ) ) & 575 & / fse3w(:,:,jk) * tmask(:,:,jk) 576 END DO 577 #if defined key_zdfddm 578 DO jk = 2, jpkm1 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 639 CASE( 1 ) !== simplified EOS ==! 640 ! 641 DO jk = 1, jpkm1 579 642 DO jj = 1, jpj 580 643 DO ji = 1, jpi 581 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 582 IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 583 rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 644 zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 645 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 646 zh = fsdept(ji,jj,jk) ! depth in meters at t-point 647 ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask 648 ! 649 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 650 pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm ! alpha 651 ! 652 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 653 pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm ! beta 654 ! 584 655 END DO 585 656 END DO 586 657 END DO 587 #endif588 END SELECT589 590 IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk )591 #if defined key_zdfddm592 IF(ln_ctl) CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk )593 #endif594 !595 IF( nn_timing == 1 ) CALL timing_stop('bn2')596 !597 END SUBROUTINE eos_bn2598 599 600 SUBROUTINE eos_alpbet( pts, palpbet, beta0 )601 !!----------------------------------------------------------------------602 !! *** ROUTINE eos_alpbet ***603 !!604 !! ** Purpose : Calculates the in situ thermal/haline expansion ratio at T-points605 !!606 !! ** Method : calculates alpha / beta ratio at T-points607 !! * nn_eos = 0 : UNESCO sea water properties608 !! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial609 !! polynomial expression of McDougall (1987).610 !! Scalar beta0 is returned = 1.611 !! * nn_eos = 1 : linear equation of state (temperature only)612 !! The ratio is undefined, so we return alpha as palpbet613 !! Scalar beta0 is returned = 0.614 !! * nn_eos = 2 : linear equation of state (temperature & salinity)615 !! The alpha/beta ratio is returned as ralpbet616 !! Scalar beta0 is returned = 1.617 !!618 !! ** Action : - palpbet : thermal/haline expansion ratio at T-points619 !! : beta0 : 1. or 0.620 !!----------------------------------------------------------------------621 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity622 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpbet ! thermal/haline expansion ratio623 REAL(wp), INTENT( out) :: beta0 ! set = 1 except with case 1 eos, rho=rho(T)624 !!625 INTEGER :: ji, jj, jk ! dummy loop indices626 REAL(wp) :: zt, zs, zh ! local scalars627 !!----------------------------------------------------------------------628 !629 IF( nn_timing == 1 ) CALL timing_start('eos_alpbet')630 !631 SELECT CASE ( nn_eos )632 !633 CASE ( 0 ) ! Jackett and McDougall (1994) formulation634 DO jk = 1, jpk635 DO jj = 1, jpj636 DO ji = 1, jpi637 zt = pts(ji,jj,jk,jp_tem) ! potential temperature638 zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35)639 zh = fsdept(ji,jj,jk) ! depth in meters640 !641 palpbet(ji,jj,jk) = &642 & ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt &643 & - 0.203814e-03_wp ) * zt &644 & + 0.170907e-01_wp ) * zt &645 & + 0.665157e-01_wp &646 & + ( - 0.678662e-05_wp * zs &647 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs &648 & + ( ( - 0.302285e-13_wp * zh &649 & - 0.251520e-11_wp * zs &650 & + 0.512857e-12_wp * zt * zt ) * zh &651 & - 0.164759e-06_wp * zs &652 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt &653 & + 0.380374e-04_wp ) * zh654 END DO655 END DO656 END DO657 beta0 = 1._wp658 !659 CASE ( 1 ) !== Linear formulation = F( temperature ) ==!660 palpbet(:,:,:) = rn_alpha661 beta0 = 0._wp662 !663 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!664 palpbet(:,:,:) = ralpbet665 beta0 = 1._wp666 658 ! 667 659 CASE DEFAULT … … 672 664 END SELECT 673 665 ! 674 IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet') 675 ! 676 END SUBROUTINE eos_alpbet 677 678 679 FUNCTION tfreez( psal, pdep ) RESULT( ptf ) 666 IF(ln_ctl) CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', & 667 & tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', ovlap=1, kdim=jpk ) 668 ! 669 IF( nn_timing == 1 ) CALL timing_stop('rab_3d') 670 ! 671 END SUBROUTINE rab_3d 672 673 SUBROUTINE rab_2d( pts, pdep, pab ) 674 !!---------------------------------------------------------------------- 675 !! *** ROUTINE rab_2d *** 676 !! 677 !! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked) 678 !! 679 !! ** Action : - pab : thermal/haline expansion ratio at T-points 680 !!---------------------------------------------------------------------- 681 REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT(in ) :: pts ! pot. temperature & salinity 682 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 683 REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT( out) :: pab ! thermal/haline expansion ratio 684 ! 685 INTEGER :: ji, jj, jk ! dummy loop indices 686 REAL(wp) :: zt , zh , zs ! local scalars 687 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 688 !!---------------------------------------------------------------------- 689 ! 690 IF( nn_timing == 1 ) CALL timing_start('rab_2d') 691 ! 692 pab(:,:,:) = 0._wp 693 ! 694 SELECT CASE ( nn_eos ) 695 ! 696 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 697 ! 698 DO jj = 1, jpjm1 699 DO ji = 1, fs_jpim1 ! vector opt. 700 ! 701 zh = pdep(ji,jj) * r1_Z0 ! depth 702 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 703 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 704 ! 705 ! alpha 706 zn3 = ALP003 707 ! 708 zn2 = ALP012*zt + ALP102*zs+ALP002 709 ! 710 zn1 = ((ALP031*zt & 711 & + ALP121*zs+ALP021)*zt & 712 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 713 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 714 ! 715 zn0 = ((((ALP050*zt & 716 & + ALP140*zs+ALP040)*zt & 717 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 718 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 719 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 720 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 721 ! 722 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 723 ! 724 pab(ji,jj,jp_tem) = zn * r1_rau0 725 ! 726 ! beta 727 zn3 = BET003 728 ! 729 zn2 = BET012*zt + BET102*zs+BET002 730 ! 731 zn1 = ((BET031*zt & 732 & + BET121*zs+BET021)*zt & 733 & + (BET211*zs+BET111)*zs+BET011)*zt & 734 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 735 ! 736 zn0 = ((((BET050*zt & 737 & + BET140*zs+BET040)*zt & 738 & + (BET230*zs+BET130)*zs+BET030)*zt & 739 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 740 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 741 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 742 ! 743 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 744 ! 745 pab(ji,jj,jp_sal) = zn / zs * r1_rau0 746 ! 747 ! 748 END DO 749 END DO 750 ! 751 CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. ) ! Lateral boundary conditions 752 CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. ) 753 ! 754 CASE( 1 ) !== simplified EOS ==! 755 ! 756 DO jj = 1, jpjm1 757 DO ji = 1, fs_jpim1 ! vector opt. 758 ! 759 zt = pts (ji,jj,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 760 zs = pts (ji,jj,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 761 zh = pdep (ji,jj) ! depth at the partial step level 762 ! 763 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 764 pab(ji,jj,jp_tem) = zn * r1_rau0 ! alpha 765 ! 766 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 767 pab(ji,jj,jp_sal) = zn * r1_rau0 ! beta 768 ! 769 END DO 770 END DO 771 ! 772 CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. ) ! Lateral boundary conditions 773 CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. ) 774 ! 775 CASE DEFAULT 776 IF(lwp) WRITE(numout,cform_err) 777 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 778 nstop = nstop + 1 779 ! 780 END SELECT 781 ! 782 IF(ln_ctl) CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', & 783 & tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' ) 784 ! 785 IF( nn_timing == 1 ) CALL timing_stop('rab_2d') 786 ! 787 END SUBROUTINE rab_2d 788 789 790 SUBROUTINE rab_0d( pts, pdep, pab ) 791 !!---------------------------------------------------------------------- 792 !! *** ROUTINE rab_0d *** 793 !! 794 !! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked) 795 !! 796 !! ** Action : - pab : thermal/haline expansion ratio at T-points 797 !!---------------------------------------------------------------------- 798 REAL(wp), DIMENSION(jpts) , INTENT(in ) :: pts ! pot. temperature & salinity 799 REAL(wp), INTENT(in ) :: pdep ! depth [m] 800 REAL(wp), DIMENSION(jpts) , INTENT( out) :: pab ! thermal/haline expansion ratio 801 ! 802 REAL(wp) :: zt , zh , zs ! local scalars 803 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 804 !!---------------------------------------------------------------------- 805 ! 806 IF( nn_timing == 1 ) CALL timing_start('rab_2d') 807 ! 808 pab(:) = 0._wp 809 ! 810 SELECT CASE ( nn_eos ) 811 ! 812 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 813 ! 814 ! 815 zh = pdep * r1_Z0 ! depth 816 zt = pts (jp_tem) * r1_T0 ! temperature 817 zs = SQRT( ABS( pts(jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 818 ! 819 ! alpha 820 zn3 = ALP003 821 ! 822 zn2 = ALP012*zt + ALP102*zs+ALP002 823 ! 824 zn1 = ((ALP031*zt & 825 & + ALP121*zs+ALP021)*zt & 826 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 827 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 828 ! 829 zn0 = ((((ALP050*zt & 830 & + ALP140*zs+ALP040)*zt & 831 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 832 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 833 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 834 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 835 ! 836 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 837 ! 838 pab(jp_tem) = zn * r1_rau0 839 ! 840 ! beta 841 zn3 = BET003 842 ! 843 zn2 = BET012*zt + BET102*zs+BET002 844 ! 845 zn1 = ((BET031*zt & 846 & + BET121*zs+BET021)*zt & 847 & + (BET211*zs+BET111)*zs+BET011)*zt & 848 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 849 ! 850 zn0 = ((((BET050*zt & 851 & + BET140*zs+BET040)*zt & 852 & + (BET230*zs+BET130)*zs+BET030)*zt & 853 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 854 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 855 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 856 ! 857 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 858 ! 859 pab(jp_sal) = zn / zs * r1_rau0 860 ! 861 ! 862 ! 863 CASE( 1 ) !== simplified EOS ==! 864 ! 865 zt = pts(jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 866 zs = pts(jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 867 zh = pdep ! depth at the partial step level 868 ! 869 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 870 pab(jp_tem) = zn * r1_rau0 ! alpha 871 ! 872 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 873 pab(jp_sal) = zn * r1_rau0 ! beta 874 ! 875 CASE DEFAULT 876 IF(lwp) WRITE(numout,cform_err) 877 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 878 nstop = nstop + 1 879 ! 880 END SELECT 881 ! 882 IF( nn_timing == 1 ) CALL timing_stop('rab_2d') 883 ! 884 END SUBROUTINE rab_0d 885 886 887 SUBROUTINE bn2( pts, pab, pn2 ) 888 !!---------------------------------------------------------------------- 889 !! *** ROUTINE bn2 *** 890 !! 891 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the 892 !! time-step of the input arguments 893 !! 894 !! ** Method : pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w 895 !! where alpha and beta are given in pab, and computed on T-points. 896 !! N.B. N^2 is set one for all to zero at jk=1 in istate module. 897 !! 898 !! ** Action : pn2 : square of the brunt-vaisala frequency at w-point 899 !! 900 !!---------------------------------------------------------------------- 901 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celcius,psu] 902 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pab ! thermal/haline expansion coef. [Celcius-1,psu-1] 903 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2] 904 ! 905 INTEGER :: ji, jj, jk ! dummy loop indices 906 REAL(wp) :: zaw, zbw, zrw ! local scalars 907 !!---------------------------------------------------------------------- 908 ! 909 IF( nn_timing == 1 ) CALL timing_start('bn2') 910 ! 911 DO jk = 2, jpkm1 ! interior points only (2=< jk =< jpkm1 ) 912 DO jj = 1, jpj ! surface and bottom value set to zero one for all in istate.F90 913 DO ji = 1, jpi 914 zrw = ( fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk) ) & 915 & / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 916 ! 917 zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 918 zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 919 ! 920 pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 921 & - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) & 922 & / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 923 END DO 924 END DO 925 END DO 926 ! 927 IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk ) 928 ! 929 IF( nn_timing == 1 ) CALL timing_stop('bn2') 930 ! 931 END SUBROUTINE bn2 932 933 934 FUNCTION eos_pt_from_ct( ctmp, psal ) RESULT( ptmp ) 935 !!---------------------------------------------------------------------- 936 !! *** ROUTINE eos_pt_from_ct *** 937 !! 938 !! ** Purpose : Compute pot.temp. from cons. temp. [Celcius] 939 !! 940 !! ** Method : rational approximation (5/3th order) of TEOS-10 algorithm 941 !! checkvalue: pt=20.02391895 Celsius for sa=35.7g/kg, ct=20degC 942 !! 943 !! Reference : TEOS-10, UNESCO 944 !! Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC) 945 !!---------------------------------------------------------------------- 946 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ctmp ! Cons. Temp [Celcius] 947 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu] 948 ! Leave result array automatic rather than making explicitly allocated 949 REAL(wp), DIMENSION(jpi,jpj) :: ptmp ! potential temperature [Celcius] 950 ! 951 INTEGER :: ji, jj ! dummy loop indices 952 REAL(wp) :: zt , zs , ztm ! local scalars 953 REAL(wp) :: zn , zd ! local scalars 954 REAL(wp) :: zdeltaS , z1_S0 , z1_T0 955 !!---------------------------------------------------------------------- 956 ! 957 IF ( nn_timing == 1 ) CALL timing_start('eos_pt_from_ct') 958 ! 959 zdeltaS = 5._wp 960 z1_S0 = 0.875_wp/35.16504_wp 961 z1_T0 = 1._wp/40._wp 962 ! 963 DO jj = 1, jpj 964 DO ji = 1, jpi 965 ! 966 zt = ctmp (ji,jj) * z1_T0 967 zs = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 ) 968 ztm = tmask(ji,jj,1) 969 ! 970 zn = ((((-2.1385727895e-01_wp*zt & 971 & - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt & 972 & + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt & 973 & + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt & 974 & + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs & 975 & +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt & 976 & + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs & 977 & -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp 978 ! 979 zd = (2.0035003456_wp*zt & 980 & -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt & 981 & + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp 982 ! 983 ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm 984 ! 985 END DO 986 END DO 987 ! 988 IF( nn_timing == 1 ) CALL timing_stop('eos_pt_from_ct') 989 ! 990 END FUNCTION eos_pt_from_ct 991 992 993 FUNCTION eos_fzp_2d( psal, pdep ) RESULT( ptf ) 994 !!---------------------------------------------------------------------- 995 !! *** ROUTINE eos_fzp *** 996 !! 997 !! ** Purpose : Compute the freezing point temperature [Celcius] 998 !! 999 !! ** Method : UNESCO freezing point (ptf) in Celcius is given by 1000 !! ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z 1001 !! checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m 1002 !! 1003 !! Reference : UNESCO tech. papers in the marine science no. 28. 1978 1004 !!---------------------------------------------------------------------- 1005 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu] 1006 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: pdep ! depth [m] 1007 REAL(wp), DIMENSION(jpi,jpj) :: ptf ! freezing temperature [Celcius] 1008 ! 1009 INTEGER :: ji, jj ! dummy loop indices 1010 REAL(wp) :: zt, zs ! local scalars 1011 !!---------------------------------------------------------------------- 1012 ! 1013 SELECT CASE ( nn_eos ) 1014 ! 1015 CASE ( -1, 1 ) !== CT,SA (TEOS-10 formulation) ==! 1016 ! 1017 DO jj = 1, jpj 1018 DO ji = 1, jpi 1019 zs= SQRT( ABS( psal(ji,jj) ) * r1_S0 ) ! square root salinity 1020 ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 1021 & - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 1022 END DO 1023 END DO 1024 ptf(:,:) = ptf(:,:) * psal(:,:) 1025 ! 1026 IF( PRESENT( pdep ) ) ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) 1027 ! 1028 CASE ( 0 ) !== PT,SP (UNESCO formulation) ==! 1029 ! 1030 ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) ) & 1031 & - 2.154996e-4_wp * psal(:,:) ) * psal(:,:) 1032 ! 1033 IF( PRESENT( pdep ) ) ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) 1034 ! 1035 CASE DEFAULT 1036 IF(lwp) WRITE(numout,cform_err) 1037 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 1038 nstop = nstop + 1 1039 ! 1040 END SELECT 1041 ! 1042 END FUNCTION eos_fzp_2d 1043 1044 FUNCTION eos_fzp_0d( psal, pdep ) RESULT( ptf ) 1045 !!---------------------------------------------------------------------- 1046 !! *** ROUTINE eos_fzp *** 1047 !! 1048 !! ** Purpose : Compute the freezing point temperature [Celcius] 1049 !! 1050 !! ** Method : UNESCO freezing point (ptf) in Celcius is given by 1051 !! ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z 1052 !! checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m 1053 !! 1054 !! Reference : UNESCO tech. papers in the marine science no. 28. 1978 1055 !!---------------------------------------------------------------------- 1056 REAL(wp), INTENT(in) :: psal ! salinity [psu] 1057 REAL(wp), INTENT(in), OPTIONAL :: pdep ! depth [m] 1058 REAL(wp) :: ptf ! freezing temperature [Celcius] 1059 ! 1060 REAL(wp) :: zs ! local scalars 1061 !!---------------------------------------------------------------------- 1062 ! 1063 SELECT CASE ( nn_eos ) 1064 ! 1065 CASE ( -1, 1 ) !== CT,SA (TEOS-10 formulation) ==! 1066 ! 1067 zs = SQRT( ABS( psal ) * r1_S0 ) ! square root salinity 1068 ptf = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 1069 & - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 1070 ptf = ptf * psal 1071 ! 1072 IF( PRESENT( pdep ) ) ptf = ptf - 7.53e-4 * pdep 1073 ! 1074 CASE ( 0 ) !== PT,SP (UNESCO formulation) ==! 1075 ! 1076 ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal ) & 1077 & - 2.154996e-4_wp * psal ) * psal 1078 ! 1079 IF( PRESENT( pdep ) ) ptf = ptf - 7.53e-4 * pdep 1080 ! 1081 CASE DEFAULT 1082 IF(lwp) WRITE(numout,cform_err) 1083 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 1084 nstop = nstop + 1 1085 ! 1086 END SELECT 1087 ! 1088 END FUNCTION eos_fzp_0d 1089 1090 1091 SUBROUTINE eos_pen( pts, pab_pe, ppen ) 1092 !!---------------------------------------------------------------------- 1093 !! *** ROUTINE eos_pen *** 1094 !! 1095 !! ** Purpose : Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points 1096 !! 1097 !! ** Method : PE is defined analytically as the vertical 1098 !! primitive of EOS times -g integrated between 0 and z>0. 1099 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd 1100 !! = 1/z * /int_0^z rd dz - rd 1101 !! where rd is the density anomaly (see eos_rhd function) 1102 !! ab_pe are partial derivatives of PE anomaly with respect to T and S: 1103 !! ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT 1104 !! ab_pe(2) = 1/(rau0 gz) * dPE/dS + drd/dS = d(pen)/dS 1105 !! 1106 !! ** Action : - pen : PE anomaly given at T-points 1107 !! : - pab_pe : given at T-points 1108 !! pab_pe(:,:,:,jp_tem) is alpha_pe 1109 !! pab_pe(:,:,:,jp_sal) is beta_pe 1110 !!---------------------------------------------------------------------- 1111 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 1112 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab_pe ! alpha_pe and beta_pe 1113 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: ppen ! potential energy anomaly 1114 ! 1115 INTEGER :: ji, jj, jk ! dummy loop indices 1116 REAL(wp) :: zt , zh , zs , ztm ! local scalars 1117 REAL(wp) :: zn , zn0, zn1, zn2 ! - - 1118 !!---------------------------------------------------------------------- 1119 ! 1120 IF( nn_timing == 1 ) CALL timing_start('eos_pen') 1121 ! 1122 SELECT CASE ( nn_eos ) 1123 ! 1124 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 1125 ! 1126 DO jk = 1, jpkm1 1127 DO jj = 1, jpj 1128 DO ji = 1, jpi 1129 ! 1130 zh = fsdept(ji,jj,jk) * r1_Z0 ! depth 1131 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 1132 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 1133 ztm = tmask(ji,jj,jk) ! tmask 1134 ! 1135 ! potential energy non-linear anomaly 1136 zn2 = (PEN012)*zt & 1137 & + PEN102*zs+PEN002 1138 ! 1139 zn1 = ((PEN021)*zt & 1140 & + PEN111*zs+PEN011)*zt & 1141 & + (PEN201*zs+PEN101)*zs+PEN001 1142 ! 1143 zn0 = ((((PEN040)*zt & 1144 & + PEN130*zs+PEN030)*zt & 1145 & + (PEN220*zs+PEN120)*zs+PEN020)*zt & 1146 & + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt & 1147 & + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000 1148 ! 1149 zn = ( zn2 * zh + zn1 ) * zh + zn0 1150 ! 1151 ppen(ji,jj,jk) = zn * zh * r1_rau0 * ztm 1152 ! 1153 ! alphaPE non-linear anomaly 1154 zn2 = APE002 1155 ! 1156 zn1 = (APE011)*zt & 1157 & + APE101*zs+APE001 1158 ! 1159 zn0 = (((APE030)*zt & 1160 & + APE120*zs+APE020)*zt & 1161 & + (APE210*zs+APE110)*zs+APE010)*zt & 1162 & + ((APE300*zs+APE200)*zs+APE100)*zs+APE000 1163 ! 1164 zn = ( zn2 * zh + zn1 ) * zh + zn0 1165 ! 1166 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm 1167 ! 1168 ! betaPE non-linear anomaly 1169 zn2 = BPE002 1170 ! 1171 zn1 = (BPE011)*zt & 1172 & + BPE101*zs+BPE001 1173 ! 1174 zn0 = (((BPE030)*zt & 1175 & + BPE120*zs+BPE020)*zt & 1176 & + (BPE210*zs+BPE110)*zs+BPE010)*zt & 1177 & + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000 1178 ! 1179 zn = ( zn2 * zh + zn1 ) * zh + zn0 1180 ! 1181 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm 1182 ! 1183 END DO 1184 END DO 1185 END DO 1186 ! 1187 CASE( 1 ) !== Vallis (2006) simplified EOS ==! 1188 ! 1189 DO jk = 1, jpkm1 1190 DO jj = 1, jpj 1191 DO ji = 1, jpi 1192 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! temperature anomaly (t-T0) 1193 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 1194 zh = fsdept(ji,jj,jk) ! depth in meters at t-point 1195 ztm = tmask(ji,jj,jk) ! tmask 1196 zn = 0.5_wp * zh * r1_rau0 * ztm 1197 ! ! Potential Energy 1198 ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 1199 ! ! alphaPE 1200 pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn 1201 pab_pe(ji,jj,jk,jp_sal) = rn_b0 * rn_mu2 * zn 1202 ! 1203 END DO 1204 END DO 1205 END DO 1206 ! 1207 CASE DEFAULT 1208 IF(lwp) WRITE(numout,cform_err) 1209 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 1210 nstop = nstop + 1 1211 ! 1212 END SELECT 1213 ! 1214 IF( nn_timing == 1 ) CALL timing_stop('eos_pen') 1215 ! 1216 END SUBROUTINE eos_pen 1217 1218 1219 SUBROUTINE eos_init 680 1220 !!---------------------------------------------------------------------- 681 1221 !! *** ROUTINE eos_init *** 682 1222 !! 683 !! ** Purpose : Compute the sea surface freezing temperature [Celcius]684 !!685 !! ** Method : UNESCO freezing point at the surface (pressure = 0???)686 !! freezing point [Celcius]=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p687 !! checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars688 !!689 !! Reference : UNESCO tech. papers in the marine science no. 28. 1978690 !!----------------------------------------------------------------------691 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]692 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: pdep ! depth [decibars]693 ! Leave result array automatic rather than making explicitly allocated694 REAL(wp), DIMENSION(jpi,jpj) :: ptf ! freezing temperature [Celcius]695 !!----------------------------------------------------------------------696 !697 ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) ) &698 & - 2.154996e-4_wp * psal(:,:) ) * psal(:,:)699 IF ( PRESENT( pdep ) ) THEN700 ptf(:,:) = ptf(:,:) - 7.53e-4_wp * pdep(:,:)701 ENDIF702 !703 END FUNCTION tfreez704 705 706 SUBROUTINE eos_init707 !!----------------------------------------------------------------------708 !! *** ROUTINE eos_init ***709 !!710 1223 !! ** Purpose : initializations for the equation of state 711 1224 !! 712 1225 !! ** Method : Read the namelist nameos and control the parameters 713 1226 !!---------------------------------------------------------------------- 714 NAMELIST/nameos/ nn_eos, rn_alpha, rn_beta 715 !!---------------------------------------------------------------------- 716 INTEGER :: ios 1227 INTEGER :: ios ! local integer 1228 !! 1229 NAMELIST/nameos/ nn_eos, ln_useCT, rn_a0, rn_b0, rn_lambda1, rn_mu1, & 1230 & rn_lambda2, rn_mu2, rn_nu 1231 !!---------------------------------------------------------------------- 717 1232 ! 718 1233 REWIND( numnam_ref ) ! Namelist nameos in reference namelist : equation of state 719 1234 READ ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 ) 720 1235 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist', lwp ) 721 1236 ! 722 1237 REWIND( numnam_cfg ) ! Namelist nameos in configuration namelist : equation of state 723 1238 READ ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 ) 724 1239 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp ) 725 1240 IF(lwm) WRITE( numond, nameos ) 1241 ! 1242 rau0 = 1026._wp !: volumic mass of reference [kg/m3] 1243 rcp = 3991.86795711963_wp !: heat capacity [J/K] 726 1244 ! 727 1245 IF(lwp) THEN ! Control print … … 731 1249 WRITE(numout,*) ' Namelist nameos : set eos parameters' 732 1250 WRITE(numout,*) ' flag for eq. of state and N^2 nn_eos = ', nn_eos 733 WRITE(numout,*) ' thermal exp. coef. (linear) rn_alpha = ', rn_alpha 734 WRITE(numout,*) ' saline exp. coef. (linear) rn_beta = ', rn_beta 1251 IF( ln_useCT ) THEN 1252 WRITE(numout,*) ' model uses Conservative Temperature' 1253 WRITE(numout,*) ' Important: model must be initialized with CT and SA fields' 1254 ELSE 1255 WRITE(numout,*) ' model does not use Conservative Temperature' 1256 ENDIF 735 1257 ENDIF 736 1258 ! 737 1259 SELECT CASE( nn_eos ) ! check option 738 1260 ! 739 CASE( 0 ) !== Jackett and McDougall (1994) formulation==!1261 CASE( -1 ) !== polynomial TEOS-10 ==! 740 1262 IF(lwp) WRITE(numout,*) 741 IF(lwp) WRITE(numout,*) ' use of Jackett & McDougall (1994) equation of state and' 742 IF(lwp) WRITE(numout,*) ' McDougall (1987) Brunt-Vaisala frequency' 743 ! 744 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 1263 IF(lwp) WRITE(numout,*) ' use of TEOS-10 equation of state (cons. temp. and abs. salinity)' 1264 ! 1265 rdeltaS = 32._wp 1266 r1_S0 = 0.875_wp/35.16504_wp 1267 r1_T0 = 1._wp/40._wp 1268 r1_Z0 = 1.e-4_wp 1269 ! 1270 EOS000 = 8.0189615746e+02_wp 1271 EOS100 = 8.6672408165e+02_wp 1272 EOS200 = -1.7864682637e+03_wp 1273 EOS300 = 2.0375295546e+03_wp 1274 EOS400 = -1.2849161071e+03_wp 1275 EOS500 = 4.3227585684e+02_wp 1276 EOS600 = -6.0579916612e+01_wp 1277 EOS010 = 2.6010145068e+01_wp 1278 EOS110 = -6.5281885265e+01_wp 1279 EOS210 = 8.1770425108e+01_wp 1280 EOS310 = -5.6888046321e+01_wp 1281 EOS410 = 1.7681814114e+01_wp 1282 EOS510 = -1.9193502195_wp 1283 EOS020 = -3.7074170417e+01_wp 1284 EOS120 = 6.1548258127e+01_wp 1285 EOS220 = -6.0362551501e+01_wp 1286 EOS320 = 2.9130021253e+01_wp 1287 EOS420 = -5.4723692739_wp 1288 EOS030 = 2.1661789529e+01_wp 1289 EOS130 = -3.3449108469e+01_wp 1290 EOS230 = 1.9717078466e+01_wp 1291 EOS330 = -3.1742946532_wp 1292 EOS040 = -8.3627885467_wp 1293 EOS140 = 1.1311538584e+01_wp 1294 EOS240 = -5.3563304045_wp 1295 EOS050 = 5.4048723791e-01_wp 1296 EOS150 = 4.8169980163e-01_wp 1297 EOS060 = -1.9083568888e-01_wp 1298 EOS001 = 1.9681925209e+01_wp 1299 EOS101 = -4.2549998214e+01_wp 1300 EOS201 = 5.0774768218e+01_wp 1301 EOS301 = -3.0938076334e+01_wp 1302 EOS401 = 6.6051753097_wp 1303 EOS011 = -1.3336301113e+01_wp 1304 EOS111 = -4.4870114575_wp 1305 EOS211 = 5.0042598061_wp 1306 EOS311 = -6.5399043664e-01_wp 1307 EOS021 = 6.7080479603_wp 1308 EOS121 = 3.5063081279_wp 1309 EOS221 = -1.8795372996_wp 1310 EOS031 = -2.4649669534_wp 1311 EOS131 = -5.5077101279e-01_wp 1312 EOS041 = 5.5927935970e-01_wp 1313 EOS002 = 2.0660924175_wp 1314 EOS102 = -4.9527603989_wp 1315 EOS202 = 2.5019633244_wp 1316 EOS012 = 2.0564311499_wp 1317 EOS112 = -2.1311365518e-01_wp 1318 EOS022 = -1.2419983026_wp 1319 EOS003 = -2.3342758797e-02_wp 1320 EOS103 = -1.8507636718e-02_wp 1321 EOS013 = 3.7969820455e-01_wp 1322 ! 1323 ALP000 = -6.5025362670e-01_wp 1324 ALP100 = 1.6320471316_wp 1325 ALP200 = -2.0442606277_wp 1326 ALP300 = 1.4222011580_wp 1327 ALP400 = -4.4204535284e-01_wp 1328 ALP500 = 4.7983755487e-02_wp 1329 ALP010 = 1.8537085209_wp 1330 ALP110 = -3.0774129064_wp 1331 ALP210 = 3.0181275751_wp 1332 ALP310 = -1.4565010626_wp 1333 ALP410 = 2.7361846370e-01_wp 1334 ALP020 = -1.6246342147_wp 1335 ALP120 = 2.5086831352_wp 1336 ALP220 = -1.4787808849_wp 1337 ALP320 = 2.3807209899e-01_wp 1338 ALP030 = 8.3627885467e-01_wp 1339 ALP130 = -1.1311538584_wp 1340 ALP230 = 5.3563304045e-01_wp 1341 ALP040 = -6.7560904739e-02_wp 1342 ALP140 = -6.0212475204e-02_wp 1343 ALP050 = 2.8625353333e-02_wp 1344 ALP001 = 3.3340752782e-01_wp 1345 ALP101 = 1.1217528644e-01_wp 1346 ALP201 = -1.2510649515e-01_wp 1347 ALP301 = 1.6349760916e-02_wp 1348 ALP011 = -3.3540239802e-01_wp 1349 ALP111 = -1.7531540640e-01_wp 1350 ALP211 = 9.3976864981e-02_wp 1351 ALP021 = 1.8487252150e-01_wp 1352 ALP121 = 4.1307825959e-02_wp 1353 ALP031 = -5.5927935970e-02_wp 1354 ALP002 = -5.1410778748e-02_wp 1355 ALP102 = 5.3278413794e-03_wp 1356 ALP012 = 6.2099915132e-02_wp 1357 ALP003 = -9.4924551138e-03_wp 1358 ! 1359 BET000 = 1.0783203594e+01_wp 1360 BET100 = -4.4452095908e+01_wp 1361 BET200 = 7.6048755820e+01_wp 1362 BET300 = -6.3944280668e+01_wp 1363 BET400 = 2.6890441098e+01_wp 1364 BET500 = -4.5221697773_wp 1365 BET010 = -8.1219372432e-01_wp 1366 BET110 = 2.0346663041_wp 1367 BET210 = -2.1232895170_wp 1368 BET310 = 8.7994140485e-01_wp 1369 BET410 = -1.1939638360e-01_wp 1370 BET020 = 7.6574242289e-01_wp 1371 BET120 = -1.5019813020_wp 1372 BET220 = 1.0872489522_wp 1373 BET320 = -2.7233429080e-01_wp 1374 BET030 = -4.1615152308e-01_wp 1375 BET130 = 4.9061350869e-01_wp 1376 BET230 = -1.1847737788e-01_wp 1377 BET040 = 1.4073062708e-01_wp 1378 BET140 = -1.3327978879e-01_wp 1379 BET050 = 5.9929880134e-03_wp 1380 BET001 = -5.2937873009e-01_wp 1381 BET101 = 1.2634116779_wp 1382 BET201 = -1.1547328025_wp 1383 BET301 = 3.2870876279e-01_wp 1384 BET011 = -5.5824407214e-02_wp 1385 BET111 = 1.2451933313e-01_wp 1386 BET211 = -2.4409539932e-02_wp 1387 BET021 = 4.3623149752e-02_wp 1388 BET121 = -4.6767901790e-02_wp 1389 BET031 = -6.8523260060e-03_wp 1390 BET002 = -6.1618945251e-02_wp 1391 BET102 = 6.2255521644e-02_wp 1392 BET012 = -2.6514181169e-03_wp 1393 BET003 = -2.3025968587e-04_wp 1394 ! 1395 PEN000 = -9.8409626043_wp 1396 PEN100 = 2.1274999107e+01_wp 1397 PEN200 = -2.5387384109e+01_wp 1398 PEN300 = 1.5469038167e+01_wp 1399 PEN400 = -3.3025876549_wp 1400 PEN010 = 6.6681505563_wp 1401 PEN110 = 2.2435057288_wp 1402 PEN210 = -2.5021299030_wp 1403 PEN310 = 3.2699521832e-01_wp 1404 PEN020 = -3.3540239802_wp 1405 PEN120 = -1.7531540640_wp 1406 PEN220 = 9.3976864981e-01_wp 1407 PEN030 = 1.2324834767_wp 1408 PEN130 = 2.7538550639e-01_wp 1409 PEN040 = -2.7963967985e-01_wp 1410 PEN001 = -1.3773949450_wp 1411 PEN101 = 3.3018402659_wp 1412 PEN201 = -1.6679755496_wp 1413 PEN011 = -1.3709540999_wp 1414 PEN111 = 1.4207577012e-01_wp 1415 PEN021 = 8.2799886843e-01_wp 1416 PEN002 = 1.7507069098e-02_wp 1417 PEN102 = 1.3880727538e-02_wp 1418 PEN012 = -2.8477365341e-01_wp 1419 ! 1420 APE000 = -1.6670376391e-01_wp 1421 APE100 = -5.6087643219e-02_wp 1422 APE200 = 6.2553247576e-02_wp 1423 APE300 = -8.1748804580e-03_wp 1424 APE010 = 1.6770119901e-01_wp 1425 APE110 = 8.7657703198e-02_wp 1426 APE210 = -4.6988432490e-02_wp 1427 APE020 = -9.2436260751e-02_wp 1428 APE120 = -2.0653912979e-02_wp 1429 APE030 = 2.7963967985e-02_wp 1430 APE001 = 3.4273852498e-02_wp 1431 APE101 = -3.5518942529e-03_wp 1432 APE011 = -4.1399943421e-02_wp 1433 APE002 = 7.1193413354e-03_wp 1434 ! 1435 BPE000 = 2.6468936504e-01_wp 1436 BPE100 = -6.3170583896e-01_wp 1437 BPE200 = 5.7736640125e-01_wp 1438 BPE300 = -1.6435438140e-01_wp 1439 BPE010 = 2.7912203607e-02_wp 1440 BPE110 = -6.2259666565e-02_wp 1441 BPE210 = 1.2204769966e-02_wp 1442 BPE020 = -2.1811574876e-02_wp 1443 BPE120 = 2.3383950895e-02_wp 1444 BPE030 = 3.4261630030e-03_wp 1445 BPE001 = 4.1079296834e-02_wp 1446 BPE101 = -4.1503681096e-02_wp 1447 BPE011 = 1.7676120780e-03_wp 1448 BPE002 = 1.7269476440e-04_wp 1449 ! 1450 CASE( 0 ) !== polynomial EOS-80 formulation ==! 1451 ! 745 1452 IF(lwp) WRITE(numout,*) 746 IF(lwp) WRITE(numout,*) ' use of linear eos rho(T) = rau0 * ( 1.0285 - rn_alpha * T )' 747 IF( lk_zdfddm ) CALL ctl_stop( ' double diffusive mixing parameterization requires', & 748 & ' that T and S are used as state variables' ) 749 ! 750 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 751 ralpbet = rn_alpha / rn_beta 752 IF(lwp) WRITE(numout,*) 753 IF(lwp) WRITE(numout,*) ' use of linear eos rho(T,S) = rau0 * ( rn_beta * S - rn_alpha * T )' 1453 IF(lwp) WRITE(numout,*) ' use of EOS-80 equation of state (pot. temp. and pract. salinity)' 1454 ! 1455 rdeltaS = 20._wp 1456 r1_S0 = 1._wp/40._wp 1457 r1_T0 = 1._wp/40._wp 1458 r1_Z0 = 1.e-4_wp 1459 ! 1460 EOS000 = 9.5356891948e+02_wp 1461 EOS100 = 1.7136499189e+02_wp 1462 EOS200 = -3.7501039454e+02_wp 1463 EOS300 = 5.1856810420e+02_wp 1464 EOS400 = -3.7264470465e+02_wp 1465 EOS500 = 1.4302533998e+02_wp 1466 EOS600 = -2.2856621162e+01_wp 1467 EOS010 = 1.0087518651e+01_wp 1468 EOS110 = -1.3647741861e+01_wp 1469 EOS210 = 8.8478359933_wp 1470 EOS310 = -7.2329388377_wp 1471 EOS410 = 1.4774410611_wp 1472 EOS510 = 2.0036720553e-01_wp 1473 EOS020 = -2.5579830599e+01_wp 1474 EOS120 = 2.4043512327e+01_wp 1475 EOS220 = -1.6807503990e+01_wp 1476 EOS320 = 8.3811577084_wp 1477 EOS420 = -1.9771060192_wp 1478 EOS030 = 1.6846451198e+01_wp 1479 EOS130 = -2.1482926901e+01_wp 1480 EOS230 = 1.0108954054e+01_wp 1481 EOS330 = -6.2675951440e-01_wp 1482 EOS040 = -8.0812310102_wp 1483 EOS140 = 1.0102374985e+01_wp 1484 EOS240 = -4.8340368631_wp 1485 EOS050 = 1.2079167803_wp 1486 EOS150 = 1.1515380987e-01_wp 1487 EOS060 = -2.4520288837e-01_wp 1488 EOS001 = 1.0748601068e+01_wp 1489 EOS101 = -1.7817043500e+01_wp 1490 EOS201 = 2.2181366768e+01_wp 1491 EOS301 = -1.6750916338e+01_wp 1492 EOS401 = 4.1202230403_wp 1493 EOS011 = -1.5852644587e+01_wp 1494 EOS111 = -7.6639383522e-01_wp 1495 EOS211 = 4.1144627302_wp 1496 EOS311 = -6.6955877448e-01_wp 1497 EOS021 = 9.9994861860_wp 1498 EOS121 = -1.9467067787e-01_wp 1499 EOS221 = -1.2177554330_wp 1500 EOS031 = -3.4866102017_wp 1501 EOS131 = 2.2229155620e-01_wp 1502 EOS041 = 5.9503008642e-01_wp 1503 EOS002 = 1.0375676547_wp 1504 EOS102 = -3.4249470629_wp 1505 EOS202 = 2.0542026429_wp 1506 EOS012 = 2.1836324814_wp 1507 EOS112 = -3.4453674320e-01_wp 1508 EOS022 = -1.2548163097_wp 1509 EOS003 = 1.8729078427e-02_wp 1510 EOS103 = -5.7238495240e-02_wp 1511 EOS013 = 3.8306136687e-01_wp 1512 ! 1513 ALP000 = -2.5218796628e-01_wp 1514 ALP100 = 3.4119354654e-01_wp 1515 ALP200 = -2.2119589983e-01_wp 1516 ALP300 = 1.8082347094e-01_wp 1517 ALP400 = -3.6936026529e-02_wp 1518 ALP500 = -5.0091801383e-03_wp 1519 ALP010 = 1.2789915300_wp 1520 ALP110 = -1.2021756164_wp 1521 ALP210 = 8.4037519952e-01_wp 1522 ALP310 = -4.1905788542e-01_wp 1523 ALP410 = 9.8855300959e-02_wp 1524 ALP020 = -1.2634838399_wp 1525 ALP120 = 1.6112195176_wp 1526 ALP220 = -7.5817155402e-01_wp 1527 ALP320 = 4.7006963580e-02_wp 1528 ALP030 = 8.0812310102e-01_wp 1529 ALP130 = -1.0102374985_wp 1530 ALP230 = 4.8340368631e-01_wp 1531 ALP040 = -1.5098959754e-01_wp 1532 ALP140 = -1.4394226233e-02_wp 1533 ALP050 = 3.6780433255e-02_wp 1534 ALP001 = 3.9631611467e-01_wp 1535 ALP101 = 1.9159845880e-02_wp 1536 ALP201 = -1.0286156825e-01_wp 1537 ALP301 = 1.6738969362e-02_wp 1538 ALP011 = -4.9997430930e-01_wp 1539 ALP111 = 9.7335338937e-03_wp 1540 ALP211 = 6.0887771651e-02_wp 1541 ALP021 = 2.6149576513e-01_wp 1542 ALP121 = -1.6671866715e-02_wp 1543 ALP031 = -5.9503008642e-02_wp 1544 ALP002 = -5.4590812035e-02_wp 1545 ALP102 = 8.6134185799e-03_wp 1546 ALP012 = 6.2740815484e-02_wp 1547 ALP003 = -9.5765341718e-03_wp 1548 ! 1549 BET000 = 2.1420623987_wp 1550 BET100 = -9.3752598635_wp 1551 BET200 = 1.9446303907e+01_wp 1552 BET300 = -1.8632235232e+01_wp 1553 BET400 = 8.9390837485_wp 1554 BET500 = -1.7142465871_wp 1555 BET010 = -1.7059677327e-01_wp 1556 BET110 = 2.2119589983e-01_wp 1557 BET210 = -2.7123520642e-01_wp 1558 BET310 = 7.3872053057e-02_wp 1559 BET410 = 1.2522950346e-02_wp 1560 BET020 = 3.0054390409e-01_wp 1561 BET120 = -4.2018759976e-01_wp 1562 BET220 = 3.1429341406e-01_wp 1563 BET320 = -9.8855300959e-02_wp 1564 BET030 = -2.6853658626e-01_wp 1565 BET130 = 2.5272385134e-01_wp 1566 BET230 = -2.3503481790e-02_wp 1567 BET040 = 1.2627968731e-01_wp 1568 BET140 = -1.2085092158e-01_wp 1569 BET050 = 1.4394226233e-03_wp 1570 BET001 = -2.2271304375e-01_wp 1571 BET101 = 5.5453416919e-01_wp 1572 BET201 = -6.2815936268e-01_wp 1573 BET301 = 2.0601115202e-01_wp 1574 BET011 = -9.5799229402e-03_wp 1575 BET111 = 1.0286156825e-01_wp 1576 BET211 = -2.5108454043e-02_wp 1577 BET021 = -2.4333834734e-03_wp 1578 BET121 = -3.0443885826e-02_wp 1579 BET031 = 2.7786444526e-03_wp 1580 BET002 = -4.2811838287e-02_wp 1581 BET102 = 5.1355066072e-02_wp 1582 BET012 = -4.3067092900e-03_wp 1583 BET003 = -7.1548119050e-04_wp 1584 ! 1585 PEN000 = -5.3743005340_wp 1586 PEN100 = 8.9085217499_wp 1587 PEN200 = -1.1090683384e+01_wp 1588 PEN300 = 8.3754581690_wp 1589 PEN400 = -2.0601115202_wp 1590 PEN010 = 7.9263222935_wp 1591 PEN110 = 3.8319691761e-01_wp 1592 PEN210 = -2.0572313651_wp 1593 PEN310 = 3.3477938724e-01_wp 1594 PEN020 = -4.9997430930_wp 1595 PEN120 = 9.7335338937e-02_wp 1596 PEN220 = 6.0887771651e-01_wp 1597 PEN030 = 1.7433051009_wp 1598 PEN130 = -1.1114577810e-01_wp 1599 PEN040 = -2.9751504321e-01_wp 1600 PEN001 = -6.9171176978e-01_wp 1601 PEN101 = 2.2832980419_wp 1602 PEN201 = -1.3694684286_wp 1603 PEN011 = -1.4557549876_wp 1604 PEN111 = 2.2969116213e-01_wp 1605 PEN021 = 8.3654420645e-01_wp 1606 PEN002 = -1.4046808820e-02_wp 1607 PEN102 = 4.2928871430e-02_wp 1608 PEN012 = -2.8729602515e-01_wp 1609 ! 1610 APE000 = -1.9815805734e-01_wp 1611 APE100 = -9.5799229402e-03_wp 1612 APE200 = 5.1430784127e-02_wp 1613 APE300 = -8.3694846809e-03_wp 1614 APE010 = 2.4998715465e-01_wp 1615 APE110 = -4.8667669469e-03_wp 1616 APE210 = -3.0443885826e-02_wp 1617 APE020 = -1.3074788257e-01_wp 1618 APE120 = 8.3359333577e-03_wp 1619 APE030 = 2.9751504321e-02_wp 1620 APE001 = 3.6393874690e-02_wp 1621 APE101 = -5.7422790533e-03_wp 1622 APE011 = -4.1827210323e-02_wp 1623 APE002 = 7.1824006288e-03_wp 1624 ! 1625 BPE000 = 1.1135652187e-01_wp 1626 BPE100 = -2.7726708459e-01_wp 1627 BPE200 = 3.1407968134e-01_wp 1628 BPE300 = -1.0300557601e-01_wp 1629 BPE010 = 4.7899614701e-03_wp 1630 BPE110 = -5.1430784127e-02_wp 1631 BPE210 = 1.2554227021e-02_wp 1632 BPE020 = 1.2166917367e-03_wp 1633 BPE120 = 1.5221942913e-02_wp 1634 BPE030 = -1.3893222263e-03_wp 1635 BPE001 = 2.8541225524e-02_wp 1636 BPE101 = -3.4236710714e-02_wp 1637 BPE011 = 2.8711395266e-03_wp 1638 BPE002 = 5.3661089288e-04_wp 1639 ! 1640 CASE( 1 ) !== Simplified EOS ==! 1641 IF(lwp) THEN 1642 WRITE(numout,*) 1643 WRITE(numout,*) ' use of simplified eos: rhd(dT=T-10,dS=S-35,Z) = ' 1644 WRITE(numout,*) ' [-a0*(1+lambda1/2*dT+mu1*Z)*dT + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS]/rau0' 1645 WRITE(numout,*) 1646 WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0 1647 WRITE(numout,*) ' saline cont. coef. rn_b0 = ', rn_b0 1648 WRITE(numout,*) ' cabbeling coef. rn_lambda1 = ', rn_lambda1 1649 WRITE(numout,*) ' cabbeling coef. rn_lambda2 = ', rn_lambda2 1650 WRITE(numout,*) ' thermobar. coef. rn_mu1 = ', rn_mu1 1651 WRITE(numout,*) ' thermobar. coef. rn_mu2 = ', rn_mu2 1652 WRITE(numout,*) ' 2nd cabbel. coef. rn_nu = ', rn_nu 1653 WRITE(numout,*) ' Caution: rn_beta0=0 incompatible with ddm parameterization ' 1654 ENDIF 754 1655 ! 755 1656 CASE DEFAULT !== ERROR in nn_eos ==! … … 759 1660 END SELECT 760 1661 ! 1662 rau0_rcp = rau0 * rcp 1663 r1_rau0 = 1._wp / rau0 1664 r1_rcp = 1._wp / rcp 1665 r1_rau0_rcp = 1._wp / rau0_rcp 1666 ! 1667 IF(lwp) WRITE(numout,*) 1668 IF(lwp) WRITE(numout,*) ' volumic mass of reference rau0 = ', rau0 , ' kg/m^3' 1669 IF(lwp) WRITE(numout,*) ' 1. / rau0 r1_rau0 = ', r1_rau0, ' m^3/kg' 1670 IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin' 1671 IF(lwp) WRITE(numout,*) ' rau0 * rcp rau0_rcp = ', rau0_rcp 1672 IF(lwp) WRITE(numout,*) ' 1. / ( rau0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp 1673 ! 761 1674 END SUBROUTINE eos_init 762 1675
Note: See TracChangeset
for help on using the changeset viewer.
