- Timestamp:
- 2015-10-06T13:40:42+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r5758 r5777 2 2 !!====================================================================== 3 3 !! *** MODULE dynldf_iso *** 4 !! Ocean dynamics: lateral viscosity trend4 !! Ocean dynamics: lateral viscosity trend (rotated laplacian operator) 5 5 !!====================================================================== 6 6 !! History : OPA ! 97-07 (G. Madec) Original code … … 8 8 !! - ! 2004-08 (C. Talandier) New trends organization 9 9 !! 2.0 ! 2005-11 (G. Madec) s-coordinate: horizontal diffusion 10 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 11 !! ! add velocity dependent coefficient and optional read in file 10 12 !!---------------------------------------------------------------------- 11 13 … … 17 19 USE oce ! ocean dynamics and tracers 18 20 USE dom_oce ! ocean space and time domain 19 USE ldfdyn _oce ! ocean dynamics lateral physics20 USE ldftra ! lateral physics: eddy diffusivity & EIV coefficients21 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 22 USE ldftra ! lateral physics: eddy diffusivity 21 23 USE zdf_oce ! ocean vertical physics 22 24 USE ldfslp ! iso-neutral slopes 23 25 ! 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link)25 26 USE in_out_manager ! I/O manager 26 27 USE lib_mpp ! MPP library 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 29 USE prtctl ! Print control 28 30 USE wrk_nemo ! Memory Allocation … … 40 42 !! * Substitutions 41 43 # include "domzgr_substitute.h90" 42 # include "ldfdyn_substitute.h90"43 44 # include "vectopt_loop_substitute.h90" 44 45 !!---------------------------------------------------------------------- … … 81 82 !! horizontal fluxes associated with the rotated lateral mixing: 82 83 !! u-component: 83 !! ziut = ( ahmt + ahmb0) e2t * e3t / e1t di[ ub ]84 !! - ahmte2t * mi-1(uslp) dk[ mi(mk(ub)) ]85 !! zjuf = ( ahmf + ahmb0) e1f * e3f / e2f dj[ ub ]86 !! - ahmfe1f * mi(vslp) dk[ mj(mk(ub)) ]84 !! ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t di[ ub ] 85 !! - ahmt e2t * mi-1(uslp) dk[ mi(mk(ub)) ] 86 !! zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f dj[ ub ] 87 !! - ahmf e1f * mi(vslp) dk[ mj(mk(ub)) ] 87 88 !! v-component: 88 !! zivf = ( ahmf + ahmb0) e2t * e3t / e1t di[ vb ]89 !! - ahmfe2t * mj(uslp) dk[ mi(mk(vb)) ]90 !! zjvt = ( ahmt + ahmb0) e1f * e3f / e2f dj[ ub ]91 !! - ahmte1f * mj-1(vslp) dk[ mj(mk(vb)) ]89 !! zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t di[ vb ] 90 !! - ahmf e2t * mj(uslp) dk[ mi(mk(vb)) ] 91 !! zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f dj[ ub ] 92 !! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vb)) ] 92 93 !! take the horizontal divergence of the fluxes: 93 94 !! diffu = 1/(e1u*e2u*e3u) { di [ ziut ] + dj-1[ zjuf ] } … … 108 109 INTEGER :: ji, jj, jk ! dummy loop indices 109 110 REAL(wp) :: zabe1, zabe2, zcof1, zcof2 ! local scalars 110 REAL(wp) :: zmskt, zmskf , zbu, zbv, zuah, zvah! - -111 REAL(wp) :: zmskt, zmskf ! - - 111 112 REAL(wp) :: zcoef0, zcoef3, zcoef4, zmkt, zmkf ! - - 112 113 REAL(wp) :: zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - - … … 127 128 ENDIF 128 129 129 ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum 130 !!gm bug is dyn_ldf_iso called before tra_ldf_iso .... <<<<<===== TO BE CHECKED 131 ! s-coordinate: Iso-level diffusion on momentum but not on tracer 130 132 IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN 131 133 ! … … 133 135 DO jj = 2, jpjm1 134 136 DO ji = 2, jpim1 135 uslp (ji,jj,jk) = - 1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk)) * umask(ji,jj,jk)136 vslp (ji,jj,jk) = - 1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk)) * vmask(ji,jj,jk)137 wslpi(ji,jj,jk) = - 1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk)) * tmask(ji,jj,jk) * 0.5138 wslpj(ji,jj,jk) = - 1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk)) * tmask(ji,jj,jk) * 0.5137 uslp (ji,jj,jk) = - ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 138 vslp (ji,jj,jk) = - ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 139 wslpi(ji,jj,jk) = - ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5 140 wslpj(ji,jj,jk) = - ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5 139 141 END DO 140 142 END DO … … 181 183 DO jj = 2, jpjm1 182 184 DO ji = fs_2, jpi ! vector opt. 183 zabe1 = ( fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) /e1t(ji,jj)184 185 zmskt = 1. /MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1)&186 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1.)185 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) * r1_e1t(ji,jj) 186 187 zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 188 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ) , 1._wp ) 187 189 188 190 zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 189 191 192 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 193 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 194 & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk) 195 END DO 196 END DO 197 ELSE ! other coordinate system (zco or sco) : e3t 198 DO jj = 2, jpjm1 199 DO ji = fs_2, jpi ! vector opt. 200 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * fse3t(ji,jj,jk) * r1_e1t(ji,jj) 201 202 zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) & 203 & + umask(ji-1,jj,jk+1) + umask(ji,jj,jk ) , 1._wp ) 204 205 zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 206 190 207 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 191 208 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & … … 193 210 END DO 194 211 END DO 195 ELSE ! other coordinate system (zco or sco) : e3t196 DO jj = 2, jpjm1197 DO ji = fs_2, jpi ! vector opt.198 zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)199 200 zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &201 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )202 203 zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )204 205 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &206 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &207 & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)208 END DO209 END DO210 212 ENDIF 211 213 … … 213 215 DO jj = 1, jpjm1 214 216 DO ji = 1, fs_jpim1 ! vector opt. 215 zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) /e2f(ji,jj)216 217 zmskf = 1. /MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1)&218 & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1.)217 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * fse3f(ji,jj,jk) * r1_e2f(ji,jj) 218 219 zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & 220 & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ) , 1._wp ) 219 221 220 222 zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) … … 234 236 DO jj = 2, jpjm1 235 237 DO ji = 1, fs_jpim1 ! vector opt. 236 zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) /e1f(ji,jj)237 238 zmskf = 1. /MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1)&239 & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1.)238 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * fse3f(ji,jj,jk) * r1_e1f(ji,jj) 239 240 zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 241 & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ) , 1._wp ) 240 242 241 243 zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 242 244 243 zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) &244 & + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)&245 & + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)245 zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & 246 & + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) & 247 & + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk) 246 248 END DO 247 249 END DO … … 251 253 DO jj = 2, jpj 252 254 DO ji = 1, fs_jpim1 ! vector opt. 253 zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) 255 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) * r1_e2t(ji,jj) 256 257 zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 258 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ) , 1._wp ) 259 260 zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 261 262 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & 263 & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) & 264 & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk) 265 END DO 266 END DO 267 ELSE ! other coordinate system (zco or sco) : e3t 268 DO jj = 2, jpj 269 DO ji = 1, fs_jpim1 ! vector opt. 270 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * fse3t(ji,jj,jk) * r1_e2t(ji,jj) 254 271 255 272 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & … … 263 280 END DO 264 281 END DO 265 ELSE ! other coordinate system (zco or sco) : e3t266 DO jj = 2, jpj267 DO ji = 1, fs_jpim1 ! vector opt.268 zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)269 270 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &271 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )272 273 zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )274 275 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &276 & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &277 & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)278 END DO279 END DO280 282 ENDIF 281 283 … … 283 285 ! Second derivative (divergence) and add to the general trend 284 286 ! ----------------------------------------------------------- 285 286 287 DO jj = 2, jpjm1 287 DO ji = 2, jpim1 !! Question vectop possible??? !!bug 288 ! volume elements 289 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 290 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 291 ! horizontal component of isopycnal momentum diffusive trends 292 zuah =( ziut (ji+1,jj) - ziut (ji,jj ) + & 293 & zjuf (ji ,jj) - zjuf (ji,jj-1) ) / zbu 294 zvah =( zivf (ji,jj ) - zivf (ji-1,jj) + & 295 & zjvt (ji,jj+1) - zjvt (ji,jj ) ) / zbv 296 ! add the trends to the general trends 297 ua (ji,jj,jk) = ua (ji,jj,jk) + zuah 298 va (ji,jj,jk) = va (ji,jj,jk) + zvah 288 DO ji = 2, jpim1 !!gm Question vectop possible??? !!bug 289 ua(ji,jj,jk) = ua(ji,jj,jk) + ( ziut(ji+1,jj) - ziut(ji,jj ) & 290 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 291 va(ji,jj,jk) = va(ji,jj,jk) + ( zivf(ji,jj ) - zivf(ji-1,jj) & 292 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 299 293 END DO 300 294 END DO … … 362 356 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) & 363 357 + tmask(ji,jj,jk )+tmask(ji+1,jj,jk ), 1. ) 364 zmkf = 1./MAX( fmask(ji,jj-1,jk-1) +fmask(ji,jj,jk-1) &365 + fmask(ji,jj-1,jk ) +fmask(ji,jj,jk ), 1. )358 zmkf = 1./MAX( fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1) & 359 + fmask(ji,jj-1,jk ) + fmask(ji,jj,jk ), 1. ) 366 360 367 361 zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi … … 409 403 DO jk = 1, jpkm1 410 404 DO ji = 2, jpim1 411 ! volume elements 412 zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk) 413 zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk) 414 ! part of the k-component of isopycnal momentum diffusive trends 415 zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu 416 zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv 417 ! add the trends to the general trends 418 ua(ji,jj,jk) = ua(ji,jj,jk) + zuav 419 va(ji,jj,jk) = va(ji,jj,jk) + zvav 405 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 406 va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 420 407 END DO 421 408 END DO
Note: See TracChangeset
for help on using the changeset viewer.
