MODULE traldf_bilapg !!============================================================================== !! *** MODULE traldf_bilapg *** !! Ocean active tracers: horizontal component of the lateral tracer mixing trend !!============================================================================== #if defined key_ldfslp || defined key_esopa !!---------------------------------------------------------------------- !! 'key_ldfslp' rotation of the lateral mixing tensor !!---------------------------------------------------------------------- !! tra_ldf_bilapg : update the tracer trend with the horizontal diffusion !! using an horizontal biharmonic operator in s-coordinate !! ldfght : ??? !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers variables USE dom_oce ! ocean space and time domain variables USE ldftra_oce ! ocean active tracers: lateral physics USE trdmod ! ocean active tracers trends USE trdmod_oce ! ocean variables trends USE in_out_manager ! I/O manager USE ldfslp ! iso-neutral slopes available USE lbclnk ! ocean lateral boundary condition (or mpp link) USE diaptr ! poleward transport diagnostics USE prtctl ! Print control IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC tra_ldf_bilapg ! routine called by step.F90 !! * Substitutions # include "domzgr_substitute.h90" # include "ldftra_substitute.h90" # include "ldfeiv_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Header$ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_ldf_bilapg( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_ldf_bilapg *** !! !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive !! trend and add it to the general trend of tracer equation. !! !! ** Method : The lateral diffusive trends is provided by 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 !! (tb,sb) and multiply it by the eddy diffusivity coefficient !! (done by a call to ldfght routine, result in (wk1,wk2) arrays). !! Applied the domain lateral boundary conditions by call to lbc_lnk !! -2- compute the geopotential harmonic operator applied to !! (wk1,wk2) by a second call to ldfght routine (result in (wk3,wk4) !! arrays). !! -3- Add this trend to the general trend (ta,sa): !! (ta,sa) = (ta,sa) + (wk3,wk4) !! !! ** Action : - Update (ta,sa) arrays with the before geopotential !! biharmonic mixing trend. !! !! 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 , wk1 => ua, & ! use ua as workspace & wk2 => va ! use va as 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) :: & wk3, wk4 ! work array used for rotated biharmonic ! ! operator on tracers and/or momentum !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' ENDIF ! 1. Laplacian of (tb,sb) * aht ! ----------------------------- ! rotated harmonic operator applied to (tb,sb) ! and multiply by aht (output in (wk1,wk2) ) CALL ldfght ( kt, tb, sb, wk1, wk2, 1 ) ! Lateral boundary conditions on (wk1,wk2) (unchanged sign) CALL lbc_lnk( wk1, 'T', 1. ) ; CALL lbc_lnk( wk2, 'T', 1. ) ! 2. Bilaplacian of (tb,sb) ! ------------------------- ! rotated harmonic operator applied to (wk1,wk2) ! (output in (wk3,wk4) ) CALL ldfght ( kt, wk1, wk2, wk3, wk4, 2 ) ! 3. Update the tracer trends (j-slab : 2, jpj-1) ! --------------------------- ! ! =============== DO jj = 2, jpjm1 ! Vertical slab ! ! =============== DO jk = 1, jpkm1 DO ji = 2, jpim1 ! add it to the general tracer trends ta(ji,jj,jk) = ta(ji,jj,jk) + wk3(ji,jj,jk) sa(ji,jj,jk) = sa(ji,jj,jk) + wk4(ji,jj,jk) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== END SUBROUTINE tra_ldf_bilapg SUBROUTINE ldfght ( kt, pt, ps, plt, pls, kaht ) !!---------------------------------------------------------------------- !! *** ROUTINE ldfght *** !! !! ** Purpose : Apply a geopotential harmonic operator to (pt,ps) and !! multiply it by the eddy diffusivity coefficient (if kaht=1). !! Routine only used in s-coordinates (l_sco=T) with bilaplacian !! operator (ln_traldf_bilap=T) acting along geopotential surfaces !! (ln_traldf_hor). !! !! ** Method : The harmonic operator rotated along geopotential !! surfaces is applied to (pt,ps) using the slopes of geopotential !! surfaces computed in inildf routine. The result is provided in !! (plt,pls) arrays. It is computed in 2 steps: !! !! First step: horizontal part of the operator. It is computed on !! ========== pt as follows (idem on ps) !! horizontal fluxes : !! zftu = e2u*e3u/e1u di[ pt ] - e2u*uslp dk[ mi(mk(pt)) ] !! zftv = e1v*e3v/e2v dj[ pt ] - e1v*vslp dk[ mj(mk(pt)) ] !! take the horizontal divergence of the fluxes (no divided by !! the volume element : !! plt = di-1[ zftu ] + dj-1[ zftv ] !! !! Second step: vertical part of the operator. It is computed on !! =========== pt as follows (idem on ps) !! vertical fluxes : !! zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2) dk-1[ pt ] !! - e2t * wslpi di[ mi(mk(pt)) ] !! - e1t * wslpj dj[ mj(mk(pt)) ] !! take the vertical divergence of the fluxes add it to the hori- !! zontal component, divide the result by the volume element and !! if kaht=1, multiply by the eddy diffusivity coefficient: !! plt = aht / (e1t*e2t*e3t) { plt + dk[ zftw ] } !! else: !! plt = 1 / (e1t*e2t*e3t) { plt + dk[ zftw ] } !! !! * Action : !! 'key_trdtra' 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 INTEGER, INTENT( in ) :: kt ! ocean time-step index REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & pt, ps ! tracer fields (before t and s for 1st call ! ! and laplacian of these fields for 2nd call. REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: & plt, pls ! partial harmonic operator applied to ! ! pt & ps components except ! ! second order vertical derivative term) INTEGER, INTENT( in ) :: & kaht ! =1 multiply the laplacian by the eddy diffusivity coeff. ! ! =2 no multiplication !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zabe1, zabe2, zmku, zmkv, & ! temporary scalars zbtr, ztah, zsah, ztav, zsav, & zcof0, zcof1, zcof2, & zcof3, zcof4 REAL(wp), DIMENSION(jpi,jpj) :: & zftu, zfsu, & ! workspace zdkt, zdk1t, & zdks, zdk1s REAL(wp), DIMENSION(jpi,jpj,jpk) :: & zftv, zfsv ! workspace (only v components for ptr) REAL(wp), DIMENSION(jpi,jpk) :: & zftw, zfsw, & ! workspace zdit, zdjt, zdj1t, & zdis, zdjs, zdj1s !!---------------------------------------------------------------------- ! ! ********** ! ! =============== DO jk = 1, jpkm1 ! First step ! ! Horizontal slab ! ! ********** ! ! =============== ! I.1 Vertical gradient of pt and ps at level jk and jk+1 ! ------------------------------------------------------- ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) zdk1t(:,:) = ( pt(:,:,jk) - pt(:,:,jk+1) ) * tmask(:,:,jk+1) zdk1s(:,:) = ( ps(:,:,jk) - ps(:,:,jk+1) ) * tmask(:,:,jk+1) IF( jk == 1 ) THEN zdkt(:,:) = zdk1t(:,:) zdks(:,:) = zdk1s(:,:) ELSE zdkt(:,:) = ( pt(:,:,jk-1) - pt(:,:,jk) ) * tmask(:,:,jk) zdks(:,:) = ( ps(:,:,jk-1) - ps(:,:,jk) ) * tmask(:,:,jk) ENDIF ! I.2 Horizontal fluxes ! --------------------- DO jj = 1, jpjm1 DO ji = 1, jpim1 zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) zmku=1./MAX( tmask(ji+1,jj,jk )+tmask(ji,jj,jk+1) & +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk ),1. ) zmkv=1./MAX( tmask(ji,jj+1,jk )+tmask(ji,jj,jk+1) & +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk ),1. ) zcof1= -e2u(ji,jj) * uslp(ji,jj,jk) * zmku zcof2= -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv zftu(ji,jj)= umask(ji,jj,jk) * & ( zabe1 *( pt(ji+1,jj,jk) - pt(ji,jj,jk) ) & + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj) & +zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) zftv(ji,jj,jk)= vmask(ji,jj,jk) * & ( zabe2 *( pt(ji,jj+1,jk) - pt(ji,jj,jk) ) & + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj) & +zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) zfsu(ji,jj)= umask(ji,jj,jk) * & ( zabe1 *( ps(ji+1,jj,jk) - ps(ji,jj,jk) ) & + zcof1 *( zdks (ji+1,jj) + zdk1s(ji,jj) & +zdk1s(ji+1,jj) + zdks (ji,jj) ) ) zfsv(ji,jj,jk)= vmask(ji,jj,jk) * & ( zabe2 *( ps(ji,jj+1,jk) - ps(ji,jj,jk) ) & + zcof2 *( zdks (ji,jj+1) + zdk1s(ji,jj) & +zdk1s(ji,jj+1) + zdks (ji,jj) ) ) END DO END DO ! I.3 Second derivative (divergence) (not divided by the volume) ! --------------------- DO jj = 2 , jpjm1 DO ji = 2 , jpim1 ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) zsah = zfsu(ji,jj) - zfsu(ji-1,jj) + zfsv(ji,jj,jk) - zfsv(ji,jj-1,jk) plt(ji,jj,jk) = ztah pls(ji,jj,jk) = zsah END DO END DO ! ! =============== END DO ! End of slab ! ! =============== !!but this should be done somewhere after ! "zonal" mean diffusive heat and salt transport IF( ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN pht_ldf(:) = ptr_vj( zftv(:,:,:) ) pst_ldf(:) = ptr_vj( zfsv(:,:,:) ) ENDIF ! ! ************ ! ! =============== DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab ! ! ************ ! ! =============== ! II.1 horizontal tracer gradient ! ------------------------------- DO jk = 1, jpk DO ji = 1, jpim1 zdit (ji,jk) = ( pt(ji+1,jj ,jk) - pt(ji,jj ,jk) ) * umask(ji,jj ,jk) zdis (ji,jk) = ( ps(ji+1,jj ,jk) - ps(ji,jj ,jk) ) * umask(ji,jj ,jk) zdjt (ji,jk) = ( pt(ji ,jj+1,jk) - pt(ji,jj ,jk) ) * vmask(ji,jj ,jk) zdjs (ji,jk) = ( ps(ji ,jj+1,jk) - ps(ji,jj ,jk) ) * vmask(ji,jj ,jk) zdj1t(ji,jk) = ( pt(ji ,jj ,jk) - pt(ji,jj-1,jk) ) * vmask(ji,jj-1,jk) zdj1s(ji,jk) = ( ps(ji ,jj ,jk) - ps(ji,jj-1,jk) ) * vmask(ji,jj-1,jk) END DO END DO ! II.2 Vertical fluxes ! -------------------- ! Surface and bottom vertical fluxes set to zero zftw(:, 1 ) = 0.e0 zfsw(:, 1 ) = 0.e0 zftw(:,jpk) = 0.e0 zfsw(:,jpk) = 0.e0 ! interior (2=