- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/eosbn2.F90
r12178 r12928 29 29 !! eos_insitu_pot: Compute the insitu and surface referenced potential volumic mass 30 30 !! eos_insitu_2d : Compute the in situ density for 2d fields 31 !! bn2 : Compute the Brunt-Vaisala frequency32 31 !! bn2 : compute the Brunt-Vaisala frequency 33 32 !! eos_pt_from_ct: compute the potential temperature from the Conservative Temperature … … 180 179 181 180 !! * Substitutions 182 # include " vectopt_loop_substitute.h90"181 # include "do_loop_substitute.h90" 183 182 !!---------------------------------------------------------------------- 184 183 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 192 191 !! *** ROUTINE eos_insitu *** 193 192 !! 194 !! ** Purpose : Compute the in situ density (ratio rho/r au0) from193 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 195 194 !! potential temperature and salinity using an equation of state 196 195 !! selected in the nameos namelist 197 196 !! 198 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - r au0 ) / rau0197 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 199 198 !! with prd in situ density anomaly no units 200 199 !! t TEOS10: CT or EOS80: PT Celsius … … 202 201 !! z depth meters 203 202 !! rho in situ density kg/m^3 204 !! r au0 reference density kg/m^3203 !! rho0 reference density kg/m^3 205 204 !! 206 205 !! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). … … 211 210 !! 212 211 !! ln_seos : simplified equation of state 213 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / r au0212 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0 214 213 !! linear case function of T only: rn_alpha<>0, other coefficients = 0 215 214 !! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 … … 238 237 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 239 238 ! 240 DO jk = 1, jpkm1 241 DO jj = 1, jpj 242 DO ji = 1, jpi 243 ! 244 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 245 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 246 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 247 ztm = tmask(ji,jj,jk) ! tmask 239 DO_3D_11_11( 1, jpkm1 ) 240 ! 241 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 242 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 243 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 244 ztm = tmask(ji,jj,jk) ! tmask 245 ! 246 zn3 = EOS013*zt & 247 & + EOS103*zs+EOS003 248 ! 249 zn2 = (EOS022*zt & 250 & + EOS112*zs+EOS012)*zt & 251 & + (EOS202*zs+EOS102)*zs+EOS002 252 ! 253 zn1 = (((EOS041*zt & 254 & + EOS131*zs+EOS031)*zt & 255 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 256 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 257 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 258 ! 259 zn0 = (((((EOS060*zt & 260 & + EOS150*zs+EOS050)*zt & 261 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 262 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 263 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 264 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 265 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 266 ! 267 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 268 ! 269 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 270 ! 271 END_3D 272 ! 273 CASE( np_seos ) !== simplified EOS ==! 274 ! 275 DO_3D_11_11( 1, jpkm1 ) 276 zt = pts (ji,jj,jk,jp_tem) - 10._wp 277 zs = pts (ji,jj,jk,jp_sal) - 35._wp 278 zh = pdep (ji,jj,jk) 279 ztm = tmask(ji,jj,jk) 280 ! 281 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt & 282 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 283 & - rn_nu * zt * zs 284 ! 285 prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked) 286 END_3D 287 ! 288 END SELECT 289 ! 290 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu : ', kdim=jpk ) 291 ! 292 IF( ln_timing ) CALL timing_stop('eos-insitu') 293 ! 294 END SUBROUTINE eos_insitu 295 296 297 SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 298 !!---------------------------------------------------------------------- 299 !! *** ROUTINE eos_insitu_pot *** 300 !! 301 !! ** Purpose : Compute the in situ density (ratio rho/rho0) and the 302 !! potential volumic mass (Kg/m3) from potential temperature and 303 !! salinity fields using an equation of state selected in the 304 !! namelist. 305 !! 306 !! ** Action : - prd , the in situ density (no units) 307 !! - prhop, the potential volumic mass (Kg/m3) 308 !! 309 !!---------------------------------------------------------------------- 310 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 311 ! ! 2 : salinity [psu] 312 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-] 313 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prhop ! potential density (surface referenced) 314 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m] 315 ! 316 INTEGER :: ji, jj, jk, jsmp ! dummy loop indices 317 INTEGER :: jdof 318 REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars 319 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 320 REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors 321 !!---------------------------------------------------------------------- 322 ! 323 IF( ln_timing ) CALL timing_start('eos-pot') 324 ! 325 SELECT CASE ( neos ) 326 ! 327 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 328 ! 329 ! Stochastic equation of state 330 IF ( ln_sto_eos ) THEN 331 ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 332 ALLOCATE(zn_sto(1:2*nn_sto_eos)) 333 ALLOCATE(zsign(1:2*nn_sto_eos)) 334 DO jsmp = 1, 2*nn_sto_eos, 2 335 zsign(jsmp) = 1._wp 336 zsign(jsmp+1) = -1._wp 337 END DO 338 ! 339 DO_3D_11_11( 1, jpkm1 ) 340 ! 341 ! compute density (2*nn_sto_eos) times: 342 ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 343 ! (2) for t-dt, s-ds (with the opposite fluctuation) 344 DO jsmp = 1, nn_sto_eos*2 345 jdof = (jsmp + 1) / 2 346 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 347 zt = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0 ! temperature 348 zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 349 zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity 350 ztm = tmask(ji,jj,jk) ! tmask 248 351 ! 249 352 zn3 = EOS013*zt & … … 260 363 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 261 364 ! 262 zn0 = (((((EOS060*zt &365 zn0_sto(jsmp) = (((((EOS060*zt & 263 366 & + EOS150*zs+EOS050)*zt & 264 367 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & … … 268 371 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 269 372 ! 270 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 373 zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 374 END DO 375 ! 376 ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 377 prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 378 DO jsmp = 1, nn_sto_eos*2 379 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 271 380 ! 272 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 273 ! 381 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rho0 - 1._wp ) ! density anomaly (masked) 274 382 END DO 275 END DO 276 END DO 277 ! 278 CASE( np_seos ) !== simplified EOS ==! 279 ! 280 DO jk = 1, jpkm1 281 DO jj = 1, jpj 282 DO ji = 1, jpi 283 zt = pts (ji,jj,jk,jp_tem) - 10._wp 284 zs = pts (ji,jj,jk,jp_sal) - 35._wp 285 zh = pdep (ji,jj,jk) 286 ztm = tmask(ji,jj,jk) 287 ! 288 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt & 289 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 290 & - rn_nu * zt * zs 291 ! 292 prd(ji,jj,jk) = zn * r1_rau0 * ztm ! density anomaly (masked) 293 END DO 294 END DO 295 END DO 296 ! 297 END SELECT 298 ! 299 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu : ', kdim=jpk ) 300 ! 301 IF( ln_timing ) CALL timing_stop('eos-insitu') 302 ! 303 END SUBROUTINE eos_insitu 304 305 306 SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 307 !!---------------------------------------------------------------------- 308 !! *** ROUTINE eos_insitu_pot *** 309 !! 310 !! ** Purpose : Compute the in situ density (ratio rho/rau0) and the 311 !! potential volumic mass (Kg/m3) from potential temperature and 312 !! salinity fields using an equation of state selected in the 313 !! namelist. 314 !! 315 !! ** Action : - prd , the in situ density (no units) 316 !! - prhop, the potential volumic mass (Kg/m3) 317 !! 318 !!---------------------------------------------------------------------- 319 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 320 ! ! 2 : salinity [psu] 321 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-] 322 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prhop ! potential density (surface referenced) 323 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m] 324 ! 325 INTEGER :: ji, jj, jk, jsmp ! dummy loop indices 326 INTEGER :: jdof 327 REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars 328 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 329 REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors 330 !!---------------------------------------------------------------------- 331 ! 332 IF( ln_timing ) CALL timing_start('eos-pot') 333 ! 334 SELECT CASE ( neos ) 335 ! 336 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 337 ! 338 ! Stochastic equation of state 339 IF ( ln_sto_eos ) THEN 340 ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 341 ALLOCATE(zn_sto(1:2*nn_sto_eos)) 342 ALLOCATE(zsign(1:2*nn_sto_eos)) 343 DO jsmp = 1, 2*nn_sto_eos, 2 344 zsign(jsmp) = 1._wp 345 zsign(jsmp+1) = -1._wp 346 END DO 347 ! 348 DO jk = 1, jpkm1 349 DO jj = 1, jpj 350 DO ji = 1, jpi 351 ! 352 ! compute density (2*nn_sto_eos) times: 353 ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 354 ! (2) for t-dt, s-ds (with the opposite fluctuation) 355 DO jsmp = 1, nn_sto_eos*2 356 jdof = (jsmp + 1) / 2 357 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 358 zt = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0 ! temperature 359 zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 360 zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity 361 ztm = tmask(ji,jj,jk) ! tmask 362 ! 363 zn3 = EOS013*zt & 364 & + EOS103*zs+EOS003 365 ! 366 zn2 = (EOS022*zt & 367 & + EOS112*zs+EOS012)*zt & 368 & + (EOS202*zs+EOS102)*zs+EOS002 369 ! 370 zn1 = (((EOS041*zt & 371 & + EOS131*zs+EOS031)*zt & 372 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 373 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 374 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 375 ! 376 zn0_sto(jsmp) = (((((EOS060*zt & 377 & + EOS150*zs+EOS050)*zt & 378 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 379 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 380 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 381 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 382 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 383 ! 384 zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 385 END DO 386 ! 387 ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 388 prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 389 DO jsmp = 1, nn_sto_eos*2 390 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 391 ! 392 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rau0 - 1._wp ) ! density anomaly (masked) 393 END DO 394 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 395 prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos 396 END DO 397 END DO 398 END DO 383 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 384 prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos 385 END_3D 399 386 DEALLOCATE(zn0_sto,zn_sto,zsign) 400 387 ! Non-stochastic equation of state 401 388 ELSE 402 DO jk = 1, jpkm1 403 DO jj = 1, jpj 404 DO ji = 1, jpi 405 ! 406 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 407 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 408 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 409 ztm = tmask(ji,jj,jk) ! tmask 410 ! 411 zn3 = EOS013*zt & 412 & + EOS103*zs+EOS003 413 ! 414 zn2 = (EOS022*zt & 415 & + EOS112*zs+EOS012)*zt & 416 & + (EOS202*zs+EOS102)*zs+EOS002 417 ! 418 zn1 = (((EOS041*zt & 419 & + EOS131*zs+EOS031)*zt & 420 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 421 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 422 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 423 ! 424 zn0 = (((((EOS060*zt & 425 & + EOS150*zs+EOS050)*zt & 426 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 427 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 428 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 429 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 430 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 431 ! 432 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 433 ! 434 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 435 ! 436 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 437 END DO 438 END DO 439 END DO 440 ENDIF 441 442 CASE( np_seos ) !== simplified EOS ==! 443 ! 444 DO jk = 1, jpkm1 445 DO jj = 1, jpj 446 DO ji = 1, jpi 447 zt = pts (ji,jj,jk,jp_tem) - 10._wp 448 zs = pts (ji,jj,jk,jp_sal) - 35._wp 449 zh = pdep (ji,jj,jk) 450 ztm = tmask(ji,jj,jk) 451 ! ! potential density referenced at the surface 452 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt & 453 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs & 454 & - rn_nu * zt * zs 455 prhop(ji,jj,jk) = ( rau0 + zn ) * ztm 456 ! ! density anomaly (masked) 457 zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 458 prd(ji,jj,jk) = zn * r1_rau0 * ztm 459 ! 460 END DO 461 END DO 462 END DO 463 ! 464 END SELECT 465 ! 466 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk ) 467 ! 468 IF( ln_timing ) CALL timing_stop('eos-pot') 469 ! 470 END SUBROUTINE eos_insitu_pot 471 472 473 SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 474 !!---------------------------------------------------------------------- 475 !! *** ROUTINE eos_insitu_2d *** 476 !! 477 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 478 !! potential temperature and salinity using an equation of state 479 !! selected in the nameos namelist. * 2D field case 480 !! 481 !! ** Action : - prd , the in situ density (no units) (unmasked) 482 !! 483 !!---------------------------------------------------------------------- 484 REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 485 ! ! 2 : salinity [psu] 486 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 487 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 488 ! 489 INTEGER :: ji, jj, jk ! dummy loop indices 490 REAL(wp) :: zt , zh , zs ! local scalars 491 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 492 !!---------------------------------------------------------------------- 493 ! 494 IF( ln_timing ) CALL timing_start('eos2d') 495 ! 496 prd(:,:) = 0._wp 497 ! 498 SELECT CASE( neos ) 499 ! 500 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 501 ! 502 DO jj = 1, jpjm1 503 DO ji = 1, fs_jpim1 ! vector opt. 504 ! 505 zh = pdep(ji,jj) * r1_Z0 ! depth 506 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 507 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 389 DO_3D_11_11( 1, jpkm1 ) 390 ! 391 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 392 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 393 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 394 ztm = tmask(ji,jj,jk) ! tmask 508 395 ! 509 396 zn3 = EOS013*zt & … … 530 417 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 531 418 ! 532 prd(ji,jj) = zn * r1_rau0 - 1._wp ! unmasked in situ density anomaly 533 ! 534 END DO 535 END DO 536 ! 537 CALL lbc_lnk( 'eosbn2', prd, 'T', 1. ) ! Lateral boundary conditions 538 ! 419 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 420 ! 421 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 422 END_3D 423 ENDIF 424 539 425 CASE( np_seos ) !== simplified EOS ==! 540 426 ! 541 DO jj = 1, jpjm1 542 DO ji = 1, fs_jpim1 ! vector opt. 543 ! 544 zt = pts (ji,jj,jp_tem) - 10._wp 545 zs = pts (ji,jj,jp_sal) - 35._wp 546 zh = pdep (ji,jj) ! depth at the partial step level 547 ! 548 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt & 549 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 550 & - rn_nu * zt * zs 551 ! 552 prd(ji,jj) = zn * r1_rau0 ! unmasked in situ density anomaly 553 ! 554 END DO 555 END DO 556 ! 557 CALL lbc_lnk( 'eosbn2', prd, 'T', 1. ) ! Lateral boundary conditions 427 DO_3D_11_11( 1, jpkm1 ) 428 zt = pts (ji,jj,jk,jp_tem) - 10._wp 429 zs = pts (ji,jj,jk,jp_sal) - 35._wp 430 zh = pdep (ji,jj,jk) 431 ztm = tmask(ji,jj,jk) 432 ! ! potential density referenced at the surface 433 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt & 434 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs & 435 & - rn_nu * zt * zs 436 prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 437 ! ! density anomaly (masked) 438 zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 439 prd(ji,jj,jk) = zn * r1_rho0 * ztm 440 ! 441 END_3D 558 442 ! 559 443 END SELECT 560 444 ! 561 IF(ln_ctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 445 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk ) 446 ! 447 IF( ln_timing ) CALL timing_stop('eos-pot') 448 ! 449 END SUBROUTINE eos_insitu_pot 450 451 452 SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 453 !!---------------------------------------------------------------------- 454 !! *** ROUTINE eos_insitu_2d *** 455 !! 456 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 457 !! potential temperature and salinity using an equation of state 458 !! selected in the nameos namelist. * 2D field case 459 !! 460 !! ** Action : - prd , the in situ density (no units) (unmasked) 461 !! 462 !!---------------------------------------------------------------------- 463 REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 464 ! ! 2 : salinity [psu] 465 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 466 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 467 ! 468 INTEGER :: ji, jj, jk ! dummy loop indices 469 REAL(wp) :: zt , zh , zs ! local scalars 470 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 471 !!---------------------------------------------------------------------- 472 ! 473 IF( ln_timing ) CALL timing_start('eos2d') 474 ! 475 prd(:,:) = 0._wp 476 ! 477 SELECT CASE( neos ) 478 ! 479 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 480 ! 481 DO_2D_11_11 482 ! 483 zh = pdep(ji,jj) * r1_Z0 ! depth 484 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 485 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 486 ! 487 zn3 = EOS013*zt & 488 & + EOS103*zs+EOS003 489 ! 490 zn2 = (EOS022*zt & 491 & + EOS112*zs+EOS012)*zt & 492 & + (EOS202*zs+EOS102)*zs+EOS002 493 ! 494 zn1 = (((EOS041*zt & 495 & + EOS131*zs+EOS031)*zt & 496 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 497 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 498 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 499 ! 500 zn0 = (((((EOS060*zt & 501 & + EOS150*zs+EOS050)*zt & 502 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 503 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 504 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 505 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 506 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 507 ! 508 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 509 ! 510 prd(ji,jj) = zn * r1_rho0 - 1._wp ! unmasked in situ density anomaly 511 ! 512 END_2D 513 ! 514 CASE( np_seos ) !== simplified EOS ==! 515 ! 516 DO_2D_11_11 517 ! 518 zt = pts (ji,jj,jp_tem) - 10._wp 519 zs = pts (ji,jj,jp_sal) - 35._wp 520 zh = pdep (ji,jj) ! depth at the partial step level 521 ! 522 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt & 523 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 524 & - rn_nu * zt * zs 525 ! 526 prd(ji,jj) = zn * r1_rho0 ! unmasked in situ density anomaly 527 ! 528 END_2D 529 ! 530 END SELECT 531 ! 532 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 562 533 ! 563 534 IF( ln_timing ) CALL timing_stop('eos2d') … … 566 537 567 538 568 SUBROUTINE rab_3d( pts, pab )539 SUBROUTINE rab_3d( pts, pab, Kmm ) 569 540 !!---------------------------------------------------------------------- 570 541 !! *** ROUTINE rab_3d *** … … 576 547 !! ** Action : - pab : thermal/haline expansion ratio at T-points 577 548 !!---------------------------------------------------------------------- 549 INTEGER , INTENT(in ) :: Kmm ! time level index 578 550 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 579 551 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab ! thermal/haline expansion ratio … … 590 562 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 591 563 ! 592 DO jk = 1, jpkm1 593 DO jj = 1, jpj 594 DO ji = 1, jpi 595 ! 596 zh = gdept_n(ji,jj,jk) * r1_Z0 ! depth 597 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 598 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 599 ztm = tmask(ji,jj,jk) ! tmask 600 ! 601 ! alpha 602 zn3 = ALP003 603 ! 604 zn2 = ALP012*zt + ALP102*zs+ALP002 605 ! 606 zn1 = ((ALP031*zt & 607 & + ALP121*zs+ALP021)*zt & 608 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 609 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 610 ! 611 zn0 = ((((ALP050*zt & 612 & + ALP140*zs+ALP040)*zt & 613 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 614 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 615 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 616 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 617 ! 618 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 619 ! 620 pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm 621 ! 622 ! beta 623 zn3 = BET003 624 ! 625 zn2 = BET012*zt + BET102*zs+BET002 626 ! 627 zn1 = ((BET031*zt & 628 & + BET121*zs+BET021)*zt & 629 & + (BET211*zs+BET111)*zs+BET011)*zt & 630 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 631 ! 632 zn0 = ((((BET050*zt & 633 & + BET140*zs+BET040)*zt & 634 & + (BET230*zs+BET130)*zs+BET030)*zt & 635 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 636 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 637 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 638 ! 639 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 640 ! 641 pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm 642 ! 643 END DO 644 END DO 645 END DO 564 DO_3D_11_11( 1, jpkm1 ) 565 ! 566 zh = gdept(ji,jj,jk,Kmm) * r1_Z0 ! depth 567 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 568 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 569 ztm = tmask(ji,jj,jk) ! tmask 570 ! 571 ! alpha 572 zn3 = ALP003 573 ! 574 zn2 = ALP012*zt + ALP102*zs+ALP002 575 ! 576 zn1 = ((ALP031*zt & 577 & + ALP121*zs+ALP021)*zt & 578 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 579 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 580 ! 581 zn0 = ((((ALP050*zt & 582 & + ALP140*zs+ALP040)*zt & 583 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 584 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 585 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 586 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 587 ! 588 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 589 ! 590 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 591 ! 592 ! beta 593 zn3 = BET003 594 ! 595 zn2 = BET012*zt + BET102*zs+BET002 596 ! 597 zn1 = ((BET031*zt & 598 & + BET121*zs+BET021)*zt & 599 & + (BET211*zs+BET111)*zs+BET011)*zt & 600 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 601 ! 602 zn0 = ((((BET050*zt & 603 & + BET140*zs+BET040)*zt & 604 & + (BET230*zs+BET130)*zs+BET030)*zt & 605 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 606 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 607 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 608 ! 609 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 610 ! 611 pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 612 ! 613 END_3D 646 614 ! 647 615 CASE( np_seos ) !== simplified EOS ==! 648 616 ! 649 DO jk = 1, jpkm1 650 DO jj = 1, jpj 651 DO ji = 1, jpi 652 zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 653 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 654 zh = gdept_n(ji,jj,jk) ! depth in meters at t-point 655 ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask 656 ! 657 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 658 pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm ! alpha 659 ! 660 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 661 pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm ! beta 662 ! 663 END DO 664 END DO 665 END DO 617 DO_3D_11_11( 1, jpkm1 ) 618 zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 619 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 620 zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point 621 ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask 622 ! 623 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 624 pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm ! alpha 625 ! 626 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 627 pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm ! beta 628 ! 629 END_3D 666 630 ! 667 631 CASE DEFAULT … … 671 635 END SELECT 672 636 ! 673 IF( ln_ctl) CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &674 & tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk )637 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', & 638 & tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk ) 675 639 ! 676 640 IF( ln_timing ) CALL timing_stop('rab_3d') … … 679 643 680 644 681 SUBROUTINE rab_2d( pts, pdep, pab )645 SUBROUTINE rab_2d( pts, pdep, pab, Kmm ) 682 646 !!---------------------------------------------------------------------- 683 647 !! *** ROUTINE rab_2d *** … … 687 651 !! ** Action : - pab : thermal/haline expansion ratio at T-points 688 652 !!---------------------------------------------------------------------- 653 INTEGER , INTENT(in ) :: Kmm ! time level index 689 654 REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT(in ) :: pts ! pot. temperature & salinity 690 655 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] … … 704 669 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 705 670 ! 706 DO jj = 1, jpjm1 707 DO ji = 1, fs_jpim1 ! vector opt. 708 ! 709 zh = pdep(ji,jj) * r1_Z0 ! depth 710 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 711 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 712 ! 713 ! alpha 714 zn3 = ALP003 715 ! 716 zn2 = ALP012*zt + ALP102*zs+ALP002 717 ! 718 zn1 = ((ALP031*zt & 719 & + ALP121*zs+ALP021)*zt & 720 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 721 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 722 ! 723 zn0 = ((((ALP050*zt & 724 & + ALP140*zs+ALP040)*zt & 725 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 726 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 727 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 728 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 729 ! 730 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 731 ! 732 pab(ji,jj,jp_tem) = zn * r1_rau0 733 ! 734 ! beta 735 zn3 = BET003 736 ! 737 zn2 = BET012*zt + BET102*zs+BET002 738 ! 739 zn1 = ((BET031*zt & 740 & + BET121*zs+BET021)*zt & 741 & + (BET211*zs+BET111)*zs+BET011)*zt & 742 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 743 ! 744 zn0 = ((((BET050*zt & 745 & + BET140*zs+BET040)*zt & 746 & + (BET230*zs+BET130)*zs+BET030)*zt & 747 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 748 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 749 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 750 ! 751 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 752 ! 753 pab(ji,jj,jp_sal) = zn / zs * r1_rau0 754 ! 755 ! 756 END DO 757 END DO 758 ! ! Lateral boundary conditions 759 CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. ) 671 DO_2D_11_11 672 ! 673 zh = pdep(ji,jj) * r1_Z0 ! depth 674 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 675 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 676 ! 677 ! alpha 678 zn3 = ALP003 679 ! 680 zn2 = ALP012*zt + ALP102*zs+ALP002 681 ! 682 zn1 = ((ALP031*zt & 683 & + ALP121*zs+ALP021)*zt & 684 & + (ALP211*zs+ALP111)*zs+ALP011)*zt & 685 & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 686 ! 687 zn0 = ((((ALP050*zt & 688 & + ALP140*zs+ALP040)*zt & 689 & + (ALP230*zs+ALP130)*zs+ALP030)*zt & 690 & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & 691 & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & 692 & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 693 ! 694 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 695 ! 696 pab(ji,jj,jp_tem) = zn * r1_rho0 697 ! 698 ! beta 699 zn3 = BET003 700 ! 701 zn2 = BET012*zt + BET102*zs+BET002 702 ! 703 zn1 = ((BET031*zt & 704 & + BET121*zs+BET021)*zt & 705 & + (BET211*zs+BET111)*zs+BET011)*zt & 706 & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 707 ! 708 zn0 = ((((BET050*zt & 709 & + BET140*zs+BET040)*zt & 710 & + (BET230*zs+BET130)*zs+BET030)*zt & 711 & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & 712 & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & 713 & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 714 ! 715 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 716 ! 717 pab(ji,jj,jp_sal) = zn / zs * r1_rho0 718 ! 719 ! 720 END_2D 760 721 ! 761 722 CASE( np_seos ) !== simplified EOS ==! 762 723 ! 763 DO jj = 1, jpjm1 764 DO ji = 1, fs_jpim1 ! vector opt. 765 ! 766 zt = pts (ji,jj,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 767 zs = pts (ji,jj,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 768 zh = pdep (ji,jj) ! depth at the partial step level 769 ! 770 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 771 pab(ji,jj,jp_tem) = zn * r1_rau0 ! alpha 772 ! 773 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 774 pab(ji,jj,jp_sal) = zn * r1_rau0 ! beta 775 ! 776 END DO 777 END DO 778 ! ! Lateral boundary conditions 779 CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. ) 724 DO_2D_11_11 725 ! 726 zt = pts (ji,jj,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 727 zs = pts (ji,jj,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 728 zh = pdep (ji,jj) ! depth at the partial step level 729 ! 730 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 731 pab(ji,jj,jp_tem) = zn * r1_rho0 ! alpha 732 ! 733 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 734 pab(ji,jj,jp_sal) = zn * r1_rho0 ! beta 735 ! 736 END_2D 780 737 ! 781 738 CASE DEFAULT … … 785 742 END SELECT 786 743 ! 787 IF( ln_ctl) CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &788 & tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )744 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', & 745 & tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' ) 789 746 ! 790 747 IF( ln_timing ) CALL timing_stop('rab_2d') … … 793 750 794 751 795 SUBROUTINE rab_0d( pts, pdep, pab )752 SUBROUTINE rab_0d( pts, pdep, pab, Kmm ) 796 753 !!---------------------------------------------------------------------- 797 754 !! *** ROUTINE rab_0d *** … … 801 758 !! ** Action : - pab : thermal/haline expansion ratio at T-points 802 759 !!---------------------------------------------------------------------- 760 INTEGER , INTENT(in ) :: Kmm ! time level index 803 761 REAL(wp), DIMENSION(jpts) , INTENT(in ) :: pts ! pot. temperature & salinity 804 762 REAL(wp), INTENT(in ) :: pdep ! depth [m] … … 841 799 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 842 800 ! 843 pab(jp_tem) = zn * r1_r au0801 pab(jp_tem) = zn * r1_rho0 844 802 ! 845 803 ! beta … … 862 820 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 863 821 ! 864 pab(jp_sal) = zn / zs * r1_r au0822 pab(jp_sal) = zn / zs * r1_rho0 865 823 ! 866 824 ! … … 873 831 ! 874 832 zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 875 pab(jp_tem) = zn * r1_r au0 ! alpha833 pab(jp_tem) = zn * r1_rho0 ! alpha 876 834 ! 877 835 zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 878 pab(jp_sal) = zn * r1_r au0 ! beta836 pab(jp_sal) = zn * r1_rho0 ! beta 879 837 ! 880 838 CASE DEFAULT … … 889 847 890 848 891 SUBROUTINE bn2( pts, pab, pn2 )849 SUBROUTINE bn2( pts, pab, pn2, Kmm ) 892 850 !!---------------------------------------------------------------------- 893 851 !! *** ROUTINE bn2 *** … … 903 861 !! 904 862 !!---------------------------------------------------------------------- 863 INTEGER , INTENT(in ) :: Kmm ! time level index 905 864 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu] 906 865 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pab ! thermal/haline expansion coef. [Celsius-1,psu-1] … … 913 872 IF( ln_timing ) CALL timing_start('bn2') 914 873 ! 915 DO jk = 2, jpkm1 ! interior points only (2=< jk =< jpkm1 ) 916 DO jj = 1, jpj ! surface and bottom value set to zero one for all in istate.F90 917 DO ji = 1, jpi 918 zrw = ( gdepw_n(ji,jj,jk ) - gdept_n(ji,jj,jk) ) & 919 & / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 920 ! 921 zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 922 zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 923 ! 924 pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 925 & - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) & 926 & / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 927 END DO 928 END DO 929 END DO 930 ! 931 IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', kdim=jpk ) 874 DO_3D_11_11( 2, jpkm1 ) 875 zrw = ( gdepw(ji,jj,jk ,Kmm) - gdept(ji,jj,jk,Kmm) ) & 876 & / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 877 ! 878 zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 879 zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 880 ! 881 pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 882 & - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) & 883 & / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 884 END_3D 885 ! 886 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', kdim=jpk ) 932 887 ! 933 888 IF( ln_timing ) CALL timing_stop('bn2') … … 965 920 z1_T0 = 1._wp/40._wp 966 921 ! 967 DO jj = 1, jpj 968 DO ji = 1, jpi 969 ! 970 zt = ctmp (ji,jj) * z1_T0 971 zs = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 ) 972 ztm = tmask(ji,jj,1) 973 ! 974 zn = ((((-2.1385727895e-01_wp*zt & 975 & - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt & 976 & + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt & 977 & + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt & 978 & + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs & 979 & +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt & 980 & + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs & 981 & -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp 982 ! 983 zd = (2.0035003456_wp*zt & 984 & -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt & 985 & + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp 986 ! 987 ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm 988 ! 989 END DO 990 END DO 922 DO_2D_11_11 923 ! 924 zt = ctmp (ji,jj) * z1_T0 925 zs = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 ) 926 ztm = tmask(ji,jj,1) 927 ! 928 zn = ((((-2.1385727895e-01_wp*zt & 929 & - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt & 930 & + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt & 931 & + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt & 932 & + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs & 933 & +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt & 934 & + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs & 935 & -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp 936 ! 937 zd = (2.0035003456_wp*zt & 938 & -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt & 939 & + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp 940 ! 941 ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm 942 ! 943 END_2D 991 944 ! 992 945 IF( ln_timing ) CALL timing_stop('eos_pt_from_ct') … … 1020 973 ! 1021 974 z1_S0 = 1._wp / 35.16504_wp 1022 DO jj = 1, jpj 1023 DO ji = 1, jpi 1024 zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 ) ! square root salinity 1025 ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 1026 & - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 1027 END DO 1028 END DO 975 DO_2D_11_11 976 zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 ) ! square root salinity 977 ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 978 & - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 979 END_2D 1029 980 ptf(:,:) = ptf(:,:) * psal(:,:) 1030 981 ! … … 1093 1044 1094 1045 1095 SUBROUTINE eos_pen( pts, pab_pe, ppen )1046 SUBROUTINE eos_pen( pts, pab_pe, ppen, Kmm ) 1096 1047 !!---------------------------------------------------------------------- 1097 1048 !! *** ROUTINE eos_pen *** … … 1101 1052 !! ** Method : PE is defined analytically as the vertical 1102 1053 !! primitive of EOS times -g integrated between 0 and z>0. 1103 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - r au0 gz ) / rau0 gz - rd1054 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd 1104 1055 !! = 1/z * /int_0^z rd dz - rd 1105 1056 !! where rd is the density anomaly (see eos_rhd function) 1106 1057 !! ab_pe are partial derivatives of PE anomaly with respect to T and S: 1107 !! ab_pe(1) = - 1/(r au0 gz) * dPE/dT + drd/dT = - d(pen)/dT1108 !! ab_pe(2) = 1/(r au0 gz) * dPE/dS + drd/dS = d(pen)/dS1058 !! ab_pe(1) = - 1/(rho0 gz) * dPE/dT + drd/dT = - d(pen)/dT 1059 !! ab_pe(2) = 1/(rho0 gz) * dPE/dS + drd/dS = d(pen)/dS 1109 1060 !! 1110 1061 !! ** Action : - pen : PE anomaly given at T-points … … 1113 1064 !! pab_pe(:,:,:,jp_sal) is beta_pe 1114 1065 !!---------------------------------------------------------------------- 1066 INTEGER , INTENT(in ) :: Kmm ! time level index 1115 1067 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 1116 1068 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab_pe ! alpha_pe and beta_pe … … 1128 1080 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 1129 1081 ! 1130 DO jk = 1, jpkm1 1131 DO jj = 1, jpj 1132 DO ji = 1, jpi 1133 ! 1134 zh = gdept_n(ji,jj,jk) * r1_Z0 ! depth 1135 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 1136 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 1137 ztm = tmask(ji,jj,jk) ! tmask 1138 ! 1139 ! potential energy non-linear anomaly 1140 zn2 = (PEN012)*zt & 1141 & + PEN102*zs+PEN002 1142 ! 1143 zn1 = ((PEN021)*zt & 1144 & + PEN111*zs+PEN011)*zt & 1145 & + (PEN201*zs+PEN101)*zs+PEN001 1146 ! 1147 zn0 = ((((PEN040)*zt & 1148 & + PEN130*zs+PEN030)*zt & 1149 & + (PEN220*zs+PEN120)*zs+PEN020)*zt & 1150 & + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt & 1151 & + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000 1152 ! 1153 zn = ( zn2 * zh + zn1 ) * zh + zn0 1154 ! 1155 ppen(ji,jj,jk) = zn * zh * r1_rau0 * ztm 1156 ! 1157 ! alphaPE non-linear anomaly 1158 zn2 = APE002 1159 ! 1160 zn1 = (APE011)*zt & 1161 & + APE101*zs+APE001 1162 ! 1163 zn0 = (((APE030)*zt & 1164 & + APE120*zs+APE020)*zt & 1165 & + (APE210*zs+APE110)*zs+APE010)*zt & 1166 & + ((APE300*zs+APE200)*zs+APE100)*zs+APE000 1167 ! 1168 zn = ( zn2 * zh + zn1 ) * zh + zn0 1169 ! 1170 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm 1171 ! 1172 ! betaPE non-linear anomaly 1173 zn2 = BPE002 1174 ! 1175 zn1 = (BPE011)*zt & 1176 & + BPE101*zs+BPE001 1177 ! 1178 zn0 = (((BPE030)*zt & 1179 & + BPE120*zs+BPE020)*zt & 1180 & + (BPE210*zs+BPE110)*zs+BPE010)*zt & 1181 & + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000 1182 ! 1183 zn = ( zn2 * zh + zn1 ) * zh + zn0 1184 ! 1185 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm 1186 ! 1187 END DO 1188 END DO 1189 END DO 1082 DO_3D_11_11( 1, jpkm1 ) 1083 ! 1084 zh = gdept(ji,jj,jk,Kmm) * r1_Z0 ! depth 1085 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 1086 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 1087 ztm = tmask(ji,jj,jk) ! tmask 1088 ! 1089 ! potential energy non-linear anomaly 1090 zn2 = (PEN012)*zt & 1091 & + PEN102*zs+PEN002 1092 ! 1093 zn1 = ((PEN021)*zt & 1094 & + PEN111*zs+PEN011)*zt & 1095 & + (PEN201*zs+PEN101)*zs+PEN001 1096 ! 1097 zn0 = ((((PEN040)*zt & 1098 & + PEN130*zs+PEN030)*zt & 1099 & + (PEN220*zs+PEN120)*zs+PEN020)*zt & 1100 & + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt & 1101 & + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000 1102 ! 1103 zn = ( zn2 * zh + zn1 ) * zh + zn0 1104 ! 1105 ppen(ji,jj,jk) = zn * zh * r1_rho0 * ztm 1106 ! 1107 ! alphaPE non-linear anomaly 1108 zn2 = APE002 1109 ! 1110 zn1 = (APE011)*zt & 1111 & + APE101*zs+APE001 1112 ! 1113 zn0 = (((APE030)*zt & 1114 & + APE120*zs+APE020)*zt & 1115 & + (APE210*zs+APE110)*zs+APE010)*zt & 1116 & + ((APE300*zs+APE200)*zs+APE100)*zs+APE000 1117 ! 1118 zn = ( zn2 * zh + zn1 ) * zh + zn0 1119 ! 1120 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 1121 ! 1122 ! betaPE non-linear anomaly 1123 zn2 = BPE002 1124 ! 1125 zn1 = (BPE011)*zt & 1126 & + BPE101*zs+BPE001 1127 ! 1128 zn0 = (((BPE030)*zt & 1129 & + BPE120*zs+BPE020)*zt & 1130 & + (BPE210*zs+BPE110)*zs+BPE010)*zt & 1131 & + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000 1132 ! 1133 zn = ( zn2 * zh + zn1 ) * zh + zn0 1134 ! 1135 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 1136 ! 1137 END_3D 1190 1138 ! 1191 1139 CASE( np_seos ) !== Vallis (2006) simplified EOS ==! 1192 1140 ! 1193 DO jk = 1, jpkm1 1194 DO jj = 1, jpj 1195 DO ji = 1, jpi 1196 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! temperature anomaly (t-T0) 1197 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 1198 zh = gdept_n(ji,jj,jk) ! depth in meters at t-point 1199 ztm = tmask(ji,jj,jk) ! tmask 1200 zn = 0.5_wp * zh * r1_rau0 * ztm 1201 ! ! Potential Energy 1202 ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 1203 ! ! alphaPE 1204 pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn 1205 pab_pe(ji,jj,jk,jp_sal) = rn_b0 * rn_mu2 * zn 1206 ! 1207 END DO 1208 END DO 1209 END DO 1141 DO_3D_11_11( 1, jpkm1 ) 1142 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! temperature anomaly (t-T0) 1143 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 1144 zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point 1145 ztm = tmask(ji,jj,jk) ! tmask 1146 zn = 0.5_wp * zh * r1_rho0 * ztm 1147 ! ! Potential Energy 1148 ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 1149 ! ! alphaPE 1150 pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn 1151 pab_pe(ji,jj,jk,jp_sal) = rn_b0 * rn_mu2 * zn 1152 ! 1153 END_3D 1210 1154 ! 1211 1155 CASE DEFAULT … … 1235 1179 !!---------------------------------------------------------------------- 1236 1180 ! 1237 REWIND( numnam_ref ) ! Namelist nameos in reference namelist : equation of state1238 1181 READ ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 ) 1239 1182 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist' ) 1240 1183 ! 1241 REWIND( numnam_cfg ) ! Namelist nameos in configuration namelist : equation of state1242 1184 READ ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 ) 1243 1185 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist' ) 1244 1186 IF(lwm) WRITE( numond, nameos ) 1245 1187 ! 1246 r au0 = 1026._wp !: volumic mass of reference [kg/m3]1188 rho0 = 1026._wp !: volumic mass of reference [kg/m3] 1247 1189 rcp = 3991.86795711963_wp !: heat capacity [J/K] 1248 1190 ! … … 1656 1598 WRITE(numout,*) ' ==>>> use of simplified eos: ' 1657 1599 WRITE(numout,*) ' rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT ' 1658 WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / r au0'1600 WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0' 1659 1601 WRITE(numout,*) ' with the following coefficients :' 1660 1602 WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0 … … 1675 1617 END SELECT 1676 1618 ! 1677 r au0_rcp = rau0 * rcp1678 r1_r au0 = 1._wp / rau01619 rho0_rcp = rho0 * rcp 1620 r1_rho0 = 1._wp / rho0 1679 1621 r1_rcp = 1._wp / rcp 1680 r1_r au0_rcp = 1._wp / rau0_rcp1622 r1_rho0_rcp = 1._wp / rho0_rcp 1681 1623 ! 1682 1624 IF(lwp) THEN … … 1693 1635 IF(lwp) WRITE(numout,*) 1694 1636 IF(lwp) WRITE(numout,*) ' Associated physical constant' 1695 IF(lwp) WRITE(numout,*) ' volumic mass of reference r au0 = ', rau0 , ' kg/m^3'1696 IF(lwp) WRITE(numout,*) ' 1. / r au0 r1_rau0 = ', r1_rau0, ' m^3/kg'1637 IF(lwp) WRITE(numout,*) ' volumic mass of reference rho0 = ', rho0 , ' kg/m^3' 1638 IF(lwp) WRITE(numout,*) ' 1. / rho0 r1_rho0 = ', r1_rho0, ' m^3/kg' 1697 1639 IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin' 1698 IF(lwp) WRITE(numout,*) ' r au0 * rcp rau0_rcp = ', rau0_rcp1699 IF(lwp) WRITE(numout,*) ' 1. / ( r au0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp1640 IF(lwp) WRITE(numout,*) ' rho0 * rcp rho0_rcp = ', rho0_rcp 1641 IF(lwp) WRITE(numout,*) ' 1. / ( rho0 * rcp ) r1_rho0_rcp = ', r1_rho0_rcp 1700 1642 ! 1701 1643 END SUBROUTINE eos_init
Note: See TracChangeset
for help on using the changeset viewer.