Changeset 11012


Ignore:
Timestamp:
2019-05-20T18:18:05+02:00 (11 months ago)
Author:
jchanut
Message:

#1791, report domvvl change in VORTEX test case. Modification of domvvl there should be removed.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/fix_vvl_ticket1791/tests/VORTEX/MY_SRC/domvvl.F90

    r10572 r11012  
    139139      CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V  
    140140      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
    141       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F 
     141      CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F 
    142142      !                                ! Vertical interpolation of e3t,u,v  
    143143      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W 
     
    629629      ! - JC - hu_b, hv_b, hur_b, hvr_b also 
    630630       
    631       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  ) 
     631      CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F'  ) 
    632632       
    633633      ! Vertical scale factor interpolations 
     
    691691      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
    692692      ! 
    693       INTEGER ::   ji, jj, jk                                       ! dummy loop indices 
     693      INTEGER ::   ji, jj, jk, jkbot                                ! dummy loop indices 
     694      INTEGER ::   nmet                                             ! horizontal interpolation method 
    694695      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F 
    695       !!---------------------------------------------------------------------- 
     696      REAL(wp) ::  zmin, zsmall, zfractap                           ! Minimum thicknesses UVF-points 
     697      REAL(wp) ::  zdo, zup                                         ! Lower and upper interfaces depths anomalies 
     698      REAL(wp), DIMENSION(jpi,jpj) :: zs                            ! Surface interface depth anomaly 
     699      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw                        ! Interface depth anomaly 
     700      !!---------------------------------------------------------------------- 
     701      ! 
     702!      nmet = 0  ! Original method (Surely wrong) 
     703      nmet = 1  ! Interface interpolation 
     704!      nmet = 2  ! Internal interfaces interpolation only, spread barotropic increment 
     705      ! 
     706      zsmall = 1.e-10    ! Minimum thickness at U or V points (m) 
     707      zfractap = 0.8_wp  ! Fraction of T-point thickness in the shallowest direction 
    696708      ! 
    697709      IF(ln_wd_il) THEN 
     
    701713      END IF 
    702714      ! 
     715      IF ( ( nmet==1 ).OR.( nmet==2 ) ) THEN 
     716         SELECT CASE ( pout ) 
     717            ! 
     718         CASE( 'U', 'V', 'F' ) 
     719            ! 
     720            ! Retrieve ssh: 
     721            zs(:,:) = 0._wp 
     722            DO jk = 1, jpkm1 
     723                zs(:,:) = zs(:,:) + pe3_in(:,:,jk)*tmask(:,:,jk) 
     724            END DO  
     725            zs(:,:) = zs(:,:) - ht_0(:,:) 
     726            ! 
     727            ! Interface depth: 
     728            zw(:,:,:) =  0._wp 
     729            DO jk=2,jpk 
     730               zw(:,:,jk) = zw(:,:,jk-1) + pe3_in(:,:,jk-1)*tmask(:,:,jk-1) 
     731            END DO  
     732            ! 
     733            DO jk=1,jpkm1 
     734               zw(:,:,jk) = zw(:,:,jk) - zs(:,:) *tmask(:,:,jk) 
     735            END DO 
     736            zw(:,:,jpk) = gdepw_0(:,:,jpk) 
     737            ! 
     738            IF ( nmet==2 ) THEN        ! Consider "internal" interfaces only 
     739               ! 
     740               DO jj = 1, jpj 
     741                  DO ji = 1, jpi 
     742                     DO jk=1,jpk 
     743                        zw(ji,jj,jk) = ((zw(ji,jj,jk) + zs(ji,jj)-gdepw_0(ji,jj,mikt(ji,jj)))              & 
     744                                     & * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - ssmask(ji,jj))  & 
     745                                     &  + gdepw_0(ji,jj,mikt(ji,jj)) )* tmask(ji,jj,jk) 
     746                     END DO 
     747                  END DO 
     748               END DO  
     749            ENDIF  
     750            ! 
     751         END SELECT 
     752      END IF 
     753      ! 
     754      pe3_out(:,:,:) = 0.0_wp 
     755      ! 
    703756      SELECT CASE ( pout )    !==  type of interpolation  ==! 
    704757         ! 
    705758      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    706          DO jk = 1, jpk 
     759         IF (nmet==0) THEN 
     760            DO jk = 1, jpk 
     761               DO jj = 1, jpjm1 
     762                  DO ji = 1, fs_jpim1   ! vector opt. 
     763                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
     764                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     765                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     766                  END DO 
     767               END DO 
     768            END DO 
     769         ELSE 
    707770            DO jj = 1, jpjm1 
    708771               DO ji = 1, fs_jpim1   ! vector opt. 
    709                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
    710                      &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    711                      &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    712                END DO 
    713             END DO 
    714          END DO 
     772                  jkbot = mbku(ji,jj) 
     773                  zdo = hu_0(ji,jj) 
     774                  DO jk=jkbot,1,-1 
     775                     zup =   0.5_wp * ( e1e2t(ji  ,jj)*(zw(ji  ,jj,jk)-gdepw_0(ji  ,jj,jk))                     & 
     776                         &            + e1e2t(ji+1,jj)*(zw(ji+1,jj,jk)-gdepw_0(ji+1,jj,jk)) ) * r1_e1e2u(ji,jj) & 
     777                         & + 0.5_wp * (gdepw_0(ji  ,jj,jk)+gdepw_0(ji+1,jj,jk)) 
     778                     ! 
     779                     ! If there is a step, taper bottom interface: 
     780                     zmin = 0._wp 
     781                     IF (zup > zdo) THEN                      
     782                        IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN 
     783                           zmin = zfractap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk)) 
     784                        ELSE 
     785                           zmin = zfractap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk)) 
     786                        ENDIF 
     787                     ENDIF        
     788                     zmin = MAX(zmin, zsmall)         
     789                     zup  = MIN(zup, zdo-zmin) 
     790                     pe3_out(ji,jj,jk) = (zdo - zup) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) &  
     791                                       &               - umask(ji,jj,jk) * e3u_0(ji,jj,jk)  
     792                     zdo = zup 
     793                  END DO 
     794               END DO 
     795            END DO 
     796         END IF 
     797         !  
     798         IF (nmet==2) THEN           ! Spread sea level anomaly 
     799            DO jj = 1, jpjm1 
     800               DO ji = 1, fs_jpim1   ! vector opt. 
     801                  DO jk=1,jpk 
     802                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                   & 
     803                           &               + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) )               &  
     804                           &               / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )              & 
     805                           &               * 0.5_wp * r1_e1e2u(ji,jj)                              & 
     806                           &               * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )        & 
     807                           &               * ( e1e2t(ji,jj)*zs(ji,jj) + e1e2t(ji+1,jj)*zs(ji+1,jj) )        
     808                  END DO 
     809               END DO 
     810            END DO 
     811            ! 
     812         ENDIF 
     813         ! 
    715814         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 
    716815         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
    717816         ! 
    718817      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    719          DO jk = 1, jpk 
     818         IF (nmet==0) THEN 
     819            DO jk = 1, jpk 
     820               DO jj = 1, jpjm1 
     821                  DO ji = 1, fs_jpim1   ! vector opt. 
     822                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
     823                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     824                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     825                  END DO 
     826               END DO 
     827            END DO 
     828         ELSE 
    720829            DO jj = 1, jpjm1 
    721830               DO ji = 1, fs_jpim1   ! vector opt. 
    722                   pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    723                      &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    724                      &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    725                END DO 
    726             END DO 
    727          END DO 
     831                  jkbot = mbkv(ji,jj) 
     832                  zdo = hv_0(ji,jj) 
     833                  DO jk=jkbot,1,-1 
     834                     zup =   0.5_wp * ( e1e2t(ji  ,jj)*(zw(ji,jj  ,jk)-gdepw_0(ji,jj  ,jk))                     & 
     835                         &            + e1e2t(ji+1,jj)*(zw(ji,jj+1,jk)-gdepw_0(ji,jj+1,jk)) ) * r1_e1e2v(ji,jj) & 
     836                         & + 0.5_wp * (gdepw_0(ji,jj  ,jk)+gdepw_0(ji,jj+1,jk)) 
     837                     ! 
     838                     ! If there is a step, taper bottom interface: 
     839                     zmin = 0._wp 
     840                     IF (zup > zdo) THEN                      
     841                        IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN 
     842                           zmin = zfractap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk)) 
     843                        ELSE 
     844                           zmin = zfractap * (zw(ji,jj  ,jk+1)-zw(ji,jj  ,jk)) 
     845                        ENDIF 
     846                     ENDIF  
     847                     zmin = MAX(zmin, zsmall)         
     848                     zup  = MIN(zup, zdo-zmin) 
     849                     pe3_out(ji,jj,jk) = (zdo - zup) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) &  
     850                                       &               - vmask(ji,jj,jk) * e3v_0(ji,jj,jk)  
     851                     zdo = zup 
     852                  END DO 
     853               END DO 
     854            END DO 
     855         END IF 
     856         ! 
     857         IF (nmet==2) THEN           ! Spread sea level anomaly 
     858            DO jj = 1, jpjm1 
     859               DO ji = 1, fs_jpim1   ! vector opt. 
     860                  DO jk=1,jpk 
     861                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                                       & 
     862                           &               + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) )                                   &  
     863                           &               / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )                                  & 
     864                           &               * 0.5_wp * r1_e1e2v(ji,jj)                                                  & 
     865                                           * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )                            & 
     866                           &               * ( e1e2t(ji,jj)*zs(ji,jj) + e1e2t(ji,jj+1)*zs(ji,jj+1) )        
     867                  END DO 
     868               END DO 
     869            END DO 
     870            ! 
     871         ENDIF 
     872         ! 
    728873         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 
    729874         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
    730875         ! 
    731       CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    732          DO jk = 1, jpk 
     876      CASE( 'F' )                   !* from T-point to F-point : hor. surface weighted mean 
     877         IF (nmet==0) THEN 
     878            DO jk=1,jpk 
     879               DO jj = 1, jpjm1 
     880                  DO ji = 1, fs_jpim1   ! vector opt. 
     881                     pe3_out(ji,jj,jk) = 0.25_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     882                           &                     *    r1_e1e2f(ji,jj)                                                  & 
     883                           &                     * (  e1e2t(ji  ,jj  ) * ( pe3_in(ji  ,jj  ,jk)-e3t_0(ji  ,jj  ,jk) )  &  
     884                           &                        + e1e2t(ji  ,jj+1) * ( pe3_in(ji  ,jj+1,jk)-e3t_0(ji  ,jj+1,jk) )  & 
     885                           &                        + e1e2t(ji+1,jj  ) * ( pe3_in(ji+1,jj  ,jk)-e3t_0(ji+1,jj  ,jk) )  &  
     886                           &                        + e1e2t(ji+1,jj+1) * ( pe3_in(ji+1,jj+1,jk)-e3t_0(ji+1,jj+1,jk) ) )                                                   
     887                  END DO 
     888               END DO 
     889            END DO 
     890         ELSE 
    733891            DO jj = 1, jpjm1 
    734892               DO ji = 1, fs_jpim1   ! vector opt. 
    735                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
    736                      &                       *    r1_e1e2f(ji,jj)                                                  & 
    737                      &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    738                      &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    739                END DO 
    740             END DO 
    741          END DO 
     893                  jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1))  
     894                  zdo = hf_0(ji,jj) 
     895                  DO jk=jkbot,1,-1 
     896                     zup =     0.25_wp * (  e1e2t(ji  ,jj  ) * (zw(ji  ,jj  ,jk)-gdepw_0(ji  ,jj  ,jk))                         &  
     897                           &              + e1e2t(ji+1,jj  ) * (zw(ji+1,jj  ,jk)-gdepw_0(ji+1,jj  ,jk))                         &  
     898                           &              + e1e2t(ji  ,jj+1) * (zw(ji  ,jj+1,jk)-gdepw_0(ji  ,jj+1,jk))                         &  
     899                           &              + e1e2t(ji+1,jj+1) * (zw(ji+1,jj+1,jk)-gdepw_0(ji+1,jj+1,jk))  ) *    r1_e1e2f(ji,jj) & 
     900                           & + 0.25_wp * ( gdepw_0(ji  ,jj  ,jk) + gdepw_0(ji+1,jj  ,jk) & 
     901                           &              +gdepw_0(ji  ,jj+1,jk) + gdepw_0(ji+1,jj+1,jk) ) 
     902                     ! 
     903                     ! If there is a step, taper bottom interface: 
     904                     zmin = 0._wp 
     905                     IF (zup > zdo) THEN  
     906                        IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN 
     907                           IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN 
     908                              zmin = zfractap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk)) 
     909                           ELSE 
     910                              zmin = zfractap * (zw(ji  ,jj+1,jk+1)-zw(ji  ,jj+1,jk)) 
     911                           ENDIF 
     912                        ELSE 
     913                           IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN 
     914                              zmin = zfractap * (zw(ji+1,jj  ,jk+1)-zw(ji+1,jj  ,jk)) 
     915                           ELSE 
     916                              zmin = zfractap * (zw(ji  ,jj  ,jk+1)-zw(ji  ,jj  ,jk)) 
     917                           ENDIF 
     918                        ENDIF 
     919                     ENDIF 
     920                     zmin = MAX(zmin, zsmall)         
     921                     zup  = MIN(zup, zdo-zmin) 
     922                     ! 
     923                     pe3_out(ji,jj,jk) = ( zdo - zup ) &  
     924                                      & *( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     925                                      &  - umask(ji,jj,jk) * umask(ji,jj+1,jk) * e3f_0(ji,jj,jk)  
     926                     zdo = zup 
     927                  END DO 
     928               END DO 
     929            END DO 
     930         END IF 
     931         ! 
     932         IF (nmet==2) THEN           ! Spread sea level anomaly 
     933            ! 
     934            DO jj = 1, jpjm1 
     935               DO ji = 1, fs_jpim1   ! vector opt. 
     936                  DO jk=1,jpk 
     937                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                           & 
     938                           &               + ( pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) )                       &  
     939                           &               / ( hf_0(ji,jj) + 1._wp - ssumask(ji,jj)*ssumask(ji,jj+1) )     & 
     940                           &               * 0.25_wp * r1_e1e2f(ji,jj)                                     &  
     941                           &               * ( umask(ji,jj,jk)*umask(ji,jj+1,jk)*(1.0_wp - zlnwd) + zlnwd )& 
     942                           &               * ( e1e2t(ji  ,jj)*zs(ji  ,jj) + e1e2t(ji  ,jj+1)*zs(ji  ,jj+1) & 
     943                           &                  +e1e2t(ji+1,jj)*zs(ji+1,jj) + e1e2t(ji+1,jj+1)*zs(ji+1,jj+1) )                
     944                  END DO 
     945               END DO 
     946            END DO 
     947         END IF 
     948         ! 
    742949         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 
    743950         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
     
    9291136               ! 
    9301137               ! usr_def_istate called here only to get sshb, that is needed to initialize e3t_b and e3t_n 
    931                CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )   
     1138               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  ) 
    9321139               ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 
    9331140               ! 
Note: See TracChangeset for help on using the changeset viewer.