MODULE trazdf_imp !!====================================================================== !! *** MODULE trazdf_imp *** !! Ocean tracers: vertical component of the tracer mixing trend !!====================================================================== !! History : OPA ! 1990-10 (B. Blanke) Original code !! 7.0 ! 1991-11 (G. Madec) !! ! 1992-06 (M. Imbard) correction on tracer trend loops !! ! 1996-01 (G. Madec) statement function for e3 !! ! 1997-05 (G. Madec) vertical component of isopycnal !! ! 1997-07 (G. Madec) geopotential diffusion in s-coord !! ! 2000-08 (G. Madec) double diffusive mixing !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module !! 2.0 ! 2006-11 (G. Madec) New step reorganisation !! 3.2 ! 2009-03 (G. Madec) heat and salt content trends !! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! tra_zdf_imp : Update the tracer trend with the diagonal vertical !! part of the mixing tensor. !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers variables USE dom_oce ! ocean space and time domain variables USE zdf_oce ! ocean vertical physics variables USE trc_oce ! share passive tracers/ocean variables USE domvvl ! variable volume USE ldftra_oce ! ocean active tracers: lateral physics USE ldftra ! lateral mixing type USE ldfslp ! lateral physics: slope of diffusion USE zdfddm ! ocean vertical physics: double diffusion USE traldf_iso_grif ! active tracers: Griffies operator USE in_out_manager ! I/O manager USE lbclnk ! ocean lateral boundary conditions (or mpp link) IMPLICIT NONE PRIVATE PUBLIC tra_zdf_imp ! routine called by step.F90 !! * Substitutions # include "domzgr_substitute.h90" # include "ldftra_substitute.h90" # include "zdfddm_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_zdf_imp( kt, cdtype, p2dt, ptb, pta, kjpt ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_zdf_imp *** !! !! ** Purpose : Compute the trend due to the vertical tracer diffusion !! including the vertical component of lateral mixing (only for 2nd !! order operator, for fourth order it is already computed and add !! to the general trend in traldf.F) and add it to the general trend !! of the tracer equations. !! !! ** Method : The vertical component of the lateral diffusive trends !! is provided by a 2nd order operator rotated along neutral or geo- !! potential surfaces to which an eddy induced advection can be !! added. It is computed using before fields (forward in time) and !! isopycnal or geopotential slopes computed in routine ldfslp. !! !! Second part: vertical trend associated with the vertical physics !! =========== (including the vertical flux proportional to dk[t] !! associated with the lateral mixing, through the !! update of avt) !! The vertical diffusion of the tracer t is given by: !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) !! It is computed using a backward time scheme (t=ta). !! Surface and bottom boundary conditions: no diffusive flux on !! both tracers (bottom, applied through the masked field avt). !! Add this trend to the general trend ta,sa : !! ta = ta + dz( avt dz(t) ) !! if lk_zdfddm=T, use avs for salinity or for passive tracers !! (sa = sa + dz( avs dz(t) ) !! !! Third part: recover avt resulting from the vertical physics !! ========== alone, for further diagnostics (for example to !! compute the turbocline depth in zdfmxl.F90). !! if lk_zdfddm=T, use avt = zavt !! (avs = zavs if lk_zdfddm=T ) !! !! ** Action : - Update (ta) with before vertical diffusion trend !!--------------------------------------------------------------------- USE oce , ONLY : zwd => ua ! ua used as workspace USE oce , ONLY : zws => va ! va - - USE wrk_nemo, ONLY: wrk_use, wrk_release USE wrk_nemo, ONLY: zwi => wrk_3d_1, zwt => wrk_3d_2 ! workspace arrays !! INTEGER , INTENT(in ) :: kt ! ocean time-step index CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) INTEGER , INTENT(in ) :: kjpt ! number of tracers REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend !! INTEGER :: ji, jj, jk, jn ! dummy loop indices REAL(wp) :: zavi, zrhs, znvvl ! local scalars REAL(wp) :: ze3tb, ze3tn, ze3ta ! variable vertical scale factors !!--------------------------------------------------------------------- IF(.NOT. wrk_use(3, 1,2))THEN CALL ctl_stop('tra_zdf_imp : requested workspace arrays unavailable.') RETURN END IF IF( kt == nit000 ) THEN IF(lwp)WRITE(numout,*) IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' zavi = 0._wp ! avoid warning at compilation phase when lk_ldfslp=F ENDIF ! ! I. Local initialization ! ----------------------- zwd(1,:, : ) = 0._wp ; zwd(jpi,:,:) = 0._wp zws(1,:, : ) = 0._wp ; zws(jpi,:,:) = 0._wp zwi(1,:, : ) = 0._wp ; zwi(jpi,:,:) = 0._wp zwt(1,:, : ) = 0._wp ; zwt(jpi,:,:) = 0._wp zwt(:,:,jpk) = 0._wp ; zwt( : ,:,1) = 0._wp ! I.1 Variable volume : to take into account vertical variable vertical scale factors ! ------------------- IF( lk_vvl ) THEN ; znvvl = 1._wp ELSE ; znvvl = 0._wp ENDIF ! II. Vertical trend associated with the vertical physics ! ======================================================= ! (including the vertical flux proportional to dk[t] associated ! with the lateral mixing, through the avt update) ! dk[ avt dk[ (t,s) ] ] diffusive trends ! ! II.0 Matrix construction ! ------------------------ DO jn = 1, kjpt ! ! Matrix construction ! ------------------------ IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN #if defined key_ldfslp IF( ln_traldf_grif ) THEN DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing zavi = ah_wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature) END DO END DO END DO ! update and save of avt (and avs if double diffusive mixing) ELSE IF( l_traldf_rot ) THEN DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature) END DO END DO END DO ELSE ! no rotation but key_ldfslp defined zwt(:,:,:) = avt(:,:,:) ENDIF #else ! No isopycnal diffusion zwt(:,:,:) = avt(:,:,:) #endif ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked) DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point ze3tn = znvvl + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) END DO END DO END DO ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 zwt(ji,jj,1) = zwd(ji,jj,1) END DO END DO DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) END DO END DO END DO ! ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN #if defined key_ldfslp IF( ln_traldf_grif ) THEN DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zavi = ah_wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature) END DO END DO END DO ELSE IF( l_traldf_rot ) THEN DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on salinity) END DO END DO END DO ELSE ! no rotation but key_ldfslp defined zwt(:,:,:) = fsavs(:,:,:) ENDIF #else ! No isopycnal diffusion zwt(:,:,:) = fsavs(:,:,:) #endif ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked) DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point ze3tn = znvvl + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) END DO END DO END DO ! Surface boudary conditions DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1) ! after scale factor at T-point zwi(ji,jj,1) = 0._wp zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) END DO END DO ! ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 zwt(ji,jj,1) = zwd(ji,jj,1) END DO END DO DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) END DO END DO END DO ! END IF ! II.1. Vertical diffusion on tracer ! --------------------------- ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1) ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1) pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn) END DO END DO DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk) ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,jk) zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn) ! zrhs=right hand side pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) END DO END DO END DO ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) END DO END DO DO jk = jpk-2, 1, -1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) & & / zwt(ji,jj,jk) * tmask(ji,jj,jk) END DO END DO END DO ! END DO ! IF(.NOT. wrk_release(3, 1,2))THEN CALL ctl_stop('tra_zdf_imp : failed to release workspace arrays.') END IF ! END SUBROUTINE tra_zdf_imp !!============================================================================== END MODULE trazdf_imp