Ignore:
Timestamp:
2018-11-15T17:31:50+01:00 (2 years ago)
Author:
clem
Message:

still playing around with nonosc limiter to get a good ice thickness at the ice edge. There is still a problem when the ice concentration is very small

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90

    r10267 r10315  
    3535   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    3636 
     37   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_a, z1_a_ups, zua_ups, zva_ups 
     38    
     39   ! alternate directions for upstream 
     40   ! clem: it gives results for Lipscomb test that are the same as "ll_upsxy=false" 
     41   ! clem: needs to be set to true in 2D when using prelimiter (otherwise "wavy solutions" are created) 
     42   LOGICAL :: ll_upsxy = .TRUE. 
     43    
     44   ! prelimiter 
     45   ! clem: use it to avoid overshoot in H 
     46   LOGICAL :: ll_prelimiter = .TRUE. 
     47   LOGICAL :: ll_prelimiter_zalesak = .TRUE.  ! from: Zalesak(1979) eq. 14 => better than Devore especially in 2D 
     48   LOGICAL :: ll_prelimiter_devore  = .FALSE. ! from: Devore eq. 11 
     49 
     50   ! iterate on the limiter (only nonosc) 
     51   ! clem: useless 
     52   LOGICAL :: ll_limiter_it2 = .FALSE. 
     53    
     54 
    3755   !! * Substitutions 
    3856#  include "vectopt_loop_substitute.h90" 
     
    7492      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    7593      REAL(wp) ::   zcfl , zdt 
    76       REAL(wp) ::   zeps = 0.1_wp           ! shift in concentration to avoid division by 0 
     94      REAL(wp) ::   zeps = 0.0_wp           ! shift in concentration to avoid division by 0 
    7795      !                                     !    must be >= 0.01 and the best seems to be 0.1 
    78       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  
     96      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho, za_ups 
     97      REAL(wp), DIMENSION(jpi,jpj) ::   zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, zh_ip  
    8098      !!---------------------------------------------------------------------- 
    8199      ! 
     
    112130      END DO 
    113131 
     132      IF(.NOT. ALLOCATED(z1_a)    )   ALLOCATE(z1_a    (jpi,jpj)) 
     133      IF(.NOT. ALLOCATED(z1_a_ups))   ALLOCATE(z1_a_ups(jpi,jpj)) 
     134      IF(.NOT. ALLOCATED(zua_ups) )   ALLOCATE(zua_ups (jpi,jpj)) 
     135      IF(.NOT. ALLOCATED(zva_ups) )   ALLOCATE(zva_ups (jpi,jpj)) 
    114136      !---------------! 
    115137      !== advection ==! 
     
    117139      DO jt = 1, icycle 
    118140 
    119          zamsk = 1._wp 
    120          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area  
     141         zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:)  
     142         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:), zua_ups, zva_ups, za_ups, zua_ho, zva_ho )   ! Open water area  
    121143 
    122144         DO jl = 1, jpl 
     
    125147            ! 
    126148            WHERE( pa_i(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_i(:,:,jl) 
    127             ELSEWHERE                        ;   z1_a(:,:) = 0.   !!;   pa_i(:,:,jl) = 0. ; pv_i(:,:,jl) = 0.  
     149            ELSEWHERE                        ;   z1_a(:,:) = 0. 
    128150            END WHERE 
    129151            ! 
    130             zamsk = 1._wp 
    131             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zfu, zfv ) ! Ice area 
    132             !!zfu = zudy ; zfv = zvdx 
    133             !!CALL adv_umx( kn_umx, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, pv_i(:,:,jl), pv_i(:,:,jl) ) 
     152            zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:)  
     153            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zua_ups, zva_ups, za_ups, zua_ho, zva_ho ) ! Ice area 
     154 
     155            ! 1/A_ups 
     156            WHERE( za_ups(:,:) > epsi20 )   ;   z1_a_ups(:,:) = 1._wp / za_ups(:,:) 
     157            ELSEWHERE                       ;   z1_a_ups(:,:) = 0. 
     158            END WHERE 
     159 
     160            !!zua_ho = zudy ; zva_ho = zvdx 
     161            !!CALL adv_umx( kn_umx, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, pv_i(:,:,jl), pv_i(:,:,jl) ) 
    134162            zamsk = 0._wp ; zh_i(:,:) = pv_i (:,:,jl) * z1_a(:,:) 
    135             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl) )              ! Ice volume 
     163            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl), zua_ups, zva_ups )              ! Ice volume 
    136164            zamsk = 0._wp ; zh_s(:,:) = pv_s (:,:,jl) * z1_a(:,:) 
    137             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl) )              ! Snw volume 
     165            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl), zua_ups, zva_ups )              ! Snw volume 
    138166            zamsk = 0._wp ; zs_i(:,:) = psv_i(:,:,jl) * z1_a(:,:) 
    139             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl) )              ! Salt content 
     167            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl), zua_ups, zva_ups )              ! Salt content 
    140168            zamsk = 0._wp ; zo_i(:,:) = poa_i(:,:,jl) * z1_a(:,:) 
    141             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl) )              ! Age content 
     169            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl), zua_ups, zva_ups )              ! Age content 
    142170            DO jk = 1, nlay_i 
    143171               zamsk = 0._wp ; ze_i(:,:) = pe_i(:,:,jk,jl) * z1_a(:,:) 
    144                CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_i(:,:), pe_i(:,:,jk,jl) )           ! Ice heat content 
     172               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, ze_i(:,:), pe_i(:,:,jk,jl), zua_ups, zva_ups )           ! Ice heat content 
    145173            END DO 
    146174            DO jk = 1, nlay_s 
    147175               zamsk = 0._wp ; ze_s(:,:) = pe_s(:,:,jk,jl) * z1_a(:,:) 
    148                CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl) )           ! Snw heat content 
    149             END DO 
    150             ! 
    151             IF ( ln_pnd_H12 ) THEN   ! melt ponds 
     176               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl), zua_ups, zva_ups )           ! Snw heat content 
     177            END DO 
     178            ! 
     179            IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_a...) 
    152180               ! to avoid a problem with the limiter nonosc when A gets close to 0 
    153181               pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1) 
    154182               ! 
    155                WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_ap(:,:) = 1._wp / pa_ip(:,:,jl) 
    156                ELSEWHERE                         ;   z1_ap(:,:) = 0. 
     183               WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_ip(:,:,jl) 
     184               ELSEWHERE                         ;   z1_a(:,:) = 0. 
    157185               END WHERE 
    158186               ! 
    159                zamsk = 1._wp 
    160                CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zfu, zfv ) ! mp fraction 
     187               zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:)  
     188               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zua_ups, zva_ups, za_ups, zua_ho, zva_ho ) ! mp fraction 
     189 
     190               ! 1/A_ups 
     191               WHERE( za_ups(:,:) > epsi20 )   ;   z1_a_ups(:,:) = 1._wp / za_ups(:,:) 
     192               ELSEWHERE                       ;   z1_a_ups(:,:) = 0. 
     193               END WHERE 
     194 
    161195               zamsk = 0._wp ; zh_ip(:,:) = pv_ip(:,:,jl) * z1_a(:,:) 
    162                CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl) )              ! mp volume 
     196               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl), zua_ups, zva_ups )              ! mp volume 
    163197            ENDIF 
    164198            ! 
     
    177211              END DO 
    178212            END DO 
    179             CALL lbc_lnk( pa_i(:,:,jl), 'T',  1. ) 
     213            CALL lbc_lnk_multi( pa_i(:,:,jl), 'T',  1., pa_ip(:,:,jl), 'T',  1. ) 
    180214            ! 
    181215         END DO 
     
    186220 
    187221    
    188    SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pfu, pfv ) 
     222   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ups, pva_ups, pa_ups, pua_ho, pva_ho ) 
    189223      !!---------------------------------------------------------------------- 
    190224      !!                  ***  ROUTINE adv_umx  *** 
     
    209243      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pt             ! tracer field 
    210244      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc            ! tracer content field 
    211       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pfu, pfv       ! high order fluxes 
     245      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pua_ups, pva_ups ! upstream u*a fluxes or u 
     246      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pa_ups         ! concentration advected upstream 
     247      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes 
    212248      ! 
    213249      INTEGER  ::   ji, jj           ! dummy loop indices   
    214250      REAL(wp) ::   ztra             ! local scalar 
    215 !!clem      REAL(wp) ::   zeps = 1.e-02        ! local scalar 
    216251      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 
    217252      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt 
     
    219254      !!---------------------------------------------------------------------- 
    220255      ! 
    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  
    228256      !  upstream (_ups) advection with initial mass fluxes 
    229257      ! --------------------------------------------------- 
    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 
     258      IF( .NOT. ll_upsxy ) THEN 
     259 
     260         ! fluxes 
     261         DO jj = 1, jpjm1 
     262            DO ji = 1, fs_jpim1 
     263               zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * pt(ji+1,jj) 
     264               zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * pt(ji,jj+1) 
     265            END DO 
     266         END DO 
     267 
     268      ELSE 
     269         ! 1 if advection of A 
     270         ! z1_a_ups already defined IF advection of other tracers 
     271         IF( pamsk == 1. )   z1_a_ups(:,:) = 1._wp 
     272         ! 
     273         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     274            ! flux in x-direction 
     275            DO jj = 1, jpjm1 
     276               DO ji = 1, fs_jpim1 
     277                  zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * pt(ji+1,jj) 
     278               END DO 
     279            END DO 
     280            ! first guess of tracer content from u-flux 
     281            DO jj = 2, jpjm1 
     282               DO ji = fs_2, fs_jpim1   ! vector opt. 
     283                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) & 
     284                     &         * tmask(ji,jj,1) * z1_a_ups(ji,jj) 
     285               END DO 
     286            END DO 
     287            CALL lbc_lnk( zpt, 'T', 1. ) 
     288            ! 
     289            ! flux in y-direction 
     290            DO jj = 1, jpjm1 
     291               DO ji = 1, fs_jpim1 
     292                  zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * zpt(ji,jj+1) 
     293               END DO 
     294            END DO 
     295            ! 
     296         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     297            ! flux in y-direction 
     298            DO jj = 1, jpjm1 
     299               DO ji = 1, fs_jpim1 
     300                  zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * pt(ji,jj+1) 
     301               END DO 
     302            END DO 
     303            ! first guess of tracer content from v-flux 
     304            DO jj = 2, jpjm1 
     305               DO ji = fs_2, fs_jpim1   ! vector opt. 
     306                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 
     307                     &         * tmask(ji,jj,1) * z1_a_ups(ji,jj)  
     308               END DO 
     309            END DO 
     310            CALL lbc_lnk( zpt, 'T', 1. ) 
     311            ! 
     312            ! flux in x-direction 
     313            DO jj = 1, jpjm1 
     314               DO ji = 1, fs_jpim1 
     315                  zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * zpt(ji+1,jj) 
     316               END DO 
     317            END DO 
     318            ! 
     319         ENDIF 
     320          
     321      ENDIF 
    237322      ! guess after content field with upstream scheme 
    238323      DO jj = 2, jpjm1 
     
    244329      END DO 
    245330      CALL lbc_lnk( zt_ups, 'T', 1. ) 
    246        
     331 
    247332      ! High order (_ho) fluxes  
    248333      ! ----------------------- 
     
    251336      CASE ( 20 )                          !== centered second order ==! 
    252337         ! 
    253          CALL cen2( kn_limiter, jt, kt, pdt, zpt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  & 
     338         CALL cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  & 
    254339            &       zt_ups, zfu_ups, zfv_ups ) 
    255340         ! 
    256341      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==! 
    257342         ! 
    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,  & 
     343         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,  & 
    259344            &        zt_ups, zfu_ups, zfv_ups ) 
    260345         ! 
    261346      END SELECT 
    262347          
    263       ! output high order fluxes 
    264       ! ------------------------ 
    265       IF( PRESENT(pfu) ) THEN  
     348      ! output upstream trend of concentration and high order fluxes 
     349      ! ------------------------------------------------------------ 
     350      IF( pamsk == 1. ) THEN 
     351         ! upstream trend of concentration 
     352         pa_ups(:,:) = zt_ups(:,:) 
     353         ! upstream and high order u*a 
    266354         DO jj = 1, jpjm1 
    267355            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) 
    270             END DO 
    271          END DO 
    272          !!CALL lbc_lnk( pfu, 'U',  -1. ) ! clem: not needed I think 
    273          !!CALL lbc_lnk( pfv, 'V',  -1. ) 
     356               pua_ups(ji,jj) = zfu_ups(ji,jj) 
     357               pva_ups(ji,jj) = zfv_ups(ji,jj) 
     358               pua_ho (ji,jj) = zfu_ho (ji,jj) 
     359               pva_ho (ji,jj) = zfv_ho (ji,jj) 
     360            END DO 
     361         END DO 
     362         !!CALL lbc_lnk( pua_ho, 'U',  -1. ) ! clem: not needed I think 
     363         !!CALL lbc_lnk( pva_ho, 'V',  -1. ) 
    274364      ENDIF 
    275365       
     
    278368      DO jj = 2, jpjm1 
    279369         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)  
     370            ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua) 
    282371               &         ) * r1_e1e2t(ji,jj)   
    283372            ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     
    430519      ! 
    431520      INTEGER  ::   ji, jj    ! dummy loop indices 
    432       REAL(wp), DIMENSION(jpi,jpj) ::   zzt 
     521      REAL(wp) ::   ztra 
     522      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zzfu_ho, zzfv_ho 
    433523      !!---------------------------------------------------------------------- 
    434524      ! 
     
    478568         ! 
    479569      ENDIF 
    480       IF( kn_limiter == 1 )   CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     570      IF( kn_limiter == 1 ) THEN 
     571         IF( .NOT. ll_limiter_it2 ) THEN 
     572            CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     573         ELSE 
     574            zzfu_ho(:,:) = pfu_ho(:,:) 
     575            zzfv_ho(:,:) = pfv_ho(:,:) 
     576            ! 1st iteration of nonosc (limit the flux with the upstream solution) 
     577            CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     578            ! guess after content field with high order 
     579            DO jj = 2, jpjm1 
     580               DO ji = fs_2, fs_jpim1 
     581                  ztra       = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     582                  zzt(ji,jj) =   ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     583               END DO 
     584            END DO 
     585            CALL lbc_lnk( zzt, 'T', 1. ) 
     586            ! 2nd iteration of nonosc (limit the flux with the limited high order solution) 
     587            CALL nonosc_2d ( pdt, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     588         ENDIF 
     589      ENDIF 
    481590      ! 
    482591   END SUBROUTINE macho 
     
    733842      
    734843 
    735    SUBROUTINE nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     844   SUBROUTINE nonosc_2d( pdt, ptc, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho ) 
    736845      !!--------------------------------------------------------------------- 
    737846      !!                    ***  ROUTINE nonosc  *** 
    738847      !!      
    739       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
     848       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
    740849      !!       scheme and the before field by a nonoscillatory algorithm  
    741850      !! 
    742851      !! **  Method  :   ... ??? 
    743       !!       warning : ptc and pt_ups must be masked, but the boundaries 
     852      !!       warning : ptc and pt_low must be masked, but the boundaries 
    744853      !!       conditions on the fluxes are not necessary zalezak (1979) 
    745854      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
     
    747856      !!---------------------------------------------------------------------- 
    748857      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 
     858      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_low      ! before field & upstream guess of after field 
     859      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pfv_low, pfu_low ! upstream flux 
    751860      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux 
    752861      ! 
    753862      INTEGER  ::   ji, jj    ! dummy loop indices 
    754863      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 
     864      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign        !   -      - 
     865      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low 
    757866      !!---------------------------------------------------------------------- 
    758867      zbig = 1.e+40_wp 
    759       zsml = 1.e-15_wp 
    760  
    761       ! test on divergence 
    762       DO jj = 2, jpjm1 
    763          DO ji = fs_2, fs_jpim1   ! vector opt.   
    764             zdiv(ji,jj) =  - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) + pfu_ho(ji,jj) - pfu_ho(ji-1,jj) )   
    765          END DO 
    766       END DO 
    767       CALL lbc_lnk( zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
     868      zsml = epsi20 
    768869 
    769870      ! antidiffusive flux : high order minus low order 
     
    771872      DO jj = 1, jpjm1 
    772873         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  
     874            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj) 
     875            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj) 
     876         END DO 
     877      END DO 
     878 
     879      ! extreme case where pfu_ho has to be zero 
     880      ! ---------------------------------------- 
     881      !                                    pfu_ho 
     882      !                           *         ---> 
     883      !                        |      |  *   |        |  
     884      !                        |      |      |    *   |     
     885      !                        |      |      |        |    * 
     886      !            t_low :       i-1     i       i+1       i+2    
     887      IF( ll_prelimiter ) THEN 
     888 
     889         DO jj = 2, jpjm1 
     890            DO ji = fs_2, fs_jpim1  
     891               zti_low(ji,jj)= pt_low(ji+1,jj  ) 
     892               ztj_low(ji,jj)= pt_low(ji  ,jj+1) 
     893            END DO 
     894         END DO 
     895         CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. ) 
     896 
     897         IF( ll_prelimiter_zalesak ) THEN 
     898 
     899            !! this does not work 
     900!!            DO jj = 2, jpjm1 
     901!!               DO ji = fs_2, fs_jpim1 
     902!!                  IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj  ) - pt_low (ji  ,jj) ) .AND.     & 
     903!!               &      SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj+1) - pt_low (ji  ,jj) )           & 
     904!!               &    ) THEN 
     905!!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji  ,jj) ) .AND.   & 
     906!!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji  ,jj) )         & 
     907!!               &       ) THEN 
     908!!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
     909!!                     ENDIF 
     910!!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj  ) ) .AND.  & 
     911!!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji  ,jj-1) )        & 
     912!!               &       )  THEN 
     913!!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
     914!!                     ENDIF 
     915!!                  ENDIF 
     916!!                END DO 
     917!!            END DO             
     918 
     919            DO jj = 2, jpjm1 
     920               DO ji = fs_2, fs_jpim1 
     921                  IF ( pfu_ho(ji,jj) * ( pt_low(ji+1,jj) - pt_low(ji,jj) ) <= 0. .AND.  & 
     922                     & pfv_ho(ji,jj) * ( pt_low(ji,jj+1) - pt_low(ji,jj) ) <= 0. ) THEN 
     923                     ! 
     924                     IF(  pfu_ho(ji,jj) * ( zti_low(ji+1,jj) - zti_low(ji,jj) ) <= 0 .AND.  & 
     925                        & pfv_ho(ji,jj) * ( ztj_low(ji,jj+1) - ztj_low(ji,jj) ) <= 0)  pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0. 
     926                     ! 
     927                     IF(  pfu_ho(ji,jj) * ( pt_low(ji  ,jj) - pt_low(ji-1,jj) ) <= 0 .AND.  & 
     928                        & pfv_ho(ji,jj) * ( pt_low(ji  ,jj) - pt_low(ji,jj-1) ) <= 0)  pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0. 
     929                     ! 
     930                  ENDIF 
     931               END DO 
     932            END DO 
     933            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
     934 
     935         ELSEIF( ll_prelimiter_devore ) THEN 
     936            z1_dt = 1._wp / pdt 
     937            DO jj = 2, jpjm1 
     938               DO ji = fs_2, fs_jpim1 
     939                  zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) ) 
     940                  pfu_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , & 
     941                     &                          zsign * ( pt_low (ji  ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji  ,jj) * z1_dt , & 
     942                     &                          zsign * ( zti_low(ji+1,jj) - zti_low(ji  ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) ) 
     943                   
     944                  zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) ) 
     945                  pfv_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , & 
     946                     &                          zsign * ( pt_low (ji,jj  ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj  ) * z1_dt , & 
     947                     &                          zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) * e1e2t(ji,jj+1) * z1_dt ) ) 
     948               END DO 
     949            END DO 
     950            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
     951             
     952         ENDIF 
     953          
     954      ENDIF 
     955       
    778956      ! Search local extrema 
    779957      ! -------------------- 
    780       ! max/min of ptc & pt_ups with large negative/positive value (-/+zbig) outside ice cover 
     958      ! max/min of ptc & pt_low with large negative/positive value (-/+zbig) outside ice cover 
    781959      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 
     960         DO ji = 1, jpi 
     961!!clem            IF    ( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) == 0._wp ) THEN 
     962            IF    ( ptc(ji,jj) <= epsi20 .AND. pt_low(ji,jj) <= epsi20 ) THEN 
    784963               zbup(ji,jj) = -zbig 
    785964               zbdo(ji,jj) =  zbig 
    786             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) ) 
     965!!clem            ELSEIF( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) /= 0._wp ) THEN 
     966            ELSEIF( ptc(ji,jj) <= epsi20 .AND. pt_low(ji,jj) > epsi20 ) THEN 
     967               zbup(ji,jj) = pt_low(ji,jj) 
     968               zbdo(ji,jj) = pt_low(ji,jj) 
     969!!clem            ELSEIF( ptc(ji,jj) /= 0._wp .AND. pt_low(ji,jj) == 0._wp ) THEN 
     970            ELSEIF( ptc(ji,jj) > epsi20 .AND. pt_low(ji,jj) <= epsi20 ) THEN 
     971               zbup(ji,jj) = ptc(ji,jj) 
     972               zbdo(ji,jj) = ptc(ji,jj) 
     973            ELSE 
     974               zbup(ji,jj) = MAX( ptc(ji,jj) , pt_low(ji,jj) ) 
     975               zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_low(ji,jj) ) 
    789976            ENDIF 
    790977         END DO 
     
    795982         DO ji = fs_2, fs_jpim1   ! vector opt. 
    796983            ! 
    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 
     984            !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 
     985            !zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) ) 
     986            ! 
     987            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 
     988               &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  ) 
     989            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), & 
     990               &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  ) 
     991            ! 
     992            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &  ! positive/negative part of the flux 
     993               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) ) 
    802994            zneg = MAX( 0., pfu_ho(ji  ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) & 
    803995               & + MAX( 0., pfv_ho(ji,jj  ) ) - MIN( 0., pfv_ho(ji,jj-1) ) 
    804996            ! 
    805997            !                                  ! 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 
     998!            zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt 
     999!            zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt 
     1000            IF( zpos >= epsi10 ) THEN ; zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt 
     1001            ELSE                      ; zbetup(ji,jj) = 0. 
    8121002            ENDIF 
    8131003            ! 
    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 
     1004            IF( zneg >= epsi10 ) THEN ; zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
     1005            ELSE                      ; zbetdo(ji,jj) = 0. 
    8181006            ENDIF 
     1007            ! 
     1008            IF( zbetdo(ji,jj) < 0._wp ) zbetdo(ji,jj)=0. 
     1009            IF( zbetup(ji,jj) < 0._wp ) zbetup(ji,jj)=0. 
     1010            ! 
    8191011            ! 
    8201012         END DO 
     
    8301022            zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) ) 
    8311023            ! 
    832             pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_ups(ji,jj) 
     1024            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_low(ji,jj) 
    8331025         END DO 
    8341026      END DO 
     
    8401032            zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) ) 
    8411033            ! 
    842             pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_ups(ji,jj) 
     1034            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_low(ji,jj) 
    8431035         END DO 
    8441036      END DO 
Note: See TracChangeset for help on using the changeset viewer.