Changeset 13618


Ignore:
Timestamp:
2020-10-16T10:07:55+02:00 (3 months ago)
Author:
clem
Message:

trunk: optimization of both Prather and UMx advection schemes. Related to ticket #2552

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

Legend:

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

    r13613 r13618  
    111111      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
    112112      END WHERE 
    113       DO jl = 1, jpl 
    114          DO_2D( 0, 0, 0, 0 ) 
    115             zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    116                &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    117                &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    118                &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    119             zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    120                &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    121                &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    122                &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    123             zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    124                &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    125                &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    126                &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    127             zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
    128                &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
    129                &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
    130                &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
    131          END_2D 
    132       END DO 
     113      CALL icemax3D( ph_i , zhi_max ) 
     114      CALL icemax3D( ph_s , zhs_max ) 
     115      CALL icemax3D( ph_ip, zhip_max) 
     116      CALL icemax3D( zs_i , zsi_max ) 
    133117      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
    134118      ! 
     
    143127         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
    144128         END WHERE 
    145       END DO 
    146       DO jl = 1, jpl 
    147          DO_3D( 0, 0, 0, 0, 1, nlay_i ) 
    148             zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), & 
    149                &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
    150                &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
    151                &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
    152          END_3D 
    153       END DO 
    154       DO jl = 1, jpl 
    155          DO_3D( 0, 0, 0, 0, 1, nlay_s ) 
    156             zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), & 
    157                &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
    158                &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
    159                &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
    160          END_3D 
    161       END DO 
    162       CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
    163       CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
     129      END DO    
     130      CALL icemax4D( ze_i , zei_max ) 
     131      CALL icemax4D( ze_s , zes_max ) 
     132      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 
     133      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 
    164134      ! 
    165135      ! 
     
    399369      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      - 
    400370      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      - 
     371      REAL(wp) ::   zpsm, zps0 
     372      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    401373      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace 
    402374      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      - 
     
    417389            DO ji = Nis0 - 1, Nie0 + 1 
    418390 
     391               zpsm  = psm (ji,jj,jl) ! optimization 
     392               zps0  = ps0 (ji,jj,jl) 
     393               zpsx  = psx (ji,jj,jl) 
     394               zpsxx = psxx(ji,jj,jl) 
     395               zpsy  = psy (ji,jj,jl) 
     396               zpsyy = psyy(ji,jj,jl) 
     397               zpsxy = psxy(ji,jj,jl) 
     398 
    419399               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
    420                psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 
    421                ! 
    422                zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     400               zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 
     401               ! 
     402               zslpmax = MAX( 0._wp, zps0 ) 
    423403               zs1max  = 1.5 * zslpmax 
    424                zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
    425                zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 
     404               zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) ) 
     405               zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 
    426406               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    427                 
    428                ps0 (ji,jj,jl) = zslpmax   
    429                psx (ji,jj,jl) = zs1new         * rswitch 
    430                psxx(ji,jj,jl) = zs2new         * rswitch 
    431                psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 
    432                psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 
    433                psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    434                 
     407 
     408               zps0 = zslpmax   
     409               zpsx  = zs1new  * rswitch 
     410               zpsxx = zs2new  * rswitch 
     411               zpsy  = zpsy    * rswitch 
     412               zpsyy = zpsyy  * rswitch 
     413               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     414 
    435415               !  Calculate fluxes and moments between boxes i<-->i+1               
    436416               !                                !  Flux from i to i+1 WHEN u GT 0  
    437417               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    438                zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
     418               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 
    439419               zalfq        =  zalf * zalf 
    440420               zalf1        =  1.0 - zalf 
    441421               zalf1q       =  zalf1 * zalf1 
    442422               ! 
    443                zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl) 
    444                zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 
    445                zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 
    446                zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq 
    447                zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    448                zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl) 
    449                zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl) 
    450                 
     423               zfm (ji,jj)  =  zalf  *   zpsm  
     424               zf0 (ji,jj)  =  zalf  * ( zps0  + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 
     425               zfx (ji,jj)  =  zalfq * ( zpsx  + 3.0 * zalf1 * zpsxx ) 
     426               zfxx(ji,jj)  =  zalf  *   zpsxx * zalfq 
     427               zfy (ji,jj)  =  zalf  * ( zpsy  + zalf1 * zpsxy ) 
     428               zfxy(ji,jj)  =  zalfq *   zpsxy 
     429               zfyy(ji,jj)  =  zalf  *   zpsyy 
     430 
    451431               !                                !  Readjust moments remaining in the box. 
    452                psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    453                ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    454                psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 
    455                psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl) 
    456                psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj) 
    457                psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj) 
    458                psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     432               zpsm  =  zpsm  - zfm(ji,jj) 
     433               zps0  =  zps0  - zf0(ji,jj) 
     434               zpsx  =  zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 
     435               zpsxx =  zalf1  * zalf1q * zpsxx 
     436               zpsy  =  zpsy  - zfy (ji,jj) 
     437               zpsyy =  zpsyy - zfyy(ji,jj) 
     438               zpsxy =  zalf1q * zpsxy 
     439               ! 
     440               psm (ji,jj,jl) = zpsm ! optimization 
     441               ps0 (ji,jj,jl) = zps0  
     442               psx (ji,jj,jl) = zpsx  
     443               psxx(ji,jj,jl) = zpsxx 
     444               psy (ji,jj,jl) = zpsy  
     445               psyy(ji,jj,jl) = zpsyy 
     446               psxy(ji,jj,jl) = zpsxy 
     447               ! 
    459448            END DO 
    460449             
     
    480469 
    481470            DO ji = Nis0, Nie0 
     471               ! 
     472               zpsm  = psm (ji,jj,jl) ! optimization 
     473               zps0  = ps0 (ji,jj,jl) 
     474               zpsx  = psx (ji,jj,jl) 
     475               zpsxx = psxx(ji,jj,jl) 
     476               zpsy  = psy (ji,jj,jl) 
     477               zpsyy = psyy(ji,jj,jl) 
     478               zpsxy = psxy(ji,jj,jl) 
    482479               !                                !  Readjust moments remaining in the box. 
    483480               zbt  =       zbet(ji-1,jj) 
    484481               zbt1 = 1.0 - zbet(ji-1,jj) 
    485482               ! 
    486                psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 
    487                ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 
    488                psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 
    489                psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 
    490                psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 
    491                psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 
    492                psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 
    493                 
     483               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 
     484               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 
     485               zpsx  = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 
     486               zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 
     487               zpsy  = zbt * zpsy  + zbt1 * ( zpsy - zfy (ji-1,jj) ) 
     488               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 
     489               zpsxy = zalg1q(ji-1,jj) * zpsxy 
     490 
    494491               !   Put the temporary moments into appropriate neighboring boxes.     
    495492               !                                !   Flux from i to i+1 IF u GT 0. 
    496                zbt  =       zbet(ji-1,jj) 
    497                zbt1 = 1.0 - zbet(ji-1,jj) 
    498                psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 
    499                zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 
    500                zalf1         = 1.0 - zalf 
    501                ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 
    502                ! 
    503                ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 
    504                psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 
    505                psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             & 
    506                   &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 
    507                   &            + zbt1 * psxx(ji,jj,jl) 
    508                psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             & 
    509                   &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   & 
    510                   &            + zbt1 * psxy(ji,jj,jl) 
    511                psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 
    512                psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 
    513                 
     493               zbt   =       zbet(ji-1,jj) 
     494               zbt1  = 1.0 - zbet(ji-1,jj) 
     495               zpsm  = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 
     496               zalf  = zbt * zfm(ji-1,jj) / zpsm 
     497               zalf1 = 1.0 - zalf 
     498               ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 
     499               ! 
     500               zps0  =  zbt  * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 
     501               zpsx  =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 
     502               zpsxx =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx                            & 
     503                  &            + 5.0 * ( zalf * zalf1 * ( zpsx  - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 
     504                  &            + zbt1 * zpsxx 
     505               zpsxy =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy            & 
     506                  &            + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * zpsy ) )  & 
     507                  &            + zbt1 * zpsxy 
     508               zpsy  =  zbt  * ( zpsy  + zfy (ji-1,jj) ) + zbt1 * zpsy  
     509               zpsyy =  zbt  * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 
     510 
    514511               !                                !  Flux from i+1 to i IF u LT 0. 
    515                zbt  =       zbet(ji,jj) 
    516                zbt1 = 1.0 - zbet(ji,jj) 
    517                psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    518                zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    519                zalf1         = 1.0 - zalf 
    520                ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    521                ! 
    522                ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 
    523                psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 
    524                psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 
    525                   &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    & 
    526                   &                                           + ( zalf1 - zalf ) * ztemp ) ) 
    527                psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    528                   &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 
    529                psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 
    530                psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 
     512               zbt   =       zbet(ji,jj) 
     513               zbt1  = 1.0 - zbet(ji,jj) 
     514               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     515               zalf  = zbt1 * zfm(ji,jj) / zpsm 
     516               zalf1 = 1.0 - zalf 
     517               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     518               ! 
     519               zps0  = zbt * zps0  + zbt1 * ( zps0 + zf0(ji,jj) ) 
     520               zpsx  = zbt * zpsx  + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 
     521               zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 
     522                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) )    & 
     523                  &                         + ( zalf1 - zalf ) * ztemp ) ) 
     524               zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     525                  &                         + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 
     526               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  + zfy (ji,jj) ) 
     527               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 
     528               ! 
     529               psm (ji,jj,jl) = zpsm  ! optimization 
     530               ps0 (ji,jj,jl) = zps0  
     531               psx (ji,jj,jl) = zpsx  
     532               psxx(ji,jj,jl) = zpsxx 
     533               psy (ji,jj,jl) = zpsy  
     534               psyy(ji,jj,jl) = zpsyy 
     535               psxy(ji,jj,jl) = zpsxy 
    531536            END DO 
    532537            ! 
     
    559564      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
    560565      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     566      REAL(wp) ::   zpsm, zps0 
     567      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    561568      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
    562569      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
     
    574581         ! Limitation of moments. 
    575582         DO_2D( 1, 1, ji0, ji0 ) 
    576             !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    577             psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  ) 
    578             ! 
    579             zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     583            ! 
     584            zpsm  = psm (ji,jj,jl) ! optimization 
     585            zps0  = ps0 (ji,jj,jl) 
     586            zpsx  = psx (ji,jj,jl) 
     587            zpsxx = psxx(ji,jj,jl) 
     588            zpsy  = psy (ji,jj,jl) 
     589            zpsyy = psyy(ji,jj,jl) 
     590            zpsxy = psxy(ji,jj,jl) 
     591            ! 
     592            !  Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 
     593            zpsm = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20  ) 
     594            ! 
     595            zslpmax = MAX( 0._wp, zps0 ) 
    580596            zs1max  = 1.5 * zslpmax 
    581             zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
    582             zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 
     597            zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) ) 
     598            zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 
    583599            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    584600            ! 
    585             ps0 (ji,jj,jl) = zslpmax   
    586             psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 
    587             psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 
    588             psy (ji,jj,jl) = zs1new         * rswitch 
    589             psyy(ji,jj,jl) = zs2new         * rswitch 
    590             psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    591   
     601            zps0 = zslpmax   
     602            zpsx  = zpsx * rswitch 
     603            zpsxx = zpsxx * rswitch 
     604            zpsy = zs1new         * rswitch 
     605            zpsyy = zs2new         * rswitch 
     606            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     607 
    592608            !  Calculate fluxes and moments between boxes j<-->j+1               
    593609            !                                !  Flux from j to j+1 WHEN v GT 0    
    594610            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    595             zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
     611            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 
    596612            zalfq        =  zalf * zalf 
    597613            zalf1        =  1.0 - zalf 
    598614            zalf1q       =  zalf1 * zalf1 
    599615            ! 
    600             zfm (ji,jj)  =  zalf  * psm(ji,jj,jl) 
    601             zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) )  
    602             zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 
    603             zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl) 
    604             zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    605             zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl) 
    606             zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl) 
     616            zfm (ji,jj)  =  zalf  * zpsm 
     617            zf0 (ji,jj)  =  zalf  * ( zps0 + zalf1 * ( zpsy  + (zalf1-zalf) * zpsyy ) )  
     618            zfy (ji,jj)  =  zalfq *( zpsy + 3.0*zalf1*zpsyy ) 
     619            zfyy(ji,jj)  =  zalf  * zalfq * zpsyy 
     620            zfx (ji,jj)  =  zalf  * ( zpsx + zalf1 * zpsxy ) 
     621            zfxy(ji,jj)  =  zalfq * zpsxy 
     622            zfxx(ji,jj)  =  zalf  * zpsxx 
    607623            ! 
    608624            !                                !  Readjust moments remaining in the box. 
    609             psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    610             ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    611             psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 
    612             psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl) 
    613             psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj) 
    614             psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj) 
    615             psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     625            zpsm   =  zpsm  - zfm(ji,jj) 
     626            zps0   =  zps0  - zf0(ji,jj) 
     627            zpsy   =  zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 
     628            zpsyy  =  zalf1 * zalf1q * zpsyy 
     629            zpsx   =  zpsx  - zfx(ji,jj) 
     630            zpsxx  =  zpsxx - zfxx(ji,jj) 
     631            zpsxy  =  zalf1q * zpsxy 
     632            ! 
     633            psm (ji,jj,jl) = zpsm ! optimization 
     634            ps0 (ji,jj,jl) = zps0  
     635            psx (ji,jj,jl) = zpsx  
     636            psxx(ji,jj,jl) = zpsxx 
     637            psy (ji,jj,jl) = zpsy  
     638            psyy(ji,jj,jl) = zpsyy 
     639            psxy(ji,jj,jl) = zpsxy 
    616640         END_2D 
    617641         ! 
     
    641665            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    642666            ! 
    643             psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 
    644             ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 
    645             psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 
    646             psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 
    647             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 
    648             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 
    649             psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 
     667            zpsm  = psm (ji,jj,jl) ! optimization 
     668            zps0  = ps0 (ji,jj,jl) 
     669            zpsx  = psx (ji,jj,jl) 
     670            zpsxx = psxx(ji,jj,jl) 
     671            zpsy  = psy (ji,jj,jl) 
     672            zpsyy = psyy(ji,jj,jl) 
     673            zpsxy = psxy(ji,jj,jl) 
     674            ! 
     675            zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 
     676            zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 
     677            zpsy  = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 
     678            zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 
     679            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  - zfx (ji,jj-1) ) 
     680            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 
     681            zpsxy = zalg1q(ji,jj-1) * zpsxy 
    650682 
    651683            !   Put the temporary moments into appropriate neighboring boxes.     
    652684            !                                !   Flux from j to j+1 IF v GT 0. 
    653             zbt  =       zbet(ji,jj-1) 
    654             zbt1 = 1.0 - zbet(ji,jj-1) 
    655             psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)  
    656             zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)  
    657             zalf1         = 1.0 - zalf 
    658             ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 
    659             ! 
    660             ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 
    661             psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  & 
    662                &             + zbt1 * psy(ji,jj,jl)   
    663             psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           & 
    664                &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
    665                &             + zbt1 * psyy(ji,jj,jl) 
    666             psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            & 
    667                &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  & 
    668                &             + zbt1 * psxy(ji,jj,jl) 
    669             psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 
    670             psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 
     685            zbt   =       zbet(ji,jj-1) 
     686            zbt1  = 1.0 - zbet(ji,jj-1) 
     687            zpsm  = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm  
     688            zalf  = zbt * zfm(ji,jj-1) / zpsm  
     689            zalf1 = 1.0 - zalf 
     690            ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 
     691            ! 
     692            zps0  =   zbt  * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 
     693            zpsy  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp )  & 
     694               &             + zbt1 * zpsy   
     695            zpsyy =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy                           & 
     696               &             + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
     697               &             + zbt1 * zpsyy 
     698            zpsxy =   zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy             & 
     699               &             + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )  & 
     700               &             + zbt1 * zpsxy 
     701            zpsx  =   zbt * ( zpsx  + zfx (ji,jj-1) ) + zbt1 * zpsx  
     702            zpsxx =   zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 
    671703 
    672704            !                                !  Flux from j+1 to j IF v LT 0. 
    673             zbt  =       zbet(ji,jj) 
    674             zbt1 = 1.0 - zbet(ji,jj) 
    675             psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    676             zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    677             zalf1         = 1.0 - zalf 
    678             ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    679             ! 
    680             ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) ) 
    681             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 
    682             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 
    683                &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    & 
    684                &                                            + ( zalf1 - zalf ) * ztemp ) ) 
    685             psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    686                &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 
    687             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 
    688             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 
     705            zbt   =       zbet(ji,jj) 
     706            zbt1  = 1.0 - zbet(ji,jj) 
     707            zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     708            zalf  = zbt1 * zfm(ji,jj) / zpsm 
     709            zalf1 = 1.0 - zalf 
     710            ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     711            ! 
     712            zps0  = zbt * zps0  + zbt1 * (  zps0 + zf0(ji,jj) ) 
     713            zpsy  = zbt * zpsy  + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 
     714            zpsyy = zbt * zpsyy + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 
     715               &                         + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) )     & 
     716               &                         + ( zalf1 - zalf ) * ztemp ) ) 
     717            zpsxy = zbt * zpsxy + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     718               &                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 
     719            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  + zfx (ji,jj) ) 
     720            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 
     721            ! 
     722            psm (ji,jj,jl) = zpsm ! optimization 
     723            ps0 (ji,jj,jl) = zps0  
     724            psx (ji,jj,jl) = zpsx  
     725            psxx(ji,jj,jl) = zpsxx 
     726            psy (ji,jj,jl) = zpsy  
     727            psyy(ji,jj,jl) = zpsyy 
     728            psxy(ji,jj,jl) = zpsxy 
    689729         END_2D 
    690730         ! 
     
    11011141   END SUBROUTINE adv_pra_rst 
    11021142 
     1143   SUBROUTINE icemax3D( pice , pmax ) 
     1144      !!--------------------------------------------------------------------- 
     1145      !!                   ***  ROUTINE icemax3D ***                      
     1146      !! ** Purpose :  compute the max of the 9 points around 
     1147      !!---------------------------------------------------------------------- 
     1148      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input 
     1149      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output 
     1150      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1151      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1152      !!---------------------------------------------------------------------- 
     1153      DO jl = 1, jpl 
     1154         DO jj = Njs0-1, Nje0+1     
     1155            DO ji = Nis0, Nie0 
     1156               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 
     1157            END DO 
     1158         END DO 
     1159         DO jj = Njs0, Nje0     
     1160            DO ji = Nis0, Nie0 
     1161               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1162            END DO 
     1163         END DO 
     1164      END DO 
     1165   END SUBROUTINE icemax3D 
     1166 
     1167   SUBROUTINE icemax4D( pice , pmax ) 
     1168      !!--------------------------------------------------------------------- 
     1169      !!                   ***  ROUTINE icemax4D ***                      
     1170      !! ** Purpose :  compute the max of the 9 points around 
     1171      !!---------------------------------------------------------------------- 
     1172      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input 
     1173      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output 
     1174      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1175      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices 
     1176      !!---------------------------------------------------------------------- 
     1177      jlay = SIZE( pice , 3 )   ! size of input arrays 
     1178      DO jl = 1, jpl 
     1179         DO jk = 1, jlay 
     1180            DO jj = Njs0-1, Nje0+1     
     1181               DO ji = Nis0, Nie0 
     1182                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 
     1183               END DO 
     1184            END DO 
     1185            DO jj = Njs0, Nje0     
     1186               DO ji = Nis0, Nie0 
     1187                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1188               END DO 
     1189            END DO 
     1190         END DO 
     1191      END DO 
     1192   END SUBROUTINE icemax4D 
     1193 
    11031194#else 
    11041195   !!---------------------------------------------------------------------- 
  • NEMO/trunk/src/ICE/icedyn_adv_umx.F90

    r13609 r13618  
    115115      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
    116116      END WHERE 
    117       DO jl = 1, jpl 
    118          DO_2D( 0, 0, 0, 0 ) 
    119             zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    120                &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    121                &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    122                &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    123             zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    124                &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    125                &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    126                &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    127             zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    128                &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    129                &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    130                &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    131             zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
    132                &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
    133                &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
    134                &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
    135          END_2D 
    136       END DO 
     117      CALL icemax3D( ph_i , zhi_max ) 
     118      CALL icemax3D( ph_s , zhs_max ) 
     119      CALL icemax3D( ph_ip, zhip_max) 
     120      CALL icemax3D( zs_i , zsi_max ) 
    137121      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
    138122      ! 
     
    147131         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
    148132         END WHERE 
    149       END DO 
    150       DO jl = 1, jpl 
    151          DO_3D( 0, 0, 0, 0, 1, nlay_i ) 
    152             zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), & 
    153                &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
    154                &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
    155                &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
    156          END_3D 
    157       END DO 
    158       DO jl = 1, jpl 
    159          DO_3D( 0, 0, 0, 0, 1, nlay_s ) 
    160             zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), & 
    161                &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
    162                &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
    163                &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
    164          END_3D 
    165       END DO 
    166       CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
    167       CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
     133      END DO    
     134      CALL icemax4D( ze_i , zei_max ) 
     135      CALL icemax4D( ze_s , zes_max ) 
     136      CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp ) 
     137      CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp ) 
    168138      ! 
    169139      ! 
     
    386356            ENDIF 
    387357         ENDIF 
     358 
     359         ! --- Lateral boundary conditions --- ! 
     360         IF    ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 
     361            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 
     362               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 
     363         ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN 
     364            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 
     365               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 
     366         ELSE 
     367            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp ) 
     368         ENDIF 
     369         CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) 
     370         CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp ) 
    388371         ! 
    389372         !== Open water area ==! 
     
    393376               &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
    394377         END_2D 
    395          CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1.0_wp ) 
     378         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1._wp ) 
    396379         ! 
    397380         ! --- diagnostics --- ! 
     
    462445      !!             work on H (and not V). It is partly related to the multi-category approach 
    463446      !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 
    464       !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 
    465       !!             since sv_i and e_i are still good. 
     447      !!             concentration is small). We also limit S and T. 
    466448      !!---------------------------------------------------------------------- 
    467449      REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     
    507489      IF( pamsk == 0._wp ) THEN 
    508490         DO jl = 1, jpl 
    509             DO_2D( 1, 0, 1, 0 ) 
     491            DO_2D( 0, 0, 1, 0 ) 
    510492               IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 
    511493                  zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj) 
     
    516498               ENDIF 
    517499               ! 
     500            END_2D 
     501            DO_2D( 1, 0, 0, 0 ) 
    518502               IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 
    519503                  zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj) 
     
    550534      IF( PRESENT( pua_ho ) ) THEN 
    551535         DO jl = 1, jpl 
    552             DO_2D( 1, 0, 1, 0 ) 
    553                pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
    554                pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
     536            DO_2D( 0, 0, 1, 0 ) 
     537               pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) 
     538               pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) 
     539            END_2D 
     540            DO_2D( 1, 0, 0, 0 ) 
     541               pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
     542               pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
    555543            END_2D 
    556544         END DO 
     
    566554         END_2D 
    567555      END DO 
    568       CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1.0_wp ) 
    569556      ! 
    570557   END SUBROUTINE adv_umx 
     
    605592            ! 
    606593            DO jl = 1, jpl              !-- flux in x-direction 
    607                DO_2D( 1, 0, 1, 0 ) 
     594               DO_2D( 1, 1, 1, 0 ) 
    608595                  pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
    609596               END_2D 
     
    611598            ! 
    612599            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    613                DO_2D( 0, 0, 0, 0 ) 
     600               DO_2D( 1, 1, 0, 0 ) 
    614601                  ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              & 
    615602                     &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     
    618605               END_2D 
    619606            END DO 
    620             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    621607            ! 
    622608            DO jl = 1, jpl              !-- flux in y-direction 
    623                DO_2D( 1, 0, 1, 0 ) 
     609               DO_2D( 1, 0, 0, 0 ) 
    624610                  pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 
    625611               END_2D 
     
    629615            ! 
    630616            DO jl = 1, jpl              !-- flux in y-direction 
    631                DO_2D( 1, 0, 1, 0 ) 
     617               DO_2D( 1, 0, 1, 1 ) 
    632618                  pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
    633619               END_2D 
     
    635621            ! 
    636622            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    637                DO_2D( 0, 0, 0, 0 ) 
     623               DO_2D( 0, 0, 1, 1 ) 
    638624                  ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  & 
    639625                     &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     
    642628               END_2D 
    643629            END DO 
    644             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    645630            ! 
    646631            DO jl = 1, jpl              !-- flux in x-direction 
    647                DO_2D( 1, 0, 1, 0 ) 
     632               DO_2D( 0, 0, 1, 0 ) 
    648633                  pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 
    649634               END_2D 
     
    694679         ! 
    695680         DO jl = 1, jpl 
    696             DO_2D( 1, 0, 1, 0 ) 
     681            DO_2D( 1, 1, 1, 0 ) 
    697682               pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) ) 
     683            END_2D 
     684            DO_2D( 1, 0, 1, 1 ) 
    698685               pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) ) 
    699686            END_2D 
     
    712699            ! 
    713700            DO jl = 1, jpl              !-- flux in x-direction 
    714                DO_2D( 1, 0, 1, 0 ) 
     701               DO_2D( 1, 1, 1, 0 ) 
    715702                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
    716703               END_2D 
     
    719706 
    720707            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    721                DO_2D( 0, 0, 0, 0 ) 
     708               DO_2D( 1, 1, 0, 0 ) 
    722709                  ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              & 
    723710                     &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     
    726713               END_2D 
    727714            END DO 
    728             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    729715 
    730716            DO jl = 1, jpl              !-- flux in y-direction 
    731                DO_2D( 1, 0, 1, 0 ) 
     717               DO_2D( 1, 0, 0, 0 ) 
    732718                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 
    733719               END_2D 
     
    738724            ! 
    739725            DO jl = 1, jpl              !-- flux in y-direction 
    740                DO_2D( 1, 0, 1, 0 ) 
     726               DO_2D( 1, 0, 1, 1 ) 
    741727                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
    742728               END_2D 
     
    745731            ! 
    746732            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    747                DO_2D( 0, 0, 0, 0 ) 
     733               DO_2D( 0, 0, 1, 1 ) 
    748734                  ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  & 
    749735                     &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     
    752738               END_2D 
    753739            END DO 
    754             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    755740            ! 
    756741            DO jl = 1, jpl              !-- flux in x-direction 
    757                DO_2D( 1, 0, 1, 0 ) 
     742               DO_2D( 0, 0, 1, 0 ) 
    758743                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 
    759744               END_2D 
     
    912897         !         
    913898         DO jl = 1, jpl 
    914             DO_2D( 1, 0, 1, 0 ) 
     899            DO_2D( 0, 0, 1, 0 ) 
    915900               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    916901                  &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     
    921906         ! 
    922907         DO jl = 1, jpl 
    923             DO_2D( 1, 0, 1, 0 ) 
     908            DO_2D( 0, 0, 1, 0 ) 
    924909               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    925910               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     
    931916         ! 
    932917         DO jl = 1, jpl 
    933             DO_2D( 1, 0, 1, 0 ) 
     918            DO_2D( 0, 0, 1, 0 ) 
    934919               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    935920               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     
    945930         ! 
    946931         DO jl = 1, jpl 
    947             DO_2D( 1, 0, 1, 0 ) 
     932            DO_2D( 0, 0, 1, 0 ) 
    948933               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    949934               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     
    959944         ! 
    960945         DO jl = 1, jpl 
    961             DO_2D( 1, 0, 1, 0 ) 
     946            DO_2D( 0, 0, 1, 0 ) 
    962947               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    963948               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     
    980965      IF( ll_neg ) THEN 
    981966         DO jl = 1, jpl 
    982             DO_2D( 1, 0, 1, 0 ) 
     967            DO_2D( 0, 0, 1, 0 ) 
    983968               IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    984969                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     
    10481033      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    10491034         DO jl = 1, jpl 
    1050             DO_2D( 1, 0, 1, 0 ) 
     1035            DO_2D( 1, 0, 0, 0 ) 
    10511036               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
    10521037                  &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     
    10561041      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    10571042         DO jl = 1, jpl 
    1058             DO_2D( 1, 0, 1, 0 ) 
     1043            DO_2D( 1, 0, 0, 0 ) 
    10591044               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10601045               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
     
    10651050      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    10661051         DO jl = 1, jpl 
    1067             DO_2D( 1, 0, 1, 0 ) 
     1052            DO_2D( 1, 0, 0, 0 ) 
    10681053               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10691054               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     
    10781063      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    10791064         DO jl = 1, jpl 
    1080             DO_2D( 1, 0, 1, 0 ) 
     1065            DO_2D( 1, 0, 0, 0 ) 
    10811066               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10821067               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     
    10911076      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    10921077         DO jl = 1, jpl 
    1093             DO_2D( 1, 0, 1, 0 ) 
     1078            DO_2D( 1, 0, 0, 0 ) 
    10941079               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10951080               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     
    11121097      IF( ll_neg ) THEN 
    11131098         DO jl = 1, jpl 
    1114             DO_2D( 1, 0, 1, 0 ) 
     1099            DO_2D( 1, 0, 0, 0 ) 
    11151100               IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    11161101                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
     
    11221107      !                                                     !-- High order flux in j-direction  --! 
    11231108      DO jl = 1, jpl 
    1124          DO_2D( 1, 0, 1, 0 ) 
     1109         DO_2D( 1, 0, 0, 0 ) 
    11251110            pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 
    11261111         END_2D 
     
    11581143      ! -------------------------------------------------- 
    11591144      DO jl = 1, jpl 
    1160          DO_2D( 1, 0, 1, 0 ) 
     1145         DO_2D( 0, 0, 1, 0 ) 
    11611146            pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 
     1147         END_2D 
     1148         DO_2D( 1, 0, 0, 0 ) 
    11621149            pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 
    11631150         END_2D 
     
    12651252      ! --------------------------------- 
    12661253      DO jl = 1, jpl 
    1267          DO_2D( 1, 0, 1, 0 ) 
     1254         DO_2D( 0, 0, 1, 0 ) 
    12681255            zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 
    12691256            zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 
     
    12761263         END_2D 
    12771264 
    1278          DO_2D( 1, 0, 1, 0 ) 
     1265         DO_2D( 1, 0, 0, 0 ) 
    12791266            zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 
    12801267            zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 
     
    16331620   END SUBROUTINE Hsnow 
    16341621 
     1622   SUBROUTINE icemax3D( pice , pmax ) 
     1623      !!--------------------------------------------------------------------- 
     1624      !!                   ***  ROUTINE icemax3D ***                      
     1625      !! ** Purpose :  compute the max of the 9 points around 
     1626      !!---------------------------------------------------------------------- 
     1627      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input 
     1628      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output 
     1629      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1630      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1631      !!---------------------------------------------------------------------- 
     1632      DO jl = 1, jpl 
     1633         DO jj = Njs0-1, Nje0+1     
     1634            DO ji = Nis0, Nie0 
     1635               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 
     1636            END DO 
     1637         END DO 
     1638         DO jj = Njs0, Nje0     
     1639            DO ji = Nis0, Nie0 
     1640               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1641            END DO 
     1642         END DO 
     1643      END DO 
     1644   END SUBROUTINE icemax3D 
     1645 
     1646   SUBROUTINE icemax4D( pice , pmax ) 
     1647      !!--------------------------------------------------------------------- 
     1648      !!                   ***  ROUTINE icemax4D ***                      
     1649      !! ** Purpose :  compute the max of the 9 points around 
     1650      !!---------------------------------------------------------------------- 
     1651      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input 
     1652      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output 
     1653      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1654      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices 
     1655      !!---------------------------------------------------------------------- 
     1656      jlay = SIZE( pice , 3 )   ! size of input arrays 
     1657      DO jl = 1, jpl 
     1658         DO jk = 1, jlay 
     1659            DO jj = Njs0-1, Nje0+1     
     1660               DO ji = Nis0, Nie0 
     1661                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 
     1662               END DO 
     1663            END DO 
     1664            DO jj = Njs0, Nje0     
     1665               DO ji = Nis0, Nie0 
     1666                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1667               END DO 
     1668            END DO 
     1669         END DO 
     1670      END DO 
     1671   END SUBROUTINE icemax4D 
    16351672 
    16361673#else 
  • NEMO/trunk/src/ICE/icedyn_rdgrft.F90

    r13495 r13618  
    349349               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN 
    350350                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  & 
    351                      &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar ) 
     351                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar     ) * z1_gstar ) 
    352352               ELSE 
    353353                  apartf(ji,jl) = 0._wp 
  • NEMO/trunk/src/ICE/iceitd.F90

    r13472 r13618  
    627627         END_2D 
    628628         ! 
    629 !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 
    630          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 
    631          CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 
    632          ! 
    633          DO ji = 1, npti 
    634             jdonor(ji,jl)  = jl  
    635             ! how much of a_i you send in cat sup is somewhat arbitrary 
    636 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    637 !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
    638 !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
    639 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    640 !!          zdaice(ji,jl)  = a_i_1d(ji) 
    641 !!          zdvice(ji,jl)  = v_i_1d(ji) 
    642 !!clem: these are from UCL and work ok 
    643             zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
    644             zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
    645          END DO 
    646          ! 
    647          IF( npti > 0 ) THEN 
     629         IF( npti > 0 ) THEN             
     630            !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 
     631            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 
     632            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 
     633            ! 
     634            DO ji = 1, npti 
     635               jdonor(ji,jl)  = jl  
     636               ! how much of a_i you send in cat sup is somewhat arbitrary 
     637               !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
     638               !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
     639               !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
     640               !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
     641               !!          zdaice(ji,jl)  = a_i_1d(ji) 
     642               !!          zdvice(ji,jl)  = v_i_1d(ji) 
     643               !!clem: these are from UCL and work ok 
     644               zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
     645               zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     646            END DO 
     647            ! 
    648648            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1 
    649649            ! Reset shift parameters 
     
    666666         END_2D 
    667667         ! 
    668          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 
    669          CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 
    670          DO ji = 1, npti 
    671             jdonor(ji,jl) = jl + 1 
    672             zdaice(ji,jl) = a_i_1d(ji)  
    673             zdvice(ji,jl) = v_i_1d(ji) 
    674          END DO 
    675          ! 
    676668         IF( npti > 0 ) THEN 
     669            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 
     670            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 
     671            DO ji = 1, npti 
     672               jdonor(ji,jl) = jl + 1 
     673               zdaice(ji,jl) = a_i_1d(ji)  
     674               zdvice(ji,jl) = v_i_1d(ji) 
     675            END DO 
     676            ! 
    677677            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl 
    678678            ! Reset shift parameters 
Note: See TracChangeset for help on using the changeset viewer.