Changeset 1110


Ignore:
Timestamp:
2008-06-13T16:10:35+02:00 (12 years ago)
Author:
ctlod
Message:

trunk: remove a bug related to the combination of key_vvl and ln_trazdf_exp options, see ticket: #202

Location:
trunk/NEMO/OPA_SRC/TRA
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/TRA/tranxt.F90

    r911 r1110  
    44   !! Ocean active tracers:  time stepping on temperature and salinity 
    55   !!====================================================================== 
    6    !! History :  7.0  !  91-11  (G. Madec)  Original code 
    7    !!                 !  93-03  (M. Guyon)  symetrical conditions 
    8    !!                 !  96-02  (G. Madec & M. Imbard)  opa release 8.0 
    9    !!            8.0  !  96-04  (A. Weaver)  Euler forward step 
    10    !!            8.2  !  99-02  (G. Madec, N. Grima)  semi-implicit pressure grad. 
    11    !!            8.5  !  02-08  (G. Madec)  F90: Free form and module 
    12    !!                 !  02-11  (C. Talandier, A-M Treguier) Open boundaries 
    13    !!                 !  05-04  (C. Deltel) Add Asselin trend in the ML budget 
    14    !!            9.0  !  06-02  (L. Debreu, C. Mazauric) Agrif implementation 
     6   !! History :  OPA  !  1991-11  (G. Madec)  Original code 
     7   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions 
     8   !!            8.0  !  1996-02  (G. Madec & M. Imbard)  opa release 8.0 
     9   !!             -   !  1996-04  (A. Weaver)  Euler forward step 
     10   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad. 
     11   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module 
     12   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries 
     13   !!             -   !  2005-04  (C. Deltel) Add Asselin trend in the ML budget 
     14   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation 
     15   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf 
    1516   !!---------------------------------------------------------------------- 
    1617 
    1718   !!---------------------------------------------------------------------- 
    18    !!   tra_nxt     : time stepping on temperature and salinity 
     19   !!   tra_nxt       : time stepping on temperature and salinity 
    1920   !!---------------------------------------------------------------------- 
    2021   USE oce             ! ocean dynamics and tracers variables 
     
    2526   USE trdmod          ! ocean active tracers trends  
    2627   USE phycst 
    27    USE domvvl          ! variable volume 
    2828   USE obctra          ! open boundary condition (obc_tra routine) 
    2929   USE bdytra          ! Unstructured open boundary condition (bdy_tra routine) 
     
    3434   USE agrif_opa_interp 
    3535 
    36  
    3736   IMPLICIT NONE 
    3837   PRIVATE 
    3938 
    40    !! * Routine accessibility 
    41    PUBLIC   tra_nxt   ! routine called by step.F90 
    42  
    43    REAL(wp) ::   vemp ! total amount of volume added or removed by E-P forcing 
     39   PUBLIC   tra_nxt    ! routine called by step.F90 
    4440 
    4541   !! * Substitutions 
    4642#  include "domzgr_substitute.h90" 
    4743   !!---------------------------------------------------------------------- 
    48    !!   OPA 9.0 , LOCEAN-IPSL (2006)  
    49    !! $Id$ 
     44   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     45   !! $Id:$ 
    5046   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5147   !!---------------------------------------------------------------------- 
     
    5753      !!                   ***  ROUTINE tranxt  *** 
    5854      !! 
    59       !! ** Purpose :   Compute the temperature and salinity fields at the  
    60       !!      next time-step from their temporal trends and swap the fields. 
     55      !! ** Purpose :   Apply the boundary condition on the after temperature   
     56      !!             and salinity fields, achieved the time stepping by adding 
     57      !!             the Asselin filter on now fields and swapping the fields. 
    6158      !!  
    62       !! ** Method  :   Apply lateral boundary conditions on (ua,va) through  
    63       !!      call to lbc_lnk routine 
    64       !!      After t and s are compute using a leap-frog scheme environment: 
    65       !!         ta = tb + 2 rdttra(k) * ta 
    66       !!         sa = sb + 2 rdttra(k) * sa 
    67       !!      Compute and save in (ta,sa) an average over three time levels 
    68       !!      (before,now and after) of temperature and salinity which is 
    69       !!      used to compute rhd in eos routine and thus the hydrostatic  
    70       !!      pressure gradient (ln_dynhpg_imp = T) 
    71       !!      Apply an Asselin time filter on now tracers (tn,sn) to avoid 
    72       !!      the divergence of two consecutive time-steps and swap tracer 
    73       !!      arrays to prepare the next time_step: 
    74       !!         (zt,zs) = (ta+2tn+tb,sa+2sn+sb)/4       (ln_dynhpg_imp = T) 
    75       !!         (zt,zs) = (0,0)                            (default option) 
    76       !!         (tb,sb) = (tn,vn) + atfp [ (tb,sb) + (ta,sa) - 2 (tn,sn) ] 
    77       !!         (tn,sn) = (ta,sa)  
    78       !!         (ta,sa) = (zt,zs)  (NB: reset to 0 after use in eos.F) 
     59      !! ** Method  :   At this stage of the computation, ta and sa are the  
     60      !!             after temperature and salinity as the time stepping has 
     61      !!             been performed in trazdf_imp or trazdf_exp module. 
    7962      !! 
    80       !! ** Action  : - update (tb,sb) and (tn,sn)  
    81       !!              - (ta,sa) time averaged (t,s)      (ln_dynhpg_imp = T) 
     63      !!              - Apply lateral boundary conditions on (ta,sa)  
     64      !!             at the local domain   boundaries through lbc_lnk call,  
     65      !!             at the radiative open boundaries (lk_obc=T),  
     66      !!             at the relaxed   open boundaries (lk_bdy=T), and 
     67      !!             at the AGRIF zoom     boundaries (lk_agrif=T) 
     68      !! 
     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) 
     80      !! 
     81      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
     82      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
    8283      !!---------------------------------------------------------------------- 
    8384      USE oce, ONLY :    ztrdt => ua   ! use ua as 3D workspace    
     
    8788      !! 
    8889      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    89       REAL(wp) ::   zt, zs, zssh1 ! temporary scalars 
    90       REAL(wp) ::   zfact         ! temporary scalar 
    91       !! Variable volume 
    92       REAL(wp), DIMENSION(jpi,jpj) ::   zssh                         ! temporary scalars 
    93       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3tb, zfse3tn, zfse3ta  ! 3D workspace 
    94  
     90      REAL(wp) ::   zt, zs, zfact ! temporary scalars 
    9591      !!---------------------------------------------------------------------- 
    9692 
    97       !! Explicit physics with thickness weighted updates 
    98       IF( lk_vvl .AND. ln_zdfexp ) THEN 
    99  
    100          ! Scale factors at before and after time step 
    101          ! ------------------------------------------- 
    102          CALL dom_vvl_sf( sshb, 'T', zfse3tb ) ;    CALL dom_vvl_sf( ssha, 'T', zfse3ta ) 
    103  
    104          ! Asselin filtered scale factor at now time step 
    105          ! ---------------------------------------------- 
    106          IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN 
    107             CALL dom_vvl_sf_ini( 'T', zfse3tn )  
    108          ELSE 
    109             zssh(:,:) = atfp * ( sshb(:,:) + ssha(:,:) ) + atfp1 * sshn(:,:) 
    110             CALL dom_vvl_sf( zssh, 'T', zfse3tn )  
    111          ENDIF 
    112  
    113          ! Thickness weighting 
    114          ! ------------------- 
    115          DO jk = 1, jpkm1 
    116             DO jj = 1, jpj 
    117                DO ji = 1, jpi 
    118                   ta(ji,jj,jk) = ta(ji,jj,jk) * fse3t(ji,jj,jk) 
    119                   sa(ji,jj,jk) = sa(ji,jj,jk) * fse3t(ji,jj,jk) 
    120  
    121                   tn(ji,jj,jk) = tn(ji,jj,jk) * fse3t(ji,jj,jk) 
    122                   sn(ji,jj,jk) = sn(ji,jj,jk) * fse3t(ji,jj,jk) 
    123  
    124                   tb(ji,jj,jk) = tb(ji,jj,jk) * zfse3tb(ji,jj,jk) 
    125                   sb(ji,jj,jk) = sb(ji,jj,jk) * zfse3tb(ji,jj,jk) 
    126                END DO 
    127             END DO 
    128          END DO 
    129  
     93      IF( kt == nit000 ) THEN 
     94         IF(lwp) WRITE(numout,*) 
     95         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 
     96         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    13097      ENDIF 
    13198 
    132       IF( l_trdtra ) THEN 
    133          ztrdt(:,:,jpk) = 0.e0 
    134          ztrds(:,:,jpk) = 0.e0 
     99      ! Update after tracer on domain lateral boundaries 
     100      !  
     101      CALL lbc_lnk( ta, 'T', 1. )      ! local domain boundaries  (T-point, unchanged sign) 
     102      CALL lbc_lnk( sa, 'T', 1. ) 
     103      ! 
     104#if defined key_obc 
     105      CALL obc_tra( kt )               ! OBC open boundaries 
     106#endif 
     107#if defined key_bdy 
     108      CALL bdy_tra( kt )               ! BDY open boundaries 
     109#endif 
     110#if defined key_agrif 
     111      CALL Agrif_tra                   ! AGRIF zoom boundaries 
     112#endif 
     113 
     114      ! trends computation initialisation 
     115      IF( l_trdtra ) THEN              ! store now fields before applying the Asselin filter 
     116         ztrdt(:,:,:) = tn(:,:,:)       
     117         ztrds(:,:,:) = sn(:,:,:) 
    135118      ENDIF 
    136119 
    137       ! 0. Lateral boundary conditions on ( ta, sa )   (T-point, unchanged sign) 
    138       ! ---------------------------------============ 
    139       CALL lbc_lnk( ta, 'T', 1. )    
    140       CALL lbc_lnk( sa, 'T', 1. ) 
    141  
    142       !                                                ! =============== 
    143       DO jk = 1, jpkm1                                 ! Horizontal slab 
    144          !                                             ! =============== 
    145  
    146          ! 1. Leap-frog scheme (only in explicit case, otherwise the  
    147          ! -------------------  time stepping is already done in trazdf) 
    148          IF( ln_zdfexp ) THEN 
    149             zfact = 2. * rdttra(jk) 
    150             IF( neuler == 0 .AND. kt == nit000 )   zfact = rdttra(jk) 
    151             ta(:,:,jk) = ( tb(:,:,jk) + zfact * ta(:,:,jk) ) * tmask(:,:,jk) 
    152             sa(:,:,jk) = ( sb(:,:,jk) + zfact * sa(:,:,jk) ) * tmask(:,:,jk) 
    153             IF(l_trdtra)   CALL ctl_stop( 'tranxt: Asselin ML trend not yet accounted for.' ) 
    154          ENDIF 
    155  
    156 #if defined key_obc 
    157          !                                             ! =============== 
    158       END DO                                           !   End of slab 
    159       !                                                ! =============== 
    160       ! Update tracers on open boundaries. 
    161       CALL obc_tra( kt ) 
    162  
    163       !                                                ! =============== 
    164       DO jk = 1, jpkm1                                 ! Horizontal slab 
    165          !                                             ! =============== 
    166 #elif defined key_bdy 
    167          !                                             ! =============== 
    168       END DO                                           !   End of slab 
    169       !                                                ! =============== 
    170  
    171       ! Update tracers on open boundaries. 
    172       CALL bdy_tra( kt ) 
    173  
    174       !                                                ! =============== 
    175       DO jk = 1, jpkm1                                 ! Horizontal slab 
    176          !                                             ! =============== 
    177 #endif 
    178 #if defined key_agrif 
    179          !                                             ! =============== 
    180       END DO                                           !   End of slab 
    181       !                                                ! =============== 
    182       ! Update tracers on open boundaries. 
    183       CALL Agrif_tra 
    184       !                                                ! =============== 
    185       DO jk = 1, jpkm1                                 ! Horizontal slab 
    186          !                                             ! =============== 
    187 #endif 
    188          ! 2. Time filter and swap of arrays 
    189          ! --------------------------------- 
    190  
    191          IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg 
    192             IF( neuler == 0 .AND. kt == nit000 ) THEN 
    193                DO jj = 1, jpj 
     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) 
    194135                  DO ji = 1, jpi 
    195136                     zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25 
    196137                     zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25 
    197                      tb(ji,jj,jk) = tn(ji,jj,jk) 
    198                      sb(ji,jj,jk) = sn(ji,jj,jk) 
     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) 
    199140                     tn(ji,jj,jk) = ta(ji,jj,jk) 
    200141                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     
    203144                  END DO 
    204145               END DO 
    205                IF( l_trdtra ) THEN 
    206                   ztrdt(:,:,jk) = 0.e0 
    207                   ztrds(:,:,jk) = 0.e0 
    208                END IF 
    209             ELSE 
     146            END DO 
     147         ELSE                                            ! explicit hpg case 
     148            DO jk = 1, jpkm1 
    210149               DO jj = 1, jpj 
    211150                  DO ji = 1, jpi 
    212                      zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25 
    213                      zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25 
    214                      tb(ji,jj,jk) = atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 
    215                      sb(ji,jj,jk) = atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 
    216                      IF( l_trdtra ) THEN ! ChD ceci est a optimiser, mais ca marche 
    217                         ztrdt(ji,jj,jk) = tb(ji,jj,jk) - tn(ji,jj,jk) 
    218                         ztrds(ji,jj,jk) = sb(ji,jj,jk) - sn(ji,jj,jk) 
    219                      END IF 
     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) 
    220153                     tn(ji,jj,jk) = ta(ji,jj,jk) 
    221154                     sn(ji,jj,jk) = sa(ji,jj,jk) 
    222                      ta(ji,jj,jk) = zt 
    223                      sa(ji,jj,jk) = zs 
    224155                  END DO 
    225156               END DO 
    226             ENDIF 
    227          ELSE                                          ! Default case 
    228             IF( neuler == 0 .AND. kt == nit000 ) THEN 
    229                IF( (lk_vvl .AND. ln_zdfexp) ) THEN                      ! Varying levels 
    230                   DO jj = 1, jpj 
    231                      DO ji = 1, jpi 
    232                         zssh1 = tmask(ji,jj,jk) / fse3t(ji,jj,jk) 
    233                         tb(ji,jj,jk) = tn(ji,jj,jk) * zssh1 * tmask(ji,jj,jk) 
    234                         sb(ji,jj,jk) = sn(ji,jj,jk) * zssh1 * tmask(ji,jj,jk) 
    235                         zssh1 = tmask(ji,jj,jk) / zfse3ta(ji,jj,jk) 
    236                         tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1 * tmask(ji,jj,jk) 
    237                         sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1 * tmask(ji,jj,jk) 
    238                      END DO 
    239                   END DO 
    240                ELSE                                                     ! Fixed levels 
    241                  DO jj = 1, jpj 
    242                      DO ji = 1, jpi 
    243                         tb(ji,jj,jk) = tn(ji,jj,jk) 
    244                         sb(ji,jj,jk) = sn(ji,jj,jk) 
    245                         tn(ji,jj,jk) = ta(ji,jj,jk) 
    246                         sn(ji,jj,jk) = sa(ji,jj,jk) 
    247                      END DO 
    248                   END DO 
    249                ENDIF 
    250                IF( l_trdtra ) THEN 
    251                   ztrdt(:,:,jk) = 0.e0 
    252                   ztrds(:,:,jk) = 0.e0 
    253                END IF 
    254             ELSE 
    255                IF( l_trdtra ) THEN 
    256                   DO jj = 1, jpj 
    257                      DO ji = 1, jpi 
    258                         ztrdt(ji,jj,jk) = atfp * ( tb(ji,jj,jk) - 2*tn(ji,jj,jk) + ta(ji,jj,jk) ) 
    259                         ztrds(ji,jj,jk) = atfp * ( sb(ji,jj,jk) - 2*sn(ji,jj,jk) + sa(ji,jj,jk) ) 
    260                      END DO 
    261                   END DO 
    262                END IF 
    263                IF( (lk_vvl .AND. ln_zdfexp) ) THEN                      ! Varying levels 
    264                   DO jj = 1, jpj 
    265                      DO ji = 1, jpi 
    266                         zssh1 = tmask(ji,jj,jk) / zfse3tn(ji,jj,jk) 
    267                         tb(ji,jj,jk) = ( atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) & 
    268                           &            + atfp1 * tn(ji,jj,jk) ) * zssh1 
    269                         sb(ji,jj,jk) = ( atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) & 
    270                           &            + atfp1 * sn(ji,jj,jk) ) * zssh1 
    271                         zssh1 = tmask(ji,jj,1) / zfse3ta(ji,jj,jk) 
    272                         tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1 
    273                         sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1 
    274                      END DO 
    275                   END DO 
    276                ELSE                                                     ! Fixed levels or first varying level 
    277                   DO jj = 1, jpj 
    278                      DO ji = 1, jpi 
    279                         tb(ji,jj,jk) = atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 
    280                         sb(ji,jj,jk) = atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 
    281                         tn(ji,jj,jk) = ta(ji,jj,jk) 
    282                         sn(ji,jj,jk) = sa(ji,jj,jk) 
    283                      END DO 
    284                   END DO 
    285                ENDIF 
    286             ENDIF 
     157            END DO 
    287158         ENDIF 
    288          !                                             ! =============== 
    289       END DO                                           !   End of slab 
    290       !                                                ! =============== 
     159         ! 
     160      ENDIF 
    291161 
    292       IF( l_trdtra )   THEN      ! Take the Asselin trend into account  
    293          ztrdt(:,:,:) = ztrdt(:,:,:) / ( 2.*rdt ) 
    294          ztrds(:,:,:) = ztrds(:,:,:) / ( 2.*rdt )          
     162#if defined key_agrif 
     163      ! Update tracer at AGRIF zoom boundaries 
     164      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only 
     165#endif       
     166 
     167      ! trends diagnostics :  Asselin filter trend : (tb filtered - tb)/2dt 
     168      IF( l_trdtra ) THEN      
     169         DO jk = 1, jpkm1 
     170            zfact = 1.e0 / ( 2.*rdttra(jk) )             ! NB: euler case, (tb filtered - tb)=0 so 2dt always OK 
     171            ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact 
     172            ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact 
     173         END DO 
    295174         CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt ) 
    296175      END IF 
    297  
     176      !                        ! control print 
    298177      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt  - Tn: ', mask1=tmask,   & 
    299178         &                       tab3d_2=sn, clinfo2=       ' Sn: ', mask2=tmask ) 
    300 #if defined key_agrif 
    301       IF (.NOT.Agrif_Root())    CALL Agrif_Update_Tra( kt ) 
    302 #endif       
    303179      ! 
    304180   END SUBROUTINE tra_nxt 
  • trunk/NEMO/OPA_SRC/TRA/trazdf.F90

    r888 r1110  
    9292      CASE ( 0 )                                       ! explicit scheme 
    9393         CALL tra_zdf_exp    ( kt, r2dt ) 
    94          IF( l_trdtra )   THEN                         ! save the vertical diffusive trends for further diagnostics 
    95             ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:) 
    96             ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:) 
    97             CALL trd_mod( ztrdt, ztrds, jptra_trd_zdf, 'TRA', kt ) 
    98          ENDIF 
    9994 
    10095      CASE ( 1 )                                       ! implicit scheme  
    10196         CALL tra_zdf_imp    ( kt, r2dt ) 
    102          IF( l_trdtra )   THEN                         ! save the vertical diffusive trends for further diagnostics 
    103             DO jk = 1, jpkm1 
    104                ztrdt(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / r2dt(jk) ) - ztrdt(:,:,jk) 
    105                ztrds(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / r2dt(jk) ) - ztrds(:,:,jk) 
    106             END DO 
    107             CALL trd_mod( ztrdt, ztrds, jptra_trd_zdf, 'TRA', kt ) 
    108          ENDIF 
    10997 
    11098      END SELECT 
    11199 
    112       !                                                ! print mean trends (used for debugging) 
     100      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics 
     101         DO jk = 1, jpkm1 
     102            ztrdt(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / r2dt(jk) ) - ztrdt(:,:,jk) 
     103            ztrds(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / r2dt(jk) ) - ztrds(:,:,jk) 
     104         END DO 
     105         CALL trd_mod( ztrdt, ztrds, jptra_trd_zdf, 'TRA', kt ) 
     106      ENDIF 
     107 
     108      !                                          ! print mean trends (used for debugging) 
    113109      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' zdf  - Ta: ', mask1=tmask,               & 
    114110         &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
  • trunk/NEMO/OPA_SRC/TRA/trazdf_exp.F90

    r789 r1110  
    33   !!                    ***  MODULE  trazdf_exp  *** 
    44   !! Ocean active tracers:  vertical component of the tracer mixing trend using 
    5    !!                        an explicit time-stepping (time spllitting scheme) 
     5   !!                        a split-explicit time-stepping  
    66   !!============================================================================== 
    7    !! History : 
    8    !!   6.0  !  90-10  (B. Blanke)  Original code 
    9    !!   7.0  !  91-11  (G. Madec) 
    10    !!        !  92-06  (M. Imbard)  correction on tracer trend loops 
    11    !!        !  96-01  (G. Madec)  statement function for e3 
    12    !!        !  97-05  (G. Madec)  vertical component of isopycnal 
    13    !!        !  97-07  (G. Madec)  geopotential diffusion in s-coord 
    14    !!        !  00-08  (G. Madec)  double diffusive mixing 
    15    !!   8.5  !  02-08  (G. Madec)  F90: Free form and module 
    16    !!   9.0  !  04-08  (C. Talandier) New trends organisation 
    17    !!        !  05-11  (G. Madec)  New organisation 
     7   !! History :  OPA  !  1990-10  (B. Blanke)  Original code 
     8   !!            7.0  !  1991-11  (G. Madec) 
     9   !!                 !  1992-06  (M. Imbard)  correction on tracer trend loops 
     10   !!                 !  1996-01  (G. Madec)  statement function for e3 
     11   !!                 !  1997-05  (G. Madec)  vertical component of isopycnal 
     12   !!                 !  1997-07  (G. Madec)  geopotential diffusion in s-coord 
     13   !!                 !  2000-08  (G. Madec)  double diffusive mixing 
     14   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module 
     15   !!             -   !  2004-08  (C. Talandier) New trends organisation 
     16   !!             -   !  2005-11  (G. Madec)  New organisation 
     17   !!            3.0  !  2008-04  (G. Madec)  leap-frog time stepping done in trazdf 
    1818   !!---------------------------------------------------------------------- 
    19    !!   tra_zdf_exp  : update the tracer trend with the vertical diffusion 
    20    !!                  using an explicit time stepping 
     19 
    2120   !!---------------------------------------------------------------------- 
    22    !! * Modules used 
     21   !!   tra_zdf_exp  : compute the tracer the vertical diffusion trend using a 
     22   !!                  split-explicit time stepping and provide the after tracer 
     23   !!---------------------------------------------------------------------- 
    2324   USE oce             ! ocean dynamics and active tracers  
    2425   USE dom_oce         ! ocean space and time domain  
     26   USE domvvl          ! variablevolume levels 
    2527   USE trdmod          ! ocean active tracers trends  
    2628   USE trdmod_oce      ! ocean variables trends 
     
    3335   PRIVATE 
    3436 
    35    !! * Routine accessibility 
    36    PUBLIC tra_zdf_exp          ! routine called by step.F90 
     37   PUBLIC   tra_zdf_exp   ! routine called by step.F90 
    3738 
    3839   !! * Substitutions 
    3940#  include "domzgr_substitute.h90" 
    4041#  include "zdfddm_substitute.h90" 
     42#  include "vectopt_loop_substitute.h90" 
    4143   !!---------------------------------------------------------------------- 
    42    !!  OPA 9.0 , LOCEAN-IPSL (2005)  
    43    !! $Header$  
    44    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
     44   !! NEMO/OPA  3.0 , LOCEAN-IPSL (2008)  
     45   !! $Id:$  
     46   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4547   !!---------------------------------------------------------------------- 
    4648 
     
    5153      !!                  ***  ROUTINE tra_zdf_exp  *** 
    5254      !!                    
    53       !! ** Purpose :   Compute the trend due to the vertical tracer mixing  
    54       !!      using an explicit time stepping and add it to the general trend  
    55       !!      of the tracer equations. 
     55      !! ** Purpose :   Compute the after tracer fields due to the vertical 
     56      !!      tracer mixing alone, and then due to the whole tracer trend. 
    5657      !! 
    57       !! ** Method  :   The vertical diffusion of tracers (t & s) is given by: 
    58       !!         difft = dz( avt dz(tb) ) = 1/e3t dk+1( avt/e3w dk(tb) ) 
    59       !!      It is evaluated with an Euler scheme, using a time splitting 
    60       !!      technique. 
    61       !!      Surface and bottom boundary conditions: no diffusive flux on 
    62       !!      both tracers (bottom, applied through the masked field avt). 
    63       !!      Add this trend to the general trend ta,sa : 
    64       !!          ta = ta + dz( avt dz(t) ) 
    65       !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm= T) 
     58      !! ** Method  : - The after tracer fields due to the vertical diffusion 
     59      !!      of tracers alone is given by: 
     60      !!                zwx = tb + p2dt difft 
     61      !!      where difft = dz( avt dz(tb) ) = 1/e3t dk+1( avt/e3w dk(tb) ) 
     62      !!           (if lk_zdfddm=T use avs on salinity instead of avt) 
     63      !!      difft is evaluated with an Euler split-explit scheme using a 
     64      !!      no flux boundary condition at both surface and bottomi boundaries. 
     65      !!      (N.B. bottom condition is applied through the masked field avt). 
     66      !!              - the after tracer fields due to the whole trend is  
     67      !!      obtained in leap-frog environment by : 
     68      !!          ta = zwx + p2dt ta 
     69      !!              - in case of variable level thickness (lk_vvl=T) the  
     70      !!     the leap-frog is applied on thickness weighted tracer. That is: 
     71      !!          ta = [ tb*e3tb + e3tn*( zwx - tb + p2dt ta ) ] / e3tn 
    6672      !! 
    67       !! ** Action : - Update (ta,sa) with the before vertical diffusion trend 
     73      !! ** Action : - after tracer fields (ta,sa)  
     74      !!--------------------------------------------------------------------- 
     75      INTEGER , INTENT(in)                 ::   kt     ! ocean time-step index 
     76      REAL(wp), INTENT(in), DIMENSION(jpk) ::   p2dt   ! vertical profile of tracer time-step 
    6877      !! 
    69       !!--------------------------------------------------------------------- 
    70       !! * Arguments 
    71       INTEGER, INTENT( in ) ::   kt           ! ocean time-step index 
    72       REAL(wp), DIMENSION(jpk), INTENT( in ) ::   & 
    73          p2dt                                 ! vertical profile of tracer time-step 
    74        
    75       !! * Local declarations 
    76       INTEGER ::   ji, jj, jk, jl             ! dummy loop indices 
    77       REAL(wp) ::   & 
    78          zlavmr,                            & ! temporary scalars 
    79          zave3r, ze3tr,                     & !    "         " 
    80          zta, zsa                             !    "         "  
    81       REAL(wp), DIMENSION(jpi,jpk) ::   & 
    82          zwx, zwy, zwz, zww 
     78      INTEGER  ::   ji, jj, jk, jl            ! dummy loop indices 
     79      REAL(wp) ::   zlavmr, zave3r, ze3tr     ! temporary scalars 
     80      REAL(wp) ::   zta, zsa, ze3tb, zcoef    ! temporary scalars 
     81      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zwz, zww   ! 3D workspace 
    8382      !!--------------------------------------------------------------------- 
    8483 
     
    8988      ENDIF 
    9089 
     90      ! Initializations 
     91      ! --------------- 
     92      zlavmr = 1. / float( n_zdfexp )                           ! Local constant 
     93      ! 
     94      zwy(:,:, 1 ) = 0.e0        ;   zww(:,:, 1 ) = 0.e0        ! surface boundary conditions: no flux 
     95      zwy(:,:,jpk) = 0.e0        ;   zww(:,:,jpk) = 0.e0        ! bottom  boundary conditions: no flux 
     96      ! 
     97      zwx(:,:,:)   = tb(:,:,:)   ;   zwz(:,:,:)   = sb(:,:,:)   ! zwx and zwz arrays set to before tracer values 
    9198 
    92       ! 0. Local constant initialization 
    93       ! -------------------------------- 
    94       zlavmr = 1. / float( n_zdfexp ) 
    95  
    96       !                                                ! =============== 
    97       DO jj = 2, jpjm1                                 !  Vertical slab 
    98          !                                             ! =============== 
    99          ! 1. Initializations 
    100          ! ------------------ 
    101  
    102          ! Surface & bottom boundary conditions: no flux 
    103          DO ji = 2, jpim1 
    104             zwy(ji, 1 ) = 0.e0 
    105             zwy(ji,jpk) = 0.e0 
    106             zww(ji, 1 ) = 0.e0 
    107             zww(ji,jpk) = 0.e0 
    108          END DO 
    109  
    110          ! zwx and zwz arrays set to before tracer values 
    111          DO jk = 1, jpk 
    112             DO ji = 2, jpim1 
    113                zwx(ji,jk) = tb(ji,jj,jk) 
    114                zwz(ji,jk) = sb(ji,jj,jk) 
    115             END DO 
    116          END DO 
    117  
    118  
    119          ! 2. Time splitting loop 
    120          ! ---------------------- 
    121  
    122          DO jl = 1, n_zdfexp 
    123  
    124             ! first vertical derivative 
    125             IF( lk_zdfddm ) THEN       ! double diffusion: avs /= avt 
    126                DO jk = 2, jpk 
    127                   DO ji = 2, jpim1 
    128                      zave3r = 1.e0 / fse3w(ji,jj,jk)  
    129                      zwy(ji,jk) =   avt(ji,jj,jk) * ( zwx(ji,jk-1) - zwx(ji,jk) ) * zave3r 
    130                      zww(ji,jk) = fsavs(ji,jj,jk) * ( zwz(ji,jk-1) - zwz(ji,jk) ) * zave3r 
    131                   END DO 
    132                END DO 
    133             ELSE                      ! default : avs = avt 
    134                DO jk = 2, jpk 
    135                   DO ji = 2, jpim1 
    136                      zave3r = avt(ji,jj,jk) / fse3w(ji,jj,jk) 
    137                      zwy(ji,jk) = zave3r *(zwx(ji,jk-1) - zwx(ji,jk) ) 
    138                      zww(ji,jk) = zave3r *(zwz(ji,jk-1) - zwz(ji,jk) ) 
    139                   END DO 
    140                END DO 
    141             ENDIF 
    142  
    143             ! trend estimation at kt+l*2*rdt/n_zdfexp 
    144             DO jk = 1, jpkm1 
    145                DO ji = 2, jpim1 
    146                   ze3tr = zlavmr / fse3t(ji,jj,jk) 
    147                   ! 2nd vertical derivative 
    148                   zta = ( zwy(ji,jk) - zwy(ji,jk+1) ) * ze3tr 
    149                   zsa = ( zww(ji,jk) - zww(ji,jk+1) ) * ze3tr 
    150                   ! update the tracer trends 
    151                   ta(ji,jj,jk) = ta(ji,jj,jk) + zta 
    152                   sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 
    153                   ! update tracer fields at kt+l*2*rdt/n_zdfexp 
    154                   zwx(ji,jk) = zwx(ji,jk) + p2dt(jk) * zta * tmask(ji,jj,jk) 
    155                   zwz(ji,jk) = zwz(ji,jk) + p2dt(jk) * zsa * tmask(ji,jj,jk) 
     99      ! Split-explicit loop  (after tracer due to the vertical diffusion alone) 
     100      ! ------------------- 
     101      ! 
     102      DO jl = 1, n_zdfexp 
     103         !                     ! first vertical derivative 
     104         DO jk = 2, jpk 
     105            DO jj = 2, jpjm1  
     106               DO ji = fs_2, fs_jpim1   ! vector opt. 
     107                  zave3r = 1.e0 / fse3w(ji,jj,jk)  
     108                  zwy(ji,jj,jk) =   avt(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) * zave3r 
     109                  zww(ji,jj,jk) = fsavs(ji,jj,jk) * ( zwz(ji,jj,jk-1) - zwz(ji,jj,jk) ) * zave3r 
    156110               END DO 
    157111            END DO 
    158112         END DO 
    159          !                                             ! =============== 
    160       END DO                                           !   End of slab 
    161       !                                                ! =============== 
     113         ! 
     114         DO jk = 1, jpkm1      ! second vertical derivative   ==> tracer at kt+l*2*rdt/n_zdfexp 
     115            DO jj = 2, jpjm1  
     116               DO ji = fs_2, fs_jpim1   ! vector opt. 
     117                  ze3tr = zlavmr / fse3t(ji,jj,jk) 
     118                  zwx(ji,jj,jk) = zwx(ji,jj,jk) + p2dt(jk) * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) * ze3tr 
     119                  zwz(ji,jj,jk) = zwz(ji,jj,jk) + p2dt(jk) * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) ) * ze3tr 
     120               END DO 
     121            END DO 
     122         END DO 
     123         ! 
     124      END DO 
    162125 
     126      ! After tracer due to all trends 
     127      ! ------------------------------ 
     128      IF( lk_vvl ) THEN          ! variable level thickness : leap-frog on tracer*e3t 
     129         DO jk = 1, jpkm1 
     130            DO jj = 2, jpjm1  
     131               DO ji = fs_2, fs_jpim1   ! vector opt. 
     132                  ze3tb = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )       ! before e3t 
     133                  zta   = zwx(ji,jj,jk) - tb(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk)       ! total trends * 2*rdt 
     134                  zsa   = zwz(ji,jj,jk) - sb(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk)      
     135                  zcoef = 1.e0 / fse3t(ji,jj,jk) * tmask(ji,jj,jk) 
     136                  ta(ji,jj,jk) = (  ze3tb * tb(ji,jj,jk) + fse3t(ji,jj,jk) * zta  ) * zcoef 
     137                  sa(ji,jj,jk) = (  ze3tb * sb(ji,jj,jk) + fse3t(ji,jj,jk) * zsa  ) * zcoef 
     138               END DO 
     139            END DO 
     140         END DO 
     141      ELSE                       ! fixed level thickness : leap-frog on tracers 
     142         DO jk = 1, jpkm1 
     143            DO jj = 2, jpjm1  
     144               DO ji = fs_2, fs_jpim1   ! vector opt. 
     145                  ta(ji,jj,jk) = ( zwx(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ) *tmask(ji,jj,jk) 
     146                  sa(ji,jj,jk) = ( zwz(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) ) *tmask(ji,jj,jk) 
     147               END DO 
     148            END DO 
     149         END DO 
     150      ENDIF 
     151      ! 
    163152   END SUBROUTINE tra_zdf_exp 
    164153 
Note: See TracChangeset for help on using the changeset viewer.