Changeset 8504


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

changes in style - part4 - clarify ice advection routines

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
Files:
3 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      ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90

    r8486 r8504  
    1212   !!   'key_lim3'                                       LIM3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    !!   ice_adv_x     : advection of sea ice on x axis 
    15    !!   ice_adv_y     : advection of sea ice on y axis 
     14   !!   ice_adv_prather     : advection of sea ice using Prather scheme 
    1615   !!---------------------------------------------------------------------- 
    1716   USE dom_oce        ! ocean domain 
    1817   USE ice            ! sea-ice variables 
     18   USE sbc_oce , ONLY : nn_fsbc   ! frequency of sea-ice call 
    1919   ! 
    2020   USE lbclnk         ! lateral boundary condition - MPP exchanges 
     
    2727   PRIVATE 
    2828 
    29    PUBLIC   ice_adv_x   ! called by iceadv 
    30    PUBLIC   ice_adv_y   ! called by iceadv 
     29   PUBLIC   ice_adv_prather   ! called by iceadv 
    3130 
    3231   !! * Substitutions 
     
    3938CONTAINS 
    4039 
    41    SUBROUTINE ice_adv_x( pdf, put , pcrh, psm , ps0 ,   & 
    42       &                  psx, psxx, psy , psyy, psxy ) 
     40   SUBROUTINE ice_adv_prather( kt, pu_ice, pv_ice,  & 
     41      &                        pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    4342      !!---------------------------------------------------------------------- 
    44       !!                **  routine ice_adv_x  ** 
     43      !!                **  routine ice_adv_prather  ** 
    4544      !!   
    4645      !! ** purpose :   Computes and adds the advection trend to sea-ice 
    47       !!              variable on i-axis 
    4846      !! 
    4947      !! ** method  :   Uses Prather second order scheme that advects tracers 
    50       !!              but also theirquadratic forms. The method preserves 
    51       !!              tracer structures by conserving second order moments. 
     48      !!                but also their quadratic forms. The method preserves 
     49      !!                tracer structures by conserving second order moments. 
    5250      !! 
    5351      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681. 
    5452      !!---------------------------------------------------------------------- 
     53      INTEGER                     , INTENT(in   ) ::   kt         ! time step 
     54      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
     55      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity 
     56      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
     57      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
     58      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume 
     59      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psmv_i     ! salt content 
     60      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content 
     61      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration 
     62      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
     63      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     64      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
     65      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     66      ! 
     67      INTEGER  ::   jk, jl, jt              ! dummy loop indices 
     68      INTEGER  ::   initad                  ! number of sub-timestep for the advection 
     69      REAL(wp) ::   zcfl , zusnit           !   -      - 
     70      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zarea 
     71      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0opw 
     72      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi 
     73      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp 
     74      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei 
     75      !!---------------------------------------------------------------------- 
     76      ! 
     77      IF( kt == nit000 .AND. lwp ) THEN 
     78         WRITE(numout,*) 
     79         WRITE(numout,*) 'ice_adv_prather: Prather advection scheme' 
     80         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
     81      ENDIF 
     82      ! 
     83      ALLOCATE( zarea(jpi,jpj)     , z0opw(jpi,jpj, 1 ) ,                                           & 
     84         &      z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , z0ai(jpi,jpj,jpl) , z0es(jpi,jpj,jpl) ,   & 
     85         &      z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap(jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   & 
     86         &      z0ei (jpi,jpj,nlay_i,jpl)                                                           ) 
     87      ! 
     88      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !         
     89      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 
     90      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
     91      IF( lk_mpp )   CALL mpp_max( zcfl ) 
     92       
     93      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
     94      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
     95      ENDIF 
     96       
     97      zarea(:,:) = e1e2t(:,:) 
     98      !------------------------- 
     99      ! transported fields                                         
     100      !------------------------- 
     101      z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:)              ! Open water area  
     102      DO jl = 1, jpl 
     103         z0snw(:,:,jl) = pv_s  (:,:,  jl) * e1e2t(:,:)     ! Snow volume 
     104         z0ice(:,:,jl) = pv_i  (:,:,  jl) * e1e2t(:,:)     ! Ice  volume 
     105         z0ai (:,:,jl) = pa_i  (:,:,  jl) * e1e2t(:,:)     ! Ice area 
     106         z0smi(:,:,jl) = psmv_i(:,:,  jl) * e1e2t(:,:)     ! Salt content 
     107         z0oi (:,:,jl) = poa_i (:,:,  jl) * e1e2t(:,:)     ! Age content 
     108         z0es (:,:,jl) = pe_s  (:,:,1,jl) * e1e2t(:,:)     ! Snow heat content 
     109         DO jk = 1, nlay_i 
     110            z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
     111         END DO 
     112         IF ( nn_pnd_scheme > 0 ) THEN 
     113            z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
     114            z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
     115         ENDIF 
     116      END DO 
     117 
     118      !                                                    !--------------------------------------------! 
     119      IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
     120         !                                                 !--------------------------------------------! 
     121         DO jt = 1, initad 
     122            CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
     123               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     124            CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
     125               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     126            DO jl = 1, jpl 
     127               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     128                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     129               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
     130                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     131               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     132                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     133               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
     134                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     135               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     136                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     137               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
     138                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     139               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---      
     140                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     141               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
     142                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     143               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations --- 
     144                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     145               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &  
     146                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     147               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents --- 
     148                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     149               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   & 
     150                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     151               DO jk = 1, nlay_i                                                               !--- ice heat contents --- 
     152                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     153                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
     154                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     155                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     156                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
     157                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     158               END DO 
     159               IF ( nn_pnd_scheme > 0 ) THEN 
     160                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction -- 
     161                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
     162                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &  
     163                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
     164                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   -- 
     165                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
     166                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &  
     167                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
     168               ENDIF 
     169            END DO 
     170         END DO 
     171      !                                                    !--------------------------------------------! 
     172      ELSE                                                 !== even ice time step:  adv_y then adv_x  ==! 
     173         !                                                 !--------------------------------------------! 
     174         DO jt = 1, initad 
     175            CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
     176               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     177            CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
     178               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     179            DO jl = 1, jpl 
     180               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     181                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     182               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
     183                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     184               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     185                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     186               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
     187                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     188               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     189                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     190               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
     191                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     192               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
     193                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     194               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
     195                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     196               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
     197                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     198               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
     199                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     200               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
     201                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     202               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   & 
     203                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     204               DO jk = 1, nlay_i                                                             !--- ice heat contents --- 
     205                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     206                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
     207                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     208                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     209                     &                                      sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
     210                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     211               END DO 
     212               IF ( nn_pnd_scheme > 0 ) THEN 
     213                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction --- 
     214                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
     215                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
     216                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
     217                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   --- 
     218                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
     219                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
     220                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
     221               ENDIF 
     222            END DO 
     223         END DO 
     224      ENDIF 
     225 
     226      !------------------------------------------- 
     227      ! Recover the properties from their contents 
     228      !------------------------------------------- 
     229      pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) 
     230      DO jl = 1, jpl 
     231         pv_i  (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) 
     232         pv_s  (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) 
     233         psmv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) 
     234         poa_i (:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) 
     235         pa_i  (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) 
     236         pe_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:) 
     237         DO jk = 1, nlay_i 
     238            pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) 
     239         END DO 
     240         ! MV MP 2016 
     241         IF ( nn_pnd_scheme > 0 ) THEN 
     242            pa_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:) 
     243            pv_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:) 
     244         ENDIF 
     245         ! END MV MP 2016 
     246      END DO 
     247      ! 
     248      DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0es , z0smi , z0oi , z0ap , z0vp , z0ei ) 
     249      ! 
     250   END SUBROUTINE ice_adv_prather 
     251    
     252   SUBROUTINE adv_x( pdf, put , pcrh, psm , ps0 ,   & 
     253      &              psx, psxx, psy , psyy, psxy ) 
     254      !!---------------------------------------------------------------------- 
     255      !!                **  routine adv_x  ** 
     256      !!   
     257      !! ** purpose :   Computes and adds the advection trend to sea-ice 
     258      !!                variable on x axis 
     259      !!---------------------------------------------------------------------- 
    55260      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
    56       REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call ice_adv_x then ice_adv_y (=1) or the opposite (=0) 
     261      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
    57262      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s] 
    58263      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area 
     
    209414 
    210415      IF(ln_ctl) THEN 
    211          CALL prt_ctl(tab2d_1=psm  , clinfo1=' ice_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ') 
    212          CALL prt_ctl(tab2d_1=psx  , clinfo1=' ice_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ') 
    213          CALL prt_ctl(tab2d_1=psy  , clinfo1=' ice_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ') 
    214          CALL prt_ctl(tab2d_1=psxy , clinfo1=' ice_adv_x: psxy :') 
     416         CALL prt_ctl(tab2d_1=psm  , clinfo1=' adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ') 
     417         CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ') 
     418         CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ') 
     419         CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_x: psxy :') 
    215420      ENDIF 
    216421      ! 
    217    END SUBROUTINE ice_adv_x 
    218  
    219  
    220    SUBROUTINE ice_adv_y( pdf, pvt , pcrh, psm , ps0 ,   & 
    221       &                  psx, psxx, psy , psyy, psxy ) 
     422   END SUBROUTINE adv_x 
     423 
     424 
     425   SUBROUTINE adv_y( pdf, pvt , pcrh, psm , ps0 ,   & 
     426      &              psx, psxx, psy , psyy, psxy ) 
    222427      !!--------------------------------------------------------------------- 
    223       !!                **  routine ice_adv_y  ** 
     428      !!                **  routine adv_y  ** 
    224429      !!             
    225430      !! ** purpose :   Computes and adds the advection trend to sea-ice  
    226       !!              variable on y axis 
    227       !! 
    228       !! ** method  :   Uses Prather second order scheme that advects tracers 
    229       !!              but also their quadratic forms. The method preserves  
    230       !!              tracer structures by conserving second order moments. 
    231       !!  
    232       !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681. 
     431      !!                variable on y axis 
    233432      !!--------------------------------------------------------------------- 
    234433      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
    235       REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call ice_adv_x then ice_adv_y (=1) or the opposite (=0) 
     434      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
    236435      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s] 
    237436      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area 
     
    389588 
    390589      IF(ln_ctl) THEN 
    391          CALL prt_ctl(tab2d_1=psm  , clinfo1=' ice_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ') 
    392          CALL prt_ctl(tab2d_1=psx  , clinfo1=' ice_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ') 
    393          CALL prt_ctl(tab2d_1=psy  , clinfo1=' ice_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ') 
    394          CALL prt_ctl(tab2d_1=psxy , clinfo1=' ice_adv_y: psxy :') 
     590         CALL prt_ctl(tab2d_1=psm  , clinfo1=' adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ') 
     591         CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ') 
     592         CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ') 
     593         CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_y: psxy :') 
    395594      ENDIF 
    396595      ! 
    397    END SUBROUTINE ice_adv_y 
     596   END SUBROUTINE adv_y 
    398597 
    399598#else 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90

    r8486 r8504  
    4343CONTAINS 
    4444 
    45    SUBROUTINE ice_adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc ) 
     45   SUBROUTINE ice_adv_umx( kt, pu_ice, pv_ice,  & 
     46      &                    pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    4647      !!---------------------------------------------------------------------- 
    4748      !!                  ***  ROUTINE ice_adv_umx  *** 
     49      !!  
     50      !! **  Purpose :   Compute the now trend due to total advection of  
     51      !!                 tracers and add it to the general trend of tracer equations 
     52      !!                 using an "Ultimate-Macho" scheme 
     53      !! 
     54      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
     55      !!---------------------------------------------------------------------- 
     56      INTEGER                     , INTENT(in   ) ::   kt         ! time step 
     57      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
     58      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity 
     59      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
     60      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
     61      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume 
     62      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psmv_i     ! salt content 
     63      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content 
     64      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration 
     65      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
     66      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     67      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
     68      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     69      ! 
     70      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
     71      INTEGER  ::   initad                  ! number of sub-timestep for the advection 
     72      REAL(wp) ::   zcfl , zusnit, zdt      !   -      - 
     73      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zudy, zvdx, zcu_box, zcv_box 
     74      !!---------------------------------------------------------------------- 
     75      ! 
     76      IF( kt == nit000 .AND. lwp ) THEN 
     77         WRITE(numout,*)  
     78         WRITE(numout,*) 'ice_adv_umx : Ultimate-MACHO advection scheme' 
     79         WRITE(numout,*) '~~~~~~~~~~~' 
     80      ENDIF 
     81      ! 
     82      ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 
     83      ! 
     84      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !         
     85      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 
     86      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
     87      IF( lk_mpp )   CALL mpp_max( zcfl ) 
     88 
     89      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
     90      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
     91      ENDIF 
     92 
     93      zdt = rdt_ice / REAL(initad) 
     94 
     95      ! --- transport --- ! 
     96      zudy(:,:) = pu_ice(:,:) * e2u(:,:) 
     97      zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 
     98 
     99      ! --- define velocity for advection: u*grad(H) --- ! 
     100      DO jj = 2, jpjm1 
     101         DO ji = fs_2, fs_jpim1 
     102            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp 
     103            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj) 
     104            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj) 
     105            ENDIF 
     106 
     107            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp 
     108            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1) 
     109            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  ) 
     110            ENDIF 
     111         END DO 
     112      END DO 
     113 
     114      !---------------! 
     115      !== advection ==! 
     116      !---------------! 
     117      DO jt = 1, initad 
     118         CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area  
     119         DO jl = 1, jpl 
     120            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area 
     121            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume 
     122            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) )       ! Salt content 
     123            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) )       ! Age content 
     124            DO jk = 1, nlay_i 
     125               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content 
     126            END DO 
     127            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume 
     128            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content 
     129            IF ( nn_pnd_scheme > 0 ) THEN 
     130               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction 
     131               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume 
     132            ENDIF 
     133         END DO 
     134      END DO 
     135      ! 
     136      DEALLOCATE( zudy, zvdx, zcu_box, zcv_box ) 
     137      ! 
     138   END SUBROUTINE ice_adv_umx 
     139    
     140   SUBROUTINE adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc ) 
     141      !!---------------------------------------------------------------------- 
     142      !!                  ***  ROUTINE adv_umx  *** 
    48143      !!  
    49144      !! **  Purpose :   Compute the now trend due to total advection of  
     
    146241      IF( nn_timing == 1 )  CALL timing_stop('ice_adv_umx') 
    147242      ! 
    148    END SUBROUTINE ice_adv_umx 
     243   END SUBROUTINE adv_umx 
    149244 
    150245 
Note: See TracChangeset for help on using the changeset viewer.