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 11407 for NEMO/trunk/src/OCE/TRA – NEMO

Ignore:
Timestamp:
2019-08-06T15:36:27+02:00 (5 years ago)
Author:
acc
Message:

Final documentation tweaks to the adaptive-implicit vertical advection section of chap_ZDF.tex and associated code changes. The code changes are mostly cosmetic or to align with the documentation but this does include a modified version of traadv_fct.F90 (thanks to Jerome) which correctly accounts for any implicit component of the vertical velocity when computing anti-diffusive fluxes for flux correction. The changes have null effect when ln_zad_Aimp is .false. and have been SETTE tested with ORCA2_ICE_PISCES with ln_zad_Aimp = .true.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/OCE/TRA/traadv_fct.F90

    r10425 r11407  
    2121   USE diaar5         ! AR5 diagnostics 
    2222   USE phycst  , ONLY : rau0_rcp 
     23   USE zdf_oce , ONLY : ln_zad_Aimp 
    2324   ! 
    2425   USE in_out_manager ! I/O manager 
    … …  
    8687      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 
    8788      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry 
     89      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup 
     90      LOGICAL  ::   ll_zAimp                                 ! flag to apply adaptive implicit vertical advection 
    8891      !!---------------------------------------------------------------------- 
    8992      ! 
    … …  
    97100      l_hst = .FALSE. 
    98101      l_ptr = .FALSE. 
     102      ll_zAimp = .FALSE. 
    99103      IF( ( cdtype =='TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )       l_trd = .TRUE. 
    100104      IF(   cdtype =='TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE.  
    … …  
    116120      ! 
    117121      zwi(:,:,:) = 0._wp         
     122      ! 
     123      ! If adaptive vertical advection, check if it is needed on this PE at this time 
     124      IF( ln_zad_Aimp ) THEN 
     125         IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE. 
     126      END IF 
     127      ! If active adaptive vertical advection, build tridiagonal matrix 
     128      IF( ll_zAimp ) THEN 
     129         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 
     130         DO jk = 1, jpkm1 
     131            DO jj = 2, jpjm1 
     132               DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.) 
     133                  zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t_a(ji,jj,jk) 
     134                  zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t_a(ji,jj,jk) 
     135                  zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t_a(ji,jj,jk) 
     136               END DO 
     137            END DO 
     138         END DO 
     139      END IF 
    118140      ! 
    119141      DO jn = 1, kjpt            !==  loop over the tracers  ==! 
    … …  
    169191            END DO 
    170192         END DO 
     193          
     194         IF ( ll_zAimp ) THEN 
     195            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) 
     196            ! 
     197            ztu(:,:,1) = 0._wp; ztu(:,:,jpk) = 0._wp 
     198            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
     199               DO jj = 1, jpj 
     200                  DO ji = 1, jpi 
     201                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     202                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     203                     ztu(ji,jj,jk) =  0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
     204                  END DO 
     205               END DO 
     206            END DO 
     207            DO jk = 1, jpkm1 
     208               DO jj = 2, jpjm1 
     209                  DO ji = fs_2, fs_jpim1   ! vector opt.   
     210                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztu(ji,jj,jk) - ztu(ji  ,jj  ,jk+1) ) & 
     211                        &                                  * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     212                  END DO 
     213               END DO 
     214            END DO 
     215            zwz(:,:,:) = zwz(:,:,:) + ztu(:,:,:) 
     216            ! 
     217         END IF 
    171218         !                 
    172219         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
    … …  
    280327         CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1.,  zwz, 'W',  1. ) 
    281328         ! 
     329         !          
     330         IF ( ll_zAimp ) THEN 
     331            DO jk = 1, jpkm1     !* trend and after field with monotonic scheme 
     332               DO jj = 2, jpjm1 
     333                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     334                     !                             ! total intermediate advective trends 
     335                     ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     336                        &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     337                        &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     338                     ztu(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
     339                  END DO 
     340               END DO 
     341            END DO 
     342            ! 
     343            CALL tridia_solver( zwdia, zwsup, zwinf, ztu, ztu , 0 ) 
     344            ! 
     345            ztu(:,:,1) = 0._wp 
     346            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
     347               DO jj = 2, jpjm1 
     348                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     349                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     350                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     351                     zwz(ji,jj,jk) =  zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztu(ji,jj,jk) + zfm_wk * ztu(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
     352                  END DO 
     353               END DO 
     354            END DO 
     355         END IF 
    282356         !        !==  monotonicity algorithm  ==! 
    283357         ! 
    … …  
    289363            DO jj = 2, jpjm1 
    290364               DO ji = fs_2, fs_jpim1   ! vector opt.   
    291                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    292                      &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    293                      &                                   + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) & 
    294                      &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    295                END DO 
    296             END DO 
    297          END DO 
     365                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     366                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     367                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     368                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) 
     369                  zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
     370               END DO 
     371            END DO 
     372         END DO 
     373         ! 
     374         IF ( ll_zAimp ) THEN 
     375            ! 
     376            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
     377               DO jj = 2, jpjm1 
     378                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     379                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     380                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     381                     zwz(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
     382                  END DO 
     383               END DO 
     384            END DO 
     385            DO jk = 1, jpkm1 
     386               DO jj = 2, jpjm1 
     387                  DO ji = fs_2, fs_jpim1   ! vector opt.   
     388                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) & 
     389                        &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     390                  END DO 
     391               END DO 
     392            END DO 
     393         END IF          
    298394         ! 
    299395         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport 
    … …  
    318414      END DO                     ! end of tracer loop 
    319415      ! 
     416      IF ( ll_zAimp ) THEN 
     417         DEALLOCATE( zwdia, zwinf, zwsup ) 
     418      ENDIF 
    320419      IF( l_trd .OR. l_hst ) THEN  
    321420         DEALLOCATE( ztrdx, ztrdy, ztrdz ) 
Note: See TracChangeset for help on using the changeset viewer.