Changeset 9976


Ignore:
Timestamp:
2018-07-19T17:58:12+02:00 (2 years ago)
Author:
acc
Message:

Branch: dev_r9956_ENHANCE05_ZAD_AIMP. First set of changes to implement an adaptive-implicit vertical advection option (see ticket #2042). This code compiles and runs but has issues when the new option is activated.

Location:
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/cfgs/SHARED/field_def_nemo-oce.xml

    r9957 r9976  
    7878        <field id="Lf_NHpf"      long_name="MLE: Lf = N H / f"   unit="m" /> 
    7979 
     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" /> 
    8084        <!-- next variables available with key_diahth --> 
    8185        <field id="mlddzt"       long_name="Thermocline Depth (depth of max dT/dz)"         standard_name="depth_at_maximum_upward_derivative_of_sea_water_potential_temperature"             unit="m"                         /> 
  • NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/cfgs/SHARED/namelist_ref

    r9957 r9976  
    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/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/dom_oce.F90

    r9957 r9976  
    135135   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0           ,   e3f_n            !: f- vert. scale factor [m] 
    136136   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0 ,   e3w_b ,   e3w_n            !: w- vert. scale factor [m] 
    137    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m] 
    138    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m] 
     137   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0 ,  e3uw_b ,  e3uw_n , e3uw_a   !: uw-vert. scale factor [m] 
     138   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0 ,  e3vw_b ,  e3vw_n , e3vw_a   !: vw-vert. scale factor [m] 
    139139 
    140140   !                                                        !  ref.   ! before  !   now   ! 
     
    267267         &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,         & 
    268268         &      e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,         &                
    269          &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) )                        
     269         &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,         &                
     270         &      e3uw_a(jpi,jpj,jpk) , e3vw_a(jpi,jpj,jpk) ,     STAT=ierr(5) )                        
    270271         ! 
    271272      ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           & 
  • NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/domvvl.F90

    r9957 r9976  
    149149 
    150150      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 
    151       e3t_a(:,:,:) = e3t_n(:,:,:) 
    152       e3u_a(:,:,:) = e3u_n(:,:,:) 
    153       e3v_a(:,:,:) = e3v_n(:,:,:) 
     151      e3t_a (:,:,:) = e3t_n (:,:,:) 
     152      e3u_a (:,:,:) = e3u_n (:,:,:) 
     153      e3v_a (:,:,:) = e3v_n (:,:,:) 
     154      e3uw_a(:,:,:) = e3uw_n(:,:,:) 
     155      e3vw_a(:,:,:) = e3vw_n(:,:,:) 
    154156      ! 
    155157      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep) 
     
    534536      END IF 
    535537 
    536       ! *********************************** ! 
    537       ! After scale factors at u- v- points ! 
    538       ! *********************************** ! 
     538      ! ******************************************* ! 
     539      ! After scale factors at u- v- uw- vw- points ! 
     540      ! ******************************************* ! 
    539541 
    540542      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 
    541543      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 
     544      CALL dom_vvl_interpol( e3u_a(:,:,:), e3uw_a(:,:,:), 'UW' ) 
     545      CALL dom_vvl_interpol( e3v_a(:,:,:), e3vw_a(:,:,:), 'VW' ) 
    542546 
    543547      ! *********************************** ! 
  • NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90

    r9957 r9976  
    7171      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars 
    7272      REAL(wp) ::   zzws, ze3va        !   -      - 
     73      REAL(wp) ::   z1_e3un, z1_e3vn   !   -      - 
     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 
     161      IF( ln_zad_Aimp ) THEN   !! 
     162         IF( ln_dynadv_vec ) THEN      !==  Vector invariant advection  ==! 
     163            SELECT CASE( nldf_dyn ) 
     164            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
     165               DO jk = 1, jpkm1 
     166                  DO jj = 2, jpjm1  
     167                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     168                        z1_e3un = 1._wp / e3u_n(ji,jj,jk) 
     169                        zzwi = ( ( avm   (ji+1,jj,jk  ) + avm   (ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     170                           &     / e3uw_a(ji  ,jj,jk  ) ) * z1_e3un * wumask(ji,jj,jk  ) 
     171                        zzws = ( ( avm   (ji+1,jj,jk+1) + avm   (ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     172                           &     / e3uw_a(ji  ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1) 
     173                        zWu  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji+1,jj,jk  )   & 
     174                           &               + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1)   ) 
     175                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un ) 
     176                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWu, 0._wp ) * z1_e3un ) 
     177                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un ) 
     178                     END DO 
     179                  END DO 
     180               END DO 
     181            CASE DEFAULT               ! iso-level lateral mixing 
     182               DO jk = 1, jpkm1 
     183                  DO jj = 2, jpjm1  
     184                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     185                        z1_e3un = 1._wp / e3u_n(ji,jj,jk) 
     186                        zzwi = ( ( avm   (ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   & 
     187                           &     / e3uw_a(ji  ,jj,jk  ) ) * z1_e3un * wumask(ji,jj,jk  ) 
     188                        zzws = ( ( avm   (ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   & 
     189                           &     / e3uw_a(ji  ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1) 
     190                        zWu  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji+1,jj,jk  )   & 
     191                           &               + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1)   ) 
     192                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un ) 
     193                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWu, 0._wp ) * z1_e3un ) 
     194                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un ) 
     195                     END DO 
     196                  END DO 
     197               END DO 
     198            END SELECT 
     199         ELSE                          !==  Flux form advection  ==! 
     200            SELECT CASE( nldf_dyn ) 
     201            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
     202               DO jk = 1, jpkm1 
     203                  DO jj = 2, jpjm1  
     204                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     205                        ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     206                        zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     207                           &         / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     208                        zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     209                           &         / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     210                        zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) 
     211                        zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 
     212                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWui, 0._wp ) )  
     213                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWus, 0._wp ) ) 
     214                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWui, 0._wp ) + MIN( zWus, 0._wp ) ) 
     215                     END DO 
     216                  END DO 
     217               END DO 
     218            CASE DEFAULT               ! iso-level lateral mixing 
     219               DO jk = 1, jpkm1 
     220                  DO jj = 2, jpjm1  
     221                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     222                        ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     223                        zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     224                        zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     225                        zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) 
     226                        zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 
     227                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWui, 0._wp ) )  
     228                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWus, 0._wp ) ) 
     229                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWui, 0._wp ) + MIN( zWus, 0._wp ) ) 
     230                     END DO 
     231                  END DO 
     232               END DO 
     233            END SELECT 
     234         ENDIF 
     235      ELSE 
     236         SELECT CASE( nldf_dyn ) 
     237         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
     238            DO jk = 1, jpkm1 
     239               DO jj = 2, jpjm1  
     240                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     241                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     242                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     243                        &         / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     244                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     245                        &         / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     246                     zwi(ji,jj,jk) = zzwi 
     247                     zws(ji,jj,jk) = zzws 
     248                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     249                  END DO 
     250               END DO 
     251            END DO 
     252         CASE DEFAULT               ! iso-level lateral mixing 
     253            DO jk = 1, jpkm1 
     254               DO jj = 2, jpjm1  
     255                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     256                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     257                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     258                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     259                     zwi(ji,jj,jk) = zzwi 
     260                     zws(ji,jj,jk) = zzws 
     261                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     262                  END DO 
     263               END DO 
     264            END DO 
     265         END SELECT 
     266      ENDIF 
    187267      ! 
    188268      DO jj = 2, jpjm1     !* Surface boundary conditions 
     
    274354      !                       !* Matrix construction 
    275355      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 
     356      IF( ln_zad_Aimp ) THEN   !! 
     357         IF( ln_dynadv_vec ) THEN      !==  Vector invariant advection  ==! 
     358            SELECT CASE( nldf_dyn ) 
     359            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv) 
     360               DO jk = 1, jpkm1 
     361                  DO jj = 2, jpjm1  
     362                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     363                        z1_e3vn = 1._wp / e3v_n(ji,jj,jk) 
     364                        zzwi = ( ( avm   (ji,jj+1,jk  ) + avm   (ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     365                           &     / e3vw_a(ji,jj  ,jk  ) ) * z1_e3vn * wvmask(ji,jj,jk  ) 
     366                        zzws = ( ( avm   (ji,jj+1,jk+1) + avm   (ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     367                           &     / e3vw_a(ji,jj  ,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1) 
     368                        zWv  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji,jj+1,jk  )   & 
     369                           &               + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1)   ) 
     370                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn ) 
     371                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn ) 
     372                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn ) 
     373                     END DO 
     374                  END DO 
     375               END DO 
     376            CASE DEFAULT               ! iso-level lateral mixing 
     377               DO jk = 1, jpkm1 
     378                  DO jj = 2, jpjm1  
     379                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     380                        z1_e3vn = 1._wp / e3v_n(ji,jj,jk) 
     381                        zzwi = ( ( avm   (ji,jj+1,jk  ) + avm(ji,jj,jk  ) )   & 
     382                           &     / e3vw_a(ji,jj  ,jk  ) ) * z1_e3vn * wvmask(ji,jj,jk  ) 
     383                        zzws = ( ( avm   (ji,jj+1,jk+1) + avm(ji,jj,jk+1) )   & 
     384                           &     / e3vw_a(ji  ,jj,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1) 
     385                        zWv  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji,jj+1,jk  )   & 
     386                           &               + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1)   ) 
     387                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn ) 
     388                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn ) 
     389                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn ) 
     390                     END DO 
     391                  END DO 
     392               END DO 
     393            END SELECT 
     394         ELSE                          !==  Flux form advection  ==! 
     395            SELECT CASE( nldf_dyn ) 
     396            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv) 
     397               DO jk = 1, jpkm1 
     398                  DO jj = 2, jpjm1  
     399                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     400                        ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at U-point 
     401                        zzwi = ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     402                           &         / ( ze3va * e3vw_a(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     403                        zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     404                           &         / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     405                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) 
     406                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) 
     407                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWvi, 0._wp ) )  
     408                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWvs, 0._wp ) ) 
     409                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     410                     END DO 
     411                  END DO 
     412               END DO 
     413            CASE DEFAULT               ! iso-level lateral mixing 
     414               DO jk = 1, jpkm1 
     415                  DO jj = 2, jpjm1  
     416                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     417                        ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at U-point 
     418                        zzwi = ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_a(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     419                        zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     420                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) 
     421                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) 
     422                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWvi, 0._wp ) )  
     423                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWvs, 0._wp ) ) 
     424                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     425                     END DO 
     426                  END DO 
     427               END DO 
     428            END SELECT 
     429         ENDIF 
     430      ELSE 
     431         SELECT CASE( nldf_dyn ) 
     432         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
     433            DO jk = 1, jpkm1 
     434               DO jj = 2, jpjm1    
     435                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     436                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
     437                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     438                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     439                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     440                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     441                     zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
     442                     zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
     443                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     444                  END DO 
     445               END DO 
     446            END DO 
     447         CASE DEFAULT               ! iso-level lateral mixing 
     448            DO jk = 1, jpkm1 
     449               DO jj = 2, jpjm1    
     450                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     451                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
     452                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     453                     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) 
     454                     zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
     455                     zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
     456                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     457                  END DO 
     458               END DO 
     459            END DO 
     460         END SELECT 
     461      ENDIF 
    306462      ! 
    307463      DO jj = 2, jpjm1        !* Surface boundary conditions 
  • NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/sshwzv.F90

    r9957 r9976  
    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)             ::   zcff_min, zcff_max                    ! local scalars 
     288      REAL(wp) , PARAMETER ::   Cu_min = 0.6_wp                       ! local parameters 
     289      REAL(wp) , PARAMETER ::   Cu_max = 0.95_wp                      ! local parameters 
     290      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters 
     291      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters 
     292      !!---------------------------------------------------------------------- 
     293      ! 
     294      IF( ln_timing )   CALL timing_start('wAimp') 
     295      ! 
     296      IF( kt == nit000 ) THEN 
     297         IF(lwp) WRITE(numout,*) 
     298         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 
     299         IF(lwp) WRITE(numout,*) '~~~~~ ' 
     300         ! 
     301         wi    (:,:,jpk) = 0._wp              ! bottom value : wi=0 (set once for all) 
     302         Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all) 
     303      ENDIF 
     304      ! 
     305      DO jk = 1, jpkm1            ! calculate Courant numbers 
     306         DO jj = 2, jpjm1 
     307            DO ji = 2, fs_jpim1   ! vector opt. 
     308               z1_e3w = 1._wp / e3w_n(ji,jj,jk) 
     309               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    & 
     310                  &                      + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   & 
     311                  &                          MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   & 
     312                  &                        / e2t(ji,jj)                                                       & 
     313                  &                      + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   & 
     314                  &                          MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   & 
     315                  &                        / e1t(ji,jj)                                                       & 
     316                  &                      ) * z1_e3w 
     317            END DO 
     318         END DO 
     319      END DO 
     320      ! 
     321      CALL iom_put("Courant",Cu_adv) 
     322      ! 
     323      wi(:,:,:) = 0._wp 
     324      zcff_min=999. ; zcff_max = -999. 
     325      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere 
     326         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition 
     327            DO jj = 2, jpjm1                            ! w where necessary 
     328               DO ji = 2, fs_jpim1   ! vector opt. 
     329                  ! 
     330                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) ) 
     331                  ! 
     332                  IF( zCu < Cu_min ) THEN               !<-- Fully explicit 
     333                     zcff = 0._wp 
     334                  ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit 
     335                     zcff = ( zCu - Cu_min )**2 
     336                     zcff = zcff / ( Fcu + zcff ) 
     337                  ELSE                                  !<-- Fully implicit 
     338                     zcff = ( zCu - Cu_max )/ zCu 
     339                  ENDIF 
     340                  ! 
     341                  !zcff = 0._wp 
     342                  wi(ji,jj,jk) =           zcff   * wn(ji,jj,jk) 
     343                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 
     344                  zcff_min = MIN( zcff, zcff_min ) 
     345                  zcff_max = MAX( zcff, zcff_max ) 
     346                  ! 
     347               END DO 
     348            END DO 
     349         END DO 
     350      ELSE 
     351         ! Fully explicit everywhere 
     352         write(*,*) 'NoWi ',kt 
     353      ENDIF 
     354         write(*,*) 'ZCFF ',kt,zcff_min,zcff_max 
     355      CALL iom_put("wimp",wi) 
     356      CALL iom_put("wexp",wn) 
     357      ! 
     358      IF( ln_timing )   CALL timing_stop('wAimp') 
     359      ! 
     360   END SUBROUTINE wAimp 
    267361   !!====================================================================== 
    268362END MODULE sshwzv 
  • NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/TRA/trazdf.F90

    r9957 r9976  
    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. 
     183                        zwi(ji,jj,jk) = - p2dt * ( zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )  + MIN( wi(ji,jj,jk  ) , 0._wp ) ) 
     184                        zws(ji,jj,jk) = - p2dt * ( zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)  - MAX( wi(ji,jj,jk+1) , 0._wp ) ) 
     185                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)   & 
     186                           &                   - p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 
     187                    END DO 
     188                  END DO 
     189               END DO 
     190            ELSE 
     191               DO jk = 1, jpkm1 
     192                  DO jj = 2, jpjm1 
     193                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     194                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk) 
     195                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 
     196                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     197                    END DO 
     198                  END DO 
     199               END DO 
     200            ENDIF 
    189201            ! 
    190202            !! Matrix inversion from the first level 
  • NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdf_oce.F90

    r9957 r9976  
    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/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdfosm.F90

    r9957 r9976  
    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/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdfphy.F90

    r9957 r9976  
    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/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/oce.F90

    r9957 r9976  
    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/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/step.F90

    r9957 r9976  
    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       
Note: See TracChangeset for help on using the changeset viewer.