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 11967 for NEMO/branches/2019/ENHANCE-02_ISF_nemo_TEST_MERGE/src/OCE/TRA/traadv_fct.F90 – NEMO

Ignore:
Timestamp:
2019-11-26T15:11:43+01:00 (4 years ago)
Author:
davestorkey
Message:

2019/ENHANCE-02_ISF_nemo_TEST_MERGE : Update to rev 11953.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/ENHANCE-02_ISF_nemo_TEST_MERGE/src/OCE/TRA/traadv_fct.F90

    r10425 r11967  
    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            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 
     198            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
     199               DO jj = 2, jpjm1 
     200                  DO ji = fs_2, fs_jpim1   ! vector opt.   
     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                     ztw(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                     zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 
     205                  END DO 
     206               END DO 
     207            END DO 
     208            DO jk = 1, jpkm1 
     209               DO jj = 2, jpjm1 
     210                  DO ji = fs_2, fs_jpim1   ! vector opt.   
     211                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
     212                        &                                  * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     213                  END DO 
     214               END DO 
     215            END DO 
     216            ! 
     217         END IF 
    171218         !                 
    172219         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
     
    277324            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked 
    278325         ENDIF 
     326         !          
     327         IF ( ll_zAimp ) THEN 
     328            DO jk = 1, jpkm1     !* trend and after field with monotonic scheme 
     329               DO jj = 2, jpjm1 
     330                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     331                     !                             ! total intermediate advective trends 
     332                     ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     333                        &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     334                        &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     335                     ztw(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
     336                  END DO 
     337               END DO 
     338            END DO 
     339            ! 
     340            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 
     341            ! 
     342            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
     343               DO jj = 2, jpjm1 
     344                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     345                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     346                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     347                     zwz(ji,jj,jk) =  zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
     348                  END DO 
     349               END DO 
     350            END DO 
     351         END IF 
    279352         ! 
    280353         CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1.,  zwz, 'W',  1. ) 
     
    289362            DO jj = 2, jpjm1 
    290363               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 
     364                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     365                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     366                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     367                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) 
     368                  zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
     369               END DO 
     370            END DO 
     371         END DO 
     372         ! 
     373         IF ( ll_zAimp ) THEN 
     374            ! 
     375            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 
     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                     ztw(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                     zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 
     383                  END DO 
     384               END DO 
     385            END DO 
     386            DO jk = 1, jpkm1 
     387               DO jj = 2, jpjm1 
     388                  DO ji = fs_2, fs_jpim1   ! vector opt.   
     389                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
     390                        &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     391                  END DO 
     392               END DO 
     393            END DO 
     394         END IF          
    298395         ! 
    299396         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport 
     
    318415      END DO                     ! end of tracer loop 
    319416      ! 
     417      IF ( ll_zAimp ) THEN 
     418         DEALLOCATE( zwdia, zwinf, zwsup ) 
     419      ENDIF 
    320420      IF( l_trd .OR. l_hst ) THEN  
    321421         DEALLOCATE( ztrdx, ztrdy, ztrdz ) 
Note: See TracChangeset for help on using the changeset viewer.