MODULE trazdf_iso_vopt !!============================================================================== !! *** MODULE trazdf_iso_vopt *** !! Ocean active tracers: vertical component of the tracer mixing trend !!============================================================================== #if defined key_ldfslp || defined key_esopa !!---------------------------------------------------------------------- !! 'key_ldfslp' rotation of the lateral mixing tensor !!---------------------------------------------------------------------- !! tra_zdf_iso_vopt : Update the tracer trend with the vertical part of !! the isopycnal or geopotential s-coord. operator and !! the vertical diffusion. vector optimization, use !! k-j-i loops. !! tra_zdf_iso : !! tra_zdf_zdf : !!---------------------------------------------------------------------- !! * Modules used 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 trdtra_oce ! tracers trends diagnostics variables USE ldfslp ! iso-neutral slopes USE zdfddm ! ocean vertical physics: double diffusion USE in_out_manager ! I/O manager USE lbclnk ! ocean lateral boundary conditions (or mpp link) IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC tra_zdf_iso_vopt ! routine called by step.F90 !! * Module variables REAL(wp), DIMENSION(jpk) :: & r2dt ! vertical profile of 2 x time-step !! * Substitutions # include "domzgr_substitute.h90" # include "ldftra_substitute.h90" # include "ldfeiv_substitute.h90" # include "zdfddm_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_zdf_iso_vopt( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_zdf_iso_vopt *** !! !! ** Purpose : !! ** Method : !! ** Action : !! !! History : !! 8.5 ! 02-08 (G. Madec) F90: Free form and module !!--------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * Local variables REAL(wp) :: zta, zsa ! temporary scalars !!--------------------------------------------------------------------- !! OPA 8.5, LODYC-IPSL (2002) !!--------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp)WRITE(numout,*) IF(lwp)WRITE(numout,*) 'tra_zdf_iso_vopt : vertical mixing computation' IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~~~~~ is iso-neutral diffusion : implicit vertical time stepping' #if defined key_diaeiv w_eiv(:,:,:) = 0.e0 #endif ENDIF ! I. vertical extra-diagonal part of the rotated tensor ! ----------------------------------------------------- CALL tra_zdf_iso IF( l_ctl .AND. lwp ) THEN ! print mean trends (used for debugging) zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) ) zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) ) WRITE(numout,*) ' zdf 1- Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl t_ctl = zta ; s_ctl = zsa ENDIF ! II. vertical diffusion (including the vertical diagonal part of the rotated tensor) ! ---------------------- CALL tra_zdf_zdf( kt ) IF( l_ctl .AND. lwp ) THEN ! print mean trends (used for debugging) zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) ) zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) ) WRITE(numout,*) ' zdf 1- Ta: ', zta, ' Sa: ', zsa ENDIF END SUBROUTINE tra_zdf_iso_vopt SUBROUTINE tra_zdf_zdf( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_zdf_zdf *** !! !! ** 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 neural 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 diamld). !! avt = zavt !! (avs = zavs if lk_zdfddm=T ) !! !! 'key_trdtra' defined: trend saved for futher diagnostics. !! !! macro-tasked on vertical slab (jj-loop) !! !! ** Action : - Update (ta,sa) with before vertical diffusion trend !! - Save the trend in (ttrd,strd) ('key_diatrends') !! !! History : !! 6.0 ! 90-10 (B. Blanke) Original code !! 7.0 ! 91-11 (G. Madec) !! ! 92-06 (M. Imbard) correction on tracer trend loops !! ! 96-01 (G. Madec) statement function for e3 !! ! 97-05 (G. Madec) vertical component of isopycnal !! ! 97-07 (G. Madec) geopotential diffusion in s-coord !! ! 00-08 (G. Madec) double diffusive mixing !! 8.5 ! 02-08 (G. Madec) F90: Free form and module !!--------------------------------------------------------------------- !! * Modules used USE oce , ONLY : zwd => ua, & ! ua, va used as zws => va ! workspace !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zavi, zrhs ! temporary scalars REAL(wp), DIMENSION(jpi,jpj,jpk) :: & zwi, zwt, zavsi ! temporary workspace arrays # if defined key_trdtra || defined key_trdmld REAL(wp) :: zta, zsa !temporary scalars # endif !!--------------------------------------------------------------------- !! OPA 8.5, LODYC-IPSL (2002) !!--------------------------------------------------------------------- ! I. Local constant initialization ! -------------------------------- ! ... time step = 2 rdttra ex IF( neuler == 0 .AND. kt == nit000 ) THEN r2dt(:) = rdttra(:) ! restarting with Euler time stepping ELSEIF( kt <= nit000 + 1) THEN r2dt(:) = 2. * rdttra(:) ! leapfrog ENDIF 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 ; zwt(:,:,jpk) = 0.e0 zavsi( 1 ,:,:)=0.e0 ; zavsi(jpi,:,:)=0.e0 ; zavsi(:,:,jpk) = 0.e0 zwt(:,:,1) = 0.e0 zavsi(:,:,1) = 0.e0 ! 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 ! ------------------------ ! update and save of avt (and avs if double diffusive mixing) 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 ! 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. zwi(ji,jj,jk) = - r2dt(jk) * zwt(ji,jj,jk ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) zws(ji,jj,jk) = - r2dt(jk) * zwt(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) zwd(ji,jj,jk) = 1. - 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. zwi(ji,jj,1) = 0.e0 zwd(ji,jj,1) = 1. - 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 ta(ji,jj,1) = tb(ji,jj,1) + r2dt(1) * ta(ji,jj,1) END DO END DO DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 zrhs = tb(ji,jj,jk) + r2dt(jk) * 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 #if defined key_trdtra || defined key_trdmld ! Compute and save the vertical diffusive temperature trends IF( l_traldf_iso ) THEN DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zta = ( ta(ji,jj,jk) - tb(ji,jj,jk) ) / r2dt(jk) ttrd(ji,jj,jk,4) = zta - ta(ji,jj,jk) + ttrd(ji,jj,jk,4) END DO END DO END DO ELSE DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zta = ( ta(ji,jj,jk) - tb(ji,jj,jk) ) / r2dt(jk) ttrd(ji,jj,jk,4) = zta - ta(ji,jj,jk) END DO END DO END DO ENDIF #endif ! 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. zwi(ji,jj,jk) = - r2dt(jk) * zavsi(ji,jj,jk ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) zws(ji,jj,jk) = - r2dt(jk) * zavsi(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) zwd(ji,jj,jk) = 1. - 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. zwi(ji,jj,1) = 0.e0 zwd(ji,jj,1) = 1. - 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 sa(ji,jj,1) = sb(ji,jj,1) + r2dt(1) * sa(ji,jj,1) END DO END DO DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 zrhs = sb(ji,jj,jk) + r2dt(jk) * 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 #if defined key_trdtra || defined key_trdmld ! Compute and save the vertical diffusive temperature trends IF( l_traldf_iso ) THEN DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zsa = ( sa(ji,jj,jk) - sb(ji,jj,jk) ) / r2dt(jk) strd(ji,jj,jk,4) = zsa - sa(ji,jj,jk) + strd(ji,jj,jk,4) END DO END DO END DO ELSE DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zsa = ( sa(ji,jj,jk) - sb(ji,jj,jk) ) / r2dt(jk) strd(ji,jj,jk,4) = zsa - sa(ji,jj,jk) END DO END DO END DO ENDIF #endif END SUBROUTINE tra_zdf_zdf SUBROUTINE tra_zdf_iso !!---------------------------------------------------------------------- !! *** ROUTINE tra_zdf_iso *** !! !! ** Purpose : !! Compute the trend due to the vertical tracer diffusion inclu- !! ding the vertical component of lateral mixing (only for second !! 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 neural or geopo- !! tential surfaces to which an eddy induced advection can be added !! It is computed using before fields (forward in time) and isopyc- !! nal or geopotential slopes computed in routine ldfslp. !! !! First part: vertical trends associated with the lateral mixing !! ========== (excluding the vertical flux proportional to dk[t] ) !! vertical fluxes associated with the rotated lateral mixing: !! zftw =-aht { e2t*wslpi di[ mi(mk(tb)) ] !! + e1t*wslpj dj[ mj(mk(tb)) ] } !! save avt coef. resulting from vertical physics alone in zavt: !! zavt = avt !! update and save in zavt the vertical eddy viscosity coefficient: !! avt = avt + wslpi^2+wslj^2 !! add vertical Eddy Induced advective fluxes (lk_traldf_eiv=T): !! zftw = zftw + { di[aht e2u mi(wslpi)] !! +dj[aht e1v mj(wslpj)] } mk(tb) !! take the horizontal divergence of the fluxes: !! difft = 1/(e1t*e2t*e3t) dk[ zftw ] !! Add this trend to the general trend (ta,sa): !! ta = ta + difft !! !! ** Action : !! Update (ta,sa) arrays with the before vertical diffusion trend !! Save in (ttrd,strd) arrays the trends if 'key_diatrends' defined !! !! History : !! 6.0 ! 90-10 (B. Blanke) Original code !! 7.0 ! 91-11 (G. Madec) !! ! 92-06 (M. Imbard) correction on tracer trend loops !! ! 96-01 (G. Madec) statement function for e3 !! ! 97-05 (G. Madec) vertical component of isopycnal !! ! 97-07 (G. Madec) geopotential diffusion in s-coord !! ! 00-08 (G. Madec) double diffusive mixing !! 8.5 ! 02-08 (G. Madec) F90: Free form and module !!--------------------------------------------------------------------- !! * Modules used USE oce , ONLY : zwx => ua, & ! use ua, va as zwy => va ! workspace arrays !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices #if defined key_partial_steps INTEGER :: iku, ikv #endif REAL(wp) :: & ztavg, zsavg, & ! temporary scalars zcoef0, zcoef3, & ! " " zcoef4, zcoeg3, & ! " " zbtr, zmku, zmkv , & ! " " zuwki, zvwki, zuwk, & ! " " zvwk, ztav, zsav REAL(wp), DIMENSION(jpi,jpj,jpk) :: & zwz, zwt, ztfw, zsfw ! temporary workspace arrays !!--------------------------------------------------------------------- !! OPA 8.5, LODYC-IPSL (2002) !!--------------------------------------------------------------------- ! 0. Local constant initialization ! -------------------------------- ztavg = 0.e0 zsavg = 0.e0 zwx( 1 ,:,:)=0.e0 ; zwx(jpi,:,:)=0.e0 zwy( 1 ,:,:)=0.e0 ; zwy(jpi,:,:)=0.e0 zwz( 1 ,:,:)=0.e0 ; zwz(jpi,:,:)=0.e0 zwt( 1 ,:,:)=0.e0 ; zwt(jpi,:,:)=0.e0 ztfw( 1 ,:,:)=0.e0 ; ztfw(jpi,:,:)=0.e0 zsfw( 1 ,:,:)=0.e0 ; zsfw(jpi,:,:)=0.e0 ! I. Vertical trends associated with lateral mixing ! ------------------------------------------------- ! (excluding the vertical flux proportional to dk[t] ) ! I.1 horizontal tracer gradient ! ------------------------------ DO jk = 1, jpkm1 DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. ! i-gradient of T and S at jj zwx (ji,jj,jk) = ( tb(ji+1,jj,jk)-tb(ji,jj,jk) ) * umask(ji,jj,jk) zwy (ji,jj,jk) = ( sb(ji+1,jj,jk)-sb(ji,jj,jk) ) * umask(ji,jj,jk) ! j-gradient of T and S at jj zwz (ji,jj,jk) = ( tb(ji,jj+1,jk)-tb(ji,jj,jk) ) * vmask(ji,jj,jk) zwt (ji,jj,jk) = ( sb(ji,jj+1,jk)-sb(ji,jj,jk) ) * vmask(ji,jj,jk) END DO END DO END DO # if defined key_partial_steps ! partial steps correction at the bottom ocean level DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. ! last ocean level iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1 ikv = MIN( mbathy(ji,jj), mbathy(ji ,jj+1) ) - 1 ! i-gradient of T and S zwx (ji,jj,iku) = gtu(ji,jj) zwy (ji,jj,iku) = gsu(ji,jj) ! j-gradient of T and S zwz (ji,jj,ikv) = gtv(ji,jj) zwt (ji,jj,ikv) = gsv(ji,jj) END DO END DO #endif ! I.2 Vertical fluxes ! ------------------- ! Surface and bottom vertical fluxes set to zero ztfw(:,:, 1 ) = 0.e0 zsfw(:,:, 1 ) = 0.e0 ztfw(:,:,jpk) = 0.e0 zsfw(:,:,jpk) = 0.e0 ! interior (2=