New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 13617 for NEMO – NEMO

Changeset 13617 for NEMO


Ignore:
Timestamp:
2020-10-16T10:07:20+02:00 (4 years ago)
Author:
clem
Message:

4.0-HEAD: optimization of both Prather and UMx advection schemes. Related to ticket #2552

Location:
NEMO/releases/r4.0/r4.0-HEAD/src/ICE
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_pra.F90

    r13589 r13617  
    111111      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
    112112      END WHERE 
    113       DO jl = 1, jpl 
    114          DO jj = 2, jpjm1 
    115             DO ji = fs_2, fs_jpim1 
    116                zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    117                   &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    118                   &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    119                   &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    120                zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    121                   &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    122                   &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    123                   &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    124                zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    125                   &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    126                   &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    127                   &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    128                zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
    129                   &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
    130                   &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
    131                   &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
    132             END DO 
    133          END DO 
    134       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 ) 
    135117      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1., zsi_max, 'T', 1. ) 
    136       ! 
     118 
    137119      ! enthalpies 
    138120      DO jk = 1, nlay_i 
    … …  
    145127         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
    146128         END WHERE 
    147       END DO 
    148       DO jl = 1, jpl 
    149          DO jk = 1, nlay_i 
    150             DO jj = 2, jpjm1 
    151                DO ji = fs_2, fs_jpim1 
    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 DO 
    157             END DO 
    158          END DO 
    159       END DO 
    160       DO jl = 1, jpl 
    161          DO jk = 1, nlay_s 
    162             DO jj = 2, jpjm1 
    163                DO ji = fs_2, fs_jpim1 
    164                   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), & 
    165                      &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
    166                      &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
    167                      &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
    168                END DO 
    169             END DO 
    170          END DO 
    171       END DO 
     129      END DO    
     130      CALL icemax4D( ze_i , zei_max ) 
     131      CALL icemax4D( ze_s , zes_max ) 
    172132      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
    173133      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
    … …  
    411371      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      - 
    412372      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      - 
     373      REAL(wp) ::   zpsm, zps0 
     374      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    413375      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace 
    414376      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      - 
    … …  
    429391             
    430392            DO ji = 1, jpi 
     393 
     394               zpsm  = psm (ji,jj,jl) ! optimization 
     395               zps0  = ps0 (ji,jj,jl) 
     396               zpsx  = psx (ji,jj,jl) 
     397               zpsxx = psxx(ji,jj,jl) 
     398               zpsy  = psy (ji,jj,jl) 
     399               zpsyy = psyy(ji,jj,jl) 
     400               zpsxy = psxy(ji,jj,jl) 
     401 
    431402               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
    432                psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 
    433                ! 
    434                zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     403               zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 
     404               ! 
     405               zslpmax = MAX( 0._wp, zps0 ) 
    435406               zs1max  = 1.5 * zslpmax 
    436                zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
    437                zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 
     407               zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) ) 
     408               zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 
    438409               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    439410 
    440                ps0 (ji,jj,jl) = zslpmax   
    441                psx (ji,jj,jl) = zs1new         * rswitch 
    442                psxx(ji,jj,jl) = zs2new         * rswitch 
    443                psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 
    444                psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 
    445                psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
     411               zps0 = zslpmax   
     412               zpsx  = zs1new  * rswitch 
     413               zpsxx = zs2new  * rswitch 
     414               zpsy  = zpsy    * rswitch 
     415               zpsyy = zpsyy  * rswitch 
     416               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
    446417 
    447418               !  Calculate fluxes and moments between boxes i<-->i+1               
    448419               !                                !  Flux from i to i+1 WHEN u GT 0  
    449420               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    450                zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
     421               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 
    451422               zalfq        =  zalf * zalf 
    452423               zalf1        =  1.0 - zalf 
    453424               zalf1q       =  zalf1 * zalf1 
    454425               ! 
    455                zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl) 
    456                zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 
    457                zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 
    458                zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq 
    459                zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    460                zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl) 
    461                zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl) 
     426               zfm (ji,jj)  =  zalf  *   zpsm  
     427               zf0 (ji,jj)  =  zalf  * ( zps0  + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 
     428               zfx (ji,jj)  =  zalfq * ( zpsx  + 3.0 * zalf1 * zpsxx ) 
     429               zfxx(ji,jj)  =  zalf  *   zpsxx * zalfq 
     430               zfy (ji,jj)  =  zalf  * ( zpsy  + zalf1 * zpsxy ) 
     431               zfxy(ji,jj)  =  zalfq *   zpsxy 
     432               zfyy(ji,jj)  =  zalf  *   zpsyy 
    462433 
    463434               !                                !  Readjust moments remaining in the box. 
    464                psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    465                ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    466                psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 
    467                psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl) 
    468                psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj) 
    469                psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj) 
    470                psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     435               zpsm  =  zpsm  - zfm(ji,jj) 
     436               zps0  =  zps0  - zf0(ji,jj) 
     437               zpsx  =  zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 
     438               zpsxx =  zalf1  * zalf1q * zpsxx 
     439               zpsy  =  zpsy  - zfy (ji,jj) 
     440               zpsyy =  zpsyy - zfyy(ji,jj) 
     441               zpsxy =  zalf1q * zpsxy 
     442               ! 
     443               psm (ji,jj,jl) = zpsm ! optimization 
     444               ps0 (ji,jj,jl) = zps0  
     445               psx (ji,jj,jl) = zpsx  
     446               psxx(ji,jj,jl) = zpsxx 
     447               psy (ji,jj,jl) = zpsy  
     448               psyy(ji,jj,jl) = zpsyy 
     449               psxy(ji,jj,jl) = zpsxy 
     450               ! 
    471451            END DO 
    472452 
    … …  
    492472 
    493473            DO ji = fs_2, fs_jpim1  
     474               ! 
     475               zpsm  = psm (ji,jj,jl) ! optimization 
     476               zps0  = ps0 (ji,jj,jl) 
     477               zpsx  = psx (ji,jj,jl) 
     478               zpsxx = psxx(ji,jj,jl) 
     479               zpsy  = psy (ji,jj,jl) 
     480               zpsyy = psyy(ji,jj,jl) 
     481               zpsxy = psxy(ji,jj,jl) 
     482               !                                !  Readjust moments remaining in the box. 
    494483               zbt  =       zbet(ji-1,jj) 
    495484               zbt1 = 1.0 - zbet(ji-1,jj) 
    496                !                                !  Readjust moments remaining in the box. 
    497                psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 
    498                ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 
    499                psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 
    500                psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 
    501                psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 
    502                psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 
    503                psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 
     485               ! 
     486               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 
     487               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 
     488               zpsx  = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 
     489               zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 
     490               zpsy  = zbt * zpsy  + zbt1 * ( zpsy - zfy (ji-1,jj) ) 
     491               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 
     492               zpsxy = zalg1q(ji-1,jj) * zpsxy 
    504493 
    505494               !   Put the temporary moments into appropriate neighboring boxes.     
    506495               !                                !   Flux from i to i+1 IF u GT 0. 
    507                zbt  =       zbet(ji-1,jj) 
    508                zbt1 = 1.0 - zbet(ji-1,jj) 
    509                psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 
    510                zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 
    511                zalf1         = 1.0 - zalf 
    512                ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 
    513                ! 
    514                ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 
    515                psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 
    516                psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             & 
    517                   &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 
    518                   &            + zbt1 * psxx(ji,jj,jl) 
    519                psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             & 
    520                   &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   & 
    521                   &            + zbt1 * psxy(ji,jj,jl) 
    522                psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 
    523                psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 
     496               zbt   =       zbet(ji-1,jj) 
     497               zbt1  = 1.0 - zbet(ji-1,jj) 
     498               zpsm  = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 
     499               zalf  = zbt * zfm(ji-1,jj) / zpsm 
     500               zalf1 = 1.0 - zalf 
     501               ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 
     502               ! 
     503               zps0  =  zbt  * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 
     504               zpsx  =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 
     505               zpsxx =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx                            & 
     506                  &            + 5.0 * ( zalf * zalf1 * ( zpsx  - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 
     507                  &            + zbt1 * zpsxx 
     508               zpsxy =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy            & 
     509                  &            + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * zpsy ) )  & 
     510                  &            + zbt1 * zpsxy 
     511               zpsy  =  zbt  * ( zpsy  + zfy (ji-1,jj) ) + zbt1 * zpsy  
     512               zpsyy =  zbt  * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 
    524513 
    525514               !                                !  Flux from i+1 to i IF u LT 0. 
    526                zbt  =       zbet(ji,jj) 
    527                zbt1 = 1.0 - zbet(ji,jj) 
    528                psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    529                zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    530                zalf1         = 1.0 - zalf 
    531                ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    532                ! 
    533                ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 
    534                psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 
    535                psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 
    536                   &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    & 
    537                   &                                           + ( zalf1 - zalf ) * ztemp ) ) 
    538                psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    539                   &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 
    540                psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 
    541                psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 
     515               zbt   =       zbet(ji,jj) 
     516               zbt1  = 1.0 - zbet(ji,jj) 
     517               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     518               zalf  = zbt1 * zfm(ji,jj) / zpsm 
     519               zalf1 = 1.0 - zalf 
     520               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     521               ! 
     522               zps0  = zbt * zps0  + zbt1 * ( zps0 + zf0(ji,jj) ) 
     523               zpsx  = zbt * zpsx  + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 
     524               zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 
     525                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) )    & 
     526                  &                         + ( zalf1 - zalf ) * ztemp ) ) 
     527               zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     528                  &                         + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 
     529               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  + zfy (ji,jj) ) 
     530               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 
     531               ! 
     532               psm (ji,jj,jl) = zpsm  ! optimization 
     533               ps0 (ji,jj,jl) = zps0  
     534               psx (ji,jj,jl) = zpsx  
     535               psxx(ji,jj,jl) = zpsxx 
     536               psy (ji,jj,jl) = zpsy  
     537               psyy(ji,jj,jl) = zpsyy 
     538               psxy(ji,jj,jl) = zpsxy 
    542539            END DO 
    543540             
    … …  
    570567      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
    571568      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     569      REAL(wp) ::   zpsm, zps0 
     570      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    572571      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
    573572      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
    … …  
    587586         DO jj = 1, jpj 
    588587            DO ji = jimin, jimax 
     588               ! 
     589               zpsm  = psm (ji,jj,jl) ! optimization 
     590               zps0  = ps0 (ji,jj,jl) 
     591               zpsx  = psx (ji,jj,jl) 
     592               zpsxx = psxx(ji,jj,jl) 
     593               zpsy  = psy (ji,jj,jl) 
     594               zpsyy = psyy(ji,jj,jl) 
     595               zpsxy = psxy(ji,jj,jl) 
     596               ! 
    589597               !  Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 
    590                psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  ) 
    591                ! 
    592                zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     598               zpsm = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20  ) 
     599               ! 
     600               zslpmax = MAX( 0._wp, zps0 ) 
    593601               zs1max  = 1.5 * zslpmax 
    594                zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
    595                zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 
     602               zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) ) 
     603               zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 
    596604               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    597605               ! 
    598                ps0 (ji,jj,jl) = zslpmax   
    599                psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 
    600                psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 
    601                psy (ji,jj,jl) = zs1new         * rswitch 
    602                psyy(ji,jj,jl) = zs2new         * rswitch 
    603                psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
     606               zps0 = zslpmax   
     607               zpsx  = zpsx * rswitch 
     608               zpsxx = zpsxx * rswitch 
     609               zpsy = zs1new         * rswitch 
     610               zpsyy = zs2new         * rswitch 
     611               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
    604612  
    605613               !  Calculate fluxes and moments between boxes j<-->j+1               
    606614               !                                !  Flux from j to j+1 WHEN v GT 0    
    607615               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    608                zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
     616               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 
    609617               zalfq        =  zalf * zalf 
    610618               zalf1        =  1.0 - zalf 
    611619               zalf1q       =  zalf1 * zalf1 
    612620               ! 
    613                zfm (ji,jj)  =  zalf  * psm(ji,jj,jl) 
    614                zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) )  
    615                zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 
    616                zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl) 
    617                zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    618                zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl) 
    619                zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl) 
     621               zfm (ji,jj)  =  zalf  * zpsm 
     622               zf0 (ji,jj)  =  zalf  * ( zps0 + zalf1 * ( zpsy  + (zalf1-zalf) * zpsyy ) )  
     623               zfy (ji,jj)  =  zalfq *( zpsy + 3.0*zalf1*zpsyy ) 
     624               zfyy(ji,jj)  =  zalf  * zalfq * zpsyy 
     625               zfx (ji,jj)  =  zalf  * ( zpsx + zalf1 * zpsxy ) 
     626               zfxy(ji,jj)  =  zalfq * zpsxy 
     627               zfxx(ji,jj)  =  zalf  * zpsxx 
    620628               ! 
    621629               !                                !  Readjust moments remaining in the box. 
    622                psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    623                ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    624                psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 
    625                psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl) 
    626                psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj) 
    627                psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj) 
    628                psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     630               zpsm   =  zpsm  - zfm(ji,jj) 
     631               zps0   =  zps0  - zf0(ji,jj) 
     632               zpsy   =  zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 
     633               zpsyy  =  zalf1 * zalf1q * zpsyy 
     634               zpsx   =  zpsx  - zfx(ji,jj) 
     635               zpsxx  =  zpsxx - zfxx(ji,jj) 
     636               zpsxy  =  zalf1q * zpsxy 
     637               ! 
     638               psm (ji,jj,jl) = zpsm ! optimization 
     639               ps0 (ji,jj,jl) = zps0  
     640               psx (ji,jj,jl) = zpsx  
     641               psxx(ji,jj,jl) = zpsxx 
     642               psy (ji,jj,jl) = zpsy  
     643               psyy(ji,jj,jl) = zpsyy 
     644               psxy(ji,jj,jl) = zpsxy 
    629645            END DO 
    630646         END DO 
    … …  
    658674               zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    659675               ! 
    660                psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 
    661                ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 
    662                psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 
    663                psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 
    664                psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 
    665                psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 
    666                psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 
     676               zpsm  = psm (ji,jj,jl) ! optimization 
     677               zps0  = ps0 (ji,jj,jl) 
     678               zpsx  = psx (ji,jj,jl) 
     679               zpsxx = psxx(ji,jj,jl) 
     680               zpsy  = psy (ji,jj,jl) 
     681               zpsyy = psyy(ji,jj,jl) 
     682               zpsxy = psxy(ji,jj,jl) 
     683               ! 
     684               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 
     685               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 
     686               zpsy  = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 
     687               zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 
     688               zpsx  = zbt * zpsx  + zbt1 * ( zpsx  - zfx (ji,jj-1) ) 
     689               zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 
     690               zpsxy = zalg1q(ji,jj-1) * zpsxy 
    667691 
    668692               !   Put the temporary moments into appropriate neighboring boxes.     
    669693               !                                !   Flux from j to j+1 IF v GT 0. 
    670                zbt  =       zbet(ji,jj-1) 
    671                zbt1 = 1.0 - zbet(ji,jj-1) 
    672                psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)  
    673                zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)  
    674                zalf1         = 1.0 - zalf 
    675                ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 
    676                ! 
    677                ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 
    678                psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  & 
    679                   &             + zbt1 * psy(ji,jj,jl)   
    680                psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           & 
    681                   &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
    682                   &             + zbt1 * psyy(ji,jj,jl) 
    683                psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            & 
    684                   &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  & 
    685                   &             + zbt1 * psxy(ji,jj,jl) 
    686                psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 
    687                psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 
     694               zbt   =       zbet(ji,jj-1) 
     695               zbt1  = 1.0 - zbet(ji,jj-1) 
     696               zpsm  = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm  
     697               zalf  = zbt * zfm(ji,jj-1) / zpsm  
     698               zalf1 = 1.0 - zalf 
     699               ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 
     700               ! 
     701               zps0  =   zbt  * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 
     702               zpsy  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp )  & 
     703                  &             + zbt1 * zpsy   
     704               zpsyy =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy                           & 
     705                  &             + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
     706                  &             + zbt1 * zpsyy 
     707               zpsxy =   zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy             & 
     708                  &             + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )  & 
     709                  &             + zbt1 * zpsxy 
     710               zpsx  =   zbt * ( zpsx  + zfx (ji,jj-1) ) + zbt1 * zpsx  
     711               zpsxx =   zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 
    688712 
    689713               !                                !  Flux from j+1 to j IF v LT 0. 
    690                zbt  =       zbet(ji,jj) 
    691                zbt1 = 1.0 - zbet(ji,jj) 
    692                psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    693                zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    694                zalf1         = 1.0 - zalf 
    695                ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    696                ! 
    697                ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) ) 
    698                psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 
    699                psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 
    700                   &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    & 
    701                   &                                            + ( zalf1 - zalf ) * ztemp ) ) 
    702                psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    703                   &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 
    704                psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 
    705                psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 
     714               zbt   =       zbet(ji,jj) 
     715               zbt1  = 1.0 - zbet(ji,jj) 
     716               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     717               zalf  = zbt1 * zfm(ji,jj) / zpsm 
     718               zalf1 = 1.0 - zalf 
     719               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     720               ! 
     721               zps0  = zbt * zps0  + zbt1 * (  zps0 + zf0(ji,jj) ) 
     722               zpsy  = zbt * zpsy  + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 
     723               zpsyy = zbt * zpsyy + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 
     724                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) )     & 
     725                  &                         + ( zalf1 - zalf ) * ztemp ) ) 
     726               zpsxy = zbt * zpsxy + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     727                  &                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 
     728               zpsx  = zbt * zpsx  + zbt1 * ( zpsx  + zfx (ji,jj) ) 
     729               zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 
     730               ! 
     731               psm (ji,jj,jl) = zpsm ! optimization 
     732               ps0 (ji,jj,jl) = zps0  
     733               psx (ji,jj,jl) = zpsx  
     734               psxx(ji,jj,jl) = zpsxx 
     735               psy (ji,jj,jl) = zpsy  
     736               psyy(ji,jj,jl) = zpsyy 
     737               psxy(ji,jj,jl) = zpsxy 
    706738            END DO 
    707739         END DO 
    … …  
    11311163   END SUBROUTINE adv_pra_rst 
    11321164 
     1165   SUBROUTINE icemax3D( pice , pmax ) 
     1166      !!--------------------------------------------------------------------- 
     1167      !!                   ***  ROUTINE icemax3D ***                      
     1168      !! ** Purpose :  compute the max of the 9 points around 
     1169      !!---------------------------------------------------------------------- 
     1170      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input 
     1171      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output 
     1172      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1173      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1174      !!---------------------------------------------------------------------- 
     1175      DO jl = 1, jpl 
     1176         DO jj = 1, jpj 
     1177            DO ji = 2, jpim1 
     1178               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 
     1179            END DO 
     1180         END DO 
     1181         DO jj = 2, jpjm1 
     1182            DO ji = 2, jpim1 
     1183               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1184            END DO 
     1185         END DO 
     1186      END DO 
     1187   END SUBROUTINE icemax3D 
     1188 
     1189   SUBROUTINE icemax4D( pice , pmax ) 
     1190      !!--------------------------------------------------------------------- 
     1191      !!                   ***  ROUTINE icemax4D ***                      
     1192      !! ** Purpose :  compute the max of the 9 points around 
     1193      !!---------------------------------------------------------------------- 
     1194      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input 
     1195      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output 
     1196      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1197      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices 
     1198      !!---------------------------------------------------------------------- 
     1199      jlay = SIZE( pice , 3 )   ! size of input arrays 
     1200      DO jl = 1, jpl 
     1201         DO jk = 1, jlay 
     1202            DO jj = 1, jpj 
     1203               DO ji = 2, jpim1 
     1204                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 
     1205               END DO 
     1206            END DO 
     1207            DO jj = 2, jpjm1 
     1208               DO ji = 2, jpim1 
     1209                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1210               END DO 
     1211            END DO 
     1212         END DO 
     1213      END DO 
     1214   END SUBROUTINE icemax4D 
     1215    
    11331216#else 
    11341217   !!---------------------------------------------------------------------- 

  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_umx.F90

    r13589 r13617  
    115115      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
    116116      END WHERE 
    117       DO jl = 1, jpl 
    118          DO jj = 2, jpjm1 
    119             DO ji = fs_2, fs_jpim1 
    120                zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    121                   &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    122                   &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    123                   &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    124                zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    125                   &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    126                   &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    127                   &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    128                zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    129                   &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    130                   &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    131                   &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    132                zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
    133                   &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
    134                   &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
    135                   &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
    136             END DO 
    137          END DO 
    138       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 ) 
    139121      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1., zsi_max, 'T', 1. ) 
    140       ! 
     122 
    141123      ! enthalpies 
    142124      DO jk = 1, nlay_i 
    … …  
    149131         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
    150132         END WHERE 
    151       END DO 
    152       DO jl = 1, jpl 
    153          DO jk = 1, nlay_i 
    154             DO jj = 2, jpjm1 
    155                DO ji = fs_2, fs_jpim1 
    156                   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), & 
    157                      &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
    158                      &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
    159                      &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
    160                END DO 
    161             END DO 
    162          END DO 
    163       END DO 
    164       DO jl = 1, jpl 
    165          DO jk = 1, nlay_s 
    166             DO jj = 2, jpjm1 
    167                DO ji = fs_2, fs_jpim1 
    168                   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), & 
    169                      &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
    170                      &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
    171                      &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
    172                END DO 
    173             END DO 
    174          END DO 
    175       END DO 
     133      END DO    
     134      CALL icemax4D( ze_i , zei_max ) 
     135      CALL icemax4D( ze_s , zes_max ) 
    176136      CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1. ) 
    177137      CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1. ) 
    178       ! 
    179138      ! 
    180139      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
    … …  
    17861745   END SUBROUTINE Hsnow 
    17871746 
     1747   SUBROUTINE icemax3D( pice , pmax ) 
     1748      !!--------------------------------------------------------------------- 
     1749      !!                   ***  ROUTINE icemax3D ***                      
     1750      !! ** Purpose :  compute the max of the 9 points around 
     1751      !!---------------------------------------------------------------------- 
     1752      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input 
     1753      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output 
     1754      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1755      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1756      !!---------------------------------------------------------------------- 
     1757      DO jl = 1, jpl 
     1758         DO jj = 1, jpj 
     1759            DO ji = 2, jpim1 
     1760               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 
     1761            END DO 
     1762         END DO 
     1763         DO jj = 2, jpjm1 
     1764            DO ji = 2, jpim1 
     1765               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1766            END DO 
     1767         END DO 
     1768      END DO 
     1769   END SUBROUTINE icemax3D 
     1770 
     1771   SUBROUTINE icemax4D( pice , pmax ) 
     1772      !!--------------------------------------------------------------------- 
     1773      !!                   ***  ROUTINE icemax4D ***                      
     1774      !! ** Purpose :  compute the max of the 9 points around 
     1775      !!---------------------------------------------------------------------- 
     1776      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input 
     1777      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output 
     1778      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1779      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices 
     1780      !!---------------------------------------------------------------------- 
     1781      jlay = SIZE( pice , 3 )   ! size of input arrays 
     1782      DO jl = 1, jpl 
     1783         DO jk = 1, jlay 
     1784            DO jj = 1, jpj 
     1785               DO ji = 2, jpim1 
     1786                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 
     1787               END DO 
     1788            END DO 
     1789            DO jj = 2, jpjm1 
     1790               DO ji = 2, jpim1 
     1791                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1792               END DO 
     1793            END DO 
     1794         END DO 
     1795      END DO 
     1796   END SUBROUTINE icemax4D 
    17881797 
    17891798#else 

  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_rdgrft.F90

    r13494 r13617  
    341341               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN 
    342342                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  & 
    343                      &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar ) 
     343                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar     ) * z1_gstar ) 
    344344               ELSE 
    345345                  apartf(ji,jl) = 0._wp 

  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceitd.F90

    r13284 r13617  
    615615         END DO 
    616616         ! 
    617 !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 
    618          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 
    619          CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 
    620          ! 
    621          DO ji = 1, npti 
    622             jdonor(ji,jl)  = jl  
    623             ! how much of a_i you send in cat sup is somewhat arbitrary 
    624 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    625 !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
    626 !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
    627 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    628 !!          zdaice(ji,jl)  = a_i_1d(ji) 
    629 !!          zdvice(ji,jl)  = v_i_1d(ji) 
    630 !!clem: these are from UCL and work ok 
    631             zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
    632             zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
    633          END DO 
    634          ! 
    635          IF( npti > 0 ) THEN 
     617         IF( npti > 0 ) THEN             
     618            !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 
     619            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 
     620            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 
     621            ! 
     622            DO ji = 1, npti 
     623               jdonor(ji,jl)  = jl  
     624               ! how much of a_i you send in cat sup is somewhat arbitrary 
     625               !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
     626               !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
     627               !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
     628               !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
     629               !!          zdaice(ji,jl)  = a_i_1d(ji) 
     630               !!          zdvice(ji,jl)  = v_i_1d(ji) 
     631               !!clem: these are from UCL and work ok 
     632               zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
     633               zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     634            END DO 
     635            ! 
    636636            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1 
    637637            ! Reset shift parameters 
    … …  
    656656         END DO 
    657657         ! 
    658          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 
    659          CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 
    660          DO ji = 1, npti 
    661             jdonor(ji,jl) = jl + 1 
    662             zdaice(ji,jl) = a_i_1d(ji)  
    663             zdvice(ji,jl) = v_i_1d(ji) 
    664          END DO 
    665          ! 
    666658         IF( npti > 0 ) THEN 
     659            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 
     660            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 
     661            ! 
     662            DO ji = 1, npti 
     663               jdonor(ji,jl) = jl + 1 
     664               zdaice(ji,jl) = a_i_1d(ji)  
     665               zdvice(ji,jl) = v_i_1d(ji) 
     666            END DO 
     667            ! 
    667668            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl 
    668669            ! Reset shift parameters 
Note: See TracChangeset for help on using the changeset viewer.