- Timestamp:
- 2010-12-04T16:20:50+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
r2399 r2450 53 53 SUBROUTINE tra_ldf_iso_grif( kt, cdtype, pgu, pgv, & 54 54 & ptb, pta, kjpt, pahtb0 ) 55 !!---------------------------------------------------------------------- 56 !! *** ROUTINE tra_ldf_iso_grif *** 57 !! 58 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 59 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 60 !! add it to the general trend of tracer equation. 61 !! 62 !! ** Method : The horizontal component of the lateral diffusive trends 63 !! is provided by a 2nd order operator rotated along neural or geopo- 64 !! tential surfaces to which an eddy induced advection can be added 65 !! It is computed using before fields (forward in time) and isopyc- 66 !! nal or geopotential slopes computed in routine ldfslp. 67 !! 68 !! 1st part : masked horizontal derivative of T ( di[ t ] ) 69 !! ======== with partial cell update if ln_zps=T. 70 !! 71 !! 2nd part : horizontal fluxes of the lateral mixing operator 72 !! ======== 73 !! zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 74 !! - aht e2u*uslp dk[ mi(mk(tb)) ] 75 !! zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ] 76 !! - aht e2u*vslp dk[ mj(mk(tb)) ] 77 !! take the horizontal divergence of the fluxes: 78 !! difft = 1/(e1t*e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] } 79 !! Add this trend to the general trend (ta,sa): 80 !! ta = ta + difft 81 !! 82 !! 3rd part: vertical trends of the lateral mixing operator 83 !! ======== (excluding the vertical flux proportional to dk[t] ) 84 !! vertical fluxes associated with the rotated lateral mixing: 85 !! zftw =-aht { e2t*wslpi di[ mi(mk(tb)) ] 86 !! + e1t*wslpj dj[ mj(mk(tb)) ] } 87 !! take the horizontal divergence of the fluxes: 88 !! difft = 1/(e1t*e2t*e3t) dk[ zftw ] 89 !! Add this trend to the general trend (ta,sa): 90 !! pta = pta + difft 91 !! 92 !! ** Action : Update pta arrays with the before rotated diffusion 93 !!---------------------------------------------------------------------- 94 USE oce, zftu => ua ! use ua as workspace 95 USE oce, zftv => va ! use va as workspace 96 !! 97 INTEGER , INTENT(in ) :: kt ! ocean time-step index 98 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 99 INTEGER , INTENT(in ) :: kjpt ! number of tracers 100 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 104 !! 105 INTEGER :: ji, jj, jk,jn ! dummy loop indices 106 INTEGER :: ip,jp,kp ! dummy loop indices 107 INTEGER :: iku, ikv ! temporary integer 108 INTEGER :: ierr ! temporary integer 109 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 110 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 111 REAL(wp) :: zcoef0, zbtr ! - - 112 REAL(wp), DIMENSION(jpi,jpj,0:1) :: zdkt ! 2D+1 workspace 113 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, ztfw ! 3D workspace 55 !!---------------------------------------------------------------------- 56 !! *** ROUTINE tra_ldf_iso_grif *** 57 !! 58 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 59 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 60 !! add it to the general trend of tracer equation. 61 !! 62 !! ** Method : The horizontal component of the lateral diffusive trends 63 !! is provided by a 2nd order operator rotated along neural or geopo- 64 !! tential surfaces to which an eddy induced advection can be added 65 !! It is computed using before fields (forward in time) and isopyc- 66 !! nal or geopotential slopes computed in routine ldfslp. 67 !! 68 !! 1st part : masked horizontal derivative of T ( di[ t ] ) 69 !! ======== with partial cell update if ln_zps=T. 70 !! 71 !! 2nd part : horizontal fluxes of the lateral mixing operator 72 !! ======== 73 !! zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 74 !! - aht e2u*uslp dk[ mi(mk(tb)) ] 75 !! zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ] 76 !! - aht e2u*vslp dk[ mj(mk(tb)) ] 77 !! take the horizontal divergence of the fluxes: 78 !! difft = 1/(e1t*e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] } 79 !! Add this trend to the general trend (ta,sa): 80 !! ta = ta + difft 81 !! 82 !! 3rd part: vertical trends of the lateral mixing operator 83 !! ======== (excluding the vertical flux proportional to dk[t] ) 84 !! vertical fluxes associated with the rotated lateral mixing: 85 !! zftw =-aht { e2t*wslpi di[ mi(mk(tb)) ] 86 !! + e1t*wslpj dj[ mj(mk(tb)) ] } 87 !! take the horizontal divergence of the fluxes: 88 !! difft = 1/(e1t*e2t*e3t) dk[ zftw ] 89 !! Add this trend to the general trend (ta,sa): 90 !! pta = pta + difft 91 !! 92 !! ** Action : Update pta arrays with the before rotated diffusion 93 !!---------------------------------------------------------------------- 94 USE oce, zftu => ua ! use ua as workspace 95 USE oce, zftv => va ! use va as workspace 96 !! 97 INTEGER , INTENT(in ) :: kt ! ocean time-step index 98 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 99 INTEGER , INTENT(in ) :: kjpt ! number of tracers 100 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 104 !! 105 INTEGER :: ji, jj, jk,jn ! dummy loop indices 106 INTEGER :: ip,jp,kp ! dummy loop indices 107 INTEGER :: ierr ! temporary integer 108 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 109 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 110 REAL(wp) :: zcoef0, zbtr ! - - 111 REAL(wp), DIMENSION(jpi,jpj,0:1) :: zdkt ! 2D+1 workspace 112 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, ztfw ! 3D workspace 114 113 ! 115 REAL(wp) :: zslope_skew,zslope_iso,zslope2, zbu, zbv116 REAL(wp) :: ze1ur,zdxt,ze2vr,ze3wr,zdyt,zdzt117 REAL(wp) :: ah,zah_slp,zaei_slp114 REAL(wp) :: zslope_skew,zslope_iso,zslope2, zbu, zbv 115 REAL(wp) :: ze1ur,zdxt,ze2vr,ze3wr,zdyt,zdzt 116 REAL(wp) :: ah,zah_slp,zaei_slp 118 117 #if defined key_diaar5 119 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace120 REAL(wp) :: zztmp ! local scalar118 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 119 REAL(wp) :: zztmp ! local scalar 121 120 #endif 122 121 !!---------------------------------------------------------------------- 123 122 124 IF( kt == nit000 ) THEN 125 IF(lwp) WRITE(numout,*) 126 IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype 127 IF(lwp) WRITE(numout,*) ' WARNING: STILL UNDER TEST, NOT RECOMMENDED. USE AT YOUR OWN PERIL' 128 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 129 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , STAT=ierr ) 130 IF( ierr > 0 ) THEN 131 CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator ah_wslp2 ' ) ; RETURN 132 ENDIF 133 IF( ln_traldf_gdia ) THEN 134 ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr ) 135 IF( ierr > 0 ) THEN 136 CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator diagnostics ' ) ; RETURN 137 ENDIF 138 ENDIF 139 ENDIF 140 141 ! 123 IF( kt == nit000 ) THEN 124 IF(lwp) WRITE(numout,*) 125 IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype 126 IF(lwp) WRITE(numout,*) ' WARNING: STILL UNDER TEST, NOT RECOMMENDED. USE AT YOUR OWN PERIL' 127 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 128 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , STAT=ierr ) 129 IF( ierr > 0 ) THEN 130 CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator ah_wslp2 ' ) ; RETURN 131 ENDIF 132 IF( ln_traldf_gdia ) THEN 133 ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr ) 134 IF( ierr > 0 ) THEN 135 CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator diagnostics ' ) ; RETURN 136 ENDIF 137 ENDIF 138 ENDIF 139 142 140 !!---------------------------------------------------------------------- 143 141 !! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv … … 146 144 !!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads 147 145 148 ah_wslp2(:,:,:) = 0.0 149 IF( ln_traldf_gdia ) THEN 150 psix_eiv(:,:,:) = 0.0 151 psiy_eiv(:,:,:) = 0.0 152 ENDIF 153 154 DO ip=0,1 155 DO kp=0,1 156 DO jk=1,jpkm1 157 DO jj=1,jpjm1 158 DO ji=1,fs_jpim1 159 ze3wr=1.0_wp/fse3w(ji+ip,jj,jk+kp) 160 zbu = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 161 ah = fsahtu(ji,jj,jk) ! fsaht(ji+ip,jj,jk) 162 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 163 zslope2=(zslope_skew & 164 & - umask(ji,jj,jk+kp)*(fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk))*ze1ur & 146 ah_wslp2(:,:,:) = 0.0 147 IF( ln_traldf_gdia ) THEN 148 psix_eiv(:,:,:) = 0.0 149 psiy_eiv(:,:,:) = 0.0 150 ENDIF 151 152 DO ip=0,1 153 DO kp=0,1 154 DO jk=1,jpkm1 155 DO jj=1,jpjm1 156 DO ji=1,fs_jpim1 157 ze3wr=1.0_wp/fse3w(ji+ip,jj,jk+kp) 158 zbu = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 159 ah = fsahtu(ji,jj,jk) ! fsaht(ji+ip,jj,jk) 160 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 161 zslope2 = ( zslope_skew & 162 & - umask(ji,jj,jk+kp)*(fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk))*ze1ur )**2 163 ah_wslp2(ji+ip,jj,jk+kp)=ah_wslp2(ji+ip,jj,jk+kp) + & 164 & ah*((zbu*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*zslope2 165 IF( ln_traldf_gdia ) THEN 166 zaei_slp = fsaeiw(ji+ip,jj,jk)*zslope_skew!fsaeit(ji+ip,jj,jk)*zslope_skew 167 psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 168 ENDIF 169 END DO 170 END DO 171 END DO 172 END DO 173 END DO 174 ! 175 DO jp=0,1 176 DO kp=0,1 177 DO jk=1,jpkm1 178 DO jj=1,jpjm1 179 DO ji=1,fs_jpim1 180 ze3wr=1.0_wp/fse3w(ji,jj+jp,jk+kp) 181 zbv = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 182 ah = fsahtu(ji,jj,jk)!fsaht(ji,jj+jp,jk) 183 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 184 zslope2=(zslope_skew - vmask(ji,jj,jk+kp)*(fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk))*ze2vr & 165 185 & )**2 166 ah_wslp2(ji+ip,jj,jk+kp)=ah_wslp2(ji+ip,jj,jk+kp) + & 167 & ah*((zbu*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*zslope2 168 IF( ln_traldf_gdia ) THEN 169 zaei_slp = fsaeiw(ji+ip,jj,jk)*zslope_skew!fsaeit(ji+ip,jj,jk)*zslope_skew 170 psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 171 ENDIF 172 END DO 173 END DO 174 END DO 175 END DO 176 END DO 177 178 DO jp=0,1 179 DO kp=0,1 180 DO jk=1,jpkm1 181 DO jj=1,jpjm1 182 DO ji=1,fs_jpim1 183 ze3wr=1.0_wp/fse3w(ji,jj+jp,jk+kp) 184 zbv = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 185 ah = fsahtu(ji,jj,jk)!fsaht(ji,jj+jp,jk) 186 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 187 zslope2=(zslope_skew - vmask(ji,jj,jk+kp)*(fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk))*ze2vr & 188 & )**2 189 ah_wslp2(ji,jj+jp,jk+kp)=ah_wslp2(ji,jj+jp,jk+kp) + & 186 ah_wslp2(ji,jj+jp,jk+kp)=ah_wslp2(ji,jj+jp,jk+kp) + & 190 187 & ah*((zbv*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*zslope2 191 IF( ln_traldf_gdia ) THEN192 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew193 psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp194 ENDIF195 END DO196 END DO197 END DO198 END DO188 IF( ln_traldf_gdia ) THEN 189 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew 190 psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 191 ENDIF 192 END DO 193 END DO 194 END DO 195 END DO 199 196 END DO 200 201 197 ! 202 198 ! ! =========== … … 224 220 DO ji = 1, jpim1 225 221 # endif 226 iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1 ! last level 227 ikv = MIN( mbathy(ji,jj), mbathy(ji ,jj+1) ) - 1 228 zdit(ji,jj,iku) = pgu(ji,jj,jn) 229 zdjt(ji,jj,ikv) = pgv(ji,jj,jn) 222 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 223 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 230 224 END DO 231 225 END DO … … 294 288 END DO 295 289 296 ! !== divergence and add to the general trend ==!297 DO jj = 2 , jpjm1298 DO ji = fs_2, fs_jpim1 ! vector opt.299 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )300 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) &301 &+ zftv(ji,jj-1,jk) - zftv(ji,jj,jk) )302 END DO303 END DO304 !305 END DO306 !307 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend308 DO jj = 2, jpjm1309 DO ji = fs_2, fs_jpim1 ! vector opt.310 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) &311 &/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )312 END DO313 END DO314 END DO315 !316 ! ! "Poleward" diffusive heat or salt transports (T-S case only)317 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN318 IF( jn == jp_tem) htr_ldf(:) = ptr_vj( zftv(:,:,:) ) ! 3.3 names319 IF( jn == jp_sal) str_ldf(:) = ptr_vj( zftv(:,:,:) )320 ENDIF290 ! !== divergence and add to the general trend ==! 291 DO jj = 2 , jpjm1 292 DO ji = fs_2, fs_jpim1 ! vector opt. 293 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 294 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) & 295 & + zftv(ji,jj-1,jk) - zftv(ji,jj,jk) ) 296 END DO 297 END DO 298 ! 299 END DO 300 ! 301 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend 302 DO jj = 2, jpjm1 303 DO ji = fs_2, fs_jpim1 ! vector opt. 304 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & 305 & / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 306 END DO 307 END DO 308 END DO 309 ! 310 ! ! "Poleward" diffusive heat or salt transports (T-S case only) 311 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 312 IF( jn == jp_tem) htr_ldf(:) = ptr_vj( zftv(:,:,:) ) ! 3.3 names 313 IF( jn == jp_sal) str_ldf(:) = ptr_vj( zftv(:,:,:) ) 314 ENDIF 321 315 322 316 #if defined key_diaar5 323 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 324 zztmp = 0.5 * rau0 * rcp 325 z2d(:,:) = 0.e0 326 DO jk = 1, jpkm1 327 DO jj = 2, jpjm1 328 DO ji = fs_2, fs_jpim1 ! vector opt. 329 z2d(ji,jj) = z2d(ji,jj) + zztmp * zftu(ji,jj,jk) & 330 & * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) * e1u(ji,jj) * fse3u(ji,jj,jk) 331 END DO 332 END DO 333 END DO 334 CALL lbc_lnk( z2d, 'U', -1. ) 335 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 336 z2d(:,:) = 0.e0 337 DO jk = 1, jpkm1 338 DO jj = 2, jpjm1 339 DO ji = fs_2, fs_jpim1 ! vector opt. 340 z2d(ji,jj) = z2d(ji,jj) + zztmp * zftv(ji,jj,jk) & 341 & * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) * e2v(ji,jj) * fse3v(ji,jj,jk) 342 END DO 343 END DO 344 END DO 345 CALL lbc_lnk( z2d, 'V', -1. ) 346 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 347 END IF 317 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 318 z2d(:,:) = 0._wp 319 zztmp = rau0 * rcp 320 DO jk = 1, jpkm1 321 DO jj = 2, jpjm1 322 DO ji = fs_2, fs_jpim1 ! vector opt. 323 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 324 END DO 325 END DO 326 END DO 327 z2d(:,:) = zztmp * z2d(:,:) 328 CALL lbc_lnk( z2d, 'U', -1. ) 329 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 330 DO jk = 1, jpkm1 331 DO jj = 2, jpjm1 332 DO ji = fs_2, fs_jpim1 ! vector opt. 333 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 334 END DO 335 END DO 336 END DO 337 z2d(:,:) = zztmp * z2d(:,:) 338 CALL lbc_lnk( z2d, 'V', -1. ) 339 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 340 END IF 348 341 #endif 349 !350 END DO351 !342 ! 343 END DO 344 ! 352 345 END SUBROUTINE tra_ldf_iso_grif 353 346
Note: See TracChangeset
for help on using the changeset viewer.
