- Timestamp:
- 2013-04-19T08:48:21+02:00 (11 years ago)
- Location:
- branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA
- Files:
-
- 16 edited
-
eosbn2.F90 (modified) (40 diffs)
-
traadv_cen2.F90 (modified) (9 diffs)
-
traadv_muscl.F90 (modified) (6 diffs)
-
traadv_muscl2.F90 (modified) (3 diffs)
-
traadv_qck.F90 (modified) (4 diffs)
-
traadv_tvd.F90 (modified) (2 diffs)
-
traadv_ubs.F90 (modified) (4 diffs)
-
trabbc.F90 (modified) (2 diffs)
-
trabbl.F90 (modified) (6 diffs)
-
tradmp.F90 (modified) (4 diffs)
-
traldf.F90 (modified) (3 diffs)
-
tranpc.F90 (modified) (3 diffs)
-
tranxt.F90 (modified) (4 diffs)
-
traqsr.F90 (modified) (8 diffs)
-
trasbc.F90 (modified) (7 diffs)
-
trazdf.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r3625 r3876 18 18 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 19 19 !! - ! 2010-10 (G. Nurser, G. Madec) add eos_alpbet used in ldfslp 20 !! 3.5 ! 2012-03 (F. Roquet, G. Madec) add pen_ddt_dds used in trdpen 21 !! - ! 2012-05 (F. Roquet) add Vallis and original JM95 equation of state 22 !! - ! 2012-05 (F. Roquet) add eos_alpha_beta 20 23 !!---------------------------------------------------------------------- 21 24 … … 28 31 !! eos_bn2 : Compute the Brunt-Vaisala frequency 29 32 !! eos_alpbet : calculates the in situ thermal/haline expansion ratio 33 !! eos_expansion_ratio : calculates the partial derivative of in situ density with respect to T and S 30 34 !! tfreez : Compute the surface freezing temperature 31 35 !! eos_init : set eos parameters (namelist) … … 51 55 END INTERFACE 52 56 53 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules 54 PUBLIC eos_init ! called by istate module 55 PUBLIC bn2 ! called by step module 56 PUBLIC eos_alpbet ! called by ldfslp module 57 PUBLIC tfreez ! called by sbcice_... modules 57 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules 58 PUBLIC eos_init ! called by istate module 59 PUBLIC bn2 ! called by step module 60 PUBLIC eos_alpbet ! called by ldfslp module 61 PUBLIC eos_alpha_beta ! used for density diagnostics in dyntra 62 PUBLIC eos_pen ! used for pe diagnostics in trdpen 63 PUBLIC tfreez ! called by sbcice_... modules 58 64 59 65 ! !!* Namelist (nameos) * 60 INTEGER , PUBLIC :: nn_eos = 0 !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.66 INTEGER , PUBLIC :: nn_eos = 0 !: = -1/0/1/2/3 type of eq. of state and associated Brunt-Vaisala frequ. 61 67 REAL(wp), PUBLIC :: rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state) 62 68 REAL(wp), PUBLIC :: rn_beta = 7.7e-4_wp !: saline expension coeff. (linear equation of state) 63 69 64 REAL(wp), PUBLIC :: ralpbet !: alpha / beta ratio 70 REAL(wp) :: vallis_ratio = 0 ! 1027/rau0 71 REAL(wp) :: vallis_diff = 0 ! (1027-rau0)/rau0 65 72 66 73 !! * Substitutions … … 82 89 !! defined through the namelist parameter nn_eos. 83 90 !! 84 !! ** Method : 3 cases: 85 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. 91 !! ** Method : 5 cases: 92 !! nn_eos = -1 : Jackett and McDougall (1995) equation of state. 93 !! Check value: rho = 1041.83267 kg/m**3 for p=3000 dbar, 94 !! t = 3 deg celcius, s=35.5 psu 95 !! nn_eos = 0 : standard NEMO equation of state, based on 96 !! Jackett and McDougall (1995) equation of state. 86 97 !! the in situ density is computed directly as a function of 87 98 !! potential temperature relative to the surface (the opa t … … 103 114 !! salinity 104 115 !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 116 !! nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 117 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 118 !! where the pressure p in decibars is approximated by the depth in meters. 105 119 !! Note that no boundary condition problem occurs in this routine 106 120 !! as pts are defined over the whole domain. … … 108 122 !! ** Action : compute prd , the in situ density (no units) 109 123 !! 110 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994111 !! ----------------------------------------------------------------------112 !! 124 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 125 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 126 !!---------------------------------------------------------------------- 113 127 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 114 128 ! ! 2 : salinity [psu] 115 129 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-] 116 ! !130 ! 117 131 INTEGER :: ji, jj, jk ! dummy loop indices 118 132 REAL(wp) :: zt , zs , zh , zsr ! local scalars … … 123 137 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 124 138 !!---------------------------------------------------------------------- 125 126 139 ! 127 140 IF( nn_timing == 1 ) CALL timing_start('eos') … … 131 144 SELECT CASE( nn_eos ) 132 145 ! 133 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 146 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==! 147 ! 134 148 !CDIR NOVERRCHK 135 149 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) … … 156 170 ! 157 171 ! add the compression terms 172 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 173 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 174 zb = zbw + ze * zs 175 ! 176 zd = -1.480266e-4_wp 177 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 178 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 179 za = ( zd*zsr + zc ) *zs + zaw 180 ! 181 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 182 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 183 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 184 zk0= ( zb1*zsr + za1 )*zs + zkw 185 ! 186 ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 187 prd(ji,jj,jk) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) & 188 & - rau0 ) * r1_rau0 * tmask(ji,jj,jk) 189 END DO 190 END DO 191 END DO 192 ! 193 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 194 ! 195 !CDIR NOVERRCHK 196 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 197 ! 198 DO jk = 1, jpkm1 199 DO jj = 1, jpj 200 DO ji = 1, jpi 201 zt = pts (ji,jj,jk,jp_tem) 202 zs = pts (ji,jj,jk,jp_sal) 203 zh = fsdept(ji,jj,jk) ! depth 204 zsr= zws (ji,jj,jk) ! square root salinity 205 ! 206 ! compute volumic mass pure water at atm pressure 207 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 208 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 209 ! seawater volumic mass atm pressure 210 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 211 & -4.0899e-3_wp ) *zt+0.824493_wp 212 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 213 zr4= 4.8314e-4_wp 214 ! 215 ! potential volumic mass (reference to the surface) 216 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 217 ! 218 ! add the compression terms 158 219 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 159 220 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp … … 187 248 END DO 188 249 ! 250 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 251 DO jk = 1, jpkm1 252 DO jj = 1, jpj 253 DO ji = 1, jpi 254 zt = pts (ji,jj,jk,jp_tem) - 10._wp 255 zs = pts (ji,jj,jk,jp_sal) - 35._wp 256 zh = fsdept(ji,jj,jk) ! depth 257 ! masked in situ density anomaly 258 prd(ji,jj,jk) = vallis_diff + vallis_ratio * ( & 259 & -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh & 260 & ) * tmask(ji,jj,jk) 261 END DO 262 END DO 263 END DO 264 ! 189 265 END SELECT 190 266 ! … … 207 283 !! namelist parameter nn_eos. 208 284 !! 209 !! ** Method : 210 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. 285 !! ** Method : 5 cases: 286 !! nn_eos = -1 : Jackett and McDougall (1995) equation of state. 287 !! nn_eos = 0 : standard NEMO equation of state, based on 288 !! Jackett and McDougall (1995) equation of state. 211 289 !! the in situ density is computed directly as a function of 212 290 !! potential temperature relative to the surface (the opa t … … 235 313 !! = rn_beta * s - rn_alpha * tn - 1. 236 314 !! rhop(t,s) = rho(t,s) 315 !! nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 316 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 317 !! rhop(t,s) = rho(t,s,0) 237 318 !! Note that no boundary condition problem occurs in this routine 238 319 !! as (tn,sn) or (ta,sa) are defined over the whole domain. … … 241 322 !! - prhop, the potential volumic mass (Kg/m3) 242 323 !! 243 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 324 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 325 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 244 326 !! Brown and Campana, Mon. Weather Rev., 1978 245 327 !!---------------------------------------------------------------------- … … 262 344 SELECT CASE ( nn_eos ) 263 345 ! 264 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!346 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==! 265 347 !CDIR NOVERRCHK 266 348 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) … … 290 372 ! 291 373 ! add the compression terms 374 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 375 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 376 zb = zbw + ze * zs 377 ! 378 zd = -1.480266e-4_wp 379 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 380 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 381 za = ( zd*zsr + zc ) *zs + zaw 382 ! 383 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 384 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 385 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 386 zk0= ( zb1*zsr + za1 )*zs + zkw 387 ! 388 ! masked in situ density anomaly 389 prd(ji,jj,jk) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) & 390 & - rau0 ) * r1_rau0 * tmask(ji,jj,jk) 391 END DO 392 END DO 393 END DO 394 ! 395 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 396 !CDIR NOVERRCHK 397 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 398 ! 399 DO jk = 1, jpkm1 400 DO jj = 1, jpj 401 DO ji = 1, jpi 402 zt = pts (ji,jj,jk,jp_tem) 403 zs = pts (ji,jj,jk,jp_sal) 404 zh = fsdept(ji,jj,jk) ! depth 405 zsr= zws (ji,jj,jk) ! square root salinity 406 ! 407 ! compute volumic mass pure water at atm pressure 408 zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt & 409 & -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 410 ! seawater volumic mass atm pressure 411 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 412 & -4.0899e-3_wp ) *zt+0.824493_wp 413 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 414 zr4= 4.8314e-4_wp 415 ! 416 ! potential volumic mass (reference to the surface) 417 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 418 ! 419 ! save potential volumic mass 420 prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk) 421 ! 422 ! add the compression terms 292 423 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 293 424 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp … … 323 454 END DO 324 455 ! 456 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 457 DO jk = 1, jpkm1 458 DO jj = 1, jpj 459 DO ji = 1, jpi 460 zt = pts (ji,jj,jk,jp_tem) - 10._wp 461 zs = pts (ji,jj,jk,jp_sal) - 35._wp 462 zh = fsdept(ji,jj,jk) ! depth 463 ! masked in situ density anomaly 464 prd(ji,jj,jk) = vallis_diff + vallis_ratio * ( & 465 & -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh & 466 & ) * tmask(ji,jj,jk) 467 ! masked in situ density anomaly 468 prhop(ji,jj,jk) = ( 1.0_wp - ( 1.67e-4_wp + 0.5e-5_wp*zt ) *zt + 0.78e-3_wp *zs ) & 469 & * 1027._wp * tmask(ji,jj,jk) 470 ! 471 END DO 472 END DO 473 END DO 474 ! 325 475 END SELECT 326 476 ! … … 342 492 !! defined through the namelist parameter nn_eos. * 2D field case 343 493 !! 344 !! ** Method : 345 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. 494 !! ** Method : 5 cases: 495 !! nn_eos = -1 : Jackett and McDougall (1995) equation of state. 496 !! nn_eos = 0 : standard NEMO equation of state, based on 497 !! Jackett and McDougall (1995) equation of state. 346 498 !! the in situ density is computed directly as a function of 347 499 !! potential temperature relative to the surface (the opa t … … 363 515 !! salinity 364 516 !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 517 !! nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 518 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 519 !! where the pressure p in decibars is approximated by the depth in meters. 365 520 !! Note that no boundary condition problem occurs in this routine 366 521 !! as pts are defined over the whole domain. … … 368 523 !! ** Action : - prd , the in situ density (no units) 369 524 !! 370 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 525 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 526 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 371 527 !!---------------------------------------------------------------------- 372 528 !! … … 391 547 SELECT CASE( nn_eos ) 392 548 ! 393 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!549 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==! 394 550 ! 395 551 !CDIR NOVERRCHK … … 403 559 DO ji = 1, fs_jpim1 ! vector opt. 404 560 zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask 405 zt = pts (ji,jj,jp_tem) ! interpolated T 406 zs = pts (ji,jj,jp_sal) ! interpolated S 561 zt = pts (ji,jj,jp_tem) ! interpolated T 562 zs = pts (ji,jj,jp_sal) ! interpolated S 563 zsr = zws (ji,jj) ! square root of interpolated S 564 zh = pdep (ji,jj) ! depth at the partial step level 565 ! 566 ! compute volumic mass pure water at atm pressure 567 zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt & 568 & -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 569 ! seawater volumic mass atm pressure 570 zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt & 571 & -4.0899e-3_wp ) *zt+0.824493_wp 572 zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 573 zr4 = 4.8314e-4_wp 574 ! 575 ! potential volumic mass (reference to the surface) 576 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 577 ! 578 ! add the compression terms 579 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 580 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 581 zb = zbw + ze * zs 582 ! 583 zd = -1.480266e-4_wp 584 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 585 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 586 za = ( zd*zsr + zc ) *zs + zaw 587 ! 588 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 589 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 590 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 591 zk0= ( zb1*zsr + za1 )*zs + zkw 592 ! 593 ! masked in situ density anomaly 594 prd(ji,jj) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) - rau0 ) / rau0 * zmask 595 END DO 596 END DO 597 ! 598 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 599 ! 600 !CDIR NOVERRCHK 601 DO jj = 1, jpjm1 602 !CDIR NOVERRCHK 603 DO ji = 1, fs_jpim1 ! vector opt. 604 zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) ) 605 END DO 606 END DO 607 DO jj = 1, jpjm1 608 DO ji = 1, fs_jpim1 ! vector opt. 609 zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask 610 zt = pts (ji,jj,jp_tem) ! interpolated T 611 zs = pts (ji,jj,jp_sal) ! interpolated S 407 612 zsr = zws (ji,jj) ! square root of interpolated S 408 613 zh = pdep (ji,jj) ! depth at the partial step level … … 455 660 END DO 456 661 ! 662 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 663 DO jj = 1, jpjm1 664 DO ji = 1, fs_jpim1 ! vector opt. 665 zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask 666 zt = pts (ji,jj,jp_tem) - 10._wp ! interpolated T 667 zs = pts (ji,jj,jp_sal) - 35._wp ! interpolated S 668 zh = pdep (ji,jj) ! depth at the partial step level 669 ! masked in situ density anomaly 670 prd(ji,jj) = vallis_diff + vallis_ratio * ( & 671 & -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh & 672 & ) * zmask 673 ! 674 END DO 675 END DO 676 ! 457 677 END SELECT 458 678 … … 473 693 !! step of the input arguments 474 694 !! 475 !! ** Method : 476 !! * nn_eos = 0 : UNESCO sea water properties 477 !! The brunt-vaisala frequency is computed using the polynomial 478 !! polynomial expression of McDougall (1987): 695 !! ** Method : 5 cases: 696 !! * nn_eos = -1 : The brunt-vaisala frequency is computed for 697 !! the Jackett and McDougall (1995) equation of state: 479 698 !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 699 !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 700 !! computed and used in zdfddm module : 701 !! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 702 !! * nn_eos = 0 : The brunt-vaisala frequency is based on the 703 !! formulation of McDougall (1987). 480 704 !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 481 705 !! computed and used in zdfddm module : … … 494 718 !! ** Action : - pn2 : the brunt-vaisala frequency 495 719 !! 496 !! References : McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987. 720 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 721 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 722 !! McDougall, J. Phys. Oceano., 1987 497 723 !!---------------------------------------------------------------------- 498 724 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] … … 501 727 !! 502 728 INTEGER :: ji, jj, jk ! dummy loop indices 503 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 729 REAL(wp) :: zgde3w, zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop ! local scalars 730 REAL(wp) :: ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM ! - - 731 REAL(wp) :: zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! - - 732 REAL(wp) :: zalpbet, zalpha, zbeta ! - - 504 733 #if defined key_zdfddm 505 734 REAL(wp) :: zds ! local scalars 506 735 #endif 736 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 507 737 !!---------------------------------------------------------------------- 508 738 509 739 ! 510 740 IF( nn_timing == 1 ) CALL timing_start('bn2') 741 ! 742 CALL wrk_alloc( jpi, jpj, jpk, zws ) 511 743 ! 512 744 ! pn2 : interior points only (2=< jk =< jpkm1 ) … … 515 747 SELECT CASE( nn_eos ) 516 748 ! 517 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 749 CASE( -1 ) !== Jackett and McDougall (1995) ==! 750 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 751 DO jk = 2, jpkm1 752 DO jj = 1, jpj 753 DO ji = 1, jpi 754 zgde3w = grav / fse3w(ji,jj,jk) 755 zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-pt 756 zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) ! salinity at w-pt 757 zsr = 0.5 * ( zws(ji,jj,jk) + zws(ji,jj,jk-1) ) ! square root of S at w-pt 758 zh = fsdepw(ji,jj,jk) ! depth in meters at w-point 759 ! 760 ! compute volumic mass pure water at atm pressure 761 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 762 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 763 ! seawater volumic mass atm pressure 764 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 765 & -4.0899e-3_wp ) *zt+0.824493_wp 766 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 767 zr4= 4.8314e-4_wp 768 ! 769 ! potential volumic mass (reference to the surface) 770 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 771 ! 772 ! add the compression terms 773 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 774 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 775 zb = zbw + ze * zs 776 ! 777 zd = -1.480266e-4_wp 778 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 779 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 780 za = ( zd*zsr + zc ) *zs + zaw 781 ! 782 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 783 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 784 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 785 zk0= ( zb1*zsr + za1 )*zs + zkw 786 ! 787 zM= ( zb*zh - za )*zh + zk0 788 zDenom= 1._wp - zh / zM 789 ! compute in-situ density 790 zrho = zrhop/zDenom 791 ! alpha 792 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 793 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 794 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 795 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 796 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 797 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 798 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 799 zdM = (zdzb*zh-zdza)*zh+zdzk0 800 zdDenom = zh * zdM / (zM*zM) 801 zalpha = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 802 ! beta 803 zdzb = ze 804 zdza = -3.0644505e-2_wp*zsr+zc 805 zdzk0 = 1.5_wp*zb1*zsr+za1 806 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 807 zdM = (zdzb*zh-zdza)*zh+zdzk0 808 zdDenom = zh * zdM / (zM*zM) 809 zbeta = ( zdrhop - zrho*zdDenom ) / zDenom / rau0 810 ! alpha/beta 811 zalpbet = zalpha / zbeta 812 ! N2 813 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 814 & * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 815 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 816 #if defined key_zdfddm 817 ! !!bug **** caution a traiter zds=dk[S]= 0 !!!! 818 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s]) 819 IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 820 rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 821 #endif 822 END DO 823 END DO 824 END DO 825 ! 826 CASE( 0 ) !== McDougall (1987) formulation ==! 518 827 DO jk = 2, jpkm1 519 828 DO jj = 1, jpj … … 524 833 zh = fsdepw(ji,jj,jk) ! depth in meters at w-point 525 834 ! 526 zal bet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt& ! ratio alpha/beta835 zalpbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & ! ratio alpha/beta 527 836 & - 0.203814e-03_wp ) * zt & 528 837 & + 0.170907e-01_wp ) * zt & … … 550 859 ! 551 860 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 552 & * ( zal bet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) &861 & * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 553 862 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 554 863 #if defined key_zdfddm … … 556 865 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s]) 557 866 IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 558 rrau(ji,jj,jk) = zal bet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds867 rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 559 868 #endif 560 869 END DO … … 564 873 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 565 874 DO jk = 2, jpkm1 566 pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 875 pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) & 876 & / fse3w(:,:,jk) * tmask(:,:,jk) 567 877 END DO 568 878 ! … … 579 889 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 580 890 IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 581 rrau(ji,jj,jk) = r alpbet* ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds891 rrau(ji,jj,jk) = rn_alpha / rn_beta * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 582 892 END DO 583 893 END DO 584 894 END DO 585 895 #endif 896 ! 897 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 898 DO jk = 2, jpkm1 899 DO jj = 1, jpj 900 DO ji = 1, jpi 901 zgde3w = grav / fse3w(ji,jj,jk) 902 zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) - 10._wp ! potential temperature at w-pt 903 zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35._wp ! salinity anomaly (s-35) at w-pt 904 zh = fsdepw(ji,jj,jk) ! depth in meters at w-point 905 ! 906 zalpha = 1027._wp * ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) / rau0 ! alpha 907 zbeta = 1027._wp * 0.78e-3_wp / rau0 ! beta 908 zalpbet = zalpha / zbeta ! ratio alpha/beta 909 ! 910 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 911 & * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 912 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 913 #if defined key_zdfddm 914 ! !!bug **** caution a traiter zds=dk[S]= 0 !!!! 915 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s]) 916 IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 917 rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 918 #endif 919 END DO 920 END DO 921 END DO 922 ! 586 923 END SELECT 587 924 … … 591 928 #endif 592 929 ! 930 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 931 ! 593 932 IF( nn_timing == 1 ) CALL timing_stop('bn2') 594 933 ! … … 603 942 !! 604 943 !! ** Method : calculates alpha / beta ratio at T-points 605 !! * nn_eos = 0 : UNESCO sea water properties 606 !! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial 607 !! polynomial expression of McDougall (1987). 944 !! * nn_eos = -1 : alpha / beta ratio is computed for 945 !! the Jackett and McDougall (1995) equation of state: 946 !! Scalar beta0 is returned = 1. 947 !! * nn_eos = 0 : alpha / beta ratio using the formulation of McDougall (1987). 608 948 !! Scalar beta0 is returned = 1. 609 949 !! * nn_eos = 1 : linear equation of state (temperature only) … … 611 951 !! Scalar beta0 is returned = 0. 612 952 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 613 !! The alpha/beta ratio is returned as ralpbet 953 !! The alpha/beta ratio is returned as palpbet 954 !! Scalar beta0 is returned = 1. 955 !! * nn_eos = 3 : Vallis equation of state (Vallis 2006) 614 956 !! Scalar beta0 is returned = 1. 615 957 !! … … 622 964 !! 623 965 INTEGER :: ji, jj, jk ! dummy loop indices 624 REAL(wp) :: zt, zs, zh ! local scalars 966 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop ! temporary scalars 967 REAL(wp) :: ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM ! - - 968 REAL(wp) :: zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 969 REAL(wp) :: zalpha, zbeta ! local scalars 970 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 625 971 !!---------------------------------------------------------------------- 626 972 ! 627 973 IF( nn_timing == 1 ) CALL timing_start('eos_alpbet') 628 974 ! 975 CALL wrk_alloc( jpi, jpj, jpk, zws ) 976 ! 629 977 SELECT CASE ( nn_eos ) 630 978 ! 631 CASE ( 0 ) ! Jackett and McDougall (1994) formulation 979 CASE ( -1 ) ! Jackett and McDougall (1995) formulation 980 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 981 DO jk = 1, jpkm1 982 DO jj = 1, jpj 983 DO ji = 1, jpi 984 zt = pts (ji,jj,jk,jp_tem) 985 zs = pts (ji,jj,jk,jp_sal) 986 zh = fsdept(ji,jj,jk) ! depth 987 zsr= zws (ji,jj,jk) ! square root salinity 988 ! 989 ! compute volumic mass pure water at atm pressure 990 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 991 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 992 ! seawater volumic mass atm pressure 993 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 994 & -4.0899e-3_wp ) *zt+0.824493_wp 995 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 996 zr4= 4.8314e-4_wp 997 ! 998 ! potential volumic mass (reference to the surface) 999 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 1000 ! 1001 ! add the compression terms 1002 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 1003 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 1004 zb = zbw + ze * zs 1005 ! 1006 zd = -1.480266e-4_wp 1007 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 1008 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 1009 za = ( zd*zsr + zc ) *zs + zaw 1010 ! 1011 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 1012 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 1013 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 1014 zk0= ( zb1*zsr + za1 )*zs + zkw 1015 ! 1016 zM= ( zb*zh - za )*zh + zk0 1017 zDenom= 1._wp - zh / zM 1018 ! compute in-situ density 1019 zrho = zrhop/zDenom 1020 ! alpha 1021 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1022 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 1023 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 1024 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 1025 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 1026 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 1027 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1028 zdM = (zdzb*zh-zdza)*zh+zdzk0 1029 zdDenom = zh * zdM / (zM*zM) 1030 zalpha = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 1031 ! beta 1032 zdzb = ze 1033 zdza = -3.0644505e-2_wp*zsr+zc 1034 zdzk0 = 1.5_wp*zb1*zsr+za1 1035 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 1036 zdM = (zdzb*zh-zdza)*zh+zdzk0 1037 zdDenom = zh * zdM / (zM*zM) 1038 zbeta = ( zdrhop - zrho*zdDenom ) / zDenom / rau0 1039 ! alpha/beta 1040 palpbet(ji,jj,jk) = zalpha / zbeta * tmask(ji,jj,jk) 1041 END DO 1042 END DO 1043 END DO 1044 beta0 = 1._wp 1045 ! 1046 CASE ( 0 ) ! McDougall (1987) formulation 632 1047 DO jk = 1, jpk 633 1048 DO jj = 1, jpj … … 660 1075 ! 661 1076 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 662 palpbet(:,:,:) = ralpbet 1077 palpbet(:,:,:) = rn_alpha / rn_beta 1078 beta0 = 1._wp 1079 ! 1080 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 1081 DO jk = 1, jpkm1 1082 DO jj = 1, jpj 1083 DO ji = 1, jpi 1084 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! potential temperature 1085 zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35) 1086 zh = fsdepw(ji,jj,jk) ! depth in meters at w-point 1087 ! 1088 zalpha = ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) ! alpha/vallis_ratio 1089 zbeta = 0.78e-3_wp ! beta/vallis_ratio 1090 ! 1091 palpbet(ji,jj,jk) = zalpha / zbeta * tmask(ji,jj,jk) 1092 END DO 1093 END DO 1094 END DO 663 1095 beta0 = 1._wp 664 1096 ! … … 670 1102 END SELECT 671 1103 ! 1104 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 1105 ! 672 1106 IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet') 673 1107 ! … … 675 1109 676 1110 1111 SUBROUTINE eos_alpha_beta( pts, palpha, pbeta ) 1112 !!---------------------------------------------------------------------- 1113 !! *** ROUTINE eos_alpha_beta *** 1114 !! 1115 !! ** Purpose : Calculates the in situ thermal/haline expansion terms at T-points 1116 !! 1117 !! ** Method : calculates alpha and beta at T-points 1118 !! * nn_eos = -1 : Jackett and McDougall (1995) equation of state 1119 !! * nn_eos = 0 : modified Jackett and McDougall (1995) equation of state 1120 !! * nn_eos = 1 : linear equation of state (temperature only) 1121 !! palpha = rn_alpha 1122 !! pbeta = 0 1123 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 1124 !! palpha = rn_alpha 1125 !! pbeta = rn_beta 1126 !! * nn_eos = 3 : Vallis equation of state (Vallis 2006) 1127 !! alpha and beta ratios are returned as 3-D arrays. 1128 !! 1129 !! ** Action : - palpha : thermal expansion ratio at T-points 1130 !! : - pbeta : haline expansion ratio at T-points 1131 !!---------------------------------------------------------------------- 1132 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 1133 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpha ! alpha 1134 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pbeta ! beta 1135 ! 1136 INTEGER :: ji, jj, jk ! dummy loop indices 1137 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop ! temporary scalars 1138 REAL(wp) :: ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM ! - - 1139 REAL(wp) :: zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 1140 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 1141 !!---------------------------------------------------------------------- 1142 ! 1143 IF ( nn_timing == 1 ) CALL timing_start('eos_alpha_beta') 1144 ! 1145 CALL wrk_alloc( jpi, jpj, jpk, zws ) 1146 ! 1147 SELECT CASE ( nn_eos ) 1148 ! 1149 CASE ( -1 ) ! Jackett and McDougall (1995) formulation 1150 ! 1151 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 1152 DO jk = 1, jpkm1 1153 DO jj = 1, jpj 1154 DO ji = 1, jpi 1155 zt = pts (ji,jj,jk,jp_tem) 1156 zs = pts (ji,jj,jk,jp_sal) 1157 zh = fsdept(ji,jj,jk) ! depth 1158 zsr= zws (ji,jj,jk) ! square root salinity 1159 ! 1160 ! compute volumic mass pure water at atm pressure 1161 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 1162 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 1163 ! seawater volumic mass atm pressure 1164 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 1165 & -4.0899e-3_wp ) *zt+0.824493_wp 1166 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 1167 zr4= 4.8314e-4_wp 1168 ! 1169 ! potential volumic mass (reference to the surface) 1170 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 1171 ! 1172 ! add the compression terms 1173 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 1174 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 1175 zb = zbw + ze * zs 1176 ! 1177 zd = -1.480266e-4_wp 1178 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 1179 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 1180 za = ( zd*zsr + zc ) *zs + zaw 1181 ! 1182 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 1183 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 1184 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 1185 zk0= ( zb1*zsr + za1 )*zs + zkw 1186 ! 1187 zM= ( zb*zh - za )*zh + zk0 1188 zDenom= 1._wp - zh / zM 1189 ! compute in-situ density 1190 zrho = zrhop/zDenom 1191 ! alpha 1192 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1193 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 1194 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 1195 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 1196 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 1197 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 1198 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1199 zdM = (zdzb*zh-zdza)*zh+zdzk0 1200 zdDenom = zh * zdM / (zM*zM) 1201 palpha(ji,jj,jk) = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 1202 ! beta 1203 zdzb = ze 1204 zdza = -3.0644505e-2_wp*zsr+zc 1205 zdzk0 = 1.5_wp*zb1*zsr+za1 1206 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 1207 zdM = (zdzb*zh-zdza)*zh+zdzk0 1208 zdDenom = zh * zdM / (zM*zM) 1209 pbeta(ji,jj,jk) = ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 1210 END DO 1211 END DO 1212 END DO 1213 ! 1214 CASE ( 0 ) ! modified Jackett and McDougall (1995) formulation 1215 ! 1216 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 1217 DO jk = 1, jpkm1 1218 DO jj = 1, jpj 1219 DO ji = 1, jpi 1220 zt = pts (ji,jj,jk,jp_tem) 1221 zs = pts (ji,jj,jk,jp_sal) 1222 zh = fsdept(ji,jj,jk) ! depth 1223 zsr= zws (ji,jj,jk) ! square root salinity 1224 ! 1225 ! compute volumic mass pure water at atm pressure 1226 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 1227 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 1228 ! seawater volumic mass atm pressure 1229 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 1230 & -4.0899e-3_wp ) *zt+0.824493_wp 1231 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 1232 zr4= 4.8314e-4_wp 1233 ! 1234 ! potential volumic mass (reference to the surface) 1235 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 1236 ! 1237 ! add the compression terms 1238 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 1239 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 1240 zb = zbw + ze * zs 1241 ! 1242 zd = -2.042967e-2_wp 1243 zc = (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 1244 zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 1245 za = ( zd*zsr + zc ) *zs + zaw 1246 ! 1247 zb1= (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp 1248 za1= ( ( 2.326469e-3_wp*zt+1.553190_wp) *zt-65.00517_wp ) *zt+1044.077_wp 1249 zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 1250 zk0= ( zb1*zsr + za1 )*zs + zkw 1251 ! 1252 zM= ( zb*zh - za )*zh + zk0 1253 zDenom= 1._wp - zh / zM 1254 ! compute in-situ density 1255 zrho = zrhop/zDenom 1256 ! alpha 1257 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1258 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 1259 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 1260 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 1261 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 1262 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 1263 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1264 zdM = (zdzb*zh-zdza)*zh+zdzk0 1265 zdDenom = zh * zdM / (zM*zM) 1266 palpha(ji,jj,jk) = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 1267 ! beta 1268 zdzb = ze 1269 zdza = -3.0644505e-2_wp*zsr+zc 1270 zdzk0 = 1.5_wp*zb1*zsr+za1 1271 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 1272 zdM = (zdzb*zh-zdza)*zh+zdzk0 1273 zdDenom = zh * zdM / (zM*zM) 1274 pbeta(ji,jj,jk) = ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 1275 END DO 1276 END DO 1277 END DO 1278 ! 1279 CASE ( 1 ) !== Linear formulation = F( temperature ) ==! 1280 palpha(:,:,:) = rn_alpha 1281 pbeta (:,:,:) = 0._wp 1282 ! 1283 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 1284 palpha(:,:,:) = rn_alpha 1285 pbeta (:,:,:) = rn_beta 1286 ! 1287 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 1288 DO jk = 1, jpkm1 1289 DO jj = 1, jpj 1290 DO ji = 1, jpi 1291 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! potential temperature (t-10) 1292 zh = fsdept(ji,jj,jk) ! depth in meters at w-point 1293 palpha(ji,jj,jk) = vallis_ratio * & 1294 & ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) * tmask(ji,jj,jk) ! alpha 1295 ! 1296 END DO 1297 END DO 1298 END DO 1299 pbeta (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:) ! beta 1300 ! 1301 CASE DEFAULT 1302 IF(lwp) WRITE(numout,cform_err) 1303 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 1304 nstop = nstop + 1 1305 ! 1306 END SELECT 1307 ! 1308 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 1309 ! 1310 IF( nn_timing == 1 ) CALL timing_stop('eos_alpha_beta') 1311 ! 1312 END SUBROUTINE eos_alpha_beta 1313 1314 1315 1316 SUBROUTINE eos_pen( pts, palphaPE, pbetaPE, pPEanom ) 1317 !!---------------------------------------------------------------------- 1318 !! *** ROUTINE eos_pen *** 1319 !! 1320 !! ** Purpose : Calculates alpha_PE, beta_PE and PE at T-points 1321 !! 1322 !! ** Method : PE is defined analytically as the vertical 1323 !! primitive of EOS times -g integrated between 0 and z>0. 1324 !! PEanom is the PE anomaly: PEanom = ( PE - rau0 gz ) / rau0 gz 1325 !! = -1/z * /int_z^0 rd , where rd density anomaly 1326 !! dalphaPE and dbetaPE are partial derivatives of PE anomaly with respect to T and S: 1327 !! alphaPE = - 1/(rau0 gz) * dPE/dT = - dPEanom/dT 1328 !! betaPE = 1/(rau0 gz) * dPE/dS = - dPEanom/dS 1329 !! 1330 !! * nn_eos = -1 : Jackett and McDougall (1995) equation of state. 1331 !! * nn_eos = 0 : modified Jackett and McDougall (1995) equation of state. 1332 !! * nn_eos = 1 : linear equation of state (temperature only) 1333 !! palpha = rau0*g*rn_alpha*z 1334 !! pbeta = 0 1335 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 1336 !! palpha = rau0*g*rn_alpha*z 1337 !! pbeta = -rau0*g*rn_beta*z 1338 !! * nn_eos = 3 : Vallis equation of state (Vallis 2006) 1339 !! 1340 !! ** Action : - pPEanom : PE anomaly given at T-points 1341 !! : - palphaPE : given at T-points 1342 !! : - pbetaPE : given at T-points 1343 !!---------------------------------------------------------------------- 1344 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 1345 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palphaPE ! partial derivative of PE anom 1346 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pbetaPE ! with respect to T and S, resp. 1347 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pPEanom ! potential energy anomaly 1348 ! 1349 INTEGER :: ji, jj, jk ! dummy loop indices 1350 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop ! temporary scalars 1351 REAL(wp) :: ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM ! - - 1352 REAL(wp) :: zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 1353 REAL(wp) :: zap1, zdelta, zsdelta, zAp, zBp, zlogBp, zCp, zatanCp ! 1354 REAL(wp) :: zddelta, zdAp, zdBp, zdlogBp, zdCp, zP, zdP, zEp ! 1355 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 1356 !!---------------------------------------------------------------------- 1357 ! 1358 IF ( nn_timing == 1 ) CALL timing_start('eos_pen') 1359 ! 1360 CALL wrk_alloc( jpi, jpj, jpk, zws ) 1361 ! 1362 SELECT CASE ( nn_eos ) 1363 ! 1364 CASE ( -1 ) ! Jackett and McDougall (1995) formulation 1365 ! 1366 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 1367 DO jk = 1, jpkm1 1368 DO jj = 1, jpj 1369 DO ji = 1, jpi 1370 zt = pts (ji,jj,jk,jp_tem) 1371 zs = pts (ji,jj,jk,jp_sal) 1372 zh = fsdept(ji,jj,jk) ! depth 1373 zsr= zws (ji,jj,jk) ! square root salinity 1374 ! 1375 ! compute volumic mass pure water at atm pressure 1376 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 1377 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 1378 ! seawater volumic mass atm pressure 1379 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 1380 & -4.0899e-3_wp ) *zt+0.824493_wp 1381 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 1382 zr4= 4.8314e-4_wp 1383 ! 1384 ! potential volumic mass (reference to the surface) 1385 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 1386 ! 1387 ! add the compression terms 1388 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 1389 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 1390 zb = zbw + ze * zs 1391 ! 1392 zd = -1.480266e-4_wp 1393 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 1394 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 1395 za = ( zd*zsr + zc ) *zs + zaw 1396 ! 1397 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 1398 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 1399 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 1400 zk0= ( zb1*zsr + za1 )*zs + zkw 1401 ! 1402 zM= ( zb*zh - za )*zh + zk0 1403 zDenom= 1._wp - zh / zM 1404 ! compute in-situ density 1405 zrho = zrhop/zDenom 1406 ! 1407 zap1 = 1._wp + za 1408 zdelta = 4._wp*zb*zk0 - zap1*zap1 1409 zsdelta = sqrt( abs( zdelta ) ) 1410 zAp = zap1 / zb / zsdelta 1411 zBp = ( zM - zh ) / zk0 1412 zlogBp = log(abs(zBp)) 1413 zCp = zh * zsdelta / (2._wp*zk0-zh*zap1) 1414 zatanCp = atan(zCp) 1415 zP = zh + zAp * zatanCp + 0.5_wp*zlogBp/zb 1416 ! compute potential energy 1417 pPEanom(ji,jj,jk) = ( ( zrhop * zP ) / rau0 / zh - 1._wp ) * tmask(ji,jj,jk) !!wrong 1418 ! 1419 ! dPEdt 1420 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1421 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 1422 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 1423 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 1424 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 1425 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 1426 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1427 zdM = (zdzb*zh-zdza)*zh+zdzk0 1428 zdDenom = zh * zdM / (zM*zM) 1429 ! 1430 zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 1431 zdAp = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 1432 zdlogBp = zdM/(zM-zh) - zdzk0/zk0 1433 zdCp = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 1434 zdP = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 1435 ! 1436 palphaPE(ji,jj,jk) = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 1437 ! 1438 ! dPEds 1439 zdzb = ze 1440 zdza = -3.0644505e-2_wp*zsr+zc 1441 zdzk0 = 1.5_wp*zb1*zsr+za1 1442 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 1443 zdM = (zdzb*zh-zdza)*zh+zdzk0 1444 zdDenom = zh * zdM / (zM*zM) 1445 ! 1446 zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 1447 zdAp = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 1448 zdlogBp = zdM/(zM-zh) - zdzk0/zk0 1449 zdCp = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 1450 zdP = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 1451 ! 1452 pbetaPE(ji,jj,jk) = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 1453 ! 1454 END DO 1455 END DO 1456 END DO 1457 ! 1458 CASE ( 0 ) ! modified Jackett and McDougall (1995) formulation 1459 ! 1460 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 1461 DO jk = 1, jpkm1 1462 DO jj = 1, jpj 1463 DO ji = 1, jpi 1464 zt = pts (ji,jj,jk,jp_tem) 1465 zs = pts (ji,jj,jk,jp_sal) 1466 zh = fsdept(ji,jj,jk) ! depth 1467 zsr= zws (ji,jj,jk) ! square root salinity 1468 ! 1469 ! compute volumic mass pure water at atm pressure 1470 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 1471 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 1472 ! seawater volumic mass atm pressure 1473 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 1474 & -4.0899e-3_wp ) *zt+0.824493_wp 1475 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 1476 zr4= 4.8314e-4_wp 1477 ! 1478 ! potential volumic mass (reference to the surface) 1479 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 1480 ! 1481 ! add the compression terms 1482 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 1483 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 1484 zb = zbw + ze * zs 1485 ! 1486 zd = -2.042967e-2_wp 1487 zc = (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 1488 zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 1489 za = ( zd*zsr + zc ) *zs + zaw 1490 ! 1491 zb1= (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp 1492 za1= ( ( 2.326469e-3_wp*zt+1.553190_wp) *zt-65.00517_wp ) *zt+1044.077_wp 1493 zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 1494 zk0= ( zb1*zsr + za1 )*zs + zkw 1495 ! 1496 zM= ( zb*zh - za )*zh + zk0 1497 zDenom= 1._wp - zh / zM 1498 ! compute in-situ density 1499 zrho = zrhop/zDenom 1500 ! 1501 zap1 = 1._wp + za 1502 zdelta = 4._wp*zb*zk0 - zap1*zap1 1503 zsdelta = sqrt( abs( zdelta ) ) 1504 zAp = zap1 / zb / zsdelta 1505 zBp = ( zM - zh ) / zk0 1506 zlogBp = log(abs(zBp)) 1507 zCp = zh * zsdelta / (2._wp*zk0-zh*zap1) 1508 zatanCp = atan(zCp) 1509 zP = zh + zAp * zatanCp + 0.5_wp*zlogBp/zb 1510 ! compute potential energy 1511 pPEanom(ji,jj,jk) = - grav * zrhop * zP * tmask(ji,jj,jk) 1512 ! 1513 ! dPEdt 1514 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1515 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 1516 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 1517 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 1518 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 1519 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 1520 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1521 zdM = (zdzb*zh-zdza)*zh+zdzk0 1522 zdDenom = zh * zdM / (zM*zM) 1523 ! 1524 zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 1525 zdAp = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 1526 zdlogBp = zdM/(zM-zh) - zdzk0/zk0 1527 zdCp = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 1528 zdP = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 1529 ! 1530 palphaPE(ji,jj,jk) = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 1531 ! 1532 ! dPEds 1533 zdzb = ze 1534 zdza = -3.0644505e-2_wp*zsr+zc 1535 zdzk0 = 1.5_wp*zb1*zsr+za1 1536 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 1537 zdM = (zdzb*zh-zdza)*zh+zdzk0 1538 zdDenom = zh * zdM / (zM*zM) 1539 ! 1540 zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 1541 zdAp = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 1542 zdlogBp = zdM/(zM-zh) - zdzk0/zk0 1543 zdCp = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 1544 zdP = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 1545 ! 1546 pbetaPE(ji,jj,jk) = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 1547 ! 1548 END DO 1549 END DO 1550 END DO 1551 ! 1552 CASE ( 1 ) !== Linear formulation = F( temperature ) ==! 1553 palphaPE(:,:,:) = rn_alpha * tmask(:,:,:) 1554 pbetaPE (:,:,:) = 0._wp 1555 DO jk = 1, jpkm1 1556 pPEanom(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 1557 END DO 1558 ! 1559 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 1560 palphaPE(:,:,:) = rn_alpha * tmask(:,:,:) 1561 pbetaPE (:,:,:) = rn_beta * tmask(:,:,:) 1562 DO jk = 1, jpkm1 1563 pPEanom(:,:,jk) = ( rn_beta * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 1564 END DO 1565 ! 1566 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 1567 DO jk = 1, jpkm1 1568 DO jj = 1, jpj 1569 DO ji = 1, jpi 1570 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! potential temperature (t-10) 1571 zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35) 1572 zh = fsdept(ji,jj,jk) ! depth in meters at w-point 1573 ! 1574 palphaPE(ji,jj,jk) = vallis_ratio * & 1575 & ( 1.67e-4_wp * ( 1._wp + 0.55e-4_wp*zh ) + 1.e-5_wp*zt ) * tmask(ji,jj,jk) ! alphaPE 1576 pPEanom(ji,jj,jk) = vallis_diff + vallis_ratio * & 1577 & ( - ( 1.67e-4_wp*(1._wp+0.55e-4_wp*zh) + 0.5e-5_wp*zt )*zt + 0.78e-3_wp*zs + 2.195e-6_wp*zh ) & 1578 & * tmask(ji,jj,jk) 1579 ! 1580 END DO 1581 END DO 1582 END DO 1583 pbetaPE (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:) ! betaPE=beta 1584 ! 1585 CASE DEFAULT 1586 IF(lwp) WRITE(numout,cform_err) 1587 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 1588 nstop = nstop + 1 1589 ! 1590 END SELECT 1591 ! 1592 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 1593 ! 1594 IF( nn_timing == 1 ) CALL timing_stop('eos_pen') 1595 ! 1596 END SUBROUTINE eos_pen 1597 1598 1599 677 1600 FUNCTION tfreez( psal ) RESULT( ptf ) 678 1601 !!---------------------------------------------------------------------- 679 !! *** ROUTINE eos_init***1602 !! *** ROUTINE tfreez *** 680 1603 !! 681 1604 !! ** Purpose : Compute the sea surface freezing temperature [Celcius] … … 724 1647 SELECT CASE( nn_eos ) ! check option 725 1648 ! 726 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!1649 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==! 727 1650 IF(lwp) WRITE(numout,*) 728 IF(lwp) WRITE(numout,*) ' use of Jackett & McDougall (1994) equation of state and' 729 IF(lwp) WRITE(numout,*) ' McDougall (1987) Brunt-Vaisala frequency' 1651 IF(lwp) WRITE(numout,*) ' use of Jackett & McDougall (1995) equation of state' 1652 ! 1653 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 1654 IF(lwp) WRITE(numout,*) 1655 IF(lwp) WRITE(numout,*) ' use of modified Jackett & McDougall (1995) equation of state' 1656 IF(lwp) WRITE(numout,*) ' and McDougall (1987) Brunt-Vaisala frequency' 730 1657 ! 731 1658 CASE( 1 ) !== Linear formulation = F( temperature ) ==! … … 736 1663 ! 737 1664 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 738 ralpbet = rn_alpha / rn_beta739 1665 IF(lwp) WRITE(numout,*) 740 1666 IF(lwp) WRITE(numout,*) ' use of linear eos rho(T,S) = rau0 * ( rn_beta * S - rn_alpha * T )' 1667 ! 1668 CASE( 3 ) !== Vallis (2006) formulation ==! 1669 IF(lwp) WRITE(numout,*) 1670 IF(lwp) WRITE(numout,*) ' use of Vallis (2006) equation of state' 1671 vallis_ratio = 1027._wp / rau0 1672 vallis_diff = (1027._wp-rau0) / rau0 741 1673 ! 742 1674 CASE DEFAULT !== ERROR in nn_eos ==! -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90
r3718 r3876 4 4 !! Ocean tracers: horizontal & vertical advective trend 5 5 !!====================================================================== 6 !! History : 8.2 ! 2001-08 (G. Madec, E. Durand)trahad+trazad=traadv7 !! 8 !! 9.0! 2004-08 (C. Talandier) New trends organization6 !! History : OPA ! 2001-08 (G. Madec, E. Durand) v8.2 trahad+trazad=traadv 7 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 8 !! - ! 2004-08 (C. Talandier) New trends organization 9 9 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 10 10 !! 2.0 ! 2006-04 (R. Benshila, G. Madec) Step reorganization … … 21 21 USE dom_oce ! ocean space and time domain 22 22 USE eosbn2 ! equation of state 23 USE trd mod_oce ! tracers trends24 USE trdtra ! tr acers trends23 USE trd_oce ! trends: ocean variables 24 USE trdtra ! trends manager: tracers 25 25 USE closea ! closed sea 26 26 USE sbcrnf ! river runoffs … … 37 37 PRIVATE 38 38 39 PUBLIC tra_adv_cen2 40 PUBLIC ups_orca_set 39 PUBLIC tra_adv_cen2 ! routine called by step.F90 40 PUBLIC ups_orca_set ! routine used by traadv_cen2_jki.F90 41 41 42 42 LOGICAL :: l_trd ! flag to compute trends … … 55 55 56 56 SUBROUTINE tra_adv_cen2( kt, kit000, cdtype, pun, pvn, pwn, & 57 & ptb, ptn, pta, kjpt )57 & ptb, ptn, pta, kjpt ) 58 58 !!---------------------------------------------------------------------- 59 59 !! *** ROUTINE tra_adv_cen2 *** … … 85 85 !! * Add this trend now to the general trend of tracer (ta,sa): 86 86 !! pta = pta + ztra 87 !! * trend diagnostic ( 'key_trdtra' defined): the trend is87 !! * trend diagnostic (l_trdtra=T or l_trctra=T): the trend is 88 88 !! saved for diagnostics. The trends saved is expressed as 89 !! Uh.gradh(T), i.e. 90 !! save trend = ztra + ptn divn 89 !! Uh.gradh(T), i.e. save trend = ztra + ptn divn 91 90 !! 92 91 !! Part II : vertical advection … … 104 103 !! Add this trend now to the general trend of tracer (ta,sa): 105 104 !! pta = pta + ztra 106 !! Trend diagnostic ( 'key_trdtra' defined): the trend is105 !! Trend diagnostic (l_trdtra=T or l_trctra=T): the trend is 107 106 !! saved for diagnostics. The trends saved is expressed as : 108 107 !! save trend = w.gradz(T) = ztra - ptn divn. … … 144 143 IF(lwp) WRITE(numout,*) 145 144 ! 146 IF 145 IF( .NOT. ALLOCATED( upsmsk ) ) THEN 147 146 ALLOCATE( upsmsk(jpi,jpj), STAT=ierr ) 148 147 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array') … … 260 259 END DO 261 260 262 ! ! trend diagnostics (contribution of upstream fluxes)261 ! ! trend diagnostics 263 262 IF( l_trd ) THEN 264 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, zwx, pun, ptn(:,:,:,jn) )265 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, zwy, pvn, ptn(:,:,:,jn) )266 CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zwz, pwn, ptn(:,:,:,jn) )263 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 264 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 265 CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 267 266 END IF 268 267 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 272 271 ENDIF 273 272 ! 274 END DO273 END DO 275 274 276 275 ! --------------------------- required in restart file to ensure restartability) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
r3718 r3876 16 16 !!---------------------------------------------------------------------- 17 17 USE oce ! ocean dynamics and active tracers 18 USE trc_oce ! share passive tracers/Ocean variables 18 19 USE dom_oce ! ocean space and time domain 19 USE trdmod_oce ! tracers trends 20 USE trdtra ! tracers trends 21 USE in_out_manager ! I/O manager 20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! tracers trends manager 22 22 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 23 USE trabbl ! tracers: bottom boundary layer 24 USE lib_mpp ! distribued memory computing 25 USE lbclnk ! ocean lateral boundary condition (or mpp link) 23 USE sbcrnf ! river runoffs 26 24 USE diaptr ! poleward transport diagnostics 27 USE trc_oce ! share passive tracers/Ocean variables25 ! 28 26 USE wrk_nemo ! Memory Allocation 29 27 USE timing ! Timing 30 28 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 31 USE eosbn2 ! equation of state 32 USE sbcrnf ! river runoffs 29 USE in_out_manager ! I/O manager 30 USE lib_mpp ! distribued memory computing 31 USE lbclnk ! ocean lateral boundary condition (or mpp link) 32 !!gm USE trabbl ! tracers: bottom boundary layer 33 !!gm USE eosbn2 ! equation of state 33 34 34 35 IMPLICIT NONE … … 191 192 zalpha = 0.5 - z0u 192 193 zu = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 193 zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * (zu * zslpx(ji+1,jj,jk))194 zzwy = ptb(ji ,jj,jk,jn) + xind(ji,jj,jk) * (zu * zslpx(ji ,jj,jk))194 zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 195 zzwy = ptb(ji ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji ,jj,jk) 195 196 zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 196 197 ! … … 198 199 zalpha = 0.5 - z0v 199 200 zv = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 200 zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * (zv * zslpy(ji,jj+1,jk))201 zzwy = ptb(ji,jj ,jk,jn) + xind(ji,jj,jk) * (zv * zslpy(ji,jj ,jk))201 zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 202 zzwy = ptb(ji,jj ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj ,jk) 202 203 zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 203 204 END DO … … 222 223 ! ! trend diagnostics (contribution of upstream fluxes) 223 224 IF( l_trd ) THEN 224 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, zwx, pun, ptb(:,:,:,jn) )225 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, zwy, pvn, ptb(:,:,:,jn) )225 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 226 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 226 227 END IF 227 228 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 273 274 zalpha = 0.5 + z0w 274 275 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr 275 zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * (zw * zslpx(ji,jj,jk+1))276 zzwy = ptb(ji,jj,jk ,jn) + xind(ji,jj,jk) * (zw * zslpx(ji,jj,jk ))276 zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 277 zzwy = ptb(ji,jj,jk ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk ) 277 278 zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 278 279 END DO … … 293 294 END DO 294 295 ! ! Save the vertical advective trends for diagnostic 295 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zwx, pwn, ptb(:,:,:,jn) )296 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 296 297 ! 297 298 ENDDO -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
r3625 r3876 13 13 !!---------------------------------------------------------------------- 14 14 USE oce ! ocean dynamics and active tracers 15 USE trc_oce ! share passive tracers/Ocean variables 15 16 USE dom_oce ! ocean space and time domain 16 USE trd mod_oce ! tracers trends17 USE trdtra ! tr acers trends17 USE trd_oce ! trends: ocean variables 18 USE trdtra ! trends manager: tracers 18 19 USE in_out_manager ! I/O manager 19 20 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 20 USE trabbl ! tracers: bottom boundary layer 21 !!gm USE trabbl ! tracers: bottom boundary layer 22 USE diaptr ! poleward transport diagnostics 23 ! 21 24 USE lib_mpp ! distribued memory computing 22 25 USE lbclnk ! ocean lateral boundary condition (or mpp link) 23 USE diaptr ! poleward transport diagnostics24 USE trc_oce ! share passive tracers/Ocean variables25 26 USE wrk_nemo ! Memory Allocation 26 27 USE timing ! Timing … … 201 202 ! ! trend diagnostics (contribution of upstream fluxes) 202 203 IF( l_trd ) THEN 203 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, zwx, pun, ptb(:,:,:,jn) )204 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, zwy, pvn, ptb(:,:,:,jn) )204 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 205 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 205 206 END IF 206 207 … … 284 285 END DO 285 286 ! ! trend diagnostics (contribution of upstream fluxes) 286 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zwx, pwn, ptb(:,:,:,jn) )287 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 287 288 ! 288 289 END DO -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
r3625 r3876 17 17 USE oce ! ocean dynamics and active tracers 18 18 USE dom_oce ! ocean space and time domain 19 USE trdmod_oce ! ocean space and time domain 20 USE trdtra ! ocean tracers trends 21 USE trabbl ! advective term in the BBL 19 USE trc_oce ! share passive tracers/Ocean variables 20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! trends manager: tracers 22 !!gm USE trabbl ! advective term in the BBL 23 USE dynspg_oce ! surface pressure gradient variables 24 USE diaptr ! poleward transport diagnostics 25 ! 22 26 USE lib_mpp ! distribued memory computing 23 27 USE lbclnk ! ocean lateral boundary condition (or mpp link) 24 USE dynspg_oce ! surface pressure gradient variables25 28 USE in_out_manager ! I/O manager 26 USE diaptr ! poleward transport diagnostics27 USE trc_oce ! share passive tracers/Ocean variables28 29 USE wrk_nemo ! Memory Allocation 29 30 USE timing ! Timing … … 233 234 END DO 234 235 ! ! trend diagnostics (contribution of upstream fluxes) 235 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, zwx, pun, ptn(:,:,:,jn) )236 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 236 237 ! 237 238 END DO … … 359 360 END DO 360 361 ! ! trend diagnostics (contribution of upstream fluxes) 361 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, zwy, pvn, ptn(:,:,:,jn) )362 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 362 363 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 363 364 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN … … 422 423 END DO 423 424 ! ! Save the vertical advective trends for diagnostic 424 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zwz, pwn, ptn(:,:,:,jn) )425 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 425 426 ! 426 427 END DO -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r3625 r3876 22 22 USE oce ! ocean dynamics and active tracers 23 23 USE dom_oce ! ocean space and time domain 24 USE trdmod_oce ! tracers trends 24 USE trc_oce ! share passive tracers/Ocean variables 25 USE trd_oce ! trends: ocean variables 25 26 USE trdtra ! tracers trends 26 USE in_out_manager ! I/O manager27 27 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 28 USE diaptr ! poleward transport diagnostics 29 ! 28 30 USE lib_mpp ! MPP library 29 31 USE lbclnk ! ocean lateral boundary condition (or mpp link) 30 USE diaptr ! poleward transport diagnostics 31 USE trc_oce ! share passive tracers/Ocean variables 32 USE in_out_manager ! I/O manager 32 33 USE wrk_nemo ! Memory Allocation 33 34 USE timing ! Timing … … 228 229 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 229 230 230 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, ztrdx, pun, ptn(:,:,:,jn) )231 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )232 CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )231 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 232 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 233 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 233 234 END IF 234 235 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r3787 r3876 14 14 USE oce ! ocean dynamics and active tracers 15 15 USE dom_oce ! ocean space and time domain 16 USE trdmod_oce ! ocean space and time domain 17 USE trdtra 18 USE lib_mpp 16 USE trc_oce ! share passive tracers/Ocean variables 17 USE trd_oce ! trends: ocean variables 18 USE trdtra ! trends manager: tracers 19 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 20 USE diaptr ! poleward transport diagnostics 21 ! 22 USE lib_mpp ! I/O library 19 23 USE lbclnk ! ocean lateral boundary condition (or mpp link) 20 24 USE in_out_manager ! I/O manager 21 USE diaptr ! poleward transport diagnostics22 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient23 USE trc_oce ! share passive tracers/Ocean variables24 25 USE wrk_nemo ! Memory Allocation 25 26 USE timing ! Timing … … 51 52 !! and add it to the general trend of passive tracer equations. 52 53 !! 53 !! ** Method : The upstream biased 3rd order scheme (UBS) is based on an54 !! ** Method : The upstream biased scheme (UBS) is based on a 3rd order 54 55 !! upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 55 56 !! It is only used in the horizontal direction. … … 182 183 ! ! trend diagnostics (contribution of upstream fluxes) 183 184 IF( l_trd ) THEN 184 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, zwx, pun, ptn(:,:,:,jn) )185 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, zwy, pvn, ptn(:,:,:,jn) )185 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 186 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 186 187 END IF 187 188 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 265 266 END DO 266 267 END DO 267 CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zltv )268 CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv ) 268 269 ENDIF 269 270 ! -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
r3625 r3876 18 18 USE dom_oce ! domain: ocean 19 19 USE phycst ! physical constants 20 USE trd mod_oce ! trends: ocean variables21 USE trdtra ! trends : activetracers20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! trends manager: tracers 22 22 USE in_out_manager ! I/O manager 23 23 USE prtctl ! Print control … … 99 99 IF( l_trdtra ) THEN ! Save the geothermal heat flux trend for diagnostics 100 100 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 101 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_bbc, ztrdt )101 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt ) 102 102 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 103 103 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
r3764 r3876 28 28 USE phycst ! physical constant 29 29 USE eosbn2 ! equation of state 30 USE trd mod_oce! trends: ocean variables31 USE trdtra ! trends : active tracers32 USE iom ! IOM server 30 USE trd_oce ! trends: ocean variables 31 USE trdtra ! trends manager: tracers 32 USE iom ! IOM server 33 33 USE in_out_manager ! I/O manager 34 34 USE lbclnk ! ocean lateral boundary conditions … … 148 148 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 149 149 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 150 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_bbl, ztrdt )151 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_bbl, ztrds )152 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 150 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt ) 151 CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds ) 152 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 153 153 ENDIF 154 154 ! … … 180 180 !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 181 181 !!---------------------------------------------------------------------- 182 !183 182 INTEGER , INTENT(in ) :: kjpt ! number of tracers 184 183 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields … … 399 398 - 0.121555e-07 ) * pfh 400 399 !!---------------------------------------------------------------------- 401 402 400 ! 403 401 IF( nn_timing == 1 ) CALL timing_start( 'bbl') … … 405 403 CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 406 404 ! 407 408 405 IF( kt == kit000 ) THEN 409 406 IF(lwp) WRITE(numout,*) … … 411 408 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 412 409 ENDIF 413 414 410 ! !* bottom temperature, salinity, velocity and depth 415 411 #if defined key_vectopt_loop -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
r3294 r3876 27 27 USE oce ! ocean: variables 28 28 USE dom_oce ! ocean: domain variables 29 USE trd mod_oce ! ocean: trendvariables30 USE trdtra ! active tracers: trends29 USE trd_oce ! trends: ocean variables 30 USE trdtra ! trends manager: tracers 31 31 USE zdf_oce ! ocean: vertical physics 32 32 USE phycst ! physical constants … … 47 47 PUBLIC dtacof_zoom ! routine called by in both tradmp.F90 and trcdmp.F90 48 48 49 ! !!* Namelist namtra_dmp : T & S newtonian damping *49 ! !!* Namelist namtra_dmp : T & S newtonian damping * 50 50 LOGICAL, PUBLIC :: ln_tradmp = .TRUE. !: internal damping flag 51 51 INTEGER :: nn_hdmp = -1 ! = 0/-1/'latitude' for damping over T and S … … 111 111 ! 112 112 CALL wrk_alloc( jpi, jpj, jpk, jpts, zts_dta ) 113 ! 113 114 ! !== input T-S data at kt ==! 114 115 CALL dta_tsd( kt, zts_dta ) ! read and interpolates T-S data at kt … … 171 172 ! 172 173 IF( l_trdtra ) THEN ! trend diagnostic 173 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_dmp, ttrdmp )174 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_dmp, strdmp )174 CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ttrdmp ) 175 CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, strdmp ) 175 176 ENDIF 176 177 ! ! Control print -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
r3294 r3876 23 23 USE traldf_iso_grif ! lateral mixing (tra_ldf_iso_grif routine) 24 24 USE traldf_lap ! lateral mixing (tra_ldf_lap routine) 25 USE trd mod_oce ! ocean space and time domain26 USE trdtra ! ocean active tracers trends25 USE trd_oce ! trends: ocean variables 26 USE trdtra ! trends manager: tracers 27 27 USE prtctl ! Print control 28 28 USE in_out_manager ! I/O manager … … 35 35 PRIVATE 36 36 37 PUBLIC tra_ldf 38 PUBLIC tra_ldf_init 37 PUBLIC tra_ldf ! called by step.F90 38 PUBLIC tra_ldf_init ! called by opa.F90 39 39 ! 40 40 INTEGER :: nldf = 0 ! type of lateral diffusion used defined from ln_traldf_... namlist logicals) … … 112 112 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 113 113 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 114 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_ldf, ztrdt )115 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_ldf, ztrds )114 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 115 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds ) 116 116 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 117 117 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
r3294 r3876 17 17 USE dom_oce ! ocean space and time domain 18 18 USE zdf_oce ! ocean vertical physics 19 USE trd mod_oce ! ocean active tracer trends20 USE trdtra ! ocean active tracer trends19 USE trd_oce ! trends: ocean variables 20 USE trdtra ! trends manager: tracers 21 21 USE eosbn2 ! equation of state (eos routine) 22 22 USE lbclnk ! lateral boundary conditions (or mpp link) … … 54 54 !! 55 55 !! ** Action : - (ta,sa) after the application od the npc scheme 56 !! - save the associated trends ( ttrd,strd) ('key_trdtra')56 !! - save the associated trends (l_trdtra=T) 57 57 !! 58 58 !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371. … … 196 196 ! ! =============== 197 197 ! 198 199 ! Lateral boundary conditions on ( ta, sa ) ( Unchanged sign) 200 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 201 198 202 IF( l_trdtra ) THEN ! save the Non penetrative mixing trends for diagnostic 199 203 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 200 204 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 201 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_npc, ztrdt )202 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_npc, ztrds )205 CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) 206 CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds ) 203 207 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 204 208 ENDIF 205 209 206 ! Lateral boundary conditions on ( ta, sa ) ( Unchanged sign) 207 ! ------------------------------============ 208 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 209 210 211 ! 2. non penetrative convective scheme statistics 212 ! ----------------------------------------------- 210 ! non penetrative convective scheme statistics 213 211 IF( nn_npcp /= 0 .AND. MOD( kt, nn_npcp ) == 0 ) THEN 214 212 IF(lwp) WRITE(numout,*)' kt=',kt, ' number of statically instable', & -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
r3294 r3876 27 27 USE dom_oce ! ocean space and time domain variables 28 28 USE sbc_oce ! surface boundary condition: ocean 29 USE zdf_oce ! ???29 USE zdf_oce ! ocean vertical mixing 30 30 USE domvvl ! variable volume 31 31 USE dynspg_oce ! surface pressure gradient variables 32 32 USE dynhpg ! hydrostatic pressure gradient 33 USE trdmod_oce ! ocean space and time domain variables 34 USE trdtra ! ocean active tracers trends 35 USE phycst 36 USE obc_oce 33 USE trd_oce ! trends: ocean variables 34 USE trdtra ! trends manager: tracers 35 USE traqsr ! penetrative solar radiation (needed for nksr) 36 USE phycst ! physical constant 37 USE ldftra_oce ! lateral physics on tracers 38 USE bdy_oce ! BDY open boundary condition variables 39 USE bdytra ! BDY open boundary condition (bdy_tra routine) 40 USE obc_oce ! open boundary condition variables 37 41 USE obctra ! open boundary condition (obc_tra routine) 38 USE bdy_oce 39 USE bdytra ! open boundary condition (bdy_tra routine) 42 ! 40 43 USE in_out_manager ! I/O manager 41 44 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 42 45 USE prtctl ! Print control 43 USE traqsr ! penetrative solar radiation (needed for nksr) 46 USE wrk_nemo ! Memory allocation 47 USE timing ! Timing 44 48 #if defined key_agrif 45 49 USE agrif_opa_update 46 50 USE agrif_opa_interp 47 51 #endif 48 USE wrk_nemo ! Memory allocation49 USE timing ! Timing50 52 51 53 IMPLICIT NONE … … 132 134 ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 133 135 ztrds(:,:,:) = tsn(:,:,:,jp_sal) 136 IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend 137 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 138 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds ) 139 ENDIF 134 140 ENDIF 135 141 … … 159 165 ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact 160 166 END DO 161 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_atf, ztrdt )162 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_atf, ztrds )167 CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt ) 168 CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds ) 163 169 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 164 170 END IF … … 168 174 & tab3d_2=tsn(:,:,:,jp_sal), clinfo2= ' Sn: ', mask2=tmask ) 169 175 ! 170 ! 171 IF( nn_timing == 1 ) CALL timing_stop('tra_nxt') 176 IF( nn_timing == 1 ) CALL timing_stop('tra_nxt') 172 177 ! 173 178 END SUBROUTINE tra_nxt -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r3680 r3876 13 13 14 14 !!---------------------------------------------------------------------- 15 !! tra_qsr : trend due to the solar radiation penetration16 !! tra_qsr_init : solar radiation penetration initialization15 !! tra_qsr : trend due to the solar radiation penetration 16 !! tra_qsr_init : solar radiation penetration initialization 17 17 !!---------------------------------------------------------------------- 18 USE oce ! ocean dynamics and active tracers 19 USE dom_oce ! ocean space and time domain 20 USE sbc_oce ! surface boundary condition: ocean 21 USE trc_oce ! share SMS/Ocean variables 22 USE trdmod_oce ! ocean variables trends 23 USE trdtra ! ocean active tracers trends 24 USE in_out_manager ! I/O manager 25 USE phycst ! physical constants 26 USE prtctl ! Print control 27 USE iom ! I/O manager 28 USE fldread ! read input fields 29 USE lib_mpp ! MPP library 18 USE oce ! ocean dynamics and active tracers 19 USE dom_oce ! ocean space and time domain 20 USE sbc_oce ! surface boundary condition: ocean 21 USE trc_oce ! share SMS/Ocean variables 22 USE trd_oce ! trends: ocean variables 23 USE trdtra ! trends manager: tracers 24 USE phycst ! physical constants 25 ! 26 USE in_out_manager ! I/O manager 27 USE prtctl ! Print control 28 USE iom ! I/O manager 29 USE fldread ! read input fields 30 USE lib_mpp ! MPP library 30 31 USE wrk_nemo ! Memory Allocation 31 32 USE timing ! Timing 32 33 33 34 34 IMPLICIT NONE … … 48 48 REAL(wp), PUBLIC :: rn_si1 = 23.0_wp !: deepest depth of extinction (water type I) (2 bands) 49 49 50 ! Module variables 51 REAL(wp) :: xsi0r !: inverse of rn_si0 52 REAL(wp) :: xsi1r !: inverse of rn_si1 50 INTEGER , PUBLIC :: nksr !: levels below which the light cannot penetrate ( depth larger than 391 m) 51 52 REAL(wp) :: xsi0r, xsi1r ! inverse of rn_si0 and rn_si1, resp. 53 REAL(wp), DIMENSION(3,61) :: rkrgb ! tabulated attenuation coefficients for RGB absorption 53 54 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read) 54 INTEGER, PUBLIC :: nksr ! levels below which the light cannot penetrate ( depth larger than 391 m)55 REAL(wp), DIMENSION(3,61) :: rkrgb !: tabulated attenuation coefficients for RGB absorption56 55 57 56 !! * Substitutions … … 87 86 !! 88 87 !! ** Action : - update ta with the penetrative solar radiation trend 89 !! - s ave the trend in ttrd ('key_trdtra')88 !! - send the trend to trdtra (l_trdtra=T) 90 89 !! 91 90 !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 92 91 !! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 93 92 !!---------------------------------------------------------------------- 94 !95 93 INTEGER, INTENT(in) :: kt ! ocean time-step 96 94 ! … … 116 114 ENDIF 117 115 118 IF( l_trdtra ) THEN ! Save t a and sa trends116 IF( l_trdtra ) THEN ! Save temperature trend 119 117 CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 120 118 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) … … 141 139 ! Compute now qsr tracer content field 142 140 ! ************************************ 143 144 141 ! ! ============================================== ! 145 142 IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! … … 168 165 IF( nn_chldta == 1 .OR. lk_vvl ) THEN !* Variable Chlorophyll or ocean volume 169 166 ! 170 IF( nn_chldta == 1 ) THEN ! *Variable Chlorophyll167 IF( nn_chldta == 1 ) THEN !- Variable Chlorophyll 171 168 ! 172 169 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step … … 184 181 END DO 185 182 END DO 186 ELSE !Variable ocean volume but constant chrlorophyll187 zchl = 0.05 ! constant chlorophyll183 ELSE !- Variable ocean volume but constant chrlorophyll 184 zchl = 0.05 ! constant chlorophyll 188 185 irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 189 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll186 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 190 187 zekg(:,:) = rkrgb(2,irgb) 191 188 zekr(:,:) = rkrgb(3,irgb) 192 189 ENDIF 193 190 ! 194 zcoef = ( 1. - rn_abs ) / 3.e0 !equi-partition in R-G-B195 ze0(:,:,1) = rn_abs 196 ze1(:,:,1) = zcoef * qsr(:,:)197 ze2(:,:,1) = zcoef * qsr(:,:)198 ze3(:,:,1) = zcoef * qsr(:,:)199 zea(:,:,1) = qsr(:,:)191 zcoef = ( 1. - rn_abs ) / 3.e0 !- equi-partition in R-G-B 192 ze0(:,:,1) = rn_abs * qsr(:,:) 193 ze1(:,:,1) = zcoef * qsr(:,:) 194 ze2(:,:,1) = zcoef * qsr(:,:) 195 ze3(:,:,1) = zcoef * qsr(:,:) 196 zea(:,:,1) = qsr(:,:) 200 197 ! 201 198 DO jk = 2, nksr+1 … … 283 280 IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics 284 281 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 285 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_qsr, ztrdt )282 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 286 283 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 287 284 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
r3764 r3876 18 18 USE dom_oce ! ocean space domain variables 19 19 USE phycst ! physical constant 20 USE sbcmod ! ln_rnf 21 USE sbcrnf ! River runoff 20 22 USE traqsr ! solar radiation penetration 21 USE trdmod_oce ! ocean trends 22 USE trdtra ! ocean trends 23 USE trd_oce ! trends: ocean variables 24 USE trdtra ! trends manager: tracers 25 ! 23 26 USE in_out_manager ! I/O manager 24 27 USE prtctl ! Print control 25 USE sbcrnf ! River runoff 26 USE sbcmod ! ln_rnf 27 USE iom 28 USE iom ! I/O library 28 29 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 30 USE wrk_nemo ! Memory Allocation … … 91 92 !! where emp, the surface freshwater budget (evaporation minus 92 93 !! precipitation minus runoff) given in kg/m2/s is divided 93 !! by rau0 = 10 20kg/m3 (density of sea water) to obtain m/s.94 !! by rau0 = 1035 kg/m3 (density of sea water) to obtain m/s. 94 95 !! Note: even though Fwe does not appear explicitly for 95 96 !! temperature in this routine, the heat carried by the water … … 107 108 !! ** Action : - Update the 1st level of (ta,sa) with the trend associated 108 109 !! with the tracer surface boundary condition 109 !! - s ave the trend it in ttrd ('key_trdtra')110 !! - send trends to trdtra module (l_trdtra=T) 110 111 !!---------------------------------------------------------------------- 111 112 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 124 125 ENDIF 125 126 126 IF( l_trdtra ) 127 IF( l_trdtra ) THEN !* Save ta and sa trends 127 128 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 128 129 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) … … 137 138 138 139 !---------------------------------------- 139 ! EMP, EMPSand QNS effects140 ! EMP, SFX and QNS effects 140 141 !---------------------------------------- 141 142 ! Set before sbc tracer content fields … … 146 147 & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN 147 148 IF(lwp) WRITE(numout,*) ' nit000-1 surface tracer content forcing fields red in the restart file' 148 zfact = 0.5 e0149 zfact = 0.5_wp 149 150 CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) ) ! before heat content sbc trend 150 151 CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) ) ! before salt content sbc trend 151 152 ELSE ! No restart or restart not found: Euler forward time stepping 152 zfact = 1. e0153 sbc_tsc_b(:,:,:) = 0. e0153 zfact = 1._wp 154 sbc_tsc_b(:,:,:) = 0._wp 154 155 ENDIF 155 156 ELSE ! Swap of forcing fields 156 157 ! ! ---------------------- 157 zfact = 0.5 e0158 zfact = 0.5_wp 158 159 sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:) 159 160 ENDIF … … 226 227 ENDIF 227 228 228 IF( l_trdtra ) THEN ! save the horizontal diffusivetrends for further diagnostics229 IF( l_trdtra ) THEN ! save the trends for further diagnostics 229 230 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 230 231 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 231 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_nsr, ztrdt )232 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_nsr, ztrds )232 CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt ) 233 CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds ) 233 234 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 234 235 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
r3294 r3876 19 19 USE sbc_oce ! surface boundary condition: ocean 20 20 USE dynspg_oce 21 22 21 USE trazdf_exp ! vertical diffusion: explicit (tra_zdf_exp routine) 23 22 USE trazdf_imp ! vertical diffusion: implicit (tra_zdf_imp routine) 24 25 23 USE ldftra_oce ! ocean active tracers: lateral physics 26 USE trdmod_oce ! ocean active tracers: lateral physics 27 USE trdtra ! ocean tracers trends 24 USE trd_oce ! trends: ocean variables 25 USE trdtra ! trends manager: tracers 26 ! 28 27 USE in_out_manager ! I/O manager 29 28 USE prtctl ! Print control … … 96 95 ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) ) - ztrds(:,:,jk) 97 96 END DO 98 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_zdf, ztrdt )99 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_zdf, ztrds )97 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 98 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 100 99 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 101 100 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
