- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r5149 r6808 4 4 !! Ocean tracers: horizontal component of the lateral tracer mixing trend 5 5 !!====================================================================== 6 !! History : OPA ! 1994-08 (G. Madec, M. Imbard) 7 !! 8.0 ! 1997-05 (G. Madec) split into traldf and trazdf 8 !! NEMO ! 2002-08 (G. Madec) Free form, F90 9 !! 1.0 ! 2005-11 (G. Madec) merge traldf and trazdf :-) 10 !! 3.3 ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC 6 !! History : OPA ! 1994-08 (G. Madec, M. Imbard) 7 !! 8.0 ! 1997-05 (G. Madec) split into traldf and trazdf 8 !! NEMO ! 2002-08 (G. Madec) Free form, F90 9 !! 1.0 ! 2005-11 (G. Madec) merge traldf and trazdf :-) 10 !! 3.3 ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC 11 !! 3.7 ! 2014-01 (G. Madec, S. Masson) restructuration/simplification of aht/aeiv specification 12 !! - ! 2014-02 (F. Lemarie, G. Madec) triad operator (Griffies) + Method of Stabilizing Correction 11 13 !!---------------------------------------------------------------------- 12 #if defined key_ldfslp || defined key_esopa 14 13 15 !!---------------------------------------------------------------------- 14 !! 'key_ldfslp' slope of the lateral diffusive direction 16 !! tra_ldf_iso : update the tracer trend with the horizontal component of a iso-neutral laplacian operator 17 !! and with the vertical part of the isopycnal or geopotential s-coord. operator 15 18 !!---------------------------------------------------------------------- 16 !! tra_ldf_iso : update the tracer trend with the horizontal 17 !! component of a iso-neutral laplacian operator 18 !! and with the vertical part of 19 !! the isopycnal or geopotential s-coord. operator 20 !!---------------------------------------------------------------------- 21 USE oce ! ocean dynamics and active tracers 22 USE dom_oce ! ocean space and time domain 23 USE trc_oce ! share passive tracers/Ocean variables 24 USE zdf_oce ! ocean vertical physics 25 USE ldftra_oce ! ocean active tracers: lateral physics 26 USE ldfslp ! iso-neutral slopes 27 USE diaptr ! poleward transport diagnostics 28 USE in_out_manager ! I/O manager 29 USE iom ! I/O library 30 USE phycst ! physical constants 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 32 USE wrk_nemo ! Memory Allocation 33 USE timing ! Timing 19 USE oce ! ocean dynamics and active tracers 20 USE dom_oce ! ocean space and time domain 21 USE trc_oce ! share passive tracers/Ocean variables 22 USE zdf_oce ! ocean vertical physics 23 USE ldftra ! lateral diffusion: tracer eddy coefficients 24 USE ldfslp ! iso-neutral slopes 25 USE diaptr ! poleward transport diagnostics 26 ! 27 USE in_out_manager ! I/O manager 28 USE iom ! I/O library 29 USE phycst ! physical constants 30 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 31 USE wrk_nemo ! Memory Allocation 32 USE timing ! Timing 34 33 35 34 IMPLICIT NONE … … 39 38 40 39 !! * Substitutions 41 # include "domzgr_substitute.h90"42 # include "ldftra_substitute.h90"43 40 # include "vectopt_loop_substitute.h90" 44 41 !!---------------------------------------------------------------------- 45 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)42 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 46 43 !! $Id$ 47 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 49 46 CONTAINS 50 47 51 SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pgu, pgv,&52 & pgui, pgvi,&53 & ptb, pta, kjpt, pahtb0)48 SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv , & 49 & pgui, pgvi, & 50 & ptb , ptbb, pta , kjpt, kpass ) 54 51 !!---------------------------------------------------------------------- 55 52 !! *** ROUTINE tra_ldf_iso *** … … 66 63 !! 67 64 !! 1st part : masked horizontal derivative of T ( di[ t ] ) 68 !! ======== with partial cell update if ln_zps=T. 65 !! ======== with partial cell update if ln_zps=T 66 !! with top cell update if ln_isfcav 69 67 !! 70 68 !! 2nd part : horizontal fluxes of the lateral mixing operator 71 69 !! ======== 72 !! zftu = (aht+ahtb0)e2u*e3u/e1u di[ tb ]73 !! - ahte2u*uslp dk[ mi(mk(tb)) ]74 !! zftv = (aht+ahtb0)e1v*e3v/e2v dj[ tb ]75 !! - ahte2u*vslp dk[ mj(mk(tb)) ]70 !! zftu = pahu e2u*e3u/e1u di[ tb ] 71 !! - pahu e2u*uslp dk[ mi(mk(tb)) ] 72 !! zftv = pahv e1v*e3v/e2v dj[ tb ] 73 !! - pahv e2u*vslp dk[ mj(mk(tb)) ] 76 74 !! take the horizontal divergence of the fluxes: 77 !! difft = 1/(e1 t*e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] }75 !! difft = 1/(e1e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] } 78 76 !! Add this trend to the general trend (ta,sa): 79 77 !! ta = ta + difft … … 82 80 !! ======== (excluding the vertical flux proportional to dk[t] ) 83 81 !! vertical fluxes associated with the rotated lateral mixing: 84 !! zftw = -aht {e2t*wslpi di[ mi(mk(tb)) ]85 !! +e1t*wslpj dj[ mj(mk(tb)) ] }82 !! zftw = - { mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ] 83 !! + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ] } 86 84 !! take the horizontal divergence of the fluxes: 87 !! difft = 1/(e1 t*e2t*e3t) dk[ zftw ]85 !! difft = 1/(e1e2t*e3t) dk[ zftw ] 88 86 !! Add this trend to the general trend (ta,sa): 89 87 !! pta = pta + difft … … 91 89 !! ** Action : Update pta arrays with the before rotated diffusion 92 90 !!---------------------------------------------------------------------- 93 USE oce , ONLY: zftu => ua , zftv => va ! (ua,va) used as workspace94 !95 91 INTEGER , INTENT(in ) :: kt ! ocean time-step index 96 INTEGER , INTENT(in ) :: kit000 92 INTEGER , INTENT(in ) :: kit000 ! first time step index 97 93 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 98 94 INTEGER , INTENT(in ) :: kjpt ! number of tracers 99 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels 100 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at pstep levels 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 103 REAL(wp) , INTENT(in ) :: pahtb0 ! background diffusion coef 95 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 96 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 97 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 98 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 99 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! tracer (kpass=1) or laplacian of tracer (kpass=2) 100 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptbb ! tracer (only used in kpass=2) 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 104 102 ! 105 103 INTEGER :: ji, jj, jk, jn ! dummy loop indices 106 104 INTEGER :: ikt 107 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 108 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 109 REAL(wp) :: zcoef0, zbtr, ztra ! - - 110 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdkt, zdk1t, zdit, zdjt, ztfw 105 INTEGER :: ierr ! local integer 106 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars 107 REAL(wp) :: zmskv, zahv_w, zabe2, zcof2, zcoef4 ! - - 108 REAL(wp) :: zcoef0, ze3w_2, zsign, z2dt, z1_2dt ! - - 109 #if defined key_diaar5 110 REAL(wp) :: zztmp ! local scalar 111 #endif 112 REAL(wp), POINTER, DIMENSION(:,:) :: zdkt, zdk1t, z2d 113 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, zftu, zftv, ztfw 112 114 !!---------------------------------------------------------------------- 113 115 ! 114 116 IF( nn_timing == 1 ) CALL timing_start('tra_ldf_iso') 115 117 ! 116 CALL wrk_alloc( jpi, jpj, z2d ) 117 CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 118 ! 119 118 CALL wrk_alloc( jpi,jpj, zdkt, zdk1t, z2d ) 119 CALL wrk_alloc( jpi,jpj,jpk, zdit, zdjt , zftu, zftv, ztfw ) 120 ! 120 121 IF( kt == kit000 ) THEN 121 122 IF(lwp) WRITE(numout,*) 122 123 IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype 123 124 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 125 ! 126 akz (:,:,:) = 0._wp 127 ah_wslp2(:,:,:) = 0._wp 128 ENDIF 129 ! ! set time step size (Euler/Leapfrog) 130 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt ! at nit000 (Euler) 131 ELSE ; z2dt = 2.* rdt ! (Leapfrog) 132 ENDIF 133 z1_2dt = 1._wp / z2dt 134 ! 135 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) 136 ELSE ; zsign = -1._wp 137 ENDIF 138 139 !!---------------------------------------------------------------------- 140 !! 0 - calculate ah_wslp2 and akz 141 !!---------------------------------------------------------------------- 142 ! 143 IF( kpass == 1 ) THEN !== first pass only ==! 144 ! 145 DO jk = 2, jpkm1 146 DO jj = 2, jpjm1 147 DO ji = fs_2, fs_jpim1 ! vector opt. 148 ! 149 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 150 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 151 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 152 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 153 ! 154 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 155 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 156 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 157 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 158 ! 159 ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 160 & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 161 END DO 162 END DO 163 END DO 164 ! 165 IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient 166 DO jk = 2, jpkm1 167 DO jj = 2, jpjm1 168 DO ji = fs_2, fs_jpim1 169 akz(ji,jj,jk) = 0.25_wp * ( & 170 & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & 171 & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & 172 & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & 173 & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) 174 END DO 175 END DO 176 END DO 177 ! 178 IF( ln_traldf_blp ) THEN ! bilaplacian operator 179 DO jk = 2, jpkm1 180 DO jj = 1, jpjm1 181 DO ji = 1, fs_jpim1 182 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 183 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) 184 END DO 185 END DO 186 END DO 187 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 188 DO jk = 2, jpkm1 189 DO jj = 1, jpjm1 190 DO ji = 1, fs_jpim1 191 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 192 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 193 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 194 END DO 195 END DO 196 END DO 197 ENDIF 198 ! 199 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 200 akz(:,:,:) = ah_wslp2(:,:,:) 201 ENDIF 124 202 ENDIF 125 203 ! … … 131 209 !! I - masked horizontal derivative 132 210 !!---------------------------------------------------------------------- 133 !!bug ajout.... why? (1,jpj,:) and (jpi,1,:) should be sufficient....134 zdit (1,:,:) = 0. e0 ; zdit (jpi,:,:) = 0.e0135 zdjt (1,:,:) = 0. e0 ; zdjt (jpi,:,:) = 0.e0211 !!gm : bug.... why (x,:,:)? (1,jpj,:) and (jpi,1,:) should be sufficient.... 212 zdit (1,:,:) = 0._wp ; zdit (jpi,:,:) = 0._wp 213 zdjt (1,:,:) = 0._wp ; zdjt (jpi,:,:) = 0._wp 136 214 !!end 137 215 … … 145 223 END DO 146 224 END DO 147 148 ! partial cell correction 149 IF( ln_zps ) THEN ! partial steps correction at the last ocean level 150 DO jj = 1, jpjm1 225 IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient 226 DO jj = 1, jpjm1 ! bottom correction (partial bottom cell) 151 227 DO ji = 1, fs_jpim1 ! vector opt. 152 ! IF useless if zpshde defines pgu everywhere153 228 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 154 229 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 155 230 END DO 156 231 END DO 232 IF( ln_isfcav ) THEN ! first wet level beneath a cavity 233 DO jj = 1, jpjm1 234 DO ji = 1, fs_jpim1 ! vector opt. 235 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 236 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 237 END DO 238 END DO 239 ENDIF 157 240 ENDIF 158 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the first wet level beneath a cavity 159 DO jj = 1, jpjm1 241 ! 242 !!---------------------------------------------------------------------- 243 !! II - horizontal trend (full) 244 !!---------------------------------------------------------------------- 245 ! 246 DO jk = 1, jpkm1 ! Horizontal slab 247 ! 248 ! !== Vertical tracer gradient 249 zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1) ! level jk+1 250 ! 251 IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2) 252 ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 253 ENDIF 254 DO jj = 1 , jpjm1 !== Horizontal fluxes 160 255 DO ji = 1, fs_jpim1 ! vector opt. 161 IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 162 IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 163 END DO 164 END DO 165 END IF 166 167 !!---------------------------------------------------------------------- 168 !! II - horizontal trend (full) 169 !!---------------------------------------------------------------------- 170 !!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t ) 171 ! 1. Vertical tracer gradient at level jk and jk+1 172 ! ------------------------------------------------ 173 ! 174 ! interior value 175 DO jk = 2, jpkm1 176 DO jj = 1, jpj 177 DO ji = 1, jpi ! vector opt. 178 zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) 179 ! 180 zdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn ) ) * wmask(ji,jj,jk) 181 END DO 182 END DO 183 END DO 184 ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 185 zdk1t(:,:,1) = ( ptb(:,:,1,jn ) - ptb(:,:,2,jn) ) * wmask(:,:,2) 186 zdkt (:,:,1) = zdk1t(:,:,1) 187 IF ( ln_isfcav ) THEN 188 DO jj = 1, jpj 189 DO ji = 1, jpi ! vector opt. 190 ikt = mikt(ji,jj) ! surface level 191 zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 192 zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 193 END DO 194 END DO 195 END IF 196 197 ! 2. Horizontal fluxes 198 ! -------------------- 199 DO jk = 1, jpkm1 200 DO jj = 1 , jpjm1 201 DO ji = 1, fs_jpim1 ! vector opt. 202 zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 203 zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 204 ! 205 zmsku = 1. / MAX( tmask(ji+1,jj,jk ) + tmask(ji,jj,jk+1) & 206 & + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk ), 1. ) 207 ! 208 zmskv = 1. / MAX( tmask(ji,jj+1,jk ) + tmask(ji,jj,jk+1) & 209 & + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk ), 1. ) 210 ! 211 zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 212 zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 256 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 257 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 258 ! 259 zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & 260 & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) 261 ! 262 zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & 263 & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) 264 ! 265 zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 266 zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 213 267 ! 214 268 zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 215 & + zcof1 * ( zdkt (ji+1,jj,jk) + zdk1t(ji,jj,jk) &216 & + zdk1t(ji+1,jj,jk) + zdkt (ji,jj,jk) ) ) * umask(ji,jj,jk)269 & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & 270 & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) 217 271 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 218 & + zcof2 * ( zdkt (ji,jj+1,jk) + zdk1t(ji,jj,jk) & 219 & + zdk1t(ji,jj+1,jk) + zdkt (ji,jj,jk) ) ) * vmask(ji,jj,jk) 220 END DO 221 END DO 222 223 ! II.4 Second derivative (divergence) and add to the general trend 224 ! ---------------------------------------------------------------- 225 DO jj = 2 , jpjm1 272 & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & 273 & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) 274 END DO 275 END DO 276 ! 277 DO jj = 2 , jpjm1 !== horizontal divergence and add to pta 226 278 DO ji = fs_2, fs_jpim1 ! vector opt. 227 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 228 ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) 229 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 230 END DO 231 END DO 232 ! ! =============== 279 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & 280 & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & 281 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 282 END DO 283 END DO 233 284 END DO ! End of slab 234 ! ! =============== 235 ! 236 ! "Poleward" diffusive heat or salt transports (T-S case only) 237 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 238 ! note sign is reversed to give down-gradient diffusive transports (#1043) 239 IF( jn == jp_tem) htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 240 IF( jn == jp_sal) str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 241 ENDIF 242 243 IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 244 ! 245 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 246 z2d(:,:) = 0._wp 247 DO jk = 1, jpkm1 248 DO jj = 2, jpjm1 249 DO ji = fs_2, fs_jpim1 ! vector opt. 250 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 251 END DO 252 END DO 253 END DO 254 z2d(:,:) = - rau0_rcp * z2d(:,:) ! note sign is reversed to give down-gradient diffusive transports (#1043) 255 CALL lbc_lnk( z2d, 'U', -1. ) 256 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 257 ! 258 z2d(:,:) = 0._wp 259 DO jk = 1, jpkm1 260 DO jj = 2, jpjm1 261 DO ji = fs_2, fs_jpim1 ! vector opt. 262 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 263 END DO 264 END DO 265 END DO 266 z2d(:,:) = - rau0_rcp * z2d(:,:) ! note sign is reversed to give down-gradient diffusive transports (#1043) 267 CALL lbc_lnk( z2d, 'V', -1. ) 268 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 269 END IF 270 ! 271 ENDIF 272 273 !!---------------------------------------------------------------------- 274 !! III - vertical trend of T & S (extra diagonal terms only) 275 !!---------------------------------------------------------------------- 276 277 ! Local constant initialization 278 ! ----------------------------- 279 ztfw(1,:,:) = 0.e0 ; ztfw(jpi,:,:) = 0.e0 280 285 286 !!---------------------------------------------------------------------- 287 !! III - vertical trend (full) 288 !!---------------------------------------------------------------------- 289 ! 290 ztfw(1,:,:) = 0._wp ; ztfw(jpi,:,:) = 0._wp 291 ! 281 292 ! Vertical fluxes 282 293 ! --------------- 294 ! ! Surface and bottom vertical fluxes set to zero 295 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 283 296 284 ! Surface and bottom vertical fluxes set to zero 285 ztfw(:,:, 1 ) = 0.e0 ; ztfw(:,:,jpk) = 0.e0 286 287 ! interior (2=<jk=<jpk-1) 288 DO jk = 2, jpkm1 297 DO jk = 2, jpkm1 ! interior (2=<jk=<jpk-1) 289 298 DO jj = 2, jpjm1 290 299 DO ji = fs_2, fs_jpim1 ! vector opt. 291 zcoef0 = - fsahtw(ji,jj,jk) * wmask(ji,jj,jk) 292 ! 293 zmsku = 1./MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 294 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk), 1. ) 295 zmskv = 1./MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 296 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk), 1. ) 297 ! 298 zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) 299 zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 300 ! 301 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 302 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 303 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 304 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 305 ! 306 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 307 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 308 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 309 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 310 ! 311 zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) !wslpi & j are already w-masked 312 zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 300 313 ! 301 314 ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & … … 306 319 END DO 307 320 END DO 308 309 310 ! I.5 Divergence of vertical fluxes added to the general tracer trend 311 ! ------------------------------------------------------------------- 312 DO jk = 1, jpkm1 321 ! !== add the vertical 33 flux ==! 322 IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz 323 DO jk = 2, jpkm1 324 DO jj = 1, jpjm1 325 DO ji = fs_2, fs_jpim1 326 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & 327 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 328 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 329 END DO 330 END DO 331 END DO 332 ! 333 ELSE ! bilaplacian 334 SELECT CASE( kpass ) 335 CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 336 DO jk = 2, jpkm1 337 DO jj = 1, jpjm1 338 DO ji = fs_2, fs_jpim1 339 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & 340 & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & 341 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 342 END DO 343 END DO 344 END DO 345 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb and ptbb gradients, resp. 346 DO jk = 2, jpkm1 347 DO jj = 1, jpjm1 348 DO ji = fs_2, fs_jpim1 349 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & 350 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 351 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) 352 END DO 353 END DO 354 END DO 355 END SELECT 356 ENDIF 357 ! 358 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pta ==! 313 359 DO jj = 2, jpjm1 314 360 DO ji = fs_2, fs_jpim1 ! vector opt. 315 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 316 ztra = ( ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1) ) * zbtr 317 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 361 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & 362 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 318 363 END DO 319 364 END DO 320 365 END DO 321 366 ! 322 END DO 323 ! 324 CALL wrk_dealloc( jpi, jpj, z2d ) 325 CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 367 IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR. & !== first pass only ( laplacian) ==! 368 ( kpass == 2 .AND. ln_traldf_blp ) ) THEN !== 2nd pass (bilaplacian) ==! 369 ! 370 ! ! "Poleward" diffusive heat or salt transports (T-S case only) 371 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 372 ! note sign is reversed to give down-gradient diffusive transports (#1043) 373 IF( jn == jp_tem) htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 374 IF( jn == jp_sal) str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 375 ENDIF 376 ! 377 IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 378 ! 379 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 380 z2d(:,:) = zftu(ji,jj,1) 381 DO jk = 2, jpkm1 382 DO jj = 2, jpjm1 383 DO ji = fs_2, fs_jpim1 ! vector opt. 384 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 385 END DO 386 END DO 387 END DO 388 !!gm CAUTION I think there is an error of sign when using BLP operator.... 389 !!gm a multiplication by zsign is required (to be checked twice !) 390 z2d(:,:) = - rau0_rcp * z2d(:,:) ! note sign is reversed to give down-gradient diffusive transports (#1043) 391 CALL lbc_lnk( z2d, 'U', -1. ) 392 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 393 ! 394 z2d(:,:) = zftv(ji,jj,1) 395 DO jk = 2, jpkm1 396 DO jj = 2, jpjm1 397 DO ji = fs_2, fs_jpim1 ! vector opt. 398 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 399 END DO 400 END DO 401 END DO 402 z2d(:,:) = - rau0_rcp * z2d(:,:) ! note sign is reversed to give down-gradient diffusive transports (#1043) 403 CALL lbc_lnk( z2d, 'V', -1. ) 404 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 405 END IF 406 ! 407 ENDIF 408 ! 409 ENDIF !== end pass selection ==! 410 ! 411 ! ! =============== 412 END DO ! end tracer loop 413 ! ! =============== 414 ! 415 CALL wrk_dealloc( jpi, jpj, zdkt, zdk1t, z2d ) 416 CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , zftu, zftv, ztfw ) 326 417 ! 327 418 IF( nn_timing == 1 ) CALL timing_stop('tra_ldf_iso') 328 419 ! 329 420 END SUBROUTINE tra_ldf_iso 330 331 #else332 !!----------------------------------------------------------------------333 !! default option : Dummy code NO rotation of the diffusive tensor334 !!----------------------------------------------------------------------335 CONTAINS336 SUBROUTINE tra_ldf_iso( kt, kit000,cdtype, pgu, pgv, pgui, pgvi, ptb, pta, kjpt, pahtb0 ) ! Empty routine337 INTEGER:: kt, kit000338 CHARACTER(len=3) :: cdtype339 REAL, DIMENSION(:,:,:) :: pgu, pgv, pgui, pgvi ! tracer gradient at pstep levels340 REAL, DIMENSION(:,:,:,:) :: ptb, pta341 WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt, kit000, cdtype, &342 & pgu(1,1,1), pgv(1,1,1), ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0343 END SUBROUTINE tra_ldf_iso344 #endif345 421 346 422 !!==============================================================================
Note: See TracChangeset
for help on using the changeset viewer.
