Changeset 5758 for branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
- Timestamp:
- 2015-09-24T08:31:40+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_bilapg.F90
r4990 r5758 1 1 MODULE dynldf_bilapg 2 !!====================================================================== 3 !! *** MODULE dynldf_bilapg *** 4 !! Ocean dynamics: lateral viscosity trend 5 !!====================================================================== 6 !! History : OPA ! 1997-07 (G. Madec) Original code 7 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 8 !! 2.0 ! 2004-08 (C. Talandier) New trends organization 9 !!---------------------------------------------------------------------- 10 #if defined key_ldfslp || defined key_esopa 11 !!---------------------------------------------------------------------- 12 !! 'key_ldfslp' Rotation of mixing tensor 13 !!---------------------------------------------------------------------- 14 !! dyn_ldf_bilapg : update the momentum trend with the horizontal part 15 !! of the horizontal s-coord. bilaplacian diffusion 16 !! ldfguv : 17 !!---------------------------------------------------------------------- 18 USE oce ! ocean dynamics and tracers 19 USE dom_oce ! ocean space and time domain 20 USE ldfdyn_oce ! ocean dynamics lateral physics 21 USE zdf_oce ! ocean vertical physics 22 USE ldfslp ! iso-neutral slopes available 23 USE ldftra_oce, ONLY: ln_traldf_iso 24 ! 25 USE in_out_manager ! I/O manager 26 USE lib_mpp ! MPP library 27 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 28 USE prtctl ! Print control 29 USE wrk_nemo ! Memory Allocation 30 USE timing ! Timing 2 !!============================================================================== 3 4 5 !! ====>>> Empty TO BE REMOVED 6 31 7 32 IMPLICIT NONE 33 PRIVATE 8 !!============================================================================== 34 9 35 PUBLIC dyn_ldf_bilapg ! called by step.F9036 37 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zfvw , zdiu, zdiv ! 2D workspace (ldfguv)38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdju, zdj1u, zdjv, zdj1v ! 2D workspace (ldfguv)39 40 !! * Substitutions41 # include "domzgr_substitute.h90"42 # include "ldfdyn_substitute.h90"43 !!----------------------------------------------------------------------44 !! NEMO/OPA 3.3 , NEMO Consortium (2010)45 !! $Id$46 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)47 !!----------------------------------------------------------------------48 CONTAINS49 50 INTEGER FUNCTION dyn_ldf_bilapg_alloc()51 !!----------------------------------------------------------------------52 !! *** ROUTINE dyn_ldf_bilapg_alloc ***53 !!----------------------------------------------------------------------54 ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) , &55 & zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc )56 !57 IF( dyn_ldf_bilapg_alloc /= 0 ) CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays')58 END FUNCTION dyn_ldf_bilapg_alloc59 60 61 SUBROUTINE dyn_ldf_bilapg( kt )62 !!----------------------------------------------------------------------63 !! *** ROUTINE dyn_ldf_bilapg ***64 !!65 !! ** Purpose : Compute the before trend of the horizontal momentum66 !! diffusion and add it to the general trend of momentum equation.67 !!68 !! ** Method : The lateral momentum diffusive trends is provided by a69 !! a 4th order operator rotated along geopotential surfaces. It is70 !! computed using before fields (forward in time) and geopotential71 !! slopes computed in routine inildf.72 !! -1- compute the geopotential harmonic operator applied to73 !! (ub,vb) and multiply it by the eddy diffusivity coefficient74 !! (done by a call to ldfgpu and ldfgpv routines) The result is in75 !! (zwk1,zwk2) arrays. Applied the domain lateral boundary conditions76 !! by call to lbc_lnk.77 !! -2- applied to (zwk1,zwk2) the geopotential harmonic operator78 !! by a second call to ldfgpu and ldfgpv routines respectively. The79 !! result is in (zwk3,zwk4) arrays.80 !! -3- Add this trend to the general trend (ta,sa):81 !! (ua,va) = (ua,va) + (zwk3,zwk4)82 !!83 !! ** Action : - Update (ua,va) arrays with the before geopotential84 !! biharmonic mixing trend.85 !!----------------------------------------------------------------------86 INTEGER, INTENT( in ) :: kt ! ocean time-step index87 !88 INTEGER :: ji, jj, jk ! dummy loop indices89 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwk1, zwk2, zwk3, zwk490 !!----------------------------------------------------------------------91 !92 IF( nn_timing == 1 ) CALL timing_start('dyn_ldf_bilapg')93 !94 CALL wrk_alloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 )95 !96 IF( kt == nit000 ) THEN97 IF(lwp) WRITE(numout,*)98 IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate'99 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'100 ! ! allocate dyn_ldf_bilapg arrays101 IF( dyn_ldf_bilapg_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays')102 ENDIF103 104 ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum105 IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN106 !107 DO jk = 1, jpk ! set the slopes of iso-level108 DO jj = 2, jpjm1109 DO ji = 2, jpim1110 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)111 vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)112 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.5113 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.5114 END DO115 END DO116 END DO117 ! Lateral boundary conditions on the slopes118 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. )119 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. )120 121 !!bug122 IF( kt == nit000 ) then123 IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)), &124 & ' wi', sqrt(MAXVAL(wslpi)) , ' wj', sqrt(MAXVAL(wslpj))125 endif126 !!end127 ENDIF128 129 zwk1(:,:,:) = 0.e0 ; zwk3(:,:,:) = 0.e0130 zwk2(:,:,:) = 0.e0 ; zwk4(:,:,:) = 0.e0131 132 ! Laplacian of (ub,vb) multiplied by ahm133 ! --------------------------------------134 CALL ldfguv( ub, vb, zwk1, zwk2, 1 ) ! rotated harmonic operator applied to (ub,vb)135 ! ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) )136 CALL lbc_lnk( zwk1, 'U', -1. ) ; CALL lbc_lnk( zwk2, 'V', -1. ) ! Lateral boundary conditions137 138 ! Bilaplacian of (ub,vb)139 ! ----------------------140 CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 ) ! rotated harmonic operator applied to (zwk1,zwk2)141 ! ! (output in (zwk3,zwk4) )142 143 ! Update the momentum trends144 ! --------------------------145 DO jj = 2, jpjm1 ! add the diffusive trend to the general momentum trends146 DO jk = 1, jpkm1147 DO ji = 2, jpim1148 ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk)149 va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk)150 END DO151 END DO152 END DO153 !154 CALL wrk_dealloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 )155 !156 IF( nn_timing == 1 ) CALL timing_stop('dyn_ldf_bilapg')157 !158 END SUBROUTINE dyn_ldf_bilapg159 160 161 SUBROUTINE ldfguv( pu, pv, plu, plv, kahm )162 !!----------------------------------------------------------------------163 !! *** ROUTINE ldfguv ***164 !!165 !! ** Purpose : Apply a geopotential harmonic operator to (pu,pv)166 !! (defined at u- and v-points) and multiply it by the eddy167 !! viscosity coefficient (if kahm=1).168 !!169 !! ** Method : The harmonic operator rotated along geopotential170 !! surfaces is applied to (pu,pv) using the slopes of geopotential171 !! surfaces computed in inildf routine. The result is provided in172 !! (plu,plv) arrays. It is computed in 2 stepv:173 !!174 !! First step: horizontal part of the operator. It is computed on175 !! ========== pu as follows (idem on pv)176 !! horizontal fluxes :177 !! zftu = e2u*e3u/e1u di[ pu ] - e2u*uslp dk[ mi(mk(pu)) ]178 !! zftv = e1v*e3v/e2v dj[ pu ] - e1v*vslp dk[ mj(mk(pu)) ]179 !! take the horizontal divergence of the fluxes (no divided by180 !! the volume element :181 !! plu = di-1[ zftu ] + dj-1[ zftv ]182 !!183 !! Second step: vertical part of the operator. It is computed on184 !! =========== pu as follows (idem on pv)185 !! vertical fluxes :186 !! zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2) dk-1[ pu ]187 !! - e2t * wslpi di[ mi(mk(pu)) ]188 !! - e1t * wslpj dj[ mj(mk(pu)) ]189 !! take the vertical divergence of the fluxes add it to the hori-190 !! zontal component, divide the result by the volume element and191 !! if kahm=1, multiply by the eddy diffusivity coefficient:192 !! plu = aht / (e1t*e2t*e3t) { plu + dk[ zftw ] }193 !! else:194 !! plu = 1 / (e1t*e2t*e3t) { plu + dk[ zftw ] }195 !!196 !! ** Action :197 !! plu, plv : partial harmonic operator applied to198 !! pu and pv (all the components except199 !! second order vertical derivative term)200 !!----------------------------------------------------------------------201 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu , pv ! 1st call: before horizontal velocity202 ! ! 2nd call: ahm x these fields203 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: plu, plv ! partial harmonic operator applied to204 ! ! pu and pv (all the components except205 ! ! second order vertical derivative term)206 INTEGER , INTENT(in ) :: kahm ! =1 1st call ; =2 2nd call207 !208 INTEGER :: ji, jj, jk ! dummy loop indices209 REAL(wp) :: zabe1 , zabe2 , zcof1 , zcof2 ! local scalar210 REAL(wp) :: zcoef0, zcoef3, zcoef4 ! - -211 REAL(wp) :: zbur, zbvr, zmkt, zmkf, zuav, zvav ! - -212 REAL(wp) :: zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - -213 !214 REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v215 !!----------------------------------------------------------------------216 !217 IF( nn_timing == 1 ) CALL timing_start('ldfguv')218 !219 CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )220 !221 ! ! ********** ! ! ===============222 DO jk = 1, jpkm1 ! First step ! ! Horizontal slab223 ! ! ********** ! ! ===============224 225 ! I.1 Vertical gradient of pu and pv at level jk and jk+1226 ! -------------------------------------------------------227 ! surface boundary condition: zdku(jk=1)=zdku(jk=2)228 ! zdkv(jk=1)=zdkv(jk=2)229 230 zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1)231 zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1)232 233 IF( jk == 1 ) THEN234 zdku(:,:) = zdk1u(:,:)235 zdkv(:,:) = zdk1v(:,:)236 ELSE237 zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk)238 zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk)239 ENDIF240 241 ! -----f-----242 ! I.2 Horizontal fluxes on U |243 ! ------------------------=== t u t244 ! |245 ! i-flux at t-point -----f-----246 DO jj = 1, jpjm1247 DO ji = 2, jpi248 zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)249 250 zmkt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &251 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )252 253 zcof1 = -e2t(ji,jj) * zmkt &254 * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )255 256 ziut(ji,jj) = tmask(ji,jj,jk) * &257 ( zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) &258 + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &259 +zdk1u(ji,jj) + zdku (ji-1,jj) ) )260 END DO261 END DO262 263 ! j-flux at f-point264 DO jj = 1, jpjm1265 DO ji = 1, jpim1266 zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)267 268 zmkf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &269 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )270 271 zcof2 = -e1f(ji,jj) * zmkf &272 * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )273 274 zjuf(ji,jj) = fmask(ji,jj,jk) * &275 ( zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) &276 + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &277 +zdk1u(ji,jj+1) + zdku (ji,jj) ) )278 END DO279 END DO280 281 ! | t |282 ! I.3 Horizontal fluxes on V | |283 ! ------------------------=== f---v---f284 ! | |285 ! i-flux at f-point | t |286 DO jj = 1, jpjm1287 DO ji = 1, jpim1288 zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)289 290 zmkf = 1./MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &291 + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )292 293 zcof1 = -e2f(ji,jj) * zmkf &294 * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )295 296 zivf(ji,jj) = fmask(ji,jj,jk) * &297 ( zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) &298 + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) &299 +zdk1u(ji,jj) + zdku (ji+1,jj) ) )300 END DO301 END DO302 303 ! j-flux at t-point304 DO jj = 2, jpj305 DO ji = 1, jpim1306 zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)307 308 zmkt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &309 + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )310 311 zcof2 = -e1t(ji,jj) * zmkt &312 * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )313 314 zjvt(ji,jj) = tmask(ji,jj,jk) * &315 ( zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) &316 + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) &317 +zdk1u(ji,jj-1) + zdku (ji,jj) ) )318 END DO319 END DO320 321 322 ! I.4 Second derivative (divergence) (not divided by the volume)323 ! ---------------------324 325 DO jj = 2, jpjm1326 DO ji = 2, jpim1327 plu(ji,jj,jk) = ziut (ji+1,jj) - ziut (ji,jj ) &328 + zjuf (ji ,jj) - zjuf (ji,jj-1)329 plv(ji,jj,jk) = zivf (ji,jj ) - zivf (ji-1,jj) &330 + zjvt (ji,jj+1) - zjvt (ji,jj )331 END DO332 END DO333 334 ! ! ===============335 END DO ! End of slab336 ! ! ===============337 338 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,339 340 ! ! ************ ! ! ===============341 DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab342 ! ! ************ ! ! ===============343 344 ! II.1 horizontal (pu,pv) gradients345 ! ---------------------------------346 347 DO jk = 1, jpk348 DO ji = 2, jpi349 ! i-gradient of u at jj350 zdiu (ji,jk) = tmask(ji,jj ,jk) * ( pu(ji,jj ,jk) - pu(ji-1,jj ,jk) )351 ! j-gradient of u and v at jj352 zdju (ji,jk) = fmask(ji,jj ,jk) * ( pu(ji,jj+1,jk) - pu(ji ,jj ,jk) )353 zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pv(ji,jj ,jk) - pv(ji ,jj-1,jk) )354 ! j-gradient of u and v at jj+1355 zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj ,jk) - pu(ji ,jj-1,jk) )356 zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji ,jj ,jk) )357 END DO358 END DO359 DO jk = 1, jpk360 DO ji = 1, jpim1361 ! i-gradient of v at jj362 zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pv(ji+1,jj,jk) - pv(ji ,jj ,jk) )363 END DO364 END DO365 366 367 ! II.2 Vertical fluxes368 ! --------------------369 370 ! Surface and bottom vertical fluxes set to zero371 372 zfuw(:, 1 ) = 0.e0373 zfvw(:, 1 ) = 0.e0374 zfuw(:,jpk) = 0.e0375 zfvw(:,jpk) = 0.e0376 377 ! interior (2=<jk=<jpk-1) on pu field378 379 DO jk = 2, jpkm1380 DO ji = 2, jpim1381 ! i- and j-slopes at uw-point382 zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )383 zuwslpj = 0.5 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )384 ! coef. for the vertical dirative385 zcoef0 = e1u(ji,jj) * e2u(ji,jj) / fse3u(ji,jj,jk) &386 * ( zuwslpi * zuwslpi + zuwslpj * zuwslpj )387 ! weights for the i-k, j-k averaging at t- and f-points, resp.388 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &389 + tmask(ji,jj,jk )+tmask(ji+1,jj,jk ), 1. )390 zmkf = 1./MAX( fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1) &391 + fmask(ji,jj-1,jk )+fmask(ji,jj,jk ), 1. )392 ! coef. for the horitontal derivative393 zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi394 zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj395 ! vertical flux on u field396 zfuw(ji,jk) = umask(ji,jj,jk) * &397 ( zcoef0 * ( pu (ji,jj,jk-1) - pu (ji,jj,jk) ) &398 + zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &399 +zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &400 + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &401 +zdj1u(ji,jk ) + zdju (ji ,jk ) ) )402 END DO403 END DO404 405 ! interior (2=<jk=<jpk-1) on pv field406 407 DO jk = 2, jpkm1408 DO ji = 2, jpim1409 ! i- and j-slopes at vw-point410 zvwslpi = 0.5 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )411 zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )412 ! coef. for the vertical derivative413 zcoef0 = e1v(ji,jj) * e2v(ji,jj) / fse3v(ji,jj,jk) &414 * ( zvwslpi * zvwslpi + zvwslpj * zvwslpj )415 ! weights for the i-k, j-k averaging at f- and t-points, resp.416 zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &417 + fmask(ji-1,jj,jk )+fmask(ji,jj,jk ), 1. )418 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &419 + tmask(ji,jj,jk )+tmask(ji,jj+1,jk ), 1. )420 ! coef. for the horizontal derivatives421 zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi422 zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj423 ! vertical flux on pv field424 zfvw(ji,jk) = vmask(ji,jj,jk) * &425 ( zcoef0 * ( pv (ji,jj,jk-1) - pv (ji,jj,jk) ) &426 + zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &427 +zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &428 + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &429 +zdjv (ji,jk ) + zdj1v(ji ,jk ) ) )430 END DO431 END DO432 433 434 ! II.3 Divergence of vertical fluxes added to the horizontal divergence435 ! ---------------------------------------------------------------------436 IF( (kahm -nkahm_smag) ==1 ) THEN437 ! multiply the laplacian by the eddy viscosity coefficient438 DO jk = 1, jpkm1439 DO ji = 2, jpim1440 ! eddy coef. divided by the volume element441 zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )442 zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )443 ! vertical divergence444 zuav = zfuw(ji,jk) - zfuw(ji,jk+1)445 zvav = zfvw(ji,jk) - zfvw(ji,jk+1)446 ! harmonic operator applied to (pu,pv) and multiply by ahm447 plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur448 plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr449 END DO450 END DO451 ELSEIF( (kahm +nkahm_smag ) == 2 ) THEN452 ! second call, no multiplication453 DO jk = 1, jpkm1454 DO ji = 2, jpim1455 ! inverse of the volume element456 zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )457 zbvr = 1. / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )458 ! vertical divergence459 zuav = zfuw(ji,jk) - zfuw(ji,jk+1)460 zvav = zfvw(ji,jk) - zfvw(ji,jk+1)461 ! harmonic operator applied to (pu,pv)462 plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur463 plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr464 END DO465 END DO466 ELSE467 IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm468 IF(lwp)WRITE(numout,*) ' We stop'469 STOP 'ldfguv'470 ENDIF471 ! ! ===============472 END DO ! End of slab473 ! ! ===============474 475 CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )476 !477 IF( nn_timing == 1 ) CALL timing_stop('ldfguv')478 !479 END SUBROUTINE ldfguv480 481 #else482 !!----------------------------------------------------------------------483 !! Dummy module : NO rotation of mixing tensor484 !!----------------------------------------------------------------------485 CONTAINS486 SUBROUTINE dyn_ldf_bilapg( kt ) ! Dummy routine487 INTEGER, INTENT(in) :: kt488 WRITE(*,*) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt489 END SUBROUTINE dyn_ldf_bilapg490 #endif491 492 !!======================================================================493 10 END MODULE dynldf_bilapg
Note: See TracChangeset
for help on using the changeset viewer.
