- Timestamp:
- 2015-07-16T11:04:29+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2_crs.F90
r5105 r5601 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 USE dom_oce! ocean space and time domain39 USE crs ! ocean space and time domain 34 40 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 41 ! 42 !USE in_out_manager ! I/O manager 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 crslbclnk ! ocean lateral boundary conditions 40 48 USE timing ! Timing 41 USE crs42 49 43 50 IMPLICIT NONE … … 46 53 ! !! * Interface 47 54 INTERFACE eos_crs 48 MODULE PROCEDURE eos_insitu_ crs, eos_insitu_pot_crs, eos_insitu_2d_crs55 MODULE PROCEDURE eos_insitu_pot , eos_insitu_2d 49 56 END INTERFACE 50 INTERFACE bn2_crs 51 MODULE PROCEDURE eos_bn2_crs 57 ! 58 INTERFACE eos_rab_crs 59 MODULE PROCEDURE rab_crs_3d, rab_crs_2d, rab_crs_0d 52 60 END INTERFACE 53 61 ! 54 62 PUBLIC eos_crs ! called by step, istate, tranpc and zpsgrd modules 63 PUBLIC bn2_crs ! called by step module 64 PUBLIC eos_rab_crs ! called by ldfslp, zdfddm, trabbl 55 65 PUBLIC eos_init_crs ! called by istate module 56 PUBLIC bn2_crs ! called by step module57 PUBLIC eos_alpbet_crs ! called by ldfslp module58 PUBLIC tfreez_crs ! called by sbcice_... modules59 66 60 67 ! !!* Namelist (nameos) * 61 68 INTEGER , PUBLIC :: nn_eos = 0 !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 62 REAL(wp), PUBLIC :: rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state) 63 REAL(wp), PUBLIC :: rn_beta = 7.7e-4_wp !: saline expension coeff. (linear equation of state) 64 65 REAL(wp), PUBLIC :: ralpbet_crs !: alpha / beta ratio 69 LOGICAL , PUBLIC :: ln_useCT = .FALSE. ! determine if eos_pt_from_ct is used to compute sst_m 70 71 ! !!! simplified eos coefficients 72 ! default value: Vallis 2006 73 REAL(wp) :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff. 74 REAL(wp) :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff. 75 REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2 76 REAL(wp) :: rn_lambda2 = 5.4914e-4_wp ! cabbeling coeff. in S^2 77 REAL(wp) :: rn_mu1 = 1.4970e-4_wp ! thermobaric coeff. in T 78 REAL(wp) :: rn_mu2 = 1.1090e-5_wp ! thermobaric coeff. in S 79 REAL(wp) :: rn_nu = 2.4341e-3_wp ! cabbeling coeff. in theta*salt 80 81 ! TEOS10/EOS80 parameters 82 REAL(wp) :: r1_S0, r1_T0, r1_Z0, rdeltaS 83 84 ! EOS parameters 85 REAL(wp) :: EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600 86 REAL(wp) :: EOS010 , EOS110 , EOS210 , EOS310 , EOS410 , EOS510 87 REAL(wp) :: EOS020 , EOS120 , EOS220 , EOS320 , EOS420 88 REAL(wp) :: EOS030 , EOS130 , EOS230 , EOS330 89 REAL(wp) :: EOS040 , EOS140 , EOS240 90 REAL(wp) :: EOS050 , EOS150 91 REAL(wp) :: EOS060 92 REAL(wp) :: EOS001 , EOS101 , EOS201 , EOS301 , EOS401 93 REAL(wp) :: EOS011 , EOS111 , EOS211 , EOS311 94 REAL(wp) :: EOS021 , EOS121 , EOS221 95 REAL(wp) :: EOS031 , EOS131 96 REAL(wp) :: EOS041 97 REAL(wp) :: EOS002 , EOS102 , EOS202 98 REAL(wp) :: EOS012 , EOS112 99 REAL(wp) :: EOS022 100 REAL(wp) :: EOS003 , EOS103 101 REAL(wp) :: EOS013 102 103 ! ALPHA parameters 104 REAL(wp) :: ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500 105 REAL(wp) :: ALP010 , ALP110 , ALP210 , ALP310 , ALP410 106 REAL(wp) :: ALP020 , ALP120 , ALP220 , ALP320 107 REAL(wp) :: ALP030 , ALP130 , ALP230 108 REAL(wp) :: ALP040 , ALP140 109 REAL(wp) :: ALP050 110 REAL(wp) :: ALP001 , ALP101 , ALP201 , ALP301 111 REAL(wp) :: ALP011 , ALP111 , ALP211 112 REAL(wp) :: ALP021 , ALP121 113 REAL(wp) :: ALP031 114 REAL(wp) :: ALP002 , ALP102 115 REAL(wp) :: ALP012 116 REAL(wp) :: ALP003 117 118 ! BETA parameters 119 REAL(wp) :: BET000 , BET100 , BET200 , BET300 , BET400 , BET500 120 REAL(wp) :: BET010 , BET110 , BET210 , BET310 , BET410 121 REAL(wp) :: BET020 , BET120 , BET220 , BET320 122 REAL(wp) :: BET030 , BET130 , BET230 123 REAL(wp) :: BET040 , BET140 124 REAL(wp) :: BET050 125 REAL(wp) :: BET001 , BET101 , BET201 , BET301 126 REAL(wp) :: BET011 , BET111 , BET211 127 REAL(wp) :: BET021 , BET121 128 REAL(wp) :: BET031 129 REAL(wp) :: BET002 , BET102 130 REAL(wp) :: BET012 131 REAL(wp) :: BET003 132 133 ! PEN parameters 134 REAL(wp) :: PEN000 , PEN100 , PEN200 , PEN300 , PEN400 135 REAL(wp) :: PEN010 , PEN110 , PEN210 , PEN310 136 REAL(wp) :: PEN020 , PEN120 , PEN220 137 REAL(wp) :: PEN030 , PEN130 138 REAL(wp) :: PEN040 139 REAL(wp) :: PEN001 , PEN101 , PEN201 140 REAL(wp) :: PEN011 , PEN111 141 REAL(wp) :: PEN021 142 REAL(wp) :: PEN002 , PEN102 143 REAL(wp) :: PEN012 144 145 ! ALPHA_PEN parameters 146 REAL(wp) :: APE000 , APE100 , APE200 , APE300 147 REAL(wp) :: APE010 , APE110 , APE210 148 REAL(wp) :: APE020 , APE120 149 REAL(wp) :: APE030 150 REAL(wp) :: APE001 , APE101 151 REAL(wp) :: APE011 152 REAL(wp) :: APE002 153 154 ! BETA_PEN parameters 155 REAL(wp) :: BPE000 , BPE100 , BPE200 , BPE300 156 REAL(wp) :: BPE010 , BPE110 , BPE210 157 REAL(wp) :: BPE020 , BPE120 158 REAL(wp) :: BPE030 159 REAL(wp) :: BPE001 , BPE101 160 REAL(wp) :: BPE011 161 REAL(wp) :: BPE002 66 162 67 163 !! * Substitutions … … 69 165 # include "vectopt_loop_substitute.h90" 70 166 !!---------------------------------------------------------------------- 71 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)72 !! $Id: eosbn2.F90 3294 2012-01-28 16:44:18Z rblod$167 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 168 !! $Id: eosbn2.F90 4990 2014-12-15 16:42:49Z timgraham $ 73 169 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 74 170 !!---------------------------------------------------------------------- 75 171 CONTAINS 76 172 77 SUBROUTINE eos_insitu_crs( pts, prd ) 78 !!---------------------------------------------------------------------- 79 !! *** ROUTINE eos_insitu *** 80 !! 81 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 82 !! potential temperature and salinity using an equation of state 83 !! defined through the namelist parameter nn_eos. 84 !! 85 !! ** Method : 3 cases: 86 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. 87 !! the in situ density is computed directly as a function of 88 !! potential temperature relative to the surface (the opa t 89 !! variable), salt and pressure (assuming no pressure variation 90 !! along geopotential surfaces, i.e. the pressure p in decibars 91 !! is approximated by the depth in meters. 92 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 93 !! with pressure p decibars 94 !! potential temperature t deg celsius 95 !! salinity s psu 96 !! reference volumic mass rau0 kg/m**3 97 !! in situ volumic mass rho kg/m**3 98 !! in situ density anomalie prd no units 99 !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, 100 !! t = 40 deg celcius, s=40 psu 101 !! nn_eos = 1 : linear equation of state function of temperature only 102 !! prd(t) = 0.0285 - rn_alpha * t 103 !! nn_eos = 2 : linear equation of state function of temperature and 104 !! salinity 105 !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 106 !! Note that no boundary condition problem occurs in this routine 107 !! as pts are defined over the whole domain. 108 !! 109 !! ** Action : compute prd , the in situ density (no units) 110 !! 111 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 112 !!---------------------------------------------------------------------- 113 !! 114 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 115 ! ! 2 : salinity [psu] 116 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-] 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) :: zrau0r ! - - 125 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 126 !!---------------------------------------------------------------------- 127 128 ! 129 IF( nn_timing == 1 ) CALL timing_start('eos') 130 ! 131 CALL wrk_alloc( jpi, jpj, jpk, zws ) 132 ! 133 SELECT CASE( nn_eos ) 134 ! 135 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 136 zrau0r = 1.e0 / rau0 137 !CDIR NOVERRCHK 138 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 139 ! 140 DO jk = 1, jpkm1 141 DO jj = 1, jpj 142 DO ji = 1, jpi 143 zt = pts (ji,jj,jk,jp_tem) 144 zs = pts (ji,jj,jk,jp_sal) 145 zh = gdept_crs(ji,jj,jk) ! depth 146 zsr= zws (ji,jj,jk) ! square root salinity 147 ! 148 ! compute volumic mass pure water at atm pressure 149 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 150 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 151 ! seawater volumic mass atm pressure 152 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 153 & -4.0899e-3_wp ) *zt+0.824493_wp 154 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 155 zr4= 4.8314e-4_wp 156 ! 157 ! potential volumic mass (reference to the surface) 158 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 159 ! 160 ! add the compression terms 161 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 162 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 163 zb = zbw + ze * zs 164 ! 165 zd = -2.042967e-2_wp 166 zc = (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 167 zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 168 za = ( zd*zsr + zc ) *zs + zaw 169 ! 170 zb1= (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp 171 za1= ( ( 2.326469e-3_wp*zt+1.553190_wp) *zt-65.00517_wp ) *zt+1044.077_wp 172 zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 173 zk0= ( zb1*zsr + za1 )*zs + zkw 174 ! 175 ! masked in situ density anomaly 176 prd(ji,jj,jk) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) & 177 & - rau0 ) * zrau0r * tmask_crs(ji,jj,jk) 178 END DO 179 END DO 180 END DO 181 ! 182 CASE( 1 ) !== Linear formulation function of temperature only ==! 183 DO jk = 1, jpkm1 184 prd(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask_crs(:,:,jk) 185 END DO 186 ! 187 CASE( 2 ) !== Linear formulation function of temperature and salinity ==! 188 DO jk = 1, jpkm1 189 prd(:,:,jk) = ( rn_beta * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask_crs(:,:,jk) 190 END DO 191 ! 192 END SELECT 193 ! 194 ! IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos : ', ovlap=1, kdim=jpk ) 195 ! 196 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 197 ! 198 IF( nn_timing == 1 ) CALL timing_stop('eos') 199 ! 200 END SUBROUTINE eos_insitu_crs 201 202 203 SUBROUTINE eos_insitu_pot_crs( pts, prd, prhop ) 173 SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 204 174 !!---------------------------------------------------------------------- 205 175 !! *** ROUTINE eos_insitu_pot *** … … 210 180 !! namelist parameter nn_eos. 211 181 !! 212 !! ** Method :213 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state.214 !! the in situ density is computed directly as a function of215 !! potential temperature relative to the surface (the opa t216 !! variable), salt and pressure (assuming no pressure variation217 !! along geopotential surfaces, i.e. the pressure p in decibars218 !! is approximated by the depth in meters.219 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0220 !! rhop(t,s) = rho(t,s,0)221 !! with pressure p decibars222 !! potential temperature t deg celsius223 !! salinity s psu224 !! reference volumic mass rau0 kg/m**3225 !! in situ volumic mass rho kg/m**3226 !! in situ density anomalie prd no units227 !!228 !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,229 !! t = 40 deg celcius, s=40 psu230 !!231 !! nn_eos = 1 : linear equation of state function of temperature only232 !! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t233 !! rhop(t,s) = rho(t,s)234 !!235 !! nn_eos = 2 : linear equation of state function of temperature and236 !! salinity237 !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0238 !! = rn_beta * s - rn_alpha * tn - 1.239 !! rhop(t,s) = rho(t,s)240 !! Note that no boundary condition problem occurs in this routine241 !! as (tn,sn) or (ta,sa) are defined over the whole domain.242 !!243 182 !! ** Action : - prd , the in situ density (no units) 244 183 !! - prhop, the potential volumic mass (Kg/m3) 245 184 !! 246 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 247 !! Brown and Campana, Mon. Weather Rev., 1978 248 !!---------------------------------------------------------------------- 249 !! 250 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 185 !!---------------------------------------------------------------------- 186 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 251 187 ! ! 2 : salinity [psu] 252 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-] 253 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prhop ! potential density (surface referenced) 254 ! 255 INTEGER :: ji, jj, jk ! dummy loop indices 256 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! local scalars 257 REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zrau0r ! - - 258 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 259 !!---------------------------------------------------------------------- 260 ! 261 IF( nn_timing == 1 ) CALL timing_start('eos-p') 262 ! 263 CALL wrk_alloc( jpi, jpj, jpk, zws ) 188 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk ), INTENT( out) :: prd ! in situ density [-] 189 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk ), INTENT( out) :: prhop ! potential density (surface referenced) 190 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk ), INTENT(in ) :: pdep ! depth [m] 191 ! 192 INTEGER :: ji, jj, jk ! dummy loop indices 193 REAL(wp) :: zt , zh , zs , ztm ! local scalars 194 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 195 !!---------------------------------------------------------------------- 196 ! 197 IF( nn_timing == 1 ) CALL timing_start('eos-pot_crs') 264 198 ! 265 199 SELECT CASE ( nn_eos ) 266 200 ! 267 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 268 zrau0r = 1.e0 / rau0 269 !CDIR NOVERRCHK 270 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 201 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 271 202 ! 272 203 DO jk = 1, jpkm1 273 DO jj = 1, jpj 274 DO ji = 1, jpi 275 zt = pts (ji,jj,jk,jp_tem) 276 zs = pts (ji,jj,jk,jp_sal) 277 zh = gdept_crs(ji,jj,jk) ! depth 278 zsr= zws (ji,jj,jk) ! square root salinity 279 ! 280 ! compute volumic mass pure water at atm pressure 281 zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt & 282 & -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 283 ! seawater volumic mass atm pressure 284 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 285 & -4.0899e-3_wp ) *zt+0.824493_wp 286 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 287 zr4= 4.8314e-4_wp 288 ! 289 ! potential volumic mass (reference to the surface) 290 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 291 ! 292 ! save potential volumic mass 293 prhop(ji,jj,jk) = zrhop * tmask_crs(ji,jj,jk) 294 ! 295 ! add the compression terms 296 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 297 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 298 zb = zbw + ze * zs 299 ! 300 zd = -2.042967e-2_wp 301 zc = (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 302 zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 303 za = ( zd*zsr + zc ) *zs + zaw 304 ! 305 zb1= ( -0.1909078_wp *zt+7.390729_wp ) *zt-55.87545_wp 306 za1= ( ( 2.326469e-3_wp*zt+1.553190_wp ) *zt-65.00517_wp ) *zt + 1044.077_wp 307 zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 308 zk0= ( zb1*zsr + za1 )*zs + zkw 309 ! 310 ! masked in situ density anomaly 311 prd(ji,jj,jk) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) & 312 & - rau0 ) * zrau0r * tmask_crs(ji,jj,jk) 204 DO jj = 1, jpj_crs 205 DO ji = 1, jpi_crs 206 ! 207 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 208 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 209 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 210 ztm = tmask_crs(ji,jj,jk) ! tmask 211 ! 212 zn3 = EOS013*zt & 213 & + EOS103*zs+EOS003 214 ! 215 zn2 = (EOS022*zt & 216 & + EOS112*zs+EOS012)*zt & 217 & + (EOS202*zs+EOS102)*zs+EOS002 218 ! 219 zn1 = (((EOS041*zt & 220 & + EOS131*zs+EOS031)*zt & 221 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 222 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 223 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 224 ! 225 zn0 = (((((EOS060*zt & 226 & + EOS150*zs+EOS050)*zt & 227 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 228 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 229 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 230 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 231 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 232 ! 233 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 234 ! 235 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 236 ! 237 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 313 238 END DO 314 239 END DO 315 240 END DO 316 241 ! 317 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 242 CASE( 1 ) !== simplified EOS ==! 243 ! 318 244 DO jk = 1, jpkm1 319 prd (:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask_crs(:,:,jk) 320 prhop(:,:,jk) = ( 1.e0_wp + prd (:,:,jk) ) * rau0 * tmask_crs(:,:,jk) 245 DO jj = 1, jpj_crs 246 DO ji = 1, jpi_crs 247 zt = pts (ji,jj,jk,jp_tem) - 10._wp 248 zs = pts (ji,jj,jk,jp_sal) - 35._wp 249 zh = pdep (ji,jj,jk) 250 ztm = tmask_crs(ji,jj,jk) 251 ! ! potential density referenced at the surface 252 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt & 253 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs & 254 & - rn_nu * zt * zs 255 prhop(ji,jj,jk) = ( rau0 + zn ) * ztm 256 ! ! density anomaly (masked) 257 zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 258 prd(ji,jj,jk) = zn * r1_rau0 * ztm 259 ! 260 END DO 261 END DO 321 262 END DO 322 263 ! 323 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==!324 DO jk = 1, jpkm1325 prd (:,:,jk) = ( rn_beta * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask_crs(:,:,jk)326 prhop(:,:,jk) = ( 1.e0_wp + prd (:,:,jk) ) * rau0 * tmask_crs(:,:,jk)327 END DO328 !329 264 END SELECT 330 265 ! 331 ! IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk ) 332 ! 333 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 334 ! 335 IF( nn_timing == 1 ) CALL timing_stop('eos-p') 336 ! 337 END SUBROUTINE eos_insitu_pot_crs 338 339 340 SUBROUTINE eos_insitu_2d_crs( pts, pdep, prd ) 266 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk ) 267 ! 268 IF( nn_timing == 1 ) CALL timing_stop('eos-pot_crs') 269 ! 270 END SUBROUTINE eos_insitu_pot 271 272 SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 341 273 !!---------------------------------------------------------------------- 342 274 !! *** ROUTINE eos_insitu_2d *** … … 346 278 !! defined through the namelist parameter nn_eos. * 2D field case 347 279 !! 348 !! ** Method : 349 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. 350 !! the in situ density is computed directly as a function of 351 !! potential temperature relative to the surface (the opa t 352 !! variable), salt and pressure (assuming no pressure variation 353 !! along geopotential surfaces, i.e. the pressure p in decibars 354 !! is approximated by the depth in meters. 355 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 356 !! with pressure p decibars 357 !! potential temperature t deg celsius 358 !! salinity s psu 359 !! reference volumic mass rau0 kg/m**3 360 !! in situ volumic mass rho kg/m**3 361 !! in situ density anomalie prd no units 362 !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, 363 !! t = 40 deg celcius, s=40 psu 364 !! nn_eos = 1 : linear equation of state function of temperature only 365 !! prd(t) = 0.0285 - rn_alpha * t 366 !! nn_eos = 2 : linear equation of state function of temperature and 367 !! salinity 368 !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 369 !! Note that no boundary condition problem occurs in this routine 370 !! as pts are defined over the whole domain. 371 !! 372 !! ** Action : - prd , the in situ density (no units) 373 !! 374 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 375 !!---------------------------------------------------------------------- 376 !! 377 REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 280 !! ** Action : - prd , the in situ density (no units) (unmasked) 281 !! 282 !!---------------------------------------------------------------------- 283 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 378 284 ! ! 2 : salinity [psu] 379 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 380 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 381 !! 382 INTEGER :: ji, jj ! dummy loop indices 383 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! temporary scalars 384 REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zmask ! - - 385 REAL(wp), POINTER, DIMENSION(:,:) :: zws 386 !!---------------------------------------------------------------------- 387 ! 388 !WRITE(numout,*) ' pts1 ' , pts(:,:,1) 389 !WRITE(numout,*) ' pts2 ' , pts(:,:,2) 390 !WRITE(numout,*) ' jpi ' , jpi 391 !WRITE(numout,*) ' fs_jpim1 ' , fs_jpim1 392 !WRITE(numout,*) ' dim ' , size(pts,1) 393 IF( nn_timing == 1 ) CALL timing_start('eos2d') 394 ! 395 CALL wrk_alloc( jpi, jpj, zws ) 396 ! 397 285 REAL(wp), DIMENSION(jpi_crs,jpj_crs) , INTENT(in ) :: pdep ! depth [m] 286 REAL(wp), DIMENSION(jpi_crs,jpj_crs) , INTENT( out) :: prd ! in situ density 287 ! 288 INTEGER :: ji, jj, jk ! dummy loop indices 289 REAL(wp) :: zt , zh , zs ! local scalars 290 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 291 !!---------------------------------------------------------------------- 292 ! 293 IF( nn_timing == 1 ) CALL timing_start('eos2d') 294 ! 398 295 prd(:,:) = 0._wp 399 296 ! 400 297 SELECT CASE( nn_eos ) 401 298 ! 402 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 403 ! 404 !CDIR NOVERRCHK 405 DO jj = 1, jpjm1 406 !CDIR NOVERRCHK 407 DO ji = 1, fs_jpim1 ! vector opt. 408 zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) ) 299 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 300 ! 301 DO jj = 1, jpj_crsm1 302 DO ji = 1, jpi_crsm1 ! vector opt. 303 ! 304 zh = pdep(ji,jj) * r1_Z0 ! depth 305 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 306 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 307 ! 308 zn3 = EOS013*zt & 309 & + EOS103*zs+EOS003 310 ! 311 zn2 = (EOS022*zt & 312 & + EOS112*zs+EOS012)*zt & 313 & + (EOS202*zs+EOS102)*zs+EOS002 314 ! 315 zn1 = (((EOS041*zt & 316 & + EOS131*zs+EOS031)*zt & 317 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 318 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 319 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 320 ! 321 zn0 = (((((EOS060*zt & 322 & + EOS150*zs+EOS050)*zt & 323 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 324 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 325 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 326 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 327 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 328 ! 329 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 330 ! 331 prd(ji,jj) = zn * r1_rau0 - 1._wp ! unmasked in situ density anomaly 332 ! 409 333 END DO 410 334 END DO 411 DO jj = 1, jpjm1 412 DO ji = 1, fs_jpim1 ! vector opt. 413 zmask = tmask_crs(ji,jj,1) ! land/sea bottom mask = surf. mask 414 zt = pts (ji,jj,jp_tem) ! interpolated T 415 zs = pts (ji,jj,jp_sal) ! interpolated S 416 zsr = zws (ji,jj) ! square root of interpolated S 417 zh = pdep (ji,jj) ! depth at the partial step level 418 ! 419 ! compute volumic mass pure water at atm pressure 420 zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt & 421 & -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 422 ! seawater volumic mass atm pressure 423 zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt & 424 & -4.0899e-3_wp ) *zt+0.824493_wp 425 zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 426 zr4 = 4.8314e-4_wp 427 ! 428 ! potential volumic mass (reference to the surface) 429 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 430 ! 431 ! add the compression terms 432 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 433 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 434 zb = zbw + ze * zs 435 ! 436 zd = -2.042967e-2_wp 437 zc = (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 438 zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt -4.721788_wp 439 za = ( zd*zsr + zc ) *zs + zaw 440 ! 441 zb1= (-0.1909078_wp *zt+7.390729_wp ) *zt-55.87545_wp 442 za1= ( ( 2.326469e-3_wp*zt+1.553190_wp ) *zt-65.00517_wp ) *zt+1044.077_wp 443 zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt & 444 & +2098.925_wp ) *zt+190925.6_wp 445 zk0= ( zb1*zsr + za1 )*zs + zkw 446 ! 447 ! masked in situ density anomaly 448 prd(ji,jj) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) - rau0 ) / rau0 * zmask 335 ! 336 CALL crs_lbc_lnk( prd, 'T', 1. ) ! Lateral boundary conditions 337 ! 338 CASE( 1 ) !== simplified EOS ==! 339 ! 340 DO jj = 1, jpj_crsm1 341 DO ji = 1, jpi_crsm1 ! vector opt. 342 ! 343 zt = pts (ji,jj,jp_tem) - 10._wp 344 zs = pts (ji,jj,jp_sal) - 35._wp 345 zh = pdep (ji,jj) ! depth at the partial step level 346 ! 347 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt & 348 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 349 & - rn_nu * zt * zs 350 ! 351 prd(ji,jj) = zn * r1_rau0 ! unmasked in situ density anomaly 352 ! 449 353 END DO 450 354 END DO 451 355 ! 452 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 453 DO jj = 1, jpjm1 454 DO ji = 1, fs_jpim1 ! vector opt. 455 prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask_crs(ji,jj,1) 456 END DO 457 END DO 458 ! 459 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 460 DO jj = 1, jpjm1 461 DO ji = 1, fs_jpim1 ! vector opt. 462 prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask_crs(ji,jj,1) 463 END DO 464 END DO 356 CALL crs_lbc_lnk( prd, 'T', 1. ) ! Lateral boundary conditions 465 357 ! 466 358 END SELECT 467 !WRITE(numout,*) ' prd ' , prd(:,:) 468 !WRITE(numout,*) ' zws ' , zws(:,:) 469 !WRITE(numout,*) ' pdep ' , pdep(:,:) 470 471 472 473 ! IF(ln_ctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 474 ! 475 CALL wrk_dealloc( jpi, jpj, zws ) 476 ! 477 IF( nn_timing == 1 ) CALL timing_stop('eos2d') 478 ! 479 END SUBROUTINE eos_insitu_2d_crs 480 481 482 SUBROUTINE eos_bn2_crs( pts, pn2 ) 483 !!---------------------------------------------------------------------- 484 !! *** ROUTINE eos_bn2 *** 485 !! 486 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the time- 487 !! step of the input arguments 488 !! 489 !! ** Method : 490 !! * nn_eos = 0 : UNESCO sea water properties 491 !! The brunt-vaisala frequency is computed using the polynomial 492 !! polynomial expression of McDougall (1987): 493 !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 494 !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 495 !! computed and used in zdfddm module : 496 !! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 497 !! * nn_eos = 1 : linear equation of state (temperature only) 498 !! N^2 = grav * rn_alpha * dk[ t ]/e3w 499 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 500 !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 501 !! The use of potential density to compute N^2 introduces e r r o r 502 !! in the sign of N^2 at great depths. We recommand the use of 503 !! nn_eos = 0, except for academical studies. 504 !! Macro-tasked on horizontal slab (jk-loop) 505 !! N.B. N^2 is set to zero at the first level (JK=1) in inidtr 506 !! and is never used at this level. 507 !! 508 !! ** Action : - pn2 : the brunt-vaisala frequency 509 !! 510 !! References : McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987. 511 !!---------------------------------------------------------------------- 512 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 513 ! ! 2 : salinity [psu] 514 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pn2 ! Brunt-Vaisala frequency [s-1] 515 !! 516 INTEGER :: ji, jj, jk ! dummy loop indices 517 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 518 #if defined key_zdfddm 519 REAL(wp) :: zds ! local scalars 520 #endif 521 !!---------------------------------------------------------------------- 522 523 ! 524 IF( nn_timing == 1 ) CALL timing_start('bn2') 525 ! 526 ! pn2 : interior points only (2=< jk =< jpkm1 ) 527 ! -------------------------- 528 ! 529 SELECT CASE( nn_eos ) 530 ! 531 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 532 DO jk = 2, jpkm1 533 DO jj = 1, jpj 534 DO ji = 1, jpi 535 zgde3w = grav / e3w_max_crs(ji,jj,jk) 536 zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-pt 537 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 538 zh = gdepw_crs(ji,jj,jk) ! depth in meters at w-point 539 ! 540 zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & ! ratio alpha/beta 541 & - 0.203814e-03_wp ) * zt & 542 & + 0.170907e-01_wp ) * zt & 543 & + 0.665157e-01_wp & 544 & + ( - 0.678662e-05_wp * zs & 545 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs & 546 & + ( ( - 0.302285e-13_wp * zh & 547 & - 0.251520e-11_wp * zs & 548 & + 0.512857e-12_wp * zt * zt ) * zh & 549 & - 0.164759e-06_wp * zs & 550 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt & 551 & + 0.380374e-04_wp ) * zh 359 ! 360 IF(ln_ctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 361 ! 362 IF( nn_timing == 1 ) CALL timing_stop('eos2d') 363 ! 364 END SUBROUTINE eos_insitu_2d 365 366 SUBROUTINE rab_crs_3d( pts, pab ) 367 !!---------------------------------------------------------------------- 368 !! *** ROUTINE rab_3d *** 369 !! 370 !! ** Purpose : Calculates thermal/haline expansion ratio at T-points 371 !! 372 !! ** Method : calculates alpha / beta at T-points 373 !! 374 !! ** Action : - pab : thermal/haline expansion ratio at T-points 375 !!---------------------------------------------------------------------- 376 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 377 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk,jpts), INTENT( out) :: pab ! thermal/haline expansion ratio 378 ! 379 INTEGER :: ji, jj, jk ! dummy loop indices 380 REAL(wp) :: zt , zh , zs , ztm ! local scalars 381 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 382 !!---------------------------------------------------------------------- 383 ! 384 IF( nn_timing == 1 ) CALL timing_start('rab_3d') 385 ! 386 SELECT CASE ( nn_eos ) 387 ! 388 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 389 ! 390 DO jk = 1, jpkm1 391 DO jj = 1, jpj_crs 392 DO ji = 1, jpi_crs 393 ! 394 zh = gdept_crs(ji,jj,jk) * r1_Z0 ! depth 395 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 396 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 397 ztm = tmask_crs(ji,jj,jk) ! tmask 398 ! 399 ! alpha 400 zn3 = ALP003 401 ! 402 zn2 = ALP012*zt + ALP102*zs+ALP002 403 ! 404 zn1 = ((ALP031*zt & 405 & + ALP121*zs+ALP021)*zt & 406 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 407 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 552 408 ! 553 zbeta = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & ! beta 554 & - 0.301985e-05_wp ) * zt & 555 & + 0.785567e-03_wp & 556 & + ( 0.515032e-08_wp * zs & 557 & + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs & 558 & + ( ( 0.121551e-17_wp * zh & 559 & - 0.602281e-15_wp * zs & 560 & - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh & 561 & + 0.408195e-10_wp * zs & 562 & + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt & 563 & - 0.121555e-07_wp ) * zh 409 zn0 = ((((ALP050*zt & 410 & + ALP140*zs+ALP040)*zt & 411 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 412 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 413 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 414 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 564 415 ! 565 !cbr zgde3w: divide by 0 566 !pn2(ji,jj,jk) = zgde3w * zbeta * tmask_crs(ji,jj,jk) & ! N^2 567 ! & * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 568 ! & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 569 pn2(ji,jj,jk) = zbeta * tmask_crs(ji,jj,jk) & ! N^2 570 & * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 571 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 572 IF( e3w_max_crs(ji,jj,jk) .NE. 0._wp ) pn2(ji,jj,jk) = zgde3w * e3w_max_crs(ji,jj,jk) 573 574 #if defined key_zdfddm 575 ! !!bug **** caution a traiter zds=dk[S]= 0 !!!! 576 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s]) 577 IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 578 rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 579 #endif 416 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 417 ! 418 pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm 419 ! 420 ! beta 421 zn3 = BET003 422 ! 423 zn2 = BET012*zt + BET102*zs+BET002 424 ! 425 zn1 = ((BET031*zt & 426 & + BET121*zs+BET021)*zt & 427 & + (BET211*zs+BET111)*zs+BET011)*zt & 428 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 429 ! 430 zn0 = ((((BET050*zt & 431 & + BET140*zs+BET040)*zt & 432 & + (BET230*zs+BET130)*zs+BET030)*zt & 433 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 434 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 435 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 436 ! 437 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 438 ! 439 pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm 440 ! 580 441 END DO 581 442 END DO 582 443 END DO 583 444 ! 584 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 585 DO jk = 2, jpkm1 586 pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / e3w_max_crs(:,:,jk) * tmask_crs(:,:,jk) 587 END DO 588 ! 589 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 590 DO jk = 2, jpkm1 591 !cbr: bug divide by 0. 592 !pn2(:,:,jk) = grav * ( rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) & 593 ! & - rn_beta * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) ) ) & 594 ! & / e3w_max_crs(:,:,jk) * tmask_crs(:,:,jk) 595 pn2(:,:,jk) = grav * ( rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) & 596 & - rn_beta * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) ) ) & 597 & * tmask_crs(:,:,jk) 598 DO jj = 1, jpj 599 DO ji = 1, jpi 600 IF( e3w_max_crs(ji,jj,jk) .NE. 0._wp ) pn2(ji,jj,jk) = pn2(ji,jj,jk) / e3w_max_crs(ji,jj,jk) 601 ENDDO 602 ENDDO 603 END DO 604 #if defined key_zdfddm 605 DO jk = 2, jpkm1 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 606 DO jj = 1, jpj 607 DO ji = 1, jpi 608 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 609 IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 610 rrau(ji,jj,jk) = ralpbet_crs * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 445 CASE( 1 ) !== simplified EOS ==! 446 ! 447 DO jk = 1, jpkm1 448 DO jj = 1, jpj_crs 449 DO ji = 1, jpi_crs 450 zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 451 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 452 zh = gdept_crs(ji,jj,jk) ! depth in meters at t-point 453 ztm = tmask_crs(ji,jj,jk) ! land/sea bottom mask = surf. mask 454 ! 455 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 456 pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm ! alpha 457 ! 458 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 459 pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm ! beta 460 ! 611 461 END DO 612 462 END DO 613 463 END DO 614 #endif615 END SELECT616 617 ! IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk )618 #if defined key_zdfddm619 ! IF(ln_ctl) CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk )620 #endif621 !622 IF( nn_timing == 1 ) CALL timing_stop('bn2')623 !624 END SUBROUTINE eos_bn2_crs625 626 627 SUBROUTINE eos_alpbet_crs( pts, palpbet, beta0 )628 !!----------------------------------------------------------------------629 !! *** ROUTINE eos_alpbet ***630 !!631 !! ** Purpose : Calculates the in situ thermal/haline expansion ratio at T-points632 !!633 !! ** Method : calculates alpha / beta ratio at T-points634 !! * nn_eos = 0 : UNESCO sea water properties635 !! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial636 !! polynomial expression of McDougall (1987).637 !! Scalar beta0 is returned = 1.638 !! * nn_eos = 1 : linear equation of state (temperature only)639 !! The ratio is undefined, so we return alpha as palpbet640 !! Scalar beta0 is returned = 0.641 !! * nn_eos = 2 : linear equation of state (temperature & salinity)642 !! The alpha/beta ratio is returned as ralpbet643 !! Scalar beta0 is returned = 1.644 !!645 !! ** Action : - palpbet : thermal/haline expansion ratio at T-points646 !! : beta0 : 1. or 0.647 !!----------------------------------------------------------------------648 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity649 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpbet ! thermal/haline expansion ratio650 REAL(wp), INTENT( out) :: beta0 ! set = 1 except with case 1 eos, rho=rho(T)651 !!652 INTEGER :: ji, jj, jk ! dummy loop indices653 REAL(wp) :: zt, zs, zh ! local scalars654 !!----------------------------------------------------------------------655 !656 IF( nn_timing == 1 ) CALL timing_start('eos_alpbet')657 !658 SELECT CASE ( nn_eos )659 !660 CASE ( 0 ) ! Jackett and McDougall (1994) formulation661 DO jk = 1, jpk662 DO jj = 1, jpj663 DO ji = 1, jpi664 zt = pts(ji,jj,jk,jp_tem) ! potential temperature665 zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35)666 zh = fsdept(ji,jj,jk) ! depth in meters667 !668 palpbet(ji,jj,jk) = &669 & ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt &670 & - 0.203814e-03_wp ) * zt &671 & + 0.170907e-01_wp ) * zt &672 & + 0.665157e-01_wp &673 & + ( - 0.678662e-05_wp * zs &674 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs &675 & + ( ( - 0.302285e-13_wp * zh &676 & - 0.251520e-11_wp * zs &677 & + 0.512857e-12_wp * zt * zt ) * zh &678 & - 0.164759e-06_wp * zs &679 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt &680 & + 0.380374e-04_wp ) * zh681 END DO682 END DO683 END DO684 beta0 = 1._wp685 !686 CASE ( 1 ) !== Linear formulation = F( temperature ) ==!687 palpbet(:,:,:) = rn_alpha688 beta0 = 0._wp689 !690 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!691 palpbet(:,:,:) = ralpbet_crs692 beta0 = 1._wp693 464 ! 694 465 CASE DEFAULT … … 699 470 END SELECT 700 471 ! 701 IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet') 702 ! 703 END SUBROUTINE eos_alpbet_crs 704 705 706 FUNCTION tfreez_crs( psal ) RESULT( ptf ) 472 IF(ln_ctl) CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', & 473 & tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', ovlap=1, kdim=jpk ) 474 ! 475 IF( nn_timing == 1 ) CALL timing_stop('rab_3d') 476 ! 477 END SUBROUTINE rab_crs_3d 478 479 SUBROUTINE rab_crs_2d( pts, pdep, pab ) 480 !!---------------------------------------------------------------------- 481 !! *** ROUTINE rab_2d *** 482 !! 483 !! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked) 484 !! 485 !! ** Action : - pab : thermal/haline expansion ratio at T-points 486 !!---------------------------------------------------------------------- 487 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpts) , INTENT(in ) :: pts ! pot. temperature & salinity 488 REAL(wp), DIMENSION(jpi_crs,jpj_crs) , INTENT(in ) :: pdep ! depth [m] 489 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpts) , INTENT( out) :: pab ! thermal/haline expansion ratio 490 ! 491 INTEGER :: ji, jj, jk ! dummy loop indices 492 REAL(wp) :: zt , zh , zs ! local scalars 493 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 494 !!---------------------------------------------------------------------- 495 ! 496 IF( nn_timing == 1 ) CALL timing_start('rab_2d') 497 ! 498 pab(:,:,:) = 0._wp 499 ! 500 SELECT CASE ( nn_eos ) 501 ! 502 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 503 ! 504 DO jj = 1, jpj_crsm1 505 DO ji = 1, jpi_crsm1 ! vector opt. 506 ! 507 zh = pdep(ji,jj) * r1_Z0 ! depth 508 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 509 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 510 ! 511 ! alpha 512 zn3 = ALP003 513 ! 514 zn2 = ALP012*zt + ALP102*zs+ALP002 515 ! 516 zn1 = ((ALP031*zt & 517 & + ALP121*zs+ALP021)*zt & 518 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 519 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 520 ! 521 zn0 = ((((ALP050*zt & 522 & + ALP140*zs+ALP040)*zt & 523 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 524 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 525 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 526 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 527 ! 528 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 529 ! 530 pab(ji,jj,jp_tem) = zn * r1_rau0 531 ! 532 ! beta 533 zn3 = BET003 534 ! 535 zn2 = BET012*zt + BET102*zs+BET002 536 ! 537 zn1 = ((BET031*zt & 538 & + BET121*zs+BET021)*zt & 539 & + (BET211*zs+BET111)*zs+BET011)*zt & 540 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 541 ! 542 zn0 = ((((BET050*zt & 543 & + BET140*zs+BET040)*zt & 544 & + (BET230*zs+BET130)*zs+BET030)*zt & 545 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 546 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 547 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 548 ! 549 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 550 ! 551 pab(ji,jj,jp_sal) = zn / zs * r1_rau0 552 ! 553 ! 554 END DO 555 END DO 556 ! 557 CALL crs_lbc_lnk( pab(:,:,jp_tem), 'T', 1. ) ! Lateral boundary conditions 558 CALL crs_lbc_lnk( pab(:,:,jp_sal), 'T', 1. ) 559 ! 560 CASE( 1 ) !== simplified EOS ==! 561 ! 562 DO jj = 1, jpj_crsm1 563 DO ji = 1, jpi_crsm1 ! vector opt. 564 ! 565 zt = pts (ji,jj,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 566 zs = pts (ji,jj,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 567 zh = pdep (ji,jj) ! depth at the partial step level 568 ! 569 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 570 pab(ji,jj,jp_tem) = zn * r1_rau0 ! alpha 571 ! 572 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 573 pab(ji,jj,jp_sal) = zn * r1_rau0 ! beta 574 ! 575 END DO 576 END DO 577 ! 578 CALL crs_lbc_lnk( pab(:,:,jp_tem), 'T', 1. ) ! Lateral boundary conditions 579 CALL crs_lbc_lnk( pab(:,:,jp_sal), 'T', 1. ) 580 ! 581 CASE DEFAULT 582 IF(lwp) WRITE(numout,cform_err) 583 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 584 nstop = nstop + 1 585 ! 586 END SELECT 587 ! 588 IF(ln_ctl) CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', & 589 & tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' ) 590 ! 591 IF( nn_timing == 1 ) CALL timing_stop('rab_2d') 592 ! 593 END SUBROUTINE rab_crs_2d 594 595 596 SUBROUTINE rab_crs_0d( pts, pdep, pab ) 597 !!---------------------------------------------------------------------- 598 !! *** ROUTINE rab_0d *** 599 !! 600 !! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked) 601 !! 602 !! ** Action : - pab : thermal/haline expansion ratio at T-points 603 !!---------------------------------------------------------------------- 604 REAL(wp), DIMENSION(jpts) , INTENT(in ) :: pts ! pot. temperature & salinity 605 REAL(wp), INTENT(in ) :: pdep ! depth [m] 606 REAL(wp), DIMENSION(jpts) , INTENT( out) :: pab ! thermal/haline expansion ratio 607 ! 608 REAL(wp) :: zt , zh , zs ! local scalars 609 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 610 !!---------------------------------------------------------------------- 611 ! 612 IF( nn_timing == 1 ) CALL timing_start('rab_2d') 613 ! 614 pab(:) = 0._wp 615 ! 616 SELECT CASE ( nn_eos ) 617 ! 618 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 619 ! 620 ! 621 zh = pdep * r1_Z0 ! depth 622 zt = pts (jp_tem) * r1_T0 ! temperature 623 zs = SQRT( ABS( pts(jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 624 ! 625 ! alpha 626 zn3 = ALP003 627 ! 628 zn2 = ALP012*zt + ALP102*zs+ALP002 629 ! 630 zn1 = ((ALP031*zt & 631 & + ALP121*zs+ALP021)*zt & 632 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 633 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 634 ! 635 zn0 = ((((ALP050*zt & 636 & + ALP140*zs+ALP040)*zt & 637 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 638 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 639 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 640 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 641 ! 642 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 643 ! 644 pab(jp_tem) = zn * r1_rau0 645 ! 646 ! beta 647 zn3 = BET003 648 ! 649 zn2 = BET012*zt + BET102*zs+BET002 650 ! 651 zn1 = ((BET031*zt & 652 & + BET121*zs+BET021)*zt & 653 & + (BET211*zs+BET111)*zs+BET011)*zt & 654 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 655 ! 656 zn0 = ((((BET050*zt & 657 & + BET140*zs+BET040)*zt & 658 & + (BET230*zs+BET130)*zs+BET030)*zt & 659 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 660 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 661 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 662 ! 663 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 664 ! 665 pab(jp_sal) = zn / zs * r1_rau0 666 ! 667 ! 668 ! 669 CASE( 1 ) !== simplified EOS ==! 670 ! 671 zt = pts(jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 672 zs = pts(jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 673 zh = pdep ! depth at the partial step level 674 ! 675 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 676 pab(jp_tem) = zn * r1_rau0 ! alpha 677 ! 678 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 679 pab(jp_sal) = zn * r1_rau0 ! beta 680 ! 681 CASE DEFAULT 682 IF(lwp) WRITE(numout,cform_err) 683 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 684 nstop = nstop + 1 685 ! 686 END SELECT 687 ! 688 IF( nn_timing == 1 ) CALL timing_stop('rab_2d') 689 ! 690 END SUBROUTINE rab_crs_0d 691 692 693 SUBROUTINE bn2_crs( pts, pab, pn2 ) 694 !!---------------------------------------------------------------------- 695 !! *** ROUTINE bn2 *** 696 !! 697 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the 698 !! time-step of the input arguments 699 !! 700 !! ** Method : pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w 701 !! where alpha and beta are given in pab, and computed on T-points. 702 !! N.B. N^2 is set one for all to zero at jk=1 in istate module. 703 !! 704 !! ** Action : pn2 : square of the brunt-vaisala frequency at w-point 705 !! 706 !!---------------------------------------------------------------------- 707 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celcius,psu] 708 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk,jpts), INTENT(in ) :: pab ! thermal/haline expansion coef. [Celcius-1,psu-1] 709 REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk ), INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2] 710 ! 711 INTEGER :: ji, jj, jk ! dummy loop indices 712 REAL(wp) :: zaw, zbw, zrw ! local scalars 713 !!---------------------------------------------------------------------- 714 ! 715 pn2(:,:,:)=0._wp 716 717 IF( nn_timing == 1 ) CALL timing_start('bn2') 718 ! 719 DO jk = 2, jpkm1 ! interior points only (2=< jk =< jpkm1 ) 720 DO jj = 1, jpj_crs ! surface and bottom value set to zero one for all in istate.F90 721 DO ji = 1, jpi_crs 722 !zrw = ( gdepw_crs(ji,jj,jk ) - gdept_crs(ji,jj,jk) ) & 723 ! & / ( gdept_crs(ji,jj,jk-1) - gdept_crs(ji,jj,jk) ) 724 zrw = gdepw_crs(ji,jj,jk ) - gdept_crs(ji,jj,jk) 725 !?IF( gdept_crs(ji,jj,jk-1) - gdept_crs(ji,jj,jk) .NE. 0._wp )THEN 726 IF( gdept_crs(ji,jj,jk-1) - gdept_crs(ji,jj,jk) .LT. 0._wp )THEN 727 zrw = zrw / ( gdept_crs(ji,jj,jk-1) - gdept_crs(ji,jj,jk) ) 728 ELSE 729 zrw = 0._wp 730 ENDIF 731 ! 732 zaw = pab(ji,jj,jk,jp_tem) * (1._wp - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 733 zbw = pab(ji,jj,jk,jp_sal) * (1._wp - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 734 ! 735 IF( e3w_max_crs(ji,jj,jk) .NE. 0._wp ) THEN 736 pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 737 & - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) & 738 & * tmask_crs(ji,jj,jk) / e3w_max_crs(ji,jj,jk) 739 ENDIF 740 END DO 741 END DO 742 END DO 743 ! 744 IF( nn_timing == 1 ) CALL timing_stop('bn2') 745 ! 746 END SUBROUTINE bn2_crs 747 748 SUBROUTINE eos_init_crs 707 749 !!---------------------------------------------------------------------- 708 750 !! *** ROUTINE eos_init *** 709 751 !! 710 !! ** Purpose : Compute the sea surface freezing temperature [Celcius]711 !!712 !! ** Method : UNESCO freezing point at the surface (pressure = 0???)713 !! freezing point [Celcius]=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p714 !! checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars715 !!716 !! Reference : UNESCO tech. papers in the marine science no. 28. 1978717 !!----------------------------------------------------------------------718 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]719 ! Leave result array automatic rather than making explicitly allocated720 REAL(wp), DIMENSION(jpi,jpj) :: ptf ! freezing temperature [Celcius]721 !!----------------------------------------------------------------------722 !723 ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) ) &724 & - 2.154996e-4_wp * psal(:,:) ) * psal(:,:)725 !726 END FUNCTION tfreez_crs727 728 729 SUBROUTINE eos_init_crs730 !!----------------------------------------------------------------------731 !! *** ROUTINE eos_init ***732 !!733 752 !! ** Purpose : initializations for the equation of state 734 753 !! 735 754 !! ** Method : Read the namelist nameos and control the parameters 736 755 !!---------------------------------------------------------------------- 737 INTEGER :: ios ! Local integer output status for namelist read 738 !! 739 NAMELIST/nameos/ nn_eos, rn_alpha, rn_beta 740 !!---------------------------------------------------------------------- 741 ! 742 REWIND( numnam_ref ) 743 READ ( numnam_ref, nameos, IOSTAT = ios, ERR = 901) 756 INTEGER :: ios ! local integer 757 !! 758 NAMELIST/nameos/ nn_eos, ln_useCT, rn_a0, rn_b0, rn_lambda1, rn_mu1, & 759 & rn_lambda2, rn_mu2, rn_nu 760 !!---------------------------------------------------------------------- 761 ! 762 REWIND( numnam_ref ) ! Namelist nameos in reference namelist : equation of state 763 READ ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 ) 744 764 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist', lwp ) 745 746 REWIND( numnam_cfg ) 765 ! 766 REWIND( numnam_cfg ) ! Namelist nameos in configuration namelist : equation of state 747 767 READ ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 ) 748 768 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp ) 749 IF(lwm) WRITE ( numond, nameos ) 750 769 IF(lwm) WRITE( numond, nameos ) 770 ! 771 rau0 = 1026._wp !: volumic mass of reference [kg/m3] 772 rcp = 3991.86795711963_wp !: heat capacity [J/K] 751 773 ! 752 774 IF(lwp) THEN ! Control print … … 756 778 WRITE(numout,*) ' Namelist nameos : set eos parameters' 757 779 WRITE(numout,*) ' flag for eq. of state and N^2 nn_eos = ', nn_eos 758 WRITE(numout,*) ' thermal exp. coef. (linear) rn_alpha = ', rn_alpha 759 WRITE(numout,*) ' saline exp. coef. (linear) rn_beta = ', rn_beta 780 IF( ln_useCT ) THEN 781 WRITE(numout,*) ' model uses Conservative Temperature' 782 WRITE(numout,*) ' Important: model must be initialized with CT and SA fields' 783 ENDIF 760 784 ENDIF 761 785 ! 762 786 SELECT CASE( nn_eos ) ! check option 763 787 ! 764 CASE( 0 ) !== Jackett and McDougall (1994) formulation==!788 CASE( -1 ) !== polynomial TEOS-10 ==! 765 789 IF(lwp) WRITE(numout,*) 766 IF(lwp) WRITE(numout,*) ' use of Jackett & McDougall (1994) equation of state and' 767 IF(lwp) WRITE(numout,*) ' McDougall (1987) Brunt-Vaisala frequency' 768 ! 769 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 790 IF(lwp) WRITE(numout,*) ' use of TEOS-10 equation of state (cons. temp. and abs. salinity)' 791 ! 792 rdeltaS = 32._wp 793 r1_S0 = 0.875_wp/35.16504_wp 794 r1_T0 = 1._wp/40._wp 795 r1_Z0 = 1.e-4_wp 796 ! 797 EOS000 = 8.0189615746e+02_wp 798 EOS100 = 8.6672408165e+02_wp 799 EOS200 = -1.7864682637e+03_wp 800 EOS300 = 2.0375295546e+03_wp 801 EOS400 = -1.2849161071e+03_wp 802 EOS500 = 4.3227585684e+02_wp 803 EOS600 = -6.0579916612e+01_wp 804 EOS010 = 2.6010145068e+01_wp 805 EOS110 = -6.5281885265e+01_wp 806 EOS210 = 8.1770425108e+01_wp 807 EOS310 = -5.6888046321e+01_wp 808 EOS410 = 1.7681814114e+01_wp 809 EOS510 = -1.9193502195_wp 810 EOS020 = -3.7074170417e+01_wp 811 EOS120 = 6.1548258127e+01_wp 812 EOS220 = -6.0362551501e+01_wp 813 EOS320 = 2.9130021253e+01_wp 814 EOS420 = -5.4723692739_wp 815 EOS030 = 2.1661789529e+01_wp 816 EOS130 = -3.3449108469e+01_wp 817 EOS230 = 1.9717078466e+01_wp 818 EOS330 = -3.1742946532_wp 819 EOS040 = -8.3627885467_wp 820 EOS140 = 1.1311538584e+01_wp 821 EOS240 = -5.3563304045_wp 822 EOS050 = 5.4048723791e-01_wp 823 EOS150 = 4.8169980163e-01_wp 824 EOS060 = -1.9083568888e-01_wp 825 EOS001 = 1.9681925209e+01_wp 826 EOS101 = -4.2549998214e+01_wp 827 EOS201 = 5.0774768218e+01_wp 828 EOS301 = -3.0938076334e+01_wp 829 EOS401 = 6.6051753097_wp 830 EOS011 = -1.3336301113e+01_wp 831 EOS111 = -4.4870114575_wp 832 EOS211 = 5.0042598061_wp 833 EOS311 = -6.5399043664e-01_wp 834 EOS021 = 6.7080479603_wp 835 EOS121 = 3.5063081279_wp 836 EOS221 = -1.8795372996_wp 837 EOS031 = -2.4649669534_wp 838 EOS131 = -5.5077101279e-01_wp 839 EOS041 = 5.5927935970e-01_wp 840 EOS002 = 2.0660924175_wp 841 EOS102 = -4.9527603989_wp 842 EOS202 = 2.5019633244_wp 843 EOS012 = 2.0564311499_wp 844 EOS112 = -2.1311365518e-01_wp 845 EOS022 = -1.2419983026_wp 846 EOS003 = -2.3342758797e-02_wp 847 EOS103 = -1.8507636718e-02_wp 848 EOS013 = 3.7969820455e-01_wp 849 ! 850 ALP000 = -6.5025362670e-01_wp 851 ALP100 = 1.6320471316_wp 852 ALP200 = -2.0442606277_wp 853 ALP300 = 1.4222011580_wp 854 ALP400 = -4.4204535284e-01_wp 855 ALP500 = 4.7983755487e-02_wp 856 ALP010 = 1.8537085209_wp 857 ALP110 = -3.0774129064_wp 858 ALP210 = 3.0181275751_wp 859 ALP310 = -1.4565010626_wp 860 ALP410 = 2.7361846370e-01_wp 861 ALP020 = -1.6246342147_wp 862 ALP120 = 2.5086831352_wp 863 ALP220 = -1.4787808849_wp 864 ALP320 = 2.3807209899e-01_wp 865 ALP030 = 8.3627885467e-01_wp 866 ALP130 = -1.1311538584_wp 867 ALP230 = 5.3563304045e-01_wp 868 ALP040 = -6.7560904739e-02_wp 869 ALP140 = -6.0212475204e-02_wp 870 ALP050 = 2.8625353333e-02_wp 871 ALP001 = 3.3340752782e-01_wp 872 ALP101 = 1.1217528644e-01_wp 873 ALP201 = -1.2510649515e-01_wp 874 ALP301 = 1.6349760916e-02_wp 875 ALP011 = -3.3540239802e-01_wp 876 ALP111 = -1.7531540640e-01_wp 877 ALP211 = 9.3976864981e-02_wp 878 ALP021 = 1.8487252150e-01_wp 879 ALP121 = 4.1307825959e-02_wp 880 ALP031 = -5.5927935970e-02_wp 881 ALP002 = -5.1410778748e-02_wp 882 ALP102 = 5.3278413794e-03_wp 883 ALP012 = 6.2099915132e-02_wp 884 ALP003 = -9.4924551138e-03_wp 885 ! 886 BET000 = 1.0783203594e+01_wp 887 BET100 = -4.4452095908e+01_wp 888 BET200 = 7.6048755820e+01_wp 889 BET300 = -6.3944280668e+01_wp 890 BET400 = 2.6890441098e+01_wp 891 BET500 = -4.5221697773_wp 892 BET010 = -8.1219372432e-01_wp 893 BET110 = 2.0346663041_wp 894 BET210 = -2.1232895170_wp 895 BET310 = 8.7994140485e-01_wp 896 BET410 = -1.1939638360e-01_wp 897 BET020 = 7.6574242289e-01_wp 898 BET120 = -1.5019813020_wp 899 BET220 = 1.0872489522_wp 900 BET320 = -2.7233429080e-01_wp 901 BET030 = -4.1615152308e-01_wp 902 BET130 = 4.9061350869e-01_wp 903 BET230 = -1.1847737788e-01_wp 904 BET040 = 1.4073062708e-01_wp 905 BET140 = -1.3327978879e-01_wp 906 BET050 = 5.9929880134e-03_wp 907 BET001 = -5.2937873009e-01_wp 908 BET101 = 1.2634116779_wp 909 BET201 = -1.1547328025_wp 910 BET301 = 3.2870876279e-01_wp 911 BET011 = -5.5824407214e-02_wp 912 BET111 = 1.2451933313e-01_wp 913 BET211 = -2.4409539932e-02_wp 914 BET021 = 4.3623149752e-02_wp 915 BET121 = -4.6767901790e-02_wp 916 BET031 = -6.8523260060e-03_wp 917 BET002 = -6.1618945251e-02_wp 918 BET102 = 6.2255521644e-02_wp 919 BET012 = -2.6514181169e-03_wp 920 BET003 = -2.3025968587e-04_wp 921 ! 922 PEN000 = -9.8409626043_wp 923 PEN100 = 2.1274999107e+01_wp 924 PEN200 = -2.5387384109e+01_wp 925 PEN300 = 1.5469038167e+01_wp 926 PEN400 = -3.3025876549_wp 927 PEN010 = 6.6681505563_wp 928 PEN110 = 2.2435057288_wp 929 PEN210 = -2.5021299030_wp 930 PEN310 = 3.2699521832e-01_wp 931 PEN020 = -3.3540239802_wp 932 PEN120 = -1.7531540640_wp 933 PEN220 = 9.3976864981e-01_wp 934 PEN030 = 1.2324834767_wp 935 PEN130 = 2.7538550639e-01_wp 936 PEN040 = -2.7963967985e-01_wp 937 PEN001 = -1.3773949450_wp 938 PEN101 = 3.3018402659_wp 939 PEN201 = -1.6679755496_wp 940 PEN011 = -1.3709540999_wp 941 PEN111 = 1.4207577012e-01_wp 942 PEN021 = 8.2799886843e-01_wp 943 PEN002 = 1.7507069098e-02_wp 944 PEN102 = 1.3880727538e-02_wp 945 PEN012 = -2.8477365341e-01_wp 946 ! 947 APE000 = -1.6670376391e-01_wp 948 APE100 = -5.6087643219e-02_wp 949 APE200 = 6.2553247576e-02_wp 950 APE300 = -8.1748804580e-03_wp 951 APE010 = 1.6770119901e-01_wp 952 APE110 = 8.7657703198e-02_wp 953 APE210 = -4.6988432490e-02_wp 954 APE020 = -9.2436260751e-02_wp 955 APE120 = -2.0653912979e-02_wp 956 APE030 = 2.7963967985e-02_wp 957 APE001 = 3.4273852498e-02_wp 958 APE101 = -3.5518942529e-03_wp 959 APE011 = -4.1399943421e-02_wp 960 APE002 = 7.1193413354e-03_wp 961 ! 962 BPE000 = 2.6468936504e-01_wp 963 BPE100 = -6.3170583896e-01_wp 964 BPE200 = 5.7736640125e-01_wp 965 BPE300 = -1.6435438140e-01_wp 966 BPE010 = 2.7912203607e-02_wp 967 BPE110 = -6.2259666565e-02_wp 968 BPE210 = 1.2204769966e-02_wp 969 BPE020 = -2.1811574876e-02_wp 970 BPE120 = 2.3383950895e-02_wp 971 BPE030 = 3.4261630030e-03_wp 972 BPE001 = 4.1079296834e-02_wp 973 BPE101 = -4.1503681096e-02_wp 974 BPE011 = 1.7676120780e-03_wp 975 BPE002 = 1.7269476440e-04_wp 976 ! 977 CASE( 0 ) !== polynomial EOS-80 formulation ==! 978 ! 770 979 IF(lwp) WRITE(numout,*) 771 IF(lwp) WRITE(numout,*) ' use of linear eos rho(T) = rau0 * ( 1.0285 - rn_alpha * T )' 772 IF( lk_zdfddm ) CALL ctl_stop( ' double diffusive mixing parameterization requires', & 773 & ' that T and S are used as state variables' ) 774 ! 775 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 776 ralpbet_crs = rn_alpha / rn_beta 777 IF(lwp) WRITE(numout,*) 778 IF(lwp) WRITE(numout,*) ' use of linear eos rho(T,S) = rau0 * ( rn_beta * S - rn_alpha * T )' 980 IF(lwp) WRITE(numout,*) ' use of EOS-80 equation of state (pot. temp. and pract. salinity)' 981 ! 982 rdeltaS = 20._wp 983 r1_S0 = 1._wp/40._wp 984 r1_T0 = 1._wp/40._wp 985 r1_Z0 = 1.e-4_wp 986 ! 987 EOS000 = 9.5356891948e+02_wp 988 EOS100 = 1.7136499189e+02_wp 989 EOS200 = -3.7501039454e+02_wp 990 EOS300 = 5.1856810420e+02_wp 991 EOS400 = -3.7264470465e+02_wp 992 EOS500 = 1.4302533998e+02_wp 993 EOS600 = -2.2856621162e+01_wp 994 EOS010 = 1.0087518651e+01_wp 995 EOS110 = -1.3647741861e+01_wp 996 EOS210 = 8.8478359933_wp 997 EOS310 = -7.2329388377_wp 998 EOS410 = 1.4774410611_wp 999 EOS510 = 2.0036720553e-01_wp 1000 EOS020 = -2.5579830599e+01_wp 1001 EOS120 = 2.4043512327e+01_wp 1002 EOS220 = -1.6807503990e+01_wp 1003 EOS320 = 8.3811577084_wp 1004 EOS420 = -1.9771060192_wp 1005 EOS030 = 1.6846451198e+01_wp 1006 EOS130 = -2.1482926901e+01_wp 1007 EOS230 = 1.0108954054e+01_wp 1008 EOS330 = -6.2675951440e-01_wp 1009 EOS040 = -8.0812310102_wp 1010 EOS140 = 1.0102374985e+01_wp 1011 EOS240 = -4.8340368631_wp 1012 EOS050 = 1.2079167803_wp 1013 EOS150 = 1.1515380987e-01_wp 1014 EOS060 = -2.4520288837e-01_wp 1015 EOS001 = 1.0748601068e+01_wp 1016 EOS101 = -1.7817043500e+01_wp 1017 EOS201 = 2.2181366768e+01_wp 1018 EOS301 = -1.6750916338e+01_wp 1019 EOS401 = 4.1202230403_wp 1020 EOS011 = -1.5852644587e+01_wp 1021 EOS111 = -7.6639383522e-01_wp 1022 EOS211 = 4.1144627302_wp 1023 EOS311 = -6.6955877448e-01_wp 1024 EOS021 = 9.9994861860_wp 1025 EOS121 = -1.9467067787e-01_wp 1026 EOS221 = -1.2177554330_wp 1027 EOS031 = -3.4866102017_wp 1028 EOS131 = 2.2229155620e-01_wp 1029 EOS041 = 5.9503008642e-01_wp 1030 EOS002 = 1.0375676547_wp 1031 EOS102 = -3.4249470629_wp 1032 EOS202 = 2.0542026429_wp 1033 EOS012 = 2.1836324814_wp 1034 EOS112 = -3.4453674320e-01_wp 1035 EOS022 = -1.2548163097_wp 1036 EOS003 = 1.8729078427e-02_wp 1037 EOS103 = -5.7238495240e-02_wp 1038 EOS013 = 3.8306136687e-01_wp 1039 ! 1040 ALP000 = -2.5218796628e-01_wp 1041 ALP100 = 3.4119354654e-01_wp 1042 ALP200 = -2.2119589983e-01_wp 1043 ALP300 = 1.8082347094e-01_wp 1044 ALP400 = -3.6936026529e-02_wp 1045 ALP500 = -5.0091801383e-03_wp 1046 ALP010 = 1.2789915300_wp 1047 ALP110 = -1.2021756164_wp 1048 ALP210 = 8.4037519952e-01_wp 1049 ALP310 = -4.1905788542e-01_wp 1050 ALP410 = 9.8855300959e-02_wp 1051 ALP020 = -1.2634838399_wp 1052 ALP120 = 1.6112195176_wp 1053 ALP220 = -7.5817155402e-01_wp 1054 ALP320 = 4.7006963580e-02_wp 1055 ALP030 = 8.0812310102e-01_wp 1056 ALP130 = -1.0102374985_wp 1057 ALP230 = 4.8340368631e-01_wp 1058 ALP040 = -1.5098959754e-01_wp 1059 ALP140 = -1.4394226233e-02_wp 1060 ALP050 = 3.6780433255e-02_wp 1061 ALP001 = 3.9631611467e-01_wp 1062 ALP101 = 1.9159845880e-02_wp 1063 ALP201 = -1.0286156825e-01_wp 1064 ALP301 = 1.6738969362e-02_wp 1065 ALP011 = -4.9997430930e-01_wp 1066 ALP111 = 9.7335338937e-03_wp 1067 ALP211 = 6.0887771651e-02_wp 1068 ALP021 = 2.6149576513e-01_wp 1069 ALP121 = -1.6671866715e-02_wp 1070 ALP031 = -5.9503008642e-02_wp 1071 ALP002 = -5.4590812035e-02_wp 1072 ALP102 = 8.6134185799e-03_wp 1073 ALP012 = 6.2740815484e-02_wp 1074 ALP003 = -9.5765341718e-03_wp 1075 ! 1076 BET000 = 2.1420623987_wp 1077 BET100 = -9.3752598635_wp 1078 BET200 = 1.9446303907e+01_wp 1079 BET300 = -1.8632235232e+01_wp 1080 BET400 = 8.9390837485_wp 1081 BET500 = -1.7142465871_wp 1082 BET010 = -1.7059677327e-01_wp 1083 BET110 = 2.2119589983e-01_wp 1084 BET210 = -2.7123520642e-01_wp 1085 BET310 = 7.3872053057e-02_wp 1086 BET410 = 1.2522950346e-02_wp 1087 BET020 = 3.0054390409e-01_wp 1088 BET120 = -4.2018759976e-01_wp 1089 BET220 = 3.1429341406e-01_wp 1090 BET320 = -9.8855300959e-02_wp 1091 BET030 = -2.6853658626e-01_wp 1092 BET130 = 2.5272385134e-01_wp 1093 BET230 = -2.3503481790e-02_wp 1094 BET040 = 1.2627968731e-01_wp 1095 BET140 = -1.2085092158e-01_wp 1096 BET050 = 1.4394226233e-03_wp 1097 BET001 = -2.2271304375e-01_wp 1098 BET101 = 5.5453416919e-01_wp 1099 BET201 = -6.2815936268e-01_wp 1100 BET301 = 2.0601115202e-01_wp 1101 BET011 = -9.5799229402e-03_wp 1102 BET111 = 1.0286156825e-01_wp 1103 BET211 = -2.5108454043e-02_wp 1104 BET021 = -2.4333834734e-03_wp 1105 BET121 = -3.0443885826e-02_wp 1106 BET031 = 2.7786444526e-03_wp 1107 BET002 = -4.2811838287e-02_wp 1108 BET102 = 5.1355066072e-02_wp 1109 BET012 = -4.3067092900e-03_wp 1110 BET003 = -7.1548119050e-04_wp 1111 ! 1112 PEN000 = -5.3743005340_wp 1113 PEN100 = 8.9085217499_wp 1114 PEN200 = -1.1090683384e+01_wp 1115 PEN300 = 8.3754581690_wp 1116 PEN400 = -2.0601115202_wp 1117 PEN010 = 7.9263222935_wp 1118 PEN110 = 3.8319691761e-01_wp 1119 PEN210 = -2.0572313651_wp 1120 PEN310 = 3.3477938724e-01_wp 1121 PEN020 = -4.9997430930_wp 1122 PEN120 = 9.7335338937e-02_wp 1123 PEN220 = 6.0887771651e-01_wp 1124 PEN030 = 1.7433051009_wp 1125 PEN130 = -1.1114577810e-01_wp 1126 PEN040 = -2.9751504321e-01_wp 1127 PEN001 = -6.9171176978e-01_wp 1128 PEN101 = 2.2832980419_wp 1129 PEN201 = -1.3694684286_wp 1130 PEN011 = -1.4557549876_wp 1131 PEN111 = 2.2969116213e-01_wp 1132 PEN021 = 8.3654420645e-01_wp 1133 PEN002 = -1.4046808820e-02_wp 1134 PEN102 = 4.2928871430e-02_wp 1135 PEN012 = -2.8729602515e-01_wp 1136 ! 1137 APE000 = -1.9815805734e-01_wp 1138 APE100 = -9.5799229402e-03_wp 1139 APE200 = 5.1430784127e-02_wp 1140 APE300 = -8.3694846809e-03_wp 1141 APE010 = 2.4998715465e-01_wp 1142 APE110 = -4.8667669469e-03_wp 1143 APE210 = -3.0443885826e-02_wp 1144 APE020 = -1.3074788257e-01_wp 1145 APE120 = 8.3359333577e-03_wp 1146 APE030 = 2.9751504321e-02_wp 1147 APE001 = 3.6393874690e-02_wp 1148 APE101 = -5.7422790533e-03_wp 1149 APE011 = -4.1827210323e-02_wp 1150 APE002 = 7.1824006288e-03_wp 1151 ! 1152 BPE000 = 1.1135652187e-01_wp 1153 BPE100 = -2.7726708459e-01_wp 1154 BPE200 = 3.1407968134e-01_wp 1155 BPE300 = -1.0300557601e-01_wp 1156 BPE010 = 4.7899614701e-03_wp 1157 BPE110 = -5.1430784127e-02_wp 1158 BPE210 = 1.2554227021e-02_wp 1159 BPE020 = 1.2166917367e-03_wp 1160 BPE120 = 1.5221942913e-02_wp 1161 BPE030 = -1.3893222263e-03_wp 1162 BPE001 = 2.8541225524e-02_wp 1163 BPE101 = -3.4236710714e-02_wp 1164 BPE011 = 2.8711395266e-03_wp 1165 BPE002 = 5.3661089288e-04_wp 1166 ! 1167 CASE( 1 ) !== Simplified EOS ==! 1168 IF(lwp) THEN 1169 WRITE(numout,*) 1170 WRITE(numout,*) ' use of simplified eos: rhd(dT=T-10,dS=S-35,Z) = ' 1171 WRITE(numout,*) ' [-a0*(1+lambda1/2*dT+mu1*Z)*dT + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS]/rau0' 1172 WRITE(numout,*) 1173 WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0 1174 WRITE(numout,*) ' saline cont. coef. rn_b0 = ', rn_b0 1175 WRITE(numout,*) ' cabbeling coef. rn_lambda1 = ', rn_lambda1 1176 WRITE(numout,*) ' cabbeling coef. rn_lambda2 = ', rn_lambda2 1177 WRITE(numout,*) ' thermobar. coef. rn_mu1 = ', rn_mu1 1178 WRITE(numout,*) ' thermobar. coef. rn_mu2 = ', rn_mu2 1179 WRITE(numout,*) ' 2nd cabbel. coef. rn_nu = ', rn_nu 1180 WRITE(numout,*) ' Caution: rn_beta0=0 incompatible with ddm parameterization ' 1181 ENDIF 779 1182 ! 780 1183 CASE DEFAULT !== ERROR in nn_eos ==! … … 784 1187 END SELECT 785 1188 ! 1189 r1_rau0 = 1._wp / rau0 1190 r1_rcp = 1._wp / rcp 1191 r1_rau0_rcp = 1._wp / ( rau0 * rcp ) 1192 ! 1193 IF(lwp) WRITE(numout,*) 1194 IF(lwp) WRITE(numout,*) ' volumic mass of reference rau0 = ', rau0 , ' kg/m^3' 1195 IF(lwp) WRITE(numout,*) ' 1. / rau0 r1_rau0 = ', r1_rau0, ' m^3/kg' 1196 IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin' 1197 IF(lwp) WRITE(numout,*) ' 1. / ( rau0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp 1198 ! 786 1199 END SUBROUTINE eos_init_crs 787 1200
Note: See TracChangeset
for help on using the changeset viewer.
