Changeset 1847


Ignore:
Timestamp:
2010-04-21T15:38:47+02:00 (11 years ago)
Author:
gm
Message:

MFL: change the time averaged in the semi-implicit hpg (tranxt) see ticket #663

File:
1 edited

Legend:

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

    r1601 r1847  
    1515   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf 
    1616   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option 
     17   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter 
    1718   !!---------------------------------------------------------------------- 
    1819 
     
    4445   PUBLIC   tra_nxt    ! routine called by step.F90 
    4546 
    46    REAL(wp), DIMENSION(jpk) ::   r2dt_t   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 
     47   REAL(wp)                 ::   rbcp, r1_2bcp   ! Brown & Campana parameters for semi-implicit hpg 
     48   REAL(wp), DIMENSION(jpk) ::   r2dt_t          ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 
    4749 
    4850   !! * Substitutions 
    4951#  include "domzgr_substitute.h90" 
    5052   !!---------------------------------------------------------------------- 
    51    !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     53   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    5254   !! $Id$ 
    5355   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    8587      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    8688      !! 
    87       INTEGER  ::   jk    ! dummy loop indices 
    88       REAL(wp) ::   zfact ! temporary scalars 
     89      INTEGER  ::   jk      ! dummy loop indices 
     90      REAL(wp) ::   zfact   ! temporary scalars 
    8991      !!---------------------------------------------------------------------- 
    9092 
     
    9395         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 
    9496         IF(lwp) WRITE(numout,*) '~~~~~~~' 
     97         ! 
     98         rbcp    = 0.25 * (1 + atfp) * (1 + atfp * atfp)       ! Brown & Campana parameter for semi-implicit hpg 
     99         r1_2bcp = 1.e0 - 2.e0 * rbcp 
    95100      ENDIF 
    96101 
     
    160165      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T) 
    161166      !!              - swap tracer fields to prepare the next time_step. 
    162       !!                This can be summurized for tempearture as: 
    163       !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T 
    164       !!             ztm = 0                   otherwise 
     167      !!                This can be summurized for temperature as: 
     168      !!             ztm = (rbcp*ta+(2-rbcp)*tn+rbcp*tb)       ln_dynhpg_imp = T 
     169      !!             ztm = 0                                   otherwise 
     170      !!                   with rbcp=(1+atfp)*(1+atfp*atfp)/4 
    165171      !!             tb  = tn + atfp*[ tb - 2 tn + ta ] 
    166172      !!             tn  = ta  
     
    189195         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
    190196            DO jk = 1, jpkm1                              ! (only swap) 
     197               tn(:,:,jk) = ta(:,:,jk)                                                          ! tb <-- tn 
     198               sn(:,:,jk) = sa(:,:,jk) 
     199            END DO 
     200         ELSE                                             ! general case (Leapfrog + Asselin filter) 
     201            DO jk = 1, jpkm1 
    191202               DO jj = 1, jpj 
    192203                  DO ji = 1, jpi 
    193                      tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tb <-- tn 
     204                     ztm = rbcp * ( ta(ji,jj,jk) + tb(ji,jj,jk) ) + r1_2bcp * tn(ji,jj,jk)      ! mean t 
     205                     zsm = rbcp * ( sa(ji,jj,jk) + sb(ji,jj,jk) ) + r1_2bcp * sn(ji,jj,jk)  
     206                     ztf = atfp * ( ta(ji,jj,jk) + tb(ji,jj,jk)   - 2.      * tn(ji,jj,jk) )    ! Asselin filter on t  
     207                     zsf = atfp * ( sa(ji,jj,jk) + sb(ji,jj,jk)   - 2.      * sn(ji,jj,jk) ) 
     208                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                          ! tb <-- filtered tn  
     209                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 
     210                     tn(ji,jj,jk) = ta(ji,jj,jk)                                                ! tn <-- ta 
    194211                     sn(ji,jj,jk) = sa(ji,jj,jk) 
    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 
     212                     ta(ji,jj,jk) = ztm                                                         ! ta <-- mean t 
    211213                     sa(ji,jj,jk) = zsm 
    212214                  END DO 
     
    220222         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
    221223            DO jk = 1, jpkm1 
     224               tn(:,:,jk) = ta(:,:,jk)                                                    ! tn <-- ta 
     225               sn(:,:,jk) = sa(:,:,jk) 
     226            END DO 
     227         ELSE                                             ! general case (Leapfrog + Asselin filter 
     228            DO jk = 1, jpkm1 
    222229               DO jj = 1, jpj 
    223230                  DO ji = 1, jpi 
    224                      tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta 
    225                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    226                   END DO 
    227                END DO 
    228             END DO 
    229          ELSE                                             ! general case (Leapfrog + Asselin filter 
    230             DO jk = 1, jpkm1 
    231                DO jj = 1, jpj 
    232                   DO ji = 1, jpi 
    233                      ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t  
     231                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )      ! Asselin filter on t  
    234232                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 
    235                      tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn  
     233                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                    ! tb <-- filtered tn  
    236234                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 
    237                      tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta 
     235                     tn(ji,jj,jk) = ta(ji,jj,jk)                                          ! tn <-- ta 
    238236                     sn(ji,jj,jk) = sa(ji,jj,jk) 
    239237                  END DO 
     
    258256      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 
    259257      !!              - swap tracer fields to prepare the next time_step. 
    260       !!                This can be summurized for tempearture as: 
    261       !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T 
    262       !!                  /(e3t_a   +2*e3t_n   +e3t_b   )     
    263       !!             ztm = 0                                otherwise 
     258      !!                This can be summurized for temperature as: 
     259      !!             ztm = ( rbcp*e3t_a*ta + (2-rbtp)*e3t_n*tn + rbcp*e3t_b*tb )   ln_dynhpg_imp = T 
     260      !!                  /( rbcp*e3t_a    + (2-rbcp)*e3t_n    + rbcp*e3t_b    )    
     261      !!             ztm = 0                                                       otherwise 
    264262      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 
    265263      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] ) 
     
    286284 
    287285      !                                              ! ----------------------- ! 
    288       IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case ! 
     286      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !  (mean tracer computation) 
    289287         !                                           ! ----------------------- ! 
    290288         ! 
    291289         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
    292             DO jk = 1, jpkm1                              ! (only swap) 
    293                DO jj = 1, jpj 
    294                   DO ji = 1, jpi 
    295                      tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta 
    296                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    297                   END DO 
    298                END DO 
    299             END DO 
    300          ELSE 
     290            DO jk = 1, jpkm1                                   ! (only swap) 
     291               tn(:,:,jk) = ta(:,:,jk)                         ! tn <-- ta 
     292               sn(:,:,jk) = sa(:,:,jk) 
     293            END DO 
     294            !                                             ! general case (Leapfrog + Asselin filter) 
     295         ELSE                                             ! apply filter on thickness weighted tracer and swap 
    301296            DO jk = 1, jpkm1 
    302297               DO jj = 1, jpj 
     
    311306                     ! 
    312307                     !                                         ! Asselin filter on thickness and tracer content 
    313                      ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 
    314                      ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   )  
    315                      zsc_f  = atfp * ( zsca   - 2.* zscn   + zscb   )  
     308                     ze3t_f = atfp * ( ze3t_a + ze3t_b - 2.* ze3t_n ) 
     309                     ztc_f  = atfp * ( ztca   + ztcb   - 2.* ztcn   )  
     310                     zsc_f  = atfp * ( zsca   + zscb   - 2.* zscn   )  
    316311                     ! 
    317312                     !                                         ! filtered tracer including the correction  
     
    320315                     zsf   = ze3fr * ( zscn   + zsc_f  ) 
    321316                     !                                         ! mean thickness and tracer 
    322                      ze3mr = 1.e0  / ( ze3t_a + 2.* ze3t_n + ze3t_b ) 
    323                      ztm   = ze3mr * ( ztca   + 2.* ztcn   + ztcb   ) 
    324                      zsm   = ze3mr * ( zsca   + 2.* zscn   + zscb   ) 
     317                     ze3mr = 1.e0  / (  rbcp * ( ze3t_a + ze3t_b ) + r1_2bcp * ze3t_n ) 
     318                     ztm   = ze3mr * (  rbcp * ( ztca   + ztcb   ) + r1_2bcp * ztcn   ) 
     319                     zsm   = ze3mr * (  rbcp * ( zsca   + zscb   ) + r1_2bcp * zscn   ) 
    325320!!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !! 
    326 !!gm                 e3t_m(ji,jj,jk) = 0.25 / ze3mr 
     321!!gm                 e3t_m(ji,jj,jk) = 1.e0 / ze3mr 
    327322                     !                                         ! swap of arrays 
    328323                     tb(ji,jj,jk) = ztf                             ! tb <-- tn + filter 
     
    337332         ENDIF 
    338333         !                                           ! ----------------------- ! 
    339       ELSE                                           !    explicit hpg case    ! 
     334      ELSE                                           !    explicit hpg case    !  (NO mean tracer computation) 
    340335         !                                           ! ----------------------- ! 
    341336         ! 
    342337         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step 
    343             DO jk = 1, jpkm1                              ! No filter nor thickness weighting computation required 
    344                DO jj = 1, jpj                             ! ONLY swap 
    345                   DO ji = 1, jpi 
    346                      tn(ji,jj,jk) = ta(ji,jj,jk)                                 ! tn <-- ta 
    347                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    348                   END DO 
    349                END DO 
     338            DO jk = 1, jpkm1                                   ! (only swap) 
     339               tn(:,:,jk) = ta(:,:,jk)                         ! tn <-- ta 
     340               sn(:,:,jk) = sa(:,:,jk) 
    350341            END DO 
    351342            !                                             ! general case (Leapfrog + Asselin filter) 
Note: See TracChangeset for help on using the changeset viewer.