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

Changeset 10364 for 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.

Location:
NEMO/trunk
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/cfgs/SHARED/field_def_nemo-oce.xml

    r9990 r10364  
    7777        <!-- variables available with MLE --> 
    7878        <field id="Lf_NHpf"      long_name="MLE: Lf = N H / f"   unit="m" /> 
     79 
     80        <!-- next variables available with ln_zad_Aimp=.true. --> 
     81        <field id="Courant"      long_name="Courant number"                                                                 unit="#"   grid_ref="grid_T_3D" /> 
     82        <field id="wimp"         long_name="Implicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" /> 
     83        <field id="wexp"         long_name="Explicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" /> 
     84        <field id="wi_cff"       long_name="Fraction of implicit vertical velocity"                                         unit="#"   grid_ref="grid_T_3D" /> 
    7985 
    8086        <!-- next variables available with key_diahth --> 

  • NEMO/trunk/cfgs/SHARED/namelist_ref

    r10190 r10364  
    972972&namzdf        !   vertical physics manager                             (default: NO selection) 
    973973!----------------------------------------------------------------------- 
     974   !                       ! adaptive-implicit vertical advection 
     975   ln_zad_Aimp = .false.      !  Courant number dependent scheme (Shchepetkin 2015) 
     976   ! 
    974977   !                       ! type of vertical closure (required) 
    975978   ln_zdfcst   = .false.      !  constant mixing 

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

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

    r10068 r10364  
    2828#endif 
    2929   ! 
     30   USE iom  
    3031   USE in_out_manager ! I/O manager 
    3132   USE restart        ! only for lrst_oce 
    … …  
    4142   PUBLIC   ssh_nxt    ! called by step.F90 
    4243   PUBLIC   wzv        ! called by step.F90 
     44   PUBLIC   wAimp      ! called by step.F90 
    4345   PUBLIC   ssh_swp    ! called by step.F90 
    4446 
    … …  
    265267   END SUBROUTINE ssh_swp 
    266268 
     269   SUBROUTINE wAimp( kt ) 
     270      !!---------------------------------------------------------------------- 
     271      !!                ***  ROUTINE wAimp  *** 
     272      !!                    
     273      !! ** Purpose :   compute the Courant number and partition vertical velocity 
     274      !!                if a proportion needs to be treated implicitly 
     275      !! 
     276      !! ** Method  : -  
     277      !! 
     278      !! ** action  :   wn      : now vertical velocity (to be handled explicitly) 
     279      !!            :   wi      : now vertical velocity (for implicit treatment) 
     280      !! 
     281      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     282      !!---------------------------------------------------------------------- 
     283      INTEGER, INTENT(in) ::   kt   ! time step 
     284      ! 
     285      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     286      REAL(wp)             ::   zCu, zcff, z1_e3w                     ! local scalars 
     287      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters 
     288      REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters 
     289      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters 
     290      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters 
     291      !!---------------------------------------------------------------------- 
     292      ! 
     293      IF( ln_timing )   CALL timing_start('wAimp') 
     294      ! 
     295      IF( kt == nit000 ) THEN 
     296         IF(lwp) WRITE(numout,*) 
     297         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 
     298         IF(lwp) WRITE(numout,*) '~~~~~ ' 
     299         ! 
     300         Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all) 
     301      ENDIF 
     302      ! 
     303      DO jk = 1, jpkm1            ! calculate Courant numbers 
     304         DO jj = 2, jpjm1 
     305            DO ji = 2, fs_jpim1   ! vector opt. 
     306               z1_e3w = 1._wp / e3w_n(ji,jj,jk) 
     307               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    & 
     308                  &                      + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   & 
     309                  &                          MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   & 
     310                  &                        * r1_e1e2t(ji,jj)                                                  & 
     311                  &                      + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   & 
     312                  &                          MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   & 
     313                  &                        * r1_e1e2t(ji,jj)                                                  & 
     314                  &                      ) * z1_e3w 
     315            END DO 
     316         END DO 
     317      END DO 
     318      ! 
     319      CALL iom_put("Courant",Cu_adv) 
     320      ! 
     321      wi(:,:,:) = 0._wp                                 ! Includes top and bottom values set to zero 
     322      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere 
     323         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition 
     324            DO jj = 2, jpjm1                            ! w where necessary 
     325               DO ji = 2, fs_jpim1   ! vector opt. 
     326                  ! 
     327                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) ) 
     328                  ! 
     329                  IF( zCu < Cu_min ) THEN               !<-- Fully explicit 
     330                     zcff = 0._wp 
     331                  ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit 
     332                     zcff = ( zCu - Cu_min )**2 
     333                     zcff = zcff / ( Fcu + zcff ) 
     334                  ELSE                                  !<-- Mostly implicit 
     335                     zcff = ( zCu - Cu_max )/ zCu 
     336                  ENDIF 
     337                  zcff = MIN(1._wp, zcff) 
     338                  ! 
     339                  wi(ji,jj,jk) =           zcff   * wn(ji,jj,jk) 
     340                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 
     341                  ! 
     342                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient 
     343               END DO 
     344            END DO 
     345         END DO 
     346      ELSE 
     347         ! Fully explicit everywhere 
     348         Cu_adv = 0.0_wp                                ! Reuse array to output coefficient 
     349      ENDIF 
     350      CALL iom_put("wimp",wi)  
     351      CALL iom_put("wi_cff",Cu_adv) 
     352      CALL iom_put("wexp",wn) 
     353      ! 
     354      IF( ln_timing )   CALL timing_stop('wAimp') 
     355      ! 
     356   END SUBROUTINE wAimp 
    267357   !!====================================================================== 
    268358END MODULE sshwzv 

  • 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 

  • NEMO/trunk/src/OCE/ZDF/zdf_oce.F90

    r10068 r10364  
    1818 
    1919   !                            !!* namelist namzdf: vertical physics * 
     20   !                             ! Adaptive-implicit vertical advection flag 
     21   LOGICAL , PUBLIC ::   ln_zad_Aimp !: adaptive (Courant number-based) implicit vertical advection 
    2022   !                             ! vertical closure scheme flags 
    2123   LOGICAL , PUBLIC ::   ln_zdfcst   !: constant coefficients 

  • NEMO/trunk/src/OCE/ZDF/zdfosm.F90

    r10069 r10364  
    3232 
    3333   !!---------------------------------------------------------------------- 
    34    !!   'key_zdfosm'                                             OSMOSIS scheme 
     34   !!   'ln_zdfosm'                                             OSMOSIS scheme 
    3535   !!---------------------------------------------------------------------- 
    3636   !!   zdf_osm       : update momentum and tracer Kz from osm scheme 

  • NEMO/trunk/src/OCE/ZDF/zdfphy.F90

    r10069 r10364  
    8080         &             ln_zdfswm,                                    &     ! surface  wave-induced mixing 
    8181         &             ln_zdfiwm,                                    &     ! internal  -      -      - 
     82         &             ln_zad_Aimp,                                  &     ! apdative-implicit vertical advection 
    8283         &             rn_avm0, rn_avt0, nn_avb, nn_havtb                  ! coefficients 
    8384      !!---------------------------------------------------------------------- 
    … …  
    101102      IF(lwp) THEN                      ! Parameter print 
    102103         WRITE(numout,*) '   Namelist namzdf : set vertical mixing mixing parameters' 
     104         WRITE(numout,*) '      adaptive-implicit vertical advection' 
     105         WRITE(numout,*) '         Courant number targeted application   ln_zad_Aimp = ', ln_zad_Aimp 
    103106         WRITE(numout,*) '      vertical closure scheme' 
    104107         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfcst = ', ln_zdfcst 
    … …  
    127130      ENDIF 
    128131 
     132      IF( ln_zad_Aimp ) THEN 
     133         IF( zdf_phy_alloc() /= 0 )   & 
     134        &       CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' ) 
     135         wi(:,:,:) = 0._wp 
     136      ENDIF 
    129137      !                          !==  Background eddy viscosity and diffusivity  ==! 
    130138      IF( nn_avb == 0 ) THEN             ! Define avmb, avtb from namelist parameter 
    … …  
    316324      ! 
    317325   END SUBROUTINE zdf_phy 
     326   INTEGER FUNCTION zdf_phy_alloc() 
     327      !!---------------------------------------------------------------------- 
     328      !!                 ***  FUNCTION zdf_phy_alloc  *** 
     329      !!---------------------------------------------------------------------- 
     330     ! Allocate wi array (declared in oce.F90) for use with the adaptive-implicit vertical velocity option 
     331     ALLOCATE(     wi(jpi,jpj,jpk), Cu_adv(jpi,jpj,jpk),  STAT= zdf_phy_alloc ) 
     332     IF( zdf_phy_alloc /= 0 )   CALL ctl_warn('zdf_phy_alloc: failed to allocate ln_zad_Aimp=T required arrays') 
     333     IF( lk_mpp             )   CALL mpp_sum ( zdf_phy_alloc ) 
     334   END FUNCTION zdf_phy_alloc 
    318335 
    319336   !!====================================================================== 

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

    r10068 r10364  
    2222   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s] 
    2323   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wn             !: vertical velocity            [m/s] 
     24   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wi             !: vertical vel. (adaptive-implicit) [m/s] 
    2425   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           hdivn          !: horizontal divergence        [s-1] 
    2526   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa    !: 4D T-S fields                  [Celsius,psu]  
    … …  
    2930   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0  [no units] 
    3031   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhop   !: potential volumic mass                           [kg/m3] 
     32   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   Cu_adv                   !: vertical Courant number (adaptive-implicit) 
    3133 
    3234   !! free surface                                      !  before  ! now    ! after  ! 

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

    r10068 r10364  
    160160 
    161161                            CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor) 
    162       IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
     162      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
    163163                            CALL wzv           ( kstp )  ! now cross-level velocity  
     164      IF( ln_zad_Aimp )     CALL wAimp         ( kstp )  ! Adaptive-implicit vertical advection partitioning 
    164165                            CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation 
    165166                             
    … …  
    198199         IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component) 
    199200                            CALL wzv        ( kstp )              ! now cross-level velocity  
     201         IF( ln_zad_Aimp )  CALL wAimp      ( kstp )  ! Adaptive-implicit vertical advection partitioning 
    200202      ENDIF 
    201203       

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

    r10068 r10364  
    2424   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2525   USE lib_mpp         ! distributed memory computing 
     26   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables 
    2627   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy 
    2728 
    … …  
    3233   PUBLIC stp_ctl           ! routine called by step.F90 
    3334 
    34    INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, istatus 
     35   INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus 
    3536   !!---------------------------------------------------------------------- 
    3637   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    6869      INTEGER , DIMENSION(3) ::   ilocu, ilocs1, ilocs2 
    6970      INTEGER , DIMENSION(2) ::   iloch 
    70       REAL(wp), DIMENSION(5) ::   zmax 
     71      REAL(wp), DIMENSION(9) ::   zmax 
    7172      CHARACTER(len=20) :: clname 
    7273      !!---------------------------------------------------------------------- 
    … …  
    9091            istatus = NF90_DEF_VAR( idrun,       's_min', NF90_DOUBLE, (/ idtime /), ids1  ) 
    9192            istatus = NF90_DEF_VAR( idrun,       's_max', NF90_DOUBLE, (/ idtime /), ids2  ) 
     93            istatus = NF90_DEF_VAR( idrun,       't_min', NF90_DOUBLE, (/ idtime /), idt1  ) 
     94            istatus = NF90_DEF_VAR( idrun,       't_max', NF90_DOUBLE, (/ idtime /), idt2  ) 
     95            IF( ln_zad_Aimp ) THEN 
     96               istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1  ) 
     97               istatus = NF90_DEF_VAR( idrun,       'Cu_max', NF90_DOUBLE, (/ idtime /), idc1  ) 
     98            ENDIF 
    9299            istatus = NF90_ENDDEF(idrun) 
    93100         ENDIF 
     101         zmax(8:9) = 0._wp    ! initialise to zero in case ln_zad_Aimp option is not in use 
    94102          
    95103      ENDIF 
    … …  
    109117      zmax(3) = MAXVAL( -tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max 
    110118      zmax(4) = MAXVAL(  tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   !       salinity max 
    111       zmax(5) = REAL( nstop , wp )                                            ! stop indicator 
     119      zmax(5) = MAXVAL( -tsn(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max 
     120      zmax(6) = MAXVAL(  tsn(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp )   !       temperature max 
     121      zmax(7) = REAL( nstop , wp )                                            ! stop indicator 
     122      IF( ln_zad_Aimp ) THEN 
     123         zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 
     124         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) !       cell Courant no. max 
     125      ENDIF 
    112126      ! 
    113127      IF( lk_mpp ) THEN 
    114          CALL mpp_max_multiple( zmax(:), 5 )    ! max over the global domain 
     128         CALL mpp_max_multiple( zmax(:), 9 )    ! max over the global domain 
    115129         ! 
    116          nstop = NINT( zmax(5) )                 ! nstop indicator sheared among all local domains 
     130         nstop = NINT( zmax(7) )                 ! nstop indicator sheared among all local domains 
    117131      ENDIF 
    118132      ! 
    … …  
    172186         istatus = NF90_PUT_VAR( idrun,  ids1, (/-zmax(3)/), (/kt/), (/1/) ) 
    173187         istatus = NF90_PUT_VAR( idrun,  ids2, (/ zmax(4)/), (/kt/), (/1/) ) 
     188         istatus = NF90_PUT_VAR( idrun,  idt1, (/-zmax(5)/), (/kt/), (/1/) ) 
     189         istatus = NF90_PUT_VAR( idrun,  idt2, (/ zmax(6)/), (/kt/), (/1/) ) 
     190         IF( ln_zad_Aimp ) THEN 
     191            istatus = NF90_PUT_VAR( idrun,  idw1, (/ zmax(8)/), (/kt/), (/1/) ) 
     192            istatus = NF90_PUT_VAR( idrun,  idc1, (/ zmax(9)/), (/kt/), (/1/) ) 
     193         ENDIF 
    174194         IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun) 
    175195         IF( kt == nitend         ) istatus = NF90_CLOSE(idrun) 
Note: See TracChangeset for help on using the changeset viewer.