MODULE dynldf_bilapg !!====================================================================== !! *** MODULE dynldf_bilapg *** !! Ocean dynamics: lateral viscosity trend !!====================================================================== #if defined key_ldfslp || defined key_esopa !!---------------------------------------------------------------------- !! 'key_ldfslp' Rotation of mixing tensor !!---------------------------------------------------------------------- !! dyn_ldf_bilapg : update the momentum trend with the horizontal part !! of the horizontal s-coord. bilaplacian diffusion !! ldfguv : !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE ldfdyn_oce ! ocean dynamics lateral physics USE zdf_oce ! ocean vertical physics USE in_out_manager ! I/O manager USE trdmod ! ocean dynamics trends USE trdmod_oce ! ocean variables trends USE ldfslp ! iso-neutral slopes available USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC dyn_ldf_bilapg ! called by step.F90 !! * Substitutions # include "domzgr_substitute.h90" # include "ldfdyn_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_ldf_bilapg( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_bilapg *** !! !! ** Purpose : Compute the before trend of the horizontal momentum !! diffusion and add it to the general trend of momentum equation. !! !! ** Method : The lateral momentum diffusive trends is provided by a !! a 4th order operator rotated along geopotential surfaces. It is !! computed using before fields (forward in time) and geopotential !! slopes computed in routine inildf. !! -1- compute the geopotential harmonic operator applied to !! (ub,vb) and multiply it by the eddy diffusivity coefficient !! (done by a call to ldfgpu and ldfgpv routines) The result is in !! (zwk1,zwk2) arrays. Applied the domain lateral boundary conditions !! by call to lbc_lnk. !! -2- applied to (zwk1,zwk2) the geopotential harmonic operator !! by a second call to ldfgpu and ldfgpv routines respectively. The !! result is in (zwk3,zwk4) arrays. !! -3- Add this trend to the general trend (ta,sa): !! (ua,va) = (ua,va) + (zwk3,zwk4) !! 'key_trddyn' defined: the trend is saved for diagnostics. !! !! ** Action : - Update (ua,va) arrays with the before geopotential !! biharmonic mixing trend. !! - save the trend in (zwk3,zwk4) ('key_trddyn') !! !! History : !! 8.0 ! 97-07 (G. Madec) Original code !! 8.5 ! 02-08 (G. Madec) F90: Free form and module !! 9.0 ! 04-08 (C. Talandier) New trends organization !!---------------------------------------------------------------------- !! * Modules used USE oce, ONLY : zwk3 => ta, & ! use ta as 3D workspace zwk4 => sa ! use sa as 3D workspace !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp), DIMENSION(jpi,jpj,jpk) :: & zwk1, zwk2 ! work array used for rotated biharmonic ! ! operator on tracers and/or momentum !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' zwk1(:,:,:) = 0.e0 ; zwk3(:,:,:) = 0.e0 zwk2(:,:,:) = 0.e0 ; zwk4(:,:,:) = 0.e0 ENDIF ! Laplacian of (ub,vb) multiplied by ahm ! -------------------------------------- ! rotated harmonic operator applied to (ub,vb) ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) ) CALL ldfguv ( ub, vb, zwk1, zwk2, 1 ) ! Lateral boundary conditions on (zwk1,zwk2) CALL lbc_lnk( zwk1, 'U', -1. ) CALL lbc_lnk( zwk2, 'V', -1. ) ! Bilaplacian of (ub,vb) ! ---------------------- ! rotated harmonic operator applied to (zwk1,zwk2) (output in (zwk3,zwk4) ) CALL ldfguv ( zwk1, zwk2, zwk3, zwk4, 2 ) ! Update the momentum trends (j-slab : 2, jpj-1) ! -------------------------- ! ! =============== DO jj = 2, jpjm1 ! Vertical slab ! ! =============== DO jk = 1, jpkm1 DO ji = 2, jpim1 ! add the diffusive trend to the general momentum trends ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk) va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== END SUBROUTINE dyn_ldf_bilapg SUBROUTINE ldfguv( pu, pv, plu, plv, kahm ) !!---------------------------------------------------------------------- !! *** ROUTINE ldfguv *** !! !! ** Purpose : Apply a geopotential harmonic operator to (pu,pv) !! (defined at u- and v-points) and multiply it by the eddy !! viscosity coefficient (if kahm=1). !! !! ** Method : The harmonic operator rotated along geopotential !! surfaces is applied to (pu,pv) using the slopes of geopotential !! surfaces computed in inildf routine. The result is provided in !! (plu,plv) arrays. It is computed in 2 stepv: !! !! First step: horizontal part of the operator. It is computed on !! ========== pu as follows (idem on pv) !! horizontal fluxes : !! zftu = e2u*e3u/e1u di[ pu ] - e2u*uslp dk[ mi(mk(pu)) ] !! zftv = e1v*e3v/e2v dj[ pu ] - e1v*vslp dk[ mj(mk(pu)) ] !! take the horizontal divergence of the fluxes (no divided by !! the volume element : !! plu = di-1[ zftu ] + dj-1[ zftv ] !! !! Second step: vertical part of the operator. It is computed on !! =========== pu as follows (idem on pv) !! vertical fluxes : !! zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2) dk-1[ pu ] !! - e2t * wslpi di[ mi(mk(pu)) ] !! - e1t * wslpj dj[ mj(mk(pu)) ] !! take the vertical divergence of the fluxes add it to the hori- !! zontal component, divide the result by the volume element and !! if kahm=1, multiply by the eddy diffusivity coefficient: !! plu = aht / (e1t*e2t*e3t) { plu + dk[ zftw ] } !! else: !! plu = 1 / (e1t*e2t*e3t) { plu + dk[ zftw ] } !! !! ** Action : !! plu, plv : partial harmonic operator applied to !! pu and pv (all the components except !! second order vertical derivative term) !! 'key_trddyn' defined: the trend is saved for diagnostics. !! !! History : !! 8.0 ! 97-07 (G. Madec) Original code !! 8.5 ! 02-08 (G. Madec) F90: Free form and module !!---------------------------------------------------------------------- !! * Arguments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & pu, pv ! momentum fields (before u and v for the 1st call, and ! ! laplacian of these fields multiplied by ahm for the 2nd REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: & plu, plv ! partial harmonic operator applied to ! ! pu and pv (all the components except ! ! second order vertical derivative term) INTEGER, INTENT( in ) :: & kahm ! =1 the laplacian is multiplied by the eddy diffusivity coef. ! ! =2 no multiplication !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zabe1, zabe2, zcof1, zcof2, & ! temporary scalars zcoef0, zcoef3, zcoef4 REAL(wp) :: & zbur, zbvr, zmkt, zmkf, zuav, zvav, & zuwslpi, zuwslpj, zvwslpi, zvwslpj REAL(wp), DIMENSION(jpi,jpj) :: & ziut, zjuf , zjvt, zivf, & ! workspace zdku, zdk1u, zdkv, zdk1v REAL(wp), DIMENSION(jpi,jpk) :: & zfuw, zfvw, zdiu, zdiv, & ! workspace zdju, zdj1u, zdjv, zdj1v !!---------------------------------------------------------------------- ! ! ********** ! ! =============== DO jk = 1, jpkm1 ! First step ! ! Horizontal slab ! ! ********** ! ! =============== ! I.1 Vertical gradient of pu and pv at level jk and jk+1 ! ------------------------------------------------------- ! surface boundary condition: zdku(jk=1)=zdku(jk=2) ! zdkv(jk=1)=zdkv(jk=2) zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1) zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1) IF( jk == 1 ) THEN zdku(:,:) = zdk1u(:,:) zdkv(:,:) = zdk1v(:,:) ELSE zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk) zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk) ENDIF ! -----f----- ! I.2 Horizontal fluxes on U | ! ------------------------=== t u t ! | ! i-flux at t-point -----f----- DO jj = 1, jpjm1 DO ji = 2, jpi zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) zmkt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. ) zcof1 = -e2t(ji,jj) * zmkt & * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ziut(ji,jj) = tmask(ji,jj,jk) * & ( zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) END DO END DO ! j-flux at f-point DO jj = 1, jpjm1 DO ji = 1, jpim1 zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) zmkf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. ) zcof2 = -e1f(ji,jj) * zmkf & * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) zjuf(ji,jj) = fmask(ji,jj,jk) * & ( zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) END DO END DO ! | t | ! I.3 Horizontal fluxes on V | | ! ------------------------=== f---v---f ! | | ! i-flux at f-point | t | DO jj = 1, jpjm1 DO ji = 1, jpim1 zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) zmkf = 1./MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. ) zcof1 = -e2f(ji,jj) * zmkf & * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) zivf(ji,jj) = fmask(ji,jj,jk) * & ( zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) END DO END DO ! j-flux at t-point DO jj = 2, jpj DO ji = 1, jpim1 zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) zmkt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) zcof2 = -e1t(ji,jj) * zmkt & * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) zjvt(ji,jj) = tmask(ji,jj,jk) * & ( zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) END DO END DO ! I.4 Second derivative (divergence) (not divided by the volume) ! --------------------- DO jj = 2, jpjm1 DO ji = 2, jpim1 plu(ji,jj,jk) = ziut (ji+1,jj) - ziut (ji,jj ) & + zjuf (ji ,jj) - zjuf (ji,jj-1) plv(ji,jj,jk) = zivf (ji,jj ) - zivf (ji-1,jj) & + zjvt (ji,jj+1) - zjvt (ji,jj ) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, ! ! ************ ! ! =============== DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab ! ! ************ ! ! =============== ! II.1 horizontal (pu,pv) gradients ! --------------------------------- DO jk = 1, jpk DO ji = 2, jpi ! i-gradient of u at jj zdiu (ji,jk) = tmask(ji,jj ,jk) * ( pu(ji,jj ,jk) - pu(ji-1,jj ,jk) ) ! j-gradient of u and v at jj zdju (ji,jk) = fmask(ji,jj ,jk) * ( pu(ji,jj+1,jk) - pu(ji ,jj ,jk) ) zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pv(ji,jj ,jk) - pv(ji ,jj-1,jk) ) ! j-gradient of u and v at jj+1 zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj ,jk) - pu(ji ,jj-1,jk) ) zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji ,jj ,jk) ) END DO END DO DO jk = 1, jpk DO ji = 1, jpim1 ! i-gradient of v at jj zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pv(ji+1,jj,jk) - pv(ji ,jj ,jk) ) END DO END DO ! II.2 Vertical fluxes ! -------------------- ! Surface and bottom vertical fluxes set to zero zfuw(:, 1 ) = 0.e0 zfvw(:, 1 ) = 0.e0 zfuw(:,jpk) = 0.e0 zfvw(:,jpk) = 0.e0 ! interior (2=