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 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) semi-implicit hpg with asselin filter + modified LF-RA !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! 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 sbc_oce ! surface boundary condition: ocean 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 traqsr ! penetrative solar radiation (needed for nksr) USE agrif_opa_update USE agrif_opa_interp IMPLICIT NONE PRIVATE PUBLIC tra_nxt ! routine called by step.F90 REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg REAL(wp), DIMENSION(jpk) :: r2dt_t ! 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) !! !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. !!---------------------------------------------------------------------- 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 !== initialisation ==! IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' IF(lwp) WRITE(numout,*) '~~~~~~~' ! rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp) ! Brown & Campana parameter for semi-implicit hpg 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 ! !== 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 IF( l_trdtra ) THEN ! trends computation: store now fields before applying the Asselin filter ztrdt(:,:,:) = tn(:,:,:) ztrds(:,:,:) = sn(:,:,:) ENDIF ! !== modifed Leap-Frog + Asselin filter (modified LF-RA) scheme ==! 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 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tra( kt ) ! Update tracer at AGRIF zoom boundaries (children only) #endif IF( l_trdtra ) THEN ! trends computation: 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 temperature as: !! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T !! ztm = 0 otherwise !! with rbcp=1/4 * (1-atfp^4) / (1-atfp) !! 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) :: zt_d, zs_d ! temporary scalars REAL(wp) :: ztn, zsn ! - - !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ENDIF ! IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step DO jk = 1, jpkm1 ! (only swap) tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta sn(:,:,jk) = sa(:,:,jk) ! sn <-- sa END DO ELSE ! General case (Leapfrog + Asselin filter) DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi IF( ln_dynhpg_imp ) THEN ! implicit hpg: keep tn, sn in memory ztn = tn(ji,jj,jk) zsn = sn(ji,jj,jk) ENDIF ! ! time laplacian on tracers ! ! used for both Asselin and Brown & Campana filters zt_d = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk) zs_d = sa(ji,jj,jk) - 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ! ! ! swap of arrays tb(ji,jj,jk) = tn(ji,jj,jk) + atfp * zt_d ! tb <-- tn filtered sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_d ! sb <-- sn filtered tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta sn(ji,jj,jk) = sa(ji,jj,jk) ! sn <-- sa ! ! semi imlicit hpg computation (Brown & Campana) IF( ln_dynhpg_imp ) THEN ta(ji,jj,jk) = ztn + rbcp * zt_d ! ta <-- Brown & Campana average sa(ji,jj,jk) = zsn + rbcp * zs_d ! sa <-- Brown & Campana average ENDIF END DO END DO END DO 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 temperature as: !! ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) ln_dynhpg_imp = T !! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] ) !! 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 :: ze3t_a, ze3t_n, ze3t_b ! temporary scalars REAL :: ztc_a, ztc_n, ztc_b ! - - REAL :: zsc_a, zsc_n, zsc_b ! - - REAL :: ztc_f, zsc_f, ztc_d, zsc_d ! - - REAL :: ze3t_f, ze3t_d ! - - REAL :: zfact1, zfact2 ! - - !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ENDIF IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step ! ! (only swap) DO jk = 1, jpkm1 ! tn <-- ta tn(:,:,jk) = ta(:,:,jk) ! sn <-- sa sn(:,:,jk) = sa(:,:,jk) END DO ! ! General case (Leapfrog + Asselin filter) ELSE ! apply filter on thickness weighted tracer and swap DO jk = 1, jpkm1 zfact1 = atfp * rdttra(jk) zfact2 = zfact1 / rau0 DO jj = 1, jpj DO ji = 1, jpi ! ! scale factors at Before, now and after ze3t_b = fse3t_b(ji,jj,jk) ze3t_n = fse3t_n(ji,jj,jk) ze3t_a = fse3t_a(ji,jj,jk) ze3t_d = fse3t_d(ji,jj,jk) ! ! tracer content at Before, now and after ztc_b = tb(ji,jj,jk) * ze3t_b ; zsc_b = sb(ji,jj,jk) * ze3t_b ztc_n = tn(ji,jj,jk) * ze3t_n ; zsc_n = sn(ji,jj,jk) * ze3t_n ztc_a = ta(ji,jj,jk) * ze3t_a ; zsc_a = sa(ji,jj,jk) * ze3t_a ! ! ! Time laplacian on tracer contents ! ! used for both Asselin and Brown & Campana filters ztc_d = ztc_a - 2. * ztc_n + ztc_b zsc_d = zsc_a - 2. * zsc_n + zsc_b ! ! Asselin Filter on thicknesses and tracer contents ze3t_f = ze3t_n + atfp * ze3t_d ztc_f = ztc_n + atfp * ztc_d zsc_f = zsc_n + atfp * zsc_d ! ! Filter correction IF( jk == 1 ) THEN ! WRITE(numout,*) 'filter correction: sbc_trd_hc_n' ze3t_f = ze3t_f - zfact2 * ( emp_b (ji,jj) - emp (ji,jj) ) ztc_f = ztc_f - zfact1 * ( sbc_trd_hc_n(ji,jj) - sbc_trd_hc_b(ji,jj) ) ENDIF IF( ln_traqsr .AND. ( jk .LE. nksr ) ) THEN ! WRITE(numout,*) 'jk =', jk ! WRITE(numout,*) 'filter correction: qsr_trd_hc_n' ztc_f = ztc_f - zfact1 * ( qsr_trd_hc_n(ji,jj,jk) - qsr_trd_hc_b(ji,jj,jk) ) ENDIF ! swap of arrays ze3t_f = 1.e0 / ze3t_f tb(ji,jj,jk) = ztc_f * ze3t_f ! tb <-- tn filtered sb(ji,jj,jk) = zsc_f * ze3t_f ! sb <-- sn filtered tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta sn(ji,jj,jk) = sa(ji,jj,jk) ! sn <-- sa ! ! semi imlicit hpg computation (Brown & Campana) IF( ln_dynhpg_imp ) THEN ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d ) ta(ji,jj,jk) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average sa(ji,jj,jk) = ze3t_d * ( zsc_n + rbcp * zsc_d ) ! sa <-- Brown & Campana average ENDIF END DO END DO END DO ENDIF ! END SUBROUTINE tra_nxt_vvl !!====================================================================== END MODULE tranxt