- Timestamp:
- 2016-11-28T17:04:10+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r5836 r7351 14 14 15 15 !!---------------------------------------------------------------------- 16 !! tra_ldf_iso : update the tracer trend with the horizontal component of a iso-neutral laplacian operator17 !! and with the vertical part of the isopycnal or geopotential s-coord. operator16 !! 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 18 18 !!---------------------------------------------------------------------- 19 USE oce 20 USE dom_oce 21 USE trc_oce 22 USE zdf_oce 23 USE ldftra 24 USE ldfslp 25 USE diaptr 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 26 ! 27 USE in_out_manager 28 USE iom 29 USE phycst 30 USE lbclnk 31 USE wrk_nemo 32 USE timing 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 33 33 34 34 IMPLICIT NONE … … 38 38 39 39 !! * Substitutions 40 # include "domzgr_substitute.h90"41 40 # include "vectopt_loop_substitute.h90" 42 41 !!---------------------------------------------------------------------- … … 103 102 ! 104 103 INTEGER :: ji, jj, jk, jn ! dummy loop indices 104 INTEGER :: ikt 105 105 INTEGER :: ierr ! local integer 106 106 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars … … 127 127 ah_wslp2(:,:,:) = 0._wp 128 128 ENDIF 129 !130 129 ! ! set time step size (Euler/Leapfrog) 131 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt tra(1)! at nit000 (Euler)132 ELSE ; z2dt = 2.* rdt tra(1)! (Leapfrog)130 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt ! at nit000 (Euler) 131 ELSE ; z2dt = 2.* rdt ! (Leapfrog) 133 132 ENDIF 134 133 z1_2dt = 1._wp / z2dt … … 138 137 ENDIF 139 138 140 141 139 !!---------------------------------------------------------------------- 142 140 !! 0 - calculate ah_wslp2 and akz … … 149 147 DO ji = fs_2, fs_jpim1 ! vector opt. 150 148 ! 151 zmsku = tmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) &149 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 152 150 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 153 zmskv = tmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) &151 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 154 152 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 155 153 ! … … 183 181 DO ji = 1, fs_jpim1 184 182 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 185 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) ) )183 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) 186 184 END DO 187 185 END DO … … 191 189 DO jj = 1, jpjm1 192 190 DO ji = 1, fs_jpim1 193 ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk)191 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 194 192 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 195 193 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt … … 228 226 DO jj = 1, jpjm1 ! bottom correction (partial bottom cell) 229 227 DO ji = 1, fs_jpim1 ! vector opt. 230 !!gm the following anonymous remark is to considered: ! IF useless if zpshde defines pgu everywhere231 228 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 232 229 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) … … 242 239 ENDIF 243 240 ENDIF 244 241 ! 245 242 !!---------------------------------------------------------------------- 246 243 !! II - horizontal trend (full) … … 255 252 ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 256 253 ENDIF 257 !!gm I don't understand why we should need this.... since wmask is used instead of tmask258 ! IF ( ln_isfcav ) THEN259 ! DO jj = 1, jpj260 ! DO ji = 1, jpi ! vector opt.261 ! ikt = mikt(ji,jj) ! surface level262 ! zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1)263 ! zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt)264 ! END DO265 ! END DO266 ! END IF267 !!gm268 269 254 DO jj = 1 , jpjm1 !== Horizontal fluxes 270 255 DO ji = 1, fs_jpim1 ! vector opt. 271 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)272 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)273 ! 274 zmsku = 1. / MAX( tmask(ji+1,jj,jk ) + tmask(ji,jj,jk+1) &275 & + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk ), 1. )276 ! 277 zmskv = 1. / MAX( tmask(ji,jj+1,jk ) + tmask(ji,jj,jk+1) &278 & + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk ), 1. )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. ) 279 264 ! 280 265 zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku … … 294 279 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & 295 280 & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & 296 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))281 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 297 282 END DO 298 283 END DO 299 284 END DO ! End of slab 300 285 301 302 286 !!---------------------------------------------------------------------- 303 287 !! III - vertical trend (full) 304 288 !!---------------------------------------------------------------------- 305 289 ! 306 290 ztfw(1,:,:) = 0._wp ; ztfw(jpi,:,:) = 0._wp 307 291 ! 308 292 ! Vertical fluxes 309 293 ! --------------- 310 311 ! Surface and bottom vertical fluxes set to zero 294 ! ! Surface and bottom vertical fluxes set to zero 312 295 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 313 296 314 ! interior (2=<jk=<jpk-1) 315 DO jk = 2, jpkm1 297 DO jk = 2, jpkm1 ! interior (2=<jk=<jpk-1) 316 298 DO jj = 2, jpjm1 317 299 DO ji = fs_2, fs_jpim1 ! vector opt. 318 300 ! 319 zmsku = tmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) &301 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 320 302 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 321 zmskv = tmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) &303 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 322 304 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 323 305 ! … … 337 319 END DO 338 320 END DO 339 !340 321 ! !== add the vertical 33 flux ==! 341 322 IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz … … 343 324 DO jj = 1, jpjm1 344 325 DO ji = fs_2, fs_jpim1 345 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &326 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & 346 327 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 347 328 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) … … 358 339 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & 359 340 & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & 360 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) / fse3w(ji,jj,jk)341 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 361 342 END DO 362 343 END DO … … 366 347 DO jj = 1, jpjm1 367 348 DO ji = fs_2, fs_jpim1 368 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &349 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & 369 350 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 370 351 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) … … 379 360 DO ji = fs_2, fs_jpim1 ! vector opt. 380 361 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & 381 & / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))362 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 382 363 END DO 383 364 END DO
Note: See TracChangeset
for help on using the changeset viewer.
