Changeset 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r4292 r6225 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 #if defined key_diaar5 31 USE phycst ! physical constants 32 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 33 #endif 34 USE wrk_nemo ! Memory Allocation 35 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 36 33 37 34 IMPLICIT NONE … … 41 38 42 39 !! * Substitutions 43 # include "domzgr_substitute.h90"44 # include "ldftra_substitute.h90"45 40 # include "vectopt_loop_substitute.h90" 46 41 !!---------------------------------------------------------------------- 47 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)42 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 48 43 !! $Id$ 49 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 51 46 CONTAINS 52 47 53 SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pgu, pgv, & 54 & 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 ) 55 51 !!---------------------------------------------------------------------- 56 52 !! *** ROUTINE tra_ldf_iso *** … … 67 63 !! 68 64 !! 1st part : masked horizontal derivative of T ( di[ t ] ) 69 !! ======== 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 70 67 !! 71 68 !! 2nd part : horizontal fluxes of the lateral mixing operator 72 69 !! ======== 73 !! zftu = (aht+ahtb0)e2u*e3u/e1u di[ tb ]74 !! - ahte2u*uslp dk[ mi(mk(tb)) ]75 !! zftv = (aht+ahtb0)e1v*e3v/e2v dj[ tb ]76 !! - 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)) ] 77 74 !! take the horizontal divergence of the fluxes: 78 !! difft = 1/(e1 t*e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] }75 !! difft = 1/(e1e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] } 79 76 !! Add this trend to the general trend (ta,sa): 80 77 !! ta = ta + difft … … 83 80 !! ======== (excluding the vertical flux proportional to dk[t] ) 84 81 !! vertical fluxes associated with the rotated lateral mixing: 85 !! zftw = -aht {e2t*wslpi di[ mi(mk(tb)) ]86 !! +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)) ] } 87 84 !! take the horizontal divergence of the fluxes: 88 !! difft = 1/(e1 t*e2t*e3t) dk[ zftw ]85 !! difft = 1/(e1e2t*e3t) dk[ zftw ] 89 86 !! Add this trend to the general trend (ta,sa): 90 87 !! pta = pta + difft … … 92 89 !! ** Action : Update pta arrays with the before rotated diffusion 93 90 !!---------------------------------------------------------------------- 94 USE oce , ONLY: zftu => ua , zftv => va ! (ua,va) used as workspace95 !96 91 INTEGER , INTENT(in ) :: kt ! ocean time-step index 97 INTEGER , INTENT(in ) :: kit000 92 INTEGER , INTENT(in ) :: kit000 ! first time step index 98 93 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 99 94 INTEGER , INTENT(in ) :: kjpt ! number of tracers 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] 100 97 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! 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 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 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 107 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 108 REAL(wp) :: zcoef0, zbtr, ztra ! - - 104 INTEGER :: ikt 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 109 #if defined key_diaar5 110 REAL(wp) :: zztmp! local scalar110 REAL(wp) :: zztmp ! local scalar 111 111 #endif 112 REAL(wp), POINTER, DIMENSION(:,: ) ::zdkt, zdk1t, z2d113 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, ztfw112 REAL(wp), POINTER, DIMENSION(:,:) :: zdkt, zdk1t, z2d 113 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, zftu, zftv, ztfw 114 114 !!---------------------------------------------------------------------- 115 115 ! 116 116 IF( nn_timing == 1 ) CALL timing_start('tra_ldf_iso') 117 117 ! 118 CALL wrk_alloc( jpi, jpj, zdkt, zdk1t, z2d ) 119 CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw ) 120 ! 121 118 CALL wrk_alloc( jpi,jpj, zdkt, zdk1t, z2d ) 119 CALL wrk_alloc( jpi,jpj,jpk, zdit, zdjt , zftu, zftv, ztfw ) 120 ! 122 121 IF( kt == kit000 ) THEN 123 122 IF(lwp) WRITE(numout,*) 124 123 IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype 125 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 126 202 ENDIF 127 203 ! … … 133 209 !! I - masked horizontal derivative 134 210 !!---------------------------------------------------------------------- 135 !!bug ajout.... why? (1,jpj,:) and (jpi,1,:) should be sufficient....136 zdit (1,:,:) = 0. e0 ; zdit (jpi,:,:) = 0.e0137 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 138 214 !!end 139 215 … … 147 223 END DO 148 224 END DO 149 IF( ln_zps ) THEN ! partial steps correction at the last ocean level150 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 228 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 153 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 154 END DO 155 END DO 229 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 230 END DO 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 156 240 ENDIF 157 241 ! 158 242 !!---------------------------------------------------------------------- 159 243 !! II - horizontal trend (full) 160 244 !!---------------------------------------------------------------------- 161 !CDIR PARALLEL DO PRIVATE( zdk1t ) 162 ! ! =============== 245 ! 163 246 DO jk = 1, jpkm1 ! Horizontal slab 164 ! ! =============== 165 ! 1. Vertical tracer gradient at level jk and jk+1 166 ! ------------------------------------------------ 167 ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 168 zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 169 ! 170 IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) 171 ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 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) 172 253 ENDIF 173 174 ! 2. Horizontal fluxes 175 ! -------------------- 176 DO jj = 1 , jpjm1 254 DO jj = 1 , jpjm1 !== Horizontal fluxes 177 255 DO ji = 1, fs_jpim1 ! vector opt. 178 zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)179 zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)180 ! 181 zmsku = 1. / MAX( tmask(ji+1,jj,jk ) + tmask(ji,jj,jk+1) &182 & + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk ), 1. )183 ! 184 zmskv = 1. / MAX( tmask(ji,jj+1,jk ) + tmask(ji,jj,jk+1) &185 & + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk ), 1. )186 ! 187 zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku188 zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv256 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 189 267 ! 190 268 zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 191 & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) &192 & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * 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) 193 271 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 194 & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & 195 & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) 196 END DO 197 END DO 198 199 ! II.4 Second derivative (divergence) and add to the general trend 200 ! ---------------------------------------------------------------- 201 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 202 278 DO ji = fs_2, fs_jpim1 ! vector opt. 203 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 204 ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) 205 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 206 END DO 207 END DO 208 ! ! =============== 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 209 284 END DO ! End of slab 210 ! ! =============== 211 ! 212 ! "Poleward" diffusive heat or salt transports (T-S case only) 213 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 214 ! note sign is reversed to give down-gradient diffusive transports (#1043) 215 IF( jn == jp_tem) htr_ldf(:) = ptr_vj( -zftv(:,:,:) ) 216 IF( jn == jp_sal) str_ldf(:) = ptr_vj( -zftv(:,:,:) ) 217 ENDIF 218 219 #if defined key_diaar5 220 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 221 z2d(:,:) = 0._wp 222 ! note sign is reversed to give down-gradient diffusive transports (#1043) 223 zztmp = -1.0_wp * rau0 * rcp 224 DO jk = 1, jpkm1 225 DO jj = 2, jpjm1 226 DO ji = fs_2, fs_jpim1 ! vector opt. 227 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 228 END DO 229 END DO 230 END DO 231 z2d(:,:) = zztmp * z2d(:,:) 232 CALL lbc_lnk( z2d, 'U', -1. ) 233 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 234 z2d(:,:) = 0._wp 235 DO jk = 1, jpkm1 236 DO jj = 2, jpjm1 237 DO ji = fs_2, fs_jpim1 ! vector opt. 238 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 239 END DO 240 END DO 241 END DO 242 z2d(:,:) = zztmp * z2d(:,:) 243 CALL lbc_lnk( z2d, 'V', -1. ) 244 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 245 END IF 246 #endif 247 248 !!---------------------------------------------------------------------- 249 !! III - vertical trend of T & S (extra diagonal terms only) 250 !!---------------------------------------------------------------------- 251 252 ! Local constant initialization 253 ! ----------------------------- 254 ztfw(1,:,:) = 0.e0 ; ztfw(jpi,:,:) = 0.e0 255 285 286 !!---------------------------------------------------------------------- 287 !! III - vertical trend (full) 288 !!---------------------------------------------------------------------- 289 ! 290 ztfw(1,:,:) = 0._wp ; ztfw(jpi,:,:) = 0._wp 291 ! 256 292 ! Vertical fluxes 257 293 ! --------------- 294 ! ! Surface and bottom vertical fluxes set to zero 295 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 258 296 259 ! Surface and bottom vertical fluxes set to zero 260 ztfw(:,:, 1 ) = 0.e0 ; ztfw(:,:,jpk) = 0.e0 261 262 ! interior (2=<jk=<jpk-1) 263 DO jk = 2, jpkm1 297 DO jk = 2, jpkm1 ! interior (2=<jk=<jpk-1) 264 298 DO jj = 2, jpjm1 265 299 DO ji = fs_2, fs_jpim1 ! vector opt. 266 zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk) 267 ! 268 zmsku = 1./MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 269 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk), 1. ) 270 zmskv = 1./MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 271 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk), 1. ) 272 ! 273 zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) 274 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) 275 313 ! 276 314 ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & … … 281 319 END DO 282 320 END DO 283 284 285 ! I.5 Divergence of vertical fluxes added to the general tracer trend 286 ! ------------------------------------------------------------------- 287 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 ==! 288 359 DO jj = 2, jpjm1 289 360 DO ji = fs_2, fs_jpim1 ! vector opt. 290 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 291 ztra = ( ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1) ) * zbtr 292 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) 293 363 END DO 294 364 END DO 295 365 END DO 296 366 ! 297 END DO 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 ! ! =============== 298 414 ! 299 415 CALL wrk_dealloc( jpi, jpj, zdkt, zdk1t, z2d ) 300 CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , ztfw )416 CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , zftu, zftv, ztfw ) 301 417 ! 302 418 IF( nn_timing == 1 ) CALL timing_stop('tra_ldf_iso') 303 419 ! 304 420 END SUBROUTINE tra_ldf_iso 305 306 #else307 !!----------------------------------------------------------------------308 !! default option : Dummy code NO rotation of the diffusive tensor309 !!----------------------------------------------------------------------310 CONTAINS311 SUBROUTINE tra_ldf_iso( kt, kit000,cdtype, pgu, pgv, ptb, pta, kjpt, pahtb0 ) ! Empty routine312 INTEGER:: kt, kit000313 CHARACTER(len=3) :: cdtype314 REAL, DIMENSION(:,:,:) :: pgu, pgv ! tracer gradient at pstep levels315 REAL, DIMENSION(:,:,:,:) :: ptb, pta316 WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt, kit000, cdtype, &317 & pgu(1,1,1), pgv(1,1,1), ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0318 END SUBROUTINE tra_ldf_iso319 #endif320 421 321 422 !!==============================================================================
Note: See TracChangeset
for help on using the changeset viewer.