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 8504 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90 – NEMO

Ignore:
Timestamp:
2017-09-06T17:47:28+02:00 (7 years ago)
Author:
clem
Message:

changes in style - part4 - clarify ice advection routines

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90

    r8500 r8504  
    3636   PUBLIC   ice_adv    ! called by icestp 
    3737 
    38    INTEGER ::   ncfl   ! number of ice time step with CFL>1/2   
    39  
    4038   !! * Substitution 
    4139#  include "vectopt_loop_substitute.h90" 
     
    5149      !!                   ***  ROUTINE ice_adv *** 
    5250      !!                     
    53       !! ** purpose : advection/diffusion process of sea ice 
     51      !! ** purpose : advection of sea ice 
    5452      !! 
    5553      !! ** method  : variables included in the process are scalar,    
     
    7270      REAL(wp), DIMENSION(jpi,jpj)           ::   zatold, zeiold, zesold, zsmvold  
    7371      REAL(wp), DIMENSION(jpi,jpj,jpl)       ::   zhimax, zviold, zvsold 
    74       ! --- ultimate macho only --- ! 
    75       REAL(wp) ::   zdt 
    76       REAL(wp), ALLOCATABLE, DIMENSION(:,:)      ::   zudy, zvdx, zcu_box, zcv_box 
    77       ! --- prather only --- ! 
    78       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zarea 
    79       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0opw 
    80       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi 
    81       ! MV MP 2016 
    82       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp 
    83       ! END MV MP 2016 
    84       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei 
    8572      !!--------------------------------------------------------------------- 
    8673      IF( nn_timing == 1 )  CALL timing_start('iceadv') 
    8774  
    8875      IF( kt == nit000 .AND. lwp ) THEN 
    89          WRITE(numout,*)'' 
    90          WRITE(numout,*)'iceadv : sea-ice advection' 
    91          WRITE(numout,*)'~~~~~~~' 
    92          ncfl = 0                ! nb of time step with CFL > 1/2 
     76         WRITE(numout,*) 
     77         WRITE(numout,*) 'iceadv: sea-ice advection' 
     78         WRITE(numout,*) '~~~~~~' 
    9379      ENDIF 
    9480       
    9581      CALL ice_var_agg( 1 ) ! integrated values + ato_i 
    96  
    97       !-------------------------------------! 
    98       !   Advection of sea ice properties   ! 
    99       !-------------------------------------! 
    10082 
    10183      ! conservation test 
     
    10385       
    10486      ! store old values for diag 
    105       zviold = v_i 
    106       zvsold = v_s 
    107       zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 ) 
    108       zeiold (:,:) = et_i 
    109       zesold (:,:) = et_s  
    110  
    111       !--- Thickness correction init. --- ! 
     87      zviold (:,:,:) = v_i(:,:,:) 
     88      zvsold (:,:,:) = v_s(:,:,:) 
     89      zsmvold(:,:)   = SUM( smv_i(:,:,:), dim=3 ) 
     90      zeiold (:,:)   = et_i(:,:) 
     91      zesold (:,:)   = et_s(:,:) 
     92 
     93      ! Thickness correction init. 
    11294      zatold(:,:) = at_i 
    11395      WHERE( a_i(:,:,:) >= epsi20 ) 
     
    119101      END WHERE 
    120102          
    121       ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- ! 
     103      ! Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick 
    122104      zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 
    123105      DO jl = 1, jpl 
     
    130112      END DO 
    131113      CALL lbc_lnk( zhimax(:,:,:), 'T', 1. ) 
    132           
    133       ! --- If ice drift field is too fast, use an appropriate time step for advection --- !         
    134       zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )    ! CFL test for stability 
    135       zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
    136       IF( lk_mpp )   CALL mpp_max( zcfl ) 
    137        
    138       IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
    139       ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
    140       ENDIF 
    141        
    142 !!      IF( zcfl > 0.5_wp .AND. lwp ) THEN 
    143 !!         ncfl = ncfl + 1 
    144 !!         IF( ncfl > 0 ) THEN    
    145 !!            WRITE(cltmp,'(i6.1)') ncfl 
    146 !!            CALL ctl_warn( 'ice_adv: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 
    147 !!         ENDIF 
    148 !!      ENDIF 
    149  
     114       
     115      !---------- 
     116      ! Advection 
     117      !---------- 
    150118      SELECT CASE ( nn_limadv ) 
    151           
    152                        !=============================! 
    153       CASE ( 0 )       !==  Ultimate-MACHO scheme  ==!                    
    154                        !=============================! 
    155        
    156          ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 
    157        
    158          IF( kt == nit000 .AND. lwp ) THEN 
    159             WRITE(numout,*)'' 
    160             WRITE(numout,*)'ice_adv_umx : Ultimate-MACHO advection scheme' 
    161             WRITE(numout,*)'~~~~~~~~~~~' 
    162          ENDIF 
    163          ! 
    164          zdt = rdt_ice / REAL(initad) 
    165           
    166          ! transport 
    167          zudy(:,:) = u_ice(:,:) * e2u(:,:) 
    168          zvdx(:,:) = v_ice(:,:) * e1v(:,:) 
    169           
    170          ! define velocity for advection: u*grad(H) 
    171          DO jj = 2, jpjm1 
    172             DO ji = fs_2, fs_jpim1 
    173                IF    ( u_ice(ji,jj) * u_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp 
    174                ELSEIF( u_ice(ji,jj)                  >  0._wp ) THEN   ;   zcu_box(ji,jj) = u_ice(ji-1,jj) 
    175                ELSE                                                    ;   zcu_box(ji,jj) = u_ice(ji  ,jj) 
    176                ENDIF 
    177                 
    178                IF    ( v_ice(ji,jj) * v_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp 
    179                ELSEIF( v_ice(ji,jj)                  >  0._wp ) THEN   ;   zcv_box(ji,jj) = v_ice(ji,jj-1) 
    180                ELSE                                                    ;   zcv_box(ji,jj) = v_ice(ji,jj  ) 
    181                ENDIF 
    182             END DO 
    183          END DO 
    184           
    185          ! advection 
    186          DO jt = 1, initad 
    187             CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, ato_i(:,:) )       ! Open water area  
    188             DO jl = 1, jpl 
    189                CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_i(:,:,jl) )      ! Ice area 
    190                CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_i(:,:,jl) )      ! Ice  volume 
    191                CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, smv_i(:,:,jl) )    ! Salt content 
    192                CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, oa_i (:,:,jl) )    ! Age content 
    193                DO jk = 1, nlay_i 
    194                   CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_i(:,:,jk,jl) )   ! Ice  heat content 
    195                END DO 
    196                CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) )      ! Snow volume 
    197                CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) )    ! Snow heat content 
    198                ! MV MP 2016 
    199                IF ( nn_pnd_scheme > 0 ) THEN 
    200                   CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_ip(:,:,jl) )  ! Melt pond fraction 
    201                   CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_ip(:,:,jl) )  ! Melt pond volume 
    202                ENDIF 
    203                ! END MV MP 2016 
    204             END DO 
    205          END DO 
    206          ! 
    207          at_i(:,:) = a_i(:,:,1)      ! total ice fraction 
    208          DO jl = 2, jpl 
    209             at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    210          END DO 
    211          ! 
    212          DEALLOCATE( zudy, zvdx, zcu_box, zcv_box ) 
    213           
    214                        !=============================! 
    215       CASE ( -1 )      !==     Prather scheme      ==!                    
    216                        !=============================! 
    217  
    218          ALLOCATE( zarea(jpi,jpj)     , z0opw(jpi,jpj, 1 ) ,                                           & 
    219             &      z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , z0ai(jpi,jpj,jpl) , z0es(jpi,jpj,jpl) ,   & 
    220             &      z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap(jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   & 
    221             &      z0ei (jpi,jpj,nlay_i,jpl)                                                           ) 
    222           
    223          IF( kt == nit000 .AND. lwp ) THEN 
    224             WRITE(numout,*)'' 
    225             WRITE(numout,*)'ice_adv_xy : Prather advection scheme' 
    226             WRITE(numout,*)'~~~~~~~~~~~' 
    227          ENDIF 
    228           
    229          zarea(:,:) = e1e2t(:,:) 
    230           
    231          !------------------------- 
    232          ! transported fields                                         
    233          !------------------------- 
    234          z0opw(:,:,1) = ato_i(:,:) * e1e2t(:,:)             ! Open water area  
    235          DO jl = 1, jpl 
    236             z0snw(:,:,jl) = v_s  (:,:,  jl) * e1e2t(:,:)  ! Snow volume 
    237             z0ice(:,:,jl) = v_i  (:,:,  jl) * e1e2t(:,:)  ! Ice  volume 
    238             z0ai (:,:,jl) = a_i  (:,:,  jl) * e1e2t(:,:)  ! Ice area 
    239             z0smi(:,:,jl) = smv_i(:,:,  jl) * e1e2t(:,:)  ! Salt content 
    240             z0oi (:,:,jl) = oa_i (:,:,  jl) * e1e2t(:,:)  ! Age content 
    241             z0es (:,:,jl) = e_s  (:,:,1,jl) * e1e2t(:,:)  ! Snow heat content 
    242             DO jk = 1, nlay_i 
    243                z0ei(:,:,jk,jl) = e_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    244             END DO 
    245             ! MV MP 2016 
    246             IF ( nn_pnd_scheme > 0 ) THEN 
    247                z0ap(:,:,jl)  = a_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 
    248                z0vp(:,:,jl)  = v_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 
    249             ENDIF 
    250             ! END MV MP 2016 
    251          END DO 
    252  
    253          IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
    254             DO jt = 1, initad 
    255                CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
    256                   &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    257                CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
    258                   &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    259                DO jl = 1, jpl 
    260                   CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    261                      &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    262                   CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
    263                      &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    264                   CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    265                      &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    266                   CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    267                      &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    268                   CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    269                      &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    270                   CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    271                      &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    272                   CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---      
    273                      &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    274                   CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
    275                      &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    276                   CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations --- 
    277                      &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    278                   CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &  
    279                      &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    280                   CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents --- 
    281                      &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    282                   CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   & 
    283                      &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    284                   DO jk = 1, nlay_i                                                                !--- ice heat contents --- 
    285                      CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    286                         &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    287                         &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    288                      CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    289                         &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    290                         &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    291                   END DO 
    292                   ! MV MP 2016 
    293                   CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction -- 
    294                      &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    295                   CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &  
    296                      &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    297                   CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   -- 
    298                      &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    299                   CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &  
    300                      &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    301                   ! END MV MP 2016 
    302                END DO 
    303             END DO 
    304          ELSE 
    305             DO jt = 1, initad 
    306                CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
    307                   &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    308                CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
    309                   &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    310                DO jl = 1, jpl 
    311                   CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    312                      &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    313                   CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
    314                      &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    315                   CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    316                      &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    317                   CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    318                      &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    319                   CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    320                      &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    321                   CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    322                      &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    323                   CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
    324                      &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    325                   CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
    326                      &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    327                   CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
    328                      &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    329                   CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
    330                      &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    331                   CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
    332                      &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    333                   CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   & 
    334                      &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    335                   DO jk = 1, nlay_i                                                           !--- ice heat contents --- 
    336                      CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    337                         &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    338                         &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    339                      CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    340                         &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    341                         &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    342                   END DO 
    343                   ! MV MP 2016 
    344                   IF ( nn_pnd_scheme > 0 ) THEN 
    345                      CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction --- 
    346                      &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    347                      CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
    348                      &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    349                      CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   --- 
    350                      &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    351                      CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
    352                      &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    353                   ENDIF 
    354                   ! END MV MP 2016 
    355                END DO 
    356             END DO 
    357          ENDIF 
    358  
    359          !------------------------------------------- 
    360          ! Recover the properties from their contents 
    361          !------------------------------------------- 
    362          ato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) 
    363          DO jl = 1, jpl 
    364             v_i  (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) 
    365             v_s  (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) 
    366             smv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) 
    367             oa_i (:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) 
    368             a_i  (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) 
    369             e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:) 
    370             DO jk = 1, nlay_i 
    371                e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) 
    372             END DO 
    373             ! MV MP 2016 
    374             IF ( nn_pnd_scheme > 0 ) THEN 
    375                a_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:) 
    376                v_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:) 
    377             ENDIF 
    378             ! END MV MP 2016 
    379          END DO 
    380          ! 
    381          at_i(:,:) = a_i(:,:,1)      ! total ice fraction 
    382          DO jl = 2, jpl 
    383             at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    384          END DO 
    385          ! 
    386          DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0es , z0smi , z0oi , z0ap , z0vp , z0ei ) 
    387          ! 
     119      CASE ( 0 )                    !--  Ultimate-MACHO scheme 
     120         CALL ice_adv_umx( kt, u_ice, v_ice,  & 
     121            &              ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     122          
     123      CASE ( -1 )                   !-- Prather scheme 
     124         CALL ice_adv_prather( kt, u_ice, v_ice,  & 
     125            &                  ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     126          
    388127      END SELECT 
    389        
    390       ! --- diags --- 
     128 
     129      ! total ice fraction 
     130      at_i(:,:) = a_i(:,:,1) 
     131      DO jl = 2, jpl 
     132         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
     133      END DO 
     134       
     135      !------------ 
     136      ! diagnostics 
     137      !------------ 
    391138      DO jj = 1, jpj 
    392139         DO ji = 1, jpi 
     
    404151      IF( iom_use('destrp') )   CALL iom_put( "destrp" , diag_trp_es         )         ! advected snw enthalpy (W/m2) 
    405152       
     153      !-------------------------------------- 
     154      ! Thickness correction in case too high 
     155      !-------------------------------------- 
    406156      IF( nn_limdyn == 2 ) THEN 
    407157         ! 
    408          CALL ice_var_zapsmall      !--- zap small areas ---! 
     158         CALL ice_var_zapsmall                       !-- zap small areas 
    409159         ! 
    410          DO jl = 1, jpl             !--- Thickness correction in case too high --- ! 
     160         DO jl = 1, jpl 
    411161            DO jj = 1, jpj 
    412162               DO ji = 1, jpi 
    413                   IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
     163                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to zhimax 
    414164                     ! 
    415165                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 
     
    428178         END DO 
    429179                   
    430          WHERE( ht_i(:,:,jpl) > hi_max(jpl) )             !--- bound ht_i to hi_max (99 m) 
     180         WHERE( ht_i(:,:,jpl) > hi_max(jpl) )        !-- bound ht_i to hi_max (99 m) 
    431181            ht_i(:,:,jpl) = hi_max(jpl) 
    432182            a_i (:,:,jpl) = v_i(:,:,jpl) / hi_max(jpl) 
    433183         END WHERE 
    434184 
    435          IF ( nn_pnd_scheme > 0 ) THEN                    !--- correct pond fraction to avoid a_ip > a_i 
     185         IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i 
    436186            WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
    437187         ENDIF 
     
    442192      ! Impose a_i < amax if no ridging/rafting or in mono-category 
    443193      !------------------------------------------------------------ 
    444       ! 
    445       IF( l_piling ) THEN ! simple conservative piling, comparable with 1-cat models 
     194      IF( l_piling ) THEN                            !-- simple conservative piling, comparable with 1-cat models 
    446195         at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    447196         DO jl = 1, jpl 
     
    452201      ENDIF 
    453202       
    454       ! --- agglomerate variables ----------------- 
     203      ! agglomerate variables  
    455204      vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 
    456205      vt_s(:,:) = SUM( v_s(:,:,:), dim=3 ) 
     
    464213      ! END MP 2016 
    465214       
    466       ! --- open water = 1 if at_i=0 -------------------------------- 
     215      ! open water = 1 if at_i=0  
    467216      WHERE( at_i == 0._wp ) ato_i = 1._wp  
    468217       
     
    470219      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    471220         
    472       ! ------------------------------------------------- 
     221      ! -------------- 
    473222      ! control prints 
    474       ! ------------------------------------------------- 
     223      ! -------------- 
    475224      IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 
    476225      ! 
Note: See TracChangeset for help on using the changeset viewer.