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