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

Changeset 10281


Ignore:
Timestamp:
2018-11-07T12:22:39+01:00 (5 years ago)
Author:
clem
Message:

advection routine containing all the tests implemented (Zalesak, Devore, etc)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn_adv_umx.F90

    r10277 r10281  
    3535   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    3636 
     37   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_a 
     38    
     39   ! alternate directions for upstream 
     40   ! clem: it gives results for Lipscomb test that are the same as "ll_upsxy=false" 
     41   LOGICAL :: ll_upsxy = .FALSE. 
     42    
     43   ! prelimiter 
     44   LOGICAL :: ll_prelimiter = .TRUE. 
     45   LOGICAL :: ll_prelimiter_zalesak = .TRUE.   ! from: Zalesak(1979) eq. 14 => better than Devore 
     46   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11 
     47 
     48   ! iterate on the limiter (only nonosc) 
     49   LOGICAL :: ll_limiter_it2 = .TRUE. 
     50    
     51 
    3752   !! * Substitutions 
    3853#  include "vectopt_loop_substitute.h90" 
     
    7792      !                                     !    must be >= 0.01 and the best seems to be 0.1 
    7893      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zfu, zfv 
    79       REAL(wp), DIMENSION(jpi,jpj) ::   z1_a, zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, z1_ap, zh_ip  
     94      REAL(wp), DIMENSION(jpi,jpj) ::   zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, zh_ip  
    8095      !!---------------------------------------------------------------------- 
    8196      ! 
     
    112127      END DO 
    113128 
     129      IF(.NOT. ALLOCATED(z1_a )) ALLOCATE(z1_a (jpi,jpj)) 
    114130      !---------------! 
    115131      !== advection ==! 
     
    149165            END DO 
    150166            ! 
    151             IF ( ln_pnd_H12 ) THEN   ! melt ponds 
     167            IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_a...) 
    152168               ! to avoid a problem with the limiter nonosc when A gets close to 0 
    153169               pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1) 
    154170               ! 
    155                WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_ap(:,:) = 1._wp / pa_ip(:,:,jl) 
    156                ELSEWHERE                         ;   z1_ap(:,:) = 0. 
     171               WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_ip(:,:,jl) 
     172               ELSEWHERE                         ;   z1_a(:,:) = 0. 
    157173               END WHERE 
    158174               ! 
     
    213229      INTEGER  ::   ji, jj           ! dummy loop indices   
    214230      REAL(wp) ::   ztra             ! local scalar 
    215 !!clem      REAL(wp) ::   zeps = 1.e-02        ! local scalar 
    216231      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 
    217       REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt 
     232      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt, zz1_a 
    218233      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc  
    219234      !!---------------------------------------------------------------------- 
    220235      ! 
    221       ! add a constant value to avoid problems with zeros 
    222       DO jj = 1, jpj 
    223          DO ji = 1, jpi 
    224             zpt(ji,jj) = pt(ji,jj) !!clem + zeps * tmask(ji,jj,1) 
    225          END DO 
    226       END DO 
    227  
    228236      !  upstream (_ups) advection with initial mass fluxes 
    229237      ! --------------------------------------------------- 
    230       ! fluxes 
    231       DO jj = 1, jpjm1 
    232          DO ji = 1, fs_jpim1 
    233             zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj) 
    234             zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1) 
    235          END DO 
    236       END DO 
     238      IF( .NOT. ll_upsxy ) THEN 
     239 
     240         ! fluxes 
     241         DO jj = 1, jpjm1 
     242            DO ji = 1, fs_jpim1 
     243               zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj) 
     244               zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1) 
     245            END DO 
     246         END DO 
     247 
     248      ELSE 
     249         ! 
     250         IF    ( pamsk == 1. ) THEN ; zz1_a(:,:) = 1._wp     ;        ! 1   if advection of A 
     251         ELSEIF( pamsk == 0. ) THEN ; zz1_a(:,:) = z1_a(:,:) ; ENDIF  ! 1/A if advection of V 
     252         ! 
     253         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     254            ! flux in x-direction 
     255            DO jj = 1, jpjm1 
     256               DO ji = 1, fs_jpim1 
     257                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj) 
     258               END DO 
     259            END DO 
     260            ! first guess of tracer content from u-flux 
     261            DO jj = 2, jpjm1 
     262               DO ji = fs_2, fs_jpim1   ! vector opt. 
     263                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * zz1_a(ji,jj) 
     264               END DO 
     265            END DO 
     266            CALL lbc_lnk( zpt, 'T', 1. ) 
     267            ! 
     268            ! flux in y-direction 
     269            DO jj = 1, jpjm1 
     270               DO ji = 1, fs_jpim1 
     271                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1) 
     272               END DO 
     273            END DO 
     274            ! 
     275         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     276            ! flux in y-direction 
     277            DO jj = 1, jpjm1 
     278               DO ji = 1, fs_jpim1 
     279                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1) 
     280               END DO 
     281            END DO 
     282            ! first guess of tracer content from v-flux 
     283            DO jj = 2, jpjm1 
     284               DO ji = fs_2, fs_jpim1   ! vector opt. 
     285                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * zz1_a(ji,jj)  
     286               END DO 
     287            END DO 
     288            CALL lbc_lnk( zpt, 'T', 1. ) 
     289            ! 
     290            ! flux in x-direction 
     291            DO jj = 1, jpjm1 
     292               DO ji = 1, fs_jpim1 
     293                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj) 
     294               END DO 
     295            END DO 
     296            ! 
     297         ENDIF 
     298          
     299      ENDIF 
    237300      ! guess after content field with upstream scheme 
    238301      DO jj = 2, jpjm1 
     
    251314      CASE ( 20 )                          !== centered second order ==! 
    252315         ! 
    253          CALL cen2( kn_limiter, jt, kt, pdt, zpt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  & 
     316         CALL cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  & 
    254317            &       zt_ups, zfu_ups, zfv_ups ) 
    255318         ! 
    256319      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==! 
    257320         ! 
    258          CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, zpt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  & 
     321         CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  & 
    259322            &        zt_ups, zfu_ups, zfv_ups ) 
    260323         ! 
     
    266329         DO jj = 1, jpjm1 
    267330            DO ji = 1, fs_jpim1 
    268                pfu(ji,jj) =  zfu_ho(ji,jj) !!clem - zeps * puc(ji,jj) 
    269                pfv(ji,jj) =  zfv_ho(ji,jj) !!clem - zeps * pvc(ji,jj) 
     331               pfu(ji,jj) =  zfu_ho(ji,jj) 
     332               pfv(ji,jj) =  zfv_ho(ji,jj) 
    270333            END DO 
    271334         END DO 
     
    278341      DO jj = 2, jpjm1 
    279342         DO ji = fs_2, fs_jpim1  
    280             ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )         &  !        Div(uaH) or        Div(ua) 
    281 !!clem               &           + ( puc   (ji,jj) - puc   (ji-1,jj) + pvc   (ji,jj) - pvc   (ji,jj-1) ) * zeps  &  ! epsi * Div(ua)  or epsi * Div(u)  
     343            ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua) 
    282344               &         ) * r1_e1e2t(ji,jj)   
    283345            ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     
    430492      ! 
    431493      INTEGER  ::   ji, jj    ! dummy loop indices 
    432       REAL(wp), DIMENSION(jpi,jpj) ::   zzt 
     494      REAL(wp) ::   ztra 
     495      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zzfu_ho, zzfv_ho 
    433496      !!---------------------------------------------------------------------- 
    434497      ! 
     
    478541         ! 
    479542      ENDIF 
    480       IF( kn_limiter == 1 )   CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     543      IF( kn_limiter == 1 ) THEN 
     544         IF( .NOT. ll_limiter_it2 ) THEN 
     545            CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     546         ELSE 
     547            zzfu_ho(:,:) = pfu_ho(:,:) 
     548            zzfv_ho(:,:) = pfv_ho(:,:) 
     549            ! 1st iteration of nonosc (limit the flux with the upstream solution) 
     550            CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     551            ! guess after content field with high order 
     552            DO jj = 2, jpjm1 
     553               DO ji = fs_2, fs_jpim1 
     554                  ztra       = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     555                  zzt(ji,jj) =   ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     556               END DO 
     557            END DO 
     558            CALL lbc_lnk( zzt, 'T', 1. ) 
     559            ! 2nd iteration of nonosc (limit the flux with the limited high order solution) 
     560            CALL nonosc_2d ( pdt, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     561         ENDIF 
     562      ENDIF 
    481563      ! 
    482564   END SUBROUTINE macho 
     
    733815      
    734816 
    735    SUBROUTINE nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     817   SUBROUTINE nonosc_2d( pdt, ptc, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho ) 
    736818      !!--------------------------------------------------------------------- 
    737819      !!                    ***  ROUTINE nonosc  *** 
    738820      !!      
    739       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
     821       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
    740822      !!       scheme and the before field by a nonoscillatory algorithm  
    741823      !! 
    742824      !! **  Method  :   ... ??? 
    743       !!       warning : ptc and pt_ups must be masked, but the boundaries 
     825      !!       warning : ptc and pt_low must be masked, but the boundaries 
    744826      !!       conditions on the fluxes are not necessary zalezak (1979) 
    745827      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
     
    747829      !!---------------------------------------------------------------------- 
    748830      REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step 
    749       REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_ups      ! before field & upstream guess of after field 
    750       REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pfv_ups, pfu_ups ! upstream flux 
     831      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_low      ! before field & upstream guess of after field 
     832      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pfv_low, pfu_low ! upstream flux 
    751833      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux 
    752834      ! 
    753835      INTEGER  ::   ji, jj    ! dummy loop indices 
    754836      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt   ! local scalars 
    755       REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo         !   -      - 
    756       REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv 
     837      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign        !   -      - 
     838      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv, zti_low, ztj_low 
    757839      !!---------------------------------------------------------------------- 
    758840      zbig = 1.e+40_wp 
     
    771853      DO jj = 1, jpjm1 
    772854         DO ji = 1, fs_jpim1   ! vector opt. 
    773             pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_ups(ji,jj) 
    774             pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_ups(ji,jj) 
    775          END DO 
    776       END DO 
    777  
     855            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj) 
     856            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj) 
     857         END DO 
     858      END DO 
     859 
     860      ! extreme case where pfu_ho has to be zero 
     861      ! ---------------------------------------- 
     862      !                                    pfu_ho 
     863      !                           *         ---> 
     864      !                        |      |  *   |        |  
     865      !                        |      |      |    *   |     
     866      !                        |      |      |        |    * 
     867      !            t_low :       i-1     i       i+1       i+2    
     868      IF( ll_prelimiter ) THEN 
     869 
     870         DO jj = 2, jpjm1 
     871            DO ji = fs_2, fs_jpim1  
     872               zti_low(ji,jj)= pt_low(ji+1,jj  ) 
     873               ztj_low(ji,jj)= pt_low(ji  ,jj+1) 
     874            END DO 
     875         END DO 
     876         CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. ) 
     877 
     878         IF( ll_prelimiter_zalesak ) THEN 
     879            DO jj = 2, jpjm1 
     880               DO ji = fs_2, fs_jpim1 
     881                  IF(    SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj) - pt_low (ji  ,jj) ) ) THEN 
     882                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj) - zti_low(ji  ,jj) ) )   pfu_ho(ji,jj) = 0. 
     883                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj) ) )   pfu_ho(ji,jj) = 0. 
     884                  ENDIF 
     885                  IF(    SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji,jj+1) - pt_low (ji,jj  ) ) ) THEN 
     886                     IF( SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) )   pfv_ho(ji,jj) = 0. 
     887                     IF( SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji,jj  ) - pt_low (ji,jj-1) ) )   pfv_ho(ji,jj) = 0. 
     888                  ENDIF 
     889               END DO 
     890            END DO 
     891            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
     892 
     893         ELSEIF( ll_prelimiter_devore ) THEN 
     894            z1_dt = 1._wp / pdt 
     895            DO jj = 2, jpjm1 
     896               DO ji = fs_2, fs_jpim1 
     897                  zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) ) 
     898                  pfu_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , & 
     899                     &                          zsign * ( pt_low (ji  ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji  ,jj) * z1_dt , & 
     900                     &                          zsign * ( zti_low(ji+1,jj) - zti_low(ji  ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) ) 
     901                   
     902                  zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) ) 
     903                  pfv_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , & 
     904                     &                          zsign * ( pt_low (ji,jj  ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj  ) * z1_dt , & 
     905                     &                          zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) * e1e2t(ji,jj+1) * z1_dt ) ) 
     906               END DO 
     907            END DO 
     908            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
     909             
     910         ENDIF 
     911          
     912      ENDIF 
     913       
    778914      ! Search local extrema 
    779915      ! -------------------- 
    780       ! max/min of ptc & pt_ups with large negative/positive value (-/+zbig) outside ice cover 
     916      ! max/min of ptc & pt_low with large negative/positive value (-/+zbig) outside ice cover 
    781917      DO jj = 1, jpj 
    782          DO ji = fs_2, fs_jpim1 
    783             IF( ptc(ji,jj) == 0._wp .AND. pt_ups(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN 
     918         DO ji = 1, jpi 
     919            IF( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN 
    784920               zbup(ji,jj) = -zbig 
    785921               zbdo(ji,jj) =  zbig 
    786922            ELSE                
    787                zbup(ji,jj) = MAX( ptc(ji,jj) , pt_ups(ji,jj) ) 
    788                zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_ups(ji,jj) ) 
     923               zbup(ji,jj) = MAX( ptc(ji,jj) , pt_low(ji,jj) ) 
     924               zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_low(ji,jj) ) 
    789925            ENDIF 
    790926         END DO 
     
    795931         DO ji = fs_2, fs_jpim1   ! vector opt. 
    796932            ! 
    797             zup  = MAX( zbup(ji,jj), zbup(ji-1,jj), zbup(ji+1,jj), zbup(ji,jj-1), zbup(ji,jj+1) )         ! search max/min in neighbourhood 
    798             zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) ) 
    799             ! 
    800             zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) & 
    801                & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) )  ! positive/negative part of the flux 
     933            !zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1) )  ! search max/min in neighbourhood 
     934            !zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) ) 
     935            ! 
     936            zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1), &  ! search max/min in neighbourhood 
     937               &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  ) 
     938            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), & 
     939               &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  ) 
     940            ! 
     941            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &  ! positive/negative part of the flux 
     942               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) ) 
    802943            zneg = MAX( 0., pfu_ho(ji  ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) & 
    803944               & + MAX( 0., pfv_ho(ji,jj  ) ) - MIN( 0., pfv_ho(ji,jj-1) ) 
    804945            ! 
    805946            !                                  ! up & down beta terms 
    806 !!clem            zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt 
    807 !!clem            zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt 
    808             IF( zpos >= epsi20 ) THEN 
    809                zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt 
    810             ELSE 
    811                zbetup(ji,jj) = zbig 
     947!!clem            zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt 
     948!!clem            zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt 
     949            IF( zpos >= epsi20 ) THEN ; zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt 
     950            ELSE                      ; zbetup(ji,jj) = zbig 
    812951            ENDIF 
    813952            ! 
    814             IF( zneg >= epsi20 ) THEN 
    815                zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
    816             ELSE 
    817                zbetdo(ji,jj) = zbig 
     953            IF( zneg >= epsi20 ) THEN ; zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
     954            ELSE                      ; zbetdo(ji,jj) = zbig 
    818955            ENDIF 
    819956            ! 
     
    830967            zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) ) 
    831968            ! 
    832             pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_ups(ji,jj) 
     969            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_low(ji,jj) 
    833970         END DO 
    834971      END DO 
     
    840977            zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) ) 
    841978            ! 
    842             pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_ups(ji,jj) 
     979            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_low(ji,jj) 
    843980         END DO 
    844981      END DO 
Note: See TracChangeset for help on using the changeset viewer.