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 13694 for NEMO/branches/2020/r12377_ticket2386/src/ICE/icedyn_adv_umx.F90 – NEMO

Ignore:
Timestamp:
2020-10-28T18:03:31+01:00 (3 years ago)
Author:
andmirek
Message:

Ticket #2386: merge with trunk rev 13688

Location:
NEMO/branches/2020/r12377_ticket2386
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/r12377_ticket2386

    • Property svn:externals
      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@13507        sette 
         10^/utils/CI/sette@13559        sette 
  • NEMO/branches/2020/r12377_ticket2386/src/ICE/icedyn_adv_umx.F90

    r13540 r13694  
    9292      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    9393      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    94       REAL(wp) ::   zdt, zvi_cen 
     94      REAL(wp) ::   zdt, z1_dt, zvi_cen 
    9595      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
    9696      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx, zcu_box, zcv_box 
     
    104104      ! 
    105105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs  
     106      !! diagnostics 
     107      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat       
    106108      !!---------------------------------------------------------------------- 
    107109      ! 
     
    113115      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
    114116      END WHERE 
    115       DO jl = 1, jpl 
    116          DO_2D( 0, 0, 0, 0 ) 
    117             zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    118                &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    119                &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    120                &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    121             zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    122                &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    123                &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    124                &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    125             zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    126                &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    127                &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    128                &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    129             zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
    130                &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
    131                &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
    132                &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
    133          END_2D 
    134       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 ) 
    135121      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 ) 
    136122      ! 
     
    145131         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
    146132         END WHERE 
    147       END DO 
    148       DO jl = 1, jpl 
    149          DO_3D( 0, 0, 0, 0, 1, nlay_i ) 
    150             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), & 
    151                &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
    152                &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
    153                &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
    154          END_3D 
    155       END DO 
    156       DO jl = 1, jpl 
    157          DO_3D( 0, 0, 0, 0, 1, nlay_s ) 
    158             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), & 
    159                &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
    160                &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
    161                &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
    162          END_3D 
    163       END DO 
    164       CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
    165       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 ) 
    166138      ! 
    167139      ! 
     
    179151      ENDIF 
    180152      zdt = rDt_ice / REAL(icycle) 
     153      z1_dt = 1._wp / zdt 
    181154 
    182155      ! --- transport --- ! 
     
    207180      !---------------! 
    208181      DO jt = 1, icycle 
     182 
     183         ! diagnostics 
     184         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     185         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
     186         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     187            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 
    209188 
    210189         ! record at_i before advection (for open water) 
     
    377356            ENDIF 
    378357         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 ) 
    379371         ! 
    380372         !== Open water area ==! 
     
    384376               &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
    385377         END_2D 
    386          CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1.0_wp ) 
    387          ! 
     378         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1._wp ) 
     379         ! 
     380         ! --- diagnostics --- ! 
     381         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     382            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
     383         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
     384            &                                        - zdiag_adv_salt(:,:) ) * z1_dt 
     385         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     386            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 
     387            &                                        - zdiag_adv_heat(:,:) ) * z1_dt 
    388388         ! 
    389389         ! --- Ensure non-negative fields and in-bound thicknesses --- ! 
     
    445445      !!             work on H (and not V). It is partly related to the multi-category approach 
    446446      !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 
    447       !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 
    448       !!             since sv_i and e_i are still good. 
     447      !!             concentration is small). We also limit S and T. 
    449448      !!---------------------------------------------------------------------- 
    450449      REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     
    490489      IF( pamsk == 0._wp ) THEN 
    491490         DO jl = 1, jpl 
    492             DO_2D( 1, 0, 1, 0 ) 
     491            DO_2D( 0, 0, 1, 0 ) 
    493492               IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 
    494493                  zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj) 
     
    499498               ENDIF 
    500499               ! 
     500            END_2D 
     501            DO_2D( 1, 0, 0, 0 ) 
    501502               IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 
    502503                  zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj) 
     
    533534      IF( PRESENT( pua_ho ) ) THEN 
    534535         DO jl = 1, jpl 
    535             DO_2D( 1, 0, 1, 0 ) 
    536                pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
    537                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) 
    538543            END_2D 
    539544         END DO 
     
    549554         END_2D 
    550555      END DO 
    551       CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1.0_wp ) 
    552556      ! 
    553557   END SUBROUTINE adv_umx 
     
    588592            ! 
    589593            DO jl = 1, jpl              !-- flux in x-direction 
    590                DO_2D( 1, 0, 1, 0 ) 
     594               DO_2D( 1, 1, 1, 0 ) 
    591595                  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) 
    592596               END_2D 
     
    594598            ! 
    595599            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    596                DO_2D( 0, 0, 0, 0 ) 
     600               DO_2D( 1, 1, 0, 0 ) 
    597601                  ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              & 
    598602                     &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     
    601605               END_2D 
    602606            END DO 
    603             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    604607            ! 
    605608            DO jl = 1, jpl              !-- flux in y-direction 
    606                DO_2D( 1, 0, 1, 0 ) 
     609               DO_2D( 1, 0, 0, 0 ) 
    607610                  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) 
    608611               END_2D 
     
    612615            ! 
    613616            DO jl = 1, jpl              !-- flux in y-direction 
    614                DO_2D( 1, 0, 1, 0 ) 
     617               DO_2D( 1, 0, 1, 1 ) 
    615618                  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) 
    616619               END_2D 
     
    618621            ! 
    619622            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    620                DO_2D( 0, 0, 0, 0 ) 
     623               DO_2D( 0, 0, 1, 1 ) 
    621624                  ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  & 
    622625                     &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     
    625628               END_2D 
    626629            END DO 
    627             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    628630            ! 
    629631            DO jl = 1, jpl              !-- flux in x-direction 
    630                DO_2D( 1, 0, 1, 0 ) 
     632               DO_2D( 0, 0, 1, 0 ) 
    631633                  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) 
    632634               END_2D 
     
    677679         ! 
    678680         DO jl = 1, jpl 
    679             DO_2D( 1, 0, 1, 0 ) 
     681            DO_2D( 1, 1, 1, 0 ) 
    680682               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 ) 
    681685               pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) ) 
    682686            END_2D 
     
    695699            ! 
    696700            DO jl = 1, jpl              !-- flux in x-direction 
    697                DO_2D( 1, 0, 1, 0 ) 
     701               DO_2D( 1, 1, 1, 0 ) 
    698702                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
    699703               END_2D 
     
    702706 
    703707            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    704                DO_2D( 0, 0, 0, 0 ) 
     708               DO_2D( 1, 1, 0, 0 ) 
    705709                  ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              & 
    706710                     &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     
    709713               END_2D 
    710714            END DO 
    711             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    712715 
    713716            DO jl = 1, jpl              !-- flux in y-direction 
    714                DO_2D( 1, 0, 1, 0 ) 
     717               DO_2D( 1, 0, 0, 0 ) 
    715718                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 
    716719               END_2D 
     
    721724            ! 
    722725            DO jl = 1, jpl              !-- flux in y-direction 
    723                DO_2D( 1, 0, 1, 0 ) 
     726               DO_2D( 1, 0, 1, 1 ) 
    724727                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
    725728               END_2D 
     
    728731            ! 
    729732            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    730                DO_2D( 0, 0, 0, 0 ) 
     733               DO_2D( 0, 0, 1, 1 ) 
    731734                  ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  & 
    732735                     &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     
    735738               END_2D 
    736739            END DO 
    737             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    738740            ! 
    739741            DO jl = 1, jpl              !-- flux in x-direction 
    740                DO_2D( 1, 0, 1, 0 ) 
     742               DO_2D( 0, 0, 1, 0 ) 
    741743                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 
    742744               END_2D 
     
    895897         !         
    896898         DO jl = 1, jpl 
    897             DO_2D( 1, 0, 1, 0 ) 
     899            DO_2D( 0, 0, 1, 0 ) 
    898900               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    899901                  &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     
    904906         ! 
    905907         DO jl = 1, jpl 
    906             DO_2D( 1, 0, 1, 0 ) 
     908            DO_2D( 0, 0, 1, 0 ) 
    907909               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    908910               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     
    914916         ! 
    915917         DO jl = 1, jpl 
    916             DO_2D( 1, 0, 1, 0 ) 
     918            DO_2D( 0, 0, 1, 0 ) 
    917919               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    918920               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     
    928930         ! 
    929931         DO jl = 1, jpl 
    930             DO_2D( 1, 0, 1, 0 ) 
     932            DO_2D( 0, 0, 1, 0 ) 
    931933               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    932934               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     
    942944         ! 
    943945         DO jl = 1, jpl 
    944             DO_2D( 1, 0, 1, 0 ) 
     946            DO_2D( 0, 0, 1, 0 ) 
    945947               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    946948               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     
    963965      IF( ll_neg ) THEN 
    964966         DO jl = 1, jpl 
    965             DO_2D( 1, 0, 1, 0 ) 
     967            DO_2D( 0, 0, 1, 0 ) 
    966968               IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    967969                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     
    973975      !                                                     !-- High order flux in i-direction  --! 
    974976      DO jl = 1, jpl 
    975          DO_2D( 1, 0, 1, 0 ) 
     977         DO_2D( 0, 0, 1, 0 ) 
    976978            pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 
    977979         END_2D 
     
    10311033      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    10321034         DO jl = 1, jpl 
    1033             DO_2D( 1, 0, 1, 0 ) 
     1035            DO_2D( 1, 0, 0, 0 ) 
    10341036               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
    10351037                  &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     
    10391041      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    10401042         DO jl = 1, jpl 
    1041             DO_2D( 1, 0, 1, 0 ) 
     1043            DO_2D( 1, 0, 0, 0 ) 
    10421044               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10431045               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
     
    10481050      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    10491051         DO jl = 1, jpl 
    1050             DO_2D( 1, 0, 1, 0 ) 
     1052            DO_2D( 1, 0, 0, 0 ) 
    10511053               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10521054               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     
    10611063      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    10621064         DO jl = 1, jpl 
    1063             DO_2D( 1, 0, 1, 0 ) 
     1065            DO_2D( 1, 0, 0, 0 ) 
    10641066               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10651067               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     
    10741076      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    10751077         DO jl = 1, jpl 
    1076             DO_2D( 1, 0, 1, 0 ) 
     1078            DO_2D( 1, 0, 0, 0 ) 
    10771079               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10781080               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     
    10951097      IF( ll_neg ) THEN 
    10961098         DO jl = 1, jpl 
    1097             DO_2D( 1, 0, 1, 0 ) 
     1099            DO_2D( 1, 0, 0, 0 ) 
    10981100               IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    10991101                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
     
    11051107      !                                                     !-- High order flux in j-direction  --! 
    11061108      DO jl = 1, jpl 
    1107          DO_2D( 1, 0, 1, 0 ) 
     1109         DO_2D( 1, 0, 0, 0 ) 
    11081110            pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 
    11091111         END_2D 
     
    11411143      ! -------------------------------------------------- 
    11421144      DO jl = 1, jpl 
    1143          DO_2D( 1, 0, 1, 0 ) 
     1145         DO_2D( 0, 0, 1, 0 ) 
    11441146            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 ) 
    11451149            pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 
    11461150         END_2D 
     
    12481252      ! --------------------------------- 
    12491253      DO jl = 1, jpl 
    1250          DO_2D( 1, 0, 1, 0 ) 
     1254         DO_2D( 0, 0, 1, 0 ) 
    12511255            zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 
    12521256            zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 
     
    12591263         END_2D 
    12601264 
    1261          DO_2D( 1, 0, 1, 0 ) 
     1265         DO_2D( 1, 0, 0, 0 ) 
    12621266            zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 
    12631267            zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 
     
    16161620   END SUBROUTINE Hsnow 
    16171621 
     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 
    16181672 
    16191673#else 
Note: See TracChangeset for help on using the changeset viewer.