Changeset 10945


Ignore:
Timestamp:
2019-05-07T19:29:25+02:00 (18 months ago)
Author:
clem
Message:

finish repairing the code when only ice dynamics is activated

Location:
NEMO/trunk/src/ICE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_adv_umx.F90

    r10930 r10945  
    3232 
    3333   PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90 
    34        
    35    REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6 
    36    REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    37     
    38    ! advection for S and T:    dVS/dt = -div( uA * uHS / u ) => ll_advS = F 
    39    !                        or dVS/dt = -div( uV * uS  / u ) => ll_advS = T 
    40    LOGICAL ::   ll_advS = .FALSE. 
    4134   ! 
    42    ! alternate directions for upstream 
    43    LOGICAL ::   ll_upsxy = .TRUE. 
     35   INTEGER, PARAMETER ::   np_advS = 1         ! advection for S and T:    dVS/dt = -div(      uVS     ) => np_advS = 1 
     36   !                                                                    or dVS/dt = -div( uA * uHS / u ) => np_advS = 2 
     37   !                                                                    or dVS/dt = -div( uV * uS  / u ) => np_advS = 3 
     38   INTEGER, PARAMETER ::   np_limiter = 1      ! limiter: 1 = nonosc 
     39   !                                                      2 = superbee 
     40   !                                                      3 = h3 
     41   LOGICAL            ::   ll_upsxy  = .TRUE.   ! alternate directions for upstream 
     42   LOGICAL            ::   ll_hoxy   = .TRUE.   ! alternate directions for high order 
     43   LOGICAL            ::   ll_neg    = .TRUE.   ! if T interpolated at u/v points is negative or v_i < 1.e-6 
     44   !                                                 then interpolate T at u/v points using the upstream scheme 
     45   LOGICAL            ::   ll_prelim = .FALSE.  ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D 
    4446   ! 
    45    ! alternate directions for high order 
    46    LOGICAL ::   ll_hoxy = .TRUE. 
     47   REAL(wp)           ::   z1_6   = 1._wp /   6._wp   ! =1/6 
     48   REAL(wp)           ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    4749   ! 
    48    ! if T interpolated at u/v points is negative or v_i < 1.e-6 
    49    !    then interpolate T at u/v points using the upstream scheme 
    50    LOGICAL ::   ll_neg = .TRUE. 
    51    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zmski_small, zmskj_small 
    52    ! 
    53    ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid) 
    54    INTEGER ::   kn_limiter = 1 
    55    ! 
    56    ! prelimiter: use it to avoid overshoot in H 
    57    LOGICAL ::   ll_prelimiter_zalesak = .FALSE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 
     50   INTEGER, ALLOCATABLE, DIMENSION(:,:,:) ::   imsk_small, jmsk_small 
    5851   ! 
    5952   !! * Substitutions 
     
    10093      REAL(wp) ::   zdt, zvi_cen 
    10194      REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication 
    102       REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box  
     95      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
    10396      REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2 
     97      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zu_cat, zv_cat 
    10498      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups 
    105       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups  
    106       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip, z1_vi, z1_vs  
    107       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar 
     99      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai , z1_aip, zhvar 
    108100      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
     101      ! 
     102      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs  
    109103      !!---------------------------------------------------------------------- 
    110104      ! 
     
    150144      zudy(:,:) = pu_ice(:,:) * e2u(:,:) 
    151145      zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 
    152  
     146      ! 
     147      ! setup transport for each ice cat  
     148      DO jl = 1, jpl 
     149         zu_cat(:,:,jl) = zudy(:,:) 
     150         zv_cat(:,:,jl) = zvdx(:,:) 
     151      END DO 
     152      ! 
    153153      ! --- define velocity for advection: u*grad(H) --- ! 
    154154      DO jj = 2, jpjm1 
     
    184184         ! setup a mask where advection will be upstream 
    185185         IF( ll_neg ) THEN 
    186             IF( .NOT. ALLOCATED(zmski_small) )   ALLOCATE( zmski_small(jpi,jpj,jpl) )  
    187             IF( .NOT. ALLOCATED(zmskj_small) )   ALLOCATE( zmskj_small(jpi,jpj,jpl) )  
     186            IF( .NOT. ALLOCATED(imsk_small) )   ALLOCATE( imsk_small(jpi,jpj,jpl) )  
     187            IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_small(jpi,jpj,jpl) )  
    188188            DO jl = 1, jpl 
    189189               DO jj = 1, jpjm1 
    190190                  DO ji = 1, jpim1 
    191191                     zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) ) 
    192                      IF( zvi_cen < epsi06) THEN ; zmski_small(ji,jj,jl) = 0. 
    193                      ELSE                       ; zmski_small(ji,jj,jl) = 1. ; ENDIF 
     192                     IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0 
     193                     ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;  ENDIF 
    194194                     zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) ) 
    195                      IF( zvi_cen < epsi06) THEN ; zmskj_small(ji,jj,jl) = 0. 
    196                      ELSE                       ; zmskj_small(ji,jj,jl) = 1. ; ENDIF 
     195                     IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0 
     196                     ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;  ENDIF 
    197197                  END DO 
    198198               END DO 
     
    204204         ! ----------------------- ! 
    205205         ! 
    206          ! set u*A=u for advection of A only  
    207          DO jl = 1, jpl 
    208             zua_ho(:,:,jl) = zudy(:,:) 
    209             zva_ho(:,:,jl) = zvdx(:,:) 
    210          END DO 
    211  
    212206         !== Ice area ==! 
    213207         zamsk = 1._wp 
    214          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     208         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, & 
    215209            &                                      pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho ) 
    216          zamsk = 0._wp 
    217  
    218          IF( .NOT. ll_advS ) THEN   !-- advection form: uA * uHS / u 
     210         ! 
     211         !                             ! --------------------------------- ! 
     212         IF( np_advS == 1 ) THEN       ! -- advection form: -div( uVS ) -- ! 
     213            !                          ! --------------------------------- ! 
     214            zamsk = 0._wp 
    219215            !== Ice volume ==! 
    220216            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
    221217            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
    222218               &                                      zhvar, pv_i, zua_ups, zva_ups ) 
    223             ! 
    224219            !== Snw volume ==!          
    225220            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     
    227222               &                                      zhvar, pv_s, zua_ups, zva_ups ) 
    228223            ! 
     224            zamsk = 1._wp 
     225            !== Salt content ==! 
     226            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     227               &                                      psv_i, psv_i ) 
     228            !== Ice heat content ==! 
     229            DO jk = 1, nlay_i 
     230               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     231                  &                                      pe_i(:,:,jk,:), pe_i(:,:,jk,:) ) 
     232            END DO 
     233            !== Snw heat content ==! 
     234            DO jk = 1, nlay_s 
     235               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     236                  &                                      pe_s(:,:,jk,:), pe_s(:,:,jk,:) ) 
     237            END DO 
     238            ! 
     239            !                          ! ------------------------------------------ ! 
     240         ELSEIF( np_advS == 2 ) THEN   ! -- advection form: -div( uA * uHS / u ) -- ! 
     241            !                          ! ------------------------------------------ ! 
     242            zamsk = 0._wp 
     243            !== Ice volume ==! 
     244            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     245            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     246               &                                      zhvar, pv_i, zua_ups, zva_ups ) 
     247            !== Snw volume ==!          
     248            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     249            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     250               &                                      zhvar, pv_s, zua_ups, zva_ups ) 
    229251            !== Salt content ==! 
    230252            zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 
    231253            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
    232254               &                                      zhvar, psv_i, zua_ups, zva_ups ) 
    233             ! 
    234255            !== Ice heat content ==! 
    235256            DO jk = 1, nlay_i 
     
    238259                  &                                      zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups ) 
    239260            END DO 
    240             ! 
    241261            !== Snw heat content ==! 
    242262            DO jk = 1, nlay_s 
     
    246266            END DO 
    247267            ! 
    248          ELSE                       !-- advection form: uV * uS / u 
     268            !                          ! ----------------------------------------- ! 
     269         ELSEIF( np_advS == 3 ) THEN   ! -- advection form: -div( uV * uS / u ) -- ! 
     270            !                          ! ----------------------------------------- ! 
     271            zamsk = 0._wp 
     272            ! 
     273            ALLOCATE( zuv_ho (jpi,jpj,jpl), zvv_ho (jpi,jpj,jpl),  & 
     274               &      zuv_ups(jpi,jpj,jpl), zvv_ups(jpi,jpj,jpl), z1_vi(jpi,jpj,jpl), z1_vs(jpi,jpj,jpl) ) 
     275            ! 
    249276            ! inverse of Vi 
    250277            WHERE( pv_i(:,:,:) >= epsi20 )   ;   z1_vi(:,:,:) = 1._wp / pv_i(:,:,:) 
     
    264291            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
    265292               &                                      zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 
    266             ! 
    267293            !== Salt content ==! 
    268294            zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:) 
    269295            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, & 
    270296               &                                      zhvar, psv_i, zuv_ups, zvv_ups ) 
    271             ! 
    272297            !== Ice heat content ==! 
    273298            DO jk = 1, nlay_i 
     
    276301                  &                                      zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups ) 
    277302            END DO 
    278             ! 
    279303            !== Snow volume ==!          
    280304            zuv_ups = zua_ups 
     
    283307            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
    284308               &                                      zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 
    285             ! 
    286309            !== Snw heat content ==! 
    287310            DO jk = 1, nlay_s 
     
    291314            END DO 
    292315            ! 
     316            DEALLOCATE( zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs ) 
     317            ! 
    293318         ENDIF 
    294319         ! 
    295320         !== Ice age ==! 
    296321         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 
    297             ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with 
    298             !       fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that 
    299             !       spread into the domain. Therefore we cheat and consider that ice age should be advected as ice concentration 
    300             !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 
    301             !!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i ) 
    302             ! set u*A=u for advection of ice age 
    303             DO jl = 1, jpl 
    304                zua_ho(:,:,jl) = zudy(:,:) 
    305                zva_ho(:,:,jl) = zvdx(:,:) 
    306             END DO 
    307322            zamsk = 1._wp 
    308             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho, zva_ho, zcu_box, zcv_box, & 
     323            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
    309324               &                                      poa_i, poa_i ) 
    310             zamsk = 0._wp 
    311325         ENDIF 
    312326         ! 
    313327         !== melt ponds ==! 
    314328         IF ( ln_pnd_H12 ) THEN 
    315             ! set u*A=u for advection of Ap only  
    316             DO jl = 1, jpl 
    317                zua_ho(:,:,jl) = zudy(:,:) 
    318                zva_ho(:,:,jl) = zvdx(:,:) 
    319             END DO 
    320329            ! fraction 
    321330            zamsk = 1._wp 
    322             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     331            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, & 
    323332               &                                      pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho ) 
     333            ! volume 
    324334            zamsk = 0._wp 
    325             ! volume 
    326335            zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:) 
    327336            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     
    341350         ! 
    342351         ! --- Ensure non-negative fields and in-bound thicknesses --- ! 
    343          ! 
    344352         ! Remove negative values (conservation is ensured) 
    345353         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
     
    474482         CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. ) 
    475483         ! 
    476          IF    ( kn_limiter == 1 ) THEN 
     484         IF    ( np_limiter == 1 ) THEN 
    477485            CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 
    478          ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN 
     486         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN 
    479487            CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho ) 
    480488            CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho ) 
     
    659667         END DO 
    660668         ! 
    661          IF    ( kn_limiter == 1 ) THEN 
     669         IF    ( np_limiter == 1 ) THEN 
    662670            CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    663          ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN 
     671         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN 
    664672            CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    665673            CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     
    677685               END DO 
    678686            END DO 
    679             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     687            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    680688 
    681689            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
     
    698706               END DO 
    699707            END DO 
    700             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     708            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    701709 
    702710         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     
    709717               END DO 
    710718            END DO 
    711             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     719            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    712720            ! 
    713721            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
     
    730738               END DO 
    731739            END DO 
    732             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     740            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    733741 
    734742         ENDIF 
    735          IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     743         IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    736744          
    737745      ENDIF 
     
    771779         CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 
    772780         !                                                        !--  limiter in x --! 
    773          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     781         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    774782         !                                                        !--  advective form update in zpt  --! 
    775783         DO jl = 1, jpl 
     
    792800         ENDIF 
    793801         !                                                        !--  limiter in y --! 
    794          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     802         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    795803         !          
    796804         ! 
     
    800808         CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 
    801809         !                                                        !--  limiter in y --! 
    802          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     810         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    803811         !                                                        !--  advective form update in zpt  --! 
    804812         DO jl = 1, jpl 
     
    821829         ENDIF 
    822830         !                                                        !--  limiter in x --! 
    823          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     831         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    824832         ! 
    825833      ENDIF 
    826834 
    827       IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     835      IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    828836      ! 
    829837   END SUBROUTINE macho 
     
    967975            DO jj = 1, jpjm1 
    968976               DO ji = 1, fs_jpim1 
    969                   IF( pt_u(ji,jj,jl) < 0._wp .OR. ( zmski_small(ji,jj,jl) == 0. .AND. pamsk == 0. ) ) THEN 
     977                  IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    970978                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    971979                        &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     
    11211129            DO jj = 1, jpjm1 
    11221130               DO ji = 1, fs_jpim1 
    1123                   IF( pt_v(ji,jj,jl) < 0._wp .OR. ( zmskj_small(ji,jj,jl) == 0. .AND. pamsk == 0. ) ) THEN 
     1131                  IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    11241132                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    11251133                        &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     
    11851193      !                        |      |      |        |    * 
    11861194      !            t_ups :       i-1     i       i+1       i+2    
    1187       IF( ll_prelimiter_zalesak ) THEN 
     1195      IF( ll_prelim ) THEN 
    11881196          
    11891197         DO jl = 1, jpl 
     
    13521360               Rjp = zslpx(ji+1,jj,jl) 
    13531361 
    1354                IF( kn_limiter == 3 ) THEN 
     1362               IF( np_limiter == 3 ) THEN 
    13551363 
    13561364                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     
    13681376                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 
    13691377 
    1370                ELSEIF( kn_limiter == 2 ) THEN 
     1378               ELSEIF( np_limiter == 2 ) THEN 
    13711379                  IF( Rj /= 0. ) THEN 
    13721380                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
     
    14471455               Rjp = zslpy(ji,jj+1,jl) 
    14481456 
    1449                IF( kn_limiter == 3 ) THEN 
     1457               IF( np_limiter == 3 ) THEN 
    14501458 
    14511459                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     
    14631471                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 
    14641472 
    1465                ELSEIF( kn_limiter == 2 ) THEN 
     1473               ELSEIF( np_limiter == 2 ) THEN 
    14661474 
    14671475                  IF( Rj /= 0. ) THEN 
  • NEMO/trunk/src/ICE/icedyn_rdgrft.F90

    r10531 r10945  
    945945         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 
    946946      ENDIF 
    947       !                              ! allocate tke arrays 
     947      ! 
     948      IF( .NOT. ln_icethd ) THEN 
     949         rn_porordg = 0._wp 
     950         rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp 
     951         rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp 
     952         IF( lwp ) THEN 
     953            WRITE(numout,*) '      ==> only ice dynamics is activated, thus some parameters must be changed' 
     954            WRITE(numout,*) '            rn_porordg   = ', rn_porordg 
     955            WRITE(numout,*) '            rn_fsnwrdg   = ', rn_fsnwrdg  
     956            WRITE(numout,*) '            rn_fpndrdg   = ', rn_fpndrdg  
     957            WRITE(numout,*) '            rn_fsnwrft   = ', rn_fsnwrft  
     958            WRITE(numout,*) '            rn_fpndrft   = ', rn_fpndrft  
     959         ENDIF 
     960      ENDIF 
     961      !                              ! allocate arrays 
    948962      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' ) 
    949963      ! 
  • NEMO/trunk/src/ICE/icevar.F90

    r10930 r10945  
    563563      DO jl = 1, jpl       !==  loop over the categories  ==! 
    564564         ! 
     565         ! make sure a_i=0 where v_i<=0 
     566         WHERE( pv_i(:,:,:) <= 0._wp )   pa_i(:,:,:) = 0._wp 
     567 
    565568         !---------------------------------------- 
    566569         ! zap ice energy and send it to the ocean 
     
    569572            DO jj = 1 , jpj 
    570573               DO ji = 1 , jpi 
    571                   IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
     574                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
    572575                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 
    573576                     pe_i(ji,jj,jk,jl) = 0._wp 
     
    580583            DO jj = 1 , jpj 
    581584               DO ji = 1 , jpi 
    582                   IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
     585                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
    583586                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 
    584587                     pe_s(ji,jj,jk,jl) = 0._wp 
     
    593596         DO jj = 1 , jpj 
    594597            DO ji = 1 , jpi 
    595                IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
     598               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
    596599                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 
    597600                  pv_i   (ji,jj,jl) = 0._wp 
    598601               ENDIF 
    599                IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
     602               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
    600603                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt 
    601604                  pv_s   (ji,jj,jl) = 0._wp 
    602605               ENDIF 
    603                IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
     606               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
    604607                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 
    605608                  psv_i  (ji,jj,jl) = 0._wp 
Note: See TracChangeset for help on using the changeset viewer.