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 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90 – NEMO

Ignore:
Timestamp:
2010-12-27T18:33:53+01:00 (13 years ago)
Author:
rblod
Message:

Update NEMOGCM from branch nemo_v3_3_beta

File:
1 edited

Legend:

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

    • Property svn:eol-style deleted
    r2468 r2528  
    22   !!====================================================================== 
    33   !!                   ***  MODULE  traldf_iso  *** 
    4    !! Ocean active tracers:  horizontal component of the lateral tracer mixing trend 
     4   !! Ocean tracers:  horizontal component of the lateral tracer mixing trend 
    55   !!====================================================================== 
    6    !! History :        !  94-08  (G. Madec, M. Imbard) 
    7    !!                  !  97-05  (G. Madec)  split into traldf and trazdf 
    8    !!             8.5  !  02-08  (G. Madec)  Free form, F90 
    9    !!             9.0  !  05-11  (G. Madec)  merge traldf and trazdf :-) 
     6   !! History :  OPA  !  1994-08  (G. Madec, M. Imbard) 
     7   !!            8.0  !  1997-05  (G. Madec)  split into traldf and trazdf 
     8   !!            NEMO !  2002-08  (G. Madec)  Free form, F90 
     9   !!            1.0  !  2005-11  (G. Madec)  merge traldf and trazdf :-) 
     10   !!            3.3  !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC 
    1011   !!---------------------------------------------------------------------- 
    1112#if   defined key_ldfslp   ||   defined key_esopa 
    1213   !!---------------------------------------------------------------------- 
    1314   !!   'key_ldfslp'               slope of the lateral diffusive direction 
    14    !!---------------------------------------------------------------------- 
    1515   !!---------------------------------------------------------------------- 
    1616   !!   tra_ldf_iso  : update the tracer trend with the horizontal  
     
    2121   USE oce             ! ocean dynamics and active tracers 
    2222   USE dom_oce         ! ocean space and time domain 
     23   USE trc_oce         ! share passive tracers/Ocean variables 
     24   USE zdf_oce         ! ocean vertical physics 
    2325   USE ldftra_oce      ! ocean active tracers: lateral physics 
    24    USE trdmod          ! ocean active tracers trends  
    25    USE trdmod_oce      ! ocean variables trends 
    26    USE zdf_oce         ! ocean vertical physics 
    27    USE in_out_manager  ! I/O manager 
    28    USE iom             ! 
    2926   USE ldfslp          ! iso-neutral slopes 
    3027   USE diaptr          ! poleward transport diagnostics 
    31    USE prtctl          ! Print control 
     28   USE in_out_manager  ! I/O manager 
     29   USE iom             ! I/O library 
    3230#if defined key_diaar5 
    3331   USE phycst          ! physical constants 
     
    4543#  include "vectopt_loop_substitute.h90" 
    4644   !!---------------------------------------------------------------------- 
    47    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
    48    !! $Id$  
    49    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    50    !!---------------------------------------------------------------------- 
    51  
     45   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     46   !! $Id$ 
     47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     48   !!---------------------------------------------------------------------- 
    5249CONTAINS 
    5350 
    54    SUBROUTINE tra_ldf_iso( kt ) 
     51   SUBROUTINE tra_ldf_iso( kt, cdtype, pgu, pgv,              & 
     52      &                                ptb, pta, kjpt, pahtb0 ) 
    5553      !!---------------------------------------------------------------------- 
    5654      !!                  ***  ROUTINE tra_ldf_iso  *** 
     
    6664      !!      nal or geopotential slopes computed in routine ldfslp. 
    6765      !! 
    68       !!      1st part :  masked horizontal derivative of T & S ( di[ t ] ) 
     66      !!      1st part :  masked horizontal derivative of T ( di[ t ] ) 
    6967      !!      ========    with partial cell update if ln_zps=T. 
    7068      !! 
     
    8886      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ] 
    8987      !!      Add this trend to the general trend (ta,sa): 
    90       !!         ta = ta + difft 
    91       !! 
    92       !! ** Action :   Update (ta,sa) arrays with the before rotated diffusion 
    93       !!            trend (except the dk[ dk[.] ] term) 
     88      !!         pta = pta + difft 
     89      !! 
     90      !! ** Action :   Update pta arrays with the before rotated diffusion 
    9491      !!---------------------------------------------------------------------- 
    95       USE oce           , zftv => ua   ! use ua as workspace 
    96       USE oce           , zfsv => va   ! use va as workspace 
    97       !! 
    98       INTEGER, INTENT( in ) ::   kt    ! ocean time-step index 
    99       !! 
    100       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    101       INTEGER  ::   iku, ikv     ! temporary integer 
    102       REAL(wp) ::   zmsku, zabe1, zcof1, zcoef3, zta   ! temporary scalars 
    103       REAL(wp) ::   zmskv, zabe2, zcof2, zcoef4, zsa   !    "         " 
    104       REAL(wp) ::   zcoef0, zbtr                       !    "         " 
    105       REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt , zdk1t         ! 2D workspace 
    106       REAL(wp), DIMENSION(jpi,jpj)     ::   zdks , zdk1s, zfsu   !  "         " 
     92      USE oce         , zftu => ua   ! use ua as workspace 
     93      USE oce         , zftv => va   ! use va as workspace 
     94      !! 
     95      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
     96      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     97      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
     98      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
     99      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
     100      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
     101      REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef 
     102      !! 
     103      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     104      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
     105      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
     106      REAL(wp) ::  zcoef0, zbtr, ztra            !   -      - 
     107      REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt, zdk1t         ! 2D workspace 
     108      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, ztfw    ! 3D workspace 
    107109#if defined key_diaar5 
    108       REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                  !  "         " 
    109       REAL(wp)                         ::   zztmp                !  "         " 
    110       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zftu                 ! 3D workspace 
    111 #else 
    112       REAL(wp), DIMENSION(jpi,jpj)     ::   zftu                 ! 2D workspace 
     110      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                 ! 2D workspace 
     111      REAL(wp)                         ::   zztmp               ! local scalar 
    113112#endif 
    114       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, ztfw     ! 3D workspace 
    115       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdis, zdjs, zsfw     !  "      " 
    116113      !!---------------------------------------------------------------------- 
    117114 
    118       IF( kt == nit000 ) THEN 
     115      IF( kt == nit000 )  THEN 
    119116         IF(lwp) WRITE(numout,*) 
    120          IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator' 
     117         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype 
    121118         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    122119      ENDIF 
    123  
    124       !!---------------------------------------------------------------------- 
    125       !!   I - masked horizontal derivative of T & S 
    126       !!---------------------------------------------------------------------- 
    127 !!bug ajout.... why?   ( 1,jpj,:) and (jpi,1,:) should be sufficient.... 
    128       zdit (1,:,:) = 0.e0     ;     zdit (jpi,:,:) = 0.e0 
    129       zdis (1,:,:) = 0.e0     ;     zdis (jpi,:,:) = 0.e0 
    130       zdjt (1,:,:) = 0.e0     ;     zdjt (jpi,:,:) = 0.e0 
    131       zdjs (1,:,:) = 0.e0     ;     zdjs (jpi,:,:) = 0.e0 
    132 !!end 
    133  
    134       ! Horizontal temperature and salinity gradient  
    135       DO jk = 1, jpkm1 
    136          DO jj = 1, jpjm1 
    137             DO ji = 1, fs_jpim1   ! vector opt. 
    138                zdit(ji,jj,jk) = ( tb(ji+1,jj  ,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk) 
    139                zdis(ji,jj,jk) = ( sb(ji+1,jj  ,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk) 
    140                zdjt(ji,jj,jk) = ( tb(ji  ,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk) 
    141                zdjs(ji,jj,jk) = ( sb(ji  ,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk) 
     120      ! 
     121      !                                                          ! =========== 
     122      DO jn = 1, kjpt                                            ! tracer loop 
     123         !                                                       ! =========== 
     124         !                                                
     125         !!---------------------------------------------------------------------- 
     126         !!   I - masked horizontal derivative  
     127         !!---------------------------------------------------------------------- 
     128         !!bug ajout.... why?   ( 1,jpj,:) and (jpi,1,:) should be sufficient.... 
     129         zdit (1,:,:) = 0.e0     ;     zdit (jpi,:,:) = 0.e0 
     130         zdjt (1,:,:) = 0.e0     ;     zdjt (jpi,:,:) = 0.e0 
     131         !!end 
     132 
     133         ! Horizontal tracer gradient  
     134         DO jk = 1, jpkm1 
     135            DO jj = 1, jpjm1 
     136               DO ji = 1, fs_jpim1   ! vector opt. 
     137                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     138                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     139               END DO 
    142140            END DO 
    143141         END DO 
    144       END DO 
    145       IF( ln_zps ) THEN      ! partial steps correction at the last level  
    146          DO jj = 1, jpjm1 
    147             DO ji = 1, fs_jpim1   ! vector opt. 
    148                ! last level 
    149                iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1 
    150                ikv = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1 
    151                zdit(ji,jj,iku) = gtu(ji,jj)  
    152                zdis(ji,jj,iku) = gsu(ji,jj)                
    153                zdjt(ji,jj,ikv) = gtv(ji,jj)  
    154                zdjs(ji,jj,ikv) = gsv(ji,jj)                
     142         IF( ln_zps ) THEN      ! partial steps correction at the last ocean level  
     143            DO jj = 1, jpjm1 
     144               DO ji = 1, fs_jpim1   ! vector opt. 
     145                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
     146                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)       
     147               END DO 
     148            END DO 
     149         ENDIF 
     150 
     151         !!---------------------------------------------------------------------- 
     152         !!   II - horizontal trend  (full) 
     153         !!---------------------------------------------------------------------- 
     154!CDIR PARALLEL DO PRIVATE( zdk1t )  
     155         !                                                ! =============== 
     156         DO jk = 1, jpkm1                                 ! Horizontal slab 
     157            !                                             ! =============== 
     158            ! 1. Vertical tracer gradient at level jk and jk+1 
     159            ! ------------------------------------------------ 
     160            ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
     161            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 
     162            ! 
     163            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:) 
     164            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 
     165            ENDIF 
     166 
     167            ! 2. Horizontal fluxes 
     168            ! --------------------    
     169            DO jj = 1 , jpjm1 
     170               DO ji = 1, fs_jpim1   ! vector opt. 
     171                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) 
     172                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) 
     173                  ! 
     174                  zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   & 
     175                     &             + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. ) 
     176                  ! 
     177                  zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   & 
     178                     &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. ) 
     179                  ! 
     180                  zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
     181                  zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
     182                  ! 
     183                  zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
     184                     &              + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
     185                     &                         + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
     186                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
     187                     &              + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
     188                     &                         + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                   
     189               END DO 
     190            END DO 
     191 
     192            ! II.4 Second derivative (divergence) and add to the general trend 
     193            ! ---------------------------------------------------------------- 
     194            DO jj = 2 , jpjm1 
     195               DO ji = fs_2, fs_jpim1   ! vector opt. 
     196                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     197                  ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  ) 
     198                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     199               END DO 
     200            END DO 
     201            !                                          ! =============== 
     202         END DO                                        !   End of slab   
     203         !                                             ! =============== 
     204         ! 
     205         ! "Poleward" diffusive heat or salt transports (T-S case only) 
     206         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
     207            IF( jn == jp_tem)   htr_ldf(:) = ptr_vj( zftv(:,:,:) ) 
     208            IF( jn == jp_sal)   str_ldf(:) = ptr_vj( zftv(:,:,:) ) 
     209         ENDIF 
     210  
     211#if defined key_diaar5 
     212         IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
     213            z2d(:,:) = 0._wp  
     214            zztmp = rau0 * rcp  
     215            DO jk = 1, jpkm1 
     216               DO jj = 2, jpjm1 
     217                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     218                     z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
     219                  END DO 
     220               END DO 
     221            END DO 
     222            z2d(:,:) = zztmp * z2d(:,:) 
     223            CALL lbc_lnk( z2d, 'U', -1. ) 
     224            CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
     225            z2d(:,:) = 0._wp  
     226            DO jk = 1, jpkm1 
     227               DO jj = 2, jpjm1 
     228                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     229                     z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
     230                  END DO 
     231               END DO 
     232            END DO 
     233            z2d(:,:) = zztmp * z2d(:,:) 
     234            CALL lbc_lnk( z2d, 'V', -1. ) 
     235            CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
     236         END IF 
     237#endif 
     238 
     239         !!---------------------------------------------------------------------- 
     240         !!   III - vertical trend of T & S (extra diagonal terms only) 
     241         !!---------------------------------------------------------------------- 
     242          
     243         ! Local constant initialization 
     244         ! ----------------------------- 
     245         ztfw(1,:,:) = 0.e0     ;     ztfw(jpi,:,:) = 0.e0 
     246          
     247         ! Vertical fluxes 
     248         ! --------------- 
     249          
     250         ! Surface and bottom vertical fluxes set to zero 
     251         ztfw(:,:, 1 ) = 0.e0      ;      ztfw(:,:,jpk) = 0.e0 
     252          
     253         ! interior (2=<jk=<jpk-1) 
     254         DO jk = 2, jpkm1 
     255            DO jj = 2, jpjm1 
     256               DO ji = fs_2, fs_jpim1   ! vector opt. 
     257                  zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk) 
     258                  ! 
     259                  zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      & 
     260                     &            + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.  ) 
     261                  zmskv = 1./MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      & 
     262                     &            + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  ) 
     263                  ! 
     264                  zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) 
     265                  zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
     266                  ! 
     267                  ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
     268                     &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
     269                     &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
     270                     &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
     271               END DO 
    155272            END DO 
    156273         END DO 
    157       ENDIF 
    158  
    159       !!---------------------------------------------------------------------- 
    160       !!   II - horizontal trend of T & S (full) 
    161       !!---------------------------------------------------------------------- 
    162        
    163 #if defined key_diaar5 
    164 !CDIR PARALLEL DO PRIVATE( zdk1t, zdk1s, zfsu )  
    165 #else 
    166 !CDIR PARALLEL DO PRIVATE( zdk1t, zdk1s, zftu, zfsu )  
    167 #endif 
    168       !                                                ! =============== 
    169       DO jk = 1, jpkm1                                 ! Horizontal slab 
    170          !                                             ! =============== 
    171          ! 1. Vertical tracer gradient at level jk and jk+1 
    172          ! ------------------------------------------------ 
    173          ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
    174  
    175          zdk1t(:,:) = ( tb(:,:,jk) - tb(:,:,jk+1) ) * tmask(:,:,jk+1) 
    176          zdk1s(:,:) = ( sb(:,:,jk) - sb(:,:,jk+1) ) * tmask(:,:,jk+1) 
    177  
    178          IF( jk == 1 ) THEN 
    179             zdkt(:,:) = zdk1t(:,:) 
    180             zdks(:,:) = zdk1s(:,:) 
    181          ELSE 
    182             zdkt(:,:) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 
    183             zdks(:,:) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 
    184          ENDIF 
    185  
    186  
    187          ! 2. Horizontal fluxes 
    188          ! -------------------- 
    189  
    190          DO jj = 1 , jpjm1 
    191             DO ji = 1, fs_jpim1   ! vector opt. 
    192                zabe1 = ( fsahtu(ji,jj,jk) + ahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) 
    193                zabe2 = ( fsahtv(ji,jj,jk) + ahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) 
    194  
    195                zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   & 
    196                   &             + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. ) 
    197  
    198                zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   & 
    199                   &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. ) 
    200  
    201                zcof1 = -fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
    202                zcof2 = -fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
    203  
    204 #if defined key_diaar5 
    205                zftu(ji,jj,jk) = (  zabe1 * zdit(ji,jj,jk)   & 
    206 #else 
    207                zftu(ji,jj   ) = (  zabe1 * zdit(ji,jj,jk)   & 
    208 #endif 
    209                   &              + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
    210                   &                         + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
    211                zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
    212                   &              + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
    213                   &                         + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk) 
    214                zfsu(ji,jj   ) = (  zabe1 * zdis(ji,jj,jk)   & 
    215                   &              + zcof1 * (  zdks (ji+1,jj) + zdk1s(ji,jj)      & 
    216                   &                         + zdk1s(ji+1,jj) + zdks (ji,jj)  )  ) * umask(ji,jj,jk) 
    217                zfsv(ji,jj,jk) = (  zabe2 * zdjs(ji,jj,jk)   & 
    218                   &              + zcof2 * (  zdks (ji,jj+1) + zdk1s(ji,jj)      & 
    219                   &                         + zdk1s(ji,jj+1) + zdks (ji,jj)  )  ) * vmask(ji,jj,jk) 
     274          
     275          
     276         ! I.5 Divergence of vertical fluxes added to the general tracer trend 
     277         ! ------------------------------------------------------------------- 
     278         DO jk = 1, jpkm1 
     279            DO jj = 2, jpjm1 
     280               DO ji = fs_2, fs_jpim1   ! vector opt. 
     281                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     282                  ztra = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr 
     283                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     284               END DO 
    220285            END DO 
    221286         END DO 
    222  
    223  
    224          ! II.4 Second derivative (divergence) and add to the general trend 
    225          ! ---------------------------------------------------------------- 
    226          DO jj = 2 , jpjm1 
    227             DO ji = fs_2, fs_jpim1   ! vector opt. 
    228                zbtr= 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) 
    229 #if defined key_diaar5 
    230                zta = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  ) 
    231 #else 
    232                zta = zbtr * ( zftu(ji,jj   ) - zftu(ji-1,jj   ) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  ) 
    233 #endif 
    234                zsa = zbtr * ( zfsu(ji,jj   ) - zfsu(ji-1,jj   ) + zfsv(ji,jj,jk) - zfsv(ji,jj-1,jk)  ) 
    235                ta (ji,jj,jk) = ta (ji,jj,jk) + zta 
    236                sa (ji,jj,jk) = sa (ji,jj,jk) + zsa 
    237             END DO 
    238          END DO 
    239          !                                          ! =============== 
    240       END DO                                        !   End of slab   
    241       !                                             ! =============== 
    242  
    243       IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN   ! Poleward diffusive heat and salt transports 
    244          pht_ldf(:) = ptr_vj( zftv(:,:,:) ) 
    245          pst_ldf(:) = ptr_vj( zfsv(:,:,:) ) 
    246       ENDIF 
    247 #if defined key_diaar5 
    248       zztmp = 0.5 * rau0 * rcp  
    249       z2d(:,:) = 0.e0  
    250       DO jk = 1, jpkm1 
    251          DO jj = 2, jpjm1 
    252             DO ji = fs_2, fs_jpim1   ! vector opt. 
    253                z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
    254             END DO 
    255          END DO 
    256       END DO 
    257       z2d(:,:) =  z2d(:,:) * zztmp 
    258       CALL lbc_lnk( z2d, 'U', -1. ) 
    259       CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
    260       z2d(:,:) = 0.e0  
    261       DO jk = 1, jpkm1 
    262          DO jj = 2, jpjm1 
    263             DO ji = fs_2, fs_jpim1   ! vector opt. 
    264                z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
    265             END DO 
    266          END DO 
    267       END DO 
    268       z2d(:,:) =  z2d(:,:) * zztmp 
    269       CALL lbc_lnk( z2d, 'V', -1. ) 
    270       CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
    271 #endif 
    272  
    273       !!---------------------------------------------------------------------- 
    274       !!   III - vertical trend of T & S (extra diagonal terms only) 
    275       !!---------------------------------------------------------------------- 
    276  
    277       ! Local constant initialization 
    278       ! ----------------------------- 
    279       ztfw(1,:,:) = 0.e0     ;     ztfw(jpi,:,:) = 0.e0 
    280       zsfw(1,:,:) = 0.e0     ;     zsfw(jpi,:,:) = 0.e0 
    281  
    282  
    283       ! Vertical fluxes 
    284       ! --------------- 
    285  
    286       ! Surface and bottom vertical fluxes set to zero 
    287       ztfw(:,:, 1 ) = 0.e0      ;      ztfw(:,:,jpk) = 0.e0 
    288       zsfw(:,:, 1 ) = 0.e0      ;      zsfw(:,:,jpk) = 0.e0 
    289  
    290       ! interior (2=<jk=<jpk-1) 
    291       DO jk = 2, jpkm1 
    292          DO jj = 2, jpjm1 
    293             DO ji = fs_2, fs_jpim1   ! vector opt. 
    294                zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk) 
    295  
    296                zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      & 
    297                   &            + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.  ) 
    298  
    299                zmskv = 1./MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      & 
    300                   &            + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  ) 
    301  
    302                zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) 
    303                zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
    304  
    305                ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
    306                   &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
    307                   &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
    308                   &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
    309  
    310                zsfw(ji,jj,jk) = zcoef3 * (   zdis(ji  ,jj  ,jk-1) + zdis(ji-1,jj  ,jk)      & 
    311                   &                        + zdis(ji-1,jj  ,jk-1) + zdis(ji  ,jj  ,jk)  )   & 
    312                   &           + zcoef4 * (   zdjs(ji  ,jj  ,jk-1) + zdjs(ji  ,jj-1,jk)      & 
    313                   &                        + zdjs(ji  ,jj-1,jk-1) + zdjs(ji  ,jj  ,jk)  ) 
    314             END DO 
    315          END DO 
    316       END DO 
    317  
    318  
    319       ! I.5 Divergence of vertical fluxes added to the general tracer trend 
    320       ! ------------------------------------------------------------------- 
    321  
    322       DO jk = 1, jpkm1 
    323          DO jj = 2, jpjm1 
    324             DO ji = fs_2, fs_jpim1   ! vector opt. 
    325                zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) 
    326                zta  = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr 
    327                zsa  = (  zsfw(ji,jj,jk) - zsfw(ji,jj,jk+1)  ) * zbtr 
    328                ta(ji,jj,jk) = ta(ji,jj,jk) + zta 
    329                sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 
    330             END DO 
    331          END DO 
     287         ! 
    332288      END DO 
    333289      ! 
     
    339295   !!---------------------------------------------------------------------- 
    340296CONTAINS 
    341    SUBROUTINE tra_ldf_iso( kt )               ! Empty routine 
    342       WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt 
     297   SUBROUTINE tra_ldf_iso( kt, cdtype, pgu, pgv, ptb, pta, kjpt, pahtb0 )      ! Empty routine 
     298      CHARACTER(len=3) ::   cdtype 
     299      REAL, DIMENSION(:,:,:) ::   pgu, pgv   ! tracer gradient at pstep levels 
     300      REAL, DIMENSION(:,:,:,:) ::   ptb, pta 
     301      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt, cdtype, pgu(1,1,1), pgv(1,1,1),   & 
     302         &                                                             ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0 
    343303   END SUBROUTINE tra_ldf_iso 
    344304#endif 
Note: See TracChangeset for help on using the changeset viewer.