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 10364 for NEMO/trunk/src/OCE/DYN/dynzdf.F90 – NEMO

Ignore:
Timestamp:
2018-11-30T18:42:51+01:00 (5 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/DYN/dynzdf.F90

    r10068 r10364  
    7171      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars 
    7272      REAL(wp) ::   zzws, ze3va        !   -      - 
     73      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      - 
     74      REAL(wp) ::   zWu , zWv          !   -      - 
     75      REAL(wp) ::   zWui, zWvi         !   -      - 
     76      REAL(wp) ::   zWus, zWvs         !   -      - 
    7377      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace  
    7478      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      - 
    … …  
    155159      !                    !* Matrix construction 
    156160      zdt = r2dt * 0.5 
    157       SELECT CASE( nldf_dyn ) 
    158       CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
    159          DO jk = 1, jpkm1 
    160             DO jj = 2, jpjm1  
    161                DO ji = fs_2, fs_jpim1   ! vector opt. 
    162                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point 
    163                   zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
    164                      &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    165                   zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    166                      &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    167                   zwi(ji,jj,jk) = zzwi 
    168                   zws(ji,jj,jk) = zzws 
    169                   zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    170                END DO 
    171             END DO 
    172          END DO 
    173       CASE DEFAULT               ! iso-level lateral mixing 
    174          DO jk = 1, jpkm1 
    175             DO jj = 2, jpjm1  
    176                DO ji = fs_2, fs_jpim1   ! vector opt. 
    177                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point 
    178                   zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    179                   zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    180                   zwi(ji,jj,jk) = zzwi 
    181                   zws(ji,jj,jk) = zzws 
    182                   zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    183                END DO 
    184             END DO 
    185          END DO 
    186       END SELECT 
    187       ! 
    188       DO jj = 2, jpjm1     !* Surface boundary conditions 
    189          DO ji = fs_2, fs_jpim1   ! vector opt. 
    190             zwi(ji,jj,1) = 0._wp 
    191             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    192          END DO 
    193       END DO 
     161      IF( ln_zad_Aimp ) THEN   !! 
     162         SELECT CASE( nldf_dyn ) 
     163         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
     164            DO jk = 1, jpkm1 
     165               DO jj = 2, jpjm1  
     166                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     167                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     168                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     169                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     170                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     171                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     172                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) 
     173                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 
     174                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )  
     175                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
     176                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     177                  END DO 
     178               END DO 
     179            END DO 
     180         CASE DEFAULT               ! iso-level lateral mixing 
     181            DO jk = 1, jpkm1 
     182               DO jj = 2, jpjm1  
     183                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     184                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     185                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     186                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     187                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) 
     188                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 
     189                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
     190                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
     191                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     192                  END DO 
     193               END DO 
     194            END DO 
     195         END SELECT 
     196         DO jj = 2, jpjm1     !* Surface boundary conditions 
     197            DO ji = fs_2, fs_jpim1   ! vector opt. 
     198               zwi(ji,jj,1) = 0._wp 
     199               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
     200               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2) 
     201               zWus = 0.5_wp * ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) 
     202               zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 
     203               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) 
     204            END DO 
     205         END DO 
     206      ELSE 
     207         SELECT CASE( nldf_dyn ) 
     208         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
     209            DO jk = 1, jpkm1 
     210               DO jj = 2, jpjm1  
     211                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     212                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     213                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     214                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     215                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     216                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     217                     zwi(ji,jj,jk) = zzwi 
     218                     zws(ji,jj,jk) = zzws 
     219                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     220                  END DO 
     221               END DO 
     222            END DO 
     223         CASE DEFAULT               ! iso-level lateral mixing 
     224            DO jk = 1, jpkm1 
     225               DO jj = 2, jpjm1  
     226                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     227                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     228                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     229                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     230                     zwi(ji,jj,jk) = zzwi 
     231                     zws(ji,jj,jk) = zzws 
     232                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     233                  END DO 
     234               END DO 
     235            END DO 
     236         END SELECT 
     237         DO jj = 2, jpjm1     !* Surface boundary conditions 
     238            DO ji = fs_2, fs_jpim1   ! vector opt. 
     239               zwi(ji,jj,1) = 0._wp 
     240               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     241            END DO 
     242         END DO 
     243      ENDIF 
     244      ! 
    194245      ! 
    195246      !              !==  Apply semi-implicit bottom friction  ==! 
    … …  
    274325      !                       !* Matrix construction 
    275326      zdt = r2dt * 0.5 
    276       SELECT CASE( nldf_dyn ) 
    277       CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
    278          DO jk = 1, jpkm1 
    279             DO jj = 2, jpjm1    
    280                DO ji = fs_2, fs_jpim1   ! vector opt. 
    281                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
    282                   zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
    283                      &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    284                   zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    285                      &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    286                   zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
    287                   zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
    288                   zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    289                END DO 
    290             END DO 
    291          END DO 
    292       CASE DEFAULT               ! iso-level lateral mixing 
    293          DO jk = 1, jpkm1 
    294             DO jj = 2, jpjm1    
    295                DO ji = fs_2, fs_jpim1   ! vector opt. 
    296                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
    297                   zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    298                   zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    299                   zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
    300                   zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
    301                   zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    302                END DO 
    303             END DO 
    304          END DO 
    305       END SELECT 
    306       ! 
    307       DO jj = 2, jpjm1        !* Surface boundary conditions 
    308          DO ji = fs_2, fs_jpim1   ! vector opt. 
    309             zwi(ji,jj,1) = 0._wp 
    310             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    311          END DO 
    312       END DO 
     327      IF( ln_zad_Aimp ) THEN   !! 
     328         SELECT CASE( nldf_dyn ) 
     329         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv) 
     330            DO jk = 1, jpkm1 
     331               DO jj = 2, jpjm1  
     332                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     333                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
     334                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     335                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     336                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     337                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     338                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  ) 
     339                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) 
     340                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 
     341                     zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 
     342                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     343                  END DO 
     344               END DO 
     345            END DO 
     346         CASE DEFAULT               ! iso-level lateral mixing 
     347            DO jk = 1, jpkm1 
     348               DO jj = 2, jpjm1  
     349                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     350                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
     351                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     352                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     353                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  ) 
     354                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) 
     355                     zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp ) 
     356                     zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp ) 
     357                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     358                  END DO 
     359               END DO 
     360            END DO 
     361         END SELECT 
     362         DO jj = 2, jpjm1     !* Surface boundary conditions 
     363            DO ji = fs_2, fs_jpim1   ! vector opt. 
     364               zwi(ji,jj,1) = 0._wp 
     365               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
     366               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2) 
     367               zWvs = 0.5_wp * ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) 
     368               zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 
     369               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) 
     370            END DO 
     371         END DO 
     372      ELSE 
     373         SELECT CASE( nldf_dyn ) 
     374         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
     375            DO jk = 1, jpkm1 
     376               DO jj = 2, jpjm1    
     377                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     378                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
     379                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     380                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     381                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     382                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     383                     zwi(ji,jj,jk) = zzwi 
     384                     zws(ji,jj,jk) = zzws 
     385                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     386                  END DO 
     387               END DO 
     388            END DO 
     389         CASE DEFAULT               ! iso-level lateral mixing 
     390            DO jk = 1, jpkm1 
     391               DO jj = 2, jpjm1    
     392                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     393                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
     394                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     395                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     396                     zwi(ji,jj,jk) = zzwi 
     397                     zws(ji,jj,jk) = zzws 
     398                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     399                  END DO 
     400               END DO 
     401            END DO 
     402         END SELECT 
     403         DO jj = 2, jpjm1        !* Surface boundary conditions 
     404            DO ji = fs_2, fs_jpim1   ! vector opt. 
     405               zwi(ji,jj,1) = 0._wp 
     406               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     407            END DO 
     408         END DO 
     409      ENDIF 
     410      ! 
    313411      !              !==  Apply semi-implicit top/bottom friction  ==! 
    314412      ! 
Note: See TracChangeset for help on using the changeset viewer.