MODULE tranxt !!====================================================================== !! *** MODULE tranxt *** !! Ocean active tracers: time stepping on temperature and salinity !!====================================================================== !! History : OPA ! 1991-11 (G. Madec) Original code !! 7.0 ! 1993-03 (M. Guyon) symetrical conditions !! 8.0 ! 1996-02 (G. Madec & M. Imbard) opa release 8.0 !! - ! 1996-04 (A. Weaver) Euler forward step !! 8.2 ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure grad. !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module !! - ! 2002-11 (C. Talandier, A-M Treguier) Open boundaries !! - ! 2005-04 (C. Deltel) Add Asselin trend in the ML budget !! 2.0 ! 2006-02 (L. Debreu, C. Mazauric) Agrif implementation !! 3.0 ! 2008-06 (G. Madec) time stepping always done in trazdf !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! tra_nxt : time stepping on temperature and salinity !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers variables USE dom_oce ! ocean space and time domain variables USE zdf_oce ! ??? USE dynspg_oce ! surface pressure gradient variables USE trdmod_oce ! ocean variables trends USE trdmod ! ocean active tracers trends USE phycst USE obctra ! open boundary condition (obc_tra routine) USE bdytra ! Unstructured open boundary condition (bdy_tra routine) USE in_out_manager ! I/O manager USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control USE agrif_opa_update USE agrif_opa_interp IMPLICIT NONE PRIVATE PUBLIC tra_nxt ! routine called by step.F90 REAL(wp), DIMENSION(jpk) :: r2dt_t ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) !! * Substitutions # include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) !! $Id$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_nxt( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE tranxt *** !! !! ** Purpose : Apply the boundary condition on the after temperature !! and salinity fields, achieved the time stepping by adding !! the Asselin filter on now fields and swapping the fields. !! !! ** Method : At this stage of the computation, ta and sa are the !! after temperature and salinity as the time stepping has !! been performed in trazdf_imp or trazdf_exp module. !! !! - Apply lateral boundary conditions on (ta,sa) !! at the local domain boundaries through lbc_lnk call, !! at the radiative open boundaries (lk_obc=T), !! at the relaxed open boundaries (lk_bdy=T), and !! at the AGRIF zoom boundaries (lk_agrif=T) !! !! - Update lateral boundary conditions on AGRIF children !! domains (lk_agrif=T) !! !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) !!---------------------------------------------------------------------- USE oce, ONLY : ztrdt => ua ! use ua as 3D workspace USE oce, ONLY : ztrds => va ! use va as 3D workspace !! INTEGER, INTENT(in) :: kt ! ocean time-step index !! INTEGER :: jk ! dummy loop indices REAL(wp) :: zfact ! temporary scalars !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' IF(lwp) WRITE(numout,*) '~~~~~~~' ENDIF ! Update after tracer on domain lateral boundaries ! CALL lbc_lnk( ta, 'T', 1. ) ! local domain boundaries (T-point, unchanged sign) CALL lbc_lnk( sa, 'T', 1. ) ! #if defined key_obc CALL obc_tra( kt ) ! OBC open boundaries #endif #if defined key_bdy CALL bdy_tra( kt ) ! BDY open boundaries #endif #if defined key_agrif CALL Agrif_tra ! AGRIF zoom boundaries #endif ! set time step size (Euler/Leapfrog) IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt_t(:) = rdttra(:) ! at nit000 (Euler) ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt_t(:) = 2.* rdttra(:) ! at nit000 or nit000+1 (Leapfrog) ENDIF ! trends computation initialisation IF( l_trdtra ) THEN ! store now fields before applying the Asselin filter ztrdt(:,:,:) = tn(:,:,:) ztrds(:,:,:) = sn(:,:,:) ENDIF ! Leap-Frog + Asselin filter time stepping IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt ) ! variable volume level (vvl) ELSE ; CALL tra_nxt_fix( kt ) ! fixed volume level ENDIF #if defined key_agrif ! Update tracer at AGRIF zoom boundaries IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tra( kt ) ! children only #endif ! trends diagnostics : Asselin filter trend : (tb filtered - tb)/2dt IF( l_trdtra ) THEN DO jk = 1, jpkm1 zfact = 1.e0 / r2dt_t(jk) ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact END DO CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt ) END IF ! ! control print IF(ln_ctl) CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt - Tn: ', mask1=tmask, & & tab3d_2=sn, clinfo2= ' Sn: ', mask2=tmask ) ! END SUBROUTINE tra_nxt SUBROUTINE tra_nxt_fix( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_nxt_fix *** !! !! ** Purpose : fixed volume: apply the Asselin time filter and !! swap the tracer fields. !! !! ** Method : - Apply a Asselin time filter on now fields. !! - save in (ta,sa) an average over the three time levels !! which will be used to compute rdn and thus the semi-implicit !! hydrostatic pressure gradient (ln_dynhpg_imp = T) !! - swap tracer fields to prepare the next time_step. !! This can be summurized for tempearture as: !! ztm = (ta+2tn+tb)/4 ln_dynhpg_imp = T !! ztm = 0 otherwise !! tb = tn + atfp*[ tb - 2 tn + ta ] !! tn = ta !! ta = ztm (NB: reset to 0 after eos_bn2 call) !! !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! ocean time-step index !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: ztm, ztf ! temporary scalars REAL(wp) :: zsm, zsf ! - - !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ENDIF ! ! ! ----------------------- ! IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! ! ! ----------------------- ! ! IF( neuler == 0 .AND. kt == nit000 ) THEN ! case of Euler time-stepping at first time-step DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ztm = 0.25 * ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) ! mean t zsm = 0.25 * ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) tb(ji,jj,jk) = tn(ji,jj,jk) ! tb <-- tn sb(ji,jj,jk) = sn(ji,jj,jk) tn(ji,jj,jk) = ta(ji,jj,jk) ! tb <-- tn sn(ji,jj,jk) = sa(ji,jj,jk) ta(ji,jj,jk) = ztm ! ta <-- mean t sa(ji,jj,jk) = zsm END DO END DO END DO ELSE ! general case (Leapfrog + Asselin filter DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ztm = 0.25 * ( ta(ji,jj,jk) + 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! mean t zsm = 0.25 * ( sa(ji,jj,jk) + 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn sb(ji,jj,jk) = sn(ji,jj,jk) + zsf tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta sn(ji,jj,jk) = sa(ji,jj,jk) ta(ji,jj,jk) = ztm ! ta <-- mean t sa(ji,jj,jk) = zsm END DO END DO END DO ENDIF ! ! ----------------------- ! ELSE ! explicit hpg case ! ! ! ----------------------- ! ! IF( neuler == 0 .AND. kt == nit000 ) THEN ! case of Euler time-stepping at first time-step DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi tb(ji,jj,jk) = tn(ji,jj,jk) ! tb <-- tn sb(ji,jj,jk) = sn(ji,jj,jk) tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta sn(ji,jj,jk) = sa(ji,jj,jk) END DO END DO END DO ELSE ! general case (Leapfrog + Asselin filter DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi !RBvvl for reproducibility ! ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t ! zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) ! tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn ! sb(ji,jj,jk) = sn(ji,jj,jk) + zsf tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta sn(ji,jj,jk) = sa(ji,jj,jk) END DO END DO END DO ENDIF ENDIF ! END SUBROUTINE tra_nxt_fix SUBROUTINE tra_nxt_vvl( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_nxt_vvl *** !! !! ** Purpose : Time varying volume: apply the Asselin time filter !! and swap the tracer fields. !! !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. !! - save in (ta,sa) a thickness weighted average over the three !! time levels which will be used to compute rdn and thus the semi- !! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) !! - swap tracer fields to prepare the next time_step. !! This can be summurized for tempearture as: !! ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb) ln_dynhpg_imp = T !! /(e3t_a +2*e3t_n +e3t_b ) !! ztm = 0 otherwise !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) !! tn = ta !! ta = zt (NB: reset to 0 after eos_bn2 call) !! !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! ocean time-step index !! ! Not yet ready WRITE(*,*) 'tra_next_vvl : you should not be there' STOP ! END SUBROUTINE tra_nxt_vvl !!====================================================================== END MODULE tranxt