Changeset 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
- Timestamp:
- 2012-01-28T17:44:18+01:00 (12 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 (modified) (29 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r2715 r3294 3 3 !! *** MODULE eosbn2 *** 4 4 !! Ocean diagnostic variable : equation of state - in situ and potential density 5 !! - Brunt-Vaisala frequency 5 !! - Brunt-Vaisala frequency 6 6 !!============================================================================== 7 7 !! History : OPA ! 1989-03 (O. Marti) Original code … … 27 27 !! eos_insitu_2d : Compute the in situ density for 2d fields 28 28 !! eos_bn2 : Compute the Brunt-Vaisala frequency 29 !! eos_alpbet : calculates the in situ thermal and haline expansion coeff.29 !! eos_alpbet : calculates the in situ thermal/haline expansion ratio 30 30 !! tfreez : Compute the surface freezing temperature 31 31 !! eos_init : set eos parameters (namelist) … … 37 37 USE lib_mpp ! MPP library 38 38 USE prtctl ! Print control 39 USE wrk_nemo ! Memory Allocation 40 USE timing ! Timing 39 41 40 42 IMPLICIT NONE 41 43 PRIVATE 42 44 43 ! !! * Interface 45 ! !! * Interface 44 46 INTERFACE eos 45 47 MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 46 END INTERFACE 48 END INTERFACE 47 49 INTERFACE bn2 48 50 MODULE PROCEDURE eos_bn2 49 END INTERFACE 51 END INTERFACE 50 52 51 53 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules … … 61 63 62 64 REAL(wp), PUBLIC :: ralpbet !: alpha / beta ratio 63 65 64 66 !! * Substitutions 65 67 # include "domzgr_substitute.h90" … … 75 77 !!---------------------------------------------------------------------- 76 78 !! *** ROUTINE eos_insitu *** 77 !! 78 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 79 !! 80 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 79 81 !! potential temperature and salinity using an equation of state 80 82 !! defined through the namelist parameter nn_eos. … … 108 110 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 109 111 !!---------------------------------------------------------------------- 110 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released111 USE wrk_nemo, ONLY: zws => wrk_3d_1 ! 3D workspace112 112 !! 113 113 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] … … 122 122 REAL(wp) :: zb1, za1, zkw, zk0 ! - - 123 123 REAL(wp) :: zrau0r ! - - 124 !!---------------------------------------------------------------------- 125 126 IF( wrk_in_use(3, 1) ) THEN 127 CALL ctl_stop('eos_insitu: requested workspace array unavailable') ; RETURN 128 ENDIF 129 124 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 125 !!---------------------------------------------------------------------- 126 127 ! 128 IF( nn_timing == 1 ) CALL timing_start('eos') 129 ! 130 CALL wrk_alloc( jpi, jpj, jpk, zws ) 131 ! 130 132 SELECT CASE( nn_eos ) 131 133 ! … … 134 136 !CDIR NOVERRCHK 135 137 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 136 ! 138 ! 137 139 DO jk = 1, jpkm1 138 140 DO jj = 1, jpj … … 191 193 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos : ', ovlap=1, kdim=jpk ) 192 194 ! 193 IF( wrk_not_released(3, 1) ) CALL ctl_stop('eos_insitu: failed to release workspace array') 195 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 196 ! 197 IF( nn_timing == 1 ) CALL timing_stop('eos') 194 198 ! 195 199 END SUBROUTINE eos_insitu … … 199 203 !!---------------------------------------------------------------------- 200 204 !! *** ROUTINE eos_insitu_pot *** 201 !! 205 !! 202 206 !! ** Purpose : Compute the in situ density (ratio rho/rau0) and the 203 207 !! potential volumic mass (Kg/m3) from potential temperature and 204 !! salinity fields using an equation of state defined through the 208 !! salinity fields using an equation of state defined through the 205 209 !! namelist parameter nn_eos. 206 210 !! … … 230 234 !! nn_eos = 2 : linear equation of state function of temperature and 231 235 !! salinity 232 !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 236 !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 233 237 !! = rn_beta * s - rn_alpha * tn - 1. 234 238 !! rhop(t,s) = rho(t,s) … … 242 246 !! Brown and Campana, Mon. Weather Rev., 1978 243 247 !!---------------------------------------------------------------------- 244 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released245 USE wrk_nemo, ONLY: zws => wrk_3d_1 ! 3D workspace246 248 !! 247 249 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] … … 253 255 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! local scalars 254 256 REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zrau0r ! - - 255 !!---------------------------------------------------------------------- 256 257 IF( wrk_in_use(3, 1) ) THEN 258 CALL ctl_stop('eos_insitu_pot: requested workspace array unavailable') ; RETURN 259 ENDIF 260 257 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 258 !!---------------------------------------------------------------------- 259 ! 260 IF( nn_timing == 1 ) CALL timing_start('eos-p') 261 ! 262 CALL wrk_alloc( jpi, jpj, jpk, zws ) 263 ! 261 264 SELECT CASE ( nn_eos ) 262 265 ! … … 265 268 !CDIR NOVERRCHK 266 269 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 267 ! 270 ! 268 271 DO jk = 1, jpkm1 269 272 DO jj = 1, jpj … … 327 330 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk ) 328 331 ! 329 IF( wrk_not_released(3, 1) ) CALL ctl_stop('eos_insitu_pot: failed to release workspace array') 332 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 333 ! 334 IF( nn_timing == 1 ) CALL timing_stop('eos-p') 330 335 ! 331 336 END SUBROUTINE eos_insitu_pot … … 336 341 !! *** ROUTINE eos_insitu_2d *** 337 342 !! 338 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 343 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 339 344 !! potential temperature and salinity using an equation of state 340 345 !! defined through the namelist parameter nn_eos. * 2D field case … … 368 373 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 369 374 !!---------------------------------------------------------------------- 370 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released371 USE wrk_nemo, ONLY: zws => wrk_2d_5 ! 2D workspace372 375 !! 373 376 REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 374 377 ! ! 2 : salinity [psu] 375 378 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 376 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 379 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 377 380 !! 378 381 INTEGER :: ji, jj ! dummy loop indices 379 382 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! temporary scalars 380 383 REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zmask ! - - 381 !!---------------------------------------------------------------------- 382 383 IF( wrk_in_use(2, 5) ) THEN 384 CALL ctl_stop('eos_insitu_2d: requested workspace array unavailable') ; RETURN 385 ENDIF 384 REAL(wp), POINTER, DIMENSION(:,:) :: zws 385 !!---------------------------------------------------------------------- 386 ! 387 IF( nn_timing == 1 ) CALL timing_start('eos2d') 388 ! 389 CALL wrk_alloc( jpi, jpj, zws ) 390 ! 386 391 387 392 prd(:,:) = 0._wp … … 449 454 DO jj = 1, jpjm1 450 455 DO ji = 1, fs_jpim1 ! vector opt. 451 prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 456 prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 452 457 END DO 453 458 END DO … … 457 462 IF(ln_ctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 458 463 ! 459 IF( wrk_not_released(2, 5) ) CALL ctl_stop('eos_insitu_2d: failed to release workspace array') 464 CALL wrk_dealloc( jpi, jpj, zws ) 465 ! 466 IF( nn_timing == 1 ) CALL timing_stop('eos2d') 460 467 ! 461 468 END SUBROUTINE eos_insitu_2d … … 468 475 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the time- 469 476 !! step of the input arguments 470 !! 477 !! 471 478 !! ** Method : 472 479 !! * nn_eos = 0 : UNESCO sea water properties … … 482 489 !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 483 490 !! The use of potential density to compute N^2 introduces e r r o r 484 !! in the sign of N^2 at great depths. We recommand the use of 491 !! in the sign of N^2 at great depths. We recommand the use of 485 492 !! nn_eos = 0, except for academical studies. 486 493 !! Macro-tasked on horizontal slab (jk-loop) … … 497 504 !! 498 505 INTEGER :: ji, jj, jk ! dummy loop indices 499 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 506 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 500 507 #if defined key_zdfddm 501 508 REAL(wp) :: zds ! local scalars … … 503 510 !!---------------------------------------------------------------------- 504 511 512 ! 513 IF( nn_timing == 1 ) CALL timing_start('bn2') 514 ! 505 515 ! pn2 : interior points only (2=< jk =< jpkm1 ) 506 ! -------------------------- 516 ! -------------------------- 507 517 ! 508 518 SELECT CASE( nn_eos ) … … 542 552 & - 0.121555e-07_wp ) * zh 543 553 ! 544 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 554 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 545 555 & * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 546 556 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) … … 565 575 & - rn_beta * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) ) ) & 566 576 & / fse3w(:,:,jk) * tmask(:,:,jk) 567 END DO 577 END DO 568 578 #if defined key_zdfddm 569 579 DO jk = 2, jpkm1 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 570 580 DO jj = 1, jpj 571 581 DO ji = 1, jpi 572 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 582 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 573 583 IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 574 584 rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds … … 584 594 #endif 585 595 ! 596 IF( nn_timing == 1 ) CALL timing_stop('bn2') 597 ! 586 598 END SUBROUTINE eos_bn2 587 599 588 600 589 SUBROUTINE eos_alpbet( pts, palp h, pbeta)590 !!---------------------------------------------------------------------- 591 !! *** ROUTINE ldf_slp_grif***592 !! 593 !! ** Purpose : Calculates the thermal and haline expansion coefficients at T-points.594 !! 595 !! ** Method : calculates alpha and beta at T-points601 SUBROUTINE eos_alpbet( pts, palpbet, beta0 ) 602 !!---------------------------------------------------------------------- 603 !! *** ROUTINE eos_alpbet *** 604 !! 605 !! ** Purpose : Calculates the in situ thermal/haline expansion ratio at T-points 606 !! 607 !! ** Method : calculates alpha / beta ratio at T-points 596 608 !! * nn_eos = 0 : UNESCO sea water properties 597 !! The brunt-vaisala frequency is computed using the polynomial 598 !! polynomial expression of McDougall (1987): 599 !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 600 !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 601 !! computed and used in zdfddm module : 602 !! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 609 !! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial 610 !! polynomial expression of McDougall (1987). 611 !! Scalar beta0 is returned = 1. 603 612 !! * nn_eos = 1 : linear equation of state (temperature only) 604 !! N^2 = grav * rn_alpha * dk[ t ]/e3w 613 !! The ratio is undefined, so we return alpha as palpbet 614 !! Scalar beta0 is returned = 0. 605 615 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 606 !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 607 !! * nn_eos = 3 : Jackett JAOT 2003 ??? 608 !! 609 !! ** Action : - palph, pbeta : thermal and haline expansion coeff. at T-point 610 !!---------------------------------------------------------------------- 611 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 612 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palph, pbeta ! thermal & haline expansion coeff. 613 ! 616 !! The alpha/beta ratio is returned as ralpbet 617 !! Scalar beta0 is returned = 1. 618 !! 619 !! ** Action : - palpbet : thermal/haline expansion ratio at T-points 620 !! : beta0 : 1. or 0. 621 !!---------------------------------------------------------------------- 622 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 623 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpbet ! thermal/haline expansion ratio 624 REAL(wp), INTENT( out) :: beta0 ! set = 1 except with case 1 eos, rho=rho(T) 625 !! 614 626 INTEGER :: ji, jj, jk ! dummy loop indices 615 REAL(wp) :: zt, zs, zh ! local scalars 616 !!---------------------------------------------------------------------- 627 REAL(wp) :: zt, zs, zh ! local scalars 628 !!---------------------------------------------------------------------- 629 ! 630 IF( nn_timing == 1 ) CALL timing_start('eos_alpbet') 617 631 ! 618 632 SELECT CASE ( nn_eos ) … … 624 638 zt = pts(ji,jj,jk,jp_tem) ! potential temperature 625 639 zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35) 626 zh = fsdept(ji,jj,jk) ! depth in meters 627 ! 628 pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & 629 & - 0.301985e-05_wp ) * zt & 630 & + 0.785567e-03_wp & 631 & + ( 0.515032e-08_wp * zs & 632 & + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs & 633 & + ( ( 0.121551e-17_wp * zh & 634 & - 0.602281e-15_wp * zs & 635 & - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh & 636 & + 0.408195e-10_wp * zs & 637 & + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt & 638 & - 0.121555e-07_wp ) * zh 639 ! 640 palph(ji,jj,jk) = - pbeta(ji,jj,jk) * & 641 & ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & 642 & - 0.203814e-03_wp ) * zt & 643 & + 0.170907e-01_wp ) * zt & 644 & + 0.665157e-01_wp & 645 & + ( - 0.678662e-05_wp * zs & 646 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs & 647 & + ( ( - 0.302285e-13_wp * zh & 648 & - 0.251520e-11_wp * zs & 649 & + 0.512857e-12_wp * zt * zt ) * zh & 650 & - 0.164759e-06_wp * zs & 651 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt & 652 & + 0.380374e-04_wp ) * zh) 640 zh = fsdept(ji,jj,jk) ! depth in meters 641 ! 642 palpbet(ji,jj,jk) = & 643 & ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & 644 & - 0.203814e-03_wp ) * zt & 645 & + 0.170907e-01_wp ) * zt & 646 & + 0.665157e-01_wp & 647 & + ( - 0.678662e-05_wp * zs & 648 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs & 649 & + ( ( - 0.302285e-13_wp * zh & 650 & - 0.251520e-11_wp * zs & 651 & + 0.512857e-12_wp * zt * zt ) * zh & 652 & - 0.164759e-06_wp * zs & 653 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt & 654 & + 0.380374e-04_wp ) * zh 653 655 END DO 654 656 END DO 655 657 END DO 656 ! 657 CASE ( 1 ) 658 palph(:,:,:) = - rn_alpha 659 pbeta(:,:,:) = 0._wp 660 ! 661 CASE ( 2 ) 662 palph(:,:,:) = - rn_alpha 663 pbeta(:,:,:) = rn_beta 658 beta0 = 1._wp 659 ! 660 CASE ( 1 ) !== Linear formulation = F( temperature ) ==! 661 palpbet(:,:,:) = rn_alpha 662 beta0 = 0._wp 663 ! 664 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 665 palpbet(:,:,:) = ralpbet 666 beta0 = 1._wp 664 667 ! 665 668 CASE DEFAULT … … 669 672 ! 670 673 END SELECT 674 ! 675 IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet') 671 676 ! 672 677 END SUBROUTINE eos_alpbet … … 747 752 748 753 !!====================================================================== 749 END MODULE eosbn2 754 END MODULE eosbn2
Note: See TracChangeset
for help on using the changeset viewer.
