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 !! 3.1 ! 2009-02 (G. Madec, R. Benshila) re-introduce the vvl option !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! tra_nxt : time stepping on temperature and salinity !! tra_nxt_fix : time stepping on temperature and salinity : fixed volume case !! tra_nxt_vvl : time stepping on temperature and salinity : variable volume case !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers variables USE dom_oce ! ocean space and time domain variables USE zdf_oce ! ??? USE domvvl ! variable volume USE dynspg_oce ! surface pressure gradient variables USE dynhpg ! hydrostatic pressure gradient USE trdmod_oce ! ocean space and time domain variables USE trdtra ! 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 traswp ! swap array USE agrif_opa_update USE agrif_opa_interp USE obc_oce IMPLICIT NONE PRIVATE PUBLIC tra_nxt ! routine called by step.F90 PUBLIC tra_nxt_fix ! to be used in trcnxt PUBLIC tra_nxt_vvl ! to be used in trcnxt REAL(wp), DIMENSION(jpk) :: r2dt ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) !! * Substitutions # include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) !! $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) !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! ocean time-step index !! INTEGER :: jk ! dummy loop indices REAL(wp) :: zfact ! temporary scalars REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds !!---------------------------------------------------------------------- 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( tsa(:,:,:,jp_tem), 'T', 1. ) ! local domain boundaries (T-point, unchanged sign) CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! #if defined key_obc || defined key_bdy || defined key_agrif CALL tra_unswap #endif #if defined key_obc IF( lk_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 #if defined key_obc || defined key_bdy || defined key_agrif CALL tra_swap #endif ! set time step size (Euler/Leapfrog) IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt(:) = rdttra(:) ! at nit000 (Euler) ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt(:) = 2.* rdttra(:) ! at nit000 or nit000+1 (Leapfrog) ENDIF ! trends computation initialisation IF( l_trdtra ) THEN !* store now fields before applying the Asselin filter ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsn(:,:,:,jp_tem) ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsn(:,:,:,jp_sal) ENDIF ! Leap-Frog + Asselin filter time stepping IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, nit000, tsb, tsn, tsa, jpts ) ! variable volume level (vvl) ELSE ; CALL tra_nxt_fix( kt, nit000, tsb, tsn, tsa, jpts ) ! fixed volume level ENDIF #if defined key_agrif CALL tra_unswap ! Update tracer at AGRIF zoom boundaries IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tra( kt ) ! children only CALL tra_swap #endif ! trends computation IF( l_trdtra ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt DO jk = 1, jpkm1 zfact = 1.e0 / r2dt(jk) ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact END DO CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_atf, ztrdt ) CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_atf, ztrds ) DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) END IF ! ! control print IF(ln_ctl) CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt - Tn: ', mask1=tmask, & & tab3d_2=tsn(:,:,:,jp_sal), clinfo2= ' Sn: ', mask2=tmask ) ! END SUBROUTINE tra_nxt SUBROUTINE tra_nxt_fix( kt, kit000, & & ptb, ptn, pta, kjpt ) !!---------------------------------------------------------------------- !! *** 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 , INTENT(in ) :: kit000 ! first time-step index INTEGER , INTENT(in ) :: kjpt ! number of tracers REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend !! INTEGER :: ji, jj, jk, jn ! dummy loop indices REAL(wp) :: ztm, ztf ! temporary scalars !!---------------------------------------------------------------------- IF( kt == kit000 ) 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 == kit000 ) THEN ! Euler time-stepping at first time-step ! ! (only swap) DO jn = 1, kjpt DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptb <-- ptn END DO END DO END DO END DO ELSE ! general case (Leapfrog + Asselin filter DO jn = 1, kjpt DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ztm = 0.25 * ( pta(ji,jj,jk,jn) + 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ) ! mean pt ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptn(ji,jj,jk,jn) ) ! Asselin filter on pt ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf ! ptb <-- filtered ptn ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta pta(ji,jj,jk,jn) = ztm ! pta <-- mean pt END DO END DO END DO END DO ENDIF ! ! ----------------------- ! ELSE ! explicit hpg case ! ! ! ----------------------- ! ! IF( neuler == 0 .AND. kt == kit000 ) THEN ! Euler time-stepping at first time-step DO jn = 1, kjpt DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta END DO END DO END DO END DO ELSE ! general case (Leapfrog + Asselin filter DO jn = 1, kjpt DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ) ! Asselin filter on t ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf ! ptb <-- filtered ptn ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta END DO END DO END DO END DO ENDIF ! ENDIF ! END SUBROUTINE tra_nxt_fix SUBROUTINE tra_nxt_vvl( kt, kit000, & & ptb, ptn, pta, kjpt ) !!---------------------------------------------------------------------- !! *** 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 INTEGER , INTENT(in ) :: kit000 ! first time-step index INTEGER , INTENT(in ) :: kjpt ! number of tracers REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend !! INTEGER :: ji, jj, jk, jn ! dummy loop indices REAL(wp) :: ztm , ztc_f , ztf , ztca, ztcn, ztcb ! temporary scalar REAL(wp) :: ze3mr, ze3fr ! - - REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f ! - - !!---------------------------------------------------------------------- IF( kt == kit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ENDIF ! ! ! ----------------------- ! IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! ! ! ----------------------- ! ! IF( neuler == 0 .AND. kt == kit000 ) THEN ! Euler time-stepping at first time-step DO jn = 1, kjpt ! (only swap) DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! tn <-- ta END DO END DO END DO END DO ELSE DO jn = 1, kjpt ! (only swap) DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ze3t_b = fse3t_b(ji,jj,jk) ze3t_n = fse3t_n(ji,jj,jk) ze3t_a = fse3t_a(ji,jj,jk) ! ! tracer content at Before, now and after ztcb = ptb(ji,jj,jk,jn) * ze3t_b ztcn = ptn(ji,jj,jk,jn) * ze3t_n ztca = pta(ji,jj,jk,jn) * ze3t_a ! ! ! Asselin filter on thickness and tracer content ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) ! ! ! filtered tracer including the correction ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) ztf = ze3fr * ( ztcn + ztc_f ) ! ! mean thickness and tracer ze3mr = 1.e0 / ( ze3t_a + 2.* ze3t_n + ze3t_b ) ztm = ze3mr * ( ztca + 2.* ztcn + ztcb ) !!gm mean e3t have to be saved and used in dynhpg or it can be recomputed in dynhpg !! !!gm e3t_m(ji,jj,jk) = 0.25 / ze3mr ! ! swap of arrays ptb(ji,jj,jk,jn) = ztf ! ptb <-- ptn + filter ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta pta(ji,jj,jk,jn) = ztm ! pta <-- mean t END DO END DO END DO END DO ENDIF ! ! ----------------------- ! ELSE ! explicit hpg case ! ! ! ----------------------- ! ! IF( neuler == 0 .AND. kt == kit000 ) THEN ! case of Euler time-stepping at first time-step DO jn = 1, kjpt ! No filter nor thickness weighting computation required DO jk = 1, jpkm1 ! ONLY swap DO jj = 1, jpj DO ji = 1, jpi ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! tn <-- ta END DO END DO END DO END DO ! ! general case (Leapfrog + Asselin filter) ELSE ! apply filter on thickness weighted tracer and swap DO jn = 1, kjpt DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ze3t_b = fse3t_b(ji,jj,jk) ze3t_n = fse3t_n(ji,jj,jk) ze3t_a = fse3t_a(ji,jj,jk) ! ! tracer content at Before, now and after ztcb = ptb(ji,jj,jk,jn) * ze3t_b ztcn = ptn(ji,jj,jk,jn) * ze3t_n ztca = pta(ji,jj,jk,jn) * ze3t_a ! ! ! Asselin filter on thickness and tracer content ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) ! ! ! filtered tracer including the correction ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) ztf = ( ztcn + ztc_f ) * ze3fr ! ! swap of arrays ptb(ji,jj,jk,jn) = ztf ! tb <-- tn filtered ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! tn <-- ta END DO END DO END DO END DO ENDIF ENDIF ! END SUBROUTINE tra_nxt_vvl !!====================================================================== END MODULE tranxt