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 9019 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90 – NEMO

Ignore:
Timestamp:
2017-12-13T15:58:53+01:00 (6 years ago)
Author:
timgraham
Message:

Merge of dev_CNRS_2017 into branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90

    r7761 r9019  
    2828   PRIVATE 
    2929 
    30    PUBLIC agrif_interp_lim3 
     30   PUBLIC   agrif_interp_lim3   ! called by agrif_user.F90 
    3131 
    3232   !!---------------------------------------------------------------------- 
     
    4646      !!  computing factor for time interpolation 
    4747      !!----------------------------------------------------------------------- 
    48       CHARACTER(len=1), INTENT( in )           :: cd_type 
    49       INTEGER         , INTENT( in ), OPTIONAL :: kiter, kitermax 
    50       !! 
    51       REAL(wp) :: zbeta 
    52       !!----------------------------------------------------------------------- 
    53       ! 
    54       IF( Agrif_Root() )  RETURN 
    55       ! 
    56       SELECT CASE(cd_type) 
     48      CHARACTER(len=1), INTENT(in   )           ::  cd_type 
     49      INTEGER         , INTENT(in   ), OPTIONAL ::  kiter, kitermax 
     50      !! 
     51      REAL(wp) ::   zbeta   ! local scalar 
     52      !!----------------------------------------------------------------------- 
     53      ! 
     54      IF( Agrif_Root() .OR. nn_ice==0 )  RETURN   ! do not interpolate if inside Parent domain or if child domain does not have ice 
     55      ! 
     56      SELECT CASE( cd_type ) 
    5757      CASE('U','V') 
    5858         IF( PRESENT( kiter ) ) THEN  ! interpolation at the child sub-time step (only for ice rheology) 
     
    6666      END SELECT 
    6767      ! 
    68       Agrif_SpecialValue=-9999. 
     68      Agrif_SpecialValue    = -9999. 
    6969      Agrif_UseSpecialValue = .TRUE. 
    70       SELECT CASE(cd_type) 
    71       CASE('U') 
    72          CALL Agrif_Bc_variable( u_ice_id  , procname=interp_u_ice  , calledweight=zbeta ) 
    73       CASE('V') 
    74          CALL Agrif_Bc_variable( v_ice_id  , procname=interp_v_ice  , calledweight=zbeta ) 
    75       CASE('T') 
    76          CALL Agrif_Bc_variable( tra_ice_id, procname=interp_tra_ice, calledweight=zbeta ) 
     70      SELECT CASE( cd_type ) 
     71      CASE('U')   ;   CALL Agrif_Bc_variable( u_ice_id  , procname=interp_u_ice  , calledweight=zbeta ) 
     72      CASE('V')   ;   CALL Agrif_Bc_variable( v_ice_id  , procname=interp_v_ice  , calledweight=zbeta ) 
     73      CASE('T')   ;   CALL Agrif_Bc_variable( tra_ice_id, procname=interp_tra_ice, calledweight=zbeta ) 
    7774      END SELECT 
    78       Agrif_SpecialValue=0. 
     75      Agrif_SpecialValue    = 0._wp 
    7976      Agrif_UseSpecialValue = .FALSE. 
    8077      ! 
    8178   END SUBROUTINE agrif_interp_lim3 
    8279 
    83    !!------------------ 
    84    !! Local subroutines 
    85    !!------------------ 
     80 
    8681   SUBROUTINE interp_u_ice( ptab, i1, i2, j1, j2, before ) 
    8782      !!----------------------------------------------------------------------- 
     
    8984      !! 
    9085      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after) 
    91       !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points, 
    92       !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around 
    93       !!----------------------------------------------------------------------- 
    94       INTEGER , INTENT(in) :: i1, i2, j1, j2 
    95       REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 
    96       LOGICAL , INTENT(in) :: before 
    97       !! 
    98       REAL(wp) :: zrhoy 
     86      !! To solve issues when parent grid is "land" masked but not all the corresponding child  
     87      !! grid points, put Agrif_SpecialValue WHERE the parent grid is masked.  
     88      !! The child solution will be found in the 9(?) points around 
     89      !!----------------------------------------------------------------------- 
     90      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     91      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     92      LOGICAL                         , INTENT(in   ) ::   before 
     93      !! 
     94      REAL(wp) ::   zrhoy   ! local scalar 
    9995      !!----------------------------------------------------------------------- 
    10096      ! 
    10197      IF( before ) THEN  ! parent grid 
    10298         ptab(:,:) = e2u(i1:i2,j1:j2) * u_ice_b(i1:i2,j1:j2) 
    103          WHERE( umask(i1:i2,j1:j2,1) == 0. )  ptab(:,:) = -9999. 
     99         WHERE( umask(i1:i2,j1:j2,1) == 0. )   ptab(i1:i2,j1:j2) = Agrif_SpecialValue 
    104100      ELSE               ! child grid 
    105101         zrhoy = Agrif_Rhoy() 
    106          u_ice(i1:i2,j1:j2) = ptab(:,:) / ( e2u(i1:i2,j1:j2) * zrhoy ) * umask(i1:i2,j1:j2,1) 
     102         u_ice(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) / ( e2u(i1:i2,j1:j2) * zrhoy ) * umask(i1:i2,j1:j2,1) 
    107103      ENDIF 
    108104      ! 
     
    115111      !! 
    116112      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after) 
    117       !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points, 
    118       !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around 
     113      !! To solve issues when parent grid is "land" masked but not all the corresponding child  
     114      !! grid points, put Agrif_SpecialValue WHERE the parent grid is masked.  
     115      !! The child solution will be found in the 9(?) points around 
    119116      !!-----------------------------------------------------------------------       
    120       INTEGER , INTENT(in) :: i1, i2, j1, j2 
    121       REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 
    122       LOGICAL , INTENT(in) :: before 
    123       !! 
    124       REAL(wp) :: zrhox 
     117      INTEGER                         , INTENT(in   ) ::  i1, i2, j1, j2 
     118      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     119      LOGICAL                         , INTENT(in   ) ::  before 
     120      !! 
     121      REAL(wp) ::   zrhox   ! local scalar 
    125122      !!----------------------------------------------------------------------- 
    126123      ! 
    127124      IF( before ) THEN  ! parent grid 
    128125         ptab(:,:) = e1v(i1:i2,j1:j2) * v_ice_b(i1:i2,j1:j2) 
    129          WHERE( vmask(i1:i2,j1:j2,1) == 0. )  ptab(:,:) = -9999. 
     126         WHERE( vmask(i1:i2,j1:j2,1) == 0. )   ptab(i1:i2,j1:j2) = Agrif_SpecialValue 
    130127      ELSE               ! child grid 
    131128         zrhox = Agrif_Rhox() 
    132          v_ice(i1:i2,j1:j2) = ptab(:,:) / ( e1v(i1:i2,j1:j2) * zrhox ) * vmask(i1:i2,j1:j2,1) 
     129         v_ice(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) / ( e1v(i1:i2,j1:j2) * zrhox ) * vmask(i1:i2,j1:j2,1) 
    133130      ENDIF 
    134131      ! 
     
    141138      !! 
    142139      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after) 
    143       !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points, 
    144       !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around 
    145       !!----------------------------------------------------------------------- 
    146       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 
    147       INTEGER , INTENT(in) :: i1, i2, j1, j2, k1, k2 
    148       LOGICAL , INTENT(in) :: before 
    149       INTEGER , INTENT(in) :: nb, ndir 
    150       !! 
    151       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztab 
     140      !! To solve issues when parent grid is "land" masked but not all the corresponding child  
     141      !! grid points, put Agrif_SpecialValue WHERE the parent grid is masked.  
     142      !! The child solution will be found in the 9(?) points around 
     143      !!----------------------------------------------------------------------- 
     144      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab 
     145      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
     146      LOGICAL                               , INTENT(in   ) ::   before 
     147      INTEGER                               , INTENT(in   ) ::   nb, ndir 
     148      !! 
    152149      INTEGER  ::   ji, jj, jk, jl, jm 
    153150      INTEGER  ::   imin, imax, jmin, jmax 
     151      LOGICAL  ::   western_side, eastern_side, northern_side, southern_side 
    154152      REAL(wp) ::   zrhox, z1, z2, z3, z4, z5, z6, z7 
    155       LOGICAL  ::   western_side, eastern_side, northern_side, southern_side 
    156  
    157       !!----------------------------------------------------------------------- 
    158       ! tracers are not multiplied by grid cell here => before: * e12t ; after: * r1_e12t / rhox / rhoy 
     153      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztab 
     154      !!----------------------------------------------------------------------- 
     155      ! tracers are not multiplied by grid cell here => before: * e1e2t ; after: * r1_e1e2t / rhox / rhoy 
    159156      ! and it is ok since we conserve tracers (same as in the ocean). 
    160       ALLOCATE( ztab(SIZE(a_i_b,1),SIZE(a_i_b,2),SIZE(ptab,3)) ) 
     157      ALLOCATE( ztab(SIZE(a_i,1),SIZE(a_i,2),SIZE(ptab,3)) ) 
    161158      
    162159      IF( before ) THEN  ! parent grid 
    163160         jm = 1 
    164161         DO jl = 1, jpl 
    165             ptab(i1:i2,j1:j2,jm) = a_i_b  (i1:i2,j1:j2,jl) ; jm = jm + 1 
    166             ptab(i1:i2,j1:j2,jm) = v_i_b  (i1:i2,j1:j2,jl) ; jm = jm + 1 
    167             ptab(i1:i2,j1:j2,jm) = v_s_b  (i1:i2,j1:j2,jl) ; jm = jm + 1 
    168             ptab(i1:i2,j1:j2,jm) = smv_i_b(i1:i2,j1:j2,jl) ; jm = jm + 1 
    169             ptab(i1:i2,j1:j2,jm) = oa_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1 
     162            ptab(i1:i2,j1:j2,jm  ) = a_i_b (i1:i2,j1:j2,jl) 
     163            ptab(i1:i2,j1:j2,jm+1) = v_i_b (i1:i2,j1:j2,jl) 
     164            ptab(i1:i2,j1:j2,jm+2) = v_s_b (i1:i2,j1:j2,jl) 
     165            ptab(i1:i2,j1:j2,jm+3) = sv_i_b(i1:i2,j1:j2,jl) 
     166            ptab(i1:i2,j1:j2,jm+4) = oa_i_b(i1:i2,j1:j2,jl) 
     167            jm = jm + 5 
    170168            DO jk = 1, nlay_s 
    171                ptab(i1:i2,j1:j2,jm) = e_s_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1 
    172             ENDDO 
     169               ptab(i1:i2,j1:j2,jm) = e_s_b(i1:i2,j1:j2,jk,jl)   ;  jm = jm + 1 
     170            END DO 
    173171            DO jk = 1, nlay_i 
    174                ptab(i1:i2,j1:j2,jm) = e_i_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1 
    175             ENDDO 
    176          ENDDO 
     172               ptab(i1:i2,j1:j2,jm) = e_i_b(i1:i2,j1:j2,jk,jl)   ;  jm = jm + 1 
     173            END DO 
     174         END DO 
    177175          
    178176         DO jk = k1, k2 
    179             WHERE( tmask(i1:i2,j1:j2,1) == 0. )  ptab(i1:i2,j1:j2,jk) = -9999. 
    180          ENDDO 
     177            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   ptab(i1:i2,j1:j2,jk) = Agrif_SpecialValue 
     178         END DO 
     179         ! 
     180      ELSE               ! child grid 
     181         ! 
     182         IF( nbghostcells > 1 ) THEN   ! ==> The easiest interpolation is used 
     183            ! 
     184            jm = 1 
     185            DO jl = 1, jpl 
     186               ! 
     187               DO jj = j1, j2 
     188                  DO ji = i1, i2 
     189                     a_i (ji,jj,jl) = ptab(ji,jj,jm  ) * tmask(ji,jj,1) 
     190                     v_i (ji,jj,jl) = ptab(ji,jj,jm+1) * tmask(ji,jj,1) 
     191                     v_s (ji,jj,jl) = ptab(ji,jj,jm+2) * tmask(ji,jj,1) 
     192                     sv_i(ji,jj,jl) = ptab(ji,jj,jm+3) * tmask(ji,jj,1) 
     193                     oa_i(ji,jj,jl) = ptab(ji,jj,jm+4) * tmask(ji,jj,1) 
     194                  END DO 
     195               END DO 
     196               jm = jm + 5 
     197               ! 
     198               DO jk = 1, nlay_s 
     199                  e_s(i1:i2,j1:j2,jk,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) 
     200                  jm = jm + 1 
     201               END DO 
     202               ! 
     203               DO jk = 1, nlay_i 
     204                  e_i(i1:i2,j1:j2,jk,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) 
     205                  jm = jm + 1 
     206               END DO 
     207               ! 
     208            END DO 
     209            ! 
     210         ELSE                          ! ==> complex interpolation (only one ghost cell available) 
     211            !! Use a more complex interpolation since we mix solutions over a couple of grid points 
     212            !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...) 
     213            !! clem: for some reason (I don't know why), the following lines do not work  
     214            !!       with mpp (or in realistic configurations?). It makes the model crash 
     215            !        I think there is an issue with Agrif_SpecialValue here (not taken into account properly) 
     216            ! record ztab 
     217            jm = 1 
     218            DO jl = 1, jpl 
     219               ztab(:,:,jm  ) = a_i  (:,:,jl) 
     220               ztab(:,:,jm+1) = v_i  (:,:,jl) 
     221               ztab(:,:,jm+2) = v_s  (:,:,jl) 
     222               ztab(:,:,jm+3) = sv_i(:,:,jl) 
     223               ztab(:,:,jm+4) = oa_i(:,:,jl) 
     224               jm = jm + 5 
     225               DO jk = 1, nlay_s 
     226                  ztab(:,:,jm) = e_s(:,:,jk,jl) 
     227                  jm = jm + 1 
     228               END DO 
     229               DO jk = 1, nlay_i 
     230                  ztab(:,:,jm) = e_i(:,:,jk,jl) 
     231                  jm = jm + 1 
     232               END DO 
     233               ! 
     234            END DO 
     235            ! 
     236            ! borders of the domain 
     237            western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2) 
     238            southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2) 
     239            ! 
     240            ! spatial smoothing 
     241            zrhox = Agrif_Rhox() 
     242            z1 =      ( zrhox - 1. ) * 0.5  
     243            z3 =      ( zrhox - 1. ) / ( zrhox + 1. ) 
     244            z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. ) 
     245            z7 =    - ( zrhox - 1. ) / ( zrhox + 3. ) 
     246            z2 = 1. - z1 
     247            z4 = 1. - z3 
     248            z5 = 1. - z6 - z7 
     249            ! 
     250            ! Remove corners 
     251            imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2 
     252            IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3 
     253            IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2 
     254            IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3 
     255            IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2 
     256 
     257            ! smoothed fields 
     258            IF( eastern_side ) THEN 
     259               ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:) 
     260               DO jj = jmin, jmax 
     261                  rswitch = 0. 
     262                  IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1. 
     263                  ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  & 
     264                     &                +      umask(nlci-2,jj,1)   *  & 
     265                     &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  & 
     266                     &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) ) 
     267                  ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1) 
     268               END DO 
     269            ENDIF 
     270            !  
     271            IF( northern_side ) THEN 
     272               ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:) 
     273               DO ji = imin, imax 
     274                  rswitch = 0. 
     275                  IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1. 
     276                  ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  & 
     277                     &                +      vmask(ji,nlcj-2,1)   *  & 
     278                     &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) & 
     279                     &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) ) 
     280                  ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1) 
     281               END DO 
     282            END IF 
     283            ! 
     284            IF( western_side) THEN 
     285               ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:) 
     286               DO jj = jmin, jmax 
     287                  rswitch = 0. 
     288                  IF( u_ice(2,jj) < 0._wp ) rswitch = 1. 
     289                  ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  & 
     290                     &           +      umask(2,jj,1)   *   & 
     291                     &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) & 
     292                     &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) ) 
     293                  ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1) 
     294               END DO 
     295            ENDIF 
     296            ! 
     297            IF( southern_side ) THEN 
     298               ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:) 
     299               DO ji = imin, imax 
     300                  rswitch = 0. 
     301                  IF( v_ice(ji,2) < 0._wp ) rswitch = 1. 
     302                  ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  & 
     303                     &           +      vmask(ji,2,1)   *  & 
     304                     &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) & 
     305                     &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) ) 
     306                  ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1) 
     307               END DO 
     308            END IF 
     309            ! 
     310            ! Treatment of corners 
     311            IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south 
     312            IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north 
     313            IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south 
     314            IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north 
     315             
     316            ! retrieve ice tracers 
     317            jm = 1 
     318            DO jl = 1, jpl 
     319               ! 
     320               DO jj = j1, j2 
     321                  DO ji = i1, i2 
     322                     a_i (ji,jj,jl) = ztab(ji,jj,jm  ) * tmask(ji,jj,1) 
     323                     v_i (ji,jj,jl) = ztab(ji,jj,jm+1) * tmask(ji,jj,1) 
     324                     v_s (ji,jj,jl) = ztab(ji,jj,jm+2) * tmask(ji,jj,1) 
     325                     sv_i(ji,jj,jl) = ztab(ji,jj,jm+3) * tmask(ji,jj,1) 
     326                     oa_i (ji,jj,jl) = ztab(ji,jj,jm+4) * tmask(ji,jj,1) 
     327                  END DO 
     328               END DO 
     329               jm = jm + 5 
     330               ! 
     331               DO jk = 1, nlay_s 
     332                  e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) 
     333                  jm = jm + 1 
     334               END DO 
     335               ! 
     336               DO jk = 1, nlay_i 
     337                  e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) 
     338                  jm = jm + 1 
     339               END DO 
     340               ! 
     341            END DO 
     342           
     343         ENDIF  ! nbghostcells=1 
    181344          
    182       ELSE               ! child grid 
    183 !! ==> The easiest interpolation is the following commented lines 
    184          jm = 1 
    185          DO jl = 1, jpl 
    186             a_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 
    187             v_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 
    188             v_s  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 
    189             smv_i(i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 
    190             oa_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 
    191             DO jk = 1, nlay_s 
    192                e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 
    193             ENDDO 
    194             DO jk = 1, nlay_i 
    195                e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 
    196             ENDDO 
    197          ENDDO 
    198  
    199 !! ==> this is a more complex interpolation since we mix solutions over a couple of grid points 
    200 !!     it is advised to use it for fields modified by high order schemes (e.g. advection UM5...) 
    201 !!        clem: for some reason (I don't know why), the following lines do not work  
    202 !!              with mpp (or in realistic configurations?). It makes the model crash 
    203 !         ! record ztab 
    204 !         jm = 1 
    205 !         DO jl = 1, jpl 
    206 !            ztab(:,:,jm) = a_i  (:,:,jl) ; jm = jm + 1 
    207 !            ztab(:,:,jm) = v_i  (:,:,jl) ; jm = jm + 1 
    208 !            ztab(:,:,jm) = v_s  (:,:,jl) ; jm = jm + 1 
    209 !            ztab(:,:,jm) = smv_i(:,:,jl) ; jm = jm + 1 
    210 !            ztab(:,:,jm) = oa_i (:,:,jl) ; jm = jm + 1 
    211 !            DO jk = 1, nlay_s 
    212 !               ztab(:,:,jm) = e_s(:,:,jk,jl) ; jm = jm + 1 
    213 !            ENDDO 
    214 !            DO jk = 1, nlay_i 
    215 !               ztab(:,:,jm) = e_i(:,:,jk,jl) ; jm = jm + 1 
    216 !            ENDDO 
    217 !         ENDDO 
    218 !         ! 
    219 !         ! borders of the domain 
    220 !         western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2) 
    221 !         southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2) 
    222 !         ! 
    223 !         ! spatial smoothing 
    224 !         zrhox = Agrif_Rhox() 
    225 !         z1 =      ( zrhox - 1. ) * 0.5  
    226 !         z3 =      ( zrhox - 1. ) / ( zrhox + 1. ) 
    227 !         z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. ) 
    228 !         z7 =    - ( zrhox - 1. ) / ( zrhox + 3. ) 
    229 !         z2 = 1. - z1 
    230 !         z4 = 1. - z3 
    231 !         z5 = 1. - z6 - z7 
    232 !         ! 
    233 !         ! Remove corners 
    234 !         imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2 
    235 !         IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3 
    236 !         IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2 
    237 !         IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3 
    238 !         IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2 
    239 ! 
    240 !         ! smoothed fields 
    241 !         IF( eastern_side ) THEN 
    242 !            ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:) 
    243 !            DO jj = jmin, jmax 
    244 !               rswitch = 0. 
    245 !               IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1. 
    246 !               ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  & 
    247 !                  &                +      umask(nlci-2,jj,1)   *  & 
    248 !                  &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  & 
    249 !                  &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) ) 
    250 !               ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1) 
    251 !            END DO 
    252 !         ENDIF 
    253 !         !  
    254 !         IF( northern_side ) THEN 
    255 !            ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:) 
    256 !            DO ji = imin, imax 
    257 !               rswitch = 0. 
    258 !               IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1. 
    259 !               ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  & 
    260 !                  &                +      vmask(ji,nlcj-2,1)   *  & 
    261 !                  &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) & 
    262 !                  &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) ) 
    263 !               ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1) 
    264 !            END DO 
    265 !         END IF 
    266 !         ! 
    267 !         IF( western_side) THEN 
    268 !            ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:) 
    269 !            DO jj = jmin, jmax 
    270 !               rswitch = 0. 
    271 !               IF( u_ice(2,jj) < 0._wp ) rswitch = 1. 
    272 !               ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  & 
    273 !                  &           +      umask(2,jj,1)   *   & 
    274 !                  &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) & 
    275 !                  &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) ) 
    276 !               ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1) 
    277 !            END DO 
    278 !         ENDIF 
    279 !         ! 
    280 !         IF( southern_side ) THEN 
    281 !            ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:) 
    282 !            DO ji = imin, imax 
    283 !               rswitch = 0. 
    284 !               IF( v_ice(ji,2) < 0._wp ) rswitch = 1. 
    285 !               ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  & 
    286 !                  &           +      vmask(ji,2,1)   *  & 
    287 !                  &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) & 
    288 !                  &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) ) 
    289 !               ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1) 
    290 !            END DO 
    291 !         END IF 
    292 !         ! 
    293 !         ! Treatment of corners 
    294 !         IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south 
    295 !         IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north 
    296 !         IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south 
    297 !         IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north 
    298 ! 
    299 !         ! retrieve ice tracers 
    300 !         jm = 1 
    301 !         DO jl = 1, jpl 
    302 !            a_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1 
    303 !            v_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1 
    304 !            v_s  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1 
    305 !            smv_i(i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1 
    306 !            oa_i (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1 
    307 !            DO jk = 1, nlay_s 
    308 !               e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1 
    309 !            ENDDO 
    310 !            DO jk = 1, nlay_i 
    311 !               e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1 
    312 !            ENDDO 
    313 !         ENDDO 
    314         
    315345         ! integrated values 
    316346         vt_i (i1:i2,j1:j2) = SUM( v_i(i1:i2,j1:j2,:), dim=3 ) 
     
    319349         et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 ) 
    320350         et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 ) 
    321  
     351         ! 
    322352      ENDIF 
    323353       
     
    327357 
    328358#else 
     359   !!---------------------------------------------------------------------- 
     360   !!   Empty module                                             no sea-ice 
     361   !!---------------------------------------------------------------------- 
    329362CONTAINS 
    330363   SUBROUTINE agrif_lim3_interp_empty 
    331       !!--------------------------------------------- 
    332       !!   *** ROUTINE agrif_lim3_interp_empty *** 
    333       !!--------------------------------------------- 
    334364      WRITE(*,*)  'agrif_lim3_interp : You should not have seen this print! error?' 
    335365   END SUBROUTINE agrif_lim3_interp_empty 
    336366#endif 
     367 
     368   !!====================================================================== 
    337369END MODULE agrif_lim3_interp 
Note: See TracChangeset for help on using the changeset viewer.