MODULE trazdf_imp !!====================================================================== !! *** MODULE trazdf_imp *** !! Ocean active 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 !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! 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 ldftra_oce ! ocean active tracers: lateral physics USE ldfslp ! lateral physics: slope of diffusion USE trdmod ! ocean active tracers trends USE trdmod_oce ! ocean variables trends USE zdfddm ! ocean vertical physics: double diffusion USE in_out_manager ! I/O manager USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control USE domvvl ! variable volume USE ldftra ! lateral mixing type 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.2 , LOCEAN-IPSL (2009) !! $Id$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_zdf_imp( kt, p2dt ) !!---------------------------------------------------------------------- !! *** 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 tracers (t & s) 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) ) !! (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T ) !! !! Third part: recover avt resulting from the vertical physics !! ========== alone, for further diagnostics (for example to !! compute the turbocline depth in zdfmxl.F90). !! avt = zavt !! (avs = zavs if lk_zdfddm=T ) !! !! ** Action : - Update (ta,sa) with before vertical diffusion trend !! !!--------------------------------------------------------------------- USE oce , ONLY : zwd => ua ! ua used as workspace USE oce , ONLY : zws => va ! va - - !! INTEGER , INTENT(in) :: kt ! ocean time-step index REAL(wp), DIMENSION(jpk), INTENT(in) :: p2dt ! vertical profile of tracer time-step !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zavi, zrhs, znvvl ! temporary scalars REAL(wp) :: ze3tb, ze3tn, ze3ta ! variable vertical scale factors REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwt, zavsi ! workspace arrays !!--------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp)WRITE(numout,*) IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing' IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' zavi = 0.e0 ! avoid warning at compilation phase when lk_ldfslp=F ENDIF ! I. Local initialization ! ----------------------- zwd (1,:, : ) = 0.e0 ; zwd (jpi,:,:) = 0.e0 zws (1,:, : ) = 0.e0 ; zws (jpi,:,:) = 0.e0 zwi (1,:, : ) = 0.e0 ; zwi (jpi,:,:) = 0.e0 zwt (1,:, : ) = 0.e0 ; zwt (jpi,:,:) = 0.e0 zavsi(1,:, : ) = 0.e0 ; zavsi(jpi,:,:) = 0.e0 zwt (:,:,jpk) = 0.e0 ; zwt ( : ,:,1) = 0.e0 zavsi(:,:,jpk) = 0.e0 ; zavsi( : ,:,1) = 0.e0 ! I.1 Variable volume : to take into account vertical variable vertical scale factors ! ------------------- IF( lk_vvl ) THEN ; znvvl = 1. ELSE ; znvvl = 0.e0 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 ! ------------------------ #if defined key_ldfslp ! update and save of avt (and avs if double diffusive mixing) 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) # if defined key_zdfddm zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! dd mixing: zavsi = total vertical mixing coef. on salinity # endif END DO END DO END DO ELSE ! no rotation but key_ldfslp defined zwt (:,:,:) = avt(:,:,:) # if defined key_zdfddm zavsi(:,:,:) = avs(:,:,:) ! avs /= avt (double diffusion mixing) # endif ENDIF #else ! No isopycnal diffusion zwt(:,:,:) = avt(:,:,:) # if defined key_zdfddm zavsi(:,:,:) = avs(:,:,:) # endif #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 ) & ! after scale factor at T-point & + znvvl * fse3t_a(ji,jj,jk) ze3tn = znvvl & ! now scale factor at T-point & + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) 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 ) & ! after scale factor at T-point & + znvvl * fse3t_a(ji,jj,1) zwi(ji,jj,1) = 0.e0 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) END DO END DO ! II.1. Vertical diffusion on t ! --------------------------- !! Matrix inversion from the first level !!---------------------------------------------------------------------- ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) ! ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) ! ( ... )( ... ) ( ... ) ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) ! ! m is decomposed in the product of an upper and lower triangular matrix ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi ! The second member is in 2d array zwy ! The solution is in 2d array zwx ! The 3d arry zwt is a work space array ! zwy is used and then used as a work space array : its value is modified! ! 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 ! 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) ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1) 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 * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk) ! zrhs=right hand side ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ta(ji,jj,jk-1) END DO END DO END DO ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk ! Save the masked temperature after in ta ! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / 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 ta(ji,jj,jk) = ( ta(ji,jj,jk) - zws(ji,jj,jk) * ta(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk) END DO END DO END DO ! II.2 Vertical diffusion on salinity ! ----------------------------------- #if defined key_zdfddm ! Rebuild the Matrix as avt /= avs ! Diagonal, inferior, superior (including the bottom boundary condition via avs masked) DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ze3ta = ( 1. - znvvl ) & ! after scale factor at T-point & + znvvl * fse3t_a(ji,jj,jk) ze3tn = znvvl & ! now scale factor at T-point & + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) zws(ji,jj,jk) = - p2dt(jk) * zavsi(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) zwi(ji,jj,1) = 0.e0 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) END DO END DO #endif !! Matrix inversion from the first level !!---------------------------------------------------------------------- ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) ! ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) ! ( ... )( ... ) ( ... ) ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) ! ! m is decomposed in the product of an upper and lower triangular ! matrix ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi ! The second member is in 2d array zwy ! The solution is in 2d array zwx ! The 3d arry zwt is a work space array ! zwy is used and then used as a work space array : its value is modified! ! 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 ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ze3tb = ( 1. - znvvl ) & ! before scale factor at T-point & + znvvl * fse3t_b(ji,jj,1) ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,1) ! now scale factor at T-point sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1) END DO END DO DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ze3tb = ( 1. - znvvl ) & ! before scale factor at T-point & + znvvl * fse3t_b(ji,jj,jk) ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,jk) ! now scale factor at T-point zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk) ! zrhs=right hand side sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1) END DO END DO END DO ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk ! Save the masked temperature after in ta ! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / 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 sa(ji,jj,jk) = ( sa(ji,jj,jk) - zws(ji,jj,jk) * sa(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk) END DO END DO END DO ! END SUBROUTINE tra_zdf_imp !!============================================================================== END MODULE trazdf_imp