- Timestamp:
- 2015-05-27T12:47:55+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5177_CNRS4_stopar/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r5147 r5296 30 30 !! eos_insitu_2d : Compute the in situ density for 2d fields 31 31 !! bn2 : Compute the Brunt-Vaisala frequency 32 !! eos_rab : generic interface of in situ thermal/haline expansion ratio 32 !! eos_rab : generic interface of in situ thermal/haline expansion ratio 33 33 !! eos_rab_3d : compute in situ thermal/haline expansion ratio 34 34 !! eos_rab_2d : compute in situ thermal/haline expansion ratio for 2d fields … … 42 42 USE in_out_manager ! I/O manager 43 43 USE lib_mpp ! MPP library 44 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 44 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 45 45 USE prtctl ! Print control 46 46 USE wrk_nemo ! Memory Allocation 47 USE lbclnk ! ocean lateral boundary conditions47 USE lbclnk ! ocean lateral boundary conditions 48 48 USE timing ! Timing 49 USE stopar ! Stochastic T/S fluctuations 50 USE stopts ! Stochastic T/S fluctuations 49 51 50 52 IMPLICIT NONE … … 60 62 END INTERFACE 61 63 ! 62 INTERFACE eos_fzp 64 INTERFACE eos_fzp 63 65 MODULE PROCEDURE eos_fzp_2d, eos_fzp_0d 64 66 END INTERFACE … … 78 80 ! !!! simplified eos coefficients 79 81 ! default value: Vallis 2006 80 REAL(wp) :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff. 81 REAL(wp) :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff. 82 REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2 83 REAL(wp) :: rn_lambda2 = 5.4914e-4_wp ! cabbeling coeff. in S^2 84 REAL(wp) :: rn_mu1 = 1.4970e-4_wp ! thermobaric coeff. in T 85 REAL(wp) :: rn_mu2 = 1.1090e-5_wp ! thermobaric coeff. in S 86 REAL(wp) :: rn_nu = 2.4341e-3_wp ! cabbeling coeff. in theta*salt 87 82 REAL(wp) :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff. 83 REAL(wp) :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff. 84 REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2 85 REAL(wp) :: rn_lambda2 = 5.4914e-4_wp ! cabbeling coeff. in S^2 86 REAL(wp) :: rn_mu1 = 1.4970e-4_wp ! thermobaric coeff. in T 87 REAL(wp) :: rn_mu2 = 1.1090e-5_wp ! thermobaric coeff. in S 88 REAL(wp) :: rn_nu = 2.4341e-3_wp ! cabbeling coeff. in theta*salt 89 88 90 ! TEOS10/EOS80 parameters 89 91 REAL(wp) :: r1_S0, r1_T0, r1_Z0, rdeltaS 90 92 91 93 ! EOS parameters 92 94 REAL(wp) :: EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600 … … 106 108 REAL(wp) :: EOS022 107 109 REAL(wp) :: EOS003 , EOS103 108 REAL(wp) :: EOS013 109 110 REAL(wp) :: EOS013 111 110 112 ! ALPHA parameters 111 113 REAL(wp) :: ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500 … … 122 124 REAL(wp) :: ALP012 123 125 REAL(wp) :: ALP003 124 126 125 127 ! BETA parameters 126 128 REAL(wp) :: BET000 , BET100 , BET200 , BET300 , BET400 , BET500 … … 149 151 REAL(wp) :: PEN002 , PEN102 150 152 REAL(wp) :: PEN012 151 153 152 154 ! ALPHA_PEN parameters 153 155 REAL(wp) :: APE000 , APE100 , APE200 , APE300 … … 279 281 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 280 282 & - rn_nu * zt * zs 281 ! 283 ! 282 284 prd(ji,jj,jk) = zn * r1_rau0 * ztm ! density anomaly (masked) 283 285 END DO … … 313 315 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m] 314 316 ! 315 INTEGER :: ji, jj, jk ! dummy loop indices 316 REAL(wp) :: zt , zh , zs , ztm ! local scalars 317 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 317 INTEGER :: ji, jj, jk, jsmp, jdof ! dummy loop indices 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 318 321 !!---------------------------------------------------------------------- 319 322 ! … … 324 327 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 325 328 ! 326 DO jk = 1, jpkm1 327 DO jj = 1, jpj 328 DO ji = 1, jpi 329 ! 330 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 331 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 332 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 333 ztm = tmask(ji,jj,jk) ! tmask 334 ! 335 zn3 = EOS013*zt & 336 & + EOS103*zs+EOS003 337 ! 338 zn2 = (EOS022*zt & 339 & + EOS112*zs+EOS012)*zt & 340 & + (EOS202*zs+EOS102)*zs+EOS002 341 ! 342 zn1 = (((EOS041*zt & 343 & + EOS131*zs+EOS031)*zt & 344 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 345 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 346 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 347 ! 348 zn0 = (((((EOS060*zt & 349 & + EOS150*zs+EOS050)*zt & 350 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 351 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 352 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 353 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 354 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 355 ! 356 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 357 ! 358 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 359 ! 360 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 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 jk = 1, jpkm1 340 DO jj = 1, jpj 341 DO ji = 1, jpi 342 ! 343 ! compute density (2*nn_sto_eos) times: 344 ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 345 ! (2) for t-dt, s-ds (with the opposite fluctuation) 346 DO jsmp = 1, nn_sto_eos*2 347 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 348 zt = pts (ji,jj,jk,jp_tem) * r1_T0 + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp) ! temperature 349 zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 350 zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity 351 ztm = tmask(ji,jj,jk) ! tmask 352 ! 353 zn3 = EOS013*zt & 354 & + EOS103*zs+EOS003 355 ! 356 zn2 = (EOS022*zt & 357 & + EOS112*zs+EOS012)*zt & 358 & + (EOS202*zs+EOS102)*zs+EOS002 359 ! 360 zn1 = (((EOS041*zt & 361 & + EOS131*zs+EOS031)*zt & 362 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 363 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 364 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 365 ! 366 zn0_sto(jsmp) = (((((EOS060*zt & 367 & + EOS150*zs+EOS050)*zt & 368 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 369 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 370 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 371 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 372 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 373 ! 374 zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 375 END DO 376 ! 377 ! 378 ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 379 prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 380 DO jsmp = 1, nn_sto_eos*2 381 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 382 ! 383 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rau0 - 1._wp ) ! density anomaly (masked) 384 END DO 385 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 386 prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos 387 END DO 361 388 END DO 362 389 END DO 363 END DO 390 DEALLOCATE(zn0_sto,zn_sto,zsign) 391 ! Non-stochastic equation of state 392 ELSE 393 DO jk = 1, jpkm1 394 DO jj = 1, jpj 395 DO ji = 1, jpi 396 ! 397 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 398 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 399 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 400 ztm = tmask(ji,jj,jk) ! tmask 401 ! 402 zn3 = EOS013*zt & 403 & + EOS103*zs+EOS003 404 ! 405 zn2 = (EOS022*zt & 406 & + EOS112*zs+EOS012)*zt & 407 & + (EOS202*zs+EOS102)*zs+EOS002 408 ! 409 zn1 = (((EOS041*zt & 410 & + EOS131*zs+EOS031)*zt & 411 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 412 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 413 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 414 ! 415 zn0 = (((((EOS060*zt & 416 & + EOS150*zs+EOS050)*zt & 417 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 418 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 419 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 420 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 421 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 422 ! 423 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 424 ! 425 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 426 ! 427 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 428 END DO 429 END DO 430 END DO 431 ENDIF 364 432 ! 365 433 CASE( 1 ) !== simplified EOS ==! … … 681 749 ! 682 750 CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. ) ! Lateral boundary conditions 683 CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. ) 751 CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. ) 684 752 ! 685 753 CASE( 1 ) !== simplified EOS ==! … … 702 770 ! 703 771 CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. ) ! Lateral boundary conditions 704 CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. ) 772 CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. ) 705 773 ! 706 774 CASE DEFAULT … … 820 888 !! *** ROUTINE bn2 *** 821 889 !! 822 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the 890 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the 823 891 !! time-step of the input arguments 824 892 !! … … 827 895 !! N.B. N^2 is set one for all to zero at jk=1 in istate module. 828 896 !! 829 !! ** Action : pn2 : square of the brunt-vaisala frequency at w-point 897 !! ** Action : pn2 : square of the brunt-vaisala frequency at w-point 830 898 !! 831 899 !!---------------------------------------------------------------------- … … 844 912 DO ji = 1, jpi 845 913 zrw = ( fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk) ) & 846 & / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 847 ! 848 zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 914 & / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 915 ! 916 zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 849 917 zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 850 918 ! … … 1026 1094 !! ** Purpose : Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points 1027 1095 !! 1028 !! ** Method : PE is defined analytically as the vertical 1096 !! ** Method : PE is defined analytically as the vertical 1029 1097 !! primitive of EOS times -g integrated between 0 and z>0. 1030 1098 !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd 1031 !! = 1/z * /int_0^z rd dz - rd 1099 !! = 1/z * /int_0^z rd dz - rd 1032 1100 !! where rd is the density anomaly (see eos_rhd function) 1033 1101 !! ab_pe are partial derivatives of PE anomaly with respect to T and S: … … 1094 1162 ! 1095 1163 zn = ( zn2 * zh + zn1 ) * zh + zn0 1096 ! 1164 ! 1097 1165 pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm 1098 1166 ! … … 1109 1177 ! 1110 1178 zn = ( zn2 * zh + zn1 ) * zh + zn0 1111 ! 1179 ! 1112 1180 pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm 1113 1181 ! … … 1589 1657 END SELECT 1590 1658 ! 1591 rau0_rcp = rau0 * rcp 1659 rau0_rcp = rau0 * rcp 1592 1660 r1_rau0 = 1._wp / rau0 1593 1661 r1_rcp = 1._wp / rcp 1594 r1_rau0_rcp = 1._wp / rau0_rcp 1662 r1_rau0_rcp = 1._wp / rau0_rcp 1595 1663 ! 1596 1664 IF(lwp) WRITE(numout,*)
Note: See TracChangeset
for help on using the changeset viewer.
