MODULE dynldf_bilapg !!====================================================================== !! *** MODULE dynldf_bilapg *** !! Ocean dynamics: lateral viscosity trend !!====================================================================== !! History : OPA ! 1997-07 (G. Madec) Original code !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module !! 2.0 ! 2004-08 (C. Talandier) New trends organization !!---------------------------------------------------------------------- #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 : !!---------------------------------------------------------------------- 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 ldfslp ! iso-neutral slopes available USE ldftra_oce, ONLY: ln_traldf_iso ! USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control USE wrk_nemo ! Memory Allocation USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC dyn_ldf_bilapg ! called by step.F90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zfvw , zdiu, zdiv ! 2D workspace (ldfguv) REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdju, zdj1u, zdjv, zdj1v ! 2D workspace (ldfguv) !! * Substitutions # include "domzgr_substitute.h90" # include "ldfdyn_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id: dynldf_bilapg.F90 4990 2014-12-15 16:42:49Z timgraham $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION dyn_ldf_bilapg_alloc() !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_bilapg_alloc *** !!---------------------------------------------------------------------- ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) , & & zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc ) ! IF( dyn_ldf_bilapg_alloc /= 0 ) CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays') END FUNCTION dyn_ldf_bilapg_alloc 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) !! !! ** Action : - Update (ua,va) arrays with the before geopotential !! biharmonic mixing trend. !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! ocean time-step index ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp), POINTER, DIMENSION(:,:,:) :: zwk1, zwk2, zwk3, zwk4 !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('dyn_ldf_bilapg') ! CALL wrk_alloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 ) ! 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,*) '~~~~~~~~~~~~~~' ! ! allocate dyn_ldf_bilapg arrays IF( dyn_ldf_bilapg_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays') ENDIF ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN ! DO jk = 1, jpk ! set the slopes of iso-level DO jj = 2, jpjm1 DO ji = 2, jpim1 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 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.5 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.5 END DO END DO END DO ! Lateral boundary conditions on the slopes CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) !!bug IF( kt == nit000 ) then IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)), & & ' wi', sqrt(MAXVAL(wslpi)) , ' wj', sqrt(MAXVAL(wslpj)) endif !!end ENDIF zwk1(:,:,:) = 0.e0 ; zwk3(:,:,:) = 0.e0 zwk2(:,:,:) = 0.e0 ; zwk4(:,:,:) = 0.e0 ! Laplacian of (ub,vb) multiplied by ahm ! -------------------------------------- CALL ldfguv( ub, vb, zwk1, zwk2, 1 ) ! rotated harmonic operator applied to (ub,vb) ! ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) ) CALL lbc_lnk( zwk1, 'U', -1. ) ; CALL lbc_lnk( zwk2, 'V', -1. ) ! Lateral boundary conditions ! Bilaplacian of (ub,vb) ! ---------------------- CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 ) ! rotated harmonic operator applied to (zwk1,zwk2) ! ! (output in (zwk3,zwk4) ) ! Update the momentum trends ! -------------------------- DO jj = 2, jpjm1 ! add the diffusive trend to the general momentum trends DO jk = 1, jpkm1 DO ji = 2, jpim1 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 ! CALL wrk_dealloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 ) ! IF( nn_timing == 1 ) CALL timing_stop('dyn_ldf_bilapg') ! 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) !!---------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu , pv ! 1st call: before horizontal velocity ! ! 2nd call: ahm x these fields 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 1st call ; =2 2nd call ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zabe1 , zabe2 , zcof1 , zcof2 ! local scalar REAL(wp) :: zcoef0, zcoef3, zcoef4 ! - - REAL(wp) :: zbur, zbvr, zmkt, zmkf, zuav, zvav ! - - REAL(wp) :: zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - - ! REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('ldfguv') ! CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) ! ! ! ********** ! ! =============== 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=