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 15610 for NEMO/releases/r4.0/r4.0-HEAD/tests/CANAL/MY_SRC – NEMO

Ignore:
Timestamp:
2021-12-17T16:09:23+01:00 (3 years ago)
Author:
jchanut
Message:

#1791: add namelist parameter nn_vvl_interp to control how scale factors are set at U-V-F points. nn_vvl_interp=0 is the old interpolation scheme that do not account for steps. nn_vvl_interp=1, corrects the bottom thicknesses, but does not ensure that they get too small for stability. nn_vvl_interp=2 is a "qco like" formulation, for which scale factors anomalies are set proportional to the scale factors at rest. Set nn_vvl_interp=2 as the default.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/releases/r4.0/r4.0-HEAD/tests/CANAL/MY_SRC/domvvl.F90

    r15563 r15610  
    4949   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate 
    5050   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer 
     51   ! 
     52   INTEGER          :: nn_vvl_interp                       ! scale factors anomaly interpolation method at U-V-F points 
     53                                                           ! =0 linear with no bottom correction over steps (old) 
     54                                                           ! =1 linear with bottom correction over steps 
     55                                                           ! =2 "qco like", i.e. proportional to thicknesses at rest 
     56   ! 
    5157   !                                                       ! conservation: not used yet 
    5258   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient 
     
    720726      INTEGER ::   ji, jj, jk                                       ! dummy loop indices 
    721727      INTEGER ::   iku, ikum1, ikv, ikvm1, ikf, ikfm1               !  
    722       REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F 
     728      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/ 
     729      REAL(wp), DIMENSION(jpi,jpj) :: zssh                          ! work array to retrieve ssh (nn_vvl_interp > 1) 
    723730      !!---------------------------------------------------------------------- 
    724731      ! 
     
    732739         ! 
    733740      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    734          DO jk = 1, jpk 
     741         SELECT CASE ( nn_vvl_interp ) 
     742         CASE ( 0 )    
     743            !    
     744            DO jk = 1, jpk 
     745               DO jj = 1, jpjm1 
     746                  DO ji = 1, fs_jpim1   ! vector opt. 
     747                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
     748                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     749                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     750                  END DO 
     751               END DO 
     752            END DO 
     753            ! 
     754         CASE ( 1 )  
     755            !    
     756            DO jk = 1, jpk 
     757               DO jj = 1, jpjm1 
     758                  DO ji = 1, fs_jpim1   ! vector opt. 
     759                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
     760                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     761                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     762                  END DO 
     763               END DO 
     764            END DO 
     765            ! 
     766            ! Bottom correction: 
    735767            DO jj = 1, jpjm1 
    736768               DO ji = 1, fs_jpim1   ! vector opt. 
    737                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
    738                      &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    739                      &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    740                END DO 
    741             END DO 
    742          END DO 
    743          ! 
    744          ! Bottom correction: 
    745          DO jj = 1, jpjm1 
    746             DO ji = 1, fs_jpim1   ! vector opt. 
    747                iku    = mbku(ji  ,jj) 
    748                ikum1  = iku - 1 
    749                pe3_out(ji,jj,iku) = ( umask(ji,jj,iku) * (1.0_wp - zlnwd) + zlnwd )                              &  
    750                   &      *(  0.5_wp * r1_e1e2u(ji,jj)                                                            & 
    751                   &      *(  e1e2t(ji  ,jj) * ( SUM(tmask(ji  ,jj,:)*(pe3_in(ji  ,jj,:) - e3t_0(ji  ,jj,:))) )   & 
    752                   &        + e1e2t(ji+1,jj) * ( SUM(tmask(ji+1,jj,:)*(pe3_in(ji+1,jj,:) - e3t_0(ji+1,jj,:))) ) ) & 
    753                   &     - SUM(pe3_out(ji,jj,1:ikum1))) 
    754             END DO 
    755          END DO 
     769                  iku    = mbku(ji  ,jj) 
     770                  ikum1  = iku - 1 
     771                  pe3_out(ji,jj,iku) = ( umask(ji,jj,iku) * (1.0_wp - zlnwd) + zlnwd )                              &  
     772                     &      *(  0.5_wp * r1_e1e2u(ji,jj)                                                            & 
     773                     &      *(  e1e2t(ji  ,jj) * ( SUM(tmask(ji  ,jj,:)*(pe3_in(ji  ,jj,:) - e3t_0(ji  ,jj,:))) )   & 
     774                     &        + e1e2t(ji+1,jj) * ( SUM(tmask(ji+1,jj,:)*(pe3_in(ji+1,jj,:) - e3t_0(ji+1,jj,:))) ) ) & 
     775                     &     - SUM(pe3_out(ji,jj,1:ikum1))) 
     776               END DO 
     777            END DO 
     778            ! 
     779         CASE ( 2 )  
     780            zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3) 
     781            DO jk = 1, jpk 
     782               DO jj = 1, jpjm1 
     783                  DO ji = 1, fs_jpim1   ! vector opt. 
     784                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)    & 
     785                        &                       * (   e1e2t(ji  ,jj) * zssh(ji  ,jj) + e1e2t(ji+1,jj) * zssh(ji+1,jj)) & 
     786                        &                       * e3u_0(ji,jj,jk) / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) ) 
     787                  END DO 
     788               END DO 
     789            END DO 
     790            !    
     791         END SELECT 
    756792         ! 
    757793         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 
     
    759795         ! 
    760796      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    761          DO jk = 1, jpk 
     797         SELECT CASE ( nn_vvl_interp ) 
     798         CASE ( 0 )    
     799            !    
     800            DO jk = 1, jpk 
     801               DO jj = 1, jpjm1 
     802                  DO ji = 1, fs_jpim1   ! vector opt. 
     803                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
     804                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     805                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     806                  END DO 
     807               END DO 
     808            END DO 
     809            !                     ! 
     810         CASE ( 1 )  
     811            !    
     812            DO jk = 1, jpk 
     813               DO jj = 1, jpjm1 
     814                  DO ji = 1, fs_jpim1   ! vector opt. 
     815                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
     816                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     817                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     818                  END DO 
     819               END DO 
     820            END DO 
     821            ! 
     822            ! Bottom correction: 
    762823            DO jj = 1, jpjm1 
    763824               DO ji = 1, fs_jpim1   ! vector opt. 
    764                   pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    765                      &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    766                      &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    767                END DO 
    768             END DO 
    769          END DO 
    770          ! 
    771          ! Bottom correction: 
    772          DO jj = 1, jpjm1 
    773             DO ji = 1, fs_jpim1   ! vector opt. 
    774                ikv    = mbkv(ji  ,jj) 
    775                ikvm1  = ikv - 1 
    776                pe3_out(ji,jj,ikv) = ( vmask(ji,jj,ikv) * (1.0_wp - zlnwd) + zlnwd )                              &  
    777                   &      *(  0.5_wp * r1_e1e2v(ji,jj)                                                            & 
    778                   &      *(  e1e2t(ji,jj  ) * ( SUM(tmask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3t_0(ji,jj  ,:))) )   & 
    779                   &        + e1e2t(ji,jj+1) * ( SUM(tmask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3t_0(ji,jj+1,:))) ) ) & 
    780                   &     - SUM(pe3_out(ji,jj,1:ikvm1))) 
    781             END DO 
    782          END DO 
     825                  ikv    = mbkv(ji  ,jj) 
     826                  ikvm1  = ikv - 1 
     827                  pe3_out(ji,jj,ikv) = ( vmask(ji,jj,ikv) * (1.0_wp - zlnwd) + zlnwd )                              &  
     828                     &      *(  0.5_wp * r1_e1e2v(ji,jj)                                                            & 
     829                     &      *(  e1e2t(ji,jj  ) * ( SUM(tmask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3t_0(ji,jj  ,:))) )   & 
     830                     &        + e1e2t(ji,jj+1) * ( SUM(tmask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3t_0(ji,jj+1,:))) ) ) & 
     831                     &     - SUM(pe3_out(ji,jj,1:ikvm1))) 
     832               END DO 
     833            END DO 
     834            ! 
     835         CASE ( 2 )  
     836            zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3) 
     837            DO jk = 1, jpk 
     838               DO jj = 1, jpjm1 
     839                  DO ji = 1, fs_jpim1   ! vector opt. 
     840                     pe3_out(ji,jj,jk) = 0.5_wp * (  vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)    & 
     841                        &                       * (   e1e2t(ji  ,jj) * zssh(ji  ,jj) + e1e2t(ji,jj+1) * zssh(ji,jj+1)) & 
     842                        &                       * e3v_0(ji,jj,jk) / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) ) 
     843                  END DO 
     844               END DO 
     845            END DO 
     846            !    
     847         END SELECT 
    783848         ! 
    784849         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 
     
    786851         ! 
    787852      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    788          DO jk = 1, jpk 
     853         SELECT CASE ( nn_vvl_interp ) 
     854         CASE ( 0 )   
     855            DO jk = 1, jpk 
     856               DO jj = 1, jpjm1 
     857                  DO ji = 1, fs_jpim1   ! vector opt. 
     858                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     859                        &                       *    r1_e1e2f(ji,jj)                                                  & 
     860                        &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
     861                        &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     862                  END DO 
     863               END DO 
     864            END DO 
     865            ! 
     866         CASE ( 1 )  
     867            ! 
     868            DO jk = 1, jpk 
     869               DO jj = 1, jpjm1 
     870                  DO ji = 1, fs_jpim1   ! vector opt. 
     871                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     872                        &                       *    r1_e1e2f(ji,jj)                                                  & 
     873                        &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
     874                        &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     875                  END DO 
     876               END DO 
     877            END DO 
     878            ! 
     879            ! Bottom correction: 
    789880            DO jj = 1, jpjm1 
    790881               DO ji = 1, fs_jpim1   ! vector opt. 
    791                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
    792                      &                       *    r1_e1e2f(ji,jj)                                                  & 
    793                      &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    794                      &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    795                END DO 
    796             END DO 
    797          END DO 
    798          ! 
    799          ! Bottom correction: 
    800          DO jj = 1, jpjm1 
    801             DO ji = 1, fs_jpim1   ! vector opt. 
    802                ikf    = MIN(mbku(ji  ,jj),mbku(ji,jj+1)) 
    803                ikfm1  = ikf - 1 
    804                pe3_out(ji,jj,ikf) = ( umask(ji,jj,ikf) * umask(ji,jj+1,ikf) * (1.0_wp - zlnwd) + zlnwd )         & 
    805                   &     * ( 0.5_wp *  r1_e1e2f(ji,jj)                                                            & 
    806                   &     * (  e1e2u(ji,jj  ) * ( SUM(umask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3u_0(ji,jj  ,:))) )   & 
    807                   &        + e1e2u(ji,jj+1) * ( SUM(umask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3u_0(ji,jj+1,:))) ) ) & 
    808                   &     - SUM(pe3_out(ji,jj,1:ikfm1))) 
    809             END DO 
    810          END DO 
     882                  ikf    = MIN(mbku(ji  ,jj),mbku(ji,jj+1)) 
     883                  ikfm1  = ikf - 1 
     884                  pe3_out(ji,jj,ikf) = ( umask(ji,jj,ikf) * umask(ji,jj+1,ikf) * (1.0_wp - zlnwd) + zlnwd )         & 
     885                     &     * ( 0.5_wp *  r1_e1e2f(ji,jj)                                                            & 
     886                     &     * (  e1e2u(ji,jj  ) * ( SUM(umask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3u_0(ji,jj  ,:))) )   & 
     887                     &        + e1e2u(ji,jj+1) * ( SUM(umask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3u_0(ji,jj+1,:))) ) ) & 
     888                     &     - SUM(pe3_out(ji,jj,1:ikfm1))) 
     889               END DO 
     890            END DO 
     891            ! 
     892         CASE ( 2 )  
     893            zssh(:,:) = SUM(umask(:,:,:)*(pe3_in(:,:,:)-e3u_0(:,:,:)), DIM=3) 
     894            DO jk = 1, jpk 
     895               DO jj = 1, jpjm1 
     896                  DO ji = 1, fs_jpim1   ! vector opt. 
     897                     pe3_out(ji,jj,jk) =  (  umask(ji,jj,jk)* umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd )   & 
     898                        &                 * 0.5_wp * r1_e1e2f(ji,jj)                                           & 
     899                        &                 * (e1e2u(ji  ,jj) * zssh(ji  ,jj) + e1e2u(ji,jj+1) * zssh(ji,jj+1))  & 
     900                        &                 * e3f_0(ji,jj,jk) / ( hf_0(ji,jj) + 1._wp - ssumask(ji,jj)*ssumask(ji,jj+1) ) 
     901                  END DO 
     902               END DO 
     903            END DO       
     904            ! 
     905         END SELECT 
    811906         ! 
    812907         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 
     
    9971092               ! 
    9981093            ELSE 
    999                ! 
     1094              ! 
    10001095               ! usr_def_istate called here only to get sshb, that is needed to 
    10011096               ! initialize e3t_b and e3t_n  
     
    10611156      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & 
    10621157         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , & 
    1063          &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe 
     1158         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg            , &! not yet implemented: ln_vvl_kepe 
     1159         &              nn_vvl_interp 
    10641160      !!----------------------------------------------------------------------  
    10651161      ! 
     
    10971193         ENDIF 
    10981194         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg 
     1195         WRITE(numout,*) '         Method to compute scale factors anomaly at U/V/F points  nn_vvl_interp   = ', nn_vvl_interp 
    10991196      ENDIF 
    11001197      ! 
     
    11081205      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 
    11091206      ! 
     1207      IF( .NOT. ln_vvl_zstar .AND. (nn_vvl_interp==2 ) )  CALL ctl_stop( 'nn_vvl_interp must be < 2 if ln_vvl_zstar=F' ) 
    11101208      IF(lwp) THEN                   ! Print the choice 
    11111209         WRITE(numout,*) 
Note: See TracChangeset for help on using the changeset viewer.