Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
- Timestamp:
- 2015-12-21T12:35:23+01:00 (8 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90 (modified) (38 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
r5930 r6140 39 39 40 40 !! * Substitutions 41 # include "domzgr_substitute.h90"42 41 # include "vectopt_loop_substitute.h90" 43 42 !!---------------------------------------------------------------------- … … 53 52 !! *** ROUTINE tra_adv_fct *** 54 53 !! 55 !! ** Purpose : Compute the now trend due to total advection of 56 !! tracersand add it to the general trend of tracer equations54 !! ** Purpose : Compute the now trend due to total advection of tracers 55 !! and add it to the general trend of tracer equations 57 56 !! 58 57 !! ** Method : - 2nd or 4th FCT scheme on the horizontal direction 59 58 !! (choice through the value of kn_fct) 60 !! - 4th order compact scheme on the vertical59 !! - on the vertical the 4th order is a compact scheme 61 60 !! - corrected flux (monotonic correction) 62 61 !! 63 !! ** Action : - update (pta) with the now advective tracer trends 64 !! - send the trends for further diagnostics 62 !! ** Action : - update pta with the now advective tracer trends 63 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 64 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 65 65 !!---------------------------------------------------------------------- 66 66 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 70 70 INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4) 71 71 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 72 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step72 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 73 73 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 74 74 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 76 76 ! 77 77 INTEGER :: ji, jj, jk, jn ! dummy loop indices 78 REAL(wp) :: z 2dtt, ztra! local scalar78 REAL(wp) :: ztra ! local scalar 79 79 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - - 80 80 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - - … … 101 101 ENDIF 102 102 ! 103 ! 104 IF( lk_vvl ) zwz(:,:, 1 ) = 0._wp ! except at the surface in linear free surface case103 ! ! surface & bottom value : flux set to zero one for all 104 zwz(:,:, 1 ) = 0._wp 105 105 zwx(:,:,jpk) = 0._wp ; zwy(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 106 106 ! 107 107 zwi(:,:,:) = 0._wp 108 ! ! =========== 109 DO jn = 1, kjpt ! tracer loop 110 ! ! =========== 108 ! 109 DO jn = 1, kjpt !== loop over the tracers ==! 111 110 ! 112 111 ! !== upstream advection with initial mass fluxes & intermediate update ==! … … 126 125 END DO 127 126 ! !* upstream tracer flux in the k direction *! 128 DO jk = 2, jpkm1 127 DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask) 129 128 DO jj = 1, jpj 130 129 DO ji = 1, jpi … … 135 134 END DO 136 135 END DO 137 ! 138 IF(.NOT.lk_vvl ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked) 136 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked) 139 137 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 140 138 DO jj = 1, jpj … … 149 147 ! 150 148 DO jk = 1, jpkm1 !* trend and after field with monotonic scheme 151 z2dtt = p2dt(jk)152 149 DO jj = 2, jpjm1 153 150 DO ji = fs_2, fs_jpim1 ! vector opt. … … 155 152 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 156 153 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 157 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))154 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 158 155 ! update and guess with monotonic sheme 159 156 !!gm why tmask added in the two following lines ??? the mask is done in tranxt ! 160 157 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk) 161 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)158 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + p2dt * ztra ) * tmask(ji,jj,jk) 162 159 END DO 163 160 END DO … … 174 171 ENDIF 175 172 ! 176 !177 173 ! !== anti-diffusive flux : high order minus low order ==! 178 174 ! 179 SELECT CASE( kn_fct_h ) 180 ! 181 CASE( 2 ) !2nd order centered175 SELECT CASE( kn_fct_h ) !* horizontal anti-diffusive fluxes 176 ! 177 CASE( 2 ) !- 2nd order centered 182 178 DO jk = 1, jpkm1 183 179 DO jj = 1, jpjm1 … … 189 185 END DO 190 186 ! 191 CASE( 4 ) !4th order centered192 zltu(:,:,jpk) = 0._wp 187 CASE( 4 ) !- 4th order centered 188 zltu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 193 189 zltv(:,:,jpk) = 0._wp 194 DO jk = 1, jpkm1 195 DO jj = 1, jpjm1 ! First derivative (gradient)190 DO jk = 1, jpkm1 ! Laplacian 191 DO jj = 1, jpjm1 ! 1st derivative (gradient) 196 192 DO ji = 1, fs_jpim1 ! vector opt. 197 193 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) … … 199 195 END DO 200 196 END DO 201 DO jj = 2, jpjm1 !197 DO jj = 2, jpjm1 ! 2nd derivative * 1/ 6 202 198 DO ji = fs_2, fs_jpim1 ! vector opt. 203 199 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 … … 206 202 END DO 207 203 END DO 208 !209 204 CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn) 210 205 ! 211 DO jk = 1, jpkm1 206 DO jk = 1, jpkm1 ! Horizontal advective fluxes 212 207 DO jj = 1, jpjm1 213 208 DO ji = 1, fs_jpim1 ! vector opt. … … 221 216 END DO 222 217 ! 223 CASE( 41 ) !4th order centered ==>> !!gm coding attempt need to be tested224 ztu(:,:,jpk) = 0._wp 218 CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested 219 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 225 220 ztv(:,:,jpk) = 0._wp 226 DO jk = 1, jpkm1 ! gradient227 DO jj = 1, jpjm1 ! First derivative (gradient)221 DO jk = 1, jpkm1 ! 1st derivative (gradient) 222 DO jj = 1, jpjm1 228 223 DO ji = 1, fs_jpim1 ! vector opt. 229 224 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) … … 234 229 CALL lbc_lnk( ztu, 'U', -1. ) ; CALL lbc_lnk( ztv, 'V', -1. ) ! Lateral boundary cond. (unchanged sgn) 235 230 ! 236 DO jk = 1, jpkm1 231 DO jk = 1, jpkm1 ! Horizontal advective fluxes 237 232 DO jj = 2, jpjm1 238 233 DO ji = 2, fs_jpim1 ! vector opt. … … 250 245 ! 251 246 END SELECT 252 ! !* vertical anti-diffusive fluxes253 SELECT CASE( kn_fct_v ) ! Interior values (w-masked)254 ! 255 CASE( 2 ) !2nd order centered247 ! 248 SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values) 249 ! 250 CASE( 2 ) !- 2nd order centered 256 251 DO jk = 2, jpkm1 257 252 DO jj = 2, jpjm1 258 253 DO ji = fs_2, fs_jpim1 259 zwz(ji,jj,jk) = ( 0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 260 - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 261 END DO 262 END DO 263 END DO 264 ! 265 CASE( 4 ) ! 4th order COMPACT 266 ! 267 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! COMPACT interpolation of T at w-point 268 ! 254 zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 255 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 256 END DO 257 END DO 258 END DO 259 ! 260 CASE( 4 ) !- 4th order COMPACT 261 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! zwt = COMPACT interpolation of T at w-point 269 262 DO jk = 2, jpkm1 270 263 DO jj = 2, jpjm1 … … 276 269 ! 277 270 END SELECT 278 ! ! top ocean value: high order = upstream ==>> zwz=0 279 zwz(:,:, 1 ) = 0._wp ! only ocean surface as interior zwz values have been w-masked 271 IF( ln_linssh ) THEN ! top ocean value: high order = upstream ==>> zwz=0 272 zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked 273 ENDIF 280 274 ! 281 275 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions 282 276 CALL lbc_lnk( zwz, 'W', 1. ) 283 277 ! 284 278 ! !== monotonicity algorithm ==! 285 279 ! 286 280 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 287 288 281 ! 289 282 ! !== final trend with corrected fluxes ==! 290 283 ! … … 295 288 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 296 289 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 297 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))298 END DO 299 END DO 300 END DO 301 ! 302 IF( l_trd ) THEN 290 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 291 END DO 292 END DO 293 END DO 294 ! 295 IF( l_trd ) THEN ! trend diagnostics (contribution of upstream fluxes) 303 296 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed 304 297 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed … … 311 304 CALL wrk_dealloc( jpi,jpj,jpk, ztrdx, ztrdy, ztrdz ) 312 305 END IF 313 ! 306 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 314 307 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 315 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:)316 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:)308 IF( jn == jp_tem ) htr_adv(:) = htr_adv(:) + ptr_sj( zwy(:,:,:) ) 309 IF( jn == jp_sal ) str_adv(:) = str_adv(:) + ptr_sj( zwy(:,:,:) ) 317 310 ENDIF 318 311 ! 319 END DO 312 END DO ! end of tracer loop 320 313 ! 321 314 CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) … … 348 341 INTEGER , INTENT(in ) :: kjpt ! number of tracers 349 342 INTEGER , INTENT(in ) :: kn_fct_zts ! number of number of vertical sub-timesteps 350 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step343 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 351 344 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 352 345 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 354 347 ! 355 348 REAL(wp), DIMENSION( jpk ) :: zts ! length of sub-timestep for vertical advection 356 REAL(wp) , DIMENSION( jpk ):: zr_p2dt ! reciprocal of tracer timestep349 REAL(wp) :: zr_p2dt ! reciprocal of tracer timestep 357 350 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 358 351 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 359 352 INTEGER :: jtaken ! toggle for collecting appropriate fluxes from sub timesteps 360 353 REAL(wp) :: z_rzts ! Fractional length of Euler forward sub-timestep for vertical advection 361 REAL(wp) :: z 2dtt, ztra! local scalar354 REAL(wp) :: ztra ! local scalar 362 355 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - - 363 356 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - - … … 390 383 zwi(:,:,:) = 0._wp 391 384 z_rzts = 1._wp / REAL( kn_fct_zts, wp ) 392 zr_p2dt(:) = 1._wp / p2dt(:) 385 zr_p2dt = 1._wp / p2dt 386 ! 387 ! surface & Bottom value : flux set to zero for all tracers 388 zwz(:,:, 1 ) = 0._wp 389 zwx(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 390 zwy(:,:,jpk) = 0._wp ; zwi(:,:,jpk) = 0._wp 393 391 ! 394 392 ! ! =========== 395 393 DO jn = 1, kjpt ! tracer loop 396 394 ! ! =========== 397 ! 1. Bottom value : flux set to zero 398 ! ---------------------------------- 399 zwx(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 400 zwy(:,:,jpk) = 0._wp ; zwi(:,:,jpk) = 0._wp 401 402 ! 2. upstream advection with initial mass fluxes & intermediate update 403 ! -------------------------------------------------------------------- 404 ! upstream tracer flux in the i and j direction 405 DO jk = 1, jpkm1 395 ! 396 ! Upstream advection with initial mass fluxes & intermediate update 397 DO jk = 1, jpkm1 ! upstream tracer flux in the i and j direction 406 398 DO jj = 1, jpjm1 407 399 DO ji = 1, fs_jpim1 ! vector opt. … … 416 408 END DO 417 409 END DO 418 419 ! upstream tracer flux in the k direction 420 DO jk = 2, jpkm1 ! Interior value 410 ! ! upstream tracer flux in the k direction 411 DO jk = 2, jpkm1 ! Interior value 421 412 DO jj = 1, jpj 422 413 DO ji = 1, jpi … … 427 418 END DO 428 419 END DO 429 ! ! top value 430 IF( lk_vvl ) THEN ! variable volume: only k=1 as zwz is multiplied by wmask 431 zwz(:,:, 1 ) = 0._wp 432 ELSE ! linear free surface 433 IF( ln_isfcav ) THEN ! ice-shelf cavities 420 IF( ln_linssh ) THEN ! top value : linear free surface case only (as zwz is multiplied by wmask) 421 IF( ln_isfcav ) THEN ! ice-shelf cavities: top value 434 422 DO jj = 1, jpj 435 423 DO ji = 1, jpi … … 437 425 END DO 438 426 END DO 439 ELSE ! standard case427 ELSE ! no cavities, surface value 440 428 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 441 429 ENDIF … … 443 431 ! 444 432 DO jk = 1, jpkm1 ! total advective trend 445 z2dtt = p2dt(jk)446 433 DO jj = 2, jpjm1 447 434 DO ji = fs_2, fs_jpim1 ! vector opt. 448 ! total intermediate advective trends435 ! ! total intermediate advective trends 449 436 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 450 437 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 451 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))452 ! update and guess with monotonic sheme438 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 439 ! ! update and guess with monotonic sheme 453 440 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 454 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)441 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + p2dt * ztra ) * tmask(ji,jj,jk) 455 442 END DO 456 443 END DO … … 497 484 END DO 498 485 END DO 499 486 ! 500 487 ! !* vertical anti-diffusive flux 501 488 zwz_sav(:,:,:) = zwz(:,:,:) 502 489 ztrs (:,:,:,1) = ptb(:,:,:,jn) 503 490 zwzts (:,:,:) = 0._wp 504 IF( lk_vvl ) zwz(:,:, 1 ) = 0._wp ! surface value set to zero in vvl case505 491 ! 506 492 DO jl = 1, kn_fct_zts ! Start of sub timestepping loop … … 508 494 IF( jl == 1 ) THEN ! Euler forward to kick things off 509 495 jtb = 1 ; jtn = 1 ; jta = 2 510 zts(:) = p2dt (:)* z_rzts496 zts(:) = p2dt * z_rzts 511 497 jtaken = MOD( kn_fct_zts + 1 , 2) ! Toggle to collect every second flux 512 498 ! ! starting at jl =1 if kn_fct_zts is odd; … … 514 500 ELSEIF( jl == 2 ) THEN ! First leapfrog step 515 501 jtb = 1 ; jtn = 2 ; jta = 3 516 zts(:) = 2._wp * p2dt (:)* z_rzts502 zts(:) = 2._wp * p2dt * z_rzts 517 503 ELSE ! Shuffle pointers for subsequent leapfrog steps 518 504 jtb = MOD(jtb,3) + 1 … … 528 514 END DO 529 515 END DO 530 IF( .NOT.lk_vvl) THEN ! top value (only in linear free surface case)516 IF( ln_linssh ) THEN ! top value (only in linear free surface case) 531 517 IF( ln_isfcav ) THEN ! ice-shelf cavities 532 518 DO jj = 1, jpj … … 535 521 END DO 536 522 END DO 537 ELSE ! standard case523 ELSE ! no ocean cavities 538 524 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 539 525 ENDIF … … 547 533 ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) & 548 534 & - zts(jk) * ( zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) & 549 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))535 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 550 536 END DO 551 537 END DO … … 557 543 DO jj = 2, jpjm1 558 544 DO ji = fs_2, fs_jpim1 559 zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt (jk)- zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk)545 zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk) 560 546 END DO 561 547 END DO … … 576 562 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 577 563 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 578 & / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))564 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 579 565 END DO 580 566 END DO … … 623 609 !! in-space based differencing for fluid 624 610 !!---------------------------------------------------------------------- 625 REAL(wp) , DIMENSION(jpk) , INTENT(in ) :: p2dt ! vertical profile oftracer time-step611 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 626 612 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field 627 613 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions … … 629 615 INTEGER :: ji, jj, jk ! dummy loop indices 630 616 INTEGER :: ikm1 ! local integer 631 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn , z2dtt! local scalars617 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 632 618 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 633 619 REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo … … 652 638 DO jk = 1, jpkm1 653 639 ikm1 = MAX(jk-1,1) 654 z2dtt = p2dt(jk)655 640 DO jj = 2, jpjm1 656 641 DO ji = fs_2, fs_jpim1 ! vector opt. … … 679 664 680 665 ! up & down beta terms 681 zbt = e1 t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt666 zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 682 667 zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 683 668 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
Note: See TracChangeset
for help on using the changeset viewer.
