Ignore:
Timestamp:
2018-11-30T18:42:51+01:00 (3 years ago)
Author:
acc
Message:

Introduce Adaptive-Implicit vertical advection option to the trunk. This is code merged from branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP (see ticket #2042). The structure for the option is complete but is currently only successful with the flux-limited advection scheme (ln_traadv_mus). The use of this scheme with flux corrected advection schemes is not recommended until improvements to the nonoscillatory algorithm have been completed (work in progress elsewhere). The scheme is activated via a new namelist switch (ln_zad_Aimp) and is off by default.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/TRA/trazdf.F90

    r10224 r10364  
    136136      ! 
    137137      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    138       REAL(wp) ::  zrhs            ! local scalars 
     138      REAL(wp) ::  zrhs, zzwi, zzws ! local scalars 
    139139      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwt, zwd, zws 
    140140      !!--------------------------------------------------------------------- 
     
    177177            ! 
    178178            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked) 
    179             DO jk = 1, jpkm1 
    180                DO jj = 2, jpjm1 
    181                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    182 !!gm BUG  I think, use e3w_a instead of e3w_n, not sure of that 
    183                      zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  ) 
    184                      zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 
    185                      zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    186                  END DO 
    187                END DO 
    188             END DO 
     179            IF( ln_zad_Aimp ) THEN         ! Adaptive implicit vertical advection 
     180               DO jk = 1, jpkm1 
     181                  DO jj = 2, jpjm1 
     182                     DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.) 
     183                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  ) 
     184                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 
     185                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zzwi - zzws   & 
     186                           &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 
     187                        zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp ) 
     188                        zws(ji,jj,jk) = zzws - p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp ) 
     189                    END DO 
     190                  END DO 
     191               END DO 
     192            ELSE 
     193               DO jk = 1, jpkm1 
     194                  DO jj = 2, jpjm1 
     195                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     196                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk) 
     197                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 
     198                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     199                    END DO 
     200                  END DO 
     201               END DO 
     202            ENDIF 
    189203            ! 
    190204            !! Matrix inversion from the first level 
Note: See TracChangeset for help on using the changeset viewer.