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