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 – NEMO

Changeset 11407 for NEMO/trunk/src/OCE


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.

Location:
NEMO/trunk/src/OCE
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/OCE/DYN/sshwzv.F90

    r11293 r11407  
    286286      REAL(wp)             ::   zCu, zcff, z1_e3t                     ! local scalars 
    287287      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters 
    288       REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters 
     288      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters 
    289289      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters 
    290290      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters 
    … …  
    346346                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 
    347347                  ! 
    348                   Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient 
     348                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl 
    349349               END DO 
    350350            END DO 
    … …  
    353353      ELSE 
    354354         ! Fully explicit everywhere 
    355          Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient 
     355         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient below and in stp_ctl 
    356356         wi    (:,:,:) = 0._wp 
    357357      ENDIF 

  • 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 ) 

  • NEMO/trunk/src/OCE/stpctl.F90

    r10570 r11407  
    9696            IF( ln_zad_Aimp ) THEN 
    9797               istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1  ) 
    98                istatus = NF90_DEF_VAR( idrun,       'Cu_max', NF90_DOUBLE, (/ idtime /), idc1  ) 
     98               istatus = NF90_DEF_VAR( idrun,       'Cf_max', NF90_DOUBLE, (/ idtime /), idc1  ) 
    9999            ENDIF 
    100100            istatus = NF90_ENDDEF(idrun) 
    … …  
    123123      IF( ln_zad_Aimp ) THEN 
    124124         zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 
    125          zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) !       cell Courant no. max 
     125         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max 
    126126      ENDIF 
    127127      ! 
Note: See TracChangeset for help on using the changeset viewer.