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 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: tranxt.F90 1601 2009-08-11 10:09:19Z ctlod $ !! 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 computation IF( l_trdtra ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 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 ! Euler time-stepping at first time-step DO jk = 1, jpkm1 ! (only swap) DO jj = 1, jpj DO ji = 1, jpi tn(ji,jj,jk) = ta(ji,jj,jk) ! tb <-- tn 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 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 ! Euler time-stepping at first time-step DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi 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 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) 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 !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: ztm , ztc_f , ztf , ztca, ztcn, ztcb ! temporary scalar REAL(wp) :: zsm , zsc_f , zsf , zsca, zscn, zscb ! - - REAL(wp) :: ze3mr, ze3fr ! - - REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f ! - - !!---------------------------------------------------------------------- IF( kt == nit000 ) 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 == nit000 ) THEN ! Euler time-stepping at first time-step DO jk = 1, jpkm1 ! (only swap) DO jj = 1, jpj DO ji = 1, jpi 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 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 = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * 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 ) zsc_f = atfp * ( zsca - 2.* zscn + zscb ) ! ! ! filtered tracer including the correction ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) ztf = ze3fr * ( ztcn + ztc_f ) zsf = ze3fr * ( zscn + zsc_f ) ! ! mean thickness and tracer ze3mr = 1.e0 / ( ze3t_a + 2.* ze3t_n + ze3t_b ) ztm = ze3mr * ( ztca + 2.* ztcn + ztcb ) zsm = ze3mr * ( zsca + 2.* zscn + zscb ) !!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 tb(ji,jj,jk) = ztf ! tb <-- tn + filter sb(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 ! No filter nor thickness weighting computation required DO jj = 1, jpj ! ONLY swap DO ji = 1, jpi tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta sn(ji,jj,jk) = sa(ji,jj,jk) END DO END DO END DO ! ! general case (Leapfrog + Asselin filter) ELSE ! apply filter on thickness weighted tracer and 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 = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * 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 ) zsc_f = atfp * ( zsca - 2.* zscn + zscb ) ! ! ! filtered tracer including the correction ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) ztf = ( ztcn + ztc_f ) * ze3fr zsf = ( zscn + zsc_f ) * ze3fr ! ! swap of arrays tb(ji,jj,jk) = ztf ! tb <-- tn filtered sb(ji,jj,jk) = zsf 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_vvl !!====================================================================== END MODULE tranxt