- Timestamp:
- 2011-10-24T11:29:46+02:00 (13 years ago)
- Location:
- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA
- Files:
-
- 2 edited
-
eosbn2.F90 (modified) (21 diffs)
-
traldf_iso_grif.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r2715 r2980 3 3 !! *** MODULE eosbn2 *** 4 4 !! Ocean diagnostic variable : equation of state - in situ and potential density 5 !! - Brunt-Vaisala frequency 5 !! - Brunt-Vaisala frequency 6 6 !!============================================================================== 7 7 !! History : OPA ! 1989-03 (O. Marti) Original code … … 27 27 !! eos_insitu_2d : Compute the in situ density for 2d fields 28 28 !! eos_bn2 : Compute the Brunt-Vaisala frequency 29 !! eos_alpbet : calculates the in situ thermal and haline expansion coeff.29 !! eos_alpbet : calculates the in situ thermal/haline expansion ratio 30 30 !! tfreez : Compute the surface freezing temperature 31 31 !! eos_init : set eos parameters (namelist) … … 41 41 PRIVATE 42 42 43 ! !! * Interface 43 ! !! * Interface 44 44 INTERFACE eos 45 45 MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 46 END INTERFACE 46 END INTERFACE 47 47 INTERFACE bn2 48 48 MODULE PROCEDURE eos_bn2 49 END INTERFACE 49 END INTERFACE 50 50 51 51 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules … … 61 61 62 62 REAL(wp), PUBLIC :: ralpbet !: alpha / beta ratio 63 63 64 64 !! * Substitutions 65 65 # include "domzgr_substitute.h90" … … 75 75 !!---------------------------------------------------------------------- 76 76 !! *** ROUTINE eos_insitu *** 77 !! 78 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 77 !! 78 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 79 79 !! potential temperature and salinity using an equation of state 80 80 !! defined through the namelist parameter nn_eos. … … 134 134 !CDIR NOVERRCHK 135 135 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 136 ! 136 ! 137 137 DO jk = 1, jpkm1 138 138 DO jj = 1, jpj … … 199 199 !!---------------------------------------------------------------------- 200 200 !! *** ROUTINE eos_insitu_pot *** 201 !! 201 !! 202 202 !! ** Purpose : Compute the in situ density (ratio rho/rau0) and the 203 203 !! potential volumic mass (Kg/m3) from potential temperature and 204 !! salinity fields using an equation of state defined through the 204 !! salinity fields using an equation of state defined through the 205 205 !! namelist parameter nn_eos. 206 206 !! … … 230 230 !! nn_eos = 2 : linear equation of state function of temperature and 231 231 !! salinity 232 !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 232 !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 233 233 !! = rn_beta * s - rn_alpha * tn - 1. 234 234 !! rhop(t,s) = rho(t,s) … … 265 265 !CDIR NOVERRCHK 266 266 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 267 ! 267 ! 268 268 DO jk = 1, jpkm1 269 269 DO jj = 1, jpj … … 336 336 !! *** ROUTINE eos_insitu_2d *** 337 337 !! 338 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 338 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 339 339 !! potential temperature and salinity using an equation of state 340 340 !! defined through the namelist parameter nn_eos. * 2D field case … … 374 374 ! ! 2 : salinity [psu] 375 375 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 376 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 376 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 377 377 !! 378 378 INTEGER :: ji, jj ! dummy loop indices … … 449 449 DO jj = 1, jpjm1 450 450 DO ji = 1, fs_jpim1 ! vector opt. 451 prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 451 prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 452 452 END DO 453 453 END DO … … 468 468 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the time- 469 469 !! step of the input arguments 470 !! 470 !! 471 471 !! ** Method : 472 472 !! * nn_eos = 0 : UNESCO sea water properties … … 482 482 !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 483 483 !! The use of potential density to compute N^2 introduces e r r o r 484 !! in the sign of N^2 at great depths. We recommand the use of 484 !! in the sign of N^2 at great depths. We recommand the use of 485 485 !! nn_eos = 0, except for academical studies. 486 486 !! Macro-tasked on horizontal slab (jk-loop) … … 497 497 !! 498 498 INTEGER :: ji, jj, jk ! dummy loop indices 499 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 499 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 500 500 #if defined key_zdfddm 501 501 REAL(wp) :: zds ! local scalars … … 504 504 505 505 ! pn2 : interior points only (2=< jk =< jpkm1 ) 506 ! -------------------------- 506 ! -------------------------- 507 507 ! 508 508 SELECT CASE( nn_eos ) … … 542 542 & - 0.121555e-07_wp ) * zh 543 543 ! 544 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 544 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 545 545 & * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 546 546 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) … … 565 565 & - rn_beta * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) ) ) & 566 566 & / fse3w(:,:,jk) * tmask(:,:,jk) 567 END DO 567 END DO 568 568 #if defined key_zdfddm 569 569 DO jk = 2, jpkm1 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 570 570 DO jj = 1, jpj 571 571 DO ji = 1, jpi 572 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 572 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 573 573 IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 574 574 rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds … … 587 587 588 588 589 SUBROUTINE eos_alpbet( pts, palp h, pbeta)590 !!---------------------------------------------------------------------- 591 !! *** ROUTINE ldf_slp_grif***592 !! 593 !! ** Purpose : Calculates the thermal and haline expansion coefficients at T-points.594 !! 595 !! ** Method : calculates alpha and beta at T-points589 SUBROUTINE eos_alpbet( pts, palpbet, beta0 ) 590 !!---------------------------------------------------------------------- 591 !! *** ROUTINE eos_alpbet *** 592 !! 593 !! ** Purpose : Calculates the in situ thermal/haline expansion ratio at T-points 594 !! 595 !! ** Method : calculates alpha / beta ratio at T-points 596 596 !! * nn_eos = 0 : UNESCO sea water properties 597 !! The brunt-vaisala frequency is computed using the polynomial 598 !! polynomial expression of McDougall (1987): 599 !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 600 !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 601 !! computed and used in zdfddm module : 602 !! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 597 !! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial 598 !! polynomial expression of McDougall (1987). 599 !! Scalar beta0 is returned = 1. 603 600 !! * nn_eos = 1 : linear equation of state (temperature only) 604 !! N^2 = grav * rn_alpha * dk[ t ]/e3w 601 !! The ratio is undefined, so we return alpha as palpbet 602 !! Scalar beta0 is returned = 0. 605 603 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 606 !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 607 !! * nn_eos = 3 : Jackett JAOT 2003 ??? 608 !! 609 !! ** Action : - palph, pbeta : thermal and haline expansion coeff. at T-point 610 !!---------------------------------------------------------------------- 611 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 612 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palph, pbeta ! thermal & haline expansion coeff. 613 ! 604 !! The alpha/beta ratio is returned as ralpbet 605 !! Scalar beta0 is returned = 1. 606 !! 607 !! ** Action : - palpbet : thermal/haline expansion ratio at T-points 608 !! : beta0 : 1. or 0. 609 !!---------------------------------------------------------------------- 610 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 611 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpbet ! thermal/haline expansion ratio 612 REAL(wp), INTENT( out) :: beta0 ! set = 1 except with case 1 eos, rho=rho(T) 613 !! 614 614 INTEGER :: ji, jj, jk ! dummy loop indices 615 REAL(wp) :: zt, zs, zh ! local scalars 615 REAL(wp) :: zt, zs, zh ! local scalars 616 616 !!---------------------------------------------------------------------- 617 617 ! … … 624 624 zt = pts(ji,jj,jk,jp_tem) ! potential temperature 625 625 zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35) 626 zh = fsdept(ji,jj,jk) ! depth in meters 627 ! 628 pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & 629 & - 0.301985e-05_wp ) * zt & 630 & + 0.785567e-03_wp & 631 & + ( 0.515032e-08_wp * zs & 632 & + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs & 633 & + ( ( 0.121551e-17_wp * zh & 634 & - 0.602281e-15_wp * zs & 635 & - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh & 636 & + 0.408195e-10_wp * zs & 637 & + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt & 638 & - 0.121555e-07_wp ) * zh 639 ! 640 palph(ji,jj,jk) = - pbeta(ji,jj,jk) * & 641 & ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & 642 & - 0.203814e-03_wp ) * zt & 643 & + 0.170907e-01_wp ) * zt & 644 & + 0.665157e-01_wp & 645 & + ( - 0.678662e-05_wp * zs & 646 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs & 647 & + ( ( - 0.302285e-13_wp * zh & 648 & - 0.251520e-11_wp * zs & 649 & + 0.512857e-12_wp * zt * zt ) * zh & 650 & - 0.164759e-06_wp * zs & 651 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt & 652 & + 0.380374e-04_wp ) * zh) 626 zh = fsdept(ji,jj,jk) ! depth in meters 627 ! 628 palpbet(ji,jj,jk) = & 629 & ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & 630 & - 0.203814e-03_wp ) * zt & 631 & + 0.170907e-01_wp ) * zt & 632 & + 0.665157e-01_wp & 633 & + ( - 0.678662e-05_wp * zs & 634 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs & 635 & + ( ( - 0.302285e-13_wp * zh & 636 & - 0.251520e-11_wp * zs & 637 & + 0.512857e-12_wp * zt * zt ) * zh & 638 & - 0.164759e-06_wp * zs & 639 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt & 640 & + 0.380374e-04_wp ) * zh 653 641 END DO 654 642 END DO 655 643 END DO 656 ! 657 CASE ( 1 ) 658 palph(:,:,:) = - rn_alpha 659 pbeta(:,:,:) = 0._wp 660 ! 661 CASE ( 2 ) 662 palph(:,:,:) = - rn_alpha 663 pbeta(:,:,:) = rn_beta 644 beta0 = 1._wp 645 ! 646 CASE ( 1 ) !== Linear formulation = F( temperature ) ==! 647 palpbet(:,:,:) = rn_alpha 648 beta0 = 0._wp 649 ! 650 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 651 palpbet(:,:,:) = ralpbet 652 beta0 = 1._wp 664 653 ! 665 654 CASE DEFAULT … … 747 736 748 737 !!====================================================================== 749 END MODULE eosbn2 738 END MODULE eosbn2 -
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
r2715 r2980 4 4 !! Ocean tracers: horizontal component of the lateral tracer mixing trend 5 5 !!====================================================================== 6 !! History : 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) 6 !! History : 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) 7 7 !! ! Griffies operator version adapted from traldf_iso.F90 8 8 !!---------------------------------------------------------------------- … … 11 11 !! 'key_ldfslp' slope of the lateral diffusive direction 12 12 !!---------------------------------------------------------------------- 13 !! tra_ldf_iso_grif : update the tracer trend with the horizontal component 14 !! of the Griffies iso-neutral laplacian operator 13 !! tra_ldf_iso_grif : update the tracer trend with the horizontal component 14 !! of the Griffies iso-neutral laplacian operator 15 15 !!---------------------------------------------------------------------- 16 16 USE oce ! ocean dynamics and active tracers … … 34 34 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: psix_eiv, psiy_eiv !: eiv stream function (diag only) 35 35 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: ah_wslp2 !: aeiv*w-slope^2 36 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt ! atypic workspace36 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt3d !: vertical tracer gradient at 2 levels 37 37 38 38 !! * Substitutions … … 53 53 !! *** ROUTINE tra_ldf_iso_grif *** 54 54 !! 55 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 56 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 55 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 56 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 57 57 !! add it to the general trend of tracer equation. 58 58 !! 59 !! ** Method : The horizontal component of the lateral diffusive trends 59 !! ** Method : The horizontal component of the lateral diffusive trends 60 60 !! is provided by a 2nd order operator rotated along neural or geopo- 61 61 !! tential surfaces to which an eddy induced advection can be added … … 67 67 !! 68 68 !! 2nd part : horizontal fluxes of the lateral mixing operator 69 !! ======== 69 !! ======== 70 70 !! zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 71 71 !! - aht e2u*uslp dk[ mi(mk(tb)) ] … … 99 99 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 100 100 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 102 102 REAL(wp) , INTENT(in ) :: pahtb0 ! background diffusion coef 103 103 ! … … 108 108 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 109 109 REAL(wp) :: zcoef0, zbtr ! - - 110 !REAL(wp), POINTER, DIMENSION(:,:,:) :: zdkt ! 2D+1 workspace111 110 ! 112 111 REAL(wp) :: zslope_skew, zslope_iso, zslope2, zbu, zbv … … 121 120 CALL ctl_stop('tra_ldf_iso_grif: requested workspace arrays unavailable.') ; RETURN 122 121 ENDIF 123 ! ARP - line below uses 'bounds re-mapping' which is only defined in 124 ! Fortran 2003 and up. We would be OK if code was written to use 125 ! zdkt(:,:,1:2) instead as then wouldn't need to re-map bounds. 126 ! As it is, we make zdkt a module array and allocate it in _alloc(). 127 !zdkt(1:jpi,1:jpj,0:1) => wrk_3d_9(:,:,1:2) 128 129 IF( kt == nit000 ) THEN 122 123 IF( kt == nit000.AND..NOT.ALLOCATED(ah_wslp2) ) THEN 130 124 IF(lwp) WRITE(numout,*) 131 125 IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype 132 IF(lwp) WRITE(numout,*) ' WARNING: STILL UNDER TEST, NOT RECOMMENDED. USE AT YOUR OWN PERIL'133 126 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 134 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , zdkt (jpi,jpj,0:1), STAT=ierr )127 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , zdkt3d(jpi,jpj,0:1), STAT=ierr ) 135 128 IF( lk_mpp ) CALL mpp_sum ( ierr ) 136 129 IF( ierr > 0 ) CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate arrays') … … 143 136 144 137 !!---------------------------------------------------------------------- 145 !! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv 138 !! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv 146 139 !!---------------------------------------------------------------------- 147 140 148 !!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads141 !!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads 149 142 150 143 ah_wslp2(:,:,:) = 0._wp … … 159 152 DO jj = 1, jpjm1 160 153 DO ji = 1, fs_jpim1 154 ze1ur = 1._wp / e1u(ji,jj) 161 155 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 162 156 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 163 zah = fsahtu(ji,jj,jk) ! 157 zah = fsahtu(ji,jj,jk) ! fsaht(ji+ip,jj,jk) 164 158 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 165 zslope2 = zslope_skew - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp) 159 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 160 ! (do this by *adding* gradient of depth) 161 zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp) 166 162 zslope2 = zslope2 *zslope2 167 163 ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) & 168 164 & + zah * ( zbu * ze3wr / ( e1t(ji+ip,jj) * e2t(ji+ip,jj) ) ) * zslope2 169 165 IF( ln_traldf_gdia ) THEN 170 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew166 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew ! fsaeit(ji+ip,jj,jk)*zslope_skew 171 167 psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp 172 168 ENDIF … … 182 178 DO jj = 1, jpjm1 183 179 DO ji=1,fs_jpim1 180 ze2vr = 1._wp / e2v(ji,jj) 184 181 ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp) 185 182 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 186 zah = fsaht u(ji,jj,jk) !fsaht(ji,jj+jp,jk)183 zah = fsahtv(ji,jj,jk) ! fsaht(ji,jj+jp,jk) 187 184 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 188 zslope2 = zslope_skew - ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp) 185 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 186 ! (do this by *adding* gradient of depth) 187 zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp) 189 188 zslope2 = zslope2 * zslope2 190 189 ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) & 191 190 & + zah * ( zbv * ze3wr / ( e1t(ji,jj+jp) * e2t(ji,jj+jp) ) ) * zslope2 192 191 IF( ln_traldf_gdia ) THEN 193 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew192 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew 194 193 psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp 195 194 ENDIF … … 207 206 zftu(:,:,:) = 0._wp 208 207 zftv(:,:,:) = 0._wp 209 ! 208 ! 210 209 DO jk = 1, jpkm1 !== before lateral T & S gradients at T-level jk ==! 211 210 DO jj = 1, jpjm1 … … 216 215 END DO 217 216 END DO 218 IF( ln_zps ) THEN! partial steps: correction at the last level217 IF( ln_zps.and.l_grad_zps ) THEN ! partial steps: correction at the last level 219 218 # if defined key_vectopt_loop 220 219 DO jj = 1, 1 … … 224 223 DO ji = 1, jpim1 225 224 # endif 226 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 227 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 225 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 226 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 228 227 END DO 229 228 END DO … … 237 236 ! 238 237 ! !== Vertical tracer gradient at level jk and jk+1 239 zdkt (:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)238 zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 240 239 ! 241 ! ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)242 IF( jk == 1 ) THEN ; zdkt (:,:,0) = zdkt(:,:,1)243 ELSE ; zdkt (:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk)240 ! ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1) 241 IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1) 242 ELSE ; zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 244 243 ENDIF 245 244 246 IF( .NOT. l_triad_iso ) THEN 247 triadi = triadi_g 248 triadj = triadj_g 249 ENDIF 250 251 DO ip = 0, 1 !== Horizontal & vertical fluxes 252 DO kp = 0, 1 253 DO jj = 1, jpjm1 254 DO ji = 1, fs_jpim1 255 ze1ur = 1._wp / e1u(ji,jj) 256 zdxt = zdit(ji,jj,jk) * ze1ur 257 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 258 zdzt = zdkt(ji+ip,jj,kp) * ze3wr 259 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 260 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 261 262 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 263 zah = fsahtu(ji,jj,jk) !*umask(ji,jj,jk+kp) !fsaht(ji+ip,jj,jk) ===>> ???? 264 zah_slp = zah * zslope_iso 265 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew 266 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 267 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 245 246 IF (ln_botmix_grif) THEN 247 DO ip = 0, 1 !== Horizontal & vertical fluxes 248 DO kp = 0, 1 249 DO jj = 1, jpjm1 250 DO ji = 1, fs_jpim1 251 ze1ur = 1._wp / e1u(ji,jj) 252 zdxt = zdit(ji,jj,jk) * ze1ur 253 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 254 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 255 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 256 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 257 258 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 259 ! ln_botmix_grif is .T. don't mask zah for bottom half cells 260 zah = fsahtu(ji,jj,jk) !*umask(ji,jj,jk+kp) !fsaht(ji+ip,jj,jk) ===>> ???? 261 zah_slp = zah * zslope_iso 262 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew 263 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 264 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 265 END DO 268 266 END DO 269 267 END DO 270 268 END DO 271 END DO 272 273 DO jp = 0, 1 274 DO kp = 0, 1 275 DO jj = 1, jpjm1 276 DO ji = 1, fs_jpim1 277 ze2vr = 1._wp / e2v(ji,jj) 278 zdyt = zdjt(ji,jj,jk) * ze2vr 279 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 280 zdzt = zdkt(ji,jj+jp,kp) * ze3wr 281 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 282 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 283 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 284 zah = fsahtv(ji,jj,jk) !*vmask(ji,jj,jk+kp) !fsaht(ji,jj+jp,jk) 285 zah_slp = zah * zslope_iso 286 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew 287 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 288 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 269 270 DO jp = 0, 1 271 DO kp = 0, 1 272 DO jj = 1, jpjm1 273 DO ji = 1, fs_jpim1 274 ze2vr = 1._wp / e2v(ji,jj) 275 zdyt = zdjt(ji,jj,jk) * ze2vr 276 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 277 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 278 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 279 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 280 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 281 ! ln_botmix_grif is .T. don't mask zah for bottom half cells 282 zah = fsahtv(ji,jj,jk) !*vmask(ji,jj,jk+kp) ! fsaht(ji,jj+jp,jk) 283 zah_slp = zah * zslope_iso 284 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew 285 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 286 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 287 END DO 289 288 END DO 290 289 END DO 291 290 END DO 292 END DO 293 294 ! !== divergence and add to the general trend ==! 291 ELSE 292 DO ip = 0, 1 !== Horizontal & vertical fluxes 293 DO kp = 0, 1 294 DO jj = 1, jpjm1 295 DO ji = 1, fs_jpim1 296 ze1ur = 1._wp / e1u(ji,jj) 297 zdxt = zdit(ji,jj,jk) * ze1ur 298 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 299 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 300 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 301 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 302 303 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 304 ! ln_botmix_grif is .F. mask zah for bottom half cells 305 zah = fsahtu(ji,jj,jk) * umask(ji,jj,jk+kp) ! fsaht(ji+ip,jj,jk) ===>> ???? 306 zah_slp = zah * zslope_iso 307 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew ! fsaeit(ji+ip,jj,jk)*zslope_skew 308 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 309 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 310 END DO 311 END DO 312 END DO 313 END DO 314 315 DO jp = 0, 1 316 DO kp = 0, 1 317 DO jj = 1, jpjm1 318 DO ji = 1, fs_jpim1 319 ze2vr = 1._wp / e2v(ji,jj) 320 zdyt = zdjt(ji,jj,jk) * ze2vr 321 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 322 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 323 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 324 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 325 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 326 ! ln_botmix_grif is .F. mask zah for bottom half cells 327 zah = fsahtv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! fsaht(ji,jj+jp,jk) 328 zah_slp = zah * zslope_iso 329 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew 330 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 331 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 332 END DO 333 END DO 334 END DO 335 END DO 336 END IF 337 ! !== divergence and add to the general trend ==! 295 338 DO jj = 2 , jpjm1 296 DO ji = fs_2, fs_jpim1 339 DO ji = fs_2, fs_jpim1 ! vector opt. 297 340 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 298 341 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) & … … 303 346 END DO 304 347 ! 305 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend348 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend 306 349 DO jj = 2, jpjm1 307 DO ji = fs_2, fs_jpim1 350 DO ji = fs_2, fs_jpim1 ! vector opt. 308 351 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & 309 352 & / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) … … 312 355 END DO 313 356 ! 314 ! ! "Poleward" diffusive heat or salt transports (T-S case only)357 ! ! "Poleward" diffusive heat or salt transports (T-S case only) 315 358 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 316 359 IF( jn == jp_tem) htr_ldf(:) = ptr_vj( zftv(:,:,:) ) ! 3.3 names … … 320 363 #if defined key_diaar5 321 364 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 322 z2d(:,:) = 0._wp 323 zztmp = rau0 * rcp 365 z2d(:,:) = 0._wp 366 zztmp = rau0 * rcp 324 367 DO jk = 1, jpkm1 325 368 DO jj = 2, jpjm1 326 369 DO ji = fs_2, fs_jpim1 ! vector opt. 327 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 370 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 328 371 END DO 329 372 END DO … … 332 375 CALL lbc_lnk( z2d, 'U', -1. ) 333 376 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 334 z2d(:,:) = 0._wp 377 z2d(:,:) = 0._wp 335 378 DO jk = 1, jpkm1 336 379 DO jj = 2, jpjm1 337 380 DO ji = fs_2, fs_jpim1 ! vector opt. 338 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 381 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 339 382 END DO 340 383 END DO … … 342 385 z2d(:,:) = zztmp * z2d(:,:) 343 386 CALL lbc_lnk( z2d, 'V', -1. ) 344 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction387 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in j-direction 345 388 END IF 346 389 #endif
Note: See TracChangeset
for help on using the changeset viewer.
