New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 1361 for branches/dev_004_VVL/NEMO/OPA_SRC/TRA/tranxt.F90 – NEMO

Ignore:
Timestamp:
2009-03-31T18:26:41+02:00 (15 years ago)
Author:
rblod
Message:

dev004_VVL: new organisation see ticket #389

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/dev_004_VVL/NEMO/OPA_SRC/TRA/tranxt.F90

    r1146 r1361  
    3939   PUBLIC   tra_nxt    ! routine called by step.F90 
    4040 
     41   REAL(wp), DIMENSION(jpk) ::   r2dt_t   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 
     42 
    4143   !! * Substitutions 
    4244#  include "domzgr_substitute.h90" 
     
    6769      !!             at the AGRIF zoom     boundaries (lk_agrif=T) 
    6870      !! 
    69       !!              - Apply the Asselin time filter on now fields, 
    70       !!             save in (ta,sa) an average over the three time levels  
    71       !!             which will be used to compute rdn and thus the semi-implicit 
    72       !!             hydrostatic pressure gradient (ln_dynhpg_imp = T), and 
    73       !!             swap tracer fields to prepare the next time_step. 
    74       !!                This can be summurized for tempearture as: 
    75       !!                    zt = (ta+2tn+tb)/4       ln_dynhpg_imp = T 
    76       !!                    zt = 0                   otherwise 
    77       !!                    tb = tn + atfp*[ tb - 2 tn + ta ] 
    78       !!                    tn = ta  
    79       !!                    ta = zt        (NB: reset to 0 after eos_bn2 call) 
     71      !!              - Update lateral boundary conditions on AGRIF children 
     72      !!             domains (lk_agrif=T) 
    8073      !! 
    8174      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
     
    8780      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    8881      !! 
    89       INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    90       REAL(wp) ::   zt, zs, zfact ! temporary scalars 
     82      INTEGER  ::   jk    ! dummy loop indices 
     83      REAL(wp) ::   zfact ! temporary scalars 
    9184      !!---------------------------------------------------------------------- 
    9285 
     
    111104      CALL Agrif_tra                   ! AGRIF zoom boundaries 
    112105#endif 
     106  
     107      ! set time step size (Euler/Leapfrog) 
     108      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler) 
     109      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog) 
     110      ENDIF 
    113111 
    114112      ! trends computation initialisation 
     
    118116      ENDIF 
    119117 
    120       ! Asselin time filter and swap of arrays 
    121       !                                              
    122       IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler 1st time step : swap only 
    123          DO jk = 1, jpkm1 
    124             tb(:,:,jk) = tn(:,:,jk)                       ! ta, sa remain at their values which 
    125             sb(:,:,jk) = sn(:,:,jk)                       ! correspond to tn, sn after the sawp 
    126             tn(:,:,jk) = ta(:,:,jk)                      
    127             sn(:,:,jk) = sa(:,:,jk) 
    128          END DO 
    129          !                                            
    130       ELSE                                           ! Leap-Frog : filter + swap 
    131          !                                             
    132          IF( ln_dynhpg_imp ) THEN                         ! semi-implicite hpg case 
    133             DO jk = 1, jpkm1                              ! (save the averaged of the 3 time steps  
    134                DO jj = 1, jpj                             !                   in the after fields) 
    135                   DO ji = 1, jpi 
    136                      zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25 
    137                      zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25 
    138                      tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 
    139                      sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 
    140                      tn(ji,jj,jk) = ta(ji,jj,jk) 
    141                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    142                      ta(ji,jj,jk) = zt 
    143                      sa(ji,jj,jk) = zs 
    144                   END DO 
    145                END DO 
    146             END DO 
    147          ELSE                                            ! explicit hpg case 
    148             DO jk = 1, jpkm1 
    149                DO jj = 1, jpj 
    150                   DO ji = 1, jpi 
    151                      tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 
    152                      sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 
    153                      tn(ji,jj,jk) = ta(ji,jj,jk) 
    154                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    155                   END DO 
    156                END DO 
    157             END DO 
    158          ENDIF 
    159          ! 
     118      ! Leap-Frog + Asselin filter time stepping 
     119      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt )      ! variable volume level (vvl) 
     120      ELSE                  ;   CALL tra_nxt_fix( kt )      ! fixed    volume level 
    160121      ENDIF 
    161122 
     
    168129      IF( l_trdtra ) THEN      
    169130         DO jk = 1, jpkm1 
    170             zfact = 1.e0 / ( 2.*rdttra(jk) )             ! NB: euler case, (tb filtered - tb)=0 so 2dt always OK 
     131            zfact = 1.e0 / r2dt_t(jk)              
    171132            ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact 
    172133            ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact 
     
    180141   END SUBROUTINE tra_nxt 
    181142 
     143 
     144   SUBROUTINE tra_nxt_fix( kt ) 
     145      !!---------------------------------------------------------------------- 
     146      !!                   ***  ROUTINE tra_nxt_fix  *** 
     147      !! 
     148      !! ** Purpose :   fixed volume: apply the Asselin time filter and  
     149      !!                swap the tracer fields. 
     150      !!  
     151      !! ** Method  : - Apply a Asselin time filter on now fields. 
     152      !!              - save in (ta,sa) an average over the three time levels  
     153      !!             which will be used to compute rdn and thus the semi-implicit 
     154      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T) 
     155      !!              - swap tracer fields to prepare the next time_step. 
     156      !!                This can be summurized for tempearture as: 
     157      !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T 
     158      !!             ztm = 0                   otherwise 
     159      !!             tb  = tn + atfp*[ tb - 2 tn + ta ] 
     160      !!             tn  = ta  
     161      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call) 
     162      !! 
     163      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
     164      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     165      !!---------------------------------------------------------------------- 
     166      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     167      !! 
     168      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     169      REAL(wp) ::   ztm, ztf     ! temporary scalars 
     170      REAL(wp) ::   zsm, zsf     !    -         - 
     171      !!---------------------------------------------------------------------- 
     172 
     173      IF( kt == nit000 ) THEN 
     174         IF(lwp) WRITE(numout,*) 
     175         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping' 
     176         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     177      ENDIF 
     178      ! 
     179      !                                              ! ----------------------- ! 
     180      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case ! 
     181         !                                           ! ----------------------- ! 
     182         ! 
     183         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step 
     184            DO jk = 1, jpkm1 
     185               DO jj = 1, jpj 
     186                  DO ji = 1, jpi 
     187                     ztm = 0.25 * ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) )      ! mean t 
     188                     zsm = 0.25 * ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) 
     189                     tb(ji,jj,jk) = tn(ji,jj,jk)                                           ! tb <-- tn 
     190                     sb(ji,jj,jk) = sn(ji,jj,jk) 
     191                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tb <-- tn 
     192                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     193                     ta(ji,jj,jk) = ztm                                                    ! ta <-- mean t 
     194                     sa(ji,jj,jk) = zsm 
     195                  END DO 
     196               END DO 
     197            END DO 
     198         ELSE                                             ! general case (Leapfrog + Asselin filter 
     199            DO jk = 1, jpkm1 
     200               DO jj = 1, jpj 
     201                  DO ji = 1, jpi 
     202                     ztm = 0.25 * ( ta(ji,jj,jk) + 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! mean t 
     203                     zsm = 0.25 * ( sa(ji,jj,jk) + 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 
     204                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t  
     205                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 
     206                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn  
     207                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 
     208                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta 
     209                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     210                     ta(ji,jj,jk) = ztm                                                    ! ta <-- mean t 
     211                     sa(ji,jj,jk) = zsm 
     212                  END DO 
     213               END DO 
     214            END DO 
     215         ENDIF 
     216         !                                           ! ----------------------- ! 
     217      ELSE                                           !    explicit hpg case    ! 
     218         !                                           ! ----------------------- ! 
     219         ! 
     220         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step 
     221            DO jk = 1, jpkm1 
     222               DO jj = 1, jpj 
     223                  DO ji = 1, jpi 
     224                     tb(ji,jj,jk) = tn(ji,jj,jk)                                           ! tb <-- tn 
     225                     sb(ji,jj,jk) = sn(ji,jj,jk) 
     226                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta 
     227                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     228                  END DO 
     229               END DO 
     230            END DO 
     231         ELSE                                             ! general case (Leapfrog + Asselin filter 
     232            DO jk = 1, jpkm1 
     233               DO jj = 1, jpj 
     234                  DO ji = 1, jpi 
     235!RBvvl for reproducibility 
     236!                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t  
     237!                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 
     238!                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn  
     239!                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 
     240                     tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 
     241                     sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 
     242                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta 
     243                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     244                  END DO 
     245               END DO 
     246            END DO 
     247         ENDIF 
     248      ENDIF 
     249      ! 
     250   END SUBROUTINE tra_nxt_fix 
     251 
     252 
     253   SUBROUTINE tra_nxt_vvl( kt ) 
     254      !!---------------------------------------------------------------------- 
     255      !!                   ***  ROUTINE tra_nxt_vvl  *** 
     256      !! 
     257      !! ** Purpose :   Time varying volume: apply the Asselin time filter   
     258      !!                and swap the tracer fields. 
     259      !!  
     260      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields. 
     261      !!              - save in (ta,sa) a thickness weighted average over the three  
     262      !!             time levels which will be used to compute rdn and thus the semi- 
     263      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 
     264      !!              - swap tracer fields to prepare the next time_step. 
     265      !!                This can be summurized for tempearture as: 
     266      !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T 
     267      !!                  /(e3t_a   +2*e3t_n   +e3t_b   )     
     268      !!             ztm = 0                                otherwise 
     269      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 
     270      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] ) 
     271      !!             tn  = ta  
     272      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call) 
     273      !! 
     274      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
     275      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     276      !!---------------------------------------------------------------------- 
     277      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     278      !!      
     279       
     280      ! Not yet ready 
     281      WRITE(*,*) 'tra_next_vvl : you should not be there' 
     282      STOP 
     283      ! 
     284   END SUBROUTINE tra_nxt_vvl 
     285 
    182286   !!====================================================================== 
    183287END MODULE tranxt 
Note: See TracChangeset for help on using the changeset viewer.