Changeset 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
- Timestamp:
- 2012-01-28T17:44:18+01:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
r2715 r3294 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 … … 26 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! Memory Allocation 29 USE timing ! Timing 30 28 31 29 32 IMPLICIT NONE … … 34 37 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: psix_eiv, psiy_eiv !: eiv stream function (diag only) 35 38 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: ah_wslp2 !: aeiv*w-slope^2 36 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt ! atypic workspace39 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt3d !: vertical tracer gradient at 2 levels 37 40 38 41 !! * Substitutions … … 48 51 CONTAINS 49 52 50 SUBROUTINE tra_ldf_iso_grif( kt, cdtype, pgu, pgv, &53 SUBROUTINE tra_ldf_iso_grif( kt, kit000, cdtype, pgu, pgv, & 51 54 & ptb, pta, kjpt, pahtb0 ) 52 55 !!---------------------------------------------------------------------- 53 56 !! *** ROUTINE tra_ldf_iso_grif *** 54 57 !! 55 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 56 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 58 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 59 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 57 60 !! add it to the general trend of tracer equation. 58 61 !! 59 !! ** Method : The horizontal component of the lateral diffusive trends 62 !! ** Method : The horizontal component of the lateral diffusive trends 60 63 !! is provided by a 2nd order operator rotated along neural or geopo- 61 64 !! tential surfaces to which an eddy induced advection can be added … … 67 70 !! 68 71 !! 2nd part : horizontal fluxes of the lateral mixing operator 69 !! ======== 72 !! ======== 70 73 !! zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 71 74 !! - aht e2u*uslp dk[ mi(mk(tb)) ] … … 89 92 !! ** Action : Update pta arrays with the before rotated diffusion 90 93 !!---------------------------------------------------------------------- 91 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released92 94 USE oce , ONLY: zftu => ua , zftv => va ! (ua,va) used as 3D workspace 93 USE wrk_nemo, ONLY: zdit => wrk_3d_6 , zdjt => wrk_3d_7 , ztfw => wrk_3d_8 ! 3D workspace94 USE wrk_nemo, ONLY: z2d => wrk_2d_1 ! 2D workspace95 95 ! 96 96 INTEGER , INTENT(in ) :: kt ! ocean time-step index 97 INTEGER , INTENT(in ) :: kit000 ! first time step index 97 98 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 98 99 INTEGER , INTENT(in ) :: kjpt ! number of tracers 99 100 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 100 101 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 102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 102 103 REAL(wp) , INTENT(in ) :: pahtb0 ! background diffusion coef 103 104 ! … … 108 109 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 109 110 REAL(wp) :: zcoef0, zbtr ! - - 110 !REAL(wp), POINTER, DIMENSION(:,:,:) :: zdkt ! 2D+1 workspace111 111 ! 112 112 REAL(wp) :: zslope_skew, zslope_iso, zslope2, zbu, zbv … … 116 116 REAL(wp) :: zztmp ! local scalar 117 117 #endif 118 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 119 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, ztfw 120 REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d ! 3D workspace 118 121 !!---------------------------------------------------------------------- 119 120 IF( wrk_in_use(3, 6,7,8) .OR. wrk_in_use(2, 1) ) THEN 121 CALL ctl_stop('tra_ldf_iso_grif: requested workspace arrays unavailable.') ; RETURN 122 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( nn_timing == 1 ) CALL timing_start('tra_ldf_iso_grif') 124 ! 125 CALL wrk_alloc( jpi, jpj, z2d ) 126 CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw ) 127 ! 128 129 IF( kt == kit000 .AND. .NOT.ALLOCATED(ah_wslp2) ) THEN 130 130 IF(lwp) WRITE(numout,*) 131 131 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 132 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 134 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , zdkt (jpi,jpj,0:1), STAT=ierr )133 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , zdkt3d(jpi,jpj,0:1), STAT=ierr ) 135 134 IF( lk_mpp ) CALL mpp_sum ( ierr ) 136 135 IF( ierr > 0 ) CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate arrays') 137 136 IF( ln_traldf_gdia ) THEN 138 ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr ) 139 IF( lk_mpp ) CALL mpp_sum ( ierr ) 140 IF( ierr > 0 ) CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate diagnostics') 137 IF (.NOT. ALLOCATED(psix_eiv))THEN 138 ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr ) 139 IF( lk_mpp ) CALL mpp_sum ( ierr ) 140 IF( ierr > 0 ) CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate diagnostics') 141 ENDIF 141 142 ENDIF 142 143 ENDIF 143 144 144 145 !!---------------------------------------------------------------------- 145 !! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv 146 !! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv 146 147 !!---------------------------------------------------------------------- 147 148 148 !!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads149 !!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads 149 150 150 151 ah_wslp2(:,:,:) = 0._wp … … 159 160 DO jj = 1, jpjm1 160 161 DO ji = 1, fs_jpim1 162 ze1ur = 1._wp / e1u(ji,jj) 161 163 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 162 164 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 163 zah = fsahtu(ji,jj,jk) ! 165 zah = fsahtu(ji,jj,jk) ! fsaht(ji+ip,jj,jk) 164 166 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) 167 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 168 ! (do this by *adding* gradient of depth) 169 zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp) 166 170 zslope2 = zslope2 *zslope2 167 171 ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) & 168 172 & + zah * ( zbu * ze3wr / ( e1t(ji+ip,jj) * e2t(ji+ip,jj) ) ) * zslope2 169 173 IF( ln_traldf_gdia ) THEN 170 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew174 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew ! fsaeit(ji+ip,jj,jk)*zslope_skew 171 175 psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp 172 176 ENDIF … … 182 186 DO jj = 1, jpjm1 183 187 DO ji=1,fs_jpim1 188 ze2vr = 1._wp / e2v(ji,jj) 184 189 ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp) 185 190 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)191 zah = fsahtv(ji,jj,jk) ! fsaht(ji,jj+jp,jk) 187 192 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) 193 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 194 ! (do this by *adding* gradient of depth) 195 zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp) 189 196 zslope2 = zslope2 * zslope2 190 197 ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) & 191 198 & + zah * ( zbv * ze3wr / ( e1t(ji,jj+jp) * e2t(ji,jj+jp) ) ) * zslope2 192 199 IF( ln_traldf_gdia ) THEN 193 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew200 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew 194 201 psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp 195 202 ENDIF … … 200 207 END DO 201 208 ! 209 #if defined key_iomput 210 IF( ln_traldf_gdia .AND. cdtype == 'TRA' ) THEN 211 CALL wrk_alloc( jpi , jpj , jpk , zw3d ) 212 DO jk=1,jpkm1 213 zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk) ! u_eiv = -dpsix/dz 214 END DO 215 zw3d(:,:,jpk) = 0._wp 216 CALL iom_put( "uoce_eiv", zw3d ) ! i-eiv current 217 218 DO jk=1,jpk-1 219 zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk) ! v_eiv = -dpsiy/dz 220 END DO 221 zw3d(:,:,jpk) = 0._wp 222 CALL iom_put( "voce_eiv", zw3d ) ! j-eiv current 223 224 DO jk=1,jpk-1 225 DO jj = 2, jpjm1 226 DO ji = fs_2, fs_jpim1 ! vector opt. 227 zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + & 228 & (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx 229 END DO 230 END DO 231 END DO 232 zw3d(:,:,jpk) = 0._wp 233 CALL iom_put( "woce_eiv", zw3d ) ! vert. eiv current 234 CALL wrk_dealloc( jpi , jpj , jpk , zw3d ) 235 ENDIF 236 #endif 202 237 ! ! =========== 203 238 DO jn = 1, kjpt ! tracer loop … … 207 242 zftu(:,:,:) = 0._wp 208 243 zftv(:,:,:) = 0._wp 209 ! 244 ! 210 245 DO jk = 1, jpkm1 !== before lateral T & S gradients at T-level jk ==! 211 246 DO jj = 1, jpjm1 … … 216 251 END DO 217 252 END DO 218 IF( ln_zps ) THEN! partial steps: correction at the last level253 IF( ln_zps.and.l_grad_zps ) THEN ! partial steps: correction at the last level 219 254 # if defined key_vectopt_loop 220 255 DO jj = 1, 1 … … 224 259 DO ji = 1, jpim1 225 260 # endif 226 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 227 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 261 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 262 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 228 263 END DO 229 264 END DO … … 237 272 ! 238 273 ! !== Vertical tracer gradient at level jk and jk+1 239 zdkt (:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)274 zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 240 275 ! 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)276 ! ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1) 277 IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1) 278 ELSE ; zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 244 279 ENDIF 245 280 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 281 282 IF (ln_botmix_grif) THEN 283 DO ip = 0, 1 !== Horizontal & vertical fluxes 284 DO kp = 0, 1 285 DO jj = 1, jpjm1 286 DO ji = 1, fs_jpim1 287 ze1ur = 1._wp / e1u(ji,jj) 288 zdxt = zdit(ji,jj,jk) * ze1ur 289 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 290 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 291 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 292 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 293 294 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 295 ! ln_botmix_grif is .T. don't mask zah for bottom half cells 296 zah = fsahtu(ji,jj,jk) !*umask(ji,jj,jk+kp) !fsaht(ji+ip,jj,jk) ===>> ???? 297 zah_slp = zah * zslope_iso 298 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew 299 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 300 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 301 END DO 268 302 END DO 269 303 END DO 270 304 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 305 306 DO jp = 0, 1 307 DO kp = 0, 1 308 DO jj = 1, jpjm1 309 DO ji = 1, fs_jpim1 310 ze2vr = 1._wp / e2v(ji,jj) 311 zdyt = zdjt(ji,jj,jk) * ze2vr 312 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 313 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 314 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 315 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 316 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 317 ! ln_botmix_grif is .T. don't mask zah for bottom half cells 318 zah = fsahtv(ji,jj,jk) !*vmask(ji,jj,jk+kp) ! fsaht(ji,jj+jp,jk) 319 zah_slp = zah * zslope_iso 320 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew 321 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 322 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 323 END DO 289 324 END DO 290 325 END DO 291 326 END DO 292 END DO 293 294 ! !== divergence and add to the general trend ==! 327 ELSE 328 DO ip = 0, 1 !== Horizontal & vertical fluxes 329 DO kp = 0, 1 330 DO jj = 1, jpjm1 331 DO ji = 1, fs_jpim1 332 ze1ur = 1._wp / e1u(ji,jj) 333 zdxt = zdit(ji,jj,jk) * ze1ur 334 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 335 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 336 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 337 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 338 339 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 340 ! ln_botmix_grif is .F. mask zah for bottom half cells 341 zah = fsahtu(ji,jj,jk) * umask(ji,jj,jk+kp) ! fsaht(ji+ip,jj,jk) ===>> ???? 342 zah_slp = zah * zslope_iso 343 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew ! fsaeit(ji+ip,jj,jk)*zslope_skew 344 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 345 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 346 END DO 347 END DO 348 END DO 349 END DO 350 351 DO jp = 0, 1 352 DO kp = 0, 1 353 DO jj = 1, jpjm1 354 DO ji = 1, fs_jpim1 355 ze2vr = 1._wp / e2v(ji,jj) 356 zdyt = zdjt(ji,jj,jk) * ze2vr 357 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 358 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 359 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 360 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 361 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 362 ! ln_botmix_grif is .F. mask zah for bottom half cells 363 zah = fsahtv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! fsaht(ji,jj+jp,jk) 364 zah_slp = zah * zslope_iso 365 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew 366 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 367 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 368 END DO 369 END DO 370 END DO 371 END DO 372 END IF 373 ! !== divergence and add to the general trend ==! 295 374 DO jj = 2 , jpjm1 296 DO ji = fs_2, fs_jpim1 375 DO ji = fs_2, fs_jpim1 ! vector opt. 297 376 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 298 377 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) & … … 303 382 END DO 304 383 ! 305 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend384 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend 306 385 DO jj = 2, jpjm1 307 DO ji = fs_2, fs_jpim1 386 DO ji = fs_2, fs_jpim1 ! vector opt. 308 387 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & 309 388 & / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) … … 312 391 END DO 313 392 ! 314 ! ! "Poleward" diffusive heat or salt transports (T-S case only)393 ! ! "Poleward" diffusive heat or salt transports (T-S case only) 315 394 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 316 395 IF( jn == jp_tem) htr_ldf(:) = ptr_vj( zftv(:,:,:) ) ! 3.3 names … … 320 399 #if defined key_diaar5 321 400 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 322 z2d(:,:) = 0._wp 323 zztmp = rau0 * rcp 401 z2d(:,:) = 0._wp 402 zztmp = rau0 * rcp 324 403 DO jk = 1, jpkm1 325 404 DO jj = 2, jpjm1 326 405 DO ji = fs_2, fs_jpim1 ! vector opt. 327 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 406 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 328 407 END DO 329 408 END DO … … 332 411 CALL lbc_lnk( z2d, 'U', -1. ) 333 412 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 334 z2d(:,:) = 0._wp 413 z2d(:,:) = 0._wp 335 414 DO jk = 1, jpkm1 336 415 DO jj = 2, jpjm1 337 416 DO ji = fs_2, fs_jpim1 ! vector opt. 338 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 417 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 339 418 END DO 340 419 END DO … … 342 421 z2d(:,:) = zztmp * z2d(:,:) 343 422 CALL lbc_lnk( z2d, 'V', -1. ) 344 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction423 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in j-direction 345 424 END IF 346 425 #endif … … 348 427 END DO 349 428 ! 350 IF( wrk_not_released(3, 6,7,8) .OR. & 351 wrk_not_released(2, 1) ) CALL ctl_stop('tra_ldf_iso_grif: failed to release workspace arrays') 429 CALL wrk_dealloc( jpi, jpj, z2d ) 430 CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw ) 431 ! 432 IF( nn_timing == 1 ) CALL timing_stop('tra_ldf_iso_grif') 352 433 ! 353 434 END SUBROUTINE tra_ldf_iso_grif … … 357 438 !! default option : Dummy code NO rotation of the diffusive tensor 358 439 !!---------------------------------------------------------------------- 440 REAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: psix_eiv, psiy_eiv !: eiv stream function (diag only) 359 441 CONTAINS 360 SUBROUTINE tra_ldf_iso_grif( kt, cdtype, pgu, pgv, ptb, pta, kjpt, pahtb0 ) ! Empty routine 442 SUBROUTINE tra_ldf_iso_grif( kt, kit000, cdtype, pgu, pgv, & 443 & ptb, pta, kjpt, pahtb0 ) 361 444 CHARACTER(len=3) :: cdtype 445 INTEGER :: kit000 ! first time step index 362 446 REAL, DIMENSION(:,:,:) :: pgu, pgv ! tracer gradient at pstep levels 363 447 REAL, DIMENSION(:,:,:,:) :: ptb, pta
Note: See TracChangeset
for help on using the changeset viewer.
