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