Changeset 2602


Ignore:
Timestamp:
2011-02-21T16:23:38+01:00 (10 years ago)
Author:
cetlod
Message:

Improvment of surface boundary condition on trazdf_imp.F90 : minor bug correction+style, see ticket #800

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r2528 r2602  
    1515   !!            3.2  !  2009-03  (G. Madec)  heat and salt content trends 
    1616   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC 
     17   !!             -   !  2011-02  (A. Coward, C. Ethe, G. Madec) improvment of surface boundary condition 
    1718   !!---------------------------------------------------------------------- 
    1819   
     
    3940   PUBLIC   tra_zdf_imp   !  routine called by step.F90 
    4041 
     42   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise  
     43 
    4144   !! * Substitutions 
    4245#  include "domzgr_substitute.h90" 
     
    5558      !!                  ***  ROUTINE tra_zdf_imp  *** 
    5659      !! 
    57       !! ** Purpose :   Compute the trend due to the vertical tracer diffusion 
    58       !!     including the vertical component of lateral mixing (only for 2nd 
    59       !!     order operator, for fourth order it is already computed and add 
    60       !!     to the general trend in traldf.F) and add it to the general trend 
    61       !!     of the tracer equations. 
    62       !! 
    63       !! ** Method  :   The vertical component of the lateral diffusive trends 
    64       !!      is provided by a 2nd order operator rotated along neutral or geo- 
    65       !!      potential surfaces to which an eddy induced advection can be  
    66       !!      added. It is computed using before fields (forward in time) and  
    67       !!      isopycnal or geopotential slopes computed in routine ldfslp. 
    68       !! 
    69       !!      Second part: vertical trend associated with the vertical physics 
    70       !!      ===========  (including the vertical flux proportional to dk[t] 
    71       !!                  associated with the lateral mixing, through the 
    72       !!                  update of avt) 
    73       !!      The vertical diffusion of the tracer t  is given by: 
    74       !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 
     60      !! ** Purpose :   Compute the after tracer through a implicit computation 
     61      !!     of the vertical tracer diffusion (including the vertical component  
     62      !!     of lateral mixing (only for 2nd order operator, for fourth order  
     63      !!     it is already computed and add to the general trend in traldf)  
     64      !! 
     65      !! ** Method  :  The vertical diffusion of the tracer t  is given by: 
     66      !!                  difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 
    7567      !!      It is computed using a backward time scheme (t=ta). 
     68      !!      If lk_zdfddm=T, use avs for salinity or for passive tracers 
    7669      !!      Surface and bottom boundary conditions: no diffusive flux on 
    7770      !!      both tracers (bottom, applied through the masked field avt). 
    78       !!      Add this trend to the general trend ta,sa : 
    79       !!         ta = ta + dz( avt dz(t) ) 
    80       !!         if lk_zdfddm=T, use avs for salinity or for passive tracers 
    81       !!         (sa = sa + dz( avs dz(t) )  
    82       !! 
    83       !!      Third part: recover avt resulting from the vertical physics 
    84       !!      ==========  alone, for further diagnostics (for example to 
    85       !!                  compute the turbocline depth in zdfmxl.F90). 
    86       !!         if lk_zdfddm=T, use avt = zavt 
    87       !!         (avs = zavs if lk_zdfddm=T ) 
    88       !! 
    89       !! ** Action  : - Update (ta) with before vertical diffusion trend 
     71      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing. 
     72      !! 
     73      !! ** Action  : - pta  becomes the after tracer 
    9074      !!--------------------------------------------------------------------- 
    9175      USE oce    , ONLY :   zwd   => ua   ! ua used as workspace 
     
    10084      !! 
    10185      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices 
    102       REAL(wp) ::  zavi, zrhs, znvvl     ! local scalars 
     86      REAL(wp) ::  zrhs                  ! local scalars 
    10387      REAL(wp) ::  ze3tb, ze3tn, ze3ta   ! variable vertical scale factors 
    10488      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi, zwt   ! workspace arrays 
     
    10993         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype 
    11094         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' 
    111          zavi = 0._wp      ! avoid warning at compilation phase when lk_ldfslp=F 
     95         ! 
     96         IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator 
     97         ELSE                ;    r_vvl = 0._wp        
     98         ENDIF 
    11299      ENDIF 
    113100      ! 
    114       ! I. Local initialization 
    115       ! ----------------------- 
    116       zwd(1,:, : ) = 0._wp     ;     zwd(jpi,:,:) = 0._wp 
    117       zws(1,:, : ) = 0._wp     ;     zws(jpi,:,:) = 0._wp 
    118       zwi(1,:, : ) = 0._wp     ;     zwi(jpi,:,:) = 0._wp 
    119       zwt(1,:, : ) = 0._wp     ;     zwt(jpi,:,:) = 0._wp 
    120       zwt(:,:,jpk) = 0._wp     ;     zwt( : ,:,1) = 0._wp 
    121  
    122       ! I.1 Variable volume : to take into account vertical variable vertical scale factors 
    123       ! ------------------- 
    124       IF( lk_vvl ) THEN   ;    znvvl = 1._wp 
    125       ELSE                ;    znvvl = 0._wp 
    126       ENDIF 
    127  
    128       ! II. Vertical trend associated with the vertical physics 
    129       ! ======================================================= 
    130       !     (including the vertical flux proportional to dk[t] associated 
    131       !      with the lateral mixing, through the avt update) 
    132       !     dk[ avt dk[ (t,s) ] ] diffusive trends 
    133  
    134       ! 
    135       ! II.0 Matrix construction 
    136       ! ------------------------ 
    137       DO jn = 1, kjpt 
     101      !                                               ! ============= ! 
     102      DO jn = 1, kjpt                                 !  tracer loop  ! 
     103         !                                            ! ============= ! 
    138104         ! 
    139105         !  Matrix construction 
    140          ! ------------------------ 
    141          IF( cdtype == 'TRA' .AND. jn == jp_tem )  THEN  
     106         ! -------------------- 
     107         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer 
     108         ! 
     109         IF(  ( cdtype == 'TRA' .AND. ( ( jn == jp_tem ) .OR. ( jn == jp_sal .AND. lk_zdfddm ) ) ) .OR. & 
     110            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN 
     111            ! 
     112            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers 
     113            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt  (:,:,2:jpk) 
     114            ELSE                                            ;   zwt(:,:,2:jpk) = fsavs(:,:,2:jpk) 
     115            ENDIF 
     116            zwt(:,:,1) = 0._wp 
     117            ! 
    142118#if defined key_ldfslp 
    143             IF( ln_traldf_grif ) THEN 
     119            ! isoneutral diffusion: add the contribution  
     120            IF( ln_traldf_grif    ) THEN     ! Griffies isoneutral diff 
    144121               DO jk = 2, jpkm1 
    145122                  DO jj = 2, jpjm1 
    146123                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    147                         ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing 
    148                         zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing 
    149                         zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature) 
     124                        zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)        
    150125                     END DO 
    151126                  END DO 
    152127               END DO 
    153             ! update and save of avt (and avs if double diffusive mixing) 
    154             ELSE IF( l_traldf_rot ) THEN 
     128            ELSE IF( l_traldf_rot ) THEN     ! standard isoneutral diff 
    155129               DO jk = 2, jpkm1 
    156130                  DO jj = 2, jpjm1 
    157131                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    158                         zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing 
    159                            & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    160                            &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  ) 
    161                         zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature) 
     132                        zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk)                       & 
     133                           &                          * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     134                           &                             + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  ) 
    162135                     END DO 
    163136                  END DO 
    164137               END DO 
    165             ELSE                         ! no rotation but key_ldfslp defined 
    166                zwt(:,:,:) = avt(:,:,:) 
    167138            ENDIF 
    168 #else 
    169             ! No isopycnal diffusion 
    170             zwt(:,:,:) = avt(:,:,:)            
    171139#endif 
    172             ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked) 
     140            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked) 
    173141            DO jk = 1, jpkm1 
    174142               DO jj = 2, jpjm1 
    175143                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    176                      ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point 
    177                      ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point 
     144                     ze3ta =  ( 1. - r_vvl ) +        r_vvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point 
     145                     ze3tn =         r_vvl   + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point 
    178146                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) ) 
    179147                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 
     
    182150               END DO 
    183151            END DO 
     152            ! 
     153            !! Matrix inversion from the first level 
     154            !!---------------------------------------------------------------------- 
     155            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk ) 
     156            ! 
     157            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 ) 
     158            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 ) 
     159            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 ) 
     160            !        (        ...               )( ...  ) ( ...  ) 
     161            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
     162            ! 
     163            !   m is decomposed in the product of an upper and lower triangular matrix. 
     164            !   The 3 diagonal terms are in 3d arrays: zwd, zws, zwi. 
     165            !   Suffices i,s and d indicate "inferior" (below diagonal), diagonal 
     166            !   and "superior" (above diagonal) components of the tridiagonal system. 
     167            !   The solution will be in the 4d array pta. 
     168            !   The 3d array zwt is used as a work space array. 
     169            !   En route to the solution pta is used a to evaluate the rhs and then  
     170            !   used as a work space array: its value is modified. 
     171            ! 
    184172            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) 
     173            ! done once for all passive tracers (so included in the IF instruction) 
    185174            DO jj = 2, jpjm1 
    186175               DO ji = fs_2, fs_jpim1 
     
    196185            END DO 
    197186            ! 
    198          ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN 
    199 #if defined key_ldfslp 
    200             IF( ln_traldf_grif ) THEN 
    201                DO jk = 2, jpkm1 
    202                   DO jj = 2, jpjm1 
    203                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    204                         zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing 
    205                         ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing 
    206                         zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature) 
    207                      END DO 
    208                   END DO 
    209                END DO 
    210             ELSE IF( l_traldf_rot ) THEN 
    211                DO jk = 2, jpkm1 
    212                   DO jj = 2, jpjm1 
    213                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    214                         zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing 
    215                            & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    216                            &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  ) 
    217                         zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on salinity) 
    218                      END DO 
    219                   END DO 
    220                END DO 
    221             ELSE                         ! no rotation but key_ldfslp defined 
    222                zwt(:,:,:) = fsavs(:,:,:) 
    223             ENDIF 
    224 #else 
    225             ! No isopycnal diffusion 
    226             zwt(:,:,:) = fsavs(:,:,:)            
    227 #endif 
    228             ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked) 
    229             DO jk = 1, jpkm1 
    230                DO jj = 2, jpjm1 
    231                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    232                      ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point 
    233                      ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point 
    234                      zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) ) 
    235                      zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 
    236                      zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    237                  END DO 
    238                END DO 
    239             END DO 
    240             ! Surface boudary conditions 
    241             DO jj = 2, jpjm1 
    242                DO ji = fs_2, fs_jpim1   ! vector opt. 
    243                  ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point 
    244                  zwi(ji,jj,1) = 0._wp 
    245                  zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 
    246                END DO 
    247             END DO 
    248             ! 
    249             ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) 
    250             DO jj = 2, jpjm1 
    251                DO ji = fs_2, fs_jpim1 
    252                   zwt(ji,jj,1) = zwd(ji,jj,1) 
    253                END DO 
    254             END DO 
    255             DO jk = 2, jpkm1 
    256                DO jj = 2, jpjm1 
    257                   DO ji = fs_2, fs_jpim1 
    258                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
    259                   END DO 
    260                END DO 
    261             END DO 
    262             ! 
    263187         END IF  
    264  
    265          ! II.1. Vertical diffusion on tracer 
    266          ! --------------------------- 
    267           
     188         !          
    268189         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    269190         DO jj = 2, jpjm1 
    270191            DO ji = fs_2, fs_jpim1 
    271                ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1) 
    272                ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1) 
     192               ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 
     193               ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 
    273194               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 
    274195            END DO 
     
    277198            DO jj = 2, jpjm1 
    278199               DO ji = fs_2, fs_jpim1 
    279                   ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk) 
    280                   ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk) 
     200                  ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 
     201                  ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t  (ji,jj,jk) 
    281202                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side  
    282203                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 
     
    285206         END DO 
    286207 
    287          ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk 
     208         ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer) 
    288209         DO jj = 2, jpjm1 
    289210            DO ji = fs_2, fs_jpim1 
     
    299220            END DO 
    300221         END DO 
    301          ! 
    302       END DO 
     222         !                                            ! ================= ! 
     223      END DO                                          !  end tracer loop  ! 
     224      !                                               ! ================= ! 
    303225      ! 
    304226   END SUBROUTINE tra_zdf_imp 
Note: See TracChangeset for help on using the changeset viewer.