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

Changeset 6989


Ignore:
Timestamp:
2016-10-05T09:43:42+02:00 (7 years ago)
Author:
clem
Message:

use a namelist parameter to choose between the different advection schemes

Location:
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM
Files:
1 added
1 deleted
9 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r6796 r6989  
    7171&namicedyn     !   Ice dynamics and transport 
    7272!------------------------------------------------------------------------------ 
     73                  ! -- limtrp & limadv -- ! 
     74   nn_limadv      =    0            !  choose the advection scheme (-1=Prather ; 0=Ultimate-Macho) 
     75   nn_limadv_ord  =    5            !  choose the order of the advection scheme (if nn_limadv=0) 
    7376                  ! -- limitd_me -- ! 
    7477   nn_icestr      =    0            !  ice strength parameteriztaion                       
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r6796 r6989  
    200200 
    201201   !                                     !!** ice-dynamics namelist (namicedyn) ** 
     202                                          ! -- limtrp & limadv -- ! 
     203   INTEGER , PUBLIC ::   nn_limadv        !: choose the advection scheme (-1=Prather ; 0=Ultimate-Macho) 
     204   INTEGER , PUBLIC ::   nn_limadv_ord    !: choose the order of the advection scheme (if Ultimate-Macho)    
    202205                                          ! -- limitd_me -- ! 
    203206   INTEGER , PUBLIC ::   nn_icestr        !: ice strength parameterization (0=Hibler79 1=Rothrock75) 
     
    281284   REAL(wp), PUBLIC ::   rdt_ice          !: ice time step 
    282285   REAL(wp), PUBLIC ::   r1_rdtice        !: = 1. / rdt_ice 
    283    REAL(wp), PUBLIC ::   usecc2           !:  = 1.0 / ( rn_ecc * rn_ecc ) 
    284    REAL(wp), PUBLIC ::   rhoco            !: = rau0 * cio 
    285286   REAL(wp), PUBLIC ::   r1_nlay_i        !: 1 / nlay_i 
    286287   REAL(wp), PUBLIC ::   r1_nlay_s        !: 1 / nlay_s  
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90

    r6860 r6989  
    2929 
    3030   PUBLIC   lim_adv_umx    ! routine called by limtrp.F90 
    31     
    32    INTEGER  ::   nn_ult_order = 4       ! order of the ULTIMATE scheme 
    33     
    34    REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6 
     31       
     32   REAL(wp) ::   z1_6   = 1._wp / 6._wp   ! =1/6 
     33   REAL(wp) ::   z1_120 = 1._wp / 120._wp ! =1/120 
    3534 
    3635   !! * Substitutions 
     
    4342CONTAINS 
    4443 
    45    SUBROUTINE lim_adv_umx( lcon, kt, pdt, pu, pv, puc, pvc, pt, ptc, pfu, pfv ) 
     44   SUBROUTINE lim_adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc ) 
    4645      !!---------------------------------------------------------------------- 
    4746      !!                  ***  ROUTINE lim_adv_umx  *** 
     
    5655      !! ** Action : - pt  the after advective tracer 
    5756      !!---------------------------------------------------------------------- 
    58       LOGICAL                     , INTENT(in   )           ::   lcon       ! "continuity" equations for a and H  
    5957      INTEGER                     , INTENT(in   )           ::   kt         ! number of iteration 
    6058      REAL(wp)                    , INTENT(in   )           ::   pdt        ! tracer time-step 
    61       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pu, pv     ! 2 ice velocity components 
    62       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc, pvc   ! 2 ice velocity components  (for a or H) 
    63                                                                             ! 2 ice transport components (for tracers q) 
    64       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pt         ! tracer field 
     59      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc, pvc   ! 2 ice velocity components => u*e2 
     60      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity 
    6561      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc        ! tracer content field 
    66       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pfu, pfv   ! advective fluxes at u- and v-points  
    6762      ! 
    6863      INTEGER  ::   ji, jj           ! dummy loop indices   
     
    7772      CALL wrk_alloc( jpi,jpj,   zt_ups, zfu_ups, zfv_ups, ztrd, zfu_ho, zfv_ho, zt_u, zt_v ) 
    7873      ! 
    79       ! 
    80 !clem      zt_ups(:,:) = 0._wp 
    8174      ! 
    8275      !  upstream advection with initial mass fluxes & intermediate update 
     
    8881            zfp_vj = pvc(ji,jj) + ABS( pvc(ji,jj) ) 
    8982            zfm_vj = pvc(ji,jj) - ABS( pvc(ji,jj) ) 
    90             zfu_ups(ji,jj) = 0.5_wp * ( zfp_ui * pt(ji,jj) + zfm_ui * pt(ji+1,jj  ) ) 
    91             zfv_ups(ji,jj) = 0.5_wp * ( zfp_vj * pt(ji,jj) + zfm_vj * pt(ji  ,jj+1) ) 
     83            zfu_ups(ji,jj) = 0.5_wp * ( zfp_ui * ptc(ji,jj) + zfm_ui * ptc(ji+1,jj  ) ) 
     84            zfv_ups(ji,jj) = 0.5_wp * ( zfp_vj * ptc(ji,jj) + zfm_vj * ptc(ji  ,jj+1) ) 
    9285         END DO 
    9386      END DO 
     
    10699      ! High order (_ho) fluxes  
    107100      ! ----------------------- 
    108       SELECT CASE( nn_ult_order ) 
     101      SELECT CASE( nn_limadv_ord ) 
    109102      CASE ( 20 )                          ! centered second order 
    110103         DO jj = 2, jpjm1 
    111104            DO ji = fs_2, fs_jpim1   ! vector opt. 
    112                zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 
    113                zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 
    114             END DO 
    115          END DO 
    116          ! 
    117       CASE ( 1 , 4 )                      ! 1st to 4th order ULTIMATE-MACHO scheme 
    118          CALL macho( lcon, kt, nn_ult_order, pdt, pt, pu, pv, zt_u, zt_v ) 
    119 !!         CALL macho( lcon, kt, nn_ult_order, pdt, ptc, pu, pv, zt_u, zt_v ) 
     105               zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) ) 
     106               zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) ) 
     107            END DO 
     108         END DO 
     109         ! 
     110      CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme 
     111         CALL macho( kt, nn_limadv_ord, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 
    120112         ! 
    121113         DO jj = 2, jpjm1 
     
    153145      CALL lbc_lnk( ptc(:,:) , 'T',  1. ) 
    154146      ! 
    155       IF( PRESENT( pfu ) ) THEN 
    156          DO jj = 2, jpjm1 
    157             DO ji = fs_2, fs_jpim1   ! vector opt. 
    158                pfu(ji,jj) = zfu_ups(ji,jj) + zfu_ho(ji,jj) 
    159                pfv(ji,jj) = zfv_ups(ji,jj) + zfv_ho(ji,jj) 
    160             END DO 
    161          END DO 
    162          CALL lbc_lnk_multi( pfu, 'U', -1., pfv, 'V', -1. )         ! Lateral bondary conditions 
    163       ENDIF 
    164147      ! 
    165148      CALL wrk_dealloc( jpi,jpj,   zt_ups, zfu_ups, zfv_ups, ztrd, zfu_ho, zfv_ho, zt_u, zt_v ) 
     
    170153 
    171154 
    172    SUBROUTINE macho( lcon, kt, k_order, pdt, pt, pu, pv, pt_u, pt_v ) 
     155   SUBROUTINE macho( kt, k_order, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 
    173156      !!--------------------------------------------------------------------- 
    174157      !!                    ***  ROUTINE ultimate_x  *** 
     
    181164      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    182165      !!---------------------------------------------------------------------- 
    183       LOGICAL                     , INTENT(in   ) ::   lcon       ! "continuity" equations for a and ah  
    184166      INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration 
    185167      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme 
    186168      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step 
    187       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv     ! 2 ice velocity components 
    188       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt         ! tracer fields 
     169      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc        ! tracer fields 
     170      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components 
     171      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    189172      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points  
    190173      ! 
     
    201184         ! 
    202185         !                                                           !--  ultimate interpolation of pt at u-point  --! 
    203          CALL ultimate_x( lcon, k_order, pdt, pt, pu, pt_u ) 
    204          ! 
    205          !                                                           !--  advective form update in zzt  --! 
    206          DO jj = 1, jpj 
    207             DO ji = fs_2, fs_jpim1   ! vector opt. 
    208                ! 
    209                IF( pu(ji,jj) * pu(ji-1,jj) <= 0._wp ) THEN   ;   zc_box = 0._wp 
    210                ELSEIF(         pu(ji  ,jj) >  0._wp ) THEN   ;   zc_box = pu(ji-1,jj) 
    211                ELSE                                          ;   zc_box = pu(ji  ,jj) 
    212                ENDIF 
    213                ! 
    214                zzt(ji,jj) = pt(ji,jj) - zc_box * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e12t(ji,jj) 
    215                IF( lcon .eqv. .TRUE. ) THEN ! clem => for a.div(u) ?? 
    216                   zzt(ji,jj) = zzt(ji,jj) - pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e12t(ji,jj) 
    217                END IF 
    218             END DO 
    219          END DO 
    220          ! 
    221          !                                                           !--  ultimate interpolation of pt at v-point  --! 
    222          CALL ultimate_y( lcon, k_order, pdt, zzt, pv, pt_v ) 
    223          ! 
    224       ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==! 
    225          ! 
    226          !                                                           !--  ultimate interpolation of pt at v-point  --! 
    227          CALL ultimate_y( lcon, k_order, pdt, pt, pv, pt_v ) 
     186         CALL ultimate_x( k_order, pdt, ptc, puc, pt_u ) 
    228187         ! 
    229188         !                                                           !--  advective form update in zzt  --! 
    230189         DO jj = 2, jpjm1 
    231             DO ji = 1, jpi 
    232                ! 
    233                IF( pv(ji,jj) * pv(ji,jj-1) <= 0._wp ) THEN   ;   zc_box = 0._wp 
    234                ELSEIF(         pv(ji,jj  ) >  0._wp ) THEN   ;   zc_box = pv(ji,jj-1) 
    235                ELSE                                          ;   zc_box = pv(ji,jj  ) 
    236                ENDIF 
    237                ! 
    238                zzt(ji,jj) = pt(ji,jj) - zc_box * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e12t(ji,jj) 
    239                IF( lcon .eqv. .TRUE. ) THEN ! clem => for a.div(u) ?? 
    240                   zzt(ji,jj) = zzt(ji,jj) - pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e12t(ji,jj) 
    241                END IF 
    242             END DO 
    243          END DO 
     190            DO ji = fs_2, fs_jpim1   ! vector opt. 
     191               zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  & 
     192                  &                    - ptc  (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e12t(ji,jj) 
     193               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
     194            END DO 
     195         END DO 
     196         CALL lbc_lnk( zzt, 'T', 1. ) 
     197         ! 
     198         !                                                           !--  ultimate interpolation of pt at v-point  --! 
     199         CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v ) 
     200         ! 
     201      ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==! 
     202         ! 
     203         !                                                           !--  ultimate interpolation of pt at v-point  --! 
     204         CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v ) 
     205         ! 
     206         !                                                           !--  advective form update in zzt  --! 
     207         DO jj = 2, jpjm1 
     208            DO ji = fs_2, fs_jpim1 
     209               zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  & 
     210                  &                    - ptc  (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e12t(ji,jj) 
     211               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
     212            END DO 
     213         END DO 
     214         CALL lbc_lnk( zzt, 'T', 1. ) 
    244215         ! 
    245216         !                                                           !--  ultimate interpolation of pt at u-point  --! 
    246          CALL ultimate_x( lcon, k_order, pdt, zzt, pu, pt_u ) 
     217         CALL ultimate_x( k_order, pdt, zzt, puc, pt_u ) 
    247218         !       
    248219      ENDIF       
     
    255226 
    256227 
    257    SUBROUTINE ultimate_x( lcon, k_order, pdt, pt, pu, pt_u ) 
     228   SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u ) 
    258229      !!--------------------------------------------------------------------- 
    259230      !!                    ***  ROUTINE ultimate_x  *** 
     
    266237      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    267238      !!---------------------------------------------------------------------- 
    268       LOGICAL                     , INTENT(in   ) ::   lcon      ! "continuity" equations for a and ah  
    269239      INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index 
    270240      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
    271       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu        ! ice i-velocity component 
     241      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity component 
    272242      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields 
    273243      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point  
    274244      ! 
    275245      INTEGER  ::   ji, jj       ! dummy loop indices 
    276       REAL(wp) ::   zcu, zdx2    !   -      - 
    277       REAL(wp), POINTER, DIMENSION(:,:) :: ztu, zltu 
     246      REAL(wp) ::   zcu, zdx2, zdx4    !   -      - 
     247      REAL(wp), POINTER, DIMENSION(:,:) :: ztu1, ztu2, ztu3, ztu4 
    278248      !!---------------------------------------------------------------------- 
    279249      ! 
    280250      IF( nn_timing == 1 )  CALL timing_start('ultimate_x') 
    281251      ! 
    282       CALL wrk_alloc( jpi,jpj,   ztu, zltu ) 
     252      CALL wrk_alloc( jpi,jpj,   ztu1, ztu2, ztu3, ztu4 ) 
     253      ! 
     254      !                                                     !--  Laplacian in i-direction  --! 
     255      DO jj = 2, jpjm1         ! First derivative (gradient) 
     256         DO ji = 1, fs_jpim1 
     257            ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
     258         END DO 
     259         !                     ! Second derivative (Laplacian) 
     260         DO ji = fs_2, fs_jpim1 
     261            ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj) 
     262         END DO 
     263      END DO 
     264      CALL lbc_lnk( ztu2, 'T', 1. ) 
     265      ! 
     266      !                                                     !--  BiLaplacian in i-direction  --! 
     267      DO jj = 2, jpjm1         ! Third derivative 
     268         DO ji = 1, fs_jpim1 
     269            ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
     270         END DO 
     271         !                     ! Fourth derivative 
     272         DO ji = fs_2, fs_jpim1 
     273            ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj) 
     274         END DO 
     275      END DO 
     276      CALL lbc_lnk( ztu4, 'T', 1. ) 
     277      ! 
    283278      ! 
    284279      SELECT CASE (k_order ) 
     
    288283         DO jj = 1, jpj 
    289284            DO ji = 1, fs_jpim1   ! vector opt. 
    290                pt_u(ji,jj) = 0.5_wp * (                             ( pt(ji+1,jj) + pt(ji,jj) )    & 
    291                   &                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) )  ) * umask(ji,jj,1) 
     285               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj) + pt(ji,jj)   & 
     286                  &                                    - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
    292287            END DO 
    293288         END DO 
    294289         ! 
    295290      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23) 
     291         ! 
    296292         DO jj = 1, jpj 
    297293            DO ji = 1, fs_jpim1   ! vector opt. 
    298                pt_u(ji,jj) = 0.5_wp * (                    ( pt(ji+1,jj) + pt(ji,jj) )                 & 
    299 !!                  &                    - pdt * pu(ji,jj) * ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj)  ) * umask(ji,jj,1) 
    300                   &                    - pdt * pu(ji,jj) * ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e12u(ji,jj) ) * umask(ji,jj,1) 
     294               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     295               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   & 
     296                  &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) )  
    301297            END DO 
    302298         END DO 
     
    305301      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24) 
    306302         ! 
    307          !                                                                 !--  Laplacian in i- and in j-directions  --! 
    308          DO jj = 2, jpjm1         ! First derivative (gradient) 
    309             DO ji = 1, fs_jpim1   ! vector opt. 
    310                ztu(ji,jj) = ( pt(ji+1,jj  ) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
    311             END DO 
    312             !                     ! Second derivative (divergence) 
    313             DO ji = fs_2, fs_jpim1   ! vector opt. 
    314                zltu(ji,jj) = (  ztu(ji,jj) - ztu(ji-1,jj)  ) * r1_e1t(ji,jj) 
    315             END DO 
    316          END DO 
    317          CALL lbc_lnk( zltu, 'T', 1. )   ! Lateral boundary cond. (unchanged sign) 
    318          !     
    319303         DO jj = 1, jpj 
    320304            DO ji = 1, fs_jpim1   ! vector opt. 
    321 !!               zcu  = pu(ji,jj) * pdt * r1_e1u(ji,jj) 
    322                zcu  = pu(ji,jj) * pdt * r1_e12u(ji,jj) 
     305               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    323306               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     307!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    324308               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                         pt  (ji+1,jj) + pt  (ji,jj)        & 
    325309                  &                                               -              zcu   * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   & 
    326                   &        - r1_6 * zdx2 * ( 1._wp - zcu*zcu ) * (                         zltu(ji+1,jj) + zltu(ji,jj)        & 
    327                   &                                               - SIGN( 1._wp, zcu ) * ( zltu(ji+1,jj) - zltu(ji,jj) )  )   ) 
     310                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj) + ztu2(ji,jj)        & 
     311                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   ) 
    328312            END DO 
    329313         END DO 
    330314         ! 
    331315      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27) 
    332          ! 
    333          !                                                                 !--  Laplacian in i- and in j-directions  --! 
    334          DO jj = 1, jpjm1         ! First derivative (gradient) 
    335             DO ji = 1, fs_jpim1   ! vector opt. 
    336                ztu(ji,jj) = ( pt(ji+1,jj  ) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
    337             END DO 
    338             !                     ! Second derivative (divergence) 
    339             DO ji = fs_2, fs_jpim1   ! vector opt. 
    340                zltu(ji,jj) = (  ztu(ji,jj) - ztu(ji-1,jj)  ) * r1_e1t(ji,jj) 
    341             END DO 
    342          END DO 
    343          CALL lbc_lnk( zltu, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
    344316         ! 
    345317         DO jj = 1, jpj 
    346318            DO ji = 1, fs_jpim1   ! vector opt. 
    347 !!               zcu  = pu(ji,jj) * pdt * r1_e1u(ji,jj) 
    348                zcu  = pu(ji,jj) * pdt * r1_e12u(ji,jj) 
     319               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    349320               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     321!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    350322               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                   pt  (ji+1,jj) + pt  (ji,jj)        & 
    351323                  &                                               -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   & 
    352                   &        - r1_6 * zdx2 * ( 1._wp - zcu*zcu ) * (                   zltu(ji+1,jj) + zltu(ji,jj)        & 
    353                   &                                               - 0.5_wp * zcu * ( zltu(ji+1,jj) - zltu(ji,jj) )  )   ) 
     324                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj) + ztu2(ji,jj)        & 
     325                  &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   ) 
     326            END DO 
     327         END DO 
     328         ! 
     329      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29) 
     330         ! 
     331         DO jj = 1, jpj 
     332            DO ji = 1, fs_jpim1   ! vector opt. 
     333               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     334               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     335!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
     336               zdx4 = zdx2 * zdx2 
     337               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (               (                   pt  (ji+1,jj) + pt  (ji,jj)       & 
     338                  &                                                     -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) ) )   & 
     339                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj) + ztu2(ji,jj)       & 
     340                  &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) )   & 
     341                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj)       & 
     342                  &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) ) 
    354343            END DO 
    355344         END DO 
     
    357346      END SELECT 
    358347      ! 
    359       CALL wrk_dealloc( jpi,jpj,   ztu, zltu ) 
     348      CALL wrk_dealloc( jpi,jpj,   ztu1, ztu2, ztu3, ztu4 ) 
    360349      ! 
    361350      IF( nn_timing == 1 )  CALL timing_stop('ultimate_x') 
     
    364353    
    365354  
    366    SUBROUTINE ultimate_y( lcon, k_order, pdt, pt, pv, pt_v ) 
     355   SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v ) 
    367356      !!--------------------------------------------------------------------- 
    368357      !!                    ***  ROUTINE ultimate_y  *** 
     
    375364      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    376365      !!---------------------------------------------------------------------- 
    377       LOGICAL                     , INTENT(in   ) ::   lcon      ! "continuity" equations for a and ah  
    378366      INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index 
    379367      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
    380       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pv        ! ice j-velocity component 
     368      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity component 
    381369      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields 
    382370      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point  
    383371      ! 
    384372      INTEGER  ::   ji, jj       ! dummy loop indices 
    385       REAL(wp) ::   zcv, zdy2    !   -      - 
    386       REAL(wp), POINTER, DIMENSION(:,:) :: ztv, zltv 
     373      REAL(wp) ::   zcv, zdy2, zdy4    !   -      - 
     374      REAL(wp), POINTER, DIMENSION(:,:) :: ztv1, ztv2, ztv3, ztv4 
    387375      !!---------------------------------------------------------------------- 
    388376      ! 
    389377      IF( nn_timing == 1 )  CALL timing_start('ultimate_y') 
    390378      ! 
    391       CALL wrk_alloc( jpi,jpj,   ztv, zltv ) 
     379      CALL wrk_alloc( jpi,jpj,   ztv1, ztv2, ztv3, ztv4 ) 
     380      ! 
     381      !                                                     !--  Laplacian in j-direction  --! 
     382      DO jj = 1, jpjm1         ! First derivative (gradient) 
     383         DO ji = fs_2, fs_jpim1 
     384            ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
     385         END DO 
     386      END DO 
     387      DO jj = 2, jpjm1         ! Second derivative (Laplacian) 
     388         DO ji = fs_2, fs_jpim1 
     389            ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj) 
     390         END DO 
     391      END DO 
     392      CALL lbc_lnk( ztv2, 'T', 1. ) 
     393      ! 
     394      !                                                     !--  BiLaplacian in j-direction  --! 
     395      DO jj = 1, jpjm1         ! First derivative 
     396         DO ji = fs_2, fs_jpim1 
     397            ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
     398         END DO 
     399      END DO 
     400      DO jj = 2, jpjm1         ! Second derivative 
     401         DO ji = fs_2, fs_jpim1 
     402            ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj) 
     403         END DO 
     404      END DO 
     405      CALL lbc_lnk( ztv4, 'T', 1. ) 
     406      ! 
    392407      ! 
    393408      SELECT CASE (k_order ) 
     
    397412         DO jj = 1, jpjm1 
    398413            DO ji = 1, jpi 
    399                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  & 
    400                   &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
     414               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1) + pt(ji,jj) )  & 
     415                  &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
    401416            END DO 
    402417         END DO 
     
    405420         DO jj = 1, jpjm1 
    406421            DO ji = 1, jpi 
    407                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                    ( pt(ji,jj+1) + pt(ji,jj) )               & 
    408 !!                  &                                     - pdt * pv(ji,jj) * ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj)  ) 
    409                   &                                     - pdt * pv(ji,jj) * ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e12v(ji,jj) ) 
     422               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     423               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  & 
     424                  &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
    410425            END DO 
    411426         END DO 
     
    413428         ! 
    414429      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24) 
    415          ! 
    416          !                                                                 !--  Laplacian in i- and in j-directions  --! 
    417          DO jj = 1, jpjm1            ! First derivative (gradient) 
    418             DO ji = fs_2, fs_jpim1   ! vector opt. 
    419                ztv(ji,jj) = ( pt(ji  ,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    420             END DO 
    421          END DO 
    422          DO jj = 2, jpjm1            ! Second derivative (divergence) 
    423             DO ji = fs_2, fs_jpim1   ! vector opt. 
    424                zltv(ji,jj) = (  ztv(ji,jj) - ztv(ji,jj-1)  ) * r1_e2t(ji,jj) 
    425             END DO 
    426          END DO 
    427          CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sign) 
    428430         ! 
    429431         DO jj = 1, jpjm1 
    430432            DO ji = 1, jpi 
    431 !!               zcv  = pv(ji,jj) * pdt * r1_e2v(ji,jj) 
    432                zcv  = pv(ji,jj) * pdt * r1_e12v(ji,jj) 
     433               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    433434               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    434                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                (  pt  (ji,jj+1) + pt  (ji,jj)        & 
    435                   &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) )  )   & 
    436                   &        - r1_6 * zdy2 * ( 1._wp - zcv*zcv ) * (                         zltv(ji,jj+1) + zltv(ji,jj)        & 
    437                   &                                               - SIGN( 1._wp, zcv ) * ( zltv(ji,jj+1) - zltv(ji,jj) )  )   ) 
     435!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
     436               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)       & 
     437                  &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   & 
     438                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1) + ztv2(ji,jj)       & 
     439                  &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) 
    438440            END DO 
    439441         END DO 
    440442         ! 
    441443      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27) 
    442          ! 
    443          !                                                                 !--  Laplacian in i- and in j-directions  --! 
    444          DO jj = 1, jpjm1            ! First derivative (gradient) 
    445             DO ji = fs_2, fs_jpim1   ! vector opt. 
    446                ztv(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    447             END DO 
    448          END DO 
    449          DO jj = 2, jpjm1            ! Second derivative (divergence) 
    450             DO ji = fs_2, fs_jpim1   ! vector opt. 
    451                zltv(ji,jj) = (  ztv(ji,jj) - ztv(ji,jj-1)  ) * r1_e2t(ji,jj) 
    452             END DO 
    453          END DO 
    454          CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
    455444         ! 
    456445         DO jj = 1, jpjm1 
    457446            DO ji = 1, jpi 
    458 !!               zcv  = pv(ji,jj) * pdt * r1_e2v(ji,jj) 
    459                zcv  = pv(ji,jj) * pdt * r1_e12v(ji,jj) 
     447               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    460448               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    461                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                          (  pt  (ji,jj+1) + pt  (ji,jj)        & 
    462                   &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) )  )   & 
    463                   &        - r1_6 * zdy2 * ( 1._wp - zcv*zcv ) * (                   zltv(ji,jj+1) + zltv(ji,jj)        & 
    464                   &                                               - 0.5_wp * zcv * ( zltv(ji,jj+1) - zltv(ji,jj) )  )   ) 
     449!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
     450               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt  (ji,jj+1) + pt  (ji,jj)       & 
     451                  &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   & 
     452                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1) + ztv2(ji,jj)       & 
     453                  &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) 
     454            END DO 
     455         END DO 
     456         ! 
     457      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29) 
     458         ! 
     459         DO jj = 1, jpjm1 
     460            DO ji = 1, jpi 
     461               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     462               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     463!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
     464               zdy4 = zdy2 * zdy2 
     465               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)      & 
     466                  &                                                     -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )  & 
     467                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1) + ztv2(ji,jj)      & 
     468                  &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) )  & 
     469                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj)      & 
     470                  &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) ) 
    465471            END DO 
    466472         END DO 
     
    468474      END SELECT 
    469475      ! 
    470       CALL wrk_dealloc( jpi,jpj,   ztv, zltv ) 
     476      CALL wrk_dealloc( jpi,jpj,   ztv1, ztv2, ztv3, ztv4 ) 
    471477      ! 
    472478      IF( nn_timing == 1 )  CALL timing_stop('ultimate_y') 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r6796 r6989  
    2626   USE prtctl           ! Print control 
    2727   USE lib_fortran      ! glob_sum 
    28    USE timing          ! Timing 
    29    USE limcons        ! conservation tests 
     28   USE timing           ! Timing 
     29   USE limcons          ! conservation tests 
    3030   USE limvar 
    3131 
     
    3333   PRIVATE 
    3434 
    35    PUBLIC   lim_dyn   ! routine called by ice_step 
     35   PUBLIC   lim_dyn        ! routine called by sbcice_lim.F90 
     36   PUBLIC   lim_dyn_init   ! routine called by sbcice_lim.F90 
    3637 
    3738   !! * Substitutions 
     
    6465 
    6566      CALL lim_var_agg(1)                      ! aggregate ice categories 
    66  
    67       IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only) 
    6867      ! 
    6968      ! conservation test 
     
    147146      !!------------------------------------------------------------------- 
    148147      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    149       NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, & 
     148      NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord,  & 
     149         &                nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, & 
    150150         &                nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr, rn_lfrelax 
    151151      !!------------------------------------------------------------------- 
     
    164164         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 
    165165         WRITE(numout,*) '~~~~~~~~~~~~' 
     166         ! limtrp 
     167         WRITE(numout,*)'    choose the advection scheme (-1=Prather, 0=Ulimate-Macho)   nn_limadv     = ', nn_limadv  
     168         WRITE(numout,*)'    choose the order of the scheme (if ultimate)                nn_limadv_ord = ', nn_limadv_ord   
     169         ! limrhg 
    166170         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)         nn_icestr     = ', nn_icestr  
    167171         WRITE(numout,*)'    Including brine volume in ice strength comp.                ln_icestr_bvf = ', ln_icestr_bvf 
     
    180184      ENDIF 
    181185      ! 
    182       usecc2 = 1._wp / ( rn_ecc * rn_ecc ) 
    183       rhoco  = rau0  * rn_cio 
    184       ! 
    185186   END SUBROUTINE lim_dyn_init 
    186187 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r6853 r6989  
    114114      CHARACTER (len=50) ::   charout 
    115115 
     116      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio 
    116117      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling 
    117118      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity 
     
    203204      ! 1) define some variables and initialize arrays 
    204205      !------------------------------------------------------------------------------! 
     206      zrhoco = rau0 * rn_cio  
     207 
    205208      ! ecc2: square of yield ellipse eccenticrity 
    206209      ecc2    = rn_ecc * rn_ecc 
     
    357360                
    358361               ! delta at T points 
    359                zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )   
     362               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
    360363 
    361364               ! P/delta at T points 
     
    424427 
    425428                  ! tau_io/(v_oce - v_ice) 
    426                   zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    427                      &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     429                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     430                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    428431 
    429432                  ! tau_bottom/v_ice 
     
    469472                                
    470473                  ! tau_io/(u_oce - u_ice) 
    471                   zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    472                      &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     474                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     475                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    473476 
    474477                  ! tau_bottom/u_ice 
     
    515518 
    516519                  ! tau_io/(u_oce - u_ice) 
    517                   zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    518                      &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     520                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     521                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    519522 
    520523                  ! tau_bottom/u_ice 
     
    559562          
    560563                  ! tau_io/(v_oce - v_ice) 
    561                   zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    562                      &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     564                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     565                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    563566 
    564567                  ! tau_bottom/v_ice 
     
    651654             
    652655            ! delta at T points 
    653             zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )   
     656            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
    654657            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
    655658            delta_i(ji,jj) = zdelta + rn_creepl * rswitch 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90

    r5341 r6989  
    5353      INTEGER, INTENT(in) ::   kt       ! number of iteration 
    5454      ! 
    55       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character 
    56       CHARACTER(LEN=50)   ::   clname   ! ice output restart file name 
     55      CHARACTER(len=20)   ::   clkt     ! ocean time-step define as a character 
     56      CHARACTER(len=50)   ::   clname   ! ice output restart file name 
    5757      CHARACTER(len=256)  ::   clpath   ! full path to ice output restart file  
    5858      !!---------------------------------------------------------------------- 
     
    107107      INTEGER ::   ji, jj, jk ,jl   ! dummy loop indices 
    108108      INTEGER ::   iter 
    109       CHARACTER(len=15) ::   znam 
    110       CHARACTER(len=1)  ::   zchar, zchar1 
     109      CHARACTER(len=25) ::   znam 
     110      CHARACTER(len=2)  ::   zchar, zchar1 
    111111      REAL(wp), POINTER, DIMENSION(:,:) :: z2d 
    112112      !!---------------------------------------------------------------------- 
     
    130130      ! Prognostic variables  
    131131      DO jl = 1, jpl  
    132          WRITE(zchar,'(I1)') jl 
     132         WRITE(zchar,'(I2.2)') jl 
    133133         znam = 'v_i'//'_htc'//zchar 
    134134         z2d(:,:) = v_i(:,:,jl) 
     
    152152 
    153153      DO jl = 1, jpl  
    154          WRITE(zchar,'(I1)') jl 
     154         WRITE(zchar,'(I2.2)') jl 
    155155         znam = 'tempt_sl1'//'_htc'//zchar 
    156156         z2d(:,:) = e_s(:,:,1,jl) 
     
    159159 
    160160      DO jl = 1, jpl  
    161          WRITE(zchar,'(I1)') jl 
     161         WRITE(zchar,'(I2.2)') jl 
    162162         DO jk = 1, nlay_i  
    163             WRITE(zchar1,'(I1)') jk 
     163            WRITE(zchar1,'(I2.2)') jk 
    164164            znam = 'tempt'//'_il'//zchar1//'_htc'//zchar 
    165165            z2d(:,:) = e_i(:,:,jk,jl) 
     
    176176      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) 
    177177 
    178       DO jl = 1, jpl  
    179          WRITE(zchar,'(I1)') jl 
    180          znam = 'sxice'//'_htc'//zchar 
    181          z2d(:,:) = sxice(:,:,jl) 
    182          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    183          znam = 'syice'//'_htc'//zchar 
    184          z2d(:,:) = syice(:,:,jl) 
    185          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    186          znam = 'sxxice'//'_htc'//zchar 
    187          z2d(:,:) = sxxice(:,:,jl) 
    188          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    189          znam = 'syyice'//'_htc'//zchar 
    190          z2d(:,:) = syyice(:,:,jl) 
    191          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    192          znam = 'sxyice'//'_htc'//zchar 
    193          z2d(:,:) = sxyice(:,:,jl) 
    194          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    195          znam = 'sxsn'//'_htc'//zchar 
    196          z2d(:,:) = sxsn(:,:,jl) 
    197          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    198          znam = 'sysn'//'_htc'//zchar 
    199          z2d(:,:) = sysn(:,:,jl) 
    200          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    201          znam = 'sxxsn'//'_htc'//zchar 
    202          z2d(:,:) = sxxsn(:,:,jl) 
    203          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    204          znam = 'syysn'//'_htc'//zchar 
    205          z2d(:,:) = syysn(:,:,jl) 
    206          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    207          znam = 'sxysn'//'_htc'//zchar 
    208          z2d(:,:) = sxysn(:,:,jl) 
    209          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    210          znam = 'sxa'//'_htc'//zchar 
    211          z2d(:,:) = sxa(:,:,jl) 
    212          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    213          znam = 'sya'//'_htc'//zchar 
    214          z2d(:,:) = sya(:,:,jl) 
    215          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    216          znam = 'sxxa'//'_htc'//zchar 
    217          z2d(:,:) = sxxa(:,:,jl) 
    218          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    219          znam = 'syya'//'_htc'//zchar 
    220          z2d(:,:) = syya(:,:,jl) 
    221          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    222          znam = 'sxya'//'_htc'//zchar 
    223          z2d(:,:) = sxya(:,:,jl) 
    224          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    225          znam = 'sxc0'//'_htc'//zchar 
    226          z2d(:,:) = sxc0(:,:,jl) 
    227          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    228          znam = 'syc0'//'_htc'//zchar 
    229          z2d(:,:) = syc0(:,:,jl) 
    230          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    231          znam = 'sxxc0'//'_htc'//zchar 
    232          z2d(:,:) = sxxc0(:,:,jl) 
    233          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    234          znam = 'syyc0'//'_htc'//zchar 
    235          z2d(:,:) = syyc0(:,:,jl) 
    236          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    237          znam = 'sxyc0'//'_htc'//zchar 
    238          z2d(:,:) = sxyc0(:,:,jl) 
    239          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    240          znam = 'sxsal'//'_htc'//zchar 
    241          z2d(:,:) = sxsal(:,:,jl) 
    242          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    243          znam = 'sysal'//'_htc'//zchar 
    244          z2d(:,:) = sysal(:,:,jl) 
    245          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    246          znam = 'sxxsal'//'_htc'//zchar 
    247          z2d(:,:) = sxxsal(:,:,jl) 
    248          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    249          znam = 'syysal'//'_htc'//zchar 
    250          z2d(:,:) = syysal(:,:,jl) 
    251          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    252          znam = 'sxysal'//'_htc'//zchar 
    253          z2d(:,:) = sxysal(:,:,jl) 
    254          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    255          znam = 'sxage'//'_htc'//zchar 
    256          z2d(:,:) = sxage(:,:,jl) 
    257          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    258          znam = 'syage'//'_htc'//zchar 
    259          z2d(:,:) = syage(:,:,jl) 
    260          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    261          znam = 'sxxage'//'_htc'//zchar 
    262          z2d(:,:) = sxxage(:,:,jl) 
    263          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    264          znam = 'syyage'//'_htc'//zchar 
    265          z2d(:,:) = syyage(:,:,jl) 
    266          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    267          znam = 'sxyage'//'_htc'//zchar 
    268          z2d(:,:) = sxyage(:,:,jl) 
    269          CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    270       END DO 
    271  
    272       CALL iom_rstput( iter, nitrst, numriw, 'sxopw ' ,  sxopw  ) 
    273       CALL iom_rstput( iter, nitrst, numriw, 'syopw ' ,  syopw  ) 
    274       CALL iom_rstput( iter, nitrst, numriw, 'sxxopw' ,  sxxopw ) 
    275       CALL iom_rstput( iter, nitrst, numriw, 'syyopw' ,  syyopw ) 
    276       CALL iom_rstput( iter, nitrst, numriw, 'sxyopw' ,  sxyopw ) 
    277  
    278       DO jl = 1, jpl  
    279          WRITE(zchar,'(I1)') jl 
    280          DO jk = 1, nlay_i  
    281             WRITE(zchar1,'(I1)') jk 
    282             znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 
    283             z2d(:,:) = sxe(:,:,jk,jl) 
    284             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    285             znam = 'sye'//'_il'//zchar1//'_htc'//zchar 
    286             z2d(:,:) = sye(:,:,jk,jl) 
    287             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    288             znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 
    289             z2d(:,:) = sxxe(:,:,jk,jl) 
    290             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    291             znam = 'syye'//'_il'//zchar1//'_htc'//zchar 
    292             z2d(:,:) = syye(:,:,jk,jl) 
    293             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    294             znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 
    295             z2d(:,:) = sxye(:,:,jk,jl) 
    296             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    297          END DO 
    298       END DO 
    299  
     178      ! In case Prather scheme is used for advection, write second order moments 
     179      ! ------------------------------------------------------------------------ 
     180      IF( nn_limadv == -1 ) THEN 
     181          
     182         DO jl = 1, jpl  
     183            WRITE(zchar,'(I2.2)') jl 
     184            znam = 'sxice'//'_htc'//zchar 
     185            z2d(:,:) = sxice(:,:,jl) 
     186            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     187            znam = 'syice'//'_htc'//zchar 
     188            z2d(:,:) = syice(:,:,jl) 
     189            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     190            znam = 'sxxice'//'_htc'//zchar 
     191            z2d(:,:) = sxxice(:,:,jl) 
     192            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     193            znam = 'syyice'//'_htc'//zchar 
     194            z2d(:,:) = syyice(:,:,jl) 
     195            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     196            znam = 'sxyice'//'_htc'//zchar 
     197            z2d(:,:) = sxyice(:,:,jl) 
     198            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     199            znam = 'sxsn'//'_htc'//zchar 
     200            z2d(:,:) = sxsn(:,:,jl) 
     201            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     202            znam = 'sysn'//'_htc'//zchar 
     203            z2d(:,:) = sysn(:,:,jl) 
     204            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     205            znam = 'sxxsn'//'_htc'//zchar 
     206            z2d(:,:) = sxxsn(:,:,jl) 
     207            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     208            znam = 'syysn'//'_htc'//zchar 
     209            z2d(:,:) = syysn(:,:,jl) 
     210            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     211            znam = 'sxysn'//'_htc'//zchar 
     212            z2d(:,:) = sxysn(:,:,jl) 
     213            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     214            znam = 'sxa'//'_htc'//zchar 
     215            z2d(:,:) = sxa(:,:,jl) 
     216            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     217            znam = 'sya'//'_htc'//zchar 
     218            z2d(:,:) = sya(:,:,jl) 
     219            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     220            znam = 'sxxa'//'_htc'//zchar 
     221            z2d(:,:) = sxxa(:,:,jl) 
     222            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     223            znam = 'syya'//'_htc'//zchar 
     224            z2d(:,:) = syya(:,:,jl) 
     225            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     226            znam = 'sxya'//'_htc'//zchar 
     227            z2d(:,:) = sxya(:,:,jl) 
     228            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     229            znam = 'sxc0'//'_htc'//zchar 
     230            z2d(:,:) = sxc0(:,:,jl) 
     231            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     232            znam = 'syc0'//'_htc'//zchar 
     233            z2d(:,:) = syc0(:,:,jl) 
     234            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     235            znam = 'sxxc0'//'_htc'//zchar 
     236            z2d(:,:) = sxxc0(:,:,jl) 
     237            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     238            znam = 'syyc0'//'_htc'//zchar 
     239            z2d(:,:) = syyc0(:,:,jl) 
     240            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     241            znam = 'sxyc0'//'_htc'//zchar 
     242            z2d(:,:) = sxyc0(:,:,jl) 
     243            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     244            znam = 'sxsal'//'_htc'//zchar 
     245            z2d(:,:) = sxsal(:,:,jl) 
     246            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     247            znam = 'sysal'//'_htc'//zchar 
     248            z2d(:,:) = sysal(:,:,jl) 
     249            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     250            znam = 'sxxsal'//'_htc'//zchar 
     251            z2d(:,:) = sxxsal(:,:,jl) 
     252            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     253            znam = 'syysal'//'_htc'//zchar 
     254            z2d(:,:) = syysal(:,:,jl) 
     255            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     256            znam = 'sxysal'//'_htc'//zchar 
     257            z2d(:,:) = sxysal(:,:,jl) 
     258            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     259            znam = 'sxage'//'_htc'//zchar 
     260            z2d(:,:) = sxage(:,:,jl) 
     261            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     262            znam = 'syage'//'_htc'//zchar 
     263            z2d(:,:) = syage(:,:,jl) 
     264            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     265            znam = 'sxxage'//'_htc'//zchar 
     266            z2d(:,:) = sxxage(:,:,jl) 
     267            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     268            znam = 'syyage'//'_htc'//zchar 
     269            z2d(:,:) = syyage(:,:,jl) 
     270            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     271            znam = 'sxyage'//'_htc'//zchar 
     272            z2d(:,:) = sxyage(:,:,jl) 
     273            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     274         END DO 
     275 
     276         CALL iom_rstput( iter, nitrst, numriw, 'sxopw ' ,  sxopw  ) 
     277         CALL iom_rstput( iter, nitrst, numriw, 'syopw ' ,  syopw  ) 
     278         CALL iom_rstput( iter, nitrst, numriw, 'sxxopw' ,  sxxopw ) 
     279         CALL iom_rstput( iter, nitrst, numriw, 'syyopw' ,  syyopw ) 
     280         CALL iom_rstput( iter, nitrst, numriw, 'sxyopw' ,  sxyopw ) 
     281          
     282         DO jl = 1, jpl  
     283            WRITE(zchar,'(I2.2)') jl 
     284            DO jk = 1, nlay_i  
     285               WRITE(zchar1,'(I2.2)') jk 
     286               znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 
     287               z2d(:,:) = sxe(:,:,jk,jl) 
     288               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     289               znam = 'sye'//'_il'//zchar1//'_htc'//zchar 
     290               z2d(:,:) = sye(:,:,jk,jl) 
     291               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     292               znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 
     293               z2d(:,:) = sxxe(:,:,jk,jl) 
     294               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     295               znam = 'syye'//'_il'//zchar1//'_htc'//zchar 
     296               z2d(:,:) = syye(:,:,jk,jl) 
     297               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     298               znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 
     299               z2d(:,:) = sxye(:,:,jk,jl) 
     300               CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     301            END DO 
     302         END DO 
     303 
     304      ENDIF 
     305       
     306      ! close restart file 
     307      ! ------------------ 
    300308      IF( iter == nitrst ) THEN 
    301          CALL iom_close( numriw )                         ! close the restart file 
     309         CALL iom_close( numriw ) 
    302310         lrst_ice = .FALSE. 
    303311      ENDIF 
     
    317325      REAL(wp) ::   zfice, ziter 
    318326      REAL(wp), POINTER, DIMENSION(:,:) ::   z2d 
    319       CHARACTER(len=15) ::   znam 
    320       CHARACTER(len=1)  ::   zchar, zchar1 
     327      CHARACTER(len=25) ::   znam 
     328      CHARACTER(len=2)  ::   zchar, zchar1 
    321329      INTEGER           ::   jlibalt = jprstlib 
    322330      LOGICAL           ::   llok 
     
    356364         &                   '   control of time parameter  nrstdt' ) 
    357365 
     366      ! Prognostic variables  
    358367      DO jl = 1, jpl  
    359          WRITE(zchar,'(I1)') jl 
     368         WRITE(zchar,'(I2.2)') jl 
    360369         znam = 'v_i'//'_htc'//zchar 
    361370         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     
    379388 
    380389      DO jl = 1, jpl  
    381          WRITE(zchar,'(I1)') jl 
     390         WRITE(zchar,'(I2.2)') jl 
    382391         znam = 'tempt_sl1'//'_htc'//zchar 
    383392         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     
    386395 
    387396      DO jl = 1, jpl  
    388          WRITE(zchar,'(I1)') jl 
     397         WRITE(zchar,'(I2.2)') jl 
    389398         DO jk = 1, nlay_i  
    390             WRITE(zchar1,'(I1)') jk 
     399            WRITE(zchar1,'(I2.2)') jk 
    391400            znam = 'tempt'//'_il'//zchar1//'_htc'//zchar 
    392401            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     
    403412      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) 
    404413 
    405       DO jl = 1, jpl  
    406          WRITE(zchar,'(I1)') jl 
    407          znam = 'sxice'//'_htc'//zchar 
    408          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    409          sxice(:,:,jl) = z2d(:,:) 
    410          znam = 'syice'//'_htc'//zchar 
    411          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    412          syice(:,:,jl) = z2d(:,:) 
    413          znam = 'sxxice'//'_htc'//zchar 
    414          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    415          sxxice(:,:,jl) = z2d(:,:) 
    416          znam = 'syyice'//'_htc'//zchar 
    417          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    418          syyice(:,:,jl) = z2d(:,:) 
    419          znam = 'sxyice'//'_htc'//zchar 
    420          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    421          sxyice(:,:,jl) = z2d(:,:) 
    422          znam = 'sxsn'//'_htc'//zchar 
    423          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    424          sxsn(:,:,jl) = z2d(:,:) 
    425          znam = 'sysn'//'_htc'//zchar 
    426          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    427          sysn(:,:,jl) = z2d(:,:) 
    428          znam = 'sxxsn'//'_htc'//zchar 
    429          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    430          sxxsn(:,:,jl) = z2d(:,:) 
    431          znam = 'syysn'//'_htc'//zchar 
    432          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    433          syysn(:,:,jl) = z2d(:,:) 
    434          znam = 'sxysn'//'_htc'//zchar 
    435          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    436          sxysn(:,:,jl) = z2d(:,:) 
    437          znam = 'sxa'//'_htc'//zchar 
    438          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    439          sxa(:,:,jl) = z2d(:,:) 
    440          znam = 'sya'//'_htc'//zchar 
    441          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    442          sya(:,:,jl) = z2d(:,:) 
    443          znam = 'sxxa'//'_htc'//zchar 
    444          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    445          sxxa(:,:,jl) = z2d(:,:) 
    446          znam = 'syya'//'_htc'//zchar 
    447          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    448          syya(:,:,jl) = z2d(:,:) 
    449          znam = 'sxya'//'_htc'//zchar 
    450          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    451          sxya(:,:,jl) = z2d(:,:) 
    452          znam = 'sxc0'//'_htc'//zchar 
    453          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    454          sxc0(:,:,jl) = z2d(:,:) 
    455          znam = 'syc0'//'_htc'//zchar 
    456          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    457          syc0(:,:,jl) = z2d(:,:) 
    458          znam = 'sxxc0'//'_htc'//zchar 
    459          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    460          sxxc0(:,:,jl) = z2d(:,:) 
    461          znam = 'syyc0'//'_htc'//zchar 
    462          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    463          syyc0(:,:,jl) = z2d(:,:) 
    464          znam = 'sxyc0'//'_htc'//zchar 
    465          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    466          sxyc0(:,:,jl) = z2d(:,:) 
    467          znam = 'sxsal'//'_htc'//zchar 
    468          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    469          sxsal(:,:,jl) = z2d(:,:) 
    470          znam = 'sysal'//'_htc'//zchar 
    471          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    472          sysal(:,:,jl) = z2d(:,:) 
    473          znam = 'sxxsal'//'_htc'//zchar 
    474          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    475          sxxsal(:,:,jl) = z2d(:,:) 
    476          znam = 'syysal'//'_htc'//zchar 
    477          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    478          syysal(:,:,jl) = z2d(:,:) 
    479          znam = 'sxysal'//'_htc'//zchar 
    480          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    481          sxysal(:,:,jl) = z2d(:,:) 
    482          znam = 'sxage'//'_htc'//zchar 
    483          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    484          sxage(:,:,jl) = z2d(:,:) 
    485          znam = 'syage'//'_htc'//zchar 
    486          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    487          syage(:,:,jl) = z2d(:,:) 
    488          znam = 'sxxage'//'_htc'//zchar 
    489          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    490          sxxage(:,:,jl) = z2d(:,:) 
    491          znam = 'syyage'//'_htc'//zchar 
    492          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    493          syyage(:,:,jl) = z2d(:,:) 
    494          znam = 'sxyage'//'_htc'//zchar 
    495          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    496          sxyage(:,:,jl)= z2d(:,:) 
    497       END DO 
    498  
    499       CALL iom_get( numrir, jpdom_autoglo, 'sxopw ' ,  sxopw  ) 
    500       CALL iom_get( numrir, jpdom_autoglo, 'syopw ' ,  syopw  ) 
    501       CALL iom_get( numrir, jpdom_autoglo, 'sxxopw' ,  sxxopw ) 
    502       CALL iom_get( numrir, jpdom_autoglo, 'syyopw' ,  syyopw ) 
    503       CALL iom_get( numrir, jpdom_autoglo, 'sxyopw' ,  sxyopw ) 
    504  
    505       DO jl = 1, jpl  
    506          WRITE(zchar,'(I1)') jl 
    507          DO jk = 1, nlay_i  
    508             WRITE(zchar1,'(I1)') jk 
    509             znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 
    510             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    511             sxe(:,:,jk,jl) = z2d(:,:) 
    512             znam = 'sye'//'_il'//zchar1//'_htc'//zchar 
    513             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    514             sye(:,:,jk,jl) = z2d(:,:) 
    515             znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 
    516             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    517             sxxe(:,:,jk,jl) = z2d(:,:) 
    518             znam = 'syye'//'_il'//zchar1//'_htc'//zchar 
    519             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    520             syye(:,:,jk,jl) = z2d(:,:) 
    521             znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 
    522             CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    523             sxye(:,:,jk,jl) = z2d(:,:) 
    524          END DO 
    525       END DO 
    526       ! 
     414      ! In case Prather scheme is used for advection, read second order moments 
     415      ! ------------------------------------------------------------------------ 
     416      IF( nn_limadv == -1 ) THEN 
     417 
     418         DO jl = 1, jpl  
     419            WRITE(zchar,'(I2.2)') jl 
     420            znam = 'sxice'//'_htc'//zchar 
     421            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     422            sxice(:,:,jl) = z2d(:,:) 
     423            znam = 'syice'//'_htc'//zchar 
     424            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     425            syice(:,:,jl) = z2d(:,:) 
     426            znam = 'sxxice'//'_htc'//zchar 
     427            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     428            sxxice(:,:,jl) = z2d(:,:) 
     429            znam = 'syyice'//'_htc'//zchar 
     430            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     431            syyice(:,:,jl) = z2d(:,:) 
     432            znam = 'sxyice'//'_htc'//zchar 
     433            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     434            sxyice(:,:,jl) = z2d(:,:) 
     435            znam = 'sxsn'//'_htc'//zchar 
     436            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     437            sxsn(:,:,jl) = z2d(:,:) 
     438            znam = 'sysn'//'_htc'//zchar 
     439            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     440            sysn(:,:,jl) = z2d(:,:) 
     441            znam = 'sxxsn'//'_htc'//zchar 
     442            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     443            sxxsn(:,:,jl) = z2d(:,:) 
     444            znam = 'syysn'//'_htc'//zchar 
     445            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     446            syysn(:,:,jl) = z2d(:,:) 
     447            znam = 'sxysn'//'_htc'//zchar 
     448            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     449            sxysn(:,:,jl) = z2d(:,:) 
     450            znam = 'sxa'//'_htc'//zchar 
     451            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     452            sxa(:,:,jl) = z2d(:,:) 
     453            znam = 'sya'//'_htc'//zchar 
     454            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     455            sya(:,:,jl) = z2d(:,:) 
     456            znam = 'sxxa'//'_htc'//zchar 
     457            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     458            sxxa(:,:,jl) = z2d(:,:) 
     459            znam = 'syya'//'_htc'//zchar 
     460            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     461            syya(:,:,jl) = z2d(:,:) 
     462            znam = 'sxya'//'_htc'//zchar 
     463            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     464            sxya(:,:,jl) = z2d(:,:) 
     465            znam = 'sxc0'//'_htc'//zchar 
     466            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     467            sxc0(:,:,jl) = z2d(:,:) 
     468            znam = 'syc0'//'_htc'//zchar 
     469            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     470            syc0(:,:,jl) = z2d(:,:) 
     471            znam = 'sxxc0'//'_htc'//zchar 
     472            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     473            sxxc0(:,:,jl) = z2d(:,:) 
     474            znam = 'syyc0'//'_htc'//zchar 
     475            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     476            syyc0(:,:,jl) = z2d(:,:) 
     477            znam = 'sxyc0'//'_htc'//zchar 
     478            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     479            sxyc0(:,:,jl) = z2d(:,:) 
     480            znam = 'sxsal'//'_htc'//zchar 
     481            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     482            sxsal(:,:,jl) = z2d(:,:) 
     483            znam = 'sysal'//'_htc'//zchar 
     484            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     485            sysal(:,:,jl) = z2d(:,:) 
     486            znam = 'sxxsal'//'_htc'//zchar 
     487            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     488            sxxsal(:,:,jl) = z2d(:,:) 
     489            znam = 'syysal'//'_htc'//zchar 
     490            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     491            syysal(:,:,jl) = z2d(:,:) 
     492            znam = 'sxysal'//'_htc'//zchar 
     493            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     494            sxysal(:,:,jl) = z2d(:,:) 
     495            znam = 'sxage'//'_htc'//zchar 
     496            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     497            sxage(:,:,jl) = z2d(:,:) 
     498            znam = 'syage'//'_htc'//zchar 
     499            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     500            syage(:,:,jl) = z2d(:,:) 
     501            znam = 'sxxage'//'_htc'//zchar 
     502            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     503            sxxage(:,:,jl) = z2d(:,:) 
     504            znam = 'syyage'//'_htc'//zchar 
     505            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     506            syyage(:,:,jl) = z2d(:,:) 
     507            znam = 'sxyage'//'_htc'//zchar 
     508            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     509            sxyage(:,:,jl)= z2d(:,:) 
     510         END DO 
     511 
     512         CALL iom_get( numrir, jpdom_autoglo, 'sxopw ' ,  sxopw  ) 
     513         CALL iom_get( numrir, jpdom_autoglo, 'syopw ' ,  syopw  ) 
     514         CALL iom_get( numrir, jpdom_autoglo, 'sxxopw' ,  sxxopw ) 
     515         CALL iom_get( numrir, jpdom_autoglo, 'syyopw' ,  syyopw ) 
     516         CALL iom_get( numrir, jpdom_autoglo, 'sxyopw' ,  sxyopw ) 
     517 
     518         DO jl = 1, jpl  
     519            WRITE(zchar,'(I2.2)') jl 
     520            DO jk = 1, nlay_i  
     521               WRITE(zchar1,'(I2.2)') jk 
     522               znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 
     523               CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     524               sxe(:,:,jk,jl) = z2d(:,:) 
     525               znam = 'sye'//'_il'//zchar1//'_htc'//zchar 
     526               CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     527               sye(:,:,jk,jl) = z2d(:,:) 
     528               znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 
     529               CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     530               sxxe(:,:,jk,jl) = z2d(:,:) 
     531               znam = 'syye'//'_il'//zchar1//'_htc'//zchar 
     532               CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     533               syye(:,:,jk,jl) = z2d(:,:) 
     534               znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 
     535               CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     536               sxye(:,:,jk,jl) = z2d(:,:) 
     537            END DO 
     538         END DO 
     539         ! 
     540      END IF 
     541       
    527542      ! clem: I do not understand why the following IF is needed 
    528543      !       I suspect something inconsistent in the main code with option nn_icesal=1 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r6970 r6989  
    260260      INTEGER  ::   ji, jj   ! dummy loop indices 
    261261      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar 
    262       REAL(wp) ::   zat_v, zvtau_ice, zv_t          !   -      - 
     262      REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      - 
    263263      !!--------------------------------------------------------------------- 
     264      zrhoco = rau0 * rn_cio 
    264265      ! 
    265266      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step) 
     
    272273               zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  ) 
    273274               !                                               ! update the ocean stress modulus 
    274                taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zmodt 
    275                tmod_io(ji,jj) = rhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point 
     275               taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * zrhoco * zmodt 
     276               tmod_io(ji,jj) = zrhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point 
    276277            END DO 
    277278         END DO 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r6746 r6989  
    1818   USE sbc_oce        ! ocean surface boundary condition 
    1919   USE ice            ! ice variables 
    20    USE limadv         ! ice advection 
    2120   USE limhdf         ! ice horizontal diffusion 
    2221   USE limvar         !  
     22   USE limadv_prather ! advection scheme (Prather) 
     23   USE limadv_umx     ! advection scheme (ultimate-macho) 
    2324   ! 
    2425   USE in_out_manager ! I/O manager 
     
    3132   USE limcons        ! conservation tests 
    3233   USE limctl         ! control prints 
    33    USE limadv_umx     ! advection scheme 
    3434 
    3535   IMPLICIT NONE 
     
    5757      !! ** method  : variables included in the process are scalar,    
    5858      !!     other values are considered as second order.  
    59       !!     For advection, a second order Prather scheme is used.   
     59      !!     For advection, one can choose between 
     60      !!     a) an Ultimate-Macho scheme (whose order is defined by nn_limadv_ord) => nn_limadv=0 
     61      !!     b) and a second order Prather scheme => nn_limadv=-1 
    6062      !! 
    6163      !! ** action : 
     
    7981      ! --- ultimate macho only --- ! 
    8082      REAL(wp)                               ::   zdt 
    81       LOGICAL                                ::   lcon 
    82       REAL(wp), POINTER, DIMENSION(:,:)      ::   ze, zu_trp, zv_trp, z1_v, zudy, zvdx 
     83      REAL(wp), POINTER, DIMENSION(:,:)      ::   zudy, zvdx, zcu_box, zcv_box 
    8384      ! --- prather only --- ! 
    8485      REAL(wp), POINTER, DIMENSION(:,:)      ::   zarea 
     
    142143      zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )    ! CFL test for stability 
    143144      zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
    144       IF(lk_mpp )   CALL mpp_max( zcfl ) 
     145      IF( lk_mpp )   CALL mpp_max( zcfl ) 
    145146       
    146147      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
     
    156157!!      ENDIF 
    157158 
    158 #if defined key_limumx 
    159       !=============================! 
    160       !==  Ultimate-MACHO scheme  ==! 
    161       !=============================! 
    162        
    163       CALL wrk_alloc( jpi,jpj, ze, zu_trp, zv_trp, z1_v, zudy, zvdx ) 
    164        
    165       IF( kt == nit000 .AND. lwp ) THEN 
    166          WRITE(numout,*)'' 
    167          WRITE(numout,*)'lim_adv_umx : Ultimate-MACHO advection scheme' 
    168          WRITE(numout,*)'~~~~~~~~~~~' 
    169       ENDIF 
    170       ! 
    171       zdt = rdt_ice / REAL(initad) 
    172        
    173       ! transport 
    174       zudy(:,:) = u_ice(:,:) * e2u(:,:) 
    175       zvdx(:,:) = v_ice(:,:) * e1v(:,:) 
    176       ! 
    177       DO jt = 1, initad 
    178          lcon = .TRUE.  
    179 !!!            lcon = .false.  
    180          CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, ato_i(:,:), ato_i(:,:) )                                       ! Open water area  
     159      SELECT CASE ( nn_limadv ) 
     160          
     161                       !=============================! 
     162      CASE ( 0 )       !==  Ultimate-MACHO scheme  ==!                    
     163                       !=============================! 
     164       
     165         CALL wrk_alloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box ) 
     166       
     167         IF( kt == nit000 .AND. lwp ) THEN 
     168            WRITE(numout,*)'' 
     169            WRITE(numout,*)'lim_adv_umx : Ultimate-MACHO advection scheme' 
     170            WRITE(numout,*)'~~~~~~~~~~~' 
     171         ENDIF 
    181172         ! 
     173         zdt = rdt_ice / REAL(initad) 
     174          
     175         ! transport 
     176         zudy(:,:) = u_ice(:,:) * e2u(:,:) 
     177         zvdx(:,:) = v_ice(:,:) * e1v(:,:) 
     178          
     179         ! define velocity for advection: u*grad(H) 
     180         DO jj = 2, jpjm1 
     181            DO ji = fs_2, fs_jpim1 
     182               IF    ( u_ice(ji,jj) * u_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp 
     183               ELSEIF( u_ice(ji,jj)                  >  0._wp ) THEN   ;   zcu_box(ji,jj) = u_ice(ji-1,jj) 
     184               ELSE                                                    ;   zcu_box(ji,jj) = u_ice(ji  ,jj) 
     185               ENDIF 
     186                
     187               IF    ( v_ice(ji,jj) * v_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp 
     188               ELSEIF( v_ice(ji,jj)                  >  0._wp ) THEN   ;   zcv_box(ji,jj) = v_ice(ji,jj-1) 
     189               ELSE                                                    ;   zcv_box(ji,jj) = v_ice(ji,jj  ) 
     190               ENDIF 
     191            END DO 
     192         END DO 
     193          
     194         ! advection 
     195         DO jt = 1, initad 
     196            CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, ato_i(:,:) )       ! Open water area  
     197            DO jl = 1, jpl 
     198               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_i(:,:,jl) )      ! Ice area 
     199               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_i(:,:,jl) )      ! Ice  volume 
     200               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, smv_i(:,:,jl) )    ! Salt content 
     201               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, oa_i (:,:,jl) )    ! Age content 
     202               DO jk = 1, nlay_i 
     203                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_i(:,:,jk,jl) )   ! Ice  heat content 
     204               END DO 
     205               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) )      ! Snow volume 
     206               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) )    ! Snow heat content 
     207            END DO 
     208         END DO 
     209         ! 
     210         at_i(:,:) = a_i(:,:,1)      ! total ice fraction 
     211         DO jl = 2, jpl 
     212            at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
     213         END DO 
     214         ! 
     215         CALL wrk_dealloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box ) 
     216          
     217                       !=============================! 
     218      CASE ( -1 )      !==     Prather scheme      ==!                    
     219                       !=============================! 
     220 
     221         CALL wrk_alloc( jpi,jpj,            zarea ) 
     222         CALL wrk_alloc( jpi,jpj,1,          z0opw ) 
     223         CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
     224         CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei ) 
     225          
     226         IF( kt == nit000 .AND. lwp ) THEN 
     227            WRITE(numout,*)'' 
     228            WRITE(numout,*)'lim_adv_xy : Prather advection scheme' 
     229            WRITE(numout,*)'~~~~~~~~~~~' 
     230         ENDIF 
     231          
     232         zarea(:,:) = e12t(:,:) 
     233          
     234         !------------------------- 
     235         ! transported fields                                         
     236         !------------------------- 
     237         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area  
    182238         DO jl = 1, jpl 
    183             WHERE( v_i(:,:,jl) /= 0._wp )   ;   z1_v(:,:) = 1._wp / v_i(:,:,jl) 
    184             ELSEWHERE                       ;   z1_v(:,:) = 0._wp 
    185             END WHERE 
    186             ! 
    187             lcon = .TRUE.  
    188 !!!               lcon = .false.  
    189             CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, a_i(:,:,jl), a_i(:,:,jl) )                      ! Ice area 
    190             CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, v_i(:,:,jl), v_i(:,:,jl), zu_trp, zv_trp )      ! Ice  volume 
    191             ! 
    192             lcon = .FALSE.  
    193             ze(:,:) = smv_i(:,:,jl) * z1_v(:,:) 
    194             CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, smv_i(:,:,jl) )                           ! Salt content 
    195             ! 
    196 !!!check that               ze(:,:) = oa_i (:,:,jl) * z1_v(:,:) 
    197 !!!check that               CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, oa_i (:,:,jl) )                           ! Age content 
    198             ! 
    199             zu_trp(:,:) = zu_trp(:,:) * r1_nlay_i 
    200             zv_trp(:,:) = zv_trp(:,:) * r1_nlay_i 
    201             z1_v  (:,:) = z1_v  (:,:) * REAL( nlay_i, wp ) 
     239            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume 
     240            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume 
     241            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area 
     242            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content 
     243            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content 
     244            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content 
    202245            DO jk = 1, nlay_i 
    203                ze (:,:) = e_i(:,:,jk,jl) * z1_v(:,:) 
    204                CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, e_i(:,:,jk,jl) )                                  ! Ice  heat content 
    205             END DO 
    206             ! 
    207             WHERE( v_s(:,:,jl) /= 0._wp )   ;   z1_v(:,:) = 1._wp / v_s(:,:,jl) 
    208             ELSEWHERE                       ;   z1_v(:,:) = 0._wp 
    209             END WHERE 
    210             ! 
    211             lcon = .TRUE.  
    212 !!!               lcon = .false.  
    213             CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, v_s(:,:,jl), v_s(:,:,jl), zu_trp, zv_trp )      ! Snow volume 
    214             ! 
    215             lcon = .FALSE.  
    216             ze (:,:) = e_s(:,:,1,jl) * z1_v(:,:) 
    217             CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, e_s(:,:,1,jl) )                                     ! Snow heat content 
    218             ! 
    219          END DO 
    220       END DO 
    221       ! 
    222       at_i(:,:) = a_i(:,:,1)      ! total ice fraction 
    223       DO jl = 2, jpl 
    224          at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    225       END DO 
    226       ! 
    227       CALL wrk_dealloc( jpi,jpj, ze, zu_trp, zv_trp, z1_v, zudy, zvdx ) 
    228          
    229 #else 
    230       !=============================! 
    231       !==      Prather scheme     ==! 
    232       !=============================! 
    233  
    234       CALL wrk_alloc( jpi,jpj,            zarea ) 
    235       CALL wrk_alloc( jpi,jpj,1,          z0opw ) 
    236       CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
    237       CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei ) 
    238        
    239       IF( kt == nit000 .AND. lwp ) THEN 
    240          WRITE(numout,*)'' 
    241          WRITE(numout,*)'lim_adv_xy : Prather advection scheme' 
    242          WRITE(numout,*)'~~~~~~~~~~~' 
    243       ENDIF 
    244  
    245       zarea(:,:) = e12t(:,:) 
    246        
    247       !------------------------- 
    248       ! transported fields                                         
    249       !------------------------- 
    250       z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area  
    251       DO jl = 1, jpl 
    252          z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume 
    253          z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume 
    254          z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area 
    255          z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content 
    256          z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content 
    257          z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content 
    258          DO jk = 1, nlay_i 
    259             z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content 
    260          END DO 
    261       END DO 
    262  
    263        
    264       IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
    265          DO jt = 1, initad 
    266             CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
    267                &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    268             CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
    269                &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    270             DO jl = 1, jpl 
    271                CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    272                   &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    273                CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
    274                   &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    275                CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    276                   &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    277                CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    278                   &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    279                CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    280                   &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    281                CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    282                   &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    283                CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---      
    284                   &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    285                CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
    286                   &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    287                CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations --- 
    288                   &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    289                CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &  
    290                   &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    291                CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents --- 
    292                   &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    293                CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   & 
    294                   &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    295                DO jk = 1, nlay_i                                                                !--- ice heat contents --- 
    296                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    297                      &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    298                      &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    299                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    300                      &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    301                      &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     246               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content 
     247            END DO 
     248         END DO 
     249          
     250          
     251         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
     252            DO jt = 1, initad 
     253               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
     254                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     255               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
     256                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     257               DO jl = 1, jpl 
     258                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     259                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     260                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
     261                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     262                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     263                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     264                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
     265                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     266                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     267                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     268                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
     269                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     270                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---      
     271                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     272                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
     273                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     274                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations --- 
     275                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     276                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &  
     277                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     278                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents --- 
     279                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     280                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   & 
     281                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     282                  DO jk = 1, nlay_i                                                                !--- ice heat contents --- 
     283                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     284                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
     285                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     286                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     287                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
     288                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     289                  END DO 
    302290               END DO 
    303291            END DO 
    304          END DO 
    305       ELSE 
    306          DO jt = 1, initad 
    307             CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
    308                &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    309             CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
    310                &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    311             DO jl = 1, jpl 
    312                CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    313                   &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    314                CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
    315                   &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    316                CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    317                   &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    318                CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    319                   &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    320                CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    321                   &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    322                CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    323                   &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    324                CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
    325                   &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    326                CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
    327                   &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    328                CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
    329                   &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    330                CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
    331                   &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    332                CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
    333                   &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    334                CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   & 
    335                   &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    336                DO jk = 1, nlay_i                                                           !--- ice heat contents --- 
    337                   CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    338                      &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    339                      &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    340                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    341                      &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    342                      &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     292         ELSE 
     293            DO jt = 1, initad 
     294               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
     295                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     296               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
     297                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     298               DO jl = 1, jpl 
     299                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     300                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     301                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
     302                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     303                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     304                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     305                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
     306                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     307                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     308                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     309                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
     310                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     311                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
     312                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     313                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
     314                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     315                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
     316                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     317                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
     318                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     319                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
     320                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     321                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   & 
     322                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     323                  DO jk = 1, nlay_i                                                           !--- ice heat contents --- 
     324                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     325                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
     326                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     327                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     328                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
     329                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     330                  END DO 
    343331               END DO 
    344332            END DO 
    345          END DO 
    346       ENDIF 
    347        
    348       !------------------------------------------- 
    349       ! Recover the properties from their contents 
    350       !------------------------------------------- 
    351       ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:) 
    352       DO jl = 1, jpl 
    353          v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:) 
    354          v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:) 
    355          smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:) 
    356          oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:) 
    357          a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:) 
    358          e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:) 
    359          DO jk = 1, nlay_i 
    360             e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:) 
    361          END DO 
    362       END DO 
    363        
    364       at_i(:,:) = a_i(:,:,1)      ! total ice fraction 
    365       DO jl = 2, jpl 
    366          at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    367       END DO 
    368  
    369       CALL wrk_dealloc( jpi,jpj,            zarea ) 
    370       CALL wrk_dealloc( jpi,jpj,1,          z0opw ) 
    371       CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
    372       CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei ) 
    373        
    374 #endif 
    375  
     333         ENDIF 
     334          
     335         !------------------------------------------- 
     336         ! Recover the properties from their contents 
     337         !------------------------------------------- 
     338         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:) 
     339         DO jl = 1, jpl 
     340            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:) 
     341            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:) 
     342            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:) 
     343            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:) 
     344            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:) 
     345            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:) 
     346            DO jk = 1, nlay_i 
     347               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:) 
     348            END DO 
     349         END DO 
     350          
     351         at_i(:,:) = a_i(:,:,1)      ! total ice fraction 
     352         DO jl = 2, jpl 
     353            at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
     354         END DO 
     355          
     356         CALL wrk_dealloc( jpi,jpj,            zarea ) 
     357         CALL wrk_dealloc( jpi,jpj,1,          z0opw ) 
     358         CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
     359         CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei ) 
     360 
     361      END SELECT 
     362       
    376363      !------------------------------! 
    377364      ! Diffusion of Ice fields                   
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r6970 r6989  
    321321         &     'use more ocean levels or less ice/snow layers/categories.' ) 
    322322      ! 
     323      CALL lim_dyn_init                ! set ice dynamics parameters 
     324      ! 
    323325      CALL lim_itd_init                ! ice thickness distribution initialization 
    324326      ! 
Note: See TracChangeset for help on using the changeset viewer.