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 1530 – NEMO

Changeset 1530


Ignore:
Timestamp:
2009-07-23T18:22:51+02:00 (15 years ago)
Author:
ctlod
Message:

only style changes in advection modules of LIM2.0 and LIM3.0

Location:
trunk/NEMO
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMO/LIM_SRC_2/limadv_2.F90

    r1529 r1530  
    44   !! LIM 2.0 sea-ice model : sea-ice advection 
    55   !!====================================================================== 
     6   !! History :  OPA  ! 2000-01 (LIM)  Original code 
     7   !!                 ! 2001-05 (G. Madec, R. Hordoir) Doctor norm 
     8   !!   NEMO     1.0  ! 2003-10 (C. Ethe) F90, module 
     9   !!             -   ! 2003-12 (R. Hordoir, G. Madec) mpp 
     10   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c. 
     11   !!-------------------------------------------------------------------- 
    612#if defined key_lim2 
    713   !!---------------------------------------------------------------------- 
    … …  
    1117   !!   lim_adv_y_2  : advection of sea ice on y axis 
    1218   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    1419   USE dom_oce 
    1520   USE dom_ice_2 
    16    USE in_out_manager  ! I/O manager 
    1721   USE ice_2 
    1822   USE lbclnk 
    19    USE prtctl          ! Print control 
     23   USE in_out_manager     ! I/O manager 
     24   USE prtctl             ! Print control 
    2025 
    2126   IMPLICIT NONE 
    2227   PRIVATE 
    2328 
    24    !! * Routine accessibility 
    25    PUBLIC lim_adv_x_2    ! called by lim_trp 
    26    PUBLIC lim_adv_y_2    ! called by lim_trp 
    27  
     29   PUBLIC   lim_adv_x_2   ! called by lim_trp 
     30   PUBLIC   lim_adv_y_2   ! called by lim_trp 
     31 
     32   REAL(wp)  ::   epsi20 = 1.e-20   ! constant values 
     33   REAL(wp)  ::   rzero  = 0.e0     !    -       - 
     34   REAL(wp)  ::   rone   = 1.e0     !    -       - 
     35    
    2836   !! * Substitutions 
    2937#  include "vectopt_loop_substitute.h90" 
    30  
    31    !! * Module variables 
    32    REAL(wp)  ::            &  ! constant values 
    33       epsi20 = 1e-20    ,  & 
    34       rzero  = 0.e0     ,  & 
    35       rone   = 1.e0 
    36    !!---------------------------------------------------------------------- 
    37    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
     38   !!---------------------------------------------------------------------- 
     39   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2009)  
    3840   !! $Id$  
    39    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     41   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4042   !!---------------------------------------------------------------------- 
    4143 
    … …  
    4850      !!   
    4951      !! ** purpose :   Computes and adds the advection trend to sea-ice 
    50       !!      variable on x axis 
     52      !!              variable on i-axis 
    5153      !! 
    52       !! ** method  :   Uses Prather second order scheme that advects 
    53       !!      tracers but also theirquadratic forms. The method preserves 
    54       !!      tracer structures by conserving second order moments. 
    55       !!      Ref.: "Numerical Advection by Conservation of Second Order 
    56       !!      Moments", JGR, VOL. 91. NO. D6. PAGES 6671-6681. MAY 20, 1986 
    57       !!      
    58       !! History : 
    59       !!        !  00-01 (LIM) 
    60       !!        !  01-05 (G. Madec, R. Hordoir) opa norm 
    61       !!        !  03-10 (C. Ethe) F90, module 
    62       !!        !  03-12 (R. Hordoir, G. Madec) mpp 
     54      !! ** method  :   Uses Prather second order scheme that advects tracers 
     55      !!              but also theirquadratic forms. The method preserves 
     56      !!              tracer structures by conserving second order moments. 
     57      !! 
     58      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681. 
    6359      !!-------------------------------------------------------------------- 
    64       !! * Arguments 
    65       REAL(wp)                    , INTENT(in)  ::  & 
    66          pdf ,       &  ! ??? 
    67          pcrh           ! = 1. : lim_adv_x is called before lim_adv_y 
    68          !              ! = 0. : lim_adv_x is called after  lim_adv_y 
    69       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  ::  & 
    70          put            ! i-direction ice velocity at ocean U-point (m/s) 
    71       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::  &  
    72          ps0 , psm , &  ! ??? 
    73          psx , psy , &  ! ???  
    74          psxx, psyy, psxy 
    75  
    76       !! * Local declarations 
    77       INTEGER ::   ji, jj      ! dummy loop indices 
    78       REAL(wp)       ::  & 
    79          zrdt, zslpmax, ztemp, zin0,     &  ! temporary scalars 
    80          zs1max, zs1new, zs2new,         &  !    "         " 
    81          zalf, zalfq, zalf1, zalf1q,     &  !    "         " 
    82          zbt , zbt1                         !    "         " 
    83       REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace 
    84          zf0 , zfx , zfy , zbet,         &  !    "           " 
    85          zfxx, zfyy, zfxy,               &  !    "           " 
    86          zfm, zalg, zalg1, zalg1q           !    "           " 
     60      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
     61      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0) 
     62      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s] 
     63      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area 
     64      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected 
     65      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments  
     66      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
     67      !!  
     68      INTEGER  ::   ji, jj                               ! dummy loop indices 
     69      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars 
     70      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
     71      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     72      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
     73      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
     74      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    8775      !--------------------------------------------------------------------- 
    8876 
    8977      ! Limitation of moments.                                            
    9078 
    91       zrdt      = rdt_ice * pdf   ! If ice drift field is too fast, use an appropriate time step for advection. 
     79      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection. 
    9280 
    9381      DO jj = 1, jpj 
    … …  
    9987               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  ) 
    10088            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
    101  
     89            ! 
    10290            ps0 (ji,jj) = zslpmax   
    10391            psx (ji,jj) = zs1new      * zin0 
    … …  
    120108            zalf1        =  1.0 - zalf 
    121109            zalf1q       =  zalf1 * zalf1 
     110            ! 
    122111            zfm (ji,jj)  =  zalf  * psm(ji,jj) 
    123112            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj)  ) ) 
    … …  
    127116            zfxy(ji,jj)  =  zalfq * psxy(ji,jj) 
    128117            zfyy(ji,jj)  =  zalf  * psyy(ji,jj) 
    129  
     118            ! 
    130119            !  Readjust moments remaining in the box. 
    131120            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj) 
    … …  
    181170            zalf1       = 1.0 - zalf 
    182171            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj) 
    183             ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji-1,jj)) + zbt1 * ps0(ji,jj) 
    184             psx(ji,jj)  = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj) 
    185             psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)   & 
     172            ! 
     173            ps0 (ji,jj) = zbt * (ps0(ji,jj) + zf0(ji-1,jj)) + zbt1 * ps0(ji,jj) 
     174            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj) 
     175            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               & 
    186176               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   & 
    187177               &        + zbt1 * psxx(ji,jj) 
    … …  
    202192            zalf1       = 1.0 - zalf 
    203193            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 
    204             ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
    205             psx(ji,jj)  = zbt  * psx(ji,jj)   & 
    206                &        + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) 
    207             psxx(ji,jj) = zbt  * psxx(ji,jj)   & 
    208                &        + zbt1 * (  zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  & 
    209                &                 + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) 
    210             psxy(ji,jj) = zbt  * psxy(ji,jj)   & 
    211                &        + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  & 
    212                &                 + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  ) 
    213             psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) ) 
    214             psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) ) 
     194            ! 
     195            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
     196            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) 
     197            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * (  zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)   & 
     198               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )        & 
     199               &                                      + ( zalf1 - zalf ) * ztemp )                                 ) 
     200            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)          & 
     201               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  ) 
     202            psy(ji,jj)  = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) + zfy (ji,jj) ) 
     203            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) ) 
    215204         END DO 
    216205      END DO 
    217206 
    218207      !-- Lateral boundary conditions 
    219       CALL lbc_lnk( psm , 'T',  1. ) 
    220       CALL lbc_lnk( ps0 , 'T',  1. ) 
    221       CALL lbc_lnk( psx , 'T', -1. ) 
    222       CALL lbc_lnk( psxx, 'T',  1. ) 
    223       CALL lbc_lnk( psy , 'T', -1. ) 
    224       CALL lbc_lnk( psyy, 'T',  1. ) 
     208      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. ) 
     209      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes 
     210      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. ) 
    225211      CALL lbc_lnk( psxy, 'T',  1. ) 
    226212 
    … …  
    231217         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :') 
    232218      ENDIF 
    233  
     219      ! 
    234220   END SUBROUTINE lim_adv_x_2 
    235221 
    … …  
    241227      !!             
    242228      !! ** purpose :   Computes and adds the advection trend to sea-ice  
    243       !!      variable on y axis 
     229      !!              variable on j-axis 
    244230      !! 
    245231      !! ** method  :   Uses Prather second order scheme that advects tracers 
    246       !!      but also their quadratic forms. The method preserves tracer 
    247       !!      structures by conserving second order moments. 
     232      !!              but also their quadratic forms. The method preserves  
     233      !!              tracer structures by conserving second order moments. 
    248234      !! 
    249       !! History : 
    250       !!   1.0  !  00-01 (LIM) 
    251       !!        !  01-05 (G. Madec, R. Hordoir) opa norm 
    252       !!   2.0  !  03-10 (C. Ethe) F90, module 
    253       !!        !  03-12 (R. Hordoir, G. Madec) mpp 
     235      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681. 
    254236      !!--------------------------------------------------------------------- 
    255       !! * Arguments 
    256       REAL(wp),                     INTENT(in)  :: & 
    257          pdf,        &  ! ??? 
    258          pcrh           ! = 1. : lim_adv_x is called before lim_adv_y 
    259          !              ! = 0. : lim_adv_x is called after  lim_adv_y 
    260       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: & 
    261          pvt            ! j-direction ice velocity at ocean V-point (m/s) 
    262       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: & 
    263          psm , ps0 , psx , psy,   & 
    264          psxx, psyy, psxy 
    265  
    266       !! * Local Variables 
    267       INTEGER  ::   ji, jj                ! dummy loop indices 
    268       REAL(wp) ::   & 
    269          zrdt, zslpmax, zin0, ztemp,  &   ! temporary scalars 
    270          zs1max, zs1new, zs2new,      &   !    "         " 
    271          zalf, zalfq, zalf1, zalf1q,  &   !    "         " 
    272          zbt , zbt1                       ! 
    273       REAL(wp), DIMENSION(jpi,jpj)  :: & 
    274          zf0 , zfx , zfy ,             &  ! temporary workspace 
    275          zfxx, zfyy, zfxy,             &  !    "           " 
    276          zfm , zbet,                   &  !    "           " 
    277          zalg, zalg1, zalg1q              !    "           " 
     237      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
     238      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0) 
     239      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s] 
     240      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area 
     241      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected 
     242      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments  
     243      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
     244      !! 
     245      INTEGER  ::   ji, jj                               ! dummy loop indices 
     246      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars 
     247      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
     248      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     249      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
     250      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
     251      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    278252      !--------------------------------------------------------------------- 
    279253 
    … …  
    282256      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 
    283257 
    284        DO jj = 1, jpj 
    285           DO ji = 1, jpi 
    286              zslpmax = MAX( rzero, ps0(ji,jj) ) 
    287              zs1max  = 1.5 * zslpmax 
    288              zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) ) 
    289              zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    290                 &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  ) 
    291              zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
    292              ps0 (ji,jj) = zslpmax   
    293              psx (ji,jj)  = psx (ji,jj) * zin0 
    294              psxx(ji,jj)  = psxx(ji,jj) * zin0 
    295              psy (ji,jj) = zs1new * zin0 
    296              psyy(ji,jj) = zs2new * zin0 
    297              psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0 
    298           END DO 
    299        END DO 
    300  
    301        !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    302        psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 
    303  
    304        !  Calculate fluxes and moments between boxes j<-->j+1               
    305        DO jj = 1, jpj 
    306           DO ji = 1, jpi 
    307              !  Flux from j to j+1 WHEN v GT 0    
    308              zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) ) 
    309              zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj) 
    310              zalfq        =  zalf * zalf 
    311              zalf1        =  1.0 - zalf 
    312              zalf1q       =  zalf1 * zalf1 
    313              zfm (ji,jj)  =  zalf  * psm(ji,jj) 
    314              zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) )  
    315              zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) ) 
    316              zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj) 
    317              zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) ) 
    318              zfxy(ji,jj)  =  zalfq * psxy(ji,jj) 
    319              zfxx(ji,jj)  =  zalf  * psxx(ji,jj) 
    320  
    321              !  Readjust moments remaining in the box. 
    322              psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj) 
    323              ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj) 
    324              psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) ) 
    325              psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj) 
    326              psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj) 
    327              psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj) 
    328              psxy(ji,jj)  =  zalf1q * psxy(ji,jj) 
    329           END DO 
    330        END DO 
    331  
    332        DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
    333           DO ji = 1, jpi 
    334              zalf          = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1)  
    335              zalg  (ji,jj) = zalf 
    336              zalfq         = zalf * zalf 
    337              zalf1         = 1.0 - zalf 
    338              zalg1 (ji,jj) = zalf1 
    339              zalf1q        = zalf1 * zalf1 
    340              zalg1q(ji,jj) = zalf1q 
    341              zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji,jj+1) 
    342              zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 
    343              zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 
    344              zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * zalfq * psyy(ji,jj+1) 
    345              zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 
    346              zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1) 
    347              zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * psxx(ji,jj+1) 
    348           END DO 
    349        END DO 
     258      DO jj = 1, jpj 
     259         DO ji = 1, jpi 
     260            zslpmax = MAX( rzero, ps0(ji,jj) ) 
     261            zs1max  = 1.5 * zslpmax 
     262            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) ) 
     263            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
     264               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  ) 
     265            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
     266            ! 
     267            ps0 (ji,jj) = zslpmax   
     268            psx (ji,jj) = psx (ji,jj) * zin0 
     269            psxx(ji,jj) = psxx(ji,jj) * zin0 
     270            psy (ji,jj) = zs1new * zin0 
     271            psyy(ji,jj) = zs2new * zin0 
     272            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0 
     273         END DO 
     274      END DO 
     275 
     276      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
     277      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  ) 
     278 
     279      !  Calculate fluxes and moments between boxes j<-->j+1               
     280      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0    
     281         DO ji = 1, jpi 
     282            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) ) 
     283            zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj) 
     284            zalfq        =  zalf * zalf 
     285            zalf1        =  1.0 - zalf 
     286            zalf1q       =  zalf1 * zalf1 
     287            zfm (ji,jj)  =  zalf  * psm(ji,jj) 
     288            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) )  
     289            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) ) 
     290            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj) 
     291            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) ) 
     292            zfxy(ji,jj)  =  zalfq * psxy(ji,jj) 
     293            zfxx(ji,jj)  =  zalf  * psxx(ji,jj) 
     294            ! 
     295            !  Readjust moments remaining in the box. 
     296            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj) 
     297            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj) 
     298            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) ) 
     299            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj) 
     300            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj) 
     301            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj) 
     302            psxy(ji,jj)  =  zalf1q * psxy(ji,jj) 
     303         END DO 
     304      END DO 
     305      ! 
     306      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
     307         DO ji = 1, jpi 
     308            zalf          = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1)  
     309            zalg  (ji,jj) = zalf 
     310            zalfq         = zalf * zalf 
     311            zalf1         = 1.0 - zalf 
     312            zalg1 (ji,jj) = zalf1 
     313            zalf1q        = zalf1 * zalf1 
     314            zalg1q(ji,jj) = zalf1q 
     315            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji,jj+1) 
     316            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 
     317            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 
     318            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * zalfq * psyy(ji,jj+1) 
     319            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 
     320            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1) 
     321            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * psxx(ji,jj+1) 
     322         END DO 
     323      END DO 
    350324  
    351        !  Readjust moments remaining in the box.  
    352        DO jj = 2, jpj 
    353           DO ji = 1, jpi 
    354              zbt  =         zbet(ji,jj-1) 
    355              zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    356              psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 
    357              ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 
    358              psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) ) 
    359              psyy(ji,jj) = zalg1 (ji,jj-1)  * zalg1q(ji,jj-1) * psyy(ji,jj) 
    360              psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 
    361              psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) 
    362              psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj) 
    363           END DO 
    364        END DO 
    365  
    366        !   Put the temporary moments into appropriate neighboring boxes.     
    367        DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0. 
    368           DO ji = 1, jpi 
    369              zbt  =         zbet(ji,jj-1) 
    370              zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    371              psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj)  
    372              zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj)  
    373              zalf1       = 1.0 - zalf 
    374              ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1) 
    375              ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj) 
    376  
    377              psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   & 
    378                 &        + zbt1 * psy(ji,jj)   
    379  
    380              psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             & 
    381                 &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   &  
    382                 &        + zbt1 * psyy(ji,jj) 
    383  
    384              psxy(ji,jj) = zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)              & 
    385                                   + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) )   & 
    386                          + zbt1 * psxy(ji,jj) 
    387              psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 
    388              psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) 
    389           END DO 
    390        END DO 
    391  
    392        DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0. 
    393           DO ji = 1, jpi 
    394              zbt  =         zbet(ji,jj) 
    395              zbt1 = ( 1.0 - zbet(ji,jj) ) 
    396              psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 
    397              zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj) 
    398              zalf1       = 1.0 - zalf 
    399              ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 
    400              ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
    401              psy(ji,jj)  = zbt  * psy(ji,jj)  & 
    402                 &        + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 
    403              psyy(ji,jj) = zbt  * psyy(ji,jj)  & 
    404                 &        + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   & 
    405                 &                 + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) 
    406              psxy(ji,jj) = zbt  * psxy(ji,jj)   & 
    407                 &        + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)   & 
    408                 &                 + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) ) 
    409              psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 
    410              psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 
    411           END DO 
    412        END DO 
     325      !  Readjust moments remaining in the box.  
     326      DO jj = 2, jpj 
     327         DO ji = 1, jpi 
     328            zbt  =         zbet(ji,jj-1) 
     329            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
     330            ! 
     331            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 
     332            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 
     333            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) ) 
     334            psyy(ji,jj) = zalg1 (ji,jj-1)  * zalg1q(ji,jj-1) * psyy(ji,jj) 
     335            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 
     336            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) 
     337            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj) 
     338         END DO 
     339      END DO 
     340 
     341      !   Put the temporary moments into appropriate neighboring boxes.     
     342      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0. 
     343         DO ji = 1, jpi 
     344            zbt  =         zbet(ji,jj-1) 
     345            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
     346            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj)  
     347            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj)  
     348            zalf1       = 1.0 - zalf 
     349            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1) 
     350            ! 
     351            ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj) 
     352            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   & 
     353               &        + zbt1 * psy(ji,jj)   
     354            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             & 
     355               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   &  
     356               &        + zbt1 * psyy(ji,jj) 
     357            psxy(ji,jj) = zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)              & 
     358               &                 + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) )   & 
     359               &        + zbt1 * psxy(ji,jj) 
     360            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 
     361            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) 
     362         END DO 
     363      END DO 
     364      ! 
     365      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0. 
     366         DO ji = 1, jpi 
     367            zbt  =         zbet(ji,jj) 
     368            zbt1 = ( 1.0 - zbet(ji,jj) ) 
     369            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 
     370            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj) 
     371            zalf1       = 1.0 - zalf 
     372            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 
     373            ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
     374            psy(ji,jj)  = zbt  * psy(ji,jj)  & 
     375               &        + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 
     376            psyy(ji,jj) = zbt  * psyy(ji,jj)                                                                         & 
     377               &        + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)                           & 
     378               &                 + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) 
     379            psxy(ji,jj) = zbt  * psxy(ji,jj)                                  & 
     380               &        + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)   & 
     381               &                 + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) ) 
     382            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 
     383            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 
     384         END DO 
     385      END DO 
    413386 
    414387      !-- Lateral boundary conditions 
    415       CALL lbc_lnk( psm , 'T',  1. ) 
    416       CALL lbc_lnk( ps0 , 'T',  1. ) 
    417       CALL lbc_lnk( psx , 'T', -1. ) 
    418       CALL lbc_lnk( psxx, 'T',  1. ) 
    419       CALL lbc_lnk( psy , 'T', -1. ) 
    420       CALL lbc_lnk( psyy, 'T',  1. ) 
     388      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. ) 
     389      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes 
     390      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. ) 
    421391      CALL lbc_lnk( psxy, 'T',  1. ) 
    422392 
    … …  
    427397         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :') 
    428398      ENDIF 
    429  
     399      ! 
    430400   END SUBROUTINE lim_adv_y_2 
    431401 
    … …  
    434404   !!   Default option            Dummy module     NO LIM 2.0 sea-ice model 
    435405   !!---------------------------------------------------------------------- 
    436 CONTAINS 
    437    SUBROUTINE lim_adv_x_2         ! Empty routine 
    438    END SUBROUTINE lim_adv_x_2 
    439    SUBROUTINE lim_adv_y_2         ! Empty routine 
    440    END SUBROUTINE lim_adv_y_2 
    441  
    442406#endif 
    443407 
     408   !!====================================================================== 
    444409END MODULE limadv_2 

  • trunk/NEMO/LIM_SRC_3/limadv.F90

    r1529 r1530  
    44   !! LIM sea-ice model : sea-ice advection 
    55   !!====================================================================== 
     6   !! History :  LIM  ! 2008-03 (M. Vancoppenolle)  LIM-3 from LIM-2 code 
     7   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c. 
     8   !!-------------------------------------------------------------------- 
    69#if defined key_lim3 
    710   !!---------------------------------------------------------------------- 
    … …  
    1114   !!   lim_adv_y  : advection of sea ice on y axis 
    1215   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    1416   USE dom_oce 
    1517   USE dom_ice 
    1618   USE ice 
     19   USE lbclnk 
    1720   USE in_out_manager  ! I/O manager 
    18    USE lbclnk 
    1921   USE prtctl          ! Print control 
    2022 
    … …  
    2224   PRIVATE 
    2325 
    24    !! * Routine accessibility 
    25    PUBLIC lim_adv_x    ! called by lim_trp 
    26    PUBLIC lim_adv_y    ! called by lim_trp 
     26   PUBLIC   lim_adv_x   ! called by lim_trp 
     27   PUBLIC   lim_adv_y   ! called by lim_trp 
     28 
     29   REAL(wp)  ::   epsi20 = 1.e-20   ! constant values 
     30   REAL(wp)  ::   rzero  = 0.e0     !    -       - 
     31   REAL(wp)  ::   rone   = 1.e0     !    -       - 
    2732 
    2833   !! * Substitutions 
    2934#  include "vectopt_loop_substitute.h90" 
    30  
    31    !! * Module variables 
    32    REAL(wp)  ::            &  ! constant values 
    33       epsi20 = 1e-20    ,  & 
    34       rzero  = 0.e0     ,  & 
    35       rone   = 1.e0 
    36    !!---------------------------------------------------------------------- 
    37    !!   LIM 3.0,  UCL-LOCEAN-IPSL (2005)  
     35   !!---------------------------------------------------------------------- 
     36   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2009)  
    3837   !! $Id$ 
    39    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     38   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4039   !!---------------------------------------------------------------------- 
    4140 
    … …  
    4847      !!   
    4948      !! ** purpose :   Computes and adds the advection trend to sea-ice 
    50       !!      variable on x axis 
     49      !!              variable on i-axis 
    5150      !! 
    52       !! ** method  :   Uses Prather second order scheme that advects 
    53       !!      tracers but also theirquadratic forms. The method preserves 
    54       !!      tracer structures by conserving second order moments. 
    55       !!      Ref.: "Numerical Advection by Conservation of Second Order 
    56       !!      Moments", JGR, VOL. 91. NO. D6. PAGES 6671-6681. MAY 20, 1986 
    57       !!      
    58       !! History : 
    59       !!        !  00-01 (LIM) 
    60       !!        !  01-05 (G. Madec, R. Hordoir) opa norm 
    61       !!        !  03-10 (C. Ethe) F90, module 
    62       !!        !  03-12 (R. Hordoir, G. Madec) mpp 
     51      !! ** method  :   Uses Prather second order scheme that advects tracers 
     52      !!              but also theirquadratic forms. The method preserves 
     53      !!              tracer structures by conserving second order moments. 
     54      !! 
     55      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681. 
    6356      !!-------------------------------------------------------------------- 
    64       !! * Arguments 
    65       REAL(wp)                    , INTENT(in)  ::  & 
    66          pdf ,       &  ! ??? 
    67          pcrh           ! = 1. : lim_adv_x is called before lim_adv_y 
    68       !              ! = 0. : lim_adv_x is called after  lim_adv_y 
    69       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  ::  & 
    70          put            ! i-direction ice velocity at ocean U-point (m/s) 
    71       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::  &  
    72          ps0 , psm , &  ! ??? 
    73          psx , psy , &  ! ???  
    74          psxx, psyy, psxy 
    75  
    76       !! * Local declarations 
    77       INTEGER ::   ji, jj      ! dummy loop indices 
    78       REAL(wp)       ::  & 
    79          zrdt, zslpmax, ztemp, zin0,     &  ! temporary scalars 
    80          zs1max, zs1new, zs2new,         &  !    "         " 
    81          zalf, zalfq, zalf1, zalf1q,     &  !    "         " 
    82          zbt , zbt1                         !    "         " 
    83       REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace 
    84          zf0 , zfx , zfy , zbet,         &  !    "           " 
    85          zfxx, zfyy, zfxy,               &  !    "           " 
    86          zfm, zalg, zalg1, zalg1q           !    "           " 
     57      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
     58      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0) 
     59      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s] 
     60      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area 
     61      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected 
     62      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments  
     63      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
     64      !!  
     65      INTEGER  ::   ji, jj                               ! dummy loop indices 
     66      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars 
     67      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
     68      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     69      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
     70      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
     71      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    8772      !--------------------------------------------------------------------- 
    8873 
    8974      ! Limitation of moments.                                            
    9075 
    91       zrdt      = rdt_ice * pdf   ! If ice drift field is too fast, use an appropriate time step for advection. 
     76      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection. 
    9277 
    9378      DO jj = 1, jpj 
    … …  
    120105            zalf1        =  1.0 - zalf 
    121106            zalf1q       =  zalf1 * zalf1 
    122             zfm (ji,jj)  =  zalf  * psm(ji,jj) 
    123             zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj)  ) ) 
    124             zfx (ji,jj)  =  zalfq * ( psx(ji,jj) + 3.0 * zalf1 * psxx(ji,jj) ) 
    125             zfxx(ji,jj)  =  zalf  * zalfq * psxx(ji,jj) 
    126             zfy (ji,jj)  =  zalf  * ( psy(ji,jj) + zalf1 * psxy(ji,jj) ) 
    127             zfxy(ji,jj)  =  zalfq * psxy(ji,jj) 
    128             zfyy(ji,jj)  =  zalf  * psyy(ji,jj) 
     107            ! 
     108            zfm (ji,jj)  =  zalf  *   psm (ji,jj) 
     109            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) )  ) 
     110            zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) ) 
     111            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq 
     112            zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) ) 
     113            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj) 
     114            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj) 
    129115 
    130116            !  Readjust moments remaining in the box. 
    … …  
    148134            zalf1q        = zalf1 * zalf1 
    149135            zalg1q(ji,jj) = zalf1q 
    150             zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji+1,jj) 
    151             zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) ) 
    152             zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx(ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) ) 
    153             zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * zalfq * psxx(ji+1,jj) 
    154             zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy(ji+1,jj) - zalf1 * psxy(ji+1,jj) ) 
    155             zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj) 
    156             zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * psyy(ji+1,jj) 
     136            ! 
     137            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji+1,jj) 
     138            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) ) 
     139            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) ) 
     140            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji+1,jj) * zalfq 
     141            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) ) 
     142            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji+1,jj) 
     143            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji+1,jj) 
    157144         END DO 
    158145      END DO 
    … …  
    162149            zbt  =       zbet(ji-1,jj) 
    163150            zbt1 = 1.0 - zbet(ji-1,jj) 
     151            ! 
    164152            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) ) 
    165153            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) ) 
    … …  
    181169            zalf1       = 1.0 - zalf 
    182170            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj) 
    183             ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji-1,jj)) + zbt1 * ps0(ji,jj) 
    184             psx(ji,jj)  = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj) 
    185             psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)   & 
     171            ! 
     172            ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj) 
     173            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj) 
     174            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               & 
    186175               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   & 
    187                &        + zbt1 * psxx(ji,jj) 
     176               &                                                + zbt1 * psxx(ji,jj) 
    188177            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             & 
    189178               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   & 
    190                &         + zbt1 * psxy(ji,jj) 
     179               &                                                + zbt1 * psxy(ji,jj) 
    191180            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj) 
    192181            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj) 
    … …  
    201190            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj) 
    202191            zalf1       = 1.0 - zalf 
    203             ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 
    204             ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
    205             psx(ji,jj)  = zbt  * psx(ji,jj)   & 
    206                &        + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) 
    207             psxx(ji,jj) = zbt  * psxx(ji,jj)   & 
    208                &        + zbt1 * (  zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  & 
    209                &                 + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) 
    210             psxy(ji,jj) = zbt  * psxy(ji,jj)   & 
    211                &        + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  & 
    212                &                 + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  ) 
     192            ztemp       = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 
     193            ! 
     194            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
     195            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) 
     196            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  & 
     197               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )      & 
     198               &                                      + ( zalf1 - zalf ) * ztemp ) ) 
     199            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  & 
     200               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  ) 
    213201            psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) ) 
    214202            psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) ) 
    … …  
    217205 
    218206      !-- Lateral boundary conditions 
    219       CALL lbc_lnk( psm , 'T',  1. ) 
    220       CALL lbc_lnk( ps0 , 'T',  1. ) 
    221       CALL lbc_lnk( psx , 'T', -1. ) 
    222       CALL lbc_lnk( psxx, 'T',  1. ) 
    223       CALL lbc_lnk( psy , 'T', -1. ) 
    224       CALL lbc_lnk( psyy, 'T',  1. ) 
     207      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. ) 
     208      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes 
     209      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. ) 
    225210      CALL lbc_lnk( psxy, 'T',  1. ) 
    226211 
    … …  
    231216         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :') 
    232217      ENDIF 
    233  
     218      ! 
    234219   END SUBROUTINE lim_adv_x 
    235220 
    … …  
    241226      !!             
    242227      !! ** purpose :   Computes and adds the advection trend to sea-ice  
    243       !!      variable on y axis 
     228      !!              variable on y axis 
    244229      !! 
    245230      !! ** method  :   Uses Prather second order scheme that advects tracers 
    246       !!      but also their quadratic forms. The method preserves tracer 
    247       !!      structures by conserving second order moments. 
     231      !!              but also their quadratic forms. The method preserves  
     232      !!              tracer structures by conserving second order moments. 
     233      !!  
     234      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681. 
     235      !!--------------------------------------------------------------------- 
     236      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
     237      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0) 
     238      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s] 
     239      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area 
     240      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected 
     241      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments  
     242      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
    248243      !! 
    249       !! History : 
    250       !!   1.0  !  00-01 (LIM) 
    251       !!        !  01-05 (G. Madec, R. Hordoir) opa norm 
    252       !!   2.0  !  03-10 (C. Ethe) F90, module 
    253       !!        !  03-12 (R. Hordoir, G. Madec) mpp 
    254       !!--------------------------------------------------------------------- 
    255       !! * Arguments 
    256       REAL(wp),                     INTENT(in)  :: & 
    257          pdf,        &  ! ??? 
    258          pcrh           ! = 1. : lim_adv_x is called before lim_adv_y 
    259       !              ! = 0. : lim_adv_x is called after  lim_adv_y 
    260       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: & 
    261          pvt            ! j-direction ice velocity at ocean V-point (m/s) 
    262       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: & 
    263          psm , ps0 , psx , psy,   & 
    264          psxx, psyy, psxy 
    265  
    266       !! * Local Variables 
    267       INTEGER  ::   ji, jj                ! dummy loop indices 
    268       REAL(wp) ::   & 
    269          zrdt, zslpmax, zin0, ztemp,  &   ! temporary scalars 
    270          zs1max, zs1new, zs2new,      &   !    "         " 
    271          zalf, zalfq, zalf1, zalf1q,  &   !    "         " 
    272          zbt , zbt1                       ! 
    273       REAL(wp), DIMENSION(jpi,jpj)  :: & 
    274          zf0 , zfx , zfy ,             &  ! temporary workspace 
    275          zfxx, zfyy, zfxy,             &  !    "           " 
    276          zfm , zbet,                   &  !    "           " 
    277          zalg, zalg1, zalg1q              !    "           " 
     244      INTEGER  ::   ji, jj                               ! dummy loop indices 
     245      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars 
     246      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
     247      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     248      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
     249      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
     250      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    278251      !--------------------------------------------------------------------- 
    279252 
    … …  
    290263               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  ) 
    291264            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
     265            ! 
    292266            ps0 (ji,jj) = zslpmax   
    293             psx (ji,jj)  = psx (ji,jj) * zin0 
    294             psxx(ji,jj)  = psxx(ji,jj) * zin0 
     267            psx (ji,jj) = psx (ji,jj) * zin0 
     268            psxx(ji,jj) = psxx(ji,jj) * zin0 
    295269            psy (ji,jj) = zs1new * zin0 
    296270            psyy(ji,jj) = zs2new * zin0 
    … …  
    300274 
    301275      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    302       psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 
     276      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 
    303277 
    304278      !  Calculate fluxes and moments between boxes j<-->j+1               
    305       DO jj = 1, jpj 
    306          DO ji = 1, jpi 
    307             !  Flux from j to j+1 WHEN v GT 0    
     279      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0    
     280         DO ji = 1, jpi 
    308281            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) ) 
    309282            zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj) 
    … …  
    311284            zalf1        =  1.0 - zalf 
    312285            zalf1q       =  zalf1 * zalf1 
     286            ! 
    313287            zfm (ji,jj)  =  zalf  * psm(ji,jj) 
    314288            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) )  
    … …  
    318292            zfxy(ji,jj)  =  zalfq * psxy(ji,jj) 
    319293            zfxx(ji,jj)  =  zalf  * psxx(ji,jj) 
    320  
     294            ! 
    321295            !  Readjust moments remaining in the box. 
    322296            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj) 
    … …  
    329303         END DO 
    330304      END DO 
    331  
     305      ! 
    332306      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
    333307         DO ji = 1, jpi 
    … …  
    339313            zalf1q        = zalf1 * zalf1 
    340314            zalg1q(ji,jj) = zalf1q 
    341             zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji,jj+1) 
    342             zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 
    343             zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 
    344             zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * zalfq * psyy(ji,jj+1) 
    345             zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 
    346             zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1) 
    347             zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * psxx(ji,jj+1) 
     315            ! 
     316            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1) 
     317            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 
     318            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 
     319            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq 
     320            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 
     321            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1) 
     322            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1) 
    348323         END DO 
    349324      END DO 
    … …  
    354329            zbt  =         zbet(ji,jj-1) 
    355330            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
     331            ! 
    356332            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 
    357333            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 
    358334            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) ) 
    359             psyy(ji,jj) = zalg1 (ji,jj-1)  * zalg1q(ji,jj-1) * psyy(ji,jj) 
     335            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj) 
    360336            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 
    361337            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) 
    … …  
    373349            zalf1       = 1.0 - zalf 
    374350            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1) 
    375             ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj) 
    376  
     351            ! 
     352            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj) 
    377353            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   & 
    378                &        + zbt1 * psy(ji,jj)   
    379  
     354               &                                               + zbt1 * psy(ji,jj)   
    380355            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             & 
    381356               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   &  
    382                &        + zbt1 * psyy(ji,jj) 
    383  
    384             psxy(ji,jj) = zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)              & 
    385                + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) )   & 
    386                + zbt1 * psxy(ji,jj) 
     357               &                                               + zbt1 * psyy(ji,jj) 
     358            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               & 
     359               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   & 
     360               &                                                + zbt1 * psxy(ji,jj) 
    387361            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 
    388362            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) 
    … …  
    397371            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj) 
    398372            zalf1       = 1.0 - zalf 
    399             ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 
    400             ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
    401             psy(ji,jj)  = zbt  * psy(ji,jj)  & 
    402                &        + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 
    403             psyy(ji,jj) = zbt  * psyy(ji,jj)  & 
    404                &        + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   & 
    405                &                 + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) 
    406             psxy(ji,jj) = zbt  * psxy(ji,jj)   & 
    407                &        + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)   & 
    408                &                 + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) ) 
    409             psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 
    410             psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 
     373            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj) 
     374            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) ) 
     375            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 
     376            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   & 
     377               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          & 
     378               &                                         + ( zalf1 - zalf ) * ztemp )                                ) 
     379            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       & 
     380               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  ) 
     381            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 
     382            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 
    411383         END DO 
    412384      END DO 
    413385 
    414386      !-- Lateral boundary conditions 
    415       CALL lbc_lnk( psm , 'T',  1. ) 
    416       CALL lbc_lnk( ps0 , 'T',  1. ) 
    417       CALL lbc_lnk( psx , 'T', -1. ) 
    418       CALL lbc_lnk( psxx, 'T',  1. ) 
    419       CALL lbc_lnk( psy , 'T', -1. ) 
    420       CALL lbc_lnk( psyy, 'T',  1. ) 
     387      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. ) 
     388      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes 
     389      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. ) 
    421390      CALL lbc_lnk( psxy, 'T',  1. ) 
    422391 
    … …  
    427396         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :') 
    428397      ENDIF 
    429  
     398      ! 
    430399   END SUBROUTINE lim_adv_y 
    431400 
     401 
    432402#else 
    433403   !!---------------------------------------------------------------------- 
    434404   !!   Default option            Dummy module         NO LIM sea-ice model 
    435405   !!---------------------------------------------------------------------- 
    436 CONTAINS 
    437    SUBROUTINE lim_adv_x         ! Empty routine 
    438    END SUBROUTINE lim_adv_x 
    439    SUBROUTINE lim_adv_y           ! Empty routine 
    440    END SUBROUTINE lim_adv_y 
    441  
    442406#endif 
    443407 
     408   !!====================================================================== 
    444409END MODULE limadv 
Note: See TracChangeset for help on using the changeset viewer.