MODULE trcadv_cen2 !!====================================================================== !! *** MODULE trcadv_cen2 *** !! Ocean passive tracers: horizontal & vertical advective tracer trend !!====================================================================== #if defined key_top !!---------------------------------------------------------------------- !! 'key_top' TOP models !!---------------------------------------------------------------------- !! trc_adv_cen2 : update the tracer trend with the horizontal !! and vertical advection trends using a 2nd order !! centered finite difference scheme !!---------------------------------------------------------------------- USE oce_trc ! ocean dynamics and active tracers variables USE trp_trc ! ocean passive tracers variables USE trcbbl ! advective passive tracers in the BBL USE prtctl_trc USE trdmld_trc USE trdmld_trc_oce ! ocean variables trends IMPLICIT NONE PRIVATE PUBLIC trc_adv_cen2 ! routine called by trcstp.F90 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: upsmsk !: mixed upstream/centered scheme near some straits ! ! and in closed seas (orca 2 and 4 configurations) !! * Substitutions # include "top_substitute.h90" !!---------------------------------------------------------------------- !! TOP 1.0 , LOCEAN-IPSL (2005) !! $Id$ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- CONTAINS !!---------------------------------------------------------------------- !! Default option : 2nd order centered scheme (k-j-i loop) !!---------------------------------------------------------------------- SUBROUTINE trc_adv_cen2( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE trc_adv_cen2 *** !! !! ** Purpose : Compute the now trend due to the advection of tracers !! and add it to the general trend of passive tracer equations. !! !! ** Method : The advection is evaluated by a second order centered !! scheme using now fields (leap-frog scheme). In specific areas !! (vicinity of major river mouths, some straits, or where tn is !! is approaching the freezing point) it is mixed with an upstream !! scheme for stability reasons. !! Part 0 : compute the upstream / centered flag !! (3D array, zind, defined at T-point (00 or <0] !! zcenu = centered flux = wn * mk(trn) !! The surface boundary condition is : !! rigid-lid (default option) : zero advective flux !! free-surf ("key_fresurf_cstvol") : wn(:,:,1) * trn(:,:,1) !! Add this trend now to the general trend of tracer tra : !! tra = tra + ztra !! Trend diagnostic ('key_trdmld_trc'): the trend is saved for !! diagnostics. The trends saved is expressed as : !! save trend = w.gradz(T) = ztra - trn divn. !! !! ** Action : - update tra with the now advective tracer trends !! - save the trends in trtrd ('key_trdmld_trc') !! !! History : !! 8.2 ! 01-08 (M-A Filiberti, and M.Levy) trahad+trazad = traadv !! 8.5 ! 02-06 (G. Madec, C. Ethe) F90: Free form and module !!---------------------------------------------------------------------- !! * Modules used USE oce_trc , zwx => ua, & ! use ua as workspace & zwy => va ! use va as workspace #if defined key_trcbbl_adv REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! temporary arrays & zun, zvn, zwn #else USE oce_trc , zun => un, & ! When no bbl, zun == un & zvn => vn, & ! When no bbl, zvn == vn & zwn => wn ! When no bbl, zwn == wn #endif !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * Local save REAL(wp), DIMENSION(jpi,jpj), SAVE :: & zbtr2 !! * Local declarations INTEGER :: ji, jj, jk, jn ! dummy loop indices REAL(wp) :: & zbtr, ztra, zfui, zfvj, & ! temporary scalars zhw, ze3tr, zcofi, zcofj, & ! " " zupsut, zupsvt, & ! " " zfp_ui, zfp_vj, zfm_ui, zfm_vj, & ! " " zcofk, zupst, zcent, & ! " " zfp_w, zfm_w, & ! " " zcenut, zcenvt ! REAL(wp), DIMENSION(jpi,jpj,jpk) :: & zind ! temporary workspace arrays REAL(wp) :: & ztai, ztaj, & ! temporary scalars zfui1, zfvj1 ! " " REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrtrd #if defined key_lim3 || defined key_lim2 REAL(wp) :: & ztfreez ! freezing point #endif CHARACTER (len=22) :: charout !!---------------------------------------------------------------------- IF( kt == nittrc000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'trc_adv_cen2 : 2nd order centered advection scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ Vector optimization case' IF(lwp) WRITE(numout,*) upsmsk(:,:) = 0.e0 ! not upstream by default IF( cp_cfg == "orca" ) CALL ups_orca_set ! set mixed Upstream/centered scheme near some straits ! ! and in closed seas (orca2 and orca4 only) zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) ENDIF IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) ) #if defined key_trcbbl_adv ! Advective bottom boundary layer ! ------------------------------- zun(:,:,:) = un(:,:,:) - u_trc_bbl(:,:,:) zvn(:,:,:) = vn(:,:,:) - v_trc_bbl(:,:,:) zwn(:,:,:) = wn(:,:,:) + w_trc_bbl(:,:,:) #endif ! Upstream / centered scheme indicator ! ------------------------------------ DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi #if defined key_lim3 || defined key_lim2 ztfreez = ( - 0.0575 + 1.710523e-3 * SQRT( sn(ji,jj,1) ) & & - 2.154996e-4 * sn(ji,jj,1) ) * sn(ji,jj,1) zind(ji,jj,jk) = MAX ( & rnfmsk(ji,jj) * rnfmsk_z(jk), & ! near runoff mouths (& closed sea outflows) upsmsk(ji,jj) & ! some of some straits ! ! below ice covered area (if tn < "freezing"+0.1 ) , MAX( 0., SIGN( 1., ztfreez + 0.1 - tn(ji,jj,jk) ) ) * tmask(ji,jj,jk) & & ) #else zind(ji,jj,jk) = MAX ( & rnfmsk(ji,jj) * rnfmsk_z(jk), & ! near runoff mouths (& closed sea outflows) upsmsk(ji,jj) & ! some of some straits & ) #endif END DO END DO END DO DO jn = 1, jptra ! I. Horizontal advective fluxes ! ------------------------------ ! Second order centered tracer flux at u and v-points ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. ! upstream indicator zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) ) zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) ) ! volume fluxes * 1/2 #if ! defined key_zco zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk) zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk) #else zfui = 0.5 * e2u(ji,jj) * zun(ji,jj,jk) zfvj = 0.5 * e1v(ji,jj) * zvn(ji,jj,jk) #endif ! upstream scheme zfp_ui = zfui + ABS( zfui ) zfp_vj = zfvj + ABS( zfvj ) zfm_ui = zfui - ABS( zfui ) zfm_vj = zfvj - ABS( zfvj ) zupsut = zfp_ui * trb(ji,jj,jk,jn) + zfm_ui * trb(ji+1,jj ,jk,jn) zupsvt = zfp_vj * trb(ji,jj,jk,jn) + zfm_vj * trb(ji ,jj+1,jk,jn) ! centered scheme zcenut = zfui * ( trn(ji,jj,jk,jn) + trn(ji+1,jj ,jk,jn) ) zcenvt = zfvj * ( trn(ji,jj,jk,jn) + trn(ji ,jj+1,jk,jn) ) ! mixed centered / upstream scheme zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt END DO END DO ! 2. Tracer flux divergence at t-point added to the general trend ! ------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #if ! defined key_zco zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk) #else zbtr = zbtr2(ji,jj) #endif ! horizontal advective trends ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) ! add it to the general tracer trends tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra #if defined key_trc_diatrd ! recompute the trends in i- and j-direction as Uh gradh(T) # if defined key_s_coord || defined key_partial_steps zfui = 0.5 * e2u(ji ,jj) * fse3u(ji, jj,jk) * zun(ji, jj,jk) zfui1= 0.5 * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zun(ji-1,jj,jk) zfvj = 0.5 * e1v(ji,jj ) * fse3v(ji,jj ,jk) * zvn(ji,jj ,jk) zfvj1= 0.5 * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * zvn(ji,jj-1,jk) # else zfui = 0.5 * e2u(ji ,jj) * zun(ji, jj,jk) zfui1= 0.5 * e2u(ji-1,jj) * zun(ji-1,jj,jk) zfvj = 0.5 * e1v(ji,jj ) * zvn(ji,jj ,jk) zfvj1= 0.5 * e1v(ji,jj-1) * zvn(ji,jj-1,jk) # endif ztai = - zbtr * ( zfui * ( trn(ji+1,jj ,jk,jn) - trn(ji, jj,jk,jn) ) & & + zfui1 * ( trn(ji, jj, jk,jn) - trn(ji-1,jj,jk,jn) ) ) ztaj = - zbtr * ( zfvj * ( trn(ji ,jj+1,jk,jn) - trn(ji,jj ,jk,jn) ) & & + zfvj1 * ( trn(ji ,jj ,jk,jn) - trn(ji,jj-1,jk,jn) ) ) ! save i- and j- advective trends computed as Uh gradh(T) IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = ztai IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = ztaj #endif END DO END DO ! ! =============== END DO ! End of slab ! ! =============== ! 3. Save the horizontal advective trends for diagnostics ! ------------------------------------------------------- !CDIR BEGIN COLLAPSE TRDTRC_XY : IF( l_trdtrc )THEN ! 3.1) Passive tracer ZONAL advection trends ztrtrd(:,:,:) = 0.e0 DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! recompute the trends in i-direction as Uh gradh(T) # if ! defined key_zco zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk) zfui = 0.5 * e2u(ji ,jj) * fse3u(ji, jj,jk) * zun(ji, jj,jk) zfui1= 0.5 * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zun(ji-1,jj,jk) # else zbtr = zbtr2(ji,jj) zfui = 0.5 * e2u(ji ,jj) * zun(ji, jj,jk) zfui1= 0.5 * e2u(ji-1,jj) * zun(ji-1,jj,jk) # endif ztai = - zbtr * ( zfui * ( trn(ji+1,jj ,jk,jn) - trn(ji, jj,jk,jn) ) & & + zfui1 * ( trn(ji, jj, jk,jn) - trn(ji-1,jj,jk,jn) ) ) ! save i- and j- advective trends computed as Uh gradh(T) ztrtrd(ji,jj,jk) = ztai END DO END DO END DO IF( luttrd(jn) ) CALL trd_mod_trc(ztrtrd, jn, jptrc_trd_xad, kt) ! handle the trend ! 3.2) Passive tracer MERIDIONAL advection trends ztrtrd(:,:,:) = 0.e0 DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! recompute the trends in j-direction as Uh gradh(T) # if ! defined key_zco zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk) zfvj = 0.5 * e1v(ji,jj ) * fse3v(ji,jj ,jk) * zvn(ji,jj ,jk) zfvj1= 0.5 * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * zvn(ji,jj-1,jk) # else zbtr = zbtr2(ji,jj) zfvj = 0.5 * e1v(ji,jj ) * zvn(ji,jj ,jk) zfvj1= 0.5 * e1v(ji,jj-1) * zvn(ji,jj-1,jk) # endif ztaj = - zbtr * ( zfvj * ( trn(ji ,jj+1,jk,jn) - trn(ji,jj ,jk,jn) ) & & + zfvj1 * ( trn(ji ,jj ,jk,jn) - trn(ji,jj-1,jk,jn) ) ) ! save i- and j- advective trends computed as Uh gradh(T) ztrtrd(ji,jj,jk) = ztaj END DO END DO END DO IF( luttrd(jn) ) CALL trd_mod_trc(ztrtrd, jn, jptrc_trd_yad, kt) ! handle the trend ENDIF TRDTRC_XY !CDIR END ! ! =========== END DO ! tracer loop ! ! =========== IF(ln_ctl) THEN ! print mean trends (used for debugging) WRITE(charout, FMT="('centered2 - had')") CALL prt_ctl_trc_info(charout) CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd') ENDIF ! II. Vertical advection ! ---------------------- DO jn = 1, jptra ! Bottom value : flux set to zero zwx(:,:,jpk) = 0.e0 ! Surface value IF ( lk_dynspg_rl ) THEN ! rigid lid : flux set to zero zwx(:,:, 1 ) = 0.e0 ELSE ! free surface-constant volume zwx(:,:, 1 ) = zwn(:,:,1) * trn(:,:,1,jn) ENDIF ! 1. Vertical advective fluxes ! ---------------------------- ! Second order centered tracer flux at w-point DO jk = 2, jpk DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! upstream indicator zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) ! velocity * 1/2 zhw = 0.5 * zwn(ji,jj,jk) ! upstream scheme zfp_w = zhw + ABS( zhw ) zfm_w = zhw - ABS( zhw ) zupst = zfp_w * trb(ji,jj,jk,jn) + zfm_w * trb(ji,jj,jk-1,jn) ! centered scheme zcent = zhw * ( trn(ji,jj,jk,jn) + trn(ji,jj,jk-1,jn) ) ! centered scheme zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent END DO END DO END DO ! 2. Tracer flux divergence at t-point added to the general trend ! ------------------------- DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ze3tr = 1. / fse3t(ji,jj,jk) ! vertical advective trends ztra = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) ! add it to the general tracer trends tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra #if defined key_trc_diatrd ! save the vertical advective trends computed as w gradz(T) IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk) #endif END DO END DO END DO ! 3. Save the vertical advective trends for diagnostic ! ---------------------------------------------------- !CDIR BEGIN COLLAPSE TRDTRC_Z : IF( l_trdtrc )THEN ztrtrd(:,:,:) = 0.e0 ! Compute T/S vertical advection trends DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ze3tr = 1. / fse3t(ji,jj,jk) ! vertical advective trends ztra = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) ! save the vertical advective trends computed as w gradz(T) ztrtrd(ji,jj,jk) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk) END DO END DO END DO IF( luttrd(jn) ) CALL trd_mod_trc(ztrtrd, jn, jptrc_trd_zad, kt) ! handle the trend ENDIF TRDTRC_Z !CDIR END ! ! =========== END DO ! tracer loop ! ! =========== IF( l_trdtrc ) DEALLOCATE( ztrtrd ) IF(ln_ctl) THEN ! print mean trends (used for debugging) WRITE(charout, FMT="('centered - zad')") CALL prt_ctl_trc_info(charout) CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd') ENDIF END SUBROUTINE trc_adv_cen2 SUBROUTINE ups_orca_set !!---------------------------------------------------------------------- !! *** ROUTINE ups_orca_set *** !! !! ** Purpose : add a portion of upstream scheme in area where the !! centered scheme generates too strong overshoot !! !! ** Method : orca (R4 and R2) confiiguration setting. Set upsmsk !! array to nozero value in some straith. !! !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca !!---------------------------------------------------------------------- INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers !!---------------------------------------------------------------------- ! mixed upstream/centered scheme near river mouths ! ------------------------------------------------ SELECT CASE ( jp_cfg ) ! ! ======================= CASE ( 4 ) ! ORCA_R4 configuration ! ! ======================= ! ! Gibraltar Strait ii0 = 70 ; ii1 = 71 ij0 = 52 ; ij1 = 53 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 ! ! ! ======================= CASE ( 2 ) ! ORCA_R2 configuration ! ! ======================= ! ! Gibraltar Strait ij0 = 102 ; ij1 = 102 ii0 = 138 ; ii1 = 138 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20 ii0 = 139 ; ii1 = 139 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40 ii0 = 140 ; ii1 = 140 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 ij0 = 101 ; ij1 = 102 ii0 = 141 ; ii1 = 141 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 ! ! Bab el Mandeb Strait ij0 = 87 ; ij1 = 88 ii0 = 164 ; ii1 = 164 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10 ij0 = 88 ; ij1 = 88 ii0 = 163 ; ii1 = 163 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 ii0 = 162 ; ii1 = 162 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40 ii0 = 160 ; ii1 = 161 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 ij0 = 89 ; ij1 = 89 ii0 = 158 ; ii1 = 160 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 ij0 = 90 ; ij1 = 90 ii0 = 160 ; ii1 = 160 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 ! ! Sound Strait ij0 = 116 ; ij1 = 116 ii0 = 145 ; ii1 = 147 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 ! END SELECT END SUBROUTINE ups_orca_set #else !!---------------------------------------------------------------------- !! Default option Empty module !!---------------------------------------------------------------------- CONTAINS SUBROUTINE trc_adv_cen2( kt ) INTEGER, INTENT(in) :: kt WRITE(*,*) 'trc_adv_cen2: You should not have seen this print! error?', kt END SUBROUTINE trc_adv_cen2 #endif !!====================================================================== END MODULE trcadv_cen2